# Numerical Study for Time Delay Multistrain Tuberculosis Model of Fractional Order

Authors: Nasser Hassan Sweilam; Seham Mahyoub Al-Mekhlafi; Taghreed Abdul Rahman Assiri
Journal: Complexity
 (2017)

---

## Abstract

A novel mathematical fractional model of multistrain tuberculosis with time delay memory is presented. The proposed model is governed by a system of fractional delay differential equations, where the fractional derivative is defined in the sense of the Grünwald–Letinkov definition. Modified parameters are introduced to account for the fractional order. The stability of the equilibrium points is investigated for any time delay. Nonstandard finite deference method is proposed to solve the resulting system of fractional-order delay differential equations. Numerical simulations show that nonstandard finite difference method can be applied to solve such fractional delay differential equations simply and effectively.


---

## Body

## 1. Introduction

It is known that tuberculosis (TB) is one of the most important infectious diseases and is considered as the second largest cause of mortality by infectious diseases and a challenging disease to control [1]. Time delays required to treatment of active TB present a major obstacle to the control of a TB epidemic [2]; it worsens the disease, increases the risk of death, and enhances tuberculosis transmission to the community [3, 4]. Both patient and the health system may be responsible for the treatment delay [3]. On the other hand, mathematical models are quite important and efficient tool to describe and investigate TB diseases; see [5–9]. In [10], Silva et al. presented TB model with time-delay memory. Herein, we consider a general model of multistrain TB diseases with time-delay memory. A discrete time delay is incorporated, in the variables of active TB infection of two and three strains, to represent the required time to commencement of treatment and diagnosis [11].The multistrain TB model incorporates three strains:extensively drug-resistant (XDR), emergingmultidrug resistant (MDR), and drug sensitive, and has been developed by Arino and Soliman [12] in 2015. Several factors of spreading TB such as the exogenous reinfection, the fast infection, and secondary infection are included in this model. Sweilam et al. introduced some numerical studies for this model in [13–16].Fractional differential equations have been the focus of many studies due to their frequent appearance in various sciences [13–20]. The general theory of differential equations with delays (DDEs) is widely developed and discussed in the literature [21–25].Delayed fractional differential equations (DFDEs) are also used to describe dynamical systems [26–28]. Recently, DFDEs begin to raise the attention of many researchers [29–33]. Relatable and efficient numerical techniques for DFDEs are very necessary and important [34].Nonstandard finite difference method (NSFDM) was firstly proposed by Mickens [35] in 1980s to solve numerically theordinary differential equations (ODEs) andpartial differential equations (PDEs) with more accuracy thanstandard finite difference method (SFDM). It is considered as a powerful numerical scheme that preserves properties of exact solutions of the differential equation [36].The main aim of work is to study numerically the solutions of fractional-order model of multistrain TB with time delay memory. The presence of fractional-order and time delays in the model can lead to a notable increase in the complexity of the observed behavior, and the solution continuously depends on all the previous states. An efficient numerical method, NSFDM, is used to numerically solve the fractional-order delay model. The rest of the paper is organized as follows: In Section2, we present a fractional order model with time delay for multistrain TB. Stability of equilibrium points is presented in Section 3. NSFD for fractional-order delay differential equations is introduced in Section 4. Some numerical simulations are given in Section 5, and conclusion in Section 6. Some definitions on fractional calculus and some properties of nonstandard discretization are given in Appendix.


## 2. Fractional Multistrain TB Model with Time Delay

In this section, a multistrain TB model of fractional-order and time delay memory is presented. The population of interest is divided into eight compartments depending on their epidemiological stages as follows: susceptible(S); latently infected with drug sensitive TB (Ls); latently infected with MDR TB (Lm); latently infected with XDR TB (Lx); sensitive drug TB infectious (Is); MDR TB infectious (Im); XDR TB infectious (Ix); recovered R. One biological meaning of the given parameters is given in Table 1. One of the main assumptions of this model is that the total population N(t), with N(t)=S(t)+Ls(t)+Lm(t)+Lx(t)+Is(t)+Im(t)+Ix(t)+R(t), is variable of the time. We introduce a discrete time delay in the state variables Im and Ix, denoted by τ, that represents the time required for diagnosis and commencement of treatment of active TB infection of two and three strains. The parameters in the modified the model are described in Table 2; see [37]. The modified system of multistrain TB model of fractional-order and time delay is(1)DtαS=bα-dαS-βsαSIsN-βmαSImN-βxαSIxN,DtαLs=λsαβsαSIsN+σsαλsαβsαRIsN+γsαIs-αssαβsαLsIsN-αsmαβmαLsImN-αsxαβxαLsIxN-dα+εsα+t1sαLs,DtαLm=λmαβmαSImN+σmαλmαβmαRImN+γmαIm+αsmαβmαλmαLsImN+1-P1αt1sαLs+1-P2αt2sαIs-αmmαβmαLmImN-αmxαβxαLmIxN-dα+εmαLm,DtαLx=λxαβxαSIxN+σxαλxαβxαRIxN+γxαIx+αsxαβxαλxαLsIxN+αmxαβxαλxαLmIxN+1-P3αt2mαIm-αxxαβxαLxIxN-dα+εxαLx,DtαIs=αssαβsαLsIsN+1-λsαβsαSIsN+σsαRIsN+εsαLs-dα+δsα+t2sα+γsαIs,DtαIm=αmmαβmαLmImN+1-λmαβmαSImN+σmαRImN+αsmαLsImN+εmαLm-dα+δmα+γmαIm-t2mαImt-τ,DtαIx=αxxαβxαLxIxN+1-λxαβxαSIxN+σxαRIxN+αsxαLsIxN+αmxαLmIxN+εxαLx-dα+δxα+γxαIx-t2xαIxt-τ,DtαR=P1αt1sαLs+P2αt2sαIs+P3αt2mαIm+t2xαIxt-τ-σsαβsαRIsN-σmαβmαRImN-σxαβxαRIxN-dαR.The initial conditions for system (1) are S(ξ)=θ1(ξ),  Ls(ξ)=θ2(ξ),Lm(ξ)=θ3(ξ),  Lx(ξ)=θ4(ξ),  Is(ξ)=θ5(ξ),  Im(ξ)=θ6(ξ),Ix(ξ)=θ7(ξ),  R(ξ)=θ8(ξ),ξ∈[-τ,0], where θ=(θ1,θ2,…,θ8)T∈C, where C is the Banach space C([0,τ],R8). From biological meaning, we further assume that θi>0 for i=1,…,8. Throughout this paper, we focus on the dynamics of the solutions of (1) in the restricted region, Ω={(S,Ls,Lm,Lx,Is,Im,Ix,R)∈R8∣S+Ls+Lm+Lx+Is+Im+Ix+R≤bα/dα}. We refer here to [24, 31], for local existence, uniqueness, and continuation results.Table 1
Interpretation of the variable states of system (1).
Variable
Interpretation
S
(
t
)
Individuals have never encountered TB.
L
s
(
t
)
The individuals infected with drug-sensitive TB but not infectious.
L
m
(
t
)
Infected with MDR-TB but not infectious.
L
x
(
t
)
Infected with XDR-TB but not infectious.
I
s
(
t
)
Able to infect others with drug sensitive strain.
I
m
(
t
)
Able to infect others with MDR strain.
I
x
(
t
)
Able to infect others with XDR strain.
R
(
t
)
Recovered by getting a successful treatment.
N
(
t
)
The variable of population size.
N
(
t
)
=
S
(
t
)
+
L
s
(
t
)
+
L
m
(
t
)
+
L
x
(
t
)
+
I
s
(
t
)
+
I
m
(
t
)
+
I
x
(
t
)
+
R
(
t
)
.Table 2
All adapted parameters and their interpretation of system (1).
Parameter
Interpretation
d
α
Natural death rate
b
α
Birth rate
λ
j
α
Rate of infected individuals move toLj with strain j∈{s,m,x}
1
-
λ
j
α
Rate of newly infected individuals progressing to active TB with strainj∈{s,m,x}
β
j
α
Transmission coefficient with strainj∈{s,m,x}
ε
j
α
Rate of endogenous reactivation ofLj
γ
j
α
Rate of natural recovery to the latent stageLj
δ
j
α
Rate of death due to TB of strainj
α
j
1
α,  αj2α
Rate of exogenous reinfection ofLj1 due to contact with Ij2
1
-
σ
j
α
Efficiency of treatment in preventing infection with strainj
P
1
α
Probability of treatment success forLs
1
-
P
1
α
Proportion of treatedLs moved to Lm due to incomplete treatment or lack of strict compliance in the use of drugs
P
2
α
Probability of treatment success forIs
1
-
P
2
α
Proportion of treatedIs moved to Lm due to incomplete treatment or lack of strict compliance in the use of drugs
P
3
α
Probability of treatment success forIm
1
-
P
3
α
proportion of treatedIm moved to Lx due to incomplete treatment or lack of strict compliance in the use of drugs
t
1
s
α
Rate of treatment forLs
t
2
j
α
Rate of treatment forIj. Note that t2x is the rate of successful treatment of Ix,  j∈{x,m,s}The unique solution(St,Lst,Lmt,Lxt,Ist,Imt,Ixt,Rt) of (1) with initial condition exists for all time t≥0. Consider the solutions of (1), for (S,Ls,Lm,Lx,Is,Im,Ix,R)∈Ω′, where Ω′ is the interior of Ω, for all ξ∈[-τ,0]. Then the solutions stay in the interior of the region for all time t≥0; that is, the region is positively invariant with respect to system (1) (see, e.g., [31]). Model (1) has a disease-free equilibrium given by E0=(bα/dα,0,0,0,0,0,0,0); see [32].


### 2.1. Basic Reproduction Number

The basic reproduction number,R0, is defined as the expected number of secondary cases produced, in a completely susceptible population, by a typical infective individual [32]. Herein, we apply the method in [32] to drive R0. The order of the infected variables is (2)I≔Ls,Lm,Lx,Is,Im,Ix′.The vector representing new infections into the infected classes F is given by (3)F≔λsαβsαSIsN+σsαλsαβsαRIsNλmαβmαSImN+σmαλmαβmαRImNλxαβxαSIxN+σxαλxαβxαRIxN1-λsαβsαSIsN+σsαRIsN1-λmαβmαSImN+σmαRImN1-λxαβxαSIxN+σxαRIxN.The vector V representing other flows within and out of the infected classes I is given by (4)V≔αssαβsαLsIsN-αsmαβmαLsImN-αsxαβxαLsIxN+γsαIs-dα+εsα+t1sαLs,+γmαIm+αsmαβmαλmαLsImN+1-P1αt1sαLs+1-P2αt2sαIs,-αmmαβmαLmImN-αmxαβxαLmIxN-dα+εmαLm,+γxαIx+αsxαβxαλxαLsIxN+αmxαβxαλxαLmIxN+1-P3αt2mαIm-αxxαβxαLxIxN-dα+εxαLx,αssαβsαLsIsN+εsαLs-dα+δsα+t2sα+γsαIs,αmmαβmαLmImN+αsmαLsImN+εmαLm-dα+δmα+γmα+t2mαIm,αxxαβxαLxIxN+1-λxαβxα+αsxαLsIxN+αmxαLmIxN+εxαLx-dα+δxα+γxα+t2xαIx,.The matrix of new infections F and the matrix of transfers between compartments V are the Jacobian matrices obtained by differentiating F and V with respect to the infected variables I and evaluating at the disease-free equilibrium. They take the form (5)F≔0A0B,V≔CDEF2, where(6)A=λsαβsα000λmαβmα000λsαβsα,B=1-λsαβsα0001-λmαβmα0001-λsαβsα,C=dα+εsα+t1sα00-1+P1αt1sαdα+εmα000dα+εxα,D=-γsα00-1+P2αt2sα-γmα00-1+P3αt2mα-γxα,F2=dα+δsα+γsα+t2sα000dα+δmα+γmα+t2mα000dα+δxα+γxα+t2xα,E=-εsα000-εmα000-εxα.Then the basic reproduction number R0 for system (1) is the spectral radius of the next generation matrix and is given by(7)R0=ρFV-1=maxR0s,R0m,R0x,where(8)R0s=βsαεsα+1-λsαdα+t1sαεsα+dα+t1sαt2sα+δsα+dα+γsαt1sα+dα,R0m=βmαεmα+1-λmαdαεmα+dαt2mα+δmα+dα+dαγmα,R0x=βxαεxα+1-λxαdαεxα+dαt2xα+δxα+dα+dαγxα.


## 2.1. Basic Reproduction Number

The basic reproduction number,R0, is defined as the expected number of secondary cases produced, in a completely susceptible population, by a typical infective individual [32]. Herein, we apply the method in [32] to drive R0. The order of the infected variables is (2)I≔Ls,Lm,Lx,Is,Im,Ix′.The vector representing new infections into the infected classes F is given by (3)F≔λsαβsαSIsN+σsαλsαβsαRIsNλmαβmαSImN+σmαλmαβmαRImNλxαβxαSIxN+σxαλxαβxαRIxN1-λsαβsαSIsN+σsαRIsN1-λmαβmαSImN+σmαRImN1-λxαβxαSIxN+σxαRIxN.The vector V representing other flows within and out of the infected classes I is given by (4)V≔αssαβsαLsIsN-αsmαβmαLsImN-αsxαβxαLsIxN+γsαIs-dα+εsα+t1sαLs,+γmαIm+αsmαβmαλmαLsImN+1-P1αt1sαLs+1-P2αt2sαIs,-αmmαβmαLmImN-αmxαβxαLmIxN-dα+εmαLm,+γxαIx+αsxαβxαλxαLsIxN+αmxαβxαλxαLmIxN+1-P3αt2mαIm-αxxαβxαLxIxN-dα+εxαLx,αssαβsαLsIsN+εsαLs-dα+δsα+t2sα+γsαIs,αmmαβmαLmImN+αsmαLsImN+εmαLm-dα+δmα+γmα+t2mαIm,αxxαβxαLxIxN+1-λxαβxα+αsxαLsIxN+αmxαLmIxN+εxαLx-dα+δxα+γxα+t2xαIx,.The matrix of new infections F and the matrix of transfers between compartments V are the Jacobian matrices obtained by differentiating F and V with respect to the infected variables I and evaluating at the disease-free equilibrium. They take the form (5)F≔0A0B,V≔CDEF2, where(6)A=λsαβsα000λmαβmα000λsαβsα,B=1-λsαβsα0001-λmαβmα0001-λsαβsα,C=dα+εsα+t1sα00-1+P1αt1sαdα+εmα000dα+εxα,D=-γsα00-1+P2αt2sα-γmα00-1+P3αt2mα-γxα,F2=dα+δsα+γsα+t2sα000dα+δmα+γmα+t2mα000dα+δxα+γxα+t2xα,E=-εsα000-εmα000-εxα.Then the basic reproduction number R0 for system (1) is the spectral radius of the next generation matrix and is given by(7)R0=ρFV-1=maxR0s,R0m,R0x,where(8)R0s=βsαεsα+1-λsαdα+t1sαεsα+dα+t1sαt2sα+δsα+dα+γsαt1sα+dα,R0m=βmαεmα+1-λmαdαεmα+dαt2mα+δmα+dα+dαγmα,R0x=βxαεxα+1-λxαdαεxα+dαt2xα+δxα+dα+dαγxα.


## 3. Equilibrium Points and Their Asymptotic Stability

To discuss the local asymptotic stability for evaluating the equilibrium points, let us consider the following [38]:(9)DtαS=DtαLs=DtαLm=DtαLx=DtαIs=DtαIm=DtαIx=DtαR=0.Then, from (A.1)(10)giS¯,Ls¯,Lm¯,Lx¯,Is¯,Im¯,Ix¯,R¯=0,i=1,2,3,…,8,where (S¯,Ls¯,Lm¯,Lx¯,Is¯,Im¯,Ix¯,R¯) denotes any equilibrium point.


### 3.1. Stability of the Disease-Free Equilibrium

IfIs(t)=Im(t)=Ix(t)=0 then Ls(t)=Lm(t)=Lx(t)=0, R(t)=0, and S(t)=bα/dα. Then the disease-free equilibrium (DFE) is E0={(bα/dα,0,0,0,0,0,0,0)}.Let us consider the coordinate transformation:s(t)=S(t)-S¯,  ls(t)=Ls(t)-Ls¯,  lm(t)=Lm(t)-Lm¯,  lx(t)=Lx(t)-Lx¯,  is(t)=Is(t)-Is¯,  im(t)=Im(t)-Im¯,  ix(t)=Ix(t)-Ix¯,  r(t)=R(t)-R¯. The corresponding characteristic equation for DFE is given as follows:(11)JE0-λI=λ-a000bcd100λ-e00f0000gλ-h0pq00000λ-r0st00u00λ-v00000w00λ-x+t2mαe-λτ00000y00λ-z+t2xαe-λτ00m00njt2xαe-λτλ-a,where a=-dα,  b=-βsα,  c=-βmα,  d1=-βxα,  e=-(dα+εsα+t1sα),  f=γsα+λsαβsα,  g=(1-p1α)t1sα,  h=-(dα+εmα),  p=(1-p2α)t2sα,  q=γmα+λmαβmα,  r=-(dα+εxα),  s=(1-p3α)t2mα,  t=γxα+λxαβxα,  u=εsα,  v=-(dα+δsα+t2sα+γsα),  w=εmα,  x=-(dα+δmα+γmα),  y=εxα,  z=-(dα+δxα+γxα),  m=p1αt1sα,  n=p2αt2sα,  j=p3αt2mα.The characteristic equation associated with above matrix is [38] (12)Δλ=JE0-λI=0,a-λ2λ2-r+z+t2xαe-λτλ-yt+z-t2xαe-λτr-λ2+h+x-t2mαe-λτλ-x+t2mαe-λτh+wq-λ2+e+vλ+uf-ve=0.Lemma 1.
IfR0<1, then the disease-free equilibrium E0 is locally asymptotically stable for τ=0.Proof.
Whenτ=0, the associated transcendental characteristic equation Δ(λ) of system (1) at E0 becomes Δ(λ)=P(λ)=0, and then the eigenvalues of the Jacobian matrix are (13)λ1,2=-d,λ3,4=r+z-t2xα±r2-2z-t2xαr+z-t2xα2+4yt2,λ5,6=x-t2mα+h±x-t2mα2-2x-t2mαh+h2+4wq2,λ7,8=v+e±v2-2ve+e2+4uf2,and by using Routh Hurwitz Theorem [28], these roots are negative or have negative real parts and all eigenvalues satisfy Matignon’s conditions [39], given by (|arg⁡λi|>απ/2) so the disease-free equilibrium E0 is locally asymptotically stable.Lemma 2.
LetR0<1, and then the disease-free equilibrium E0 is locally asymptotically stable for τ>0.Proof.
Let us considerτ>0, and we noted that second and third factor of the characteristic equation (12), which are (λ2-(r+z+t2xαe(-λτ))λ-yt+(z-t2xαe(-λτ))r) and (-λ2+(h+x-t2mαe(-λτ))λ-(x+t2mαe(-λτ))h+wq), have no pure imaginary roots for any value of the delay τ, if R0<1. Hence all the roots of the characteristic equation have negative real parts and we get that DFE is locally asymptotically stable regardless of the value of the delay and all eigenvalues satisfy Matignon’s conditions [39], given by (arg⁡λi>απ/2) so the disease-free equilibrium E0 is locally asymptotically stable.


### 3.2. Stability of the Endemic Equilibrium

System (1) has an endemic equilibrium if at least one of the infected variables is not zero. The expression “analytic” is complexity and not useful for our purposes. Consider the values of parameters from Table 3. Then the basic reproduction number is R0>1. The endemic equilibrium S=338.2,  Ls=0,  Lm=0,  Lx=2233.8,  Is=0,  Im=0,  Ix=4820.6,  R=62.0. The matrices A1 and A2 associated with the linearized system at the endemic equilibrium are computed as (14)A1=-9.02260.41070.41070.4107-0.2244-0.2244-0.22440.41070-2.2127000.632100000.8800-1.332700.8000.6321004.44750.14720.1472-1.4118-0.0791-0.07500.34321.052500.500000-1.272900000.01460.500-0.4415004.1762-0.1241-0.31810.6023-0.3504-0.3504-0.56750.78130.01880.01880.01880.0188-0.01030.01960.0237-2.6245,A2=000000000000000000000000000000000000000000000-t2mα00000000-t2xα0000000t2xα0. The transcendental characteristic equation Δλ=(λI-A1-e-τλA2) is given by(15)λ8+18.8862-0.0680e-τλλ7+1.2702×102-1.2255e-τλ+0.0012e-2τλλ6+4.2273×102-7.5306e-τλ+0.0198e-2τλλ5+7.7354×102-22.1163e-τλ+0.10986e-2τλλ4+7.9327×102-34.3423e-τλ+0.2751e-2τλλ3+4.3522×102-28.4367e-τλ+0.3354e-2τλλ2+1.1156×102-11.6693e-τλ+0.18817e-2τλλ9.6913-1.8329e-τλ+0.0367e-2τλ=0,when τ=0, and we have the following characteristic equation: (16)λ8+18.8182λ7+125.7957λ6+415.2192λ5+751.5336λ4+759.2028λ3+407.1187λ2+100.0789λ+7.8951=0.The roots of (16) are -0.14208, -0.41505, -1.57465, -1.40036, -1.01044, -9.15198, -2.56181+0.05148i, and -2.56181+0.05148i, and these roots are negative or have negative real parts and all eigenvalues satisfy Matignon’s conditions [39], given by (arg⁡λi>απ/2)  i=1,2,…,8 so the endemic equilibrium is locally asymptotically stable.Table 3
Parameter values of system (1).
Parameter
Value
bα
319
0
α
1
y
e
a
r
α
dα
0.3
8
α
1
y
e
a
r
α
βsα=βmα=βxα
1
4
α
1
y
e
a
r
α
λ
s
α
=
λ
m
α
=
λ
x
α
0
.
5
α
1
y
e
a
r
α
εsα=εmα=εxα
0
.
5
α
1
y
e
a
r
α
αr1,r2α
0.0
5
α
1
y
e
a
r
α
γsα=γmα=γxα
0
.
3
α
1
y
e
a
r
α
t
1
s
α
0.8
8
α
1
y
e
a
r
α
t2rα:r∈(s,m,x)
t
2
s
α
=
0.8
8
α;
t
2
m
α
=
t
2
x
α
=
0.03
4
α
1
y
e
a
r
α
σ
r
α
0.2
5
α
1
y
e
a
r
α
Prα
0.8
8
α
1
y
e
a
r
α
δ
r
α
0.04
5
α
1
y
e
a
r
αConsider now the caseτ>0; we noted that the roots of the characteristic equation (15) have no pure imaginary roots for any value of the delay τ, if R0>1. Hence all the roots of the characteristic equation have negative real parts and all eigenvalues satisfy Matignon’s conditions [39]. Therefore, the endemic equilibrium is locally asymptotically.


## 3.1. Stability of the Disease-Free Equilibrium

IfIs(t)=Im(t)=Ix(t)=0 then Ls(t)=Lm(t)=Lx(t)=0, R(t)=0, and S(t)=bα/dα. Then the disease-free equilibrium (DFE) is E0={(bα/dα,0,0,0,0,0,0,0)}.Let us consider the coordinate transformation:s(t)=S(t)-S¯,  ls(t)=Ls(t)-Ls¯,  lm(t)=Lm(t)-Lm¯,  lx(t)=Lx(t)-Lx¯,  is(t)=Is(t)-Is¯,  im(t)=Im(t)-Im¯,  ix(t)=Ix(t)-Ix¯,  r(t)=R(t)-R¯. The corresponding characteristic equation for DFE is given as follows:(11)JE0-λI=λ-a000bcd100λ-e00f0000gλ-h0pq00000λ-r0st00u00λ-v00000w00λ-x+t2mαe-λτ00000y00λ-z+t2xαe-λτ00m00njt2xαe-λτλ-a,where a=-dα,  b=-βsα,  c=-βmα,  d1=-βxα,  e=-(dα+εsα+t1sα),  f=γsα+λsαβsα,  g=(1-p1α)t1sα,  h=-(dα+εmα),  p=(1-p2α)t2sα,  q=γmα+λmαβmα,  r=-(dα+εxα),  s=(1-p3α)t2mα,  t=γxα+λxαβxα,  u=εsα,  v=-(dα+δsα+t2sα+γsα),  w=εmα,  x=-(dα+δmα+γmα),  y=εxα,  z=-(dα+δxα+γxα),  m=p1αt1sα,  n=p2αt2sα,  j=p3αt2mα.The characteristic equation associated with above matrix is [38] (12)Δλ=JE0-λI=0,a-λ2λ2-r+z+t2xαe-λτλ-yt+z-t2xαe-λτr-λ2+h+x-t2mαe-λτλ-x+t2mαe-λτh+wq-λ2+e+vλ+uf-ve=0.Lemma 1.
IfR0<1, then the disease-free equilibrium E0 is locally asymptotically stable for τ=0.Proof.
Whenτ=0, the associated transcendental characteristic equation Δ(λ) of system (1) at E0 becomes Δ(λ)=P(λ)=0, and then the eigenvalues of the Jacobian matrix are (13)λ1,2=-d,λ3,4=r+z-t2xα±r2-2z-t2xαr+z-t2xα2+4yt2,λ5,6=x-t2mα+h±x-t2mα2-2x-t2mαh+h2+4wq2,λ7,8=v+e±v2-2ve+e2+4uf2,and by using Routh Hurwitz Theorem [28], these roots are negative or have negative real parts and all eigenvalues satisfy Matignon’s conditions [39], given by (|arg⁡λi|>απ/2) so the disease-free equilibrium E0 is locally asymptotically stable.Lemma 2.
LetR0<1, and then the disease-free equilibrium E0 is locally asymptotically stable for τ>0.Proof.
Let us considerτ>0, and we noted that second and third factor of the characteristic equation (12), which are (λ2-(r+z+t2xαe(-λτ))λ-yt+(z-t2xαe(-λτ))r) and (-λ2+(h+x-t2mαe(-λτ))λ-(x+t2mαe(-λτ))h+wq), have no pure imaginary roots for any value of the delay τ, if R0<1. Hence all the roots of the characteristic equation have negative real parts and we get that DFE is locally asymptotically stable regardless of the value of the delay and all eigenvalues satisfy Matignon’s conditions [39], given by (arg⁡λi>απ/2) so the disease-free equilibrium E0 is locally asymptotically stable.


## 3.2. Stability of the Endemic Equilibrium

System (1) has an endemic equilibrium if at least one of the infected variables is not zero. The expression “analytic” is complexity and not useful for our purposes. Consider the values of parameters from Table 3. Then the basic reproduction number is R0>1. The endemic equilibrium S=338.2,  Ls=0,  Lm=0,  Lx=2233.8,  Is=0,  Im=0,  Ix=4820.6,  R=62.0. The matrices A1 and A2 associated with the linearized system at the endemic equilibrium are computed as (14)A1=-9.02260.41070.41070.4107-0.2244-0.2244-0.22440.41070-2.2127000.632100000.8800-1.332700.8000.6321004.44750.14720.1472-1.4118-0.0791-0.07500.34321.052500.500000-1.272900000.01460.500-0.4415004.1762-0.1241-0.31810.6023-0.3504-0.3504-0.56750.78130.01880.01880.01880.0188-0.01030.01960.0237-2.6245,A2=000000000000000000000000000000000000000000000-t2mα00000000-t2xα0000000t2xα0. The transcendental characteristic equation Δλ=(λI-A1-e-τλA2) is given by(15)λ8+18.8862-0.0680e-τλλ7+1.2702×102-1.2255e-τλ+0.0012e-2τλλ6+4.2273×102-7.5306e-τλ+0.0198e-2τλλ5+7.7354×102-22.1163e-τλ+0.10986e-2τλλ4+7.9327×102-34.3423e-τλ+0.2751e-2τλλ3+4.3522×102-28.4367e-τλ+0.3354e-2τλλ2+1.1156×102-11.6693e-τλ+0.18817e-2τλλ9.6913-1.8329e-τλ+0.0367e-2τλ=0,when τ=0, and we have the following characteristic equation: (16)λ8+18.8182λ7+125.7957λ6+415.2192λ5+751.5336λ4+759.2028λ3+407.1187λ2+100.0789λ+7.8951=0.The roots of (16) are -0.14208, -0.41505, -1.57465, -1.40036, -1.01044, -9.15198, -2.56181+0.05148i, and -2.56181+0.05148i, and these roots are negative or have negative real parts and all eigenvalues satisfy Matignon’s conditions [39], given by (arg⁡λi>απ/2)  i=1,2,…,8 so the endemic equilibrium is locally asymptotically stable.Table 3
Parameter values of system (1).
Parameter
Value
bα
319
0
α
1
y
e
a
r
α
dα
0.3
8
α
1
y
e
a
r
α
βsα=βmα=βxα
1
4
α
1
y
e
a
r
α
λ
s
α
=
λ
m
α
=
λ
x
α
0
.
5
α
1
y
e
a
r
α
εsα=εmα=εxα
0
.
5
α
1
y
e
a
r
α
αr1,r2α
0.0
5
α
1
y
e
a
r
α
γsα=γmα=γxα
0
.
3
α
1
y
e
a
r
α
t
1
s
α
0.8
8
α
1
y
e
a
r
α
t2rα:r∈(s,m,x)
t
2
s
α
=
0.8
8
α;
t
2
m
α
=
t
2
x
α
=
0.03
4
α
1
y
e
a
r
α
σ
r
α
0.2
5
α
1
y
e
a
r
α
Prα
0.8
8
α
1
y
e
a
r
α
δ
r
α
0.04
5
α
1
y
e
a
r
αConsider now the caseτ>0; we noted that the roots of the characteristic equation (15) have no pure imaginary roots for any value of the delay τ, if R0>1. Hence all the roots of the characteristic equation have negative real parts and all eigenvalues satisfy Matignon’s conditions [39]. Therefore, the endemic equilibrium is locally asymptotically.


## 4. NSFD for Fractional Delay Differential Equations

In this section, we apply NSFD method with GLFDs to obtain the discretization of the delay fractional multistrain TB model (1), which will yield the following equations: (17)∑j=0n+1ωjαSn+1-j=bα-dαSn+1-βsαSn+1IsnNn-βmαSn+1ImnNn-βxαSn+1IxnNn,∑j=0n+1ωjαLsn+1-j=λsαβsαSn+1IsnNn+σsαλsαβsαRn+1IsnNn+γsαIsn-αssαβsLsn+1IsnNn-αsxαβxαLsn+1IxnNn-dα+εsα+t1sαLsn+1-αsmαβmαLsn+1ImnNn,∑j=0n+1ωjαLmn+1-j=λmαβmαSn+1ImnNn+σmαλmαβmαRn+1ImnNn+λmααsmαβmαLsn+1ImnNn+γmαImn+t1sαLsn+1-P1αt1sαLsn+1+t2sαIsn-P2αt2sαIsn-αmmαβmαLmn+1ImnNn-αmxαβxαLmn+1IxnNn-dα+εmαLmn+1,∑j=0n+1ωjαLxn+1-j=λxαβxαSn+1IxnNn+σxαλxαβxαRn+1IxnNn+λxααsxαβxαLsn+1IxnNn+γxαIxn+λxααmxαβxαLmn+1IxnNn+t2mαImn-P3αt2mαImn-αxxαβxαLxn+1IxnNn-dα+εxαLxn+1,∑j=0n+1ωjαIsn+1-j=αssαβsαLsn+1IsnNn+1-λsαβsαSn+1IsnNn+σsαRn+1IsnNn+εsαLsn+1-dα+δsαIsn+1-γsα+t2sαIsn,∑j=0n+1ωjαImn+1-j=αmmαβmαLmn+1ImnNn+1-λmαβmαSn+1ImnNn+σmαRn+1ImnNn+αsmαLsn+1ImnNn+εmαLmn+1-dα+δmαImn+1-γmαIm-t2mαImn-κ,∑j=0n+1ωjαIxn+1-j=αxxαβxαLxn+1IxnNn+1-λxαβmSn+1IxnNn+σxRn+1IxnNn+αmxLxn+1ImnNn+εxαLxn+1-dα+δxαIxn+1-γxαIx-t2xαIxn-κ,∑j=0n+1ωjαRn+1-j=P1αt1sαLsn+1+P2αt2sαIsn+P3αt2mαImn+t2xαIxn-κ-dαRn+1-σsαβsαRn+1IsnNn-σmαβmαRn+1ImnNn-σxαβxαRn+1IxnNn,where (18)Nn=Sn+Lsn+Lmn+Lxn+Isn+Imn+Ixn+Rn,and ω0α=(φi(h))-α,  i=1,2,…,8,  n=-κ,-κ+1,…,0,1, where the nonlocal approximations are used for the nonlinear terms and the following functions of denominator:(19)φ1h=edαh-1dα,φ2h=edα+εsα+t1sαh-1dα+εsα+t1sα,φ3h=edα+εmαh-1dα+εmα,φ4h=edα+εxαh-1dα+εxα,φ5h=1-e-dα+δsαhγsα+t2sα,φ6h=1-e-dα+δmαhγmα+t2mα,φ7h=1-e-dα+δxαhγxα+t2xα,φ8h=edαh-1dα.Then we obtain(20)Sn+1=bα-∑j=1n+1ωjαSn+1-jφ1h-α+dα+βsαIsn+βmαImn+βxαIxn/Nn,Lsn+1=βsαIsn/NnλsαSn+1+σsαRn+1+γsαIsn-∑j=1n+1ωjαLsn+1-jφ2h-α+dα+t1sα+εsα+1/NnαssαβsαIsn+αsmαβmαImn+αsxαβxαIxn,Lmn+1=βmαλmαImn/NnSn+1+σmαRn+1+αsmαLsn+1+γmαImn+t1sαLsn+11-P1αφ3h-α+dα+εmα+1/NnαmmαβmαImn+αmxαβxαIxn+t2sαIsn1-P2α-∑j=1n+1ωjαLmn+1-jφ3h-α+dα+εmα+1/NnαmmαβmαImn+αmxαβxαIxn,Lxn+1=βxαλxαIxn/NnSn+1+σxαRn+1+αsxαLsn+1+αmxαLmn+1+t2mαImn1-P3αφ4h-α+dα+εxα+1/NnαxxαβxαIxn+γxαIxn-∑j=1n+1ωjαLxn+1-jφ4h-α+dα+εxα+1/NnαxxαβxαIxn,Isn+1=φ5hβsαIsn/NnαssαLsn+1+1-λsαSn+1+σsαRn+1φ5h-α+dα+δsα-γsα+t2sαIsn+εsαLsn+1-∑j=1n+1ωjαIsn+1-jφ5h-α+dα+δsα,Imn+1=βmαImn/NnαmmαLmn+1+1-λmαSn+1+σmαRn+1+αsmαLsn+1φ6h-α+dα+δmα-γmImn-t2mαImn-κ+εmαLmn+1-∑j=1n+1ωjαImn+1-jφ6h-α+dα+δmα,Ixn+1=βxαIxn/NnαxxαLxn+1+1-λxαSn+1+σxαRn+1+αsxαLsn+1+αmxαLmn+1φ7h-α+dα+δxα-γxαIxn-t2xαIxn-κ+εxαLxn+1-∑j=1n+1ωjαIxn+1-jφ7h-α+dα+δxα,Rn+1=t1sαP1αLsn+1+P2αt2sαIsn+t2mαP3αImn+t2xαIxn-κ-∑j=1n+1ωjαRn+1-jφ8h-α+dα+1/NnσsαβsαIsn+σmαβmαImn+σxαβxαIxn.


## 5. Numerical Results and Simulations

In this section, we show the effectiveness of the numerical technique for delay fractional differential equations. Throughout this section, all simulations are performed with initial conditionsS0,Ls0,Lm0,Lx0,Is0,Im0,Ix0,R0=(5000,50,50,50,30,30,30,60), with the parameters in Table 3. The approximate solutions of the proposed system are given in Figures 1–12 at different values of τ and α. Figure 1 shows the behavior of the approximate solutions of R(t) in two cases with and without delay using dde23 at α=1 and τ=0.3. In Figure 2, we use the same data in Figure 1 and use NSFDM; we noted that the number of individuals R(t) increases in the case of nondelay, that is, the delays in diagnosis and commencement of treatment to the individuals of Im and Ix causing a shortage in the number of individuals of R(t). Figure 3 shows the relationship between Im(t) and Im(t-τ) and chaotic attractors at τ=0.1 and α=1. Figure 4, shows the relationship between Ix(t) and Ix(t-τ) at τ=0.1 in case of integer order. Figures 5 and 6 show the relationship between Im(t) and Im(t-τ) and Ix(t) with Ix(t-τ), respectively, in case of fraction order where α=0.95. Figures 7 and 8 show the approximate solutions for Im(t-τ) and Ix(t-τ) at τ=2,  α=0.98 by using NSFDM. Figures 9 and 10 show the approximate solutions of different τ in both fraction and integer cases; we noted that increasing the value of τ causes decreasing the values of R(t). Figures 11 and 12 show the behavior of the approximate solutions with different value of α, which are given to demonstrate how the fractional model is a generalization of the integer order model.Figure 1
The approximate solution ofR(t) with τ=0.3, using dde23.Figure 2
The approximate solution ofR(t) with τ=0.3, using NSFDM, α=1.Figure 3
The relationship betweenIm(t),  Im(t-τ) with τ=0.1,  α=1 using NSFDM.Figure 4
The relationship betweenIx(t),  Ix(t-τ) with τ=0.1,  α=1 using NSFDM.Figure 5
The relationship betweenIm(t),  Im(t-τ) with τ=0.1,  α=0.95 using NSFDM.Figure 6
The relationship betweenIx(t),  Ix(t-τ) with τ=0.1,  α=0.95 using NSFDM.Figure 7
The approximate solutionsIx(t-τ) with τ=2,  α=0.98 using NSFDM.Figure 8
The approximate solutionsIm(t-τ) with τ=2,  α=0.98 using NSFDM.Figure 9
The approximate solutions with differentτ and α=0.9 using NSFDM.Figure 10
The approximate solution ofR(t) with different τ, using NSFDM, α=1.Figure 11
The approximate solutions withτ=0.4 and different α using NSFDM.Figure 12
The approximate solutions withτ=0.4 and different α using NSFDM.


## 6. Conclusion

Fractional models have the potential to describe more complex dynamics than the integer models and can include easily the memory effect present in many real world phenomena. In this paper, multistrain TB model of fractional order derivatives with time delay memory is presented. A nonstandard numerical scheme is introduced to numerically study the approximate solution of proposed model problem. The obtained results show that the delays in diagnosis and commencement of treatment to the individuals ofIm and Ix cause a shortage in the number of individuals of R(t). The approximate solution of proposed model changes when τ and α take different values. Some figures are given to demonstrate how the fractional delay model is a generalization of the integer order model. It is concluded that NSFDM can be applied to solve such fractional delay differential equations simply and effectively. All results are obtained by using MATLAB programming.



---

*Source: 1047384-2017-07-31.xml*