- Email: [email protected]
Modeling and analysis of Tuberculosis (TB) in Khyber Pakhtunkhwa, Pakistan
Accepted Manuscript Modeling and analysis of Tuberculosis (TB) in Khyber Pakhtunkhwa, Pakistan Saif Ullah, Muhammad Altaf Khan, Muhammad Farooq, Taza Gul
PII: DOI: Reference:
S0378-4754(19)30108-9 https://doi.org/10.1016/j.matcom.2019.03.012 MATCOM 4768
To appear in:
Mathematics and Computers in Simulation
Received date : 18 May 2018 Revised date : 27 January 2019 Accepted date : 26 March 2019 Please cite this article as: S. Ullah, M.A. Khan, M. Farooq et al., Modeling and analysis of Tuberculosis (TB) in Khyber Pakhtunkhwa, Pakistan, Mathematics and Computers in Simulation (2019), https://doi.org/10.1016/j.matcom.2019.03.012 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Manuscript Click here to view linked References
Modeling and analysis of Tuberculosis (TB) in Khyber Pakhtunkhwa, Pakistan January 28, 2019 Saif Ullah† , Muhammad Altaf Khan†† ,1 Muhammad Farooq † ,Taza Gul
††
† Department
of Mathematics, University of Peshawar, Peshawar, Khyber Pakhtunkhwa, Pakistan. †† Department of Mathematics, City University of Science and Information Technology, Peshawar, Khyber Pakhtunkhwa, Pakistan. Abstract
1
In this paper, we present a mathematical model on the dynamics of Tuberculosis (TB) in Khyber Pakhtunkhwa, province of Pakistan. The model is parameterized using TB infected cases reported from the year 2002 to 2017 in Khyber Pakhtunkhwa. The basic reproduction number R0 for the given period is estimated as R0 ≈ 1.3419. Mathematical analysis of the model in terms of stability analysis is derived and discussed in detail. The model at the disease free case is locally as well as globally asymptotically stable when the basic reproduction number R0 < 1. Further, when the basic reproduction number R0 > 1, then the model is stable locally as well as globally. Moreover, the sensitivity analysis of the model parameters are performed and their corresponding graphical results are presented. Finally, the numerical results for the sensitive parameters are plotted and their effects are shown on the disease elimination.
2 3 4 5 6 7 8 9 10 11 12 13
14 15
Key words: Tuberculosis (TB) Model, Parameter estimation, Lyapunov function, Stability.
16
1
17 18 19 20 21 22
Introduction
Tuberculosis (TB) is a bacterial infectious disease caused by bacillus Mycobacterium tuberculosis (MTB) and is one of the top 10 causes of death worldwide. TB infection mainly affects the lungs (known as pulmonary TB). It may also invade other organs including brain, spine, kidneys, central nervous system or lymphatic system (known as extrapulmonary TB). TB is also known as air-born infection because it is transmitted from TB active people thorough air when they cough, spit, sneeze or speak. The main symptoms of 1
Corresponding author. Email:[email protected]
1
23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67
TB infected person are chronic cough, blood-containing sputum, night sweats, fever and weight loss. According to World Health Organization (WHO) recent reports, abut (1/4)th of the whole world’s population are exposed to TB infection (which means the people infected by TB but not yet propagate the infection). More than 10.4 million people were infected with TB and about 1.7 million expire due to this infection in 2016. The person once infected with TB have 5-15 percent lifetime risk of falling ill with TB [1]. Mostly, the TB related deaths occur in the middle income countries such as India, Indonesia, China, Philippines, Pakistan, Nigeria, and South Africa. These countries accounted more than 60% of the whole TB burden in the world [1]. Mathematical modeling is one of the important tools to explore the dynamics of a disease and provide useful techniques to control it. In the last few decades, several biologists and mathematicians have been developed different epidemic models for the transmission dynamic of TB in different areas around the globe. Waaler et al. [17] introduced the first mathematical model of TB in 1962, by dividing the whole population into three classes. The transmission model of TB infection depending on the proportion of its prevalence is develop in [13]. A two strain, an age structure and time delay TB transmission models were proposed by Castillo et al. [2, 5]. Yang et al. [21] analyzed the global stability of TB models with incomplete treatment. A similar stability analysis of TB mathematical model was discussed by Liu et al. in [9]. A fractional TB model with Caputo derivative is investigated in [15]. The most powerful method to show the global dynamics of epidemic models is that of the Lyapunov which can be found in recent papers such as [19, 20]. Most of the models are based on the qualitative behavior of TB infection. However there are few articles presenting the quantitative analysis of TB data. In [10], Liu and co-authors investigated the TB dynamics with seasonality and use the epidemiological data to estimate the model parameters. Trauer et al. [14] construct a realistic model to study the dynamics of TB in different regions of Asia-pacific. Zhang et al. [22] proposed the TB mathematical model with hospitalized and non hospitalized infective population and use the TB data of China to simulate his model. Gerberry examined the TB transmission model and explored the practical aspects of backward bifurcation [6]. A TB model with relapse and reinfection is analyzed by Robert [18]. Recently, mathematical model with optimal control strategies for TB dynamics in the Philippines has been studied by Kim et al. [7]. In Pakistan, TB is a leading public health problem and one of the major causes of morbidity and mortality. Pakistan is one of the highest countries with TB infection and ranked fifth among the twenty two countries having high burden of this infection. Annually more than 0.5 million new TB cases including 15000 children are reported and more than seventy thousand Pakistanis died because of TB infection each year. Due to high global prevalence of multidrug-resistant tuberculosis (MDR-TB), Pakistan ranked fourth in the world [12, 4, 1]. TB is also a major health problem and leading cause of death in the province of Pakistan named Khyber Pakhtunkhwa. According to national TB control program (NTP) about 462920 TB infected cases are registered and treated from 2002 till 2017 in Khyber Pakhtunkhwa [11]. In this paper we construct a mathematical model to investigate the dynamics of TB infection in Khyber Pakhtunkhwa. We use the new TB infected cases from 2002 to 2017 obtained from the NTP [11], to estimate the biological parameters of the proposed model. The structure of the paper is organized as follows: The complete analysis of the model 2
75
and the parameters estimation from availed data are presented in the next section “Methods and materials”. The fundamental properties of the model are presented in section 3. The stability analysis of the disease-free case is given in section 4. The existence and stability analysis of endemic equilibrium point is provided in section 5. Also, in this section the analysis of bifurcation is discussed. In section 6, the sensitivity analysis is presented. Numerical simulations for the model based on fitted values of parameters are obtained to illustrate the theoretical results in section 7. Finally, the conclusion is summarized in section 8.
76
2
77
2.1
68 69 70 71 72 73 74
78 79 80 81 82 83
Methods and materials Model formulation
In this section we present the proposed mathematical model of TB infection. To develop the model, total human population is divided into five epidemiological sub-compartments denoted by susceptible S(t), Exposed L(t), TB active I(t), under treatment T (t), and recovered individuals R(t). The transmission dynamics of TB infection is given by the following system of nonlinear differential equations: dS = Λ− dt dL = dt
βSI N
βSI N
− µS,
− (µ + ǫ)L + (1 − η)δT,
dI = ǫL + ηδT − (µ + γ + σ1 )I, dt
(1)
dT = γI − (µ + δ + σ2 + α)T, dt dR = αT − µR. dt 84 85 86 87 88 89 90 91 92 93 94 95
Where N(t) = S(t) + L(t) + I(t) + T (t) + R(t). Recruitment rate is Λ and the natural death rate is µ. The number of new infection per unit time due to contact with TB infective individuals is denoted by parameter β. The progression rate from class L to I is ǫ. The infective individuals are treated at the rate γ. The treated individuals leaves the class T at a rate δ. The term (1 − η)δT , where the parameter η reflects the treatment failure, further by considering η = 0, will means that all the individuals under treatment become latent, and when η = 1, means the failure of treatment and all the treated individuals will still be infectious. The treated individuals enter to either latent class due to the remainder of Mycobacterium tuberculosis or infective class I due to the failure of treatment at the rate δ. The disease mortality rates in the infective class I and treated class T are σ1 and σ2 respectively. The parameter α represents recovery rate of the treated individuals. All these parameters and their flow rate to each class is shown in Figure 1 while the model 3
97
variables are given in Table 1. Since the last equation of (1) is independent on the rest, so we ignore it and consider the first four equations of the model.
98
2.2
96
99 100 101 102 103 104 105 106 107 108 109 110 111
Estimation of model parameters
In this subsection we estimate the model parameters based on the actual TB incidence data taken from national TB control program, Pakistan [11]. In order to parameterize TB model (1), we obtained some of the parameter values from literature (see Table 2), others were estimated or fitted from the data. For instance, the demographic parameter, µ, is estimated as µ = 1/67.7 per year, where 67.7 years is the average lifespan in Pakistan [1]. The other demographic parameter, Λ, is then estimated as follows. Since the total population of Khyber Pakhtunkhwa province as at 2017 was 30,523,371 [4], we assumed that Λ/µ, which is the limiting total human population in the absence of the disease, is 30, 523, 371, so that Λ = 450, 862.20088626 per year. The rest of model parameters are estimated using least-square fitting technique which provides a better fit of the model solution to the real data as shown in Fig. 2. The best fitted parameter values are obtained by minimizing the error between the actual TB incidence data and the solution of the proposed model (1). The objective function used in the parameter estimation is as follows ˆ = argmin Θ
n X j=1
112 113 114 115 116 117 118
2
(Itj − I tj ) ,
where I tj denotes the actual cumulative TB infected case and Itj are the corresponding model solution at time tj , while n is the number of available actual data points. In Fig. 2, the actual incidence data is shown by red circles while the model fitted curve is shown by blue solid line. The associated state variables and parameters of the model (1) are tabulated in Table:1 and Table:2 respectively. Consequently, using these parameter estimates, we obtained, in this study, that the value of R0 for the 2002-2017 TB cases in Khyber Pakhtunkhwa Pakistan is R0 ≈ 1.3419. Variable Description S(t) L(t) I(t) T (t) R(t)
The number of individuals which attract the disease at any time t The number of individuals which is not yet infectious The number of TB active individuals t The number of treated individuals at any time t Recovered individuals at any time t Table 1: Description of state variables of the model (1).
4
Parameter Λ β α γ µ σ1 σ2 δ η ǫ
Description Recruitment rate Disease contact rate Progression rate from T to R Transmission rate from I to T Natural mortality rate Disease related death rate of infected individuals Disease related death rate in T Rate at which treated individuals leave the T Treatment failure rate Rate of moving from L to I
Baseline value 450,862.20088626 0.5433 0.3968 0.2873 1/67.7 0.2202 0.0550 1.1996 0.1500 0.2007
Reference Estimated Fitted Fitted Fitted [12] Fitted Fitted Fitted Fitted Fitted
Table 2: Numerical values for the parameters of the TB model (1) estimated from Khyber Pakhtunkhwa population and TB incidence data.
(µ + σ1 )I I µS Λ
S
ǫL
µL βSI N
ηδT
L
γI
(1 − η)δT
R
µR
αT T
(µ + σ2 )T
Figure 1: Flow chart.
119
3
Analysis of the model
121
In this section we prove the basic properties of the proposed TB model (1). We proceed as follow:
122
3.1
120
123 124
Positivity and boundedness of the solution
To show that the TB infection model (1) is epidemiologically meaningful, it is necessary to show that all its state variables are non-negative for all time. In other words, solutions of 5
5
×10 4 Data Model Fit
Tuberculosis Incidence Data
4.5 4 3.5 3 2.5 2 1.5 1 0.5 2002
2004
2006
2008
2010
2012
2014
2016
2018
Years
Figure 2: Data fitting (blue solid curve) of the number of TB cases in Khyber Pakhtunkhwa, using model (1). 125 126
127 128
the model system (1) with non-negative initial data will remain non-negative for all time t > 0. Lemma 1. For all t > 0 and the initial data P (0) ≥ 0, where P (t) = (S(t), L(t), I(t), T (t)), the solution of TB model (1) are non-negative if they exist. Furthermore, lim supN(t) ≤
t→∞
129 130
Λ . µ
Proof. Let χ(t) = βI and t1 = sup{t > 0 : P (t) > 0 ∈ [0, t]}, then the from first equation N of system (1) we have dS = Λ − χ(t)S(t) − µS, dt = Λ − (χ(t) + µ)S(t).
131 132
(2)
The assumption that the solution of the system (1) exists in the interval I, the solution of (2) given as below Z t1 Z t1 S(t1 ) = S(0)exp − µt1 + χ(x)dx + exp − µt1 + χ1 (x)dx 0 0 Z t1 Z p × Λ exp µp + χ(y)dy dp > 0. 0
0
6
133 134 135
Similarly, it can be shown that P > 0 for all for all time t > 0. For the second part of the proof, note that 0 < S(0) ≤ N(t), 0 < L(0) ≤ N(t), 0 < I(0) ≤ N(t), 0 < T (0) ≤ N(t). Adding components of the model (1) gives dN(t) = Λ − µN(t) − σ1 I − σ2 T ≤ Λ − µN(t). dt
136
Hence lim supN(t) ≤
t→∞
Λ . µ
137
138
139 140
141
3.2
Invariant regions
The TB model (1) will be studied in a biologically feasible region as given below. Consider the feasible region Ω ⊂ R4+ , with
Ω= 142 143
144
(S(t), L(t), I(t), T (t)) ∈
R4+
Λ . : N(t) ≤ µ
Lemma 2. The region Ω ⊂ R4+ , is positively invariant for the model (1) with non-negative initial conditions in R4+ . Proof. It follows from summing equations of the model (1) that dN(t) = Λ − µN(t) − σ1 I − σ2 L ≤ Λ − µN(t), dt
145 146 147 148
149
hence dNdt(t) ≤ 0, if N(0) ≥ Λµ . Thus N(t) ≤ N(0)eµt + Λµ (1 − eµt ). In particular, N(t) ≤ N(0)eµt . Thus, the region Ω is positively invariant. Furthermore, if N(0) > Λµ , then either the solutions enters Ω in finite time, or N(t) approaches Λµ asymptotically. Hence, the region Ω attracts all solutions in R4+ .
4
Stability analysis of disease free equilibrium
The model (1) has a disease free equilibrium (DFE), Z0 , given by Z0 = (S 0 , 0, 0, 0), 150
where S0 =
151
Λ . µ
Let denote k1 = (µ + ǫ), k2 = (µ + γ + σ1 ), k3 = (µ + δ + σ2 + α). 7
152 153 154 155 156 157 158 159 160 161 162
163
164
The basic reproduction number has the key role in the disease modeling and is defined as the average number of new infected cases produced by a single infectious individual introduced in a totally susceptible population. Usually it explores the global dynamics of the disease whether the disease is controllable or not or it can be spread or not. When the value of the basic reproduction number falls less than 1, then the disease will not be spread in the community, whereas, when its values cross the 1, then the disease can be spread in the community i.e., the disease becomes endemic and will produce deaths and sometimes outbreak and epidemics. Therefore, the basic reproduction number is a threshold parameter for a disease and its biological significance of is obvious. The basic reproduction number R0 of the model is obtain by using the next generation method given in [16]. The associated matrices for the computations of R0 is given by βSI k L − (1 − η)δT 1 N F = 0 , V = k2 I − ǫL − ηδT . k3 T − γI 0 Then the linearization of system (1) at DFE is given by 0 β 0 F = 0 0 0 , 0 0 0 k1 0 −(1 − η)δ . −ηδ V = −ǫ k2 0 −γ k3
It follows that the basic reproduction number, denoted by R0 = ρ(F V −1 ) is given by R0 = R01 + R02 + R03 ,
165
with R01 =
166 167
168 169
170 171 172
173
174
βǫ γδǫ(1 − η) γδη , R02 = , R03 = . k1 k2 k1 k2 k3 k2 k3
In order to explore that R0 plays a key role to determine the global dynamics of the model (1), we state the following result (see Theorem 2, [16]). Theorem 1. The DFE of the system (1) about an equilibrium point Z0 , is locally asymptotically stable if and only if R0 < 1, otherwise unstable. Proof. To proof the result, we follow the hypothesis given in [16], namely (A1 ) − (A5 ). The hypothesis (A1 ) − (A4 ) are easy to verify, while to check A5 one need to show that the eigenvalues of thefollowing square matrix J|E0 have negative real parts M 0 , J|E0 = J3 J4 −k1 β (1 − η)δ −k1 0 (1 − η)δ , J4 = ǫ −k2 . ηδ ηδ where, J3 = −F , M = F −V = ǫ −k2 0 γ −k3 0 γ −k3 8
175 176 177
Define s(M) = max{Reλ : λ is an eigenvalue of M}, where s(M) is the simple eigenvalue of matrix M having a positive eigenvector, then from [16], we have if R0 < 1 ⇔ s(M) < 0. Further, calculate the eigenvalues of J4 , s(J4 ) = max(−k1 , −k2 , −k3 ) < 0.
179
therefore, if R0 < 1, then s(J|E0 ) < 0, the DFE point Z0 of system (1) is locally asymptotically stable.
180
4.1
178
181 182
183 184
Global stability analysis of DFE
Theorem 2. For R0 <1, the DFE Z0 of the system (1), is stable globally asymptotically on Ω and unstable for R0 >1. Proof. To establish the global stability of the disease-free equilibrium, we construct the following Lyapunov function: L(t) = A1 L + A2 I + A3 T,
185 186
where Ai , for i = 1,2,3, are some positive constants to be chosen later. Calculating the time derivative of L(t) along the solutions of system (1), we obtain: βSI dL(t) − k1 L + (1 − η)δT + A2 ǫL + ηδT − k2 I + A3 γI − k3 T = A1 dt N ≤ A1 βI − k1 L + (1 − η)δT + A2 ǫL + ηδT − k2 I + A3 γI − k3 T , S ≤ N = A1 β + A3 γ − A2 k2 I + [A2 ǫ − A1 k1 ]L + [A1 (1 − η)δ + A2 ηδ − A3 k3 ]T A1 β + A3 γ ≤ A2 k2 − 1 I + [A2 ǫ − A1 k1 ]L + [A1 (1 − η)δ + A2 ηδ − A3 k3 ]T. A2 k2 Now choosing A1 = ǫ, A2 = k1 and k3 =
187
ǫ(1−η)δ+k1 ηδ k3
and after simplification we get
dL(t) ≤ k1 k2 R0 − 1 I. dt
(3)
190
Hence if R0 < 1 then dL(t) is negative. Therefore the largest compact invariant set in dt Ω is the singleton set Z0 . Hence LaSalle’s invariant principle [8] then implies that Z0 is globally asymptotically stable in Ω.
191
5
192
Lemma 3. A unique positive endemic equilibrium (EE) exists when R0 exceeds 1.
188 189
Existence and stability of endemic equilibrium
9
193 194
Proof. The EE of (1) at E1 = (S ∗ , L∗ , I ∗ , T ∗) is obtained by setting the right hand side of the system (1) equal to zero, we get ∗ ∗ S = N , R 0 ((µ+σ1 )+γ(µ+σ2 +α)+γδ(1−η))N ∗ (R0 −1) ∗ L = R0 {(µ+σ1 )+γ(µ+σ2 +α)+γδ(1−η)+ǫ(γ+k3 )} , I∗ = ∗ T =
ǫk3 N ∗ (R0 −1) , R0 {(µ+σ1 )+γ(µ+σ2 +α)+γδ(1−η)+ǫ(γ+k3 )} ǫγN ∗ (R0 −1) . R0 {(µ+σ1 )+γ(µ+σ2 +α)+γδ(1−η)+ǫ(γ+k3 )}
195
Clearly, from above we conclude that a unique EE exists whenever R0 > 1.
196
5.1
Local stability of EE
197
198 199
Theorem 3. For R0 > 1, then the EE of the model (1) at E1 is locally asymptotically stable. Proof. The Jacobian matrix JE1 of the system (1) around the EE E1 is obtained by
J(E1 ) = 200
−I
∗ (N ∗ −S ∗ )β
−µ
N ∗2 I ∗ (N ∗ −S ∗ )β N ∗2
0 0
I ∗S∗ β N ∗2 I ∗S∗β − N ∗2 −
k1
− (N
ǫ 0
∗ −I ∗ )S ∗ β
N ∗2 (N ∗ −I ∗ )S ∗ β N ∗2 −k2
γ
I ∗ S∗β N ∗2
δ(1 − η) − δη −k3
I ∗S∗β N ∗2
.
The characteristic equation associated to J(E1 ) is λ4 + b1 λ3 + b2 λ2 + b3 λ + b4 = 0,
201
(4)
where the coefficients are
I ∗ βS ∗ b1 = k1 + k2 + k3 + µ + , N ∗2 1 b2 = k3 (µN ∗ 2 + IβS ∗ ) + QN ∗ 2 + k2 (µN ∗ 2 + IβS ∗ ) + k1 ((k2 + k3 + µ)N ∗ 2 + I ∗ β(N ∗ − S ∗ )) N ∗2 ∗ ∗ ∗ ∗ +βS (I µ + (N − I )ǫ) , 1 b3 = k1 (k3 (µN ∗ 2 + I ∗ β(N ∗ − S ∗ )) + QN ∗ 2 + k2 (µN ∗ 2 + I ∗ β(N ∗ − S ∗ ))) + k2 (k3 (µN ∗ 2 + I ∗ βS ∗ ) N2 +I ∗ βµS ∗ ) + βk3S ∗ I ∗ + βS ∗ I ∗ (γ + µ(1 + I ∗ N ∗ )ǫ) − (βk3 S(N ∗ I ∗ )ǫ + βγδη ∗2 +γδN ((1 − η)ǫ + ηµ)) , 10
b4
202
1 ∗2 = + I ∗ β(N ∗ − S ∗ )) + βS ∗ I ∗ (k2 k3 + ǫγµ) − (βk3µS ∗ ǫ(N ∗ − I ∗ ) 2 k1 Q(µN ∗ N ∗2 ∗ ∗ ∗ ∗ ∗ +γµ(δ(1 − η)ǫN + βδηS I ) + βγδI (1 − η)ǫ(N − S )) ,
where Q = {(µ + σ1 )k3 + γ(µ + σ2 + α) + γδ(1 − η)}.
203 204 205
206
207 208 209
The the eigenvalues of the characteristic polynomial (4) will be negative if it fulfil the Routh-Hurtwiz conditions that are bi > 0 for i = 1, 2, 3, 4 and b1 b2 b3 > b23 + b21 b4 . Hence then the EE of the model (1) will be locally asymptotically stable if R0 > 1.
5.2
Backward bifurcation analysis
Here, we discuss the existence of backward bifurcation of the system (1). To analyze this, we follow the center manifold theory described in [3]. If we take β as a bifurcation parameter, then at R0 = 1, we get β∗ = β =
210 211 212
k1 k2 k3 − γδ(1 − η)ǫ − γδηk1 . k3 ǫ
The following variations are made in the variables of the system (1) so that S = x1 , L = x2 , I = x3 , and T = x4 . Further using vector notation x = (x1 , x2 , x3 , x4 )T the model (1) can then be written in the form dx = F , with F = (f1 , f2 , f3 , f4 )T as shown below dt dx1 =Λ− dt dx2 = dt
βx1 x3 N
βx1 x3 N
− µx2 ,
− k1 x2 + (1 − η)δx4 ,
dx3 = ǫx2 + ηδx4 − k2 x3 , dt dx4 = γx3 − k3 x4 , dt 213
(5)
where N = Σ4i=1 xi . The Jacobian matrix evaluated at disease-free equilibrium Z0 with β ∗ is 1 +k1 k2 k3 −µ 0 − −γδǫ(1−η)−γδηk 0 ǫk3 1 +k1 k2 k3 0 −k1 −γδǫ(1−η)−γδηk δ(1 − η) ǫk 3 J= . 0 ǫ −k2 δη 0 0 γ −k3 11
(6)
214 215 216
The Jacobian matrix J has a simple zero eigenvalue calculated at β ∗ . The right and left eigenvectors denoted by W = (w1 , w2 , w3 , w4 ) and V = (v1 , v2 , v3 , v4 ) respectively and obtained w1 = −
w3 (γδηǫ − γδǫ − γδηk1 + k1 k2 k3 ) , k3 µǫ
217
w2 = 218
w3 (k2 k3 − γδη) γw3 > 0, w3 > 0, w4 = , k3 ǫ k3
and
k1 v2 δv2 (ηk1 + ǫ(1 − η)) , v4 = . ǫ k3 ǫ Computation of a: To evaluate a, we use the following partial derivatives v1 = 0, v2 > 0, v3 =
219
−2βµ ∂ 2 f2 ∂ 2 f2 −βµ ∂ 2 f2 = , = = , 2 ∂x3 Λ ∂x2 ∂x3 ∂x3 ∂x4 Λ 220
we get −2βµ 2 a= v2 w2 w3 + v2 w3 + v2 w3 w4 . Λ
221
(7)
Computation of b: We find the following partial derivatives to evaluate b ∂ 2 f2 = 1, ∂x3 ∂β
222
we obtain b = v2 w3 > 0.
(8)
224
Since all the terms in the expressions of b are positive, therefore, b > 0. Hence, the presence of the bifurcation at β = β ∗ depends solely on the sign of a.
225
5.3
223
226 227 228
229
Global stability of endemic equilibrium (EE)
In this subsection we prove the global asymptotic stability property of the EE of the proposed model (1). At the EE, we have from system (1) by using the approach given in [19, 20] ∗ ∗ Λ = βSN ∗I + µS ∗ , ∗ ∗ k1 L∗ = βSN ∗I + (1 − η)δT ∗ , (9) k2 I ∗ = ǫL∗ + ηδT ∗ , ∗ γI = k3 T ∗ . Theorem 4. If R0 > 1, then the EE of the model (1) is globally asymptotically stable. 12
230
231
232
Proof. To show the result, we define the following Lyapunove function S(t) I(t) ∗ ∗ ∗ ∗ K(t) = ǫ S(t) − S − S ln ∗ + k1 I(t) − I − I ln ∗ S I ∗ T (t) δT (k1 η + ǫ(1 − η)) + T (t) − T ∗ − T ∗ ln ∗ . ∗ γI T Taking time derivative of (10) along the model solution, we get S∗ ′ L∗ ′ I∗ ′ dK(t) = ǫ (1 − )S + (1 − )L + k1 (1 − )I dt S L I ∗ ∗ T δT (k1 η + ǫ(1 − η)) + (1 − )T ′ . ∗ γI T
(10)
(11)
It follows from direct calculation: i S ∗ h βSI S∗ ′ S = ǫ 1− Λ− − µS ǫ 1− S S N i S ∗ h βS ∗ I ∗ βSI ∗ = ǫ 1− + µS − − µS S N∗ N S∗ S ǫβS ∗ I ∗ S ∗ IN ∗ ∗ = ǫµS 2 − 1− − ∗ + + ∗ S S N∗ S I N ǫβSI . − N
(12)
i L∗ ′ L∗ h βSI ǫ 1− L = ǫ 1− − k1 L + (1 − η)δT L L N ǫβSI ǫβSIL∗ T L∗ = − − k1 ǫL + k1 ǫL∗ + (1 − η)ǫδT − (1 − η)δǫ N NL L ∗ ∗ ∗ ǫβSI ǫβSIL ǫβS I = − − k1 ǫL + + (1 − η)δǫT ∗ + (1 − η)ǫδT N NL N∗ T L∗ −(1 − η)δǫ . (13) L I∗ ′ I ∗ h k1 1 − I = k1 1 − ǫL + ηδT − k2 I I I I∗ I∗ = k1 ǫL − k1 ǫL + k1 ηδT − k1 ηδT − k1 k2 I + k1 k2 I ∗ I I ∗ I L I ∗ ǫβS ∗ I ∗ = k1 ǫL − k1 ǫL∗ ∗ + k1 ηδT − k1 ηδT + + k1 ηδT ∗ + ǫ(1 − η)δT ∗ IL I N∗ ǫβS ∗ I ∗ I I I − − (1 − η)ǫδT ∗ ∗ − k1 ηδT ∗ ∗ ∗ ∗ N I I I
13
∗ I ∗ ǫβS ∗ I ∗ ǫβS ∗ I ∗ LI ∗ ∗ LI − (1 − η)ǫδT + k ηδT − k ηδT + + k1 ηδT ∗ + 1 1 N ∗ L∗ I L∗ I I N∗ ǫβS ∗ I ∗ I ∗ I ∗ I ǫ(1 − η)δT ∗ − − (1 − η)ǫδT − k ηδT . (14) 1 N∗ I∗ I∗ I∗
= k1 ǫL −
233
T∗ δT ∗ (k1 η + ǫ(1 − η)) T∗ I∗ δT ∗ (k1 η + ǫ(1 − η)) 1 − [γI − k T ] = 1 − [I − T] 3 γI ∗ T I∗ T T∗ I I IT ∗ = δηk1 T ∗ ∗ + δ(1 − η)ǫT ∗ ∗ − δηk1 T ∗ ∗ I I I T ∗ IT −δ(1 − η)ǫT ∗ ∗ − δηk1 T − δ(1 − η)ǫT I T +δηk1 T ∗ + δ(1 − η)ǫT ∗ . 234
235
236 237 238
239
240 241 242 243 244
(15)
Using equations (12) to (15) in (11) we get dK(t) S∗ S ǫβS ∗ I ∗ S∗ I LI ∗ SIL∗ N ∗ LS ∗ ∗ = ǫµS 2 − − ∗ + 3 − − − − (1 − ) dt S S N∗ S I ∗ L∗ I S ∗ I ∗ NL L∗ S I ∗ L L∗ T T ∗I ∗ +(1 − η)ǫδT 3 − ∗ − − IL LT ∗ T I ∗ I ∗T IT ∗ ∗ − . (16) +δηk1 T 2 − IT ∗ I ∗ T Using properties of the arithmetic mean and geometric mean in equation (16) we have S S∗ 2 − − ≤ 0, S∗ S S∗ I LI ∗ SIL∗ N ∗ LS ∗ 3 − − − − (1 − ) ≤ 0, S I∗ L∗ I S∗I ∗N L L∗ S ∗ ∗ ∗ 3 − IILL∗ − LLTT∗ − TT I I∗ ≤ 0, I∗T IT ∗ 2 − IT ∗ − I ∗ T ≤ 0.
≤ 0 when R0 >1. As, all the parameters are non-negative, therefore it follows that dK(t) dt ∗ Therefore, by the LaSalles Invariance Principle [8], (S, L, I, T ) → (S , L∗ , I ∗ , T ∗ ) as t → ∞.
6
Sensitivity analysis of R0
Here, we analyzed the effect of the various model parameters to define the relative importance of these parameters to disease transmission. For this purpose, we calculate the partial derivatives of threshold number R0 with respect to model parameters γ, δ, β and 0 η. Parameter η is the failure rate of treatment and ∂R > 0. Therefore, the total number ∂η of infected population can be decreased by reducing the treatment failure rate η. Also, 14
245 246 247 248
∂R0 ∂δ
0 > 0, therefore, the parameter δ have positive effect. Since ∂R > 0 which implies that ∂β by decreasing the transmission rate β, the infection can be reduced. On the other hand, 0 < 0. It means that by increasing the parameter γ, the since the partial derivatives ∂R ∂γ TB infection can be controlled.
Parameter β α γ δ η ǫ µ σ1 σ2
Description Disease contact rate Progression rate from T to R Transmission rate from I to T Rate at which treated individuals leave the T Treatment failure rate Rate of moving from L to I Natural death rate Death rate in I Death rate in T
Sensitivity +0.7221 -0.066192 -0.272156 +0.0778308 +0.253903 +0.0655175 -0.0962639 -0.42162 -0.0091748
Table 3: Sensitivity indices of R0 with respect to model parameters. 249 250
Definition 1. The normalized sensitivity index of R0 depending on the differentiability on a parameter Ψ is defined as follows:
260
Ψ ∂R0 . (17) R0 ∂Ψ We calculate the sensitivity indices of the basic reproduction number R0 with respect to model parameters using above formula given in (17) and are tabulated in Table 3. From Table 3 we can see that parameters β, δ, η and ǫ have a positive sign which means that R0 increases with the parameter. While the rest of parameters have negative sign means that R0 decreases for higher values of the parameters. Further, this implies that by increasing (or decreasing) of these parameters by 10% will increase (or decrease) the basic reproductive number by 7.221%, 0.7783%, 2.539% and 0.6551% respectively. On the other hand the parameters α, γ, µ, σ1 , and σ2 have negative influence hence, increasing these rates by 10% will decrease the basic reproduction number by 0.6819%, 2.721 %, 0.9621%, 4.2162% and 0.0917%.
261
7
SΨR0 =
251 252 253 254 255 256 257 258 259
262 263 264 265 266 267 268
Numerical results and Discussion
To illustrate the numerical results, the following values were taken for the initial conditions to mimic the Khyber Pakhtunkhwa of Pakistan. According to the Pakistan Bureau of Statistics the population of Khyber Pakhtunkhwa in the 2017 census is estimated as 30,523,371 [4]. We take the total initial population as N(0) = 30, 523, 371 . We assume that initial population of exposed individuals L(0) = 83000. The initial infected is as given in the data is I(0) = 8010. We assume there are no treated and recovered individuals initially, so T (0) = 0 and R(0) = 0 respectively. Hence, the initial susceptible population 15
269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285
286
287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310
is given as S(0) = N(0) − L(0) − I(0) − T (0) − R(0) = 30432361. The numerical values parameters fitted from real data are presented in the Table 2. Fig. 3, shows the simulation of the model (1), when R0 = 0.8068 < 1. It is clear that all solutions converged to the DFE, Zo . In Fig. 4, the simulation of the model (1) is shown when R0 = 1.3419 > 1. Clearly all solutions converged to the EE, E1 . Fig. 5, shows that by taking different initial values all infected individuals vanishes when R0 = 0.8068 < 1, β = 0.2433. In Fig. 6, it is clear that when R0 = 1.3419 > 1 and β = 0.5433 and by considering diffract initial conditions then all infected classes converged to the EE, E1 . Fig. 7, explains the effect of treatment rate γ on the number of cumulative new cases of infection. Clearly from Fig. 7, by increasing treatment parameter γ the total TB infective individuals decreases. In Fig. 8, the effect of treatment failure parameter η is shown. The total infected individuals decreases when η decreases. The effect of parameter δ is simulated in Fig. 9, describing that by decreasing values of δ the total infected population decreases. The behavior of TB infected cases with and without treatment when R0 = 0.8068 < 1 and β = 0.2433, is given in Fig. 10. Finally, in Fig. 11, we describes the dynamics of infected individuals with and without treatment for R0 = 1.3419 > 1 and β = 0.5433. Figs. 13 to 18 show the behavior of the sensitive parameters to the basic reproduction number R0 .
8
Conclusion
In this paper, we formulated and analyzed the transmission dynamics of TB in Khyber Pakhtunkhwa province in Pakistan. The model was parameterized by using the data since 2002-2017 obtained from NTP. The stability results for the model are obtained and found stable locally as well globally when the basic reproduction number R0 less or greater than unity respectively. The numerical results for various values were presented. The government of Pakistan initiated various programs to eliminate TB infection in the last two decades. We developed the TB incidence model to simulate the TB data of Khyber Pakhtunkhwa reported by NTP and WHO [11, 1], from 2002-2017. The value of the basic reproduction number R0 and biological parameters of the model were estimated on the basis of available data and are tabulated in Table 2. To eliminate the TB epidemics it is necessary to ensure that R0 < 1. From numerical simulations, it is clear that the basic reproduction number R0 depends on the transmission rate β, treatment rate γ, treatment failure parameter η and on the proportion δ, at which the treated individuals leave the treatment class. Therefore, to control and eradicate the TB infection, it is important to consider the following strategies. The first and most important strategy is to minimize the contact with TB infected individuals by decreasing the value of β. Second important strategy is to increase the treatment rate γ of infective individuals. Also, to decrease the TB infective we need to reduce the treatment failure cases. Although the most effective strategy is to isolate infective by decreasing the value of transmission coefficient β. However, in developing countries like Pakistan it is very difficult to isolate the infectious individuals because of the high cost of the long term treatment. Therefore, we proposed the most feasible and optimal strategy to eliminate the TB infection is to increase the treatment rate by decrease the treatment cost so that the poor people can get the treatment measure, and by proper follow up of under treatment population so that the failure of treatment 16
9
×10 4 L I T R
8
Behavior of individuas
7 6 5 4 3 2 1 0 0
50
100
150
200
250
300
350
400
450
Time(Years)
Figure 3: The plot shows the behavior of individuals when R0 = 0.8068 < 1. 311
cases decreases.
312
Conflict of interests
313
The author’s declare that no competing interests exists.
314
Funding
315
No source of funding for this study.
17
3.5
×10 7 S L I T R
3
Endemic Equiibrium
2.5
2
1.5
1
0.5
0 0
100
200
300
400
500
600
700
800
900
1000
Time(Years)
Figure 4: The plot shows the behavior of individuals when R0 = 1.3419 > 1.
10
×10 4
Tota number of infected cases
9 8 7 6 5 4 3 2 1 0 0
50
100
150
200
250
300
350
400
450
Time(Years)
Figure 5: The plot shows the total number of infected individuals with different initial conditions with R0 = 0.8068 < 1 and β = 0.2433. 18
8
×10 6
Total number of infected cases
7
6
5
4
3
2
1
0 0
50
100
150
200
250
300
Time(Years)
Figure 6: The plot shows the total number of infected individuals with different initial conditions with R0 = 1.3419 > 1 and β = 0.5433.
8
×10 6 γ=0.2873 γ=0.3873 γ=0.4873 γ=0.5873 γ=0.6873 γ=0.7873
Tota number of infected cases
7
6
5
4
3
2
1
0 0
50
100
150
200
250
300
Time(Years)
Figure 7: The plot shows the total number of infected individuals by changing the value of γ.
19
8
×10 6 η=0.15 η=0.1 η=0.01
Total number of infected cases
7
6
5
4
3
2
1
0 0
50
100
150
200
250
300
Time(Years)
Figure 8: The plot shows the total number of infected individuals for different values of η.
8
×10 6 δ=1.1996 δ=1.0 δ=0.8 δ=0.6 δ=0.4 δ=0.2
Total number of infected cases
7
6
5
4
3
2
1
0 0
50
100
150
200
250
300
Time(Years)
Figure 9: The plot shows the total number of infected individuals for different values of δ.
20
10
×10 4 With treatment Without treatment
Total number of infected cases
9 8 7 6 5 4 3 2 1 0 0
20
40
60
80
100
120
140
160
180
200
Time(Years)
Figure 10: The plot shows the comparison of the total number of infected individuals with and without treatment with R0 = 0.8068 < 1 and β = 0.2433.
10
×10 6 With treatment Without treatment
Total number of infected cases
9 8 7 6 5 4 3 2 1 0 0
50
100
150
200
250
300
Time(Years)
Figure 11: The plot shows the comparison of the total number of infected individuals with and without treatment with R0 = 1.3419 > 1 and β = 0.5433. 21
1
2.5
1
1.2
1.4
2 1.8
1.6
2
0.9
1.8
0.8 0.7
1.5
0.6
1
0.5
δ
R0
2
1.6
1.4
0.8
1.2
1.4
2 1.8 1.6
0.3
0 1
1.2
1
0.4
0.5
1
0.2 1 0.8
0.5
8 0.
0.1
0.6
0.8
0.6
0.4
δ
0
0
0.2 0
0.6
0
γ
0.2
0.4
0.6
0.8
1
γ
(a)
(b)
Figure 12: Effects of parameters γ and δ on the basic reproduction number R0 .
1
1
1.2
1.4
1.9
0.8
1.8
0.7
1.7
0.6
1.6
0.5
1.5
0.4
1.4
1
1.2
1.4
0.3
0.5 1
1.6
1
2 1.8
1.5
η
R0
2
2 1.8
0.9 2.5
1.6
2
1.3
0.2
1.2
0.1
1.1
1 0.8
0.5
0.6 0.4
η
0
0
0.2 0
1
0
γ
0.2
0.4
0.6
0.8
1
γ
(a)
(b)
Figure 13: Effects of parameters γ and η on the basic reproduction number R0 .
22
1
0.8
1.8
0.7
1.7
0.6
1.6
1.5
α
R0
1
1.9
1.2
1.4
2
1.6
2.5
2
1.8 2
0.9
0.5
1.5
1 1
0.4
1.2
1.4
1.6
0.5 1
1.4
1.8 2
0.3
1.3
0.2
1.2
0.1
1.1
1 0.8
0.5
0.6 0.4
α
0
0
0.2 0
1
0
γ
0.2
0.4
0.6
0.8
1
γ
(a)
(b)
Figure 14: Sensitivity analysis of the basic reproduction number R0 verses parameters γ and α.
1
2
0.9 2.5
1.8
0.8 1.6
η
R0
2
1.8
1.6
1.4
1.2
1
1
0.6
0.8
1.5
0.6
0.7
0.4
2
1.4
0.5
1.2
0.4
0.5
1
0.3
0 1
0.8
0.2 0.6
2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.1
0.6
0.4
1 0.8
0.5 0.4
η
0
0
0.2 0
0.4
0
β
0.2
0.4
0.6
0.8
1
β
(a)
(b)
Figure 15: Effects of parameters β and η on the basic reproduction number R0 .
23
2
3.5
1.8
3.5
3.5
3
1.6 3
3
1.4
2.5
2.5
δ
1.2 1
2.5 2
2
2
0.8 1.5
0.6 0.4
1.5
1.5
1
1
1
0.2 0.5
0.5
0 0
0.5
1
0.5
1.5
2
β
(a)
(b)
Figure 16: Effects of parameters β and δ on the basic reproduction number R0 .
1
3.5
0.9 3
4
0.8 0.7
2.5
1
3
2
1
2
5
0.5
1.
γ
0.4
0.5
0.2 1 0.8
1
0.1
0.6
1
0.5 0.4
γ
0
0
0
β
2.5
2
0
0.2
1.5
2
0.3
0 1
1. 5
R0
0.6
0.2
0.4
0.6
3 0.8
0.5
1
β
(a)
(b)
Figure 17: Effects of parameters β and γ on the basic reproduction number R0 . 316
317 318 319
320 321
322 323
References [1] World Health Organization Media Centre. Available 2018. https://www.who.int/en/news-room/fact-sheets/detail/tuberculosis. accessed 2018. [2] Carlos Castillo-Chavez and Zhilan Feng. To treat or not to treat: the case of tuberculosis. Journal of mathematical biology, 35(6):629–656, 1997. [3] Carlos Castillo-Chavez and Baojun Song. Dynamical models of tuberculosis and their applications. Mathematical biosciences and engineering, 1(2):361–404, 2004. 24
324 325
326 327
328 329
[4] World Health Organization. WHO country cooperation strategic. http://apps.who.int/iris/bitstream/10665/136607/1/ccsbrief pak en.pdf, 2016. [5] Z Feng and C Castillo-Chavez. Mathematical models for the disease dynamics of tubeculosis, 1998. [6] David J Gerberry. Practical aspects of backward bifurcation in a mathematical model for tuberculosis. Journal of theoretical biology, 388:15–36, 2016.
332
[7] Soyoung Kim, A Aurelio, and Eunok Jung. Mathematical model and intervention strategies for mitigating tuberculosis in the philippines. Journal of theoretical biology, 443:100–112, 2018.
333
[8] Joseph P LaSalle. The stability of dynamical systems, volume 25. Siam, 1976.
330 331
334 335
336 337
338 339
340 341
342 343 344
345 346 347
348 349
350 351 352
353 354 355
356 357
[9] Junli Liu and Tailei Zhang. Global stability for a tuberculosis model. Mathematical and Computer Modelling, 54(1-2):836–845, 2011. [10] Luju Liu, Xiao-Qiang Zhao, and Yicang Zhou. A tuberculosis model with seasonality. Bulletin of Mathematical Biology, 72(4):931–952, 2010. [11] National TB Control Program http://www.ntp.gov.pk/webdatabase.php.
Pakistan
(NTP).
[12] Pakistan Bureau of Statistics. Pakistans 6th census: Population of major cities 583 census. 584. http://www.pbscensus.gov.pk/. [13] Charles S Revelle, Walter R Lynn, and Floyd Feldmann. Mathematical models for the economic allocation of tuberculosis control activities in developing nations. American review of respiratory disease, 96(5):893–909, 1967. [14] James M Trauer, Justin T Denholm, and Emma S McBryde. Construction of a mathematical model for tuberculosis transmission in highly endemic regions of the asia-pacific. Journal of theoretical biology, 358:74–84, 2014. [15] Saif Ullah, Muhammad Altaf Khan, and Muhammad Farooq. A fractional model for the dynamics of tb virus. Chaos, Solitons & Fractals, 116:63–71, 2018. [16] Pauline Van den Driessche and James Watmough. Reproduction numbers and subthreshold endemic equilibria for compartmental models of disease transmission. Mathematical biosciences, 180(1-2):29–48, 2002. [17] Hans Waaler, Anton Geser, and Stig Andersen. The use of mathematical models in the study of the epidemiology of tuberculosis. American Journal of Public Health and the Nations Health, 52(6):1002–1013, 1962. [18] Robert S Wallis. Mathematical models of tuberculosis reactivation and relapse. Frontiers in microbiology, 7:669, 2016.
25
358 359
360 361
362 363 364
365 366
[19] Yi Wang and Jinde Cao. Global dynamics of multi-group SEI animal disease models with indirect transmission. Chaos, Solitons & Fractals, 69:81–89, 2014. [20] Yi Wang and Jinde Cao. Global stability of general cholera models with nonlinear incidence and removal rates. Journal of the Franklin Institute, 352(6):2464–2485, 2015. [21] Yali Yang, Jianquan Li, Zhien Ma, and Luju Liu. Global stability of two models with incomplete treatment for tuberculosis. Chaos, Solitons & Fractals, 43(1-12):79–85, 2010. [22] Jinhui Zhang, Yong Li, and Xinan Zhang. Mathematical modeling of tuberculosis data of china. Journal of theoretical biology, 365:159–163, 2015.
26
PII: DOI: Reference:
S0378-4754(19)30108-9 https://doi.org/10.1016/j.matcom.2019.03.012 MATCOM 4768
To appear in:
Mathematics and Computers in Simulation
Received date : 18 May 2018 Revised date : 27 January 2019 Accepted date : 26 March 2019 Please cite this article as: S. Ullah, M.A. Khan, M. Farooq et al., Modeling and analysis of Tuberculosis (TB) in Khyber Pakhtunkhwa, Pakistan, Mathematics and Computers in Simulation (2019), https://doi.org/10.1016/j.matcom.2019.03.012 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Manuscript Click here to view linked References
Modeling and analysis of Tuberculosis (TB) in Khyber Pakhtunkhwa, Pakistan January 28, 2019 Saif Ullah† , Muhammad Altaf Khan†† ,1 Muhammad Farooq † ,Taza Gul
††
† Department
of Mathematics, University of Peshawar, Peshawar, Khyber Pakhtunkhwa, Pakistan. †† Department of Mathematics, City University of Science and Information Technology, Peshawar, Khyber Pakhtunkhwa, Pakistan. Abstract
1
In this paper, we present a mathematical model on the dynamics of Tuberculosis (TB) in Khyber Pakhtunkhwa, province of Pakistan. The model is parameterized using TB infected cases reported from the year 2002 to 2017 in Khyber Pakhtunkhwa. The basic reproduction number R0 for the given period is estimated as R0 ≈ 1.3419. Mathematical analysis of the model in terms of stability analysis is derived and discussed in detail. The model at the disease free case is locally as well as globally asymptotically stable when the basic reproduction number R0 < 1. Further, when the basic reproduction number R0 > 1, then the model is stable locally as well as globally. Moreover, the sensitivity analysis of the model parameters are performed and their corresponding graphical results are presented. Finally, the numerical results for the sensitive parameters are plotted and their effects are shown on the disease elimination.
2 3 4 5 6 7 8 9 10 11 12 13
14 15
Key words: Tuberculosis (TB) Model, Parameter estimation, Lyapunov function, Stability.
16
1
17 18 19 20 21 22
Introduction
Tuberculosis (TB) is a bacterial infectious disease caused by bacillus Mycobacterium tuberculosis (MTB) and is one of the top 10 causes of death worldwide. TB infection mainly affects the lungs (known as pulmonary TB). It may also invade other organs including brain, spine, kidneys, central nervous system or lymphatic system (known as extrapulmonary TB). TB is also known as air-born infection because it is transmitted from TB active people thorough air when they cough, spit, sneeze or speak. The main symptoms of 1
Corresponding author. Email:[email protected]
1
23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67
TB infected person are chronic cough, blood-containing sputum, night sweats, fever and weight loss. According to World Health Organization (WHO) recent reports, abut (1/4)th of the whole world’s population are exposed to TB infection (which means the people infected by TB but not yet propagate the infection). More than 10.4 million people were infected with TB and about 1.7 million expire due to this infection in 2016. The person once infected with TB have 5-15 percent lifetime risk of falling ill with TB [1]. Mostly, the TB related deaths occur in the middle income countries such as India, Indonesia, China, Philippines, Pakistan, Nigeria, and South Africa. These countries accounted more than 60% of the whole TB burden in the world [1]. Mathematical modeling is one of the important tools to explore the dynamics of a disease and provide useful techniques to control it. In the last few decades, several biologists and mathematicians have been developed different epidemic models for the transmission dynamic of TB in different areas around the globe. Waaler et al. [17] introduced the first mathematical model of TB in 1962, by dividing the whole population into three classes. The transmission model of TB infection depending on the proportion of its prevalence is develop in [13]. A two strain, an age structure and time delay TB transmission models were proposed by Castillo et al. [2, 5]. Yang et al. [21] analyzed the global stability of TB models with incomplete treatment. A similar stability analysis of TB mathematical model was discussed by Liu et al. in [9]. A fractional TB model with Caputo derivative is investigated in [15]. The most powerful method to show the global dynamics of epidemic models is that of the Lyapunov which can be found in recent papers such as [19, 20]. Most of the models are based on the qualitative behavior of TB infection. However there are few articles presenting the quantitative analysis of TB data. In [10], Liu and co-authors investigated the TB dynamics with seasonality and use the epidemiological data to estimate the model parameters. Trauer et al. [14] construct a realistic model to study the dynamics of TB in different regions of Asia-pacific. Zhang et al. [22] proposed the TB mathematical model with hospitalized and non hospitalized infective population and use the TB data of China to simulate his model. Gerberry examined the TB transmission model and explored the practical aspects of backward bifurcation [6]. A TB model with relapse and reinfection is analyzed by Robert [18]. Recently, mathematical model with optimal control strategies for TB dynamics in the Philippines has been studied by Kim et al. [7]. In Pakistan, TB is a leading public health problem and one of the major causes of morbidity and mortality. Pakistan is one of the highest countries with TB infection and ranked fifth among the twenty two countries having high burden of this infection. Annually more than 0.5 million new TB cases including 15000 children are reported and more than seventy thousand Pakistanis died because of TB infection each year. Due to high global prevalence of multidrug-resistant tuberculosis (MDR-TB), Pakistan ranked fourth in the world [12, 4, 1]. TB is also a major health problem and leading cause of death in the province of Pakistan named Khyber Pakhtunkhwa. According to national TB control program (NTP) about 462920 TB infected cases are registered and treated from 2002 till 2017 in Khyber Pakhtunkhwa [11]. In this paper we construct a mathematical model to investigate the dynamics of TB infection in Khyber Pakhtunkhwa. We use the new TB infected cases from 2002 to 2017 obtained from the NTP [11], to estimate the biological parameters of the proposed model. The structure of the paper is organized as follows: The complete analysis of the model 2
75
and the parameters estimation from availed data are presented in the next section “Methods and materials”. The fundamental properties of the model are presented in section 3. The stability analysis of the disease-free case is given in section 4. The existence and stability analysis of endemic equilibrium point is provided in section 5. Also, in this section the analysis of bifurcation is discussed. In section 6, the sensitivity analysis is presented. Numerical simulations for the model based on fitted values of parameters are obtained to illustrate the theoretical results in section 7. Finally, the conclusion is summarized in section 8.
76
2
77
2.1
68 69 70 71 72 73 74
78 79 80 81 82 83
Methods and materials Model formulation
In this section we present the proposed mathematical model of TB infection. To develop the model, total human population is divided into five epidemiological sub-compartments denoted by susceptible S(t), Exposed L(t), TB active I(t), under treatment T (t), and recovered individuals R(t). The transmission dynamics of TB infection is given by the following system of nonlinear differential equations: dS = Λ− dt dL = dt
βSI N
βSI N
− µS,
− (µ + ǫ)L + (1 − η)δT,
dI = ǫL + ηδT − (µ + γ + σ1 )I, dt
(1)
dT = γI − (µ + δ + σ2 + α)T, dt dR = αT − µR. dt 84 85 86 87 88 89 90 91 92 93 94 95
Where N(t) = S(t) + L(t) + I(t) + T (t) + R(t). Recruitment rate is Λ and the natural death rate is µ. The number of new infection per unit time due to contact with TB infective individuals is denoted by parameter β. The progression rate from class L to I is ǫ. The infective individuals are treated at the rate γ. The treated individuals leaves the class T at a rate δ. The term (1 − η)δT , where the parameter η reflects the treatment failure, further by considering η = 0, will means that all the individuals under treatment become latent, and when η = 1, means the failure of treatment and all the treated individuals will still be infectious. The treated individuals enter to either latent class due to the remainder of Mycobacterium tuberculosis or infective class I due to the failure of treatment at the rate δ. The disease mortality rates in the infective class I and treated class T are σ1 and σ2 respectively. The parameter α represents recovery rate of the treated individuals. All these parameters and their flow rate to each class is shown in Figure 1 while the model 3
97
variables are given in Table 1. Since the last equation of (1) is independent on the rest, so we ignore it and consider the first four equations of the model.
98
2.2
96
99 100 101 102 103 104 105 106 107 108 109 110 111
Estimation of model parameters
In this subsection we estimate the model parameters based on the actual TB incidence data taken from national TB control program, Pakistan [11]. In order to parameterize TB model (1), we obtained some of the parameter values from literature (see Table 2), others were estimated or fitted from the data. For instance, the demographic parameter, µ, is estimated as µ = 1/67.7 per year, where 67.7 years is the average lifespan in Pakistan [1]. The other demographic parameter, Λ, is then estimated as follows. Since the total population of Khyber Pakhtunkhwa province as at 2017 was 30,523,371 [4], we assumed that Λ/µ, which is the limiting total human population in the absence of the disease, is 30, 523, 371, so that Λ = 450, 862.20088626 per year. The rest of model parameters are estimated using least-square fitting technique which provides a better fit of the model solution to the real data as shown in Fig. 2. The best fitted parameter values are obtained by minimizing the error between the actual TB incidence data and the solution of the proposed model (1). The objective function used in the parameter estimation is as follows ˆ = argmin Θ
n X j=1
112 113 114 115 116 117 118
2
(Itj − I tj ) ,
where I tj denotes the actual cumulative TB infected case and Itj are the corresponding model solution at time tj , while n is the number of available actual data points. In Fig. 2, the actual incidence data is shown by red circles while the model fitted curve is shown by blue solid line. The associated state variables and parameters of the model (1) are tabulated in Table:1 and Table:2 respectively. Consequently, using these parameter estimates, we obtained, in this study, that the value of R0 for the 2002-2017 TB cases in Khyber Pakhtunkhwa Pakistan is R0 ≈ 1.3419. Variable Description S(t) L(t) I(t) T (t) R(t)
The number of individuals which attract the disease at any time t The number of individuals which is not yet infectious The number of TB active individuals t The number of treated individuals at any time t Recovered individuals at any time t Table 1: Description of state variables of the model (1).
4
Parameter Λ β α γ µ σ1 σ2 δ η ǫ
Description Recruitment rate Disease contact rate Progression rate from T to R Transmission rate from I to T Natural mortality rate Disease related death rate of infected individuals Disease related death rate in T Rate at which treated individuals leave the T Treatment failure rate Rate of moving from L to I
Baseline value 450,862.20088626 0.5433 0.3968 0.2873 1/67.7 0.2202 0.0550 1.1996 0.1500 0.2007
Reference Estimated Fitted Fitted Fitted [12] Fitted Fitted Fitted Fitted Fitted
Table 2: Numerical values for the parameters of the TB model (1) estimated from Khyber Pakhtunkhwa population and TB incidence data.
(µ + σ1 )I I µS Λ
S
ǫL
µL βSI N
ηδT
L
γI
(1 − η)δT
R
µR
αT T
(µ + σ2 )T
Figure 1: Flow chart.
119
3
Analysis of the model
121
In this section we prove the basic properties of the proposed TB model (1). We proceed as follow:
122
3.1
120
123 124
Positivity and boundedness of the solution
To show that the TB infection model (1) is epidemiologically meaningful, it is necessary to show that all its state variables are non-negative for all time. In other words, solutions of 5
5
×10 4 Data Model Fit
Tuberculosis Incidence Data
4.5 4 3.5 3 2.5 2 1.5 1 0.5 2002
2004
2006
2008
2010
2012
2014
2016
2018
Years
Figure 2: Data fitting (blue solid curve) of the number of TB cases in Khyber Pakhtunkhwa, using model (1). 125 126
127 128
the model system (1) with non-negative initial data will remain non-negative for all time t > 0. Lemma 1. For all t > 0 and the initial data P (0) ≥ 0, where P (t) = (S(t), L(t), I(t), T (t)), the solution of TB model (1) are non-negative if they exist. Furthermore, lim supN(t) ≤
t→∞
129 130
Λ . µ
Proof. Let χ(t) = βI and t1 = sup{t > 0 : P (t) > 0 ∈ [0, t]}, then the from first equation N of system (1) we have dS = Λ − χ(t)S(t) − µS, dt = Λ − (χ(t) + µ)S(t).
131 132
(2)
The assumption that the solution of the system (1) exists in the interval I, the solution of (2) given as below Z t1 Z t1 S(t1 ) = S(0)exp − µt1 + χ(x)dx + exp − µt1 + χ1 (x)dx 0 0 Z t1 Z p × Λ exp µp + χ(y)dy dp > 0. 0
0
6
133 134 135
Similarly, it can be shown that P > 0 for all for all time t > 0. For the second part of the proof, note that 0 < S(0) ≤ N(t), 0 < L(0) ≤ N(t), 0 < I(0) ≤ N(t), 0 < T (0) ≤ N(t). Adding components of the model (1) gives dN(t) = Λ − µN(t) − σ1 I − σ2 T ≤ Λ − µN(t). dt
136
Hence lim supN(t) ≤
t→∞
Λ . µ
137
138
139 140
141
3.2
Invariant regions
The TB model (1) will be studied in a biologically feasible region as given below. Consider the feasible region Ω ⊂ R4+ , with
Ω= 142 143
144
(S(t), L(t), I(t), T (t)) ∈
R4+
Λ . : N(t) ≤ µ
Lemma 2. The region Ω ⊂ R4+ , is positively invariant for the model (1) with non-negative initial conditions in R4+ . Proof. It follows from summing equations of the model (1) that dN(t) = Λ − µN(t) − σ1 I − σ2 L ≤ Λ − µN(t), dt
145 146 147 148
149
hence dNdt(t) ≤ 0, if N(0) ≥ Λµ . Thus N(t) ≤ N(0)eµt + Λµ (1 − eµt ). In particular, N(t) ≤ N(0)eµt . Thus, the region Ω is positively invariant. Furthermore, if N(0) > Λµ , then either the solutions enters Ω in finite time, or N(t) approaches Λµ asymptotically. Hence, the region Ω attracts all solutions in R4+ .
4
Stability analysis of disease free equilibrium
The model (1) has a disease free equilibrium (DFE), Z0 , given by Z0 = (S 0 , 0, 0, 0), 150
where S0 =
151
Λ . µ
Let denote k1 = (µ + ǫ), k2 = (µ + γ + σ1 ), k3 = (µ + δ + σ2 + α). 7
152 153 154 155 156 157 158 159 160 161 162
163
164
The basic reproduction number has the key role in the disease modeling and is defined as the average number of new infected cases produced by a single infectious individual introduced in a totally susceptible population. Usually it explores the global dynamics of the disease whether the disease is controllable or not or it can be spread or not. When the value of the basic reproduction number falls less than 1, then the disease will not be spread in the community, whereas, when its values cross the 1, then the disease can be spread in the community i.e., the disease becomes endemic and will produce deaths and sometimes outbreak and epidemics. Therefore, the basic reproduction number is a threshold parameter for a disease and its biological significance of is obvious. The basic reproduction number R0 of the model is obtain by using the next generation method given in [16]. The associated matrices for the computations of R0 is given by βSI k L − (1 − η)δT 1 N F = 0 , V = k2 I − ǫL − ηδT . k3 T − γI 0 Then the linearization of system (1) at DFE is given by 0 β 0 F = 0 0 0 , 0 0 0 k1 0 −(1 − η)δ . −ηδ V = −ǫ k2 0 −γ k3
It follows that the basic reproduction number, denoted by R0 = ρ(F V −1 ) is given by R0 = R01 + R02 + R03 ,
165
with R01 =
166 167
168 169
170 171 172
173
174
βǫ γδǫ(1 − η) γδη , R02 = , R03 = . k1 k2 k1 k2 k3 k2 k3
In order to explore that R0 plays a key role to determine the global dynamics of the model (1), we state the following result (see Theorem 2, [16]). Theorem 1. The DFE of the system (1) about an equilibrium point Z0 , is locally asymptotically stable if and only if R0 < 1, otherwise unstable. Proof. To proof the result, we follow the hypothesis given in [16], namely (A1 ) − (A5 ). The hypothesis (A1 ) − (A4 ) are easy to verify, while to check A5 one need to show that the eigenvalues of thefollowing square matrix J|E0 have negative real parts M 0 , J|E0 = J3 J4 −k1 β (1 − η)δ −k1 0 (1 − η)δ , J4 = ǫ −k2 . ηδ ηδ where, J3 = −F , M = F −V = ǫ −k2 0 γ −k3 0 γ −k3 8
175 176 177
Define s(M) = max{Reλ : λ is an eigenvalue of M}, where s(M) is the simple eigenvalue of matrix M having a positive eigenvector, then from [16], we have if R0 < 1 ⇔ s(M) < 0. Further, calculate the eigenvalues of J4 , s(J4 ) = max(−k1 , −k2 , −k3 ) < 0.
179
therefore, if R0 < 1, then s(J|E0 ) < 0, the DFE point Z0 of system (1) is locally asymptotically stable.
180
4.1
178
181 182
183 184
Global stability analysis of DFE
Theorem 2. For R0 <1, the DFE Z0 of the system (1), is stable globally asymptotically on Ω and unstable for R0 >1. Proof. To establish the global stability of the disease-free equilibrium, we construct the following Lyapunov function: L(t) = A1 L + A2 I + A3 T,
185 186
where Ai , for i = 1,2,3, are some positive constants to be chosen later. Calculating the time derivative of L(t) along the solutions of system (1), we obtain: βSI dL(t) − k1 L + (1 − η)δT + A2 ǫL + ηδT − k2 I + A3 γI − k3 T = A1 dt N ≤ A1 βI − k1 L + (1 − η)δT + A2 ǫL + ηδT − k2 I + A3 γI − k3 T , S ≤ N = A1 β + A3 γ − A2 k2 I + [A2 ǫ − A1 k1 ]L + [A1 (1 − η)δ + A2 ηδ − A3 k3 ]T A1 β + A3 γ ≤ A2 k2 − 1 I + [A2 ǫ − A1 k1 ]L + [A1 (1 − η)δ + A2 ηδ − A3 k3 ]T. A2 k2 Now choosing A1 = ǫ, A2 = k1 and k3 =
187
ǫ(1−η)δ+k1 ηδ k3
and after simplification we get
dL(t) ≤ k1 k2 R0 − 1 I. dt
(3)
190
Hence if R0 < 1 then dL(t) is negative. Therefore the largest compact invariant set in dt Ω is the singleton set Z0 . Hence LaSalle’s invariant principle [8] then implies that Z0 is globally asymptotically stable in Ω.
191
5
192
Lemma 3. A unique positive endemic equilibrium (EE) exists when R0 exceeds 1.
188 189
Existence and stability of endemic equilibrium
9
193 194
Proof. The EE of (1) at E1 = (S ∗ , L∗ , I ∗ , T ∗) is obtained by setting the right hand side of the system (1) equal to zero, we get ∗ ∗ S = N , R 0 ((µ+σ1 )+γ(µ+σ2 +α)+γδ(1−η))N ∗ (R0 −1) ∗ L = R0 {(µ+σ1 )+γ(µ+σ2 +α)+γδ(1−η)+ǫ(γ+k3 )} , I∗ = ∗ T =
ǫk3 N ∗ (R0 −1) , R0 {(µ+σ1 )+γ(µ+σ2 +α)+γδ(1−η)+ǫ(γ+k3 )} ǫγN ∗ (R0 −1) . R0 {(µ+σ1 )+γ(µ+σ2 +α)+γδ(1−η)+ǫ(γ+k3 )}
195
Clearly, from above we conclude that a unique EE exists whenever R0 > 1.
196
5.1
Local stability of EE
197
198 199
Theorem 3. For R0 > 1, then the EE of the model (1) at E1 is locally asymptotically stable. Proof. The Jacobian matrix JE1 of the system (1) around the EE E1 is obtained by
J(E1 ) = 200
−I
∗ (N ∗ −S ∗ )β
−µ
N ∗2 I ∗ (N ∗ −S ∗ )β N ∗2
0 0
I ∗S∗ β N ∗2 I ∗S∗β − N ∗2 −
k1
− (N
ǫ 0
∗ −I ∗ )S ∗ β
N ∗2 (N ∗ −I ∗ )S ∗ β N ∗2 −k2
γ
I ∗ S∗β N ∗2
δ(1 − η) − δη −k3
I ∗S∗β N ∗2
.
The characteristic equation associated to J(E1 ) is λ4 + b1 λ3 + b2 λ2 + b3 λ + b4 = 0,
201
(4)
where the coefficients are
I ∗ βS ∗ b1 = k1 + k2 + k3 + µ + , N ∗2 1 b2 = k3 (µN ∗ 2 + IβS ∗ ) + QN ∗ 2 + k2 (µN ∗ 2 + IβS ∗ ) + k1 ((k2 + k3 + µ)N ∗ 2 + I ∗ β(N ∗ − S ∗ )) N ∗2 ∗ ∗ ∗ ∗ +βS (I µ + (N − I )ǫ) , 1 b3 = k1 (k3 (µN ∗ 2 + I ∗ β(N ∗ − S ∗ )) + QN ∗ 2 + k2 (µN ∗ 2 + I ∗ β(N ∗ − S ∗ ))) + k2 (k3 (µN ∗ 2 + I ∗ βS ∗ ) N2 +I ∗ βµS ∗ ) + βk3S ∗ I ∗ + βS ∗ I ∗ (γ + µ(1 + I ∗ N ∗ )ǫ) − (βk3 S(N ∗ I ∗ )ǫ + βγδη ∗2 +γδN ((1 − η)ǫ + ηµ)) , 10
b4
202
1 ∗2 = + I ∗ β(N ∗ − S ∗ )) + βS ∗ I ∗ (k2 k3 + ǫγµ) − (βk3µS ∗ ǫ(N ∗ − I ∗ ) 2 k1 Q(µN ∗ N ∗2 ∗ ∗ ∗ ∗ ∗ +γµ(δ(1 − η)ǫN + βδηS I ) + βγδI (1 − η)ǫ(N − S )) ,
where Q = {(µ + σ1 )k3 + γ(µ + σ2 + α) + γδ(1 − η)}.
203 204 205
206
207 208 209
The the eigenvalues of the characteristic polynomial (4) will be negative if it fulfil the Routh-Hurtwiz conditions that are bi > 0 for i = 1, 2, 3, 4 and b1 b2 b3 > b23 + b21 b4 . Hence then the EE of the model (1) will be locally asymptotically stable if R0 > 1.
5.2
Backward bifurcation analysis
Here, we discuss the existence of backward bifurcation of the system (1). To analyze this, we follow the center manifold theory described in [3]. If we take β as a bifurcation parameter, then at R0 = 1, we get β∗ = β =
210 211 212
k1 k2 k3 − γδ(1 − η)ǫ − γδηk1 . k3 ǫ
The following variations are made in the variables of the system (1) so that S = x1 , L = x2 , I = x3 , and T = x4 . Further using vector notation x = (x1 , x2 , x3 , x4 )T the model (1) can then be written in the form dx = F , with F = (f1 , f2 , f3 , f4 )T as shown below dt dx1 =Λ− dt dx2 = dt
βx1 x3 N
βx1 x3 N
− µx2 ,
− k1 x2 + (1 − η)δx4 ,
dx3 = ǫx2 + ηδx4 − k2 x3 , dt dx4 = γx3 − k3 x4 , dt 213
(5)
where N = Σ4i=1 xi . The Jacobian matrix evaluated at disease-free equilibrium Z0 with β ∗ is 1 +k1 k2 k3 −µ 0 − −γδǫ(1−η)−γδηk 0 ǫk3 1 +k1 k2 k3 0 −k1 −γδǫ(1−η)−γδηk δ(1 − η) ǫk 3 J= . 0 ǫ −k2 δη 0 0 γ −k3 11
(6)
214 215 216
The Jacobian matrix J has a simple zero eigenvalue calculated at β ∗ . The right and left eigenvectors denoted by W = (w1 , w2 , w3 , w4 ) and V = (v1 , v2 , v3 , v4 ) respectively and obtained w1 = −
w3 (γδηǫ − γδǫ − γδηk1 + k1 k2 k3 ) , k3 µǫ
217
w2 = 218
w3 (k2 k3 − γδη) γw3 > 0, w3 > 0, w4 = , k3 ǫ k3
and
k1 v2 δv2 (ηk1 + ǫ(1 − η)) , v4 = . ǫ k3 ǫ Computation of a: To evaluate a, we use the following partial derivatives v1 = 0, v2 > 0, v3 =
219
−2βµ ∂ 2 f2 ∂ 2 f2 −βµ ∂ 2 f2 = , = = , 2 ∂x3 Λ ∂x2 ∂x3 ∂x3 ∂x4 Λ 220
we get −2βµ 2 a= v2 w2 w3 + v2 w3 + v2 w3 w4 . Λ
221
(7)
Computation of b: We find the following partial derivatives to evaluate b ∂ 2 f2 = 1, ∂x3 ∂β
222
we obtain b = v2 w3 > 0.
(8)
224
Since all the terms in the expressions of b are positive, therefore, b > 0. Hence, the presence of the bifurcation at β = β ∗ depends solely on the sign of a.
225
5.3
223
226 227 228
229
Global stability of endemic equilibrium (EE)
In this subsection we prove the global asymptotic stability property of the EE of the proposed model (1). At the EE, we have from system (1) by using the approach given in [19, 20] ∗ ∗ Λ = βSN ∗I + µS ∗ , ∗ ∗ k1 L∗ = βSN ∗I + (1 − η)δT ∗ , (9) k2 I ∗ = ǫL∗ + ηδT ∗ , ∗ γI = k3 T ∗ . Theorem 4. If R0 > 1, then the EE of the model (1) is globally asymptotically stable. 12
230
231
232
Proof. To show the result, we define the following Lyapunove function S(t) I(t) ∗ ∗ ∗ ∗ K(t) = ǫ S(t) − S − S ln ∗ + k1 I(t) − I − I ln ∗ S I ∗ T (t) δT (k1 η + ǫ(1 − η)) + T (t) − T ∗ − T ∗ ln ∗ . ∗ γI T Taking time derivative of (10) along the model solution, we get S∗ ′ L∗ ′ I∗ ′ dK(t) = ǫ (1 − )S + (1 − )L + k1 (1 − )I dt S L I ∗ ∗ T δT (k1 η + ǫ(1 − η)) + (1 − )T ′ . ∗ γI T
(10)
(11)
It follows from direct calculation: i S ∗ h βSI S∗ ′ S = ǫ 1− Λ− − µS ǫ 1− S S N i S ∗ h βS ∗ I ∗ βSI ∗ = ǫ 1− + µS − − µS S N∗ N S∗ S ǫβS ∗ I ∗ S ∗ IN ∗ ∗ = ǫµS 2 − 1− − ∗ + + ∗ S S N∗ S I N ǫβSI . − N
(12)
i L∗ ′ L∗ h βSI ǫ 1− L = ǫ 1− − k1 L + (1 − η)δT L L N ǫβSI ǫβSIL∗ T L∗ = − − k1 ǫL + k1 ǫL∗ + (1 − η)ǫδT − (1 − η)δǫ N NL L ∗ ∗ ∗ ǫβSI ǫβSIL ǫβS I = − − k1 ǫL + + (1 − η)δǫT ∗ + (1 − η)ǫδT N NL N∗ T L∗ −(1 − η)δǫ . (13) L I∗ ′ I ∗ h k1 1 − I = k1 1 − ǫL + ηδT − k2 I I I I∗ I∗ = k1 ǫL − k1 ǫL + k1 ηδT − k1 ηδT − k1 k2 I + k1 k2 I ∗ I I ∗ I L I ∗ ǫβS ∗ I ∗ = k1 ǫL − k1 ǫL∗ ∗ + k1 ηδT − k1 ηδT + + k1 ηδT ∗ + ǫ(1 − η)δT ∗ IL I N∗ ǫβS ∗ I ∗ I I I − − (1 − η)ǫδT ∗ ∗ − k1 ηδT ∗ ∗ ∗ ∗ N I I I
13
∗ I ∗ ǫβS ∗ I ∗ ǫβS ∗ I ∗ LI ∗ ∗ LI − (1 − η)ǫδT + k ηδT − k ηδT + + k1 ηδT ∗ + 1 1 N ∗ L∗ I L∗ I I N∗ ǫβS ∗ I ∗ I ∗ I ∗ I ǫ(1 − η)δT ∗ − − (1 − η)ǫδT − k ηδT . (14) 1 N∗ I∗ I∗ I∗
= k1 ǫL −
233
T∗ δT ∗ (k1 η + ǫ(1 − η)) T∗ I∗ δT ∗ (k1 η + ǫ(1 − η)) 1 − [γI − k T ] = 1 − [I − T] 3 γI ∗ T I∗ T T∗ I I IT ∗ = δηk1 T ∗ ∗ + δ(1 − η)ǫT ∗ ∗ − δηk1 T ∗ ∗ I I I T ∗ IT −δ(1 − η)ǫT ∗ ∗ − δηk1 T − δ(1 − η)ǫT I T +δηk1 T ∗ + δ(1 − η)ǫT ∗ . 234
235
236 237 238
239
240 241 242 243 244
(15)
Using equations (12) to (15) in (11) we get dK(t) S∗ S ǫβS ∗ I ∗ S∗ I LI ∗ SIL∗ N ∗ LS ∗ ∗ = ǫµS 2 − − ∗ + 3 − − − − (1 − ) dt S S N∗ S I ∗ L∗ I S ∗ I ∗ NL L∗ S I ∗ L L∗ T T ∗I ∗ +(1 − η)ǫδT 3 − ∗ − − IL LT ∗ T I ∗ I ∗T IT ∗ ∗ − . (16) +δηk1 T 2 − IT ∗ I ∗ T Using properties of the arithmetic mean and geometric mean in equation (16) we have S S∗ 2 − − ≤ 0, S∗ S S∗ I LI ∗ SIL∗ N ∗ LS ∗ 3 − − − − (1 − ) ≤ 0, S I∗ L∗ I S∗I ∗N L L∗ S ∗ ∗ ∗ 3 − IILL∗ − LLTT∗ − TT I I∗ ≤ 0, I∗T IT ∗ 2 − IT ∗ − I ∗ T ≤ 0.
≤ 0 when R0 >1. As, all the parameters are non-negative, therefore it follows that dK(t) dt ∗ Therefore, by the LaSalles Invariance Principle [8], (S, L, I, T ) → (S , L∗ , I ∗ , T ∗ ) as t → ∞.
6
Sensitivity analysis of R0
Here, we analyzed the effect of the various model parameters to define the relative importance of these parameters to disease transmission. For this purpose, we calculate the partial derivatives of threshold number R0 with respect to model parameters γ, δ, β and 0 η. Parameter η is the failure rate of treatment and ∂R > 0. Therefore, the total number ∂η of infected population can be decreased by reducing the treatment failure rate η. Also, 14
245 246 247 248
∂R0 ∂δ
0 > 0, therefore, the parameter δ have positive effect. Since ∂R > 0 which implies that ∂β by decreasing the transmission rate β, the infection can be reduced. On the other hand, 0 < 0. It means that by increasing the parameter γ, the since the partial derivatives ∂R ∂γ TB infection can be controlled.
Parameter β α γ δ η ǫ µ σ1 σ2
Description Disease contact rate Progression rate from T to R Transmission rate from I to T Rate at which treated individuals leave the T Treatment failure rate Rate of moving from L to I Natural death rate Death rate in I Death rate in T
Sensitivity +0.7221 -0.066192 -0.272156 +0.0778308 +0.253903 +0.0655175 -0.0962639 -0.42162 -0.0091748
Table 3: Sensitivity indices of R0 with respect to model parameters. 249 250
Definition 1. The normalized sensitivity index of R0 depending on the differentiability on a parameter Ψ is defined as follows:
260
Ψ ∂R0 . (17) R0 ∂Ψ We calculate the sensitivity indices of the basic reproduction number R0 with respect to model parameters using above formula given in (17) and are tabulated in Table 3. From Table 3 we can see that parameters β, δ, η and ǫ have a positive sign which means that R0 increases with the parameter. While the rest of parameters have negative sign means that R0 decreases for higher values of the parameters. Further, this implies that by increasing (or decreasing) of these parameters by 10% will increase (or decrease) the basic reproductive number by 7.221%, 0.7783%, 2.539% and 0.6551% respectively. On the other hand the parameters α, γ, µ, σ1 , and σ2 have negative influence hence, increasing these rates by 10% will decrease the basic reproduction number by 0.6819%, 2.721 %, 0.9621%, 4.2162% and 0.0917%.
261
7
SΨR0 =
251 252 253 254 255 256 257 258 259
262 263 264 265 266 267 268
Numerical results and Discussion
To illustrate the numerical results, the following values were taken for the initial conditions to mimic the Khyber Pakhtunkhwa of Pakistan. According to the Pakistan Bureau of Statistics the population of Khyber Pakhtunkhwa in the 2017 census is estimated as 30,523,371 [4]. We take the total initial population as N(0) = 30, 523, 371 . We assume that initial population of exposed individuals L(0) = 83000. The initial infected is as given in the data is I(0) = 8010. We assume there are no treated and recovered individuals initially, so T (0) = 0 and R(0) = 0 respectively. Hence, the initial susceptible population 15
269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285
286
287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310
is given as S(0) = N(0) − L(0) − I(0) − T (0) − R(0) = 30432361. The numerical values parameters fitted from real data are presented in the Table 2. Fig. 3, shows the simulation of the model (1), when R0 = 0.8068 < 1. It is clear that all solutions converged to the DFE, Zo . In Fig. 4, the simulation of the model (1) is shown when R0 = 1.3419 > 1. Clearly all solutions converged to the EE, E1 . Fig. 5, shows that by taking different initial values all infected individuals vanishes when R0 = 0.8068 < 1, β = 0.2433. In Fig. 6, it is clear that when R0 = 1.3419 > 1 and β = 0.5433 and by considering diffract initial conditions then all infected classes converged to the EE, E1 . Fig. 7, explains the effect of treatment rate γ on the number of cumulative new cases of infection. Clearly from Fig. 7, by increasing treatment parameter γ the total TB infective individuals decreases. In Fig. 8, the effect of treatment failure parameter η is shown. The total infected individuals decreases when η decreases. The effect of parameter δ is simulated in Fig. 9, describing that by decreasing values of δ the total infected population decreases. The behavior of TB infected cases with and without treatment when R0 = 0.8068 < 1 and β = 0.2433, is given in Fig. 10. Finally, in Fig. 11, we describes the dynamics of infected individuals with and without treatment for R0 = 1.3419 > 1 and β = 0.5433. Figs. 13 to 18 show the behavior of the sensitive parameters to the basic reproduction number R0 .
8
Conclusion
In this paper, we formulated and analyzed the transmission dynamics of TB in Khyber Pakhtunkhwa province in Pakistan. The model was parameterized by using the data since 2002-2017 obtained from NTP. The stability results for the model are obtained and found stable locally as well globally when the basic reproduction number R0 less or greater than unity respectively. The numerical results for various values were presented. The government of Pakistan initiated various programs to eliminate TB infection in the last two decades. We developed the TB incidence model to simulate the TB data of Khyber Pakhtunkhwa reported by NTP and WHO [11, 1], from 2002-2017. The value of the basic reproduction number R0 and biological parameters of the model were estimated on the basis of available data and are tabulated in Table 2. To eliminate the TB epidemics it is necessary to ensure that R0 < 1. From numerical simulations, it is clear that the basic reproduction number R0 depends on the transmission rate β, treatment rate γ, treatment failure parameter η and on the proportion δ, at which the treated individuals leave the treatment class. Therefore, to control and eradicate the TB infection, it is important to consider the following strategies. The first and most important strategy is to minimize the contact with TB infected individuals by decreasing the value of β. Second important strategy is to increase the treatment rate γ of infective individuals. Also, to decrease the TB infective we need to reduce the treatment failure cases. Although the most effective strategy is to isolate infective by decreasing the value of transmission coefficient β. However, in developing countries like Pakistan it is very difficult to isolate the infectious individuals because of the high cost of the long term treatment. Therefore, we proposed the most feasible and optimal strategy to eliminate the TB infection is to increase the treatment rate by decrease the treatment cost so that the poor people can get the treatment measure, and by proper follow up of under treatment population so that the failure of treatment 16
9
×10 4 L I T R
8
Behavior of individuas
7 6 5 4 3 2 1 0 0
50
100
150
200
250
300
350
400
450
Time(Years)
Figure 3: The plot shows the behavior of individuals when R0 = 0.8068 < 1. 311
cases decreases.
312
Conflict of interests
313
The author’s declare that no competing interests exists.
314
Funding
315
No source of funding for this study.
17
3.5
×10 7 S L I T R
3
Endemic Equiibrium
2.5
2
1.5
1
0.5
0 0
100
200
300
400
500
600
700
800
900
1000
Time(Years)
Figure 4: The plot shows the behavior of individuals when R0 = 1.3419 > 1.
10
×10 4
Tota number of infected cases
9 8 7 6 5 4 3 2 1 0 0
50
100
150
200
250
300
350
400
450
Time(Years)
Figure 5: The plot shows the total number of infected individuals with different initial conditions with R0 = 0.8068 < 1 and β = 0.2433. 18
8
×10 6
Total number of infected cases
7
6
5
4
3
2
1
0 0
50
100
150
200
250
300
Time(Years)
Figure 6: The plot shows the total number of infected individuals with different initial conditions with R0 = 1.3419 > 1 and β = 0.5433.
8
×10 6 γ=0.2873 γ=0.3873 γ=0.4873 γ=0.5873 γ=0.6873 γ=0.7873
Tota number of infected cases
7
6
5
4
3
2
1
0 0
50
100
150
200
250
300
Time(Years)
Figure 7: The plot shows the total number of infected individuals by changing the value of γ.
19
8
×10 6 η=0.15 η=0.1 η=0.01
Total number of infected cases
7
6
5
4
3
2
1
0 0
50
100
150
200
250
300
Time(Years)
Figure 8: The plot shows the total number of infected individuals for different values of η.
8
×10 6 δ=1.1996 δ=1.0 δ=0.8 δ=0.6 δ=0.4 δ=0.2
Total number of infected cases
7
6
5
4
3
2
1
0 0
50
100
150
200
250
300
Time(Years)
Figure 9: The plot shows the total number of infected individuals for different values of δ.
20
10
×10 4 With treatment Without treatment
Total number of infected cases
9 8 7 6 5 4 3 2 1 0 0
20
40
60
80
100
120
140
160
180
200
Time(Years)
Figure 10: The plot shows the comparison of the total number of infected individuals with and without treatment with R0 = 0.8068 < 1 and β = 0.2433.
10
×10 6 With treatment Without treatment
Total number of infected cases
9 8 7 6 5 4 3 2 1 0 0
50
100
150
200
250
300
Time(Years)
Figure 11: The plot shows the comparison of the total number of infected individuals with and without treatment with R0 = 1.3419 > 1 and β = 0.5433. 21
1
2.5
1
1.2
1.4
2 1.8
1.6
2
0.9
1.8
0.8 0.7
1.5
0.6
1
0.5
δ
R0
2
1.6
1.4
0.8
1.2
1.4
2 1.8 1.6
0.3
0 1
1.2
1
0.4
0.5
1
0.2 1 0.8
0.5
8 0.
0.1
0.6
0.8
0.6
0.4
δ
0
0
0.2 0
0.6
0
γ
0.2
0.4
0.6
0.8
1
γ
(a)
(b)
Figure 12: Effects of parameters γ and δ on the basic reproduction number R0 .
1
1
1.2
1.4
1.9
0.8
1.8
0.7
1.7
0.6
1.6
0.5
1.5
0.4
1.4
1
1.2
1.4
0.3
0.5 1
1.6
1
2 1.8
1.5
η
R0
2
2 1.8
0.9 2.5
1.6
2
1.3
0.2
1.2
0.1
1.1
1 0.8
0.5
0.6 0.4
η
0
0
0.2 0
1
0
γ
0.2
0.4
0.6
0.8
1
γ
(a)
(b)
Figure 13: Effects of parameters γ and η on the basic reproduction number R0 .
22
1
0.8
1.8
0.7
1.7
0.6
1.6
1.5
α
R0
1
1.9
1.2
1.4
2
1.6
2.5
2
1.8 2
0.9
0.5
1.5
1 1
0.4
1.2
1.4
1.6
0.5 1
1.4
1.8 2
0.3
1.3
0.2
1.2
0.1
1.1
1 0.8
0.5
0.6 0.4
α
0
0
0.2 0
1
0
γ
0.2
0.4
0.6
0.8
1
γ
(a)
(b)
Figure 14: Sensitivity analysis of the basic reproduction number R0 verses parameters γ and α.
1
2
0.9 2.5
1.8
0.8 1.6
η
R0
2
1.8
1.6
1.4
1.2
1
1
0.6
0.8
1.5
0.6
0.7
0.4
2
1.4
0.5
1.2
0.4
0.5
1
0.3
0 1
0.8
0.2 0.6
2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.1
0.6
0.4
1 0.8
0.5 0.4
η
0
0
0.2 0
0.4
0
β
0.2
0.4
0.6
0.8
1
β
(a)
(b)
Figure 15: Effects of parameters β and η on the basic reproduction number R0 .
23
2
3.5
1.8
3.5
3.5
3
1.6 3
3
1.4
2.5
2.5
δ
1.2 1
2.5 2
2
2
0.8 1.5
0.6 0.4
1.5
1.5
1
1
1
0.2 0.5
0.5
0 0
0.5
1
0.5
1.5
2
β
(a)
(b)
Figure 16: Effects of parameters β and δ on the basic reproduction number R0 .
1
3.5
0.9 3
4
0.8 0.7
2.5
1
3
2
1
2
5
0.5
1.
γ
0.4
0.5
0.2 1 0.8
1
0.1
0.6
1
0.5 0.4
γ
0
0
0
β
2.5
2
0
0.2
1.5
2
0.3
0 1
1. 5
R0
0.6
0.2
0.4
0.6
3 0.8
0.5
1
β
(a)
(b)
Figure 17: Effects of parameters β and γ on the basic reproduction number R0 . 316
317 318 319
320 321
322 323
References [1] World Health Organization Media Centre. Available 2018. https://www.who.int/en/news-room/fact-sheets/detail/tuberculosis. accessed 2018. [2] Carlos Castillo-Chavez and Zhilan Feng. To treat or not to treat: the case of tuberculosis. Journal of mathematical biology, 35(6):629–656, 1997. [3] Carlos Castillo-Chavez and Baojun Song. Dynamical models of tuberculosis and their applications. Mathematical biosciences and engineering, 1(2):361–404, 2004. 24
324 325
326 327
328 329
[4] World Health Organization. WHO country cooperation strategic. http://apps.who.int/iris/bitstream/10665/136607/1/ccsbrief pak en.pdf, 2016. [5] Z Feng and C Castillo-Chavez. Mathematical models for the disease dynamics of tubeculosis, 1998. [6] David J Gerberry. Practical aspects of backward bifurcation in a mathematical model for tuberculosis. Journal of theoretical biology, 388:15–36, 2016.
332
[7] Soyoung Kim, A Aurelio, and Eunok Jung. Mathematical model and intervention strategies for mitigating tuberculosis in the philippines. Journal of theoretical biology, 443:100–112, 2018.
333
[8] Joseph P LaSalle. The stability of dynamical systems, volume 25. Siam, 1976.
330 331
334 335
336 337
338 339
340 341
342 343 344
345 346 347
348 349
350 351 352
353 354 355
356 357
[9] Junli Liu and Tailei Zhang. Global stability for a tuberculosis model. Mathematical and Computer Modelling, 54(1-2):836–845, 2011. [10] Luju Liu, Xiao-Qiang Zhao, and Yicang Zhou. A tuberculosis model with seasonality. Bulletin of Mathematical Biology, 72(4):931–952, 2010. [11] National TB Control Program http://www.ntp.gov.pk/webdatabase.php.
Pakistan
(NTP).
[12] Pakistan Bureau of Statistics. Pakistans 6th census: Population of major cities 583 census. 584. http://www.pbscensus.gov.pk/. [13] Charles S Revelle, Walter R Lynn, and Floyd Feldmann. Mathematical models for the economic allocation of tuberculosis control activities in developing nations. American review of respiratory disease, 96(5):893–909, 1967. [14] James M Trauer, Justin T Denholm, and Emma S McBryde. Construction of a mathematical model for tuberculosis transmission in highly endemic regions of the asia-pacific. Journal of theoretical biology, 358:74–84, 2014. [15] Saif Ullah, Muhammad Altaf Khan, and Muhammad Farooq. A fractional model for the dynamics of tb virus. Chaos, Solitons & Fractals, 116:63–71, 2018. [16] Pauline Van den Driessche and James Watmough. Reproduction numbers and subthreshold endemic equilibria for compartmental models of disease transmission. Mathematical biosciences, 180(1-2):29–48, 2002. [17] Hans Waaler, Anton Geser, and Stig Andersen. The use of mathematical models in the study of the epidemiology of tuberculosis. American Journal of Public Health and the Nations Health, 52(6):1002–1013, 1962. [18] Robert S Wallis. Mathematical models of tuberculosis reactivation and relapse. Frontiers in microbiology, 7:669, 2016.
25
358 359
360 361
362 363 364
365 366
[19] Yi Wang and Jinde Cao. Global dynamics of multi-group SEI animal disease models with indirect transmission. Chaos, Solitons & Fractals, 69:81–89, 2014. [20] Yi Wang and Jinde Cao. Global stability of general cholera models with nonlinear incidence and removal rates. Journal of the Franklin Institute, 352(6):2464–2485, 2015. [21] Yali Yang, Jianquan Li, Zhien Ma, and Luju Liu. Global stability of two models with incomplete treatment for tuberculosis. Chaos, Solitons & Fractals, 43(1-12):79–85, 2010. [22] Jinhui Zhang, Yong Li, and Xinan Zhang. Mathematical modeling of tuberculosis data of china. Journal of theoretical biology, 365:159–163, 2015.
26