Numerical study for multi-strain tuberculosis (TB) model of variable-order fractional derivatives.

In this paper, we presented a novel multi-strain TB model of variable-order fractional derivatives, which incorporates three strains: drug-sensitive, emerging multi-drug resistant (MDR) and extensively drug-resistant (XDR), as an extension for multi-strain TB model of nonlinear ordinary differential equations which developed in 2014 by Arino and Soliman [1]. Numerical simulations for this variable-order fractional model are the main aim of this work, where the variable-order fractional derivative is defined in the sense of Grünwald-Letnikov definition. Two numerical methods are presented for this model, the standard finite difference method (SFDM) and nonstandard finite difference method (NSFDM). Numerical comparison between SFDM and NSFDM is presented. It is concluded that, NSFDM preserves the positivity of the solutions and numerically stable in large regions than SFDM.

Introduction

Variable-order fractional calculus (i.e., the fractional differentiation and integration of variable order) is the generalization of classical calculus and fractional calculus, which were invented by Newton and Leibnitz hundreds of years ago. Now the study on it becomes a hot pot in recent ten years [2], [3], [4], [5], [6], [7]. It has turned out that many problems in physics, biology, engineering, and finance can be described excellently by models using mathematical tools from variable-order fractional calculus, such as mechanical applications [2], diffusion process [5], multifractional Gaussian noise [8], and FIR filters [9]. For more details, see [7], [10] and references therein. Understanding the transmission characteristics of infectious diseases in communities, regions and countries can lead to better approaches to decrease the transmission of these diseases [11]. Variable-order fractional derivative is good at depicting the memory property which changes with time or spatial location [3], [5].

TB is growing more resistant to treatment worldwide according to study released in August 2012 in the journal Lancet, a finding that suggests the potentially fatal disease is becoming more difficult and costly to treat [12]. In this article we focused our attention in Egypt.

We consider in this work a model developed by Arino and Soliman for TB [1]. The model incorporates three strains, drug-sensitive, MDR and XDR. Several papers considered modeling TB such as [13], [14], but the model we consider here includes several factors of spreading TB such as the fast infection, the exogenous reinfection and secondary infection along with the resistance factor. The main aim of this paper was to study numerically the multi-strain TB model of variable-order fractional derivatives which incorporates three strains: drug-sensitive, MDR and XDR. We develop a special class of numerical method, known as NSFDM for solving this model. This technique, developed by Mickens (1980) [15], [16], [17], [18], [19], [20], [21], [22], [23] has brought a creation of new numerical schemes preserving the physical properties, especially the stability properties of equilibria, of the approximated system. Numerical comparison between NSFDM and SFDM is presented. When the secondary infection generated by an infected individual exceeds the unity, there are no analytical results proved for the model, such as the existence and stability of the endemic equilibrium . In this case we use the developed NSFD numerical scheme to approximate the endemic solution numerically and investigate its stability. Furthermore, with the help of the NSFDM, we answer the following question: Given the data provided by the World Health Organization (2012) on the current parameters corresponding to the propagation of the TB in Egypt, what would be the required rate of treatment to achieve in order to control the disease? The proposed method showed its superiority in preserving the positivity (compared to the numerical standard method considered in this work) of the state variables of the systems under study. This is an essential requirement when simulating systems especially those arising in biology. This paper is organized as follows: In Section ‘Mathematical model’, Mathematical model is presented. Preliminaries and notations on variable-order fractional differential equations are given, in Section ‘Preliminaries and notations’. Equilibrium points and their asymptotic stability are presented in Section ‘Variable-order fractional derivatives for multi-strain TB model’. Variable-order fractional derivatives for the multi-strain TB model are presented; moreover, the construction of the proposed nonstandard numerical scheme is carried out in Section ‘Equilibrium points and their asymptotic stability’. In Section ‘Numerical results and simulations’, Numerical results and simulation are discussed. Finally, in Section ‘Conclusions’ we presented the conclusions.

Mathematical model

The multistrain TB-model given in [1] can be formulated as follows:All variables in above system and their definition are in Table 1. Also, all parameters and their interpretation are in Table 2.

Table 1

All variables of the system (1), (2), (3), (4), (5), (6), (7), (8) and their interpretation.

Variable	Definition
S(t)	The susceptible population individuals who have never encountered TB
L_s(t)	The individuals infected with the drug-sensitive TB strain but who are in a latent stage, i.e., who are neither showing symptoms nor infecting others
L_m(t)	Individuals latently infected with MDR-TB
L_x(t)	Individuals latently infected with XDR-TB
I_s(t)	Individuals infected with the drug-sensitive TB strain who are infectious to others (and most likely, showing symptoms as well)
I_m(t)	Those individuals who are infectious with the MDR-TB strain
I_x(t)	Individuals who infectious with the XDR-TB strain
R(t)	Those individuals for whom treatment was successful
N(t)	The total population
	N=S+L_s+L_m+L_x+I_s+I_m+I_x+R

Table 2

All parameters of the system (1), (2), (3), (4), (5), (6), (7), (8) and their interpretation.

Parameter	Interpretation
b	Birth/recruitment rate
d	Per capita natural death rate
	Disease dynamics
β_r	Transmission coefficient for strain r
λ_r	Proportion of newly infected individuals developing LTBI with strain r
1-λ_r	Proportion of newly infected individuals progressing to active TB with strain r due to fast infection
ε_r	Per capita rate of endogenous reactivation of L_r
α_r1,α_r2	Proportion of exogenous reinfection of L_r1 due to contact with I_r2
γ_r	Per capita rate of natural recovery to the latent stage L_r
δ_r	Per capita rate of death due to TB of strain r
	Treatment related
t_1s	Per capita rate of treatment for L_s
t_2r	Per capita rate of treatment for I_r. Note that t_2x is the rate of successful treatment of I_x,r∈{x,m,s}
1-σ_r	Efficiency of treatment in preventing infection with strain r
P_1	Probability of treatment success for L_s
1-P_1	Proportion of treated L_s moved to L_m due to incomplete treatment or lack of strict compliance in the use of drugs
P_2	Probability of treatment success for I_s
1-P_2	Proportion of treated I_s moved to L_m due to incomplete treatment or lack of strict compliance in the use of drugs
P_3	Probability of treatment success for I_m
1-P_3	Proportion of treated I_m moved to L_x due to incomplete treatment or lack of strict compliance in the use of drugs

The basic reproduction number

The basic reproduction number for system (1), (2), (3), (4), (5), (6), (7), (8) is given by [1]where

Theorem [1] assumes thatThen the disease free equilibrium is globally asymptotically stable when and endemic equilibria are locally asymptotically stable when .

Preliminaries and notations

In this section, some basic definitions and properties in the theory of the variable-order fractional calculus are presented.

Grünwald–Letnikov approximation

We will begin with the signal variable-order fractional differentialwhere , and denotes the variable fractional order derivative, defined bywhere , and is the nth-order Riemann–Liouville integral operator defined aswhere is the gamma function.

To apply Miken’s scheme, we have chosen this Grünwald–Letnikov approximation variable-order fractional derivative as follows [15]:where denotes the integer part of t and h is the step size; therefore, Eq. (16) is discretized aswhere , and , are the Grünwald–Letnikov coefficients defined as.

Variable-order fractional derivatives for multi-strain TB model

In the following, we introduce the multi-strain TB model of variable-order fractional derivatives which is the integer order given in system (1), (2), (3), (4), (5), (6), (7), (8), and the new system is described by variable-order fractional differential equations as follows:where is the Caputo variable fractional order derivative. Because model (18), (19), (20), (21), (22), (23), (24), (25) monitors the dynamics of human populations, all the parameters are assumed to be nonnegative.

Equilibrium points and their asymptotic stability

Let and consider the system (18), (19), (20), (21), (22), (23), (24), (25)With the initial values .

To evaluate the equilibrium points letFrom which we can get the equilibrium points .

To evaluate the asymptotic stability letSo the equilibrium point is locally asymptotically stable if all eigenvalues of Jacobian evaluated at the equilibrium point satisfy [16]To evaluate the equilibrium points, letNow, if , and .

Then the disease free equilibrium (DFE) is .

We calculate the Jacobian matrix of the system (18), (19), (20), (21), (22), (23), (24), (25) at the disease free equilibrium point as follows:where , , . The characteristic equation associated with above matrix is . Then the eigenvalues of Jacobian matrix are , , by using Theorem (Routh Hurwitz criteria) [17], these roots are negative or have negative real parts and DFE is locally asymptotically stable if all eigenvalues of the Jacobian matrix satisfies , , . For simplicity, we will determine the stability of the DFE numerically by using Table 3 and put . Then eigenvalues are , , , , , , , . So, if , the DFE is locally asymptotically stable since , , .

Table 3

All parameters in the system (18), (19), (20), (21), (22), (23), (24), (25) and the reference of the parameters.

Parameter	Value	Reference
b	3190	Assumed
d	0.38	[26]
β_s=β_m=β_x	14	[26]
λ_s=λ_m=λ_x	0.5	Assumed
ε_s=ε_m=ε_x	0.5	Assumed
α_r1,r2	0.05	Assumed
γ_s=γ_m=γ_x	0.3	Assumed
t_1s	0.88	[26]
t_2r:r∈(s,m,x)	t_2s=0.88; t_2m=t_2x=0.034	[26]
σ_r	0.25	[26]
P_r	0.88	[26]
δ_r	0.045	[26]

If at least one of the infected variables is non-zero, then the solutions for model (18), (19), (20), (21), (22), (23), (24), (25) are the endemic equilibrium [1]. This system is highly nonlinear in , and , and hence explicit solutions are not obtainable. So we solved the system (18), (19), (20), (21), (22), (23), (24), (25) numerically to obtain endemic fixed point using NSFDM.

SFD discretization

SFD methods are simple numerical methods for approximating the solutions of differential equations using finite differences to approximate the derivatives.

The forward Euler method is one of these methods, in this method the derivative term is replaced by , where h is the step size, for more details see [18].

NSFD discretization

The nonstandard finite difference schemes were introduced by Mickens in the 1980s as a powerful numerical method that preserves significant properties of exact solutions of the involved differential equation [19]. The concept of the nonstandard finite difference method is discussed in [20].

A numerical scheme is called NSFD discretization if at least one of the following conditions is satisfied [18]:

Nonlocal approximation is used.

The discretization of derivative is not traditional and uses a nonnegative function [19], [20].

To describe the main aspects of NSFD schemes, we consider an ODE in the formwhere λ is a possibly vector, parameter. Given a mesh-grid that just for simplicity we assume to be equispaced with step-size 0, NSFD schemes are constructed by the following two main steps: 1- the derivative at the left-hand side of (27) is replaced by a discrete representation in the formwhere is an approximation of , 2-the nonlinear term in (27) is replaced by a nonlocal discrete representation depending on some of the previous approximations.

For example, if there are nonlinear terms such as in the differential equation, these are replaced by or .

Let us denote by , , , , , , and the values of the approximations of , , , , , , and respectively, for and h is the timestep of the scheme. The sequences , , , , , , and should be nonnegative in order to be consistent with the biological nature of the model [21].

NSFDM has many advantages than SFDM, for more details see [20], [21], [22], [23], [24]. Generally speaking, we can say that NSFDM is more efficient and accurate than SFDM [15], [25].

NSFD for variable-order fractional derivatives system

The system (18), (19), (20), (21), (22), (23), (24), (25) can be discretized as follows:where the discretization for is given asAnd , where the nonlocal approximations are used for the nonlinear terms and the following denominator functions are used:We obtain,

Numerical results and simulations

Since most of the variable-order fractional differential equations do not have exact analytic solutions, so approximation and numerical techniques must be used. Several analytical and numerical methods have been proposed to solve variable-order fractional differential equations. For numerical solutions of the system (18), (19), (20), (21), (22), (23), (24), (25) one can use NSFDM, the approximate solution , , , , , , , is displayed in Fig. 1, when and in Fig. 2, when , in each figures, and three different values of , , are considered. The approximate solutions are displayed in Fig. 2 that, the equilibrium point of NSFDM is locally asymptotically stable when , , where the eigenvalues are given as , , , , , , . By applying the relationship (26) we obtained that, , . When , system (18), (19), (20), (21), (22), (23), (24), (25) is the classical integer-order system. Moreover, we observed that, the integer order derivative can be used to characterize the short memory of systems, and the variable-order fractional derivative can be employed to depict the variable memory of systems. In Fig. 3, we presented the result obtained by NSFDM and SFDM and ode45 schemes with step size and , and we observed that, all numerical methods converge almost to the equilibrium point when . In Table 4, we reported the convergence behavior of numerical methods to the disease free equilibrium, and in Table 5, we reported the convergence behavior of numerical methods to the equilibrium point .

Fig. 1

Profiles obtained by using NSFDM for solving variable-order fraction model with different , , an .

Fig. 2

Profiles obtained by using NSFDM for solving variable-order fraction model with different , , and .

Fig. 3

Profiles obtained by using different methods with , , , and .

Table 4

Result obtained by SFDM and NSFDM for , , and initial conditions as with different time step size.

h	SFDM	NSFDM
0.01	Convergent	Convergent
0.1	Convergent	Convergent
1	Convergent	Convergent
20	Divergent	Convergent
100	Divergent	Convergent

Table 5

Result obtained by SFDM and NSFDM for , , and initial conditions as with different time step size.

h	SFDM	NSFDM
0.01	Convergent	Convergent
0.1	Convergent	Convergent
1	Divergent	Convergent
20	Divergent	Convergent
100	Divergent	Convergent

From Table 4, we can conclude that NSFDM unconditionally converges to the correct disease free equilibria for large h, while the SFDM converges only when h is small.

From Table 5, we can conclude that NSFD scheme unconditionally converges to the equilibrium point for large h, while the SFD scheme converges only when h is small. Moreover, the system (28), (29), (30), (31), (32), (33), (34), (35) is unconditionally locally asymptotically stable.

Previous Fig. 4(a)–(d), illustrates propagation of TB along the time when as follows:

Fig. 4

Illustrate propagation of multi-strain TB along the time , , , and , by using NSFDM.

In Fig. 4(a), the relationship between and illustrates that, there are individuals succeeded treatment with them and may exposed to infection again by contagious members of the first strain. At the beginning of the period of the time the number of members increases and the number of members decreases, then after time steps the curves intersect again, will be responsible to treatment and their numbers will be decreased.

In Fig. 4(b), the relationship between and , describes the spread of infection from the members of the third strain to healthy people, then the number of infectious people increases and the number of healthy people decreases with proper time.

In Fig. 4(c), the relationship between and , describes the spread of contagious from the members of the second strain to healthy people, then the number of infectious people increases and the number of healthy people decreases with proper time.

In Fig. 4(d), the relationship between and , describes the spread of contagious from the members of the first strain to individuals who carry the disease latent of the first strain , after time steps the curves intersect again then will be responsible to treatment and the number of them decreases.

In Fig. 5, we presented the result obtained by NSFDM and SFDM schemes with step size and , . We can clearly see, all schemes converge to correct equilibrium point when .

Fig. 5

Profiles obtained by using NSFDM and SFDM with , , , and .

In Fig. 6, we presented the results obtained by NSFD and SFD schemes with step size and . As we can clearly see, the SFD scheme is unstable and the solutions are divergent, so we cannot use this scheme to solve the system when step size is large.

Fig. 6

Profiles obtained by using NSFDM and SFDM with , , , and .

From these numerical results obtained in this work we can control the disease and turn the endemic point to the disease free point as follows:

Let us consider:Then,So, we derive the rate of treatment required for achieving control of the disease.

For example, if we choose the following elements which belong to such as , , and , we obtained the disease free point (see Fig. 7).

Fig. 7

Profiles obtained by using NSFDM for , , , and , .

Conclusions

In this article, a novel multi-strain TB model of variable-order fractional derivatives which incorporates three strains: drug-sensitive, MDR and XDR, is studied. It can be concluded from the numerical results presented in this paper, that the variable-order fractional TB model given here is a general model than the integer and fractional order models. Furthermore, the integer order model can be used to characterize the short memory of systems, and the variable-order fractional model can be employed to depict the variable memory of systems. Moreover, we can conclude that NSFDM is more efficient for solving variable-order fractional mathematical model for multi-strain TB, than the SFDM, because it preserves the positivity of the solution and the stability regions using it are bigger than the SFDM stability regions. All results in this paper are obtained using MATLAB (R2013a), on a computer machine with intel (R) core i3-3110M @ 2.40 GHz and 4 GB RAM.

Conflict of Interest

The authors have declared no conflict of interest.

Compliance with Ethics Requirements

This article does not contain any studies with human or animal subjects.