# Mathematical Analysis of a General Two-Patch Model of Tuberculosis Disease with Lost Sight Individuals

Authors: Abdias Laohombé; Isabelle Ngningone Eya; Jean Jules Tewa; Alassane Bah; Samuel Bowong; Suares Clovis Oukouomi Noutchie
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

A two-patch model,SEi1,…,EinIiLi,  i=1,2, is used to analyze the spread of tuberculosis, with an arbitrary number n of latently infected compartments in each patch. A fraction of infectious individuals that begun their treatment will not return to the hospital for the examination of sputum. This fact usually occurs in sub-Saharan Africa, due to many reasons. The model incorporates migrations from one patch to another. The existence and uniqueness of the associated equilibria are discussed. A Lyapunov function is used to show that when the basic reproduction ratio is less than one, the disease-free equilibrium is globally and asymptotically stable. When it is greater than one, there exists at least one endemic equilibrium. The local stability of endemic equilibria can be illustrated using numerical simulations. Numerical simulation results are provided to illustrate the theoretical results and analyze the influence of lost sight individuals.


---

## Body

## 1. Introduction

For a given system, the focus in qualitative mathematical epidemiology is the long-term dynamics. The simplest possible attractor is a globally and asymptotically stable equilibrium. Equilibrium can be shown to be globally and asymptotically stable, using Poincaré-Bendixson theory [1], Bendixson’s negative criterion [2, 3], or the generalized version of Dulac [4]. Another method of Li and Muldowney [5–7] for demonstrating global stability in n dimensions has been developed more recently, with applications in three [8, 9] and four dimensions [10, 11]. For higher-dimensional systems, the theory of quadratic forms [12] or Lyapunov’s method can be used [13, 14]. Lyapunov’s method requires finding a function V such that the flow always crosses the level sets from higher values of V to lower values. When such a function can be found, then any isolated minimum of the function is a stable equilibrium of the flow. In this paper, the stability of a 2n+6-dimensions system is investigated using Lyapunov-LaSalle functions and quadratic forms. The function V=∑i=1nai(xi-xi*ln⁡xi) has a long story in epidemiology [15–20]. Volterra himself originally discovered this function, although he did not use the vocabulary and the theory of Lyapunov functions. Since epidemic models are Lotka-Volterra-like models, the pertinence of this function is not surprising.The issue of modeling tuberculosis motivates the model studied in this paper. It is an extension with two patches andn latent classes of the SEIL model in [21], where L denotes the lost sight individuals. These lost sight individuals usually occur in sub-Saharan Africa. For example, according to the Cameroonian National Program of Fight against Tuberculosis, about 10% of infectious individuals that begun their therapy treatment lost sight. Therefore, this fact cannot be neglected in the modeling of TB. By allowing for an arbitrary number of latently infected compartments, the model allows for the approximation of a wide class of distributions of latency durations. This is of particular importance for tuberculosis since latency may last for years or even decades.In the model given here, there are migrations between people of two patches, who have the same epidemiological characteristics. We introduce a direct transfer from the class of susceptible individuals toward the compartment of infectious individuals in each patch. The reason is that there are some individuals with bad immune effector protection, due to many reasons. We also incorporate a transfer from the infectious class to the exposed class to take into account the fact that some infectious individuals apparently recover but actually harbor TB bacteria. Let us give now the outlines of the paper. In Section2, the model is constructed; the variables and parameters of the model are explained. In Section 3, the mathematical properties of the model are given. We present the computation thresholds as basic reproduction numbers, which are bifurcation parameters R0i, using the same method as in [22]. The equilibria of the model are computed: equilibrium without disease, endemic equilibria, or coexistence equilibria. The stability of each equilibrium is investigated, using the bifurcation parameters R0i. In Section 4, numerical simulations are done to illustrate the results. In Section 5, we give the conclusion.


## 2. Model Construction

The model consisted of two patches. Each patch denotes a given population and the two populations are very much closed. Based on epidemiological status, the population of a given patchi(i=1,2) is divided into n+3 classes: susceptible (Si), latently infected (Ei), infectious (Ii), and lost sight (Li). All recruitments in a given patch i are into the susceptible class and occur at a constant rate Λi. The rate constant for nondisease related death is μi; thus 1/μi is the average lifetime. Infectious and lost sight classes of a patch i have addition death rates due to the disease with rates constants di and li, respectively. Since we do not know if lost sight class is recovered, died, or is still infectious, we assume that a fraction δi of them is still infectious and can transmit disease to susceptible class. Transmission ofMycobacterium tuberculosis occurs following adequate contacts between susceptible and infectious or lost sight classes that continue to have disease. We assume that infected individuals are not infectious and thus are not capable of transmitting bacteria. We use the universal incidence expressions βiSiIi/Ni and βiδiSiLi/Ni to indicate successful transmission ofM. tuberculosis due to nonlinear contacts dynamics in the population by infectious and lost sight classes, respectively. A fraction pi of the newly infected individuals is assumed to undergo fast progression directly to the infectious class, while the remainder are latently infected and enter the latent class. Once latently infected withM. tuberculosis, an individual will remain so for life unless reactivation occurs. To account for treatment, we define rijEi, with i=1,2 and j=1,…,n, as the fraction of infected individuals receiving effective chemoprophylaxis, and γi as the rate of effective per capita therapy. We assume that chemoprophylaxis of latently infected individuals Eij reduces their reactivation at rate rij. Thus, a fraction (1-rij)Ei of infected individuals who do not receive effective chemoprophylaxis becomes infectious with a rate constant ki, so that 1/ki represents the average latent period. Thus, individuals leave the class Eij to Ii at rate ki(1-rij). After receiving an effective therapy, infectious individuals can spontaneously recover from the disease with a rate constant γi, entering the infected class Ein. A fraction si(1-rij)Ii of infectious individuals that began their treatment will not return to the hospital for the examination of sputum. After some time, some of them will return to the hospital with the disease at a constant rate αi. This can be the situation in many African countries or refugees camps in Africa or elsewhere.The transfer diagram of the model is given by Figure1.Figure 1
Flow chart of a two-patch model of tuberculosis with lost sight individuals in sub-Saharan Africa, with migrations at rateai between susceptible individuals of two populations, at rate bij between latently infected individuals of two populations, at rate ci between infectious individuals in care centre, and at rate fi for infectious lost sight individuals. An individual of a given population (patch) has contacts only with individuals of the same population.This yields the following set of differential equations:(1)dS1dt=Λ1-(μ1+a1)S1-β1(I1+δ1L1)S1N1+a2S2,dE11dt=β1(1-p1)(I1+δ1L1)S1N1+b21E21-(μ1+b11+(1-r11)k11)E11,dE12dt=(1-r11)k11E11+b22E22-(μ1+b12+(1-r12)k12)E12,⋮dE1ndt=(1-r1,n-1)k1,n-1E1,n-1+γ1I1+b2nE2n-(μ1+b1n+(1-r1n)k1n)E1n,dI1dt=β1p1(I1+δ1L1)S1N1+(1-r1n)k1nE1n+α1L1+c2I2-(μ1+d1+γ1+(1-v1)s1+c1)I1,dL1dt=(1-v1)s1I1+f2L2-(μ1+l1+α1+f1)L1,dS2dt=Λ2-(μ2+a2)S2-β2(I2+δ2L2)S2N2+a1S1,dE21dt=β2(1-p2)(I2+δ2L2)S2N2+b11E11-(μ2+b21+(1-r21)k21)E21,dE22dt=(1-r21)k21E21+b12E12-(μ2+b22+(1-r22)k22)E22,⋮dE2ndt=(1-r2,n-1)k2,n-1E2,n-1+γ2I2+b1nE1n-(μ2+b2n+(1-r2n)k2n)E2n,dI2dt=β2p2(I2+δ2L2)S2N2+(1-r2n)k2nE2n+α2L2+c1I1-(μ2+d2+γ2+(1-v2)s2+c2)I2,dL2dt=(1-v2)s2I2+f1L1-(μ2+l2+α2+f2)L2,
where N1=S1+E11+⋯+E1n+I1+L1 and N2=S2+E21+⋯+E2n+I2+L2.System (1) can also be written as
(2)dx1dt=Λ1-(μ1+a1)x1-β1x1N1〈en+1n+2+δ1en+2n+2Y1〉,dx1dt=Λ2-(μ2+a2)x2-β2x2N2〈en+1n+2+δ2en+2n+2Y2〉,dY1dt=β1x1N1〈en+1n+2+en+2n+2Y1〉Q1+A11Y1+A12Y2,dY2dt=β2x2N2〈en+1n+2+en+2n+2Y1〉Q2+A21Y1+A22Y2,
where 〈·/·〉 is the scalar usual product in IR+n+2, Q1=((1-p1),0,…,0,p1)T, Q2=((1-p2),0,…,0,p2)T, x=(x1,x2)T, Y1=(E11,…,E1n,I1,L1)T=(y11,…,y1n,y1,n+1,y1,n+2)T, Y1=(E21,…,E2n,I2,L2)T = (y21,…,y2n,y2,n+1,y2,n+2)T, (ein+2) is the canonical basis of IR+n+2, the matrices A11 and A22 are Metzler stable [23], and A12 and A21 are Metzler stable and given, respectively, by(3)A11=(-m110⋯0000(1-r11)k11-m120⋯000⋯⋯⋯⋯⋯⋯⋯0⋯0(1-r1,n-1)k1,n-1-m1nγ1000⋯0(1-r1,n)k1,n-mI1α1000⋯0(1-v1)s1-mL1),A22=(-m210⋯0000(1-r21)k21-m220⋯000⋯⋯⋯⋯⋯⋯⋯0⋯0(1-r2,n-1)k2,n-1-m2nγ2000⋯0(1-r2,n)k2,n-mI2α2000⋯0(1-v2)s2-mL2),A12=(b2100…000b220…000…0b2n000…00c200…000f2),A21=(b1100…000b120…000…0b1n000…00c100…000f1),where mi1=(μi+bi1+(1-ri1)ki1), min=(μi+bin+(1-rin)kin), mIi=(μi+di+γi+(1-vi)si+ci), and mLi=(μi+li+αi+fi), i=1,2.Let us give another expression of system (1). System (2) can be written as
(4)dx(t)dt=Λ+Dxx-∑i=12〈Biy(t)〉〈ei2/x(t)〉Niei2,dy(t)dt=∑i=12〈Biy(t)〉〈ei2/x(t)〉NiKi+Ayy,
where 〈·/·〉 is now the usual scalar product in IR2n+4, x=(x1,x2)T, e12=(1,0)T, e22=(0,1)T, Λ=(Λ1,Λ2)T, y=(E11,…,E1n,I1,L1,E21,…,E2n,I2,L2)T, Dx=(-(μ1+a1)a2a1-(μ2+a2)), B1=(0,…,0︸n,β1,β1δ1,0,…,0︸n+2)T, B2=(0,…,0︸n+2,0,…,0︸n,β2,β2δ2)T, K1=(1-p1,0,…,0︸n-1,p1,0,0,…,0︸n+2)T, K2=(0,…,0︸n+2,1-p2,0,…,0︸n-1,p2,0)T, and Ay=(A11A12A21A22).


## 3. Mathematical Properties

### 3.1. Positivity of the Solutions

Since the variables considered here are nonnegative quantities, we have to be sure that their values are always nonnegative.Theorem 1.
The nonnegative orthantIR+2n+6 is positively invariant by (1). This means that every trajectory, which begins in the positive orthant, will stay inside.Proof.
System (4) can be written in the following form:
(5)dx(t)dt=Λ+Dxx-∑i=12〈Biy(t)〉〈ei2/x(t)〉Niei2,dy(t)dt=(∑i=12BiT〈ei2/x(t)〉NiKi+Ay)y(t).
Since x(t)≥0 and Ni(t)≥0, the matrix (∑i=12BiT(〈ei2/x(t)〉/Ni)Ki+Ay) is Metzler. It is well known that linear Metzler system lets the nonnegative orthant invariant [23]. Moreover dx(t)/dt≤Λ+Dxx and Dx is a Metzler matrix. These prove the positive invariance of the nonnegative orthant IR+2n+6 by (1).


### 3.2. Boundedness of the Trajectories

Lemma 2.
The simplex{(S1,E11,…,E1n,I1,L1,S2,E21,…,E2n,I2,L2)∈IR+2n+6/S1+E11+⋯E1n+I1 + L1+S2+E21+⋯E2n+I2+L2≤((Λ1+Λ2)/μ)+ε}, is a compact forward and absorbing set for (1).Proof.
From system (1), if N(t)=N1(t)+N2(t) denotes the entire population of (1), one has N(t)=S1+E11+⋯+E1n+I1+L1+S2+E21+⋯+E2n+I2+L2.
Then, the dynamics ofN(t) at any time t is given by
(6)dN(t)dt=Λ1+Λ2-μ1N1(t)-μ2N2(t)-d1I1-l1L1-d2I2-l2L2.
Thus, dN(t)/dt≤Λ1+Λ2-μN(t), where μ=min⁡(μ1,μ2). It follows that limSupt→+∞N(t)=(Λ1+Λ2)/μ.In the simplexΓε, (1) is mathematically well posed. The following lemma also holds.Lemma 3.
The simplexγε={(x,y)∈Γε/x≤x*} is a compact forward invariant set for (1).


### 3.3. Local Stability of the Disease-Free Equilibrium

Many epidemiological models have a threshold condition which can determine whether an infection will be eliminated from the population or become endemic. The basic reproduction number,R0, is defined as the average number of secondary infections produced by an infected individual in a completely susceptible population [16]. As discussed in [24, 25], R0 is a simply normalized bifurcation (transcritical) parameter for epidemiological models, such that R0 implies that the endemic steady state is stable (i.e., the infection persists), and R0 implies that the uninfected steady state is stable (i.e., the infection can be eliminated from the population).Equation (1) has a disease-free equilibrium given by K2=(S1*,0,…,0︸n+2,S2*,0,…,0︸n+2)T which always exists in the nonnegative orthant IR+2n+6. The explicit expressions of S1* and S2* are
(7)S1*=(μ2+a2)Λ1+a2Λ2μ1μ2+μ1a2+μ2a1,S2*=(μ1+a1)Λ2+a1Λ1μ1μ2+μ1a2+μ2a1.
The disease-free equilibrium can also be denoted by (x*,0) since it is the solution of equation Λ+Dxx*=0 of the compact system (5).Lemma 4.
Using the same method as in [17], the basic reproduction ratio of (1) is R0=ρ(K1B1T(-Ay)-1+K2B2T(-Ay)-1), where ρ denotes the spectral radius.Proof.
The expressions, which are coming from the other compartment, due to contamination, are given by the following matrix:(8)F=∂∂y(∑i=12〈Biy(t)〉〈ei2/x(t)〉NiKi)(x*,0)=∑i=12〈ei2/x*〉Ni*KiBiT.
The expressions, which are coming from the other compartment, due to reasons different from contamination, are given by V=-Ay. The next generation matrix, since 〈ei2/x*〉=Ni*, is
(9)FV-1=(∑i=12〈ei2/x*〉Ni*KiBiT)(-Ay)-1=(∑i=12KiBiT)(-Ay)-1.
We can observe that, since R0 is the largest eigenvalue of the next generation matrix,
(10)R0=ρ(FV-1)=ρ(K1B1T(-Ay)-1+K2B2T(-Ay)-1),
where ρ denotes the spectral radius.Consequently, from Theorem 2 of [22], the following lemma holds.Lemma 5.
The disease-free equilibrium of (1) is locally and asymptotically stable whenever R0≤1 and unstable if R0>1.Remark 6.
This lemma shows that ifR0≤1, a small flow of infectious individuals will not generate large outbreaks of the disease. To control the disease independently of the initial total number of infectious individuals, a global asymptotic stability property has to be established for the DFE when R0≤1.


### 3.4. Global Stability of the Disease-Free Equilibrium

The following result helps to determine the stability and is related to LaSalle’s principle [13]. Consider the differential equation
(11)dx(t)dt=X(x(t)),X(x0)=0,
where X is a C1 function defined on an open set of IRn containing the closure Ω- of a positively invariant set Ω such that the equilibrium x0 is in Ω-. The following lemma holds.Lemma 7.
We assume that system (11) is point dissipative [1] on Ω. In other words there exists a compact set K⊂Ω such that, for any y∈Ω, there exists a time t(y0)>0 such that, for any time t>t(y0), the trajectory with initial condition y0 is in the interior of K. If there exists a C1 function V≥0 defined on Ω-, then(1)
dV(t)/dt=〈δV(x)/X(x)〉≤0, for all x∈Ω;(2)
the greatest invariant setL contained in E={x∈Ω-,dV(t)/dt=0} is contained in a positively invariant set on which the restriction of system (11) is globally and asymptotically stable on L at x0.Then, x0 is a globally and asymptotically stable equilibrium of system (11) on Ω.Theorem 8.
WhenR0≤1 (this implies R01≤1 and R02≤1, where R0i is the basic reproduction number for the patch i), the disease-free equilibrium (DFE), when it is unique, is globally and asymptotically stable in Γε, since it is the unique equilibrium. This implies the global asymptotic stability of the DFE on the nonnegative orthant IR+2n+6; that is, the disease dies out in both two populations.Remark 9.
At least one endemic equilibrium can exist and coexist with the disease-free equilibrium whenR0≤1. In this case the DFE cannot be globally asymptotically stable.Proof.
Following the same method as in [23], system (4) can be written in a pseudotriangular form as
(12)dx(t)dt=A1(x,y)(x-x*)+A12(x,y)y=f(x,y),dy(t)dt=A2(x,y)y=g(x,y),
where f(x*,0)=Λ+Dxx*=0, g(x*,0)=0, A1(x,y)=Dx∈IR2x2, A12(x,y) = -∑i=12(〈ei2/x(t)〉/Ni)ei2BiT∈IR2x(2n+4), and A2(x,y)=∑i=12(〈ei2/x(t)〉/Ni)KiBiT+Ay∈IR(2n+4)x(2n+4). The Jacobian matrix of system (4) at disease-free equilibrium (x*,0) is A2(x,y)=∑i=12(〈ei2/x(t)〉/Ni)KiBiT+Ay∈IR(2n+4)x(2n+4). The Jacobian matrix of system (4) at disease-free equilibrium (x*,0) is J=(A1(x*,0)A12(x*,0)0A2(x*,0))=(Dx-∑i=12ei2BiT0∑i=12KiBiT+Ay). The matrix Dx is clearly a Metzler stable matrix and, using a result in [8], there exists a vector u=(u1,u2)>0 such that Dxu≤0. The matrix A-2=∑i=12KiBiT+Ay is a Metzler and irreducible matrix, which is stable if R0≤1. With this condition satisfied and using the same previous result in [22], the stability modulus α(A-2) of the matrix is in the spectrum of A-2 (Perron-Frobenius theorem, proof of Lemma 2) and there exists a vector c>0∈IR2n+4 such that cTA-2=α(A-2)cT, with α(A-2)≤0. Let us show the global stability of the DFE.
System (12) is in pseudotriangular form since the matrix A2(x,y) depends on (x,y). There are many results for the stability of triangular systems [21–25]. Using LaSalle’s principle [13] one can obtain attractivity of equilibrium. But, for nonlinear systems, attractivity does not generally implies stability. We can use now the result of Lemma 7. Let us consider the following candidate Lyapunov function:
(13)VDFE(x(t),y(t))=〈(0,c)(x(t),y(t))T〉.
The function VDFE(x(t),y(t)) is positive and it is time derivative along the trajectories of (4) which gives
(14)dV(x(t),y(t))dt=(∑i=12〈Biy(t)〉〈ei2/x(t)〉NiKi+Ayy)c≤cT(∑i=12KiBiT+Ay)y≤cTα(A-2)y.
Then, R0≤1 ensures that dVDFE(x(t),y(t))/dt≤0 for all y≥0 and that dVDFE(x(t),y(t))/dt=0 holds when (x(t),y(t))=(x*,0). This proves the global asymptotic stability on IR+2n+6 [14]. This achieves the proof that the DFE is globally and asymptotically stable.


### 3.5. Existence of Endemic Equilibrium

Definition 10.
An equilibrium for a multipatch model as (1) is called endemic equilibrium when the two populations coexist (the density of each compartment is different from zero) at this equilibrium.Lemma 11.
An endemic equilibrium(x-,y-) of (1) can be determined, using an equation F(x-)=1.Proof.
An endemic equilibrium(x-,y-) of (1) is obtained by setting the right-hand side of (4) which equals zero, giving
(15)Λ+Dxx-=∑i=12〈Biy-(t)〉〈ei2/x-(t)〉N-iei2,∑i=12〈Biy-(t)〉〈ei2/x-(t)〉N-iKi+Ayy-=0.
Multiplying the second equation of (15) by (-Ay)-1 gives
(16)y-=〈B1y-(t)〉〈e12/x-(t)〉N-1(-Ay)-1K1+〈B2y-(t)〉〈e22/x-(t)〉N-2(-Ay)-1K2.
From (16),
(17)〈B1y-〉=g11〈B1y-〉〈e12/x-(t)〉N-1+g12〈B2y-〉〈e22/x-(t)〉N-2,〈B2y-〉=g21〈B1y-〉〈e12/x-(t)〉N-1+g22〈B2y-〉〈e22/x-(t)〉N-2,
where
(18)g11=〈B1(-Ay)-1K1〉,g12=〈B1(-Ay)-1K2〉,g21=〈B2(-Ay)-1K1〉,g22=〈B2(-Ay)-1K2〉.
From the second equation of (17),
(19)〈B2y-〉=g21(〈e12/x-(t)〉/N-1)1-g22(〈e22/x-(t)〉/N-2)〈B1y-〉.
Using the expression of 〈B2/y-〉 in (16) gives
(20)y-=〈B1y-〉〈e12/x-(t)〉N-1(-Ay)-1K1+g21(〈e12/x-(t)〉/N-1)1-g22(〈e22/x-(t)〉/N-2)〈B1y-〉×〈e22/x-(t)〉N-2(-Ay)-1K2.
Then
(21)〈B1y-〉=〈B1y-〉〈e12/x-(t)〉N-1g11+g21(〈e12/x-(t)〉/N-1)1-g22(〈e22/x-(t)〉/N-2)〈B1y-〉〈e22/x-(t)〉N-2g12,
where g11=〈B1/(-Ay)-1K1〉 and g12=〈B1/(-Ay)-1K2〉.If〈B1/y-〉=0, we obtain the disease-free equilibrium again. If not, the following equation holds:
(22)1=〈e12/x-(t)〉N-1×[g11+g21g121-g22(〈e22/x-(t)〉/N-2)〈e22/x-(t)〉N-2],
which can also be written as
(23)〈e22x-(t)〉=N-1N-2-g11N-2〈e12/x-(t)〉g12g21〈e12/x-(t)〉+g22N-1-g11g22〈e12/x-(t)〉.
This expression of 〈e22/x-(t)〉 in (23) is positive when one of these conditions is satisfied:
(24)N-1-g11〈e12x-(t)〉>0,N-1-g11〈e12x-(t)〉<0,g12g21〈e12x-(t)〉+g22N-1<g11g22〈e12x-(t)〉.Case  1. N-1-g11〈e12/x-(t)〉>0. This yields the condition
(25)〈e12x-(t)〉<N-1g11.Case  2. N-1-g11〈e12/x-(t)〉<0 and g12g21〈e12/x-(t)〉+g22N-1<g11g22〈e12/x-(t)〉. This yields the condition
(26)〈e12x-(t)〉>max⁡(N-1g11,g22N-1g11g22-g12g21)mmmmmlmlsinceg11g22-g12g21>0.
In sum, 〈e22/x-(t)〉 exists when 〈e12/x-(t)〉 satisfies condition (25) or condition (26). Now, using the first equation of system (15) gives
(27)Λ+Dxx-=∑i=12〈Biy-(t)〉〈ei2/x-(t)〉N-iei2.
Then,
(28)〈B1/y-〉N-1=Λ1-(μ1+a1)x-1+a2x-2x-1=〈Λ+Dxx-/e12〉〈e12/x-〉,〈B2/y-〉N-2=Λ2-(μ2+a2)x-2+a1x-1x-2=〈Λ+Dxx-/e22〉〈e22/x-〉.
Putting (28) into the first equation of (17) gives
(29)N-1=〈e12x-(t)〉×[g11〈Λ+Dxx-e12〉+g12〈Λ+Dxx-e22〉]×(〈Λ+Dxx-e12〉)-1.
Putting (28) into the second equation of (17) gives
(30)N-2=〈e22x-(t)〉×[g21〈Λ+Dxx-e12〉+g22〈Λ+Dxx-e22〉]×(〈Λ+Dxx-e22〉)-1.
Then, in the compact form,
(31)N-1〈e12/x-(t)〉=g11〈Λ+Dxx-/e12〉+g12〈Λ+Dxx-/e22〉〈Λ+Dxx-/e12〉,N-2〈e22/x-(t)〉=g21〈Λ+Dxx-/e12〉+g22〈Λ+Dxx-/e22〉〈Λ+Dxx-/e22〉.
From (1),
(32)N˙1+N˙2=Λ1+Λ2-μ1N1-μ2N2-d1I1-l1L1-d2I2-l2L2,
which can be written in semicompact form as
(33)N˙1+N˙2=Λ1+Λ2-μ1N1-μ2N2-〈B3y〉-〈B4y〉,
where B3=(0,…,0,d1,l1,0,…,0,0,0)T, B4=(0,…,0,0,0,0,…,0,d2,l2)T, and N=N1+N2.At equilibrium,(34)Λ1+Λ2=μ1N-1+μ2N-2+〈B3y-〉+〈B4y-〉.
Using (16) in (34),
(35)Λ1+Λ2=μ1N-1+μ2N-2+〈B1y-(t)〉〈e12/x-(t)〉N-1h31+〈B2y-(t)〉〈e22/x-(t)〉N-2h32+〈B1y-(t)〉〈e12/x-(t)〉N-1h41+〈B2y-(t)〉〈e22/x-(t)〉N-2h42,
where h31=〈B3/(-Ay)-1K1〉, h32=〈B3/(-Ay)-1K2〉, h41=〈B4/(-Ay)-1K1〉, and h42=〈B4/(-Ay)-1K2〉. Using (28) in (35) gives
(36)Λ1+Λ2=μ1N-1+μ2N-2+〈Λ+Dxx-e12〉h31+〈Λ+Dxx-e22〉h32+〈Λ+Dxx-e12〉h41+〈Λ+Dxx-e22〉h42.
Using (31) in (36) gives
(37)1=μ1Λ1+Λ2〈e12x-〉×g11〈Λ+Dxx-/e12〉+g12〈Λ+Dxx-/e22〉〈Λ+Dxx-/e12〉+μ2Λ1+Λ2〈e22x-〉×g21〈Λ+Dxx-/e12〉+g22〈Λ+Dxx-/e22〉〈Λ+Dxx-/e22〉+h31Λ1+Λ2〈Λ+Dxx-e12〉+h32Λ1+Λ2〈Λ+Dxx-e22〉+h41Λ1+Λ2〈Λ+Dxx-e12〉+h42Λ1+Λ2〈Λ+Dxx-e22〉.
Let us set
(38)F(x-)=μ1Λ1+Λ2〈e12x-〉×g11〈Λ+Dxx-/e12〉+g12〈Λ+Dxx-/e22〉〈Λ+Dxx-/e12〉+μ2Λ1+Λ2〈e22x-〉×g21〈Λ+Dxx-/e12〉+g22〈Λ+Dxx-/e22〉〈Λ+Dxx-/e22〉+h31Λ1+Λ2〈Λ+Dxx-e12〉+h32Λ1+Λ2×〈Λ+Dxx-e22〉+h41Λ1+Λ2〈Λ+Dxx-e12〉+h42Λ1+Λ2×〈Λ+Dxx-e22〉.One can observe that it is useful to analyze how many times the functionF(x-) intersects the line F(x-)=1 in the plane (x-,F(x-)). The monotony of the function F depends on the parameters of system (1). Then, the following lemma and theorem hold.Lemma 12.
The functionF satisfies the following properties.(1)
The limit whenx-→0 is finite and is given by
(39)lim⁡x-→0F(x-)=h31Λ1+Λ2〈Λe12〉+h32Λ1+Λ2〈Λe22〉+h41Λ1+Λ2〈Λe12〉+h42Λ1+Λ2〈Λe22〉=s0.(2)
The limit whenx-→+∞ is infinite: lim⁡x-→+∞F(x-)=+∞.Theorem 13.
These properties hold for system (1).(1)
There is no boundary equilibrium. This means that the disease cannot persist in one patch while disappearing in the other one.(2)
Ifs0≤R0≤1, there exists at least one solution of F(x-)=1, depending on the monotony of the function F, and therefore there exists at least one endemic equilibrium (x-,y-) for system (5) which coexists with the disease-free equilibrium. This situation corresponds to a backward bifurcation.(3)
IfR0≤s0≤1, there exists at least one solution of F(x-)=1, depending on the monotony of the function F, and therefore there exists at least one endemic equilibrium (x-,y-) for system (5) which coexists with the disease-free equilibrium. There is no backward bifurcation in this case since the endemic equilibrium no longer exists (it disappears as the disease-free equilibrium).(4)
IfR0≤1≤s0, at least one solution of F(x-)=1 can exist, depending on parameters of system (1), since the monotony of the function F also depends on these parameters.(5)
Ifs0≤1≤R0, at least one solution of F(x-)=1 can exist, depending on the monotony of the function F, and therefore there exists at least one endemic equilibrium (x-,y-) for system (5), which coexists with the disease-free equilibrium.(6)
If1≤R0≤s0, at least one solution of F(x-)=1 can exist, depending on parameters of system (1), since the monotony of the function F also depends on these parameters.(7)
If1≤s0≤R0, at least one solution of F(x-)=1 can exist, depending on parameters of system (1), since the monotony of the function F also depends on these parameters.Proof.
When an endemic equilibrium(x-,y-) for system (5) exists, x- is solution of F(x-)=1. N-1 and N-2 are given by (31) as
(40)N-1〈e12/x-(t)〉=g11〈Λ+Dxx-/e12〉+g12〈Λ+Dxx-/e22〉〈Λ+Dxx-/e12〉,N-2〈e22/x-(t)〉=g21〈Λ+Dxx-/e12〉+g22〈Λ+Dxx-/e22〉〈Λ+Dxx-/e22〉.〈B1/y-〉 and 〈B2/y-〉 are given by (16) as 〈Bi/y-〉/N-i=〈Λ+Dxx-/ei2〉/〈ei2/x-〉, i=1,2.
The expression ofy- is given by (16) as
(41)y-=〈B1y-〉〈e12/x-(t)〉N-1(-Ay)-1K1+〈B2y-〉〈e22/x-(t)〉N-2(-Ay)-1K2.
Since the function F is not strictly monotone, the number of solutions of F(x-)=1 depends on the parameters of system (1). Therefore, there can be more than one endemic equilibrium.Remark 14.
In the second case of the theorem, we observe that at least one endemic equilibrium coexists with the disease-free equilibrium. Then, there can be a backward bifurcation in this case.


### 3.6. Stability of Endemic Equilibrium

The stability of endemic equilibrium is always a big challenge in epidemiology. The problem is more difficult here since we have2n+6 equations. For multipatch and universal incidence law, results concerning the global stability of endemic equilibrium are limited. Next we will illustrate some results concerning our model.


## 3.1. Positivity of the Solutions

Since the variables considered here are nonnegative quantities, we have to be sure that their values are always nonnegative.Theorem 1.
The nonnegative orthantIR+2n+6 is positively invariant by (1). This means that every trajectory, which begins in the positive orthant, will stay inside.Proof.
System (4) can be written in the following form:
(5)dx(t)dt=Λ+Dxx-∑i=12〈Biy(t)〉〈ei2/x(t)〉Niei2,dy(t)dt=(∑i=12BiT〈ei2/x(t)〉NiKi+Ay)y(t).
Since x(t)≥0 and Ni(t)≥0, the matrix (∑i=12BiT(〈ei2/x(t)〉/Ni)Ki+Ay) is Metzler. It is well known that linear Metzler system lets the nonnegative orthant invariant [23]. Moreover dx(t)/dt≤Λ+Dxx and Dx is a Metzler matrix. These prove the positive invariance of the nonnegative orthant IR+2n+6 by (1).


## 3.2. Boundedness of the Trajectories

Lemma 2.
The simplex{(S1,E11,…,E1n,I1,L1,S2,E21,…,E2n,I2,L2)∈IR+2n+6/S1+E11+⋯E1n+I1 + L1+S2+E21+⋯E2n+I2+L2≤((Λ1+Λ2)/μ)+ε}, is a compact forward and absorbing set for (1).Proof.
From system (1), if N(t)=N1(t)+N2(t) denotes the entire population of (1), one has N(t)=S1+E11+⋯+E1n+I1+L1+S2+E21+⋯+E2n+I2+L2.
Then, the dynamics ofN(t) at any time t is given by
(6)dN(t)dt=Λ1+Λ2-μ1N1(t)-μ2N2(t)-d1I1-l1L1-d2I2-l2L2.
Thus, dN(t)/dt≤Λ1+Λ2-μN(t), where μ=min⁡(μ1,μ2). It follows that limSupt→+∞N(t)=(Λ1+Λ2)/μ.In the simplexΓε, (1) is mathematically well posed. The following lemma also holds.Lemma 3.
The simplexγε={(x,y)∈Γε/x≤x*} is a compact forward invariant set for (1).


## 3.3. Local Stability of the Disease-Free Equilibrium

Many epidemiological models have a threshold condition which can determine whether an infection will be eliminated from the population or become endemic. The basic reproduction number,R0, is defined as the average number of secondary infections produced by an infected individual in a completely susceptible population [16]. As discussed in [24, 25], R0 is a simply normalized bifurcation (transcritical) parameter for epidemiological models, such that R0 implies that the endemic steady state is stable (i.e., the infection persists), and R0 implies that the uninfected steady state is stable (i.e., the infection can be eliminated from the population).Equation (1) has a disease-free equilibrium given by K2=(S1*,0,…,0︸n+2,S2*,0,…,0︸n+2)T which always exists in the nonnegative orthant IR+2n+6. The explicit expressions of S1* and S2* are
(7)S1*=(μ2+a2)Λ1+a2Λ2μ1μ2+μ1a2+μ2a1,S2*=(μ1+a1)Λ2+a1Λ1μ1μ2+μ1a2+μ2a1.
The disease-free equilibrium can also be denoted by (x*,0) since it is the solution of equation Λ+Dxx*=0 of the compact system (5).Lemma 4.
Using the same method as in [17], the basic reproduction ratio of (1) is R0=ρ(K1B1T(-Ay)-1+K2B2T(-Ay)-1), where ρ denotes the spectral radius.Proof.
The expressions, which are coming from the other compartment, due to contamination, are given by the following matrix:(8)F=∂∂y(∑i=12〈Biy(t)〉〈ei2/x(t)〉NiKi)(x*,0)=∑i=12〈ei2/x*〉Ni*KiBiT.
The expressions, which are coming from the other compartment, due to reasons different from contamination, are given by V=-Ay. The next generation matrix, since 〈ei2/x*〉=Ni*, is
(9)FV-1=(∑i=12〈ei2/x*〉Ni*KiBiT)(-Ay)-1=(∑i=12KiBiT)(-Ay)-1.
We can observe that, since R0 is the largest eigenvalue of the next generation matrix,
(10)R0=ρ(FV-1)=ρ(K1B1T(-Ay)-1+K2B2T(-Ay)-1),
where ρ denotes the spectral radius.Consequently, from Theorem 2 of [22], the following lemma holds.Lemma 5.
The disease-free equilibrium of (1) is locally and asymptotically stable whenever R0≤1 and unstable if R0>1.Remark 6.
This lemma shows that ifR0≤1, a small flow of infectious individuals will not generate large outbreaks of the disease. To control the disease independently of the initial total number of infectious individuals, a global asymptotic stability property has to be established for the DFE when R0≤1.


## 3.4. Global Stability of the Disease-Free Equilibrium

The following result helps to determine the stability and is related to LaSalle’s principle [13]. Consider the differential equation
(11)dx(t)dt=X(x(t)),X(x0)=0,
where X is a C1 function defined on an open set of IRn containing the closure Ω- of a positively invariant set Ω such that the equilibrium x0 is in Ω-. The following lemma holds.Lemma 7.
We assume that system (11) is point dissipative [1] on Ω. In other words there exists a compact set K⊂Ω such that, for any y∈Ω, there exists a time t(y0)>0 such that, for any time t>t(y0), the trajectory with initial condition y0 is in the interior of K. If there exists a C1 function V≥0 defined on Ω-, then(1)
dV(t)/dt=〈δV(x)/X(x)〉≤0, for all x∈Ω;(2)
the greatest invariant setL contained in E={x∈Ω-,dV(t)/dt=0} is contained in a positively invariant set on which the restriction of system (11) is globally and asymptotically stable on L at x0.Then, x0 is a globally and asymptotically stable equilibrium of system (11) on Ω.Theorem 8.
WhenR0≤1 (this implies R01≤1 and R02≤1, where R0i is the basic reproduction number for the patch i), the disease-free equilibrium (DFE), when it is unique, is globally and asymptotically stable in Γε, since it is the unique equilibrium. This implies the global asymptotic stability of the DFE on the nonnegative orthant IR+2n+6; that is, the disease dies out in both two populations.Remark 9.
At least one endemic equilibrium can exist and coexist with the disease-free equilibrium whenR0≤1. In this case the DFE cannot be globally asymptotically stable.Proof.
Following the same method as in [23], system (4) can be written in a pseudotriangular form as
(12)dx(t)dt=A1(x,y)(x-x*)+A12(x,y)y=f(x,y),dy(t)dt=A2(x,y)y=g(x,y),
where f(x*,0)=Λ+Dxx*=0, g(x*,0)=0, A1(x,y)=Dx∈IR2x2, A12(x,y) = -∑i=12(〈ei2/x(t)〉/Ni)ei2BiT∈IR2x(2n+4), and A2(x,y)=∑i=12(〈ei2/x(t)〉/Ni)KiBiT+Ay∈IR(2n+4)x(2n+4). The Jacobian matrix of system (4) at disease-free equilibrium (x*,0) is A2(x,y)=∑i=12(〈ei2/x(t)〉/Ni)KiBiT+Ay∈IR(2n+4)x(2n+4). The Jacobian matrix of system (4) at disease-free equilibrium (x*,0) is J=(A1(x*,0)A12(x*,0)0A2(x*,0))=(Dx-∑i=12ei2BiT0∑i=12KiBiT+Ay). The matrix Dx is clearly a Metzler stable matrix and, using a result in [8], there exists a vector u=(u1,u2)>0 such that Dxu≤0. The matrix A-2=∑i=12KiBiT+Ay is a Metzler and irreducible matrix, which is stable if R0≤1. With this condition satisfied and using the same previous result in [22], the stability modulus α(A-2) of the matrix is in the spectrum of A-2 (Perron-Frobenius theorem, proof of Lemma 2) and there exists a vector c>0∈IR2n+4 such that cTA-2=α(A-2)cT, with α(A-2)≤0. Let us show the global stability of the DFE.
System (12) is in pseudotriangular form since the matrix A2(x,y) depends on (x,y). There are many results for the stability of triangular systems [21–25]. Using LaSalle’s principle [13] one can obtain attractivity of equilibrium. But, for nonlinear systems, attractivity does not generally implies stability. We can use now the result of Lemma 7. Let us consider the following candidate Lyapunov function:
(13)VDFE(x(t),y(t))=〈(0,c)(x(t),y(t))T〉.
The function VDFE(x(t),y(t)) is positive and it is time derivative along the trajectories of (4) which gives
(14)dV(x(t),y(t))dt=(∑i=12〈Biy(t)〉〈ei2/x(t)〉NiKi+Ayy)c≤cT(∑i=12KiBiT+Ay)y≤cTα(A-2)y.
Then, R0≤1 ensures that dVDFE(x(t),y(t))/dt≤0 for all y≥0 and that dVDFE(x(t),y(t))/dt=0 holds when (x(t),y(t))=(x*,0). This proves the global asymptotic stability on IR+2n+6 [14]. This achieves the proof that the DFE is globally and asymptotically stable.


## 3.5. Existence of Endemic Equilibrium

Definition 10.
An equilibrium for a multipatch model as (1) is called endemic equilibrium when the two populations coexist (the density of each compartment is different from zero) at this equilibrium.Lemma 11.
An endemic equilibrium(x-,y-) of (1) can be determined, using an equation F(x-)=1.Proof.
An endemic equilibrium(x-,y-) of (1) is obtained by setting the right-hand side of (4) which equals zero, giving
(15)Λ+Dxx-=∑i=12〈Biy-(t)〉〈ei2/x-(t)〉N-iei2,∑i=12〈Biy-(t)〉〈ei2/x-(t)〉N-iKi+Ayy-=0.
Multiplying the second equation of (15) by (-Ay)-1 gives
(16)y-=〈B1y-(t)〉〈e12/x-(t)〉N-1(-Ay)-1K1+〈B2y-(t)〉〈e22/x-(t)〉N-2(-Ay)-1K2.
From (16),
(17)〈B1y-〉=g11〈B1y-〉〈e12/x-(t)〉N-1+g12〈B2y-〉〈e22/x-(t)〉N-2,〈B2y-〉=g21〈B1y-〉〈e12/x-(t)〉N-1+g22〈B2y-〉〈e22/x-(t)〉N-2,
where
(18)g11=〈B1(-Ay)-1K1〉,g12=〈B1(-Ay)-1K2〉,g21=〈B2(-Ay)-1K1〉,g22=〈B2(-Ay)-1K2〉.
From the second equation of (17),
(19)〈B2y-〉=g21(〈e12/x-(t)〉/N-1)1-g22(〈e22/x-(t)〉/N-2)〈B1y-〉.
Using the expression of 〈B2/y-〉 in (16) gives
(20)y-=〈B1y-〉〈e12/x-(t)〉N-1(-Ay)-1K1+g21(〈e12/x-(t)〉/N-1)1-g22(〈e22/x-(t)〉/N-2)〈B1y-〉×〈e22/x-(t)〉N-2(-Ay)-1K2.
Then
(21)〈B1y-〉=〈B1y-〉〈e12/x-(t)〉N-1g11+g21(〈e12/x-(t)〉/N-1)1-g22(〈e22/x-(t)〉/N-2)〈B1y-〉〈e22/x-(t)〉N-2g12,
where g11=〈B1/(-Ay)-1K1〉 and g12=〈B1/(-Ay)-1K2〉.If〈B1/y-〉=0, we obtain the disease-free equilibrium again. If not, the following equation holds:
(22)1=〈e12/x-(t)〉N-1×[g11+g21g121-g22(〈e22/x-(t)〉/N-2)〈e22/x-(t)〉N-2],
which can also be written as
(23)〈e22x-(t)〉=N-1N-2-g11N-2〈e12/x-(t)〉g12g21〈e12/x-(t)〉+g22N-1-g11g22〈e12/x-(t)〉.
This expression of 〈e22/x-(t)〉 in (23) is positive when one of these conditions is satisfied:
(24)N-1-g11〈e12x-(t)〉>0,N-1-g11〈e12x-(t)〉<0,g12g21〈e12x-(t)〉+g22N-1<g11g22〈e12x-(t)〉.Case  1. N-1-g11〈e12/x-(t)〉>0. This yields the condition
(25)〈e12x-(t)〉<N-1g11.Case  2. N-1-g11〈e12/x-(t)〉<0 and g12g21〈e12/x-(t)〉+g22N-1<g11g22〈e12/x-(t)〉. This yields the condition
(26)〈e12x-(t)〉>max⁡(N-1g11,g22N-1g11g22-g12g21)mmmmmlmlsinceg11g22-g12g21>0.
In sum, 〈e22/x-(t)〉 exists when 〈e12/x-(t)〉 satisfies condition (25) or condition (26). Now, using the first equation of system (15) gives
(27)Λ+Dxx-=∑i=12〈Biy-(t)〉〈ei2/x-(t)〉N-iei2.
Then,
(28)〈B1/y-〉N-1=Λ1-(μ1+a1)x-1+a2x-2x-1=〈Λ+Dxx-/e12〉〈e12/x-〉,〈B2/y-〉N-2=Λ2-(μ2+a2)x-2+a1x-1x-2=〈Λ+Dxx-/e22〉〈e22/x-〉.
Putting (28) into the first equation of (17) gives
(29)N-1=〈e12x-(t)〉×[g11〈Λ+Dxx-e12〉+g12〈Λ+Dxx-e22〉]×(〈Λ+Dxx-e12〉)-1.
Putting (28) into the second equation of (17) gives
(30)N-2=〈e22x-(t)〉×[g21〈Λ+Dxx-e12〉+g22〈Λ+Dxx-e22〉]×(〈Λ+Dxx-e22〉)-1.
Then, in the compact form,
(31)N-1〈e12/x-(t)〉=g11〈Λ+Dxx-/e12〉+g12〈Λ+Dxx-/e22〉〈Λ+Dxx-/e12〉,N-2〈e22/x-(t)〉=g21〈Λ+Dxx-/e12〉+g22〈Λ+Dxx-/e22〉〈Λ+Dxx-/e22〉.
From (1),
(32)N˙1+N˙2=Λ1+Λ2-μ1N1-μ2N2-d1I1-l1L1-d2I2-l2L2,
which can be written in semicompact form as
(33)N˙1+N˙2=Λ1+Λ2-μ1N1-μ2N2-〈B3y〉-〈B4y〉,
where B3=(0,…,0,d1,l1,0,…,0,0,0)T, B4=(0,…,0,0,0,0,…,0,d2,l2)T, and N=N1+N2.At equilibrium,(34)Λ1+Λ2=μ1N-1+μ2N-2+〈B3y-〉+〈B4y-〉.
Using (16) in (34),
(35)Λ1+Λ2=μ1N-1+μ2N-2+〈B1y-(t)〉〈e12/x-(t)〉N-1h31+〈B2y-(t)〉〈e22/x-(t)〉N-2h32+〈B1y-(t)〉〈e12/x-(t)〉N-1h41+〈B2y-(t)〉〈e22/x-(t)〉N-2h42,
where h31=〈B3/(-Ay)-1K1〉, h32=〈B3/(-Ay)-1K2〉, h41=〈B4/(-Ay)-1K1〉, and h42=〈B4/(-Ay)-1K2〉. Using (28) in (35) gives
(36)Λ1+Λ2=μ1N-1+μ2N-2+〈Λ+Dxx-e12〉h31+〈Λ+Dxx-e22〉h32+〈Λ+Dxx-e12〉h41+〈Λ+Dxx-e22〉h42.
Using (31) in (36) gives
(37)1=μ1Λ1+Λ2〈e12x-〉×g11〈Λ+Dxx-/e12〉+g12〈Λ+Dxx-/e22〉〈Λ+Dxx-/e12〉+μ2Λ1+Λ2〈e22x-〉×g21〈Λ+Dxx-/e12〉+g22〈Λ+Dxx-/e22〉〈Λ+Dxx-/e22〉+h31Λ1+Λ2〈Λ+Dxx-e12〉+h32Λ1+Λ2〈Λ+Dxx-e22〉+h41Λ1+Λ2〈Λ+Dxx-e12〉+h42Λ1+Λ2〈Λ+Dxx-e22〉.
Let us set
(38)F(x-)=μ1Λ1+Λ2〈e12x-〉×g11〈Λ+Dxx-/e12〉+g12〈Λ+Dxx-/e22〉〈Λ+Dxx-/e12〉+μ2Λ1+Λ2〈e22x-〉×g21〈Λ+Dxx-/e12〉+g22〈Λ+Dxx-/e22〉〈Λ+Dxx-/e22〉+h31Λ1+Λ2〈Λ+Dxx-e12〉+h32Λ1+Λ2×〈Λ+Dxx-e22〉+h41Λ1+Λ2〈Λ+Dxx-e12〉+h42Λ1+Λ2×〈Λ+Dxx-e22〉.One can observe that it is useful to analyze how many times the functionF(x-) intersects the line F(x-)=1 in the plane (x-,F(x-)). The monotony of the function F depends on the parameters of system (1). Then, the following lemma and theorem hold.Lemma 12.
The functionF satisfies the following properties.(1)
The limit whenx-→0 is finite and is given by
(39)lim⁡x-→0F(x-)=h31Λ1+Λ2〈Λe12〉+h32Λ1+Λ2〈Λe22〉+h41Λ1+Λ2〈Λe12〉+h42Λ1+Λ2〈Λe22〉=s0.(2)
The limit whenx-→+∞ is infinite: lim⁡x-→+∞F(x-)=+∞.Theorem 13.
These properties hold for system (1).(1)
There is no boundary equilibrium. This means that the disease cannot persist in one patch while disappearing in the other one.(2)
Ifs0≤R0≤1, there exists at least one solution of F(x-)=1, depending on the monotony of the function F, and therefore there exists at least one endemic equilibrium (x-,y-) for system (5) which coexists with the disease-free equilibrium. This situation corresponds to a backward bifurcation.(3)
IfR0≤s0≤1, there exists at least one solution of F(x-)=1, depending on the monotony of the function F, and therefore there exists at least one endemic equilibrium (x-,y-) for system (5) which coexists with the disease-free equilibrium. There is no backward bifurcation in this case since the endemic equilibrium no longer exists (it disappears as the disease-free equilibrium).(4)
IfR0≤1≤s0, at least one solution of F(x-)=1 can exist, depending on parameters of system (1), since the monotony of the function F also depends on these parameters.(5)
Ifs0≤1≤R0, at least one solution of F(x-)=1 can exist, depending on the monotony of the function F, and therefore there exists at least one endemic equilibrium (x-,y-) for system (5), which coexists with the disease-free equilibrium.(6)
If1≤R0≤s0, at least one solution of F(x-)=1 can exist, depending on parameters of system (1), since the monotony of the function F also depends on these parameters.(7)
If1≤s0≤R0, at least one solution of F(x-)=1 can exist, depending on parameters of system (1), since the monotony of the function F also depends on these parameters.Proof.
When an endemic equilibrium(x-,y-) for system (5) exists, x- is solution of F(x-)=1. N-1 and N-2 are given by (31) as
(40)N-1〈e12/x-(t)〉=g11〈Λ+Dxx-/e12〉+g12〈Λ+Dxx-/e22〉〈Λ+Dxx-/e12〉,N-2〈e22/x-(t)〉=g21〈Λ+Dxx-/e12〉+g22〈Λ+Dxx-/e22〉〈Λ+Dxx-/e22〉.〈B1/y-〉 and 〈B2/y-〉 are given by (16) as 〈Bi/y-〉/N-i=〈Λ+Dxx-/ei2〉/〈ei2/x-〉, i=1,2.
The expression ofy- is given by (16) as
(41)y-=〈B1y-〉〈e12/x-(t)〉N-1(-Ay)-1K1+〈B2y-〉〈e22/x-(t)〉N-2(-Ay)-1K2.
Since the function F is not strictly monotone, the number of solutions of F(x-)=1 depends on the parameters of system (1). Therefore, there can be more than one endemic equilibrium.Remark 14.
In the second case of the theorem, we observe that at least one endemic equilibrium coexists with the disease-free equilibrium. Then, there can be a backward bifurcation in this case.


## 3.6. Stability of Endemic Equilibrium

The stability of endemic equilibrium is always a big challenge in epidemiology. The problem is more difficult here since we have2n+6 equations. For multipatch and universal incidence law, results concerning the global stability of endemic equilibrium are limited. Next we will illustrate some results concerning our model.


## 4. Numerical Simulations

System (1) is simulated with parameter values presented in Table 1 using real data of the cities of Yaounde and Bafia [21].Table 1
Numerical values for the parameters of the model.
Parameters
Description
Estimated value/range
Reference
Λ
1
Recruitment rate into theS1 class
100/yr
Assumed
Λ
2
Recruitment rate into theS2 class
110/yr
Assumed
β
1
Transmission coefficient in patch 1
Variable
Assumed
β
2
Transmission coefficient in patch 2
Variable
Assumed
μ
1
Per capita death rate in patch 1
0.019896/yr
Estimated
μ
2
Per capita death rate in patch 2
0.019897/yr
Estimated
k
1
i
Progression rate fromE1i to next class in patch 1
0.00013/yr
Assumed
k
2
i
Progression rate fromE2i to next class in patch 2
0.00023/yr
Assumed
r
1
i
Effective chemoprophylaxis in theE1i class
0/yr
Estimated
r
2
i
Effective chemoprophylaxis in theE2i class
0/yr
Estimated
δ
1
Lost sight individuals still able to transmit disease in patch 1
0.46/yr
Estimated
δ
2
Lost sight individuals still able to transmit disease in patch 2
0.48/yr
Estimated
γ
1
Effective therapy inI1 class
0.8182/yr
Estimated
γ
2
Effective therapy inI2 class
0.8183/yr
Assumed
α
1
Return rate in care centre for patch 1
0.23/yr
Estimated
α
2
Return rate in care centre for patch 2
0.23/yr
Estimated
v
1
Infectious individuals in care centre of patch 1 who can stay in care centre
0.59/yr
Estimated
v
2
Infectious individuals in care centre of patch 2 who can stay in care centre
0.5903/yr
Estimated
s
1
Infectious individuals in care centre of patch 1 who disappear
0.019/yr
Estimated
s
2
Infectious individuals in care centre of patch 1 who disappear
0.01905/yr
Estimated
a
1
Migration rate of susceptible individuals  from patch 1 to patch 2
0.07/yr
Assumed
a
2
Migration rate of susceptible individuals  from patch 2 to patch 1
0.0701/yr
Assumed
b
1
i
Migration rate from latent classE1i to latent class E2i
0.05/yr
Assumed
b
2
i
Migration rate from latent classE2i to latent class E1i
0.0503/yr
Assumed
c
1
Migration rate of infectious individuals from patch 1 to patch 2
0.02/yr
Assumed
c
2
Migration rate of infectious individuals from patch 2 to patch 1
0.0201/yr
Assumed
d
1
Additional death rate inI1
0.0575/yr
Assumed
d
2
Additional death rate inI2
0.05751/yr
Assumed
l
1
Additional death rate inL1
0.0565/yr
Assumed
l
2
Additional death rate inL2
0.05651/yr
Assumed
p
1
Fast progression to theI1 class
0.015/yr
Assumed
p
2
Fast progression to theI2 class
0.016
AssumedNext we use the values in Table1 to generate several curves (Figures 2, 3, 4, 5, 6, 7, 8, and 9) in order to illustrate our theoretical results.Figure 2
Trajectories of susceptible individuals of (1) for different initial conditions when β1=0.002 and β2=0.003 (such that R01=0.3691 and R02=0.4980). We can observe the global stability of DFE when R0<1. All other parameters are as in Table 1.Figure 3
Trajectories of latently infected individuals of (1) for different initial conditions when β1=0.002 and β2=0.003 (such that R01=0.3691 and R02=0.4980). The latent classes disappear at the end and only susceptible classes persist when R0<1. All other parameters are as in Table 1.Trajectories of infectious individuals in care center (a) and lost sight individuals (b) of (1) for different initial conditions when β1=0.002 and β2=0.003 (such that R01=0.3691 and R02=0.4980). The infectious and lost sight classes disappear at the end and only susceptible classes persist when R0<1. All other parameters are as in Table 1.
(a)
(b)Figure 5
Trajectories of susceptible individuals of (1) for different initial conditions when β1=0.002 and β2=0.003 (such that R01=0.3691 and R02=0.4980). The DFE coexists with endemic equilibrium and is locally stable when R0<1. All other parameters are as in Table 1.Figure 6
Trajectories of latently infected individuals of (1) for different initial conditions when β1=0.822 and β2=0.943 (such that R01=2.3601 and R02=2.4910). The endemic equilibrium is locally and asymptotically stable when R0>1. All other parameters are as in Table 1.Figure 7
Trajectories of latently infected individuals of (1) for different initial conditions when β1=2.082 and β2=2.173 (such that R01=3.1492 and R02=4.1085). The endemic equilibrium is locally and asymptotically stable when R0>1. In this case, the DFE is unstable. All other parameters are as in Table 1.(a, b) Trajectories of infectious individuals in care center of (1) for different initial conditions. The endemic equilibrium exists and is stable when R0>1; the DFE is unstable in this case.
(a)
(b)(a, b) Trajectories of lost sight individuals of (1) for different initial conditions. The endemic equilibrium exists and is stable when R0>1; the DFE is unstable in this case.
(a)
(b)


## 5. Conclusion

The long-term dynamics of our system has been completely investigated. The model exhibits rich dynamics, depending on the values of the bifurcation parametersR0i, i=1,2, and R0. The influence of parameter R0i is significant on the spread of tuberculosis since it quantifies the intensity of pathogens transmission. The stability of equilibria depends on these parameters. We have transcritical bifurcation parameters and backward bifurcation. When the basic reproduction number is less than unity, tuberculosis can be controlled in each population if the DFE is the unique equilibrium. It can be more difficult if the DFE coexists with at least one endemic equilibrium (this is the situation of backward bifurcation). In this case, the disease-free equilibrium can be locally and asymptotically stable, as well as the endemic equilibrium. When the basic reproduction number R0 is greater than unity, tuberculosis is endemic and can be difficult to control in the population. The disease in our model cannot persist in one population while disappearing in the other one.



---

*Source: 263780-2014-07-20.xml*