Stability preserving non-standard finite difference schemes for diabetes with tuberculosis infectious model

In this article, we study the non-standard finite difference (NSFD) schemes for the diabetes mellitus with tuberculosis infection. We have classified diabetesmellitus into two categories such as diabetes without complications and with complications and studied diabetes with tuberculosis infection. We have modelled diabetes with tuberculosis infection with linear and nonlinear model and studied the local andglobal asymptotic stability analysis using theNSFD scheme. ARTICLE HISTORY Received 26 October 2018 Accepted 6 May 2019


Introduction
Non-standard finite difference method (NSFD) was introduced by Mickens (1994). This method offers several advantages over finite difference methods. Whenever continuous dynamical systems have been converted into a discrete system, the properties of the continuous system are not transferred fully to the discrete systems specially in the case of large step size in the discrete systems. The traditional numerical methods like Euler methods and Runge-Kutta (RK) methods are found to be useful only for some values of step sizes and shown to be a different kind of behaviour not behaving exactly as its continuous counter part (Mickens, 1994). But if we use the NSFD scheme, the properties of the continuous dynamic systems can be preserved into its discrete counterpart (Raja Sekhara Rao, Venkata Ratnam, & Sita Rama Murthy, 2018). This scheme is very much useful in analysing the models arising from real life problems. Recently, a variety of NSFD schemes applied to the metapopulation model and studied for stability analysis by Dang and Hoang (2016, 2018a, 2018b. In our present study, we are using this NSFD scheme to analyse the stability of a biological model of Diabetes Tuberculosis association. Diabetes mellitus (DM) is a metabolic disorder and a noncommunicable disease developing many serious complications in the host body. DM can depress the ability of the immune system and thereby increasing the susceptibility to many infectious diseases including tuberculosis causing a serious threat to public health globally (Martinez & Kornfeld, 2014). Developing tuberculosis infection in diabetic individual is three times higher than non-diabetic individual. One third world population is actively or latently infected by tuberculosis infection (World Health Organization, 2013). A large number of people are also getting affected by diabetes worldwide, thus creating a joint burden of an increased chance of tuberculosis infection in diabetic people.

Pathophysiology of tuberculosis
Tuberculosis is an airborne disease caused by Mycobacterium tuberculosis, an acid-fast, aerobic and pathogenic bacillus, measuring 0.5 µm by 3 µm. This bacillus is protected with a unique cell structure which helps it in surviving under adverse conditions within the host body. Tuberculosis bacillus takes entry through an aerosol droplet to the alveolar passage of the exposed host (Smith, 2003). Later, this bacillus is phagocytosed by host macrophage (Schlesinger, 1993) and then resides within it. Under an immunocompromised condition in the host, tuberculosis bacillus comes out from the macrophage and establishes infection.

Diabetes with associated complications
Diabetes Mellitus is noncommunicable disease caused by metabolic disorder. Type I (insulin-dependent diabetes mellitus) and Type II (non insulin-dependent diabetes mellitus) are the main broad category of diabetes mellitus. Complications associated with diabetes may vary from accute to chronic form. The acute complications are mainly associated with hypoglycaemia, hyperosmolar hyperglycaemic nonketotic coma, diabetic ketoacidosis (Loewen & Haas, 1991). The chronic complications lead to mainly microvascular and macrovascular diseases. Autoimmune diseases are also very common with chronic complications due to weak immune systems. Diabetic nephropathy, neuropathy, and retinopathy have been categorized under microvascular complications (Fowler, 2008;Vithian & Hurel, 2010). Whereas macrovascular complications include mainly coronary artery disease, peripheral artery disease, stroke (Fowler, 2008). In Diabetic people with complications, the long term hyperglycaemic conditions lead to the debilitated immune system by weakening neutrophil function, antioxidant system. It also affects the humoral immunity (Casqueiro, Casqueiro, & Alves, 2012). Therefore, the immunocompromised state of diabetic people predisposes them to many infectious diseases such as respiratory tract infections, urinary tract infection, skin and soft tissue infection, gastrointestinal and liver infection (Casqueiro et al., 2012;Muller et al., 2005). In our present study, we analyse the development of active tuberculosis (respiratory tract infection) in diabetic people as it has been reported that diabetes increases the risk of developing active tuberculosis three times higher.

Diabetes as risk factor of tuberculosis infection
With the changing life style, less physical activities and consumption of foods with high calories have accelerated the incidence of diabetes worldwide. Almost 425 million people are affected with diabetes globally, whereas India alone had more than 72 million cases of diabetes in 2017 (International Diabetes Federation, 2018). New Tuberculosis incidence every year in India is also very high and estimated to be almost 2.79 million in 2016 (World Health Organization, 2017). Though association between diabetes and tuberculosis is still not clearly understood, but there are several studies and data which strongly supports this theory. There are many cellular factors and mechanism which predispose a diabetic individual to tuberculosis infection. As for example, phagocytes and T-cells are very much useful to protect the host against infection. In an experiment with a diabetic mice having two weeks tuberculosis infection, the phagocytosis of tuberculosis bacillus was delayed due to the late response of innate immunity (Pal, Moiz, Saif, & Zeeshan, 2016), thereby posing a higher risk of tuberculosis infection in diabetic people. In an another experiment with a mice, having streptozotocin-induced diabetes and exposed to tuberculosis infection, the diminished concentration of INF − γ and a lower concentration of inducible NO synthase by macrophage were reported. Again with a high concentration of glucose, INF − γ production was further reduced (Yamashiro et al., 2005). INF − γ gives immunity to host against tuberculosis infections. Therefore, diabetic host with reduced INF − γ concentration is more prone to tuberculosis infections. Also, diabetic people tend to have higher percentage of glycated haemoglobin which promotes rapid progression of tuberculosis with higher rate of mortality (Samadi, Safavi, & Mahmoodi, 2011). Besides chemotaxis, activation and antigen presentation by phagocytes play a crucial role as a defense mechanism against tuberculosis but these activities are also impaired in diabetic people (Dooley & Chaisson, 2009).
The number of diabetic people is increasing globally, specially in asian countries like China, India. Hence, the chance of developing tuberculosis infection is also getting increased. Here we are analysing the dynamics of tuberculosis infection in diabetic people using NSFD scheme.

Principles of NSFD scheme
Continuous dynamical systems are converted to discrete systems under NSFD scheme. Mickens has provided five principles to define these NSFD schemes. We have taken these rules as it is from Egbelowo, Harley, and Jacobs, (2017), the rules are as follows: (1) The orders of the discrete representation of the derivative must be equal to the orders of the corresponding derivatives appearing in the differential equations.
(2) Denominator functions for the discrete representations for derivatives must, in general, be expressed in terms of more complicated functions of the step sizes than those conventionally used. (3) Nonlinear terms must, in general, be modelled by nonlocal discrete representations. (4) All the special conditions that correspond to either the differential equation and/or its solutions should also correspond to the difference equation and/or its solutions. (5) The discrete scheme should not introduce extraneous or spurious solutions.
Our model is a three compartmental model and has been constructed using the first order ODE system. Then it has been transformed into the discrete form using NSFD principles and NSFD scheme is very much useful in inter compartmental dynamics. Lyapunov functions have been used to study the global stability of the systems. The paper has been organized in the following manner. Introduction has been articulated in Section 1. Section 2 describes the details of the methodology of the linear model and stability analysis. Section 3 deals with our nonlinear model and stability analysis. Section 4 has the conclusion of our study.

Methodology
In this study, we have used system of first order ordinary differential equations to convert biological expressions. In our study, we have pursued Michaelis and Menten equations of enzymatic reaction (Boutayeb, Chetouani, Achouyab, & Twizell, 2006;Boutayeb, Twizell, Achouayb, & Chetouani, 2004;Dhara, Baths, & Danumjaya, 2018). Then we have used NSFD schemes to study the local and global asymptotic stability analysis. The parameter values and data have been taken from the published articles (Boutayeb et al., 2006;Pealing et al., 2015). Below, we describe the procedure in the following manner. We assume The following parameters have been used in our model.
• I: Incidence of diabetes mellitus.
• μ 1 : Natural mortality rate due to diabetes mellitus without complications.
• μ 2 : Natural mortality rate due to complications in diabetes mellitus.
• μ 3 : Natural mortality rate in diabetes mellitus with complications due to tuberculosis infection. • λ 1 : Rate of developing a complication in diabetic individuals.
• λ 2 : Rate of developing a tuberculosis infection with complicated diabetic individuals.
• γ 1 : Rate at which complications in diabetes are cured.
• γ 2 : Rate at which complications in tuberculosis are cured.

Linear model
In this model, we analyse the dynamics of diabetes mellitus with tuberculosis infection. By observing the schematic diagram of the model, we derive the governing equations of the diabetes mellitus with tuberculosis infection disease as follows: ( 1 ) with the following initial data Below, we show the schematic diagram in Figure 1. Below, we discuss the different cases of λ i , i = 1, 2. Let λ 1 be a constant rate of developing complication in diabetic individual and λ 2 be the constant rate at which diabetic individual acquiring tuberculosis infection. In this case, the proposed model (1)-(3) is linear.
To obtain the equilibrium point for the system (1)- (3), we equate the right-hand side part of the system (1)-(3) equal to zero. Thus, we obtain the equilibrium point for the system (1)-(3) is as follows: where For the stability of the continuous model (1)-(3), we refer to Dhara et al. (2018). Below, we concentrate on the stability of the discrete model.

Local asymptotic stability
In this sub-section, we discuss the local asymptotic stability of the system (8)-(10).
The homogeneous system of (8)-(10) can be expressed as The characteristic equation corresponding to (14) is given by where Lemma 2.1: The roots of the cubic equation (15) satisfy | i | < 1, i = 1, 2, 3 if and only if the following conditions hold: The equilibrium solution of a difference system (14) is locally asymptotically stable if all the characteristic roots | i | < 1, i = 1, 2, 3.
The following parameter values listed in Table 1 are used in verifying the conditions of Lemma 2.1 We now compute To implement NSFD scheme (8)-(10) by using MATLAB software, we first divide the given time domain [t 0 , t n ] into finite number of equal subintervals with step size h where h = t i+1 − t i , i = 0, 1, 2, . . . , n − 1. Below, in Figure 2 we show the grid formulation: In Figures 3 and 4 we show the profiles of D(t), C(t) and T(t) for the time step size h = 0.25 using the NSFD scheme and fourth-order Runge-Kutta (RK-4) scheme.
We observe that the profiles of D(t), C(t) and T(t) in NSFD scheme are approaching the equilibrium point (16) as time increases but RK-4 scheme fails to converge. Thus, we conclude that the NSFD scheme has preserved the stability properties of of the system (1)-(3). Discussion: From the above Figure 3, we can observe that all the three profiles are converging to equilibrium solutions. D(t) profile converges within 50 days to equilibrium solution while C(t) and T(t) profiles have converged within 150 days and 250 days respectively. The converging rate for D(t) profile was faster than C(t) and T(t) profiles as expected. Treatment of tuberculosis infection in complicated diabetics is very difficult and thereby increasing the mortality rate in T(t) profiles. As a result, converging to the equilibrium solution was delayed in T(t) profile compared to D(t) and C(t) profiles.

Global stability
In this sub-section, we discuss the global stability of the system (8)-(10).
To discuss the global stability, we need to find a Lyapunov function. Below, we state the Lemma without proof.

Lemma 2.2:
The function V(x k ) is said to be a Lyapunov function for the discrete system (8)-(10), if the following conditions are hold (Nandakumaran, Datti, & George, 2017): The equilibrium solution x * of discrete system (8)-(10) is said to be globally asymptotically stable if there exists a Lyapunov function V that satisfy the conditions of Lemma 2.2 In this regard, we express (8)-(10) using equilibrium value as follows: We construct the following Lyapunov functional which satisfies all the conditions of Lemma 2.2 as Substituting (17)- (19) in (21), we obtain Using triangle inequality in (22), we arrive at After simplifying, we rewrite the equation (23) as Using the expressions of α i , i = 1, 2, 3 in (24), we obtain Since, h, μ 1 , μ 2 and μ 3 are positive, we express the equation (25) as We observe that (20) is the required Lyapunov function which satisfies all the requirements of Lemma 2.2 and hence we obtain the global stability for all values of the step size h.

Nonlinear model
In this section, we discuss the nonlinear model and corresponding local asymptotic stability and global stability. We now express the rate of developing a complication in diabetic individuals λ 1 as and the rate of developing a tuberculosis infection in complicated diabetic individuals λ 2 as Substituting (26)- (27) in (1)- (3), we obtain the following system with For the convenience, we let β 1 = (γ 1 + μ 2 ) and β 2 = (μ 3 + γ 2 ). The equilibrium points for the system (28)-(30) are (316.56, 648.32, 4928.92).
We note that while finding the equilibrium point (33), we obtain the following nonlinear equations To solve the nonlinear equations for finding the exact expression of D(t), C(t) and T(t) is impossible and we have computed the equilibrium point (33) numerically using MATHE-MATICA with the help of the parameter values listed in Table 1. We consider the following NSFDS for the system (28)-(30) as follows: We rewrite the above system (34)-(36) as The equilibrium solutions of (37)-(39) are given by Below, we prove that by solving (40)-(42), we obtain the equilibrium solutions those are equivalent to (31)-(33). We consider the equations (40)-(42) and after simplifying them we arrive at Now we observe that the system of equations (43)-(45) are exactly same as the stationary equations of (28)-(30) and hence the system (43)-(45) has same equilibrium solutions those are equivalent to (31)-(33).

Proof:
In case of the equilibrium point (31), the matrix in (46) implies that ⎛ We observe that the matrix (47) is an upper triangular matrix and hence the eigenvalues are diagonal entries only. Thus, we have We note that | i | < 1, i = 1, 3 and while | 2 | < 1 provided β 1 > 1. Therefore, the condition for the local asymptotic stability of the equilibrium solution (31) is Theorem 3.2: The equilibrium point (32) is locally asymptotically stable for the system (37)-(39) if (β 1 + β 2 ) > 1.

Proof:
In case of the equilibrium point (32), the matrix in (46) implies that The eigenvalues of the matrix (48) are The above equation can be rewritten as We note that the roots of the quadratic equation (49) satisfy | i | < 1, i = 2, 3 if and only if the following conditions hold: Using the equilibrium point (32) and the parameters values listed in Table 1 with h = 0.25, we compute (i) 1 + P 1 + Q 1 = 0.0010401 > 0, (ii) 1 − P 1 + Q 1 = 3.5482 > 0, and Thus, all the conditions (i)-(iii) of (49) are satisfied and hence | i | < 1, i = 2, 3.

Proof:
We now study the local asymptotic stability of the equilibrium point (33). Substituting (33) and using the parameters values listed in Table 1 To show the local asymptotic stability of the equilibrium point (33) computationally, we choose the initial data: D 0 = 274, C 0 = 137 and T 0 = 74.
To implement the NSFD scheme (37)-(39), we again use the data mentioned in Table 1 and MATLAB software.
Below, in Figure 5, we show the profile of D(t), C(t) and T(t) for the time step size h = 0.25.
We observe that the profile of D(t), C(t) and T(t) are approaching the equilibrium point (33) as time increases. Thus, we conclude that the NSFD scheme has preserved the stability properties of of the system (28)-(30).
Discussion: From the above Figure 5, we can observe that all the three profiles (D(t), C(t), T(t)) are converging to the equilibrium solutions indicating the stability of our nonlinear discrete scheme. The D(t) profile has reached equilibrium solution earlier than C(t) and T(t). T(t) profile takes longer time to reach equilibrium solution because of tuberculosis infections. From the plot, we can observe that D(t) profile has taken less than 50 days only while C(t) profile has taken less than 100 days to reach equilibrium solutions. But T(t) profile has taken more than 300 days. This indicates the impact of tuberculosis infection to reach equilibrium in complicated diabetic (C(t)) people. Besides, T(t) has covered much high number of populations compared to D(t), C(t) profiles. From this plot, we may say that the tuberculosis infection spreads faster in the complicated diabetic people infecting a huge number of populations. As T(t) profile reach equilibrium solution late, so T(t) profile has increased chance of infecting a high number of diabetic people. Henceforth the association of diabetic people and tuberculosis infection turns out to be very critical.

Conclusion
In this article, we have studied the association of diabetes mellitus developing complication and its association with tuberculosis infection. We have studied both linear and nonlinear models using first oder ODE system. We have performed local and global stability for both models using the NSFD scheme. Lyapunov function has been used to derive global stability. Then we have simulated our model using MATLAB software. From the plots (Figures 3  and 5), we have observed that the nonlinear model has performed well compared to our linear one. The threefold higher risk of tuberculosis infection in diabetes mellitus has been reflected with higher accuracy in the nonlinear model. Both models have reached equilibrium solutions as the time increases indicating that stabilities have been preserved for our models. From the plots, we may infer that tuberculosis infections have an adverse effect on diabetic people and it spreads faster among them. The T(t) profile attains stability later and also covers a large number of populations than D(t), C(t) profiles. Therefore tuberculosis infection imparts a double burden on patients with diabetes. Hence necessary actions are urgently required to combat this joint epidemics.