- Email: [email protected]
Optimal control strategies for dengue transmission in pakistan
Accepted Manuscript
Optimal Control Strategies for Dengue Transmission in Pakistan F.B. Agusto, M.A. Khan PII: DOI: Reference:
S0025-5564(18)30453-X https://doi.org/10.1016/j.mbs.2018.09.007 MBS 8125
To appear in:
Mathematical Biosciences
Received date: Revised date: Accepted date:
29 July 2018 5 September 2018 6 September 2018
Please cite this article as: F.B. Agusto, M.A. Khan, Optimal Control Strategies for Dengue Transmission in Pakistan, Mathematical Biosciences (2018), doi: https://doi.org/10.1016/j.mbs.2018.09.007
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.
ACCEPTED MANUSCRIPT
Highlights • A mathematical model of dengue outbreak in Peshawar, Pakistan is considered. • Model formulations and its basic properties are presented.
CR IP T
• Sensitivity analysis and optimal control problem is formulated and discussed
AC
CE
PT
ED
M
AN US
• Control strategies are presented for disease elimination.
1
ACCEPTED MANUSCRIPT
CR IP T
Optimal Control Strategies for Dengue Transmission in Pakistan F. B. Agusto1 , M.A. Khan2∗ 2 1
Department of Ecology and Evolutionary Biology, University of Kansas, USA.
Department of Mathematics, City University of Science and Information Technology, 25000, KP, Pakistan.
6 7 8 9 10 11 12 13 14
15 16
17
18 19 20 21 22 23
M
5
ED
4
PT
3
This paper presents a deterministic model for dengue virus transmission. The model is parameterized using data from the 2017 dengue outbreak in Pakistan. We estimated the basic reproduction number (R0 ) without any interventions for the 2017 dengue outbreak in Peshawar district of Pakistan as R0 ≈ 2.64, the distribution of the reproduction number lies in the range R0 ∈ [1.21, 5.24] (with a mean R0 ≈ 2.64). Optimal control theory is then applied to investigate the optimal strategy for curtailing the spread of the disease using two time-dependent control variables determined from sensitivity analysis. These control variables are insecticide use and vaccination. The results show that the two controls avert the same number of infections in the district regardless of the weights on the costs this is due to the reciprocal relationship between the cost of insecticide use and vaccination. A strong reciprocal relationship exists between the use of insecticide and vaccination; as the cost of insecticide increases the use of vaccination increases. The use of insecticide on the other hand slightly increases when vaccination level decreases due to increase in cost.
CE
2
Abstract
Key Words: Dengue model, Vaccination, Basic reproduction number, Parameters estimation, Sensitivity analysis, optimal control strategies, numerical results.
AC
1
AN US
September 12, 2018
1
Introduction
Dengue fever is a mosquito-borne disease caused by the dengue virus. It has four different serotypes, i.e., DEN1, DEN2, DEN3, and DEN4, of the genus Flavivirus. The dengue is endemic in almost of 100 countries of the world, such as Africa, Eastern Mediterranean, the Americas, the subtropical regions, and especially the South-East Asia and Western Pacific regions are the most seriously affected [23, 55, 65]. After malaria, dengue has been recognized as the next deadliest mosquito-borne disease, with more than 100 millions ∗
Corresponding Author email: [email protected]
ACCEPTED MANUSCRIPT
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
CR IP T
30
The dengue virus transmission to the human is by the bites of mosquitoes (female mosquitoes, of type Aedes aegypti) [10], acquiring infection by taking the blood meal from an infected individual, then the infected mosquitoes pass the virus to healthy individuals (susceptible individuals). The recovery of an individual from one serotype leads to the permanently immune to it and temporarily or partially-immune to the others serotypes [63].
AN US
29
There is no specific effective treatment for the dengue fever until now but only the fluid replacement therapy at the beginning of infection diagnosis, and some traditional medicine (Brazilian traditional medical) in the dengue is treated with cat’s claw herb, but no such treatment observed [22, 24]. Besides treatment, there is no such effective and safe vaccine available for the dengue fever in the market. A report WHO shows some recent advancement about the vaccine development for the dengue fever in [64]. The first dengue vaccine, Dengvaxia (CYD-TDV) by Sanofi Pasteur, was first registered in Mexico in December, 2015. CYD-TDV is a live recombinant tetravalent dengue vaccine that has been evaluated as a 3-dose series on a 0/6/12 month schedule in Phase III clinical studies. It has been registered for use in individuals 9-45 years of age living in endemic areas [64]. Besides this, undergoing vaccine development includes, including live attenuated mono- and tetra-valent formulation, inactivated whole virus vaccines, and recombinant subunit vaccines[14, 13]. Any future dengue vaccine that perfectly protect against all serotype may or may not be possible, thus it is expected that any future dengue vaccine would be imperfect [10].
M
28
ED
27
PT
26
infections and some reports estimates about 390 million [12] and 20,000 deaths [55]. A study from 2012, shows that 3.9 billion people in 128 countries are at risk of infection with dengue viruses [17]. It is spread by many types of mosquitoes of the type Aedes, principally A. aegypti. Break bone fever or classical dengue fever relatively causes both mild morbidity and mortality, and the sufferers recover within a short span of time of about one or two weeks from the onset of fever [35, 36], some people develop hemorrhagic fever (DHF) or dengue shock syndrome (DSS) [37]. The disease severity drastically increased with mortality ranging from (5-15%) [35, 36, 38]. World Health Organization (WHO) report shows the occurrence of high number of cases of DHF annually worldwide [13, 14].
A number of mathematical models on the transmission dynamics of dengue fever have been proposed and analyzed, see [11, 23, 25, 26, 29, 30, 31, 42, 50, 49, 52, 57, 58, 60, 70] and the references therein. For instant, Rodrigues et al. [50] formulated and analyzed a vaccination models and optimal control strategies to dengue. First they studied the vaccination process by incorporating a new compartment in the model, and investigated the different ways of distributing the vaccines: pediatric and random mass vaccines, with distinct levels of efficacy and durability. They found that dengue eradication success depends on the type of vaccine as well as on the vaccination coverage. Simulating both epidemic and endemic scenarios, they showed using vaccine as a control variable, that vaccine when available can be a promising strategy to fight the disease. Rodrigues et al. [49] applied optimal control theory to study Dengue epidemics. Their proposed cost functional depends on the costs of medical treatment of the infected people, and the costs related to educational and sanitation campaigns. They showed efficient ways to decrease the number of infected mosquitoes and individuals in less time and with lower costs using current computational tools. Shim [57] studied the dynamics of dengue and cost-effectiveness of vaccine in the
CE
25
AC
24
3
ACCEPTED MANUSCRIPT
68 69 70
71 72 73 74 75 76
Philippines using an age-structured model. Their results indicate that the higher the cost of vaccination, the less cost-effective the dengue vaccination program. Furthermore, their results demonstrated that at a relatively low vaccine efficacy an age-targeted vaccination may be cost-effective provided the cost is sufficiently low vaccination. Motivated by the above studies, we develop a mathematical model to investigate the transmission dynamics of Dengue fever in Pakistan. We aim to study the impact of an imperfect vaccine in the bid to control dengue in Pakistan. Therefore, incorporated into the model is an imperfect vaccine. We parameterize the model using the 2017 dengue outbreak from Peshawar district of Khyber Pakhtunhwa Province in Pakistan. We then apply optimal control theory to study the best control strategies to curtail the disease in the community.
CR IP T
67
84
2
85 86 87 88 89 90 91 92
Model Formulation
M
82
The model features both human and mosquito populations. The human population is divided into five different classes, namely, susceptible individuals, SH (t), vaccinated, VH (t), exposed, EH (t), infected individuals, IH (t), and recovered individuals, RH (t), at any time t. Thus, the total human population denoted by NH (t) is given as NH (t) = SH (t) + VH (t) + EH (t) + IH (t) + RH (t). The population of mosquitoes is similarly divided into three different classes, namely, susceptible mosquitoes, SM (t), exposed mosquitoes, EM (t), and infected mosquitoes, IM (t), at any time t. So that, the total population of vector (mosquitoes) denoted by NM (t) is given by NM (t) = SM (t) + EM (t) + IM (t).
ED
81
PT
80
CE
79
AC
78
AN US
83
The paper is organized as follows: A model formulation and the basic mathematical properties of the model are presented in Section 2. The parameterization of the model is presented in Section 3. This is followed by a sensitivity analysis in Section 4. Information glean from the sensitivity analysis are used to identify the control variables applied in Section 5 on optimal control problem. Numerical results with detailed discussions are presented in Section 6. Finally, a detailed discussion and conclusion is presented in Section 7.
77
4
ACCEPTED MANUSCRIPT
The model that govern the system of differential equations is given by,
94
where λH =
= νH SH − (1 − ε)λH VH − (ωH + µH )VH , = λH [SH + (1 − ε)VH ] − (σH + µH )EH , = σH EH − (γH + µH + δH )IH , = γH IH − µH RH , = ΛM − λM SM − µM SM , = λM SM − (σM + µM )EM , = σM EM − µM IM ,
bM (αH EM + βH IM ) bM (αM EH + βM IH ) , and λM = . NH NH
AC
CE
PT
ED
M
95
= ΛH + ωH VH − λH SH − (νH + µH )SH ,
CR IP T
dSH dt dVH dt dEH dt dIH dt dRH dt dSM dt dEM dt dIM dt
AN US
93
Figure 1: Flow diagram of the dengue model (1). 96
The variables and their parameters are defined in Table 1.
5
(1)
ACCEPTED MANUSCRIPT
97
Analysis of the model
98
Basic qualitative properties
99
Positivity and boundedness of solutions
102
103 104 105
CR IP T
101
For the dengue model (1) to be epidemiologically meaningful, it is important to prove that all its state variables are non-negative for all time. In other words, solution of the system (1) with non-negative initial data will remain non-negative for all time t > 0. Lemma 1. Let the initial data F (0) ≥ 0 , where F (t) = (SH (t), VH (t), EH (t), IH (t), RH (t), SM (t), EM (t), IM (t)). Then the solutions F (t) of the dengue model (1) are non-negative for all t > 0. Furthermore lim sup NH (t) ≤ t→∞
ΛM ΛH , and lim sup NM (t) ≤ µH µM t→∞
with
AN US
100
NH (t) = SH (t) + VH (t) + EH (t) + IH (t) + RH (t), and NM = SM (t) + EM (t) + IM (t).
ED
Description Population of susceptible individuals Population of vaccinated individuals Population of exposed individuals Population of infected individuals Population of recovered individuals Population of susceptible mosquitoes Population of exposed mosquitoes Population of infected mosquitoes Description Vaccination rate Vaccine efficacy Vaccine waning rate Human recruitment rate Transmission probability from EM to SH . Transmission probability from IM to SH Natural mortality rate in humans Disease progression rate of infectious of exposed humans Mosquito biting rate Recovery rate of infected individuals in human compartments IH Disease related death rate Mosquito recruitment rate Transmission probability from EH to SM Transmission probability from IH to SM Natural death rate of mosquitoes Disease progression rate of infectious of exposed mosquitoes
CE
PT
Variable SH VH EH IH RH SM EM IM Parameter νH ε ωH ΛH αH βH µH σH bM γH δH ΛM αM βM µM σM
M
The proof of Lemma 1 is given in Appendix A.
AC
106
Table 1: Description of the parameters and variables for the dengue model (1). 6
ACCEPTED MANUSCRIPT
107
108 109
Invariant regions The dengue model (1) will be analyzed in a biologically-feasible region as follows. Consider the feasible region Ω = ΩH × ΩM ⊂ R5+ × R3+ with,
ΛH 4 , ΩH = (SH (t), VH (t), EH (t), IH (t), RH (t)) ∈ R+ : NH (t) ≤ µH
and
ΩM
112
113
ΛM 3 = (SM (t), EM (t), IM (t)) ∈ R+ : NM (t) ≤ . µM
Lemma 2. The region Ω ⊂ R8+ is positively-invariant for the model (1) with non-negative initial conditions in R8+ .
AN US
111
CR IP T
110
The prove of Lemma 2 is given in Appendix B.
116
Stability of disease-free equilibrium (DFE)
M
115
In the next section, the conditions for the existence and stability of the equilibria of the model (1) are stated.
114
ED
The dengue model (1) has a disease free equilibrium (DFE). The DFE is obtained by setting the right-hand sides of the equations in the model (1) to zero, which is given by
120 121
k2 ΛH Λ H νH ΛM , , 0, 0, 0, , 0, 0 . k1 k2 − νH ωH k1 k2 − νH ωH µM
CE
118 119
where k1 = νH + µH and k2 = ωH + µH . The stability of E0 can be established using the next generation operator method on system (1). Taking EH (t), IH (t), EM (t), and IM (t) as the infected compartments and then using the notation in [61], the Jacobian F and V matrices for new infectious terms and the remaining transfer terms, respectively, are defined as:
AC
117
PT
0 0 E0 = (SH , VH0 , 0, 0, 0, SM , 0, 0) =
0 0 F = 0 bM αM SM N0 H 0
0 0 0 bM βM SM NH0 0
0 bM αH [SH + (1 − ε)VH0 ] NH0 0
0 0
7
0 bM βH [SH + (1 − ε)VH0 ] NH0 0 , 0 0
ACCEPTED MANUSCRIPT
(σH + µH ) 0 0 0 −σH (γH + µH + δH ) 0 0 . V = 0 0 (σM + µM ) 0 0 0 −σM µM
129
130 131 132 133 134 135 136
137 138 139 140 141 142 143
CR IP T
AN US
128
Lemma 3. The disease-free equilibrium (DFE) of the dengue model (1) is locally asymptotically stable (LAS) if R0 < 1 and unstable if R0 > 1.
Backward Bifurcation Analysis
Disease transmission models typically undergo a simple transcritical bifurcation at R0 = 1 where there is an exchange of model’s stability as the model equilibria moves from the DFE to an endemic equilibrium. Some models such as vaccination models, on the otherhand exhibit backward bifurcation, where the stable DFE co-exists with a stable endemic equilibrium when the reproduction number is less than unity (see [7, 4, 8, 18, 27, 28, 34, 54, 53]). In a backward bifurcation setting, disease control is only feasible if R0 is reduced further to values below another sub-threshold less than unity.
M
127
ED
126
The expression R0 is the number of secondary infections in completely susceptible population due to infections from one introduced infectious individual with dengue. Further, using Theorem 2 in [61], the following result is established.
PT
124 125
where, ρ is the spectral radius and k3 = σH + µH , k4 = γH + µH + δH , k5 = σM + µM .
The important implication of this phenomenon on public health is that the classical requirement of having the reproduction number less than unity, although necessary, is no longer sufficient for disease control. This implies that effective disease control is dependent on the initial sizes of the sub-populations of the model. It is instructive, therefore, to explore whether or not the model (1) exhibits the phenomenon of backward bifurcation. In determining this possibility in the model (1), we use the Centre Manifold theory [19], as described in Theorem 4.1 by Castillo-Chavez and Song [21].
CE
123
Therefore, the reproduction number is s 0 0 b2M SM (αH µM + βH σM )(αM k2 + βM σH )[SH + (1 − ε)VH0 ] , R0 = ρ(F V −1 ) = k3 k4 k5 µM NH0
AC
122
145
Theorem 1. The dengue model (1) undergoes a backward bifurcation at R0 = 1 whenever the inequality a > 0 holds.
146
The prove of Theorem 1 is given in Appendix C.
144
8
ACCEPTED MANUSCRIPT
147
3
148
3.1
156 157 158 159 160 161
162 163 164 165 166 167 168
169 170 171 172
CR IP T
155
Between January - September 2017, the Ministry of National Health Services, Regulations and Coordination in Pakistan started reporting a high number of dengue fever cases in seven provinces of the country. From July - August 2017, a total of 2,199 cases had been reported by the ministry, Khyber Pakhtunkhwa (KPK) province had the highest number of reported laboratory confirmed cases (1,279 cases), this was followed by Sindh province with 793 cases, AJK province had the least number of reported cases (4 cases) [68]. However, by November 2017, a total of 24,807 laboratory confirmed cases with 69 associated deaths were reported in Khyber Pakhtunkhwa province [69].
AN US
154
A demographic break down of the 1,279 number of cases reported in July - August 2017 from the KPK province [68] show the 25-64 age group having the highest number of reported cases (668 cases, about 52.22% of the total), this was followed by the 15-24 age group with 371 reported cases (29.04%). A total of 199 cases were reported amongst the 0-14 age group (15.53%); the least affected age group was the 65+ group, with 41 reported cases (3.2%). 62.9% of the cases (804) are male while females represented 37.1% (475) of the reported cases.
M
153
ED
152
In order to parametrize our model, we obtained from the dengue response unit in Peshawar the 2017 dengue outbreak data from various hospitals in the Peshawar district, namely, Lady Reading Hospital (LRH), Hayatabad Medical Complex (HMC), and Khyber Teaching Hospital (KTH), see Figure 2.
PT
151
The dengue is a frequent epidemic in Pakistan. In 2011 and 2013 dengue outbreak occurred in Lahore, Punjab and in Swat Khyber Pakhtunkha district respectively. The 2011 dengue outbreak in Lahore, Punjab lead to 300 death while the 2013 outbreak in Swat Khyber Pakhtunkha district caused a total of 57 death cases (see, [39, 40] and the references therein).
CE
150
The 2017 dengue outbreak data
AC
149
Parameter Estimation
9
ACCEPTED MANUSCRIPT
5000
0
Jul 24 17
Sep 12 17 Time (days)
CR IP T
10000
Nov 01 17
AN US
Reported numberof cases
15000
Figure 2: Number of new cases in the 2017 dengue outbreak in Pakistan. The data were obtained from Peshawar district hospitals.
175
In order to parametrize the dengue model (1), first, we study the worst-case scenario where public health interventions such as vaccination are not implemented in the community. In the absence of such interventions, the model (1) reduces to the following basic model: dSH dt dEH dt dIH dt dRH dt dSM dt dEM dt dIM dt
bM (αH EM + βH IM )SH − µH SH , NH bM (αH EM + βH IM )SH = − (σH + µH )EH , NH = ΛH −
AC
CE
PT
176
Pre-intervention model
M
174
3.2
ED
173
177
= σH EH − (γH + µH + δH )IH , = γH IH − µH RH , bM (αM EH + βM IH )SM − µM SM , NH bM (αM EH + βM IH )SM = − (σM + µM )EM , NH = ΛM −
= σM EM − µM IM .
The dengue model (2) is bounded and positively invariant in the region Ω = ΩH × ΩM ⊂ R4+ × R3+ 10
(2)
ACCEPTED MANUSCRIPT
178
with, ΩH = and
180 181 182 183
184 185 186
ΛH : NH (t) ≤ , µH
ΛM 3 = (SM (t), EM (t), IM (t)) ∈ R+ : NM (t) ≤ . µM
In Appendix D, we state the conditions for the stability of the disease-free equilibrium ˜ 0 . Further analysis about backward of model (2) based on the reproduction number, R bifurcation of model (2) are given in Appendix E. Backward bifurcation is the phenomenon in which two stable equilibria (in this case disease-free and endemic equilibria) co-exist when the reproduction number is less than one.
AN US
179
(SH (t), EH (t), IH (t), RH (t)) ∈
R4+
CR IP T
ΩM
A number of vector-borne diseases models [4, 5, 15, 59] have shown that in the absence of disease related death, the phenomenon of backward bifurcation can be lost and the disease-free equilibrium will become globally stable.
192
3.3
194 195 196 197 198 199 200
201 202 203 204 205
ED
193
Data fitting
In order to parametrize dengue model (2), we obtained some of the parameter values from literature (see Table 2), others were estimated or fitted to the Peshawar dengue outbreak data obtained from Peshawar district hospital using the classical least-squares method (see Figure 3). For instance, the demographic parameter, µH , is estimated as µH = 1/67.7 per year, where 67.7 years is the average lifespan in Pakistan [62]. The other demographic parameter, ΛH , is then estimated as follows. Since the total urban population of Peshawar district as at 2017 was 1, 970, 042 [47], we assumed that ΛH /µH = 1, 970, 042 is the limiting total human population in the absence of the disease, so that ΛH = 29, 099.59 per year.
PT
190
CE
189
AC
188
M
191
There is no specific treatment for dengue or severe dengue [67]; however with early detection and access to proper medical care fatality rates can be lower to rate below 1% [67] or may vary between 2%-5% [56]. But, when left untreated, the mortality rate can be as high as 20% [56]. Hence, in Appendix F, we ignore the disease induce death rate and show that the disease-free equilibrium corresponding to model (2) is globally stable.
187
Thus, we obtained after fitting the model (2) to the 2017 dengue outbreak data the following estimates for bM = 0.8759, αH = 0.8728, βH = 0.8535, σH = 0.01298, γH = 0.01837, αM = 0.8593, βM = 0.7646, σM = 0.005728, δH = 0.000015857. Consequently, ˜ 0 (stated in using these parameter estimates, we obtained, in this study, that the value of R ¯ 0 ≈ 2.65. Appendix D) for the 2017 dengue outbreak in Peshawar district of Pakistan is R
11
15000
Model Data
CR IP T
10000
5000
0
AN US
Cumulative number of infected cases
ACCEPTED MANUSCRIPT
Jul 24 17
Sep 12 17 Time (days)
Nov 01 17
M
Figure 3: Data fitting of the number of cases using model (1), for the 2017 dengue outbreaks in Pakistan. The parameters fitted are given as bM = 0.8759, αH = 0.8728, βH = 0.8535, σH = 0.01298, γH = 0.01837, αM = 0.8593, βM = 0.7646, σM = 0.005728, δH = 0.000015857. Baseline value 29,099.59 0.8728 0.8535 1/67.7 0.01298 0.8759 0.01837 0.000015857 0.3 0.8593 0.7646 1/42 0.005728
AC
CE
PT
ED
Parameter ΛH αH βH µH σH bM γH δH ΛM αM βM µM σM
Reference Estimated Fitted Fitted [62] Fitted Fitted Fitted Fitted [43] Fitted Fitted [43] Fitted
Table 2: Numerical values for the parameters of the dengue model (1). 206
207 208
3.4
Quantify the expected disease burden
In order to quantify the expected burden of the disease in the country (under the worst˜ 0 is generated, using the parameter values and case scenario), the distribution profile of R 12
ACCEPTED MANUSCRIPT
211 212 213 214
ranges in Table 2. The results obtained, depicted in Figure 4, show the distribution of the ˜ 0 ∈ [1.21, 5.24] (with a mean R ˜ 0 ≈ 2.71, suggesting reproduction number in the range R the potential for a larger dengue outbreaks, where one infected case infects, on average, about 2.71 others). Furthermore, it is worth noting that for the worse-case scenario, the ˜ 0 values are concentrated in the interval (2.2036, 3.1518), central 50% of the generated R ˜ with the median value of R0 = 2.64 been close to the mean.
CR IP T
209 210
5 4.5
3.5 3 2.5 2 1.5 200
300
400
500
600
700
800
900
1000
M
100
AN US
R0
4
Runs
ED
˜ 0 , for the models (2). Parameter Figure 4: The Scatterplots of the reproduction number, R values (baseline) and ranges used are as given in Table 2.
218
219 220
221
222 223 224 225 226 227 228
PT
217
CE
216
˜ 0 for the 2017 Hence, to summarize, we have obtained in this study, that the value of R ˜ dengue outbreak in Peshawar district of Pakistan is R0 ≈ 2.65 and the distribution of the ˜ 0 ∈ [1.2127, 5.2426] (with a mean R ˜ 0 ≈ 2.7052) for reproduction number is in the range R the estimated worse case scenario. In the next section, we apply optimal control theory to determine the optimal control strategy that will effectively curtail the spread of dengue in the community.
AC
215
4
Sensitivity analysis
The outputs of deterministic models are governed by the model input parameters, which may exhibit some uncertainty in their determination or selection. We employed a global sensitivity analysis to assess the impact of uncertainty and the sensitivity of the outcomes of the numerical simulations to variations in each parameter of the model (1) using Latin Hypercube Sampling (LHS) and partial rank correlation coefficients (PRCC). LHS is a stratified sampling without replacement technique which allows for an efficient analysis of parameter variations across simultaneous uncertainty ranges in each parameter [16, 44, 45, 13
ACCEPTED MANUSCRIPT
237 238 239 240 241 242 243
244 245 246 247 248 249 250 251
252 253 254 255 256 257 258
CR IP T
236
Figure 5 depict the PRCC values for each parameter of the model using the reproduction ˜ 0 and R0 , as the response functions. Parameters with the highest PRCC values numbers, R have the largest impact on the reproduction numbers. Therefore, the key parameters ˜ 0 and R0 ) are separated into those that decrease influencing there production numbers (R ˜ 0 and R0 when increased (bars extending to the left for negative PRCC values) and those R ˜ 0 and R0 to increase when increased (bars extending to the right for positive that causes R PRCC values).
AN US
235
From Figures 5(a) and 5(b), it follows that the parameters that have the most negative influence on dengue transmission dynamics are human recruitment rate (ΛH ), disease progression rate in humans (σH ), human recovery rate (γH ), mosquito natural death rate ˜0 (µM ), and vaccine efficacy (ε). While the parameters with the most positive impact on R and R0 are the transmission probability per contact of infectious mosquitoes with susceptible humans (αH and βH ), mosquito recruitment rate (ΛM ), mosquito biting rate (bM ), transmission probability per contact of susceptible mosquitoes with infectious humans (αM and βM ), disease progression rate in mosquitoes (σM ), and vaccine waning rate (ωH ).
M
233 234
ED
232
Identification of these key parameters is important to the formulation of effective control strategies necessary for combating the spread of disease in the community. In particular, the results of this sensitivity analysis suggest that a strategy that reduces the parameters with positive PRCC values (i.e., βH , αH , ΛM , bM , βM , αM , σH , ωH ) will adequately reduce the spread of dengue in the community. Furthermore, a strategy that increases the parameters with negative PRCC values (i.e., ΛH , σH , γH , µM , ε) will be effective in curtailing the spread of dengue in the community.
PT
231
51]. PRCC measures the strength of the relationship between the model outcome and the parameters, stating the degree of the effect that each parameter has on the outcome [16, 44, 45, 51]. Thus, sensitivity analysis determines the parameters with the most significant impact on the outcome of the numerical simulations of the model [16, 44, 46]. To generate the LHS matrices, we assume that all the model parameters are uniformly distributed. Then a total of 1,000 simulations of the models per LHS run were carried out, using the baseline values tabulated in Table 2 and the ranges as 20% from the baseline values (in either direction).
CE
230
AC
229
14
ǫ ωH ν µMH σM β αMM b ΛMM δ µHH γ σHH β αHH ΛH
µM σM βM αM bM ΛM δH µH γH σH βH αH ΛH
-1
-0.5
0
0.5
1
-1
-0.5
0
0.5
1
(b)
AN US
(a)
CR IP T
ACCEPTED MANUSCRIPT
Figure 5: PRCC values for the dengue models (2) and (1), using as response functions: (a) The
˜ 0 . (b) The reproduction number R0 . Parameter values (ranges) used are reproduction number R as given in Table 2.
263 264 265 266 267 268
269 270 271 272 273 274 275
276 277 278 279 280
M
ED
262
The result of the sensitivity analysis suggest that a control strategy that reduces the transmission probability per contact of infectious mosquitoes with susceptible humans (αH and βH ), mosquito recruitment rate (ΛM ), mosquito biting rate (bM ), transmission probability per contact of susceptible mosquitoes with infectious humans (αM and βM ), disease progression rate in mosquitoes (σM ), and vaccine waning rate (ωH ) will adequately reduce the spread of dengue in the community. Furthermore, a strategy that increases the human recruitment rate (ΛH ), disease progression rate in humans (σH ), human recovery rate (γH ), mosquito natural death rate (µM ), and vaccine efficacy (ε) will be effective in curtailing the spread of dengue in the community.
PT
261
The optimal control problems
CE
260
5
Therefore, we introduce into the transmission model (1) a time dependent control variable uV (t), representing the vaccination rate which was previously taken as a constant parameter (νH ). Although, this parameter does not significantly impact the reproduction number (R0 ), we nevertheless consider this parameter as time dependent variable since the vaccine efficacy (ε) have a significant impact of the (R0 ). We also consider as time dependent control variable uM (t) that impact the mosquito population since mosquito death rate µM has a strong significant impact on R0 .
AC
259
The mosquito biting rate bM and transmission probabilities αH , βH , αM and βM , are another set of candidate parameters to use as a time dependent controls. Using them as control measures however, will require the use of personal protection which is hard to implement since the mosquito species causing dengue only bite during the day. So, our control efforts using control variables uV (t) and uM (t) will indirectly impact these 15
ACCEPTED MANUSCRIPT
parameters. Thus, model (1) with the time dependent control variables uV (t) and uM (t) becomes:
283
bM (αH EM + βH IM )SH − (uV (t) + µH )SH NH (1 − ε)bM (αH EM + βH IM )VH = uV (t)SH − − (ωH + µH )VH NH bM (αH EM + βH IM )[SH + (1 − ε)VH ] = − (σH + µH )EH NH = ΛH + ωH VH −
= σH EH − (γH + µH + δH )IH = γH IH − µH RH ,
bM (αM EH + βM IH )SM − µM SM − uM (t)SM NH bM (αM EH + βM IH )SM = − (σM + µM )EM − uM (t)EM NH = ΛM −
= σM EM − µM IM − uM (t)IM .
Our goal is to minimize the cost function defined as Z
0
289 290 291 292
293
+
A4 u2M (t)
dt
And we seek to find optimal controls u∗V , u∗M , such that J(u∗V , u∗M ) = min{J(uV , uM )} U
294
(4)
subject to the differential equations (3), where tf is the final time. This performance specification involves minimizing the numbers of exposed and infected individuals, along with the cost of applying the controls (uV (t), uM (t)). In this paper, a quadratic objective functional is implemented for measuring the control cost, since the costs of the intervention are nonlinear. This assumption is based on the fact that there are no linear relationship between the effects of intervention and the cost of intervention of the infective populations, such quadratic costs have been frequently used [1, 2, 9]. The coefficients, Ai , i = 1, · · · , 4, are balancing cost factors. The controls (uV (t), uM (t)) is a bounded Lebesgue integrable functions [1, 2, 9].
PT
288
A1 EH + A2 IH +
A3 u2V (t)
CE
287
AC
285 286
tf
ED
J(uV , uM ) = 284
(3)
CR IP T
dSH dt dVH dt dEH dt dIH dt dRH dt dSM dt dEM dt dIM dt
AN US
282
M
281
where the control set, U = {(uV , uM ) : [0, tf ] → [0, 1], (uV (t), uM (t)), is Lebesgue measurable}.
16
ACCEPTED MANUSCRIPT
295
296 297 298 299
The necessary conditions that an optimal control quintuple must satisfy come from the Pontryagin’s Maximum Principle [48]. This principle converts (3) and (4) into a problem of minimizing pointwise a Hamiltonian H, with respect to the controls (uV (t) and uM (t)). First we formulate the Hamiltonian from the cost functional (4) and the governing dynamics (3) to obtain the optimality conditions.
CR IP T
300
Characterization of optimal controls
305 306 307
308
PT
304
∗ ∗ ∗ Theorem 2. Given an optimal control quintuple (u∗V , u∗M ) and solutions SH , VH∗ , EH , IH , ∗ ∗ ∗ ∗ ∗ ∗ RH , SM , EM , IM of the corresponding state system (3) that minimizes J(uV , uM ) over U. Then there exists adjoint variables λSH , λVH , λEH , λIH , λRH , λSM , λEM , λIM , satisfying
CE
303
where λSH , λVH , λEH , λIH , λRH , λSM , λEM , λIM are the associated adjoints for the states SH , VH , EH , IH , RH , SM , EM , IM . The system of adjoint equations is found by taking the appropriate partial derivatives of the Hamiltonian (5) with respect to the associated state and control variables.
AC
301 302
ED
M
AN US
H = A1 EH + A2 IH + A3 u2V + A4 u2M i h bM (αH EM + βH IM )SH − (uV + µH )SH + λSH ΛH + ωH VH − NH n o (1 − ε)bM (αH EM + βH IM )VH + λVH uV SH − − (ωH + µH )VH NH o n b (α E + β I )[S + (1 − ε)V ] M H M H M H H − (σH + µH )EH + λEH NH h i h i + λIH σH EH − (γH + µH + δH )IH + λRH γH IH − µH RH h i bM (αM EH + βM IH )SM + λSM ΛM − − µM SM − uM (t)SM NH h b (α E + β I )S i M M H M H M + λEM − (σM + µM )EM − uM (t)EM NH h i + λIM σM EM − µM IM − uM (t)IM
−
∂H dλi = dt ∂i
and with transversality conditions λi (tf ) = 0, where i = SH , VH , EH , IH , RH , SM , EM , IM .
309
The optimality conditions is given as ∂H = 0, ∂uj 310
(5)
j = 1, 2.
Furthermore, the control (u∗V , u∗M ) is given as 17
(6)
ACCEPTED MANUSCRIPT
u∗M
313 314
−
dλSH dt
dλRH − dt −
dλSM dt
dλI − M dt
=
∂H , ∂RH
λRH (tf ) = 0,
=
∂H , ∂SM
λSM (tf ) = 0,
∂H , ∂IM
λIM (tf ) = 0,
···
··· =
ED
= (λSH − λVH )uV + µH λSH + λ∗H +λ∗M
dλVH dt
PT
dλSH dt
∗ SM (λEM − λSM ), NH∗
CE
316
λSH (tf ) = 0,
evaluated at the optimal controls and corresponding state variables, results in the stated adjoint system (5) and (6) given by
+(λVH − λEH )(1 − ε)λ∗H
dλEH dt
∗ ∗ ) (NH∗ − SH ∗ VH (λ − λ ) + (1 − ε)λ (λEH − λVH ) SH EH H NH∗ NH∗
= ωH (λVH − λSH ) + µH λVH + (λEH − λSH )λ∗H
AC
315
∂H , ∂SH
=
AN US
312
Proof. The existence of an optimal control is guaranteed using the result by Fleming and Rishel [33]. Thus, the differential equations governing the adjoint variables are obtained by the differentiation of the Hamiltonian function, evaluated at the optimal controls. Thus, the adjoint system can be written as:
M
311
(7)
CR IP T
u∗V
" # ∗ SH (λ∗SH − λ∗VH ) , = min 1, max 0, 2A3 " # ∗ ∗ ∗ λ∗SM SM + λ∗EM EM + λ∗IM IM = min 1, max 0, . 2A4
= (λEH − λSH )
∗ ∗ SH ∗ SM + (λ − λ )λ E S M M M NH∗ NH∗
(NH∗ − VH∗ ) , NH∗
∗ ∗ ∗ ∗ λ∗H SH ∗ VH ∗ (NH − EH ) + (λ − λ )(1 − ε)λ + (λ − λ )b α S EH VH SM EM M M M H NH∗ NH∗ NH∗ 2
+σH (λEH − λIH ) + µH λEH − A1 ,
dλIH dt
= (λEH − λSH )
∗ ∗ ∗ λ∗H SH λ∗H VH∗ ∗ (NH − IH ) + (λ − λ )(1 − ε) + (λ − λ )b β S EH VH SM EM M M M NH∗ NH∗ NH∗ 2
+γH (λIH − λRH ) + (µH + δH )λIH − A2 , dλRH dt
= (λEH − λSH )
∗ λ∗H ∗ λ∗H ∗ ∗ SM ) )λ S + (λ − λ V (1 − ε) + (λ − λ + µH λRH , EH VH EM SM M NH∗ H NH∗ H NH∗
18
ACCEPTED MANUSCRIPT
dλSM dt dλEM dt
= (λSM − λEM )λ∗M + (µM + uM )λSM , = (λSH − λEH )bM
αH αH ∗ S + (1 − ε)(λVH − λEH )bM ∗ VH∗ + σM (λEM − λIM ) NH∗ H NH
+(µM + uM )λEM ,
319
Furthermore, differentiating the Hamiltonian function with respect to the control variables in the interior of the control set and then solving for controls (u∗V , u∗M ) result in the optimality conditions given as ∂H ∂uV ∂H ∂uM
320
∗ = 2A3 uV − SH (λ∗SH − λ∗VH ) = 0,
∗ ∗ ∗ = 2A4 uM − λ∗SM SM + λ∗EM EM + λ∗IM IM = 0.
Solving for u∗V and u∗M , we have,
u∗M
326 327 328 329 330
331 332
PT
Remark 5.1. Due to the a priori boundedness of the state and adjoint functions and the resulting Lipschitz structure of the ODE’s, the uniqueness of the optimal control for small time (tf ) was obtained. The uniqueness of the optimal control quintuple follows from the uniqueness of the optimality system, which consists of (3) and (5), (6) with characterization (7). The restriction on the length of the time interval is to guarantee the uniqueness of the optimality system, the smallness in the length of time is due to the opposite time orientations of (3), (5), and (6); the state problem has initial values, and the adjoint problem has final values. This restriction is very common in control problems (see [1, 2, 9]).
CE
324 325
Using the bounds on the controls, the characterization (7) can be derived.
AC
322 323
∗ (λ∗SH − λ∗VH ) SH , 2A3 ∗ ∗ ∗ + λ∗IM IM + λ∗EM EM λ∗SM SM . = 2A4
=
ED
u∗V
321
AN US
318
bM β H ∗ bM β H SH + (λVH − λEH ) ∗ VH∗ (1 − ε) + (µM + uM )λIM . ∗ NH NH
M
317
= (λSH − λEH )
CR IP T
dλIM dt
Next, we discuss the numerical solutions of the optimality system, the corresponding optimal control and the interpretations from various cases.
19
ACCEPTED MANUSCRIPT
340 341
342 343 344 345 346 347 348 349 350 351 352 353 354 355 356
357 358 359 360 361 362 363
CR IP T
339
Thus, to illustrate the optimal control strategies, the following values were taken for the initial conditions to mimic the Peshawar District of Pakistan. According to the Pakistan Bureau of Statistics the population of Peshawar District in the 2017 census is estimated as 1,970,042 [47]. We take the total initial population as NH (0) = 1970 per 1000 of the total population. We assume that initial population of vaccinated individuals VH (0) = 10; the initial population of exposed is given as EH (0) = 10. The initial infected is as given in the data is I(0) = 1. We assume there are no recovered individuals initially, so R(0) = 0. Hence, the initial susceptible population is given as S(0) = NH (0) − V (0) − EH (0) − IH (0) − RH (0) = 1959. We assume the initial population for the mosquitoes are given as SM (0) = 1338, EM (0) = 10, IM (0) = 152. For the weight factors we choose A1 = 1.0, A2 = 1.0, A3 = 0.0001, and A4 = 0.01. It should be pointed out that the weights in the simulations here are only of theoretical sense to illustrate the control strategies proposed in this paper. The reproduction number in the absence of intervention is given ˜ 0 ≈ 2.65, this was obtained from using the estimated and fitted parameter values in as R Table 2. This value indicate that the disease is epidemic in the district.
AN US
338
M
337
ED
336
The following algorithm was used to compute the optimal controls and state values using a Runge-Kutta method of the fourth order. First, an initial estimate for the control quintuple is made. Then the state variables are solved forward in time using the dynamics (3). The results obtained for the state variables are plugged into the adjoint equations (5). These adjoint equations with given final conditions (7) are then solved backward in time, employing the backward fourth order Runge-Kutta method. Both the state and adjoint values are then used to update the control, and the process is repeated until the current state, adjoint, and controls values converge sufficiently [41].
The results of the optimal control simulations of the dengue model (3) are depicted in Figure 6. With optimal control, most of the susceptible populace are either protected or vaccinated from the virus (see figures 6(a) and 6(b)) thereby leading to fewer individuals at the risk of infection (see Figure 6(c)) and eventually resulting in fewer individuals recovering from the infection as shown in Figure 6(d). This is not the case in the absence of the control strategy. For the mosquitoes, there are fewer of them with the control than in the absence of control as observed in Figures 6(e) and 6(f).
PT
335
Numerical Illustration
CE
334
6
AC
333
20
1500
1000
500
0
20
40
60
80
1500
without control with control
1000
500
100
0
20
Time (days)
60
80
ED
1000
PT
800 600 400 200 0
20
40
60
80
100
80
100
400
200
0
20
40
60
Time (days)
(d)
1200
without control with control
1000 800 600 400 200 0
100
20
40
60
80
100
Time (days)
Time (days)
AC
80
600
0
without control with control
CE
Susceptible mosquitoes
(c)
1200
without control with control
M
Time (days)
100
Recovered individuals
500
800
AN US
1000
Infectious mosquitoes
Infectious individuals
1500
40
60
(b)
without control with control
20
40
Time (days)
(a)
0
CR IP T
without control with control
Vaccinated individuals
Susceptible individuals
ACCEPTED MANUSCRIPT
(e)
(f)
Figure 6: Simulation results of model (3) with and without controls. (a) Susceptible individuals; (b) Vaccinated individuals; (c) Infectious individuals (exposed and Infected); and (d) Recovered individuals; (e) Susceptible mosquitoes; and (f) Infectious mosquitoes (exposed and Infected). Using A1 = 1.0, A2 = 1.0, A3 = 0.0001, A4 = 0.01.
364 365 366
The corresponding time-dependent controls (uV (t) and uM (t)) are depicted in Figure 7. The time-dependent controls uV (t) is at the upper bound of unity for about 10 days before gradually decreasing to the lower bound at the end of the simulation period; while the 21
ACCEPTED MANUSCRIPT
368 369 370
control uM (t) starts at the upper bound for over 60 days and slowly decrease to the lower bound at the end of the simulation period. These results suggest that to prevent an outbreak, vaccination control should start at a relatively high levels and while the vector control is maintain at high levels for a long time.
CR IP T
367
uV
1
uM
0.6
AN US
control
0.8
0.4
0
0
20
M
0.2
40
60
80
100
ED
Time (days)
373 374
375
CE
372
Varying cost weight on controls uV (t) and uM (t) In this section, we investigate the effect of varying the weights A3 , and A4 on the controls uV (t) and uM (t). These coefficients, Ai , i = 3, 4, are balancing cost factors on the cost of implementing the control strategies.
AC
371
PT
Figure 7: Simulation results of the control profile of model (3). Using A1 = 1.0, A2 = 1.0, A3 = 0.0001, A4 = 0.01.
Varying the weight A3 on control uV (t)
381
First we vary the weight A3 on the control uV (t) representing the use of vaccination. We use the weights A3 = 0.001, A3 = 0.01, A3 = 5.0 and A3 = 10.0 and keep constant A1 = A2 = 10.0 and A4 = 1.00. These costs correspond to the case with very cheap, cheap, expensive and very expensive vaccination strategies. The solution profile are depicted in Figure 8. The solution profile for the other variables are similar to what we observed in Figures 6 and are therefore not show here.
382
Figure 8 show that as the weight A3 increases from low to high (i.e., from cheap to very
376 377 378 379 380
22
ACCEPTED MANUSCRIPT
384 385 386 387
expensive cost), control uV (t) decreases from the higher levels to very low levels. The insecticide control uM (t) remain high, starting at the upper bound for about 10 days. The control uM slightly increases as weights A3 increases. This slight adjustment in control uM (t) help to keep the community protected and the infection low as observed in the sub-populations depicted in Figure 6.
uV
1
CR IP T
383
uV
1
uM
uM
0.8
control
control
0.8 0.6 0.4 0.2
20
40
60
80
100
0
uV
M
0.8
ED
0.6
20
PT
0
40
40
60
80
100
60
80
uV
1
uM
0.2
20
(b)
1
0.4
0
Time (days)
(a)
control
AN US
0
Time (days)
0
0.4 0.2
uM
0.8
control
0
0.6
0.6 0.4 0.2 0
100
Time (days)
0
20
40
60
80
100
Time (days) (d)
CE
(c)
Figure 8: Simulation results of the control profile of model (3) using A1 = A2 = 10.0, A4 = 1.00
388
AC
with cost weights (a) A3 = 0.001; (b) A3 = 0.01; (c) A3 = 5.0; and (d) A3 = 10.0.
Varying the weight A4 on control uM (t)
394
Next, we vary the weight A4 on control uM (t) which corresponds to insecticide use. We keep constant the weights A1 = A2 = 10.0, A3 = 0.001 and use the weights A4 = 0.01, A4 = 0.10, A4 = 5.0, and A4 = 10.0. These costs as with vaccination correspond to the case with very cheap, cheap, expensive and very expensive insecticide use. The control solution profile are depicted in Figure 9. The solution profile for the other variables are similar to what we observed in Figures 6 and are also not show here.
395
At low A4 values (i.e., A4 = 0.01) which corresponds to really cheap vaccine, the control
389 390 391 392 393
23
ACCEPTED MANUSCRIPT
397 398 399
400 401 402 403 404
uM (t) is usually at the upper bound and remain at the upper bound longer than control uV (t), the vaccination rate. However, as A4 gets bigger, the level of the control uM (t) falls below the control uV (t) as the level of control uV starts to rise. Indicating the existence of a reciprocal type relation between the weight A4 and control uV (t) [32]. This reciprocal relation between the weight A4 and control uV (t) is strongerthan the reciprocal relationship between A3 and control uM (t) (see Figure 8). Note, this relationship help to keep the community protected and the infection low as observed in Figure 6. It is possible to get a reciprocal relationship with a different strength depending on the parameter set used.
1
control
uM
0.4 0.2
0.6 0.4 0.2
0
20
40
60
80
100
0
ED
(a)
1 0.8
0.2
CE
0.4
AC
40
60
80
60
80
100
80
100
(b)
0.6 0.4 0.2
uM
20
40
Time (days)
0.8
uV
0
20
1
PT
0.6
0
M
Time (days)
control
AN US
0.8
control
control
uM
uV
0.6
0
uV
1
0.8
0
CR IP T
396
0
100
Time (days)
uV uM
0
20
40
60
Time (days)
(c)
(d)
Figure 9: Simulation results of the control profile of model (3) using A3 = 0.0001 with cost weights (a) A4 = 0.001; (b) A4 = 0.01; (c) A4 = 10.0; and (d) A4 = 13.0.
405 406 407
Next, we compare in Figure 10 the impart of these weights on the infection averted ratio (IAR), average cost-effectiveness ratio (ACER), infection averted, and total cost discussed in details in the next section below.
24
ACCEPTED MANUSCRIPT
408
Benefits gained from varying controls uV (t) and uM (t)
413
Infection averted ratio
414
The infection averted ratio (IAR) is stated as:
410 411
IAR =
CR IP T
412
Furthermore, we will investigate the associated benefits gained from implementing these controls while varying Ai , i = 3, 4, using cost-effectiveness analysis [2, 6]. To achieve this, we will consider three approaches, the infection averted ratio (IAR), and the total cost obtain from the objective functional in (4).
409
Number of infection averted . Number of recovered
(8)
418
Average cost-effectiveness ratio (ACER)
416
419 420
Next, we consider the average cost-effectiveness ratio (ACER) which deals with a single intervention, evaluating it against the no intervention baseline option. ACER is calculated as Total cost produced by the intervention . (9) ACER = Total number of infection averted
425 426
427 428 429 430 431
PT
424
The results depicted in Figures 10(c) and 10(d) are expected, lower weights are more cost-effective than higher weights, surprisingly, insecticide led to lower cost compare to vaccination at higher weights. This is due to the strong reciprocal effect we observed in Figures 9(c) and 9(d) when A4 was been varied and the control uV adjusted to compensate for the low levels of uM use.
CE
423
The two weights averted approximately the same number of infection regardless of the cost (see Figure 10(a), and 10(b)). This is due to the reciprocal effect we observed in Figures 8 and 9, as the two controls uV and uM adjust to compensate for the low levels in one when the other is high. Thus, the two controls exhibit this sea-saw type behavior, thus helping to keep the level of infection the same regardless of the cost.
AC
422
ED
M
421
AN US
417
The number of infection averted is given as the difference between the total infectious individuals without control and the total infectious individuals with control. The strategy with the highest ratio is the most cost-effective.
415
25
ACCEPTED MANUSCRIPT
5
A3
Infection Averted
A4
10
5
0
0.001
0.01
5.0
A4
4 3 2 1 0
10.0
A3
0.001
Weights
10.0
0.6
0.2
0.01
5.0
Weights
M
0.4
A3 A4
10
Total Cost
0.8
×105
AN US
A4
12
10.0
ED
Average Cost-Effectiveness Ratio (ACER)
A3
(c)
5.0
(b)
1
0.001
0.01
Weights
(a)
0
CR IP T
×105
Infection Averted Ratio (IAR)
15
8 6 4 2 0
0.001
0.01
5.0
10.0
Weights
(d)
PT
Figure 10: Simulation results of the control profile of model (3) using cost weights A3 = A4 =
433 434 435 436 437 438 439 440
441 442 443
Control within interval of uncertainty for bM and µM In carrying out the sensitivity analysis test we set the values of the parameters to be within ± 20% range of their baseline values and found the parameters bM , and µM among ˜ 0 and R0 . However, the time others to have a high impact on the reproduction numbers R dependent controls uV (t) and uM (t) incorporated into system (3) did not directly act on these parameters. Hence, we vary first, the parameter bM using a range of ± 50% from the baseline value and determine its impact on the system and the control variables. We have used the ± 50% range for emphasis purpose, the result is expected to be the same in the ± 20% range.
AC
432
CE
0.0001, A3 = A4 = 0.001, A3 = A4 = 10.0, and A3 = A4 = 13 (a) Infection averted; (b) Infection averted ratio (IAR); (c) Average cost-effective ratio (ACER); (d) Total cost.
The results of the optimal control simulations of the dengue model (3) in the uncertainty interval for the parameter bM are depicted in Figures 11(a) - 11(b). The parameter bM is the mosquitoes biting rate, since Aedes egypyti bite during the day, knowledge of the 26
ACCEPTED MANUSCRIPT
impact of the control in the uncertainty range of this parameter is vital. 1400
Lower bound Baseline values Upper bound
1500
1000
500
0
0
20
40
60
80
Lower bound Baseline values Upper bound
1200 1000 800 600 400 200 0
100
0
20
Time (days)
CR IP T
Infectious individuals
2000
Infectious mosquitoes
444
40
60
80
100
Time (days) (b)
AN US
(a)
Figure 11: Simulation of the dengue model (3) as a function of time without control (solid lines) and with optimal control (dash lines) for: (a). Total number of infectious humans; (b). Total number of infectious mosquitoes. The parameter bM is set within ± 50% interval of its baseline values given in Table 2 and A1 = 1.0, A2 = 1.0, A3 = 0.0001, A4 = 0.01.
448
1
PT
0.6
0.8
CE
control uV
0.8
0.4
AC
0.2 0
1
Lower bound Baseline values Upper bound
control uM
447
M
446
We observed that the upper bound led to the most number of infectious in humans (EH and IH ), and mosquitoes. However, due to the actions of the control variables depicted in Figures 12(a)-12(b), the number of the infectious humans was drastically reduced while the infectious mosquitoes remain the same within the period of simulation.
ED
445
0
20
40
0.6 0.4 Lower bound Baseline values Upper bound
0.2
60
80
0
100
Time (days)
0
20
40
60
80
100
Time (days)
(a)
(b)
Figure 12: The optimal controls of the dengue model (3) in the ± 50% uncertainty interval of parameter bM for: (a). control uV ; (b). control uM .
449 450
Hence, regardless of the parameter values chosen within the uncertainty interval of bM the control profile for the mosquitoes remained the same. Although, the duration of the 27
ACCEPTED MANUSCRIPT
455
456 457 458 459 460 461 462
Next, we investigate the optimal control of system (3) for different values of the mosquito life span; we set µM = 1/14, 1/21, 1/42 [43]. The results are depicted in Figures 13(a)13(b). We observed the same profile for infectious humans regardless of the mosquito life span either with or without control. The control led to fewer infectious humans. For the mosquitoes, there were more infected mosquitoes as expected when the mosquitoes live longer (µM = 1/42) and fewer infected ones when they live for a shorter period of time (µM = 1/14).
CR IP T
453 454
controls at the maximum level are longer when parameter bM is at the upper bound and shorter at the lower bound. The duration of the control uV is marginally different at the upper bound. Note, that we have used very low weights, the profile is the same even for higher weights and tighter intervals as the ± 20% uncertainty interval used for the sensitivity analysis.
µ = 1/42 M
µ = 1/21 M
µM = 1/14
1500
1000
0
M
500
1400
0
20
40
60
80
100
ED
Infectious individuals
2000
AN US
452
µ = 1/42 M
1200
µ = 1/21 M
µ = 1/14
1000
M
800 600 400 200 0
0
20
40
60
80
100
Time (days)
Time (days) (a)
Infectious mosquitoes
451
(b)
PT
Figure 13: The optimal controls of the dengue model (3) when µM = 1/14, 1/21, 1/42 for : (a).
464 465 466
For the optimal controls uV (t) and uM (t), there are marginal differences in the control profiles in Figures 14(a)-14(b). The implication of these results are that, regardless of the mosquito life span the control efforts are the same and they will adequately reduce the infectious populations both in humans and mosquitoes.
AC
463
CE
Total number of infectious humans; (b). Total number of infectious mosquitoes. The parameter values used are given in Table 2 and A1 = 1.0, A2 = 1.0, A3 = 0.0001, A4 = 0.01.
28
ACCEPTED MANUSCRIPT
1
1 µM = 1/42 µ = 1/21
0.8
0.6 0.4 0.2 0
µ = 1/42 M µ = 1/21
0.6
M
µ = 1/14 M
0.4 0.2
0
20
40
60
80
0
100
Time (days)
0
20
CR IP T
M
µM = 1/14
control uM
control uV
0.8
40
60
80
100
Time (days)
(a)
(b)
471 472 473 474
475 476 477 478 479 480 481
482 483 484 485 486 487
488 489
M
ED
470
Dengue is a viral infection transmitted via infected female Aedes mosquito bites resulting in either a mild dengue fever or haemorrhagic fever [65]. There are four distinct dengue virus serotypes (namely DEN-1, DEN-2, DEN-3 and DEN-4) [65]. Since 2010, three provinces in Pakistan namely Khyber Pakhtunkhwa, Punjab and Sindh has been experiencing dengue fever epidemic [66]. The first confirmed dengue fever outbreak in the country was in reported in 1994; however, the first annual epidemic trend and sudden rise in cases occurred in Karachi in November 2005 [66].
PT
469
Discussion and Conclusion
The government of Pakistan solicited technical assistance from the World health Organization (WHO) to combat the epidemic by emphasizing on vector control, case management and community awareness [66]. Unfortunately, these strategies were not effective, as at 2010, there was another sudden rise in cases with more reported cases than the preceding five years combined [66]. The country has since faced a number of repeated outbreaks in different provinces [69]. In 2017, the country experienced another outbreak starting from the province of Khyber Pakhtunkhwa spreading to Sindh province, and AJK province [69].
CE
468
7
AC
467
AN US
Figure 14: The optimal controls of the dengue model (3) when µM = 1/14, 1/21, 1/42 for: (a). Control uV ; (b). Control uM .
In this paper, we applied optimal control theory and investigate the impact of using an imperfect vaccine and mosquito reduction strategy in the bid to control dengue in Pakistan. For our investigation we use the dengue model (1) proposed in Section 2. The analysis shows that the disease-free equilibrium of the model is locally asymptotically stable whenever the associated reproduction number (R0 ), is less than unity and unstable otherwise. ˜0 To identify the parameters with the strongest impact on the reproduction numbers (R and R0 ), the model outcome, we used Latin Hypercube Sampling (LHS) and partial rank 29
ACCEPTED MANUSCRIPT
498 499 500 501 502 503 504 505
506 507 508 509 510 511 512 513 514
515 516 517 518 519 520 521 522 523 524 525 526 527 528
529 530
CR IP T
496 497
In particular, the results of this sensitivity analysis suggest that a strategy that reduces the transmission probability per contact of infectious mosquitoes with susceptible humans (αH and βH ), mosquito recruitment rate (ΛM ), mosquito biting rate (bM ), transmission probability per contact of susceptible mosquitoes with infectious humans (αM and βM ), disease progression rate in mosquitoes (σM ), and vaccine waning rate (ωH ). Furthermore, a strategy that increases the human recruitment rate (ΛH ), disease progression rate in humans (σH ), human recovery rate (γH ), mosquito natural death rate (µM ), and vaccine efficacy (ε) will be effective in curtailing the spread of dengue in the community.
AN US
495
Hence, we introduced into the model two time-dependent control variables representing vaccination and the use of insecticide. We found that the total number infectious (exposed and infected) individuals can be reduced in the community by the application of time dependent controls. Our results are in line with the results obtained in [34,35,39] which respectively showed that vaccine when available can be a promising tool to fight the disease in a community. Our study used a simple model and strategy to reach the same conclusion unlike Rodrigues et al. [35] which investigated using optimal control strategies the different ways of distributing vaccines: pediatric and random mass vaccines, with distinct levels of efficacy and durability.
M
494
ED
493
PT
492
correlation coefficients (PRCC) methods. The results of the sensitivity analysis of the model show that the parameters that have the most influence the reproduction number, R0 , are the human recruitment rate (ΛH ), disease progression rate in humans (σH ), human recovery rate (γH ), mosquito natural death rate (µM ), and vaccine efficacy (ε), the transmission probability per contact of infectious mosquitoes with susceptible humans (αH and βH ), mosquito recruitment rate (ΛM ), mosquito biting rate (bM ), transmission probability per contact of susceptible mosquitoes with infectious humans (αM and βM ), disease progression rate in mosquitoes (σM ), and vaccine waning rate (ωH ).
The application of optimal control theory requires an objective/cost functional. The optimal control objective or cost functional we utilized, minimized the numbers of exposed and infected individuals along with the cost of applying the controls. We choose a quadratic objective functionalin this study, since the costs of the intervention are nonlinear. This assumption is based on the fact that there are no linear relationship between the effects of intervention and the cost of intervention of the infective populations. Rodrigues et al. [34] used a cost functional which depends on the costs of medical treatment of the infected people, and the costs related to educational and sanitation campaigns. They then went on to show efficient ways to decrease the number of infected mosquitoes and individuals in less time and with lower costs using current computational tools. Although, that was not our approach in this study, we nevertheless investigated the associated benefits of the control strategies using cost-effectiveness analysis. This was achieved, using three approaches, the infection averted ratio (IAR), the average cost-effectiveness ratio (ACER) and the objective functional.
CE
491
AC
490
We observed while investigating the effect of the costs on the cost weights a reciprocal relation exist between the cost of insecticde use and vaccination; as the cost of insecticide 30
ACCEPTED MANUSCRIPT
538 539 540 541 542 543
544 545 546 547
548 549 550 551 552
553
554 555 556 557
558 559 560 561 562 563 564 565
566 567 568
569
CR IP T
537
AN US
536
Furthermore, our study showed that lower weights on the control are more cost-effective than higher weights. That is controls with lower weights are cheaper strategies to use compare to controls with higher weights. Insecticide led to lower cost compare to vaccination at higher weights; therefore, insecticide is more cost-effect that vaccination at lower weight. Shim [39] in their study of the dynamics of dengue and cost-effectiveness of vaccine in the Philippines using an age-structured model showed that the higher the cost of vaccination, the less cost-effective the dengue vaccination program. Furthermore, their results demonstrated that at a relatively low vaccine efficacy an age-targeted vaccination may be cost-effective provided the cost is sufficiently low. Although, we did not investigate the effect of vaccine efficacy on cost-effectiveness and overall performance; it is easy to see that due to the reciprocal relationship between insecticide use and vaccination the reduce performance in vaccination due to reduced efficacy will be off-set by increase in insecticide use. To conclude, in this study, we have presented a deterministic model of ordinary differential equations for dengue infection transmission dynamics. And we have studied using optimal control theory the use of insecticide, and vaccination as effective control measures against the epidemics using data obtained from hospital from the Peshawar District in Pakistan. The following results were observed from our analysis and numerical simulations:
M
535
ED
534
(i) The model has a DFE that is locally asymptotically stable if R0 < 1; (ii) Hence, we estimated in this study the basic reproduction number (R0 ) without any interventions for the 2017 dengue outbreak in Peshawar district of Pakistan as R0 ≈ 2.64, the distribution of the reproduction number lies in the range R0 ∈ [1.21, 5.24] (with a mean R0 ≈ 2.64);
PT
533
increases the use of vaccination increases. Only a slight increase in insecticide use is observed when vaccination level decreased due to increase in cost. The two controls averts about the same number of infections in the district regardless of the weights on the costs this is due to the reciprocal relation existing between the cost of insecticde use and vaccination.
CE
532
(iii) The sensitivity analysis of the model shows that the dominant parameters influencing ˜ 0 and R0 ) are the human recruitment rate (ΛH ), disthe reproduction numbers s (R ease progression rate in humans (σH ), human recovery rate (γH ), mosquito natural death rate (µM ), and vaccine efficacy (ε), the transmission probability per contact of infectious mosquitoes with susceptible humans (αH and βH ), mosquito recruitment rate (ΛM ), mosquito biting rate (bM ), transmission probability per contact of susceptible mosquitoes with infectious humans (αM and βM ), disease progression rate in mosquitoes (σM ), and vaccine waning rate (ωH ).
AC
531
(iv) The application of two time-dependent controls uV (t) and uM (t) obtained from the sensitivity analysis can reduce the total number of infected humans and mosquitoes in the district; (v) Investigating the effect of the costs on the weights A3 , and A4 we observed that 31
ACCEPTED MANUSCRIPT
571 572 573 574 575 576 577 578
579
(a) A reciprocal relation exist between the cost of insecticde use and vaccination; as the cost of insecticide increases the use of vaccination increases. Only a slight increase in insecticide use is observed when vaccination level decreased due to increase in cost. (b) The two controls averts about the same number of infections in the district regardless of the weights on the costs this is due to the reciprocal relation exist between the cost of insecticde use and vaccination.
CR IP T
570
(c) Lower weights are more cost-effective than higher weights with insecticide leading to cheaper cost than vaccination.
Acknowledgments
582
References
588
589 590 591
592 593 594
595 596
597 598
599 600
601 602 603
M
ED
587
[3] F.B. Agusto, S. Bewick, and W.F. Fagan. Mathematical model of zika virus with vertical transmission. Infectious Disease Modelling, 2(2):244–267, 2017. [4] F.B. Agusto, S.Y. Del Valle, K.W. Blayneh, C.N. Ngonghala, M.J. Goncalves, N. Li, R. Zhao, and H. Gong. The impact of bed-net use on malaria prevalence. Journal of theoretical biology, 320:58–65, 2013.
PT
586
[2] F.B. Agusto. Optimal isolation control strategies and cost-effectiveness analysis of a two-strain avian influenza model. BioSystems, 113(3):155–164, 2013.
[5] F.B. Agusto, S. Easley, K. Freeman, and M. Thomas. Mathematical model of three age-structured transmission dynamics of chikungunya virus. Computational and mathematical methods in medicine, 2016, 2016.
CE
585
[1] F. Agusto and S. Lenhart. Optimal control of the spread of malaria superinfectivity. Journal of Biological Systems, 21(04):1340002, 2013.
AC
583 584
AN US
581
One of the authors, FBA, was supported by the Strategic Environmental Research and Development Program under grant RC-2639.
580
[6] F.B. Agusto and I.M. ELmojtaba. Optimal control and cost-effective analysis of malaria/visceral leishmaniasis co-infection. PloS one, 12(2):e0171102, 2017.
[7] F.B. Agusto and A.B. Gumel. Theoretical assessment of avian influenza vaccine. DCDS Series B, 13(1):1–25, 2010. [8] F.B. Agusto and A.B. Gumel. Qualitative dynamics of lowly-and highly-pathogenic avian influenza strains. Mathematical biosciences, 243(2):147–162, 2013.
[9] F.B. Agusto, N. Marcus, and K.O. Okosun. Application of optimal control to the epidemiology of malaria. Electronic Journal of Differential Equations, 2012(81):1–22, 2012. 32
ACCEPTED MANUSCRIPT
605
[10] T.L Bancroft. On the aetiology of dengue fever. Australian Medical Gazette, 25:17–18, 1906.
606
[11] G. B´eraud. Mathematical models and vaccination strategies. Vaccine, 2017.
611 612 613
614 615 616 617
618 619 620
621 622 623
624 625 626 627
CR IP T
610
[13] J.E. Blaney, J.M. Matro, B.R. Murphy, and S.S. Whitehead. Recombinant, liveattenuated tetravalent dengue virus vaccine formulations induce a balanced, broad, and protective neutralizing antibody response against each of the four serotypes in rhesus monkeys. Journal of virology, 79(9):5516–5528, 2005. [14] J.E. Blaney, N.S. Sathe, C.T. Hanson, C.Y. Firestone, B.R. Murphy, and S.S. Whitehead. Vaccine candidates for dengue virus type 1 (den1) generated by replacement of the structural genes of rden4 and rden4δ30 with those of den1. Virology journal, 4(1):23, 2007.
AN US
609
[12] S. Bhatt, P.W. Gething, O.J. Brady, J.P. Messina, A.W. Farlow, C.L. Moyes, J.M. Drake, J.S. Brownstein, A.G. Hoen, O. Sankoh, et al. The global distribution and burden of dengue. Nature, 496(7446):504, 2013.
[15] K.W. Blayneh, A.B. Gumel, S. Lenhart, and T. Clayton. Backward bifurcation and optimal control in transmission dynamics of west nile virus. Bulletin of mathematical biology, 72(4):1006–1028, 2010. [16] S.M. Blower and H.I. Dowlatabadi. Sensitivity and uncertainty analysis of complex models of disease transmission: an hiv model, as an example. International Statistical Review/Revue Internationale de Statistique, pages 229–243, 1994.
M
608
ED
607
[17] O.J. Brady, P.W. Gething, S. Bhatt, J.P. Messina, J.S. Brownstein, A.G. Hoen, C.L. Moyes, A.W. Farlow, T.W. Scott, and S.I. Hay. Refining the global spatial limits of dengue virus transmission by evidence-based consensus. PLoS neglected tropical diseases, 6(8):e1760, 2012.
PT
604
630
[19] J Carr. Applications of Centre Manifold Theory. Springer-Verlag, New York, 1981.
632 633
634 635
636 637
638 639 640
[20] C. Castillo-Chavez, S. Blower, P. van den Driessche, D. Kirschner, and A.A. Yakubu. Mathematical approaches for emerging and reemerging infectious diseases: an introduction, volume 1. Springer Science & Business Media, 2002.
AC
631
CE
629
[18] F. Brauer. Backward bifurcations in simple vaccination models. Journal of Mathematical Analysis and Applications, 298(2):418–431, 2004.
628
[21] C. Castillo-Chavez and B. Song. Dynamical models of tuberculosis and their applications. Mathematical biosciences and engineering, 1(2):361–404, 2004.
[22] Centres for Disease Control and Prevention. Dengue. https: // www. cdc. gov/ dengue/ index. html , 2016. [23] G. Chowell, P. Diaz-Duenas, J.C. Miller, A. Alcazar-Velazco, J.M. Hyman, P.W. Fenimore, and C. Castillo-Chavez. Estimation of the reproduction number of dengue fever from spatial epidemic data. Mathematical biosciences, 208(2):571–589, 2007. 33
ACCEPTED MANUSCRIPT
649
650 651 652
653 654
655 656
657 658
659 660
661 662
663 664
665 666
667 668 669
670 671 672
673 674 675 676
CR IP T
648
[27] J. Dushoff, W. Huang, and C. Castillo-Chavez. Backwards bifurcations and catastrophe in simple models of fatal diseases. Journal of mathematical biology, 36(3):227–248, 1998. [28] E.H. Elbasha and A.B. Gumel. Theoretical assessment of public health impact of imperfect prophylactic hiv-1 vaccines with therapeutic benefits. Bulletin of mathematical biology, 68(3):577, 2006.
AN US
647
[29] C. Esteva, L.and Vargas. Influence of vertical and mechanical transmission on the dynamics of dengue disease. Mathematical biosciences, 167(1):51–64, 2000. [30] L. Esteva and C. Vargas. Analysis of a dengue disease transmission model. Mathematical biosciences, 150(2):131–151, 1998. [31] L. Esteva and C. Vargas. A model for dengue disease with variable human population. Journal of mathematical biology, 38(3):220–240, 1999.
M
646
[26] I. Dorigatti, C. McCormack, G. Nedjati-Gilani, and N.M. Ferguson. Using wolbachia for dengue control: Insights from modelling. Trends in parasitology, 2017.
[32] F.B. Agusto. Optimal control and cost-effective analysis of the 2017 meningitis outbreak in nigeria. Submitted manuscript, 2018.
ED
645
[33] W.H. Fleming and R.W. Rishel. Deterministic and stochastic optimal control, volume 1. Springer Science & Business Media, 2012.
PT
644
[25] M. Derouich and A. Boutayeb. Dengue fever: Mathematical modelling and computer simulation. Applied Mathematics and Computation, 177(2):528–544, 2006.
[34] S.M. Garba, A.B. Gumel, and M.R.A. Bakar. Backward bifurcations in dengue transmission dynamics. Mathematical biosciences, 215(1):11–25, 2008.
CE
643
[24] Dengue Virus Ne. Treatment of dengue. http: // www. denguevirusnet. com/ treatment. html , 2018.
[35] D.J. Gubler. Dengue and dengue hemorrhagic fever. Clinical microbiology reviews, 11(3):480–496, 1998. [36] D.J. Gubler and G. Kuno. Dengue and dengue hemorrhagic fever: its history and resurgence as a global public health problem. Dengue and dengue hemorrhagic fever. CAB international, London, United Kingdom, pages 1–22, 1997.
AC
641 642
[37] S.B. Halstead, S. Nimmannitya, and S.N. Cohen. Observations related to pathogenesis of dengue hemorrhagic fever. iv. relation of disease severity to antibody response and virus recovered. The Yale journal of biology and medicine, 42(5):311, 1970. [38] I. Kawaguchi, A. Sasaki, and M. Boots. Why are dengue virus serotypes so distantly related? enhancement and limiting serotype similarity between dengue virus strains. Proceedings of the Royal Society of London B: Biological Sciences, 270(1530):2241– 2247, 2003. 34
ACCEPTED MANUSCRIPT
683
684 685 686
687 688 689
690 691 692
693 694 695
696 697 698
[41] S. Lenhart and J.T. Workman. Optimal control applied to biological models. CRC Press, 2007.
CR IP T
682
[42] S.B. Maier, X. Huang, E. Massad, M. Amaku, M.N. Burattini, and D. Greenhalgh. Analysis of the optimal vaccination age for dengue in brazil with a tetravalent dengue vaccine. Mathematical biosciences, 294:15–32, 2017. [43] C.A. Manore, K.S. Hickmann, S. Xu, H.J. Wearing, and J.M. Hyman. Comparing dengue and chikungunya emergence and endemic transmission in a. aegypti and a. albopictus. Journal of theoretical biology, 356:174–191, 2014.
AN US
681
[44] S. Marino, I.B. Hogue, C.J. Ray, and D. E. Kirschner. A methodology for performing global uncertainty and sensitivity analysis in systems biology. Journal of theoretical biology, 254(1):178–196, 2008. [45] M.D. McKay, R.J. Beckman, and W.J. Conover. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics, 42(1):55–61, 2000.
M
680
[40] M.R. Khanani, A. Arif, and R. Shaikh. Dengue in pakistan: Journey from a disease free to a hyper endemic nation. Journal of the Dow University of Health Sciences Karachi, 5(3):81–84, 2011.
[46] R.G. McLeod, J.F. Brewster, A.B. Gumel, and A. Slonowsky. Sensitivity and uncertainty analyses for a sars model with time-varying inputs and outputs. Mathematical Biosciences and Engineering, 3(3):527, 2006.
ED
679
[39] J. Khan, A. Ghaffar, and S.A. Khan. The changing epidemiological pattern of dengue in swat, khyber pakhtunkhwa. PloS one, 13(4):e0195706, 2018.
PT
677 678
703
[48] L.S. Pontryagin. Mathematical theory of optimal processes. CRC Press, 1987.
701
704 705 706
707 708
709 710 711
AC
700
CE
702
[47] Pakistan Bureau of Statistics. Pakistan’s 6th census: Population of major cities census. http: // www. pbscensus. gov. pk/ sites/ default/ files/ population_ of_ major_ cities_ census_ 2017% 20_ 0. pdf , 2017.
699
[49] H.S. Rodrigues, M.T.T. Monteiro, and D.F.M. Torres. Dynamics of dengue epidemics when using optimal control. Mathematical and Computer Modelling, 52(9-10):1667– 1673, 2010.
[50] H.S. Rodrigues, M.T.T. Monteiro, and D.F.M. Torres. Vaccination models and optimal control strategies to dengue. Mathematical biosciences, 247:1–12, 2014. [51] M.A. Sanchez and S.M. Blower. Uncertainty and sensitivity analysis of the basic reproductive rate: tuberculosis as an example. American Journal of Epidemiology, 145(12):1127–1137, 1997.
35
ACCEPTED MANUSCRIPT
714
715 716 717
718 719 720
[52] T. Sardar, S. Rana, and J. Chattopadhyay. A mathematical model of dengue transmission with memory. Communications in Nonlinear Science and Numerical Simulation, 22(1-3):511–525, 2015. [53] O. Sharomi, C. Podder, A. Gumel, and B. Song. Mathematical analysis of the transmission dynamics of hiv/tb coinfection in the presence of treatment. Mathematical Biosciences and Engineering, 5(1):145, 2008.
CR IP T
712 713
[54] O. Sharomi, C.N. Podder, A.B. Gumel, E.H. Elbasha, and J. Watmough. Role of incidence function in vaccine-induced backward bifurcation in some hiv models. Mathematical Biosciences, 210(2):436–463, 2007. [55] C. Shekhar. Deadly dengue: new vaccines promise to tackle this escalating global menace. Chemistry & biology, 14(8):871–872, 2007.
723
[56] S.M. Shepherd. Dengue. Medscape Logo, 2017.
731
732 733
734 735 736
737 738 739
740 741
742 743 744 745
746 747
M
730
[59] M.I. Teboh-Ewungkem, C.N. Podder, and A.B. Gumel. Mathematical study of the role of gametocytes and an imperfect vaccine on malaria transmission dynamics. Bulletin of mathematical biology, 72(1):63–93, 2010.
ED
729
[60] J.J. Tewa, J.L. Dimi, and S. Bowong. Lyapunov functions for a dengue disease transmission model. Chaos, Solitons & Fractals, 39(2):936–941, 2009.
PT
727 728
[58] C.J. Struchiner, P.M. Luz, C.T. Codeco, F.C. Coelho, and E. Massad. Current research issues in mosquito-borne diseases modeling. Contemporary Mathematics, 410:349–366, 2006.
[61] P. Van den Driessche and J. Watmough. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical biosciences, 180(1):29–48, 2002.
CE
726
[57] E. Shim. Dengue dynamics and vaccine cost-effectiveness analysis in the philippines. The American journal of tropical medicine and hygiene, 95(5):1137–1147, 2016.
[62] WHO Health Organization (WHO). Who country cooperation strategic. http: // apps. who. int/ iris/ bitstream/ 10665/ 136607/ 1/ ccsbrief_ pak_ en. pdf , 2016.
AC
724 725
AN US
722
721
[63] World Health Organization (WHO). Dengue and severe dengue. http: // www. who. int/ mediacentre/ factsheets/ fs117/ en/ , 2017.
[64] World Health Organization (WHO). Dengue vaccine research. Immunization, Vaccines and Biologicals http: // www. who. int/ immunization/ research/ development/ dengue_ vaccines/ en/ , 2017. [65] World Health Organization (WHO). Epidemiology. http: // www. who. int/ denguecontrol/ epidemiology/ en/ , 2018. 36
ACCEPTED MANUSCRIPT
757
758 759 760 761
762 763
CR IP T
756
[68] World Health Organization (WHO), Eastern Mediterranean Regional Office (EMRO). Dengue fever in pakistan. Weekly Epidemiological Monitor https: // reliefweb. int/ sites/ reliefweb. int/ files/ resources/ Epi_ Monitor_ 2017_ 10_ 36. pdf , 10(36), 03 September 2017. Accessed March 16, 2018. [69] World Health Organization (WHO), Eastern Mediterranean Regional Office (EMRO). Deaths from dengue in pakistan. Weekly Epidemiological Monitor. https: // reliefweb. int/ sites/ reliefweb. int/ files/ resources/ Epi_ Monitor_ 2017_ 10_ 46. pdf , 10(46), 12 November 2017. Accessed March 16, 2018.
AN US
754 755
[70] H.M. Yang and C.P. Ferreira. Assessing the effects of vector control on dengue transmission. Applied Mathematics and Computation, 198(1):401–413, 2008.
M
753
[67] World Health Organization (WHO). Dengue and severe dengue. http: // www. who. int/ news-room/ fact-sheets/ detail/ dengue-and-severe-dengue , 2018. Accessed September 3, 2018.
ED
752
PT
751
CE
750
[66] World Health Organization (WHO). Pakistan. http: // www. emro. who. int/ pak/ programmes/ dengue-fever. html , 2018. Accessed July 24, 2018.
AC
748 749
37
ACCEPTED MANUSCRIPT
765 766 767
A
Appendix A: Proof of Lemma 1
Lemma 1: Let the initial data F (0) ≥ 0 , where F (t) = (SH (t), VH (t), EH (t), IH (t), RH (t), SM (t), EM (t), IM (t)). Then the solutions F (t) of the model (1) are non-negative for all t > 0. Furthermore lim sup NH (t) ≤ t→∞
ΛH ΛM , and lim sup NM (t) ≤ . µH µM t→∞
CR IP T
764
where
NH (t) = SH (t) + VH (t) + EH (t) + IH (t) + RH (t), and NM (t) = SM (t) + EM (t) + IM (t). 768
772
773
AN US
M
771
dSH bM (αH EM + βH IM )SH = ΛH + ωH VH − − (νH + µH )SH . (A-1) dt NH bM (αH EM + βH IM ) Let λH (t) = , then (A-1) becomes NH dSH = ΛH + ωH VH − λH (t)SH − (νH + µH )SH dt which can be re-written as ( Z t1 ) Z t1 d SH (t) exp λH (ζ)dζ + (νH + µH )t = (ΛH + ωH VH ) exp λH (ζ)dζ(νH + µH )t . dt 0 0
ED
770
Proof. Let t1 = sup{t > 0 : F (t) > 0 ∈ [0, t]}. Thus, t1 > 0. It follows from the first equation of the system (1), that
Hence, t1
λH (ζ)dζ + (νH + µH )t1
CE
SH (t1 ) exp
Z
PT
769
0
so that,
"
SH (t1 ) = SH (0) exp − "
+ exp − ×
exp
> 0.
− SH (0) = ×
AC
774
Z
" Z
Z
Z
t1
(ΛH + ωH VH (p)) Z p exp λH (ζ)dζ + (νH + µH )p dp 0
0
t1
λH (ζ)dζ + (νH + µH )t1
0
t1
λH (ζ)dζ + (νH + µH )t1
0
p
λH (ζ)dζ + (νH + µH )p
0
38
#
# Z
0
dp
#
t1
(ΛH + ωH VH (p))
ACCEPTED MANUSCRIPT
776 777 778
779
Similarly, it can be shown that F > 0 for all t > 0. For the second part of the proof, note that 0 < SH (0) ≤ NH (t), 0 ≤ VH (0) ≤ NH (t), 0 ≤ EH (0) ≤ NH (t), 0 ≤ (IH 0) ≤ NH (t), 0 ≤ RH (0) ≤ NH (t), 0 < SM (0) ≤ NM (t), 0 ≤ EM (0) ≤ NM (t), 0 ≤ IM (0) ≤ NM (t). Adding the human and mosquito components of model (1) gives
CR IP T
775
dNH (t) = ΛH − µH NH − δH IH dt ≤ ΛH − µH NH and
AN US
780
(A-2)
dNM (t) = ΛM − µM NM dt Hence, lim sup NH (t) ≤ t→∞
784
Lemma 2: The region Ω = ΩH × ΩM ⊂ R5+ × R3+ is positively-invariant for the model (1) with non-negative initial conditions in R8+ . Proof. It follows from summing up the human and mosquito equations of model (1) that
786
787
788
789
AC
CE
785
Appendix B: Proof of Lemma 2
ED
783
B
PT
782
M
as required.
781
ΛH ΛM , and lim sup NM (t) ≤ µH µM t→∞
dNH (t) = ΛH − µH NH − δH IH dt ≤ ΛH − µH NH
and
dNM (t) = ΛM − µM NM dt
dNH (t) ΛH dNM (t) ΛM ≤ 0, if NH (0) ≥ and ≤ 0, if NM (0) ≥ . dt µH dt µM ΛH ΛM Thus, NH (t) ≤ NH (0)e−µH t + (1 − e−µH t ) and NM (t) ≤ NM (0)e−µM t + (1 − e−µM t ). µH µM ΛH ΛM In particular, NH (t) ≤ and NM (t) ≤ . µH µM
Hence,
39
ACCEPTED MANUSCRIPT
794 795 796 797 798 799
CR IP T
793
C
Appendix C: Proof of Theorem 1
The phenomenon of backward bifurcation is shown using the concept of centre manifold theory on system (1). First, we make the following changes in variables by setting SH = x1 , VH = x2 , EH = x3 , IH = x4 , RH = x5 , SM = x6 , EM = x7 and − IM = x8 Furthermore, employing the vector notation → x = (x1 , x2 , x3 , x4 , x5 , x6 , x7 , x8 ), − the dengue model (1) with vaccination can be reformulated in the form dx/dt = F → x , with F = (f1 , f2 , f3 , f4 , f5 , f6 , f7 , f8 )T as given below = ΛH + ωH x2 − λH x1 − (νH + µH )x1 ,
(C-1)
= νH x1 − (1 − ε)λH x2 − (ωH + µH )x2 , = λH [x1 + (1 − ε)x2 ] − (σH + µH )x3 , = σH x3 − (γH + µH + δH )x4 , = γH x4 − µH x5 , = ΛM − λM x6 − µM x6 ,
where,
= λM x6 − (σM + µM )x7 , = σ M x7 − µ M x8 ,
AC
800
CE
PT
dx1 dt dx2 dt dx3 dt dx4 dt dx5 dt dx6 dt dx7 dt dx8 dt
AN US
792
M
791
ΛM ΛH and NM (0) > , µH µM ΛH then either the solutions enters Ω in finite time, or NH (t) approaches and NM (t) µH ΛM asymptotically. Hence, the region Ω attracts all solutions in R8+ . approaches µM Thus, the region Ω is positively-invariant. Furthermore, if NH (0) >
ED
790
λH =
801
bM (x7 αH + βH x8 ) bM (x3 αM + x4 βM ) , λM = , NH NH
and NH = x1 + x2 + x3 + x4 + x5 .
The Jacobian of system (C-1) evaluated at the disease-free equilibrium, E0 , is given by
40
ACCEPTED MANUSCRIPT
νH 0 0 J1 = 0 0 0 0
807
∗ −bM βH 0 β∗ (1−ε)bM VH H 0 SH 0 +(1−ε)V 0 β ∗ bM (SH H) H 0 SH
0
0
0
0
0
−k3
0
0
0
0 +(1−ε)V 0 α bM (SH H) H 0 SH
0
σH
−k4
0
0
0
0
γH
−µH
0
0
−µM
0
0
0
0
0
−
0 α bM SM M 0 SH
−
0 β bM SM M 0 SH
0
0 α bM SM M 0 SH
0 β bM SM M 0 SH
0
0
0
0
0
0
0
0
0
−k5
0
σM
−µM
Consider the case R0 = 1. Suppose that βH is chosen as a bifurcation parameter. Setting R0 = 1 and solving for βH using the expression for R0 gives o n 0 0 02 + (1 − ε)VH0 ] (σH βM + k4 αM ) [SH − b2M αH SM µM k3 k4 k5 SH 0 0 b2M σM SM (SH + (1 − ε)VH0 ) (σH βM + k4 αM )
.
PT
∗ The Jacobian matrix J1 evaluated at βH = βH has a simple zero eigenvalue with the other eigenvalues having negative real parts. Hence the centre manifold theorem can be applied ∗ . to analyze the dynamics of (C-1) near βH = βH
− Computing the right and the left eigenvector of J1 , denoted respectively by → w = [w1 , w2 , w3 , w4 , → − T w5 , w6 , w7 , w8 ] and v = [v1 , v2 , v3 , v4 , v5 , v6 , v7 , v8 ], we have 0 k3 k4 k5 w7 SH (k2 SH + (1 − ε)VH0 ωH ) , 0 0 bM S M (k1 k2 − νH ωH ) (SH + (1 − ε)VH0 ) (σH βM + k4 αM ) 0 0 k3 k4 k5 w7 SH (k1 (1 − ε)VH0 + νH SH ) , = 0 0 0 bM SM (k1 k2 − νH ωH ) (SH + (1 − ε)VH ) (σH βM + k4 αM ) 0 0 0 k4 k5 w7 SH k5 w7 σH SH k5 w7 γH σH SH , w = , w = , = 4 5 0 0 0 bM SM (σH βM + k4 αM ) bM SM (σH βM + k4 αM ) bM µ H S M (σH βM + k4 αM ) k5 w 7 w7 σM = − , w7 = w7 > 0, w8 = . µM µM
w1 = −
w2 w3 w6 809
−bM αH
−k2
AC
808
0
CE
805 806
0
0α (1−ε)bM VH H − 0 SH
∗ βH =
804
0
M
803
0
ED
802
ωH
CR IP T
−k1
AN US
and 41
ACCEPTED MANUSCRIPT
0 k3 k4 v3 SH k3 v3 βM , v7 = , v1 = v2 = v5 = v6 = 0, v3 = v3 > 0, v4 = 0 σH βM + k4 αM bM S M (σH βM + k4 αM ) 0 0 bM αH (SH − εVH + VH ) k3 k4 k5 SH . v8 = −v3 − 0 0 SH σM bM σM SM (σH βM + k4 αM )
811
Using the result from Theorem 4.1 by Castillo-Chavez and Song [21], we calculate the coefficients a and b as
CR IP T
810
n 2bM 0 0 v7 (w3 αM + w4 βM ) (w1 + w2 + w3 + w4 + w5 ) SM − w6 SH + VH0 0 0 2 (SH + VH ) o 0 (w2 ε + w3 + w4 + w5 ) . −v3 (w7 αH + w8 βH ) VH0 [w1 ε + (w3 + w4 + w5 ) (ε − 1)] − SH
The coefficient b is given as follows:
812
b=
AN US
a = −
v3 w7 bM σM ((1 − ε)νH + k2 ) . µM (νH + k2 )
815
D
M
814
The coefficient b is clearly positive; the presence of backward bifurcation in model (1) is determined by the sign of coefficient a.
813
Appendix D: Stability of model (2) DFE
818 819
The stability of E0 can be established using the next generation operator method on system (2). Taking EH (t), IH (t), EM (t), and IM (t) as the infected compartments and then using the notation in [61], the Jacobian F and V matrices for new infectious terms and the remaining transfer terms, respectively, are defined as: 0 0 bM α H S H bM βH SH 0 0 NH0 NH0 0 0 0 0 , F = 0 0 bM α M S M bM βM SM 0 0 0 S0 S H H 0 0 0 0 (µH + σH ) 0 0 0 −σH (γH + µH + δH ) 0 0 . V = 0 0 (µM + σM ) 0 0 0 −σM µM
CE
817
AC
816
PT
ED
The dengue model (2) has a disease free equilibrium (DFE). The DFE is obtained by setting the right-hand sides of the equations in the model (2) to zero, which is given by Λ ΛM H 0 0 , 0, 0, 0, , 0, 0 . E0 = (SH , 0, 0, 0, SM , 0, 0) = µH µM
42
ACCEPTED MANUSCRIPT
822 823 824
825 826
827
where, g1 = σH + µH , g2 = γH + µH + δH , g3 = σM + µM . ˜ 0 is the number of secondary infections in completely susceptible popuThe expression R lation due to infections from one introduced infectious individual with dengue. Further, using Theorem 2 in [61], the following result is established.
CR IP T
821
Therefore, the reproduction number is s 2 0 ˜ 0 = bM (αH µM + βH σM )(αM g2 + βM σH )SM , R 0 g1 g2 g3 µM SH
Lemma 4. The disease-free equilibrium (DFE) of the dengue model (2) is locally asymp˜ 0 > 1. totically stable (LAS) if R˜0 < 1 and unstable if R
E
828
Appendix E: Backward bifurcation analysis of the pre-intervention model (2)
AN US
820
Let,
∗ ∗ ∗ ∗ E1 = (SH , EH , IH , RH , SV∗ , EV∗ , IV∗ )
denote the endemic equilibrium of the dengue pre-intervention model (2) where λ∗H ΛH λ∗H ΛH σH λ∗H γH ΛH σH ΛH ∗ ∗ ∗ , E = , I = , R = , H H H µH + λ∗1 g1 (µH + λ∗1 ) g1 g2 (µH + λ1 ) g1 g2 µH (µH + λH ) ΛM λ∗M ΛM λ∗M ΛM σM ∗ ∗ = , E = , I = ; (E-1) M λ∗M + µM g3 (λ∗M + µM ) M g3 µM (λ∗2 + µM )
M
829
∗ SM
with λ∗1 =
833
Substituting the expression in (E-1) into (E-2), we obtain after some simplifications the quadratic polynomial
CE
832
(E-2)
z1 λ∗1 2 + z2 λ∗1 + z3 = 0,
(E-3)
AC
831
∗ ∗ ∗ ∗ ) βM ) αH + βH IM αM + IH bM (EM bM (EH ∗ , λ = . 2 ∗ ∗ NH NH
PT
830
ED
∗ SH =
where
z1 = g3 ΛH µM [g2 µH + σH (γH + µH )] [g2 µH (bM αM + µM ) + σH (bM µH βM + µM (γH + µH ))] , z2 = g1 g2 bM µ2H (g2 αM + σH βM ) [g3 ΛH µM − bM ΛM (αH µM + βH σM )] +2g1 g2 g3 ΛH µH µ2M [g2 µH + σH (γH + µH )] ,
834 835
e 2 ). z3 = g12 g22 g3 ΛH µ2H µ2M (1 − R 0
(E-4)
e 0 < 1 and negative when R e 0 > 1. Here, z1 is obviously positive and z3 positive whenever R Thus, the following result is established: 43
ACCEPTED MANUSCRIPT
836
837
838
Theorem 3. The dengue pre-intervention model (2) has: e 0 > 1; (i) a unique endemic equilibrium exists if and only if z3 < 0 and R
(ii) a unique endemic equilibrium exists if z2 < 0, and z3 = 0 or z22 − 4z1 z3 = 0; (iii) two endemic equilibria exists if z3 > 0, z2 < 0 and z22 − 4z1 z3 > 0;
840
(iv) otherwise no endemic equilibrium exists.
842 843 844 845
In Theorem 3, case (i) obviously shows that the model (2) has a unique endemic equilibrium e 0 > 1. Furthermore, case (iii) of the aforementioned Theorem 3, shows the whenever R e 0 < 1. To show this possibility, we set the possibility of backward bifurcation whenever R 2 e 0 , denoted by R e c and is polynomial z2 − 4z1 z3 = 0, and solve for the critical values of R given as ec = R
M
0.009 0.008
ED
0.01
PT
849
force of infection λ *H
848
z22 . 4z1 g12 g22 g3 ΛH µ2H µ2M
0.007 0.006
CE
847
1−
(E-5)
e 0 such that R ec < R e 0 < 1. Using So, backward bifurcation can occur for the values of R some of the parameter values given in Table 2 except for δH = 0.4097, µH = 0.0288, ΛH = 39, µM = 0.15 and ΛM = 30, we illustrate the possibility of backward bifurcation in Figure 15 for the dengue pre-intervention model (2).
AC
846
s
AN US
841
CR IP T
839
0.005 0.004 0.003 0.002 0.001 0
0.6
0.8
Rc
1
1.2
1.4
1.6
R0
Figure 15: Backward bifurcation curve for the dengue model (2).
44
ACCEPTED MANUSCRIPT
F
851
852 853 854 855
Global Asymptotic Stability of the DFE of the preintervention model (2) when δH = 0
Next, we consider the case when the disease-induced mortality rate is ignored. With proper medical care dengue fatality rates is below 1% [67]. In this case, we can either ignore the disease-induced death rate(δH = 0). Thus, we consider the case when δH = 0 in order to study the global stability of the disease-free equilibrium.
CR IP T
850
For dengue elimination to be independent of the initial sizes of the sub-populations of the model, the global asymptotic stability of the DFE must be established. This is considered next. Consider the feasible region ∗ ∗ Γ1 = {X ∈ Γ : SH ≤ SH , SM ≤ SM },
where , X = SH , EH , IH , RH , SM , EM , IM .
857
Lemma 5. The region Γ1 is positively invariant for model (2).
858
∗ Proof. It follows from the first equation of dengue model (2), (where SH =
AN US
856
ΛH ), so that µH
M
bM (αH EM (t) + βH IM (t))SH (t) dSH (t) = ΛH − − µH SH SB (t) dt NH (t) ≤ ΛH − µH SH (t)
ΛH − SH (t) µH
(F-1)
ED
≤ µH
PT
Thus,
CE
859
∗ = µH [SH − SH (t)].
∗ ∗ SH (t) ≤ SH − [SH − SH (0)]e−µH t .
∗ ∗ Thus, if SH (0) ≤ SH for all t ≥ 0, then SH (t) ≤ SH for all t ≥ 0.
861
∗ Similarly, it follows from the fifth equation of the dengue model (2) (where SM =
862
AC
860
that
ΛM ), µM
dSM (t) bM (αM EH (t) + βM IH (t))SM (t) = ΛM − − µM SM (t), dt NH (t) ≤ ΛM − µM SM (t), ≤ µM
ΛM − SM (t) , µM
∗ = µM [SM − SM (t)].
45
(F-2)
ACCEPTED MANUSCRIPT
867 868 869 870 871
872 873
874 875
876 877
∗ ∗ and SM (0) ≤ SM for all t ≥ 0, then SM (t) ≤ SM for all t ≥ 0.
Thus, in summary, it has been shown that the region Γ1 is positively invariant and attracts all solutions in R7+ for the dengue model (2).
CR IP T
866
ΛM µM
One of the sufficient conditions for the global asymptotic stability (GAS) of the diseasefree equilibrium (DFE) is a constant population level. This is typically the case when the disease-induced death rate is ignored (δH = 0). In this case the GAS result could be established using the second generation approach given in Castillo-Chavez et al. (2002) [20].
AN US
865
∗ Hence, if SM =
∗ ∗ SM (t) ≤ SM − [SM − SM (0)]e−µM t .
Theorem 4. The DFE, E0 , of model (2), is globally asymptotically stable (GAS) in Γ1 ˜ 0 ≤ 1. whenever R The above result shows that dengue will be eliminated from the community if the threshold ˜ 0 can be brought to a value less than unity. quantity R
M
864
Thus,
Theorem 4: The DFE, E0 , of model (2), is globally asymptotically stable (GAS) in Γ1 ˜ 0 ≤ 1 and δH = 0. whenever R
ED
863
880
Let X = (SH , RH , SM ) and Z = (EH , IH , EM , IM ) and group the system (2) into
CE
PT
879
Proof. To prove the global stability of the disease-free equilibrium, we will following the approach in [4, 3].
878
882
883
= F (X, 0)
dZ dt
= G(X, Z),
(F-3)
where F (X, 0) is the right hand side of S˙H , R˙H , and S˙M with EH = IB = IH = EM = ˙ . IM = 0 and G(X, Z) is the right hand side of E˙H , I˙H , E˙M and IM
AC
881
dX dt
Next, consider the reduced system: dSH dt dRH dt dSM dt
dX = F (X, 0) given as: dt
= ΛH − µH SH (t) = −µH RH (t)
(F-4)
= ΛM − µM SM (t). 46
ACCEPTED MANUSCRIPT
Let ∗
X =
∗ ∗ ∗ (SH , RH , SM )
=
ΛH ΛM , 0, . µH µM
886
To do this, solve the first and second equations of (F-4), this gives
SB (t) =
Z
t µH u
e
0
RH (t) =
Z
t µH u
e
887
Integrating SH (t) in (F-5) gives:
888
du + RH (0) e−µH t .
ΛH [1 − eµH t ] + SH (0)e−µH t . µH
ΛH and lim RH (t) = 0. t→∞ µH
ED
t=→∞
Similarly, solving for SM (t) in (F-4) gives SM (t) =
890
converges to
894
895
896
897 898
ΛM ΛM + e−µM t [SM (0) − ] which µM µM
CE
ΛM , as t → ∞. µM
These asymptotic dynamics are independent of initial conditions in Γ. Hence, the convergence of solutions of (F-4) is global in Γ1 . Next, following [20], we require G(X, Z) to satisfy the two stated conditions:
AC
893
PT
889
892
(F-6)
Taking the limit of SH (t) in (F-6) and RH (t) in (F-5) as t → ∞, we have lim SH (t) =
891
(F-5)
M
SH (t) =
ΛH du + SH (0) e−µH t
AN US
0
CR IP T
885
be an equilibrium of the reduced system (F-4), we show that X ∗ is a globally stable equilibrium in Γ1 .
884
(i). G(X, 0) = 0 and
b b (ii). G(X, Z) = DZ G(X ∗ , 0)Z − G(X, Z), G(X, Z) ≥ 0, ΛH ΛM ∗ where (X , 0) = , 0, , 0, 0, 0, 0 and DZ G(X ∗ , 0) is the Jacobian of G(X, Z) taken µH µM with respect to (EH , IH , EM , IM ) and evaluated at (X ∗ , 0), which is an M-matrix (the diagonal elements are nonnegative). Thus,
47
ACCEPTED MANUSCRIPT
−k1
0
∗ bM αH SH ∗ NH 0
σH −k2 DZ G(X ∗ , 0) = ∗ ∗ bM αM SM bM βM SM −k3 ∗ ∗ NH NH 0 0 σM ∗ SH (bM αH EM + bM βH IM ) 1− ∗ NH 0 b G(X, Z) = ∗ SM (b α E + b β M M H M M IH ) 1− ∗ NH 0
904
SH NH ∗ NH SH
∗ SM NH ∗ NH SM
AN US
Therefore, the disease-free equilibrium is globally asymptotically stable by the theorem in [20] (page 246).
M
903
ED
902
PT
901
−µM ∗
,
NH∗ SH Hence, if the human population is at equilibrium level, we have that 1 − ∗ > 0, SH NH ∗ SM NM b and 1 − > 0; thus, G(X, Z) ≥ 0. ∗ NM SM
CE
900
0
ΛH ΛH ΛM ∗ ∗ ∗ , NH∗ = and SM = . We have in Γ1 that, SH ≤ SH , and SM ≤ SM . µH µH µM Therefore, it follows that SH ≤ NH . ∗ where, SH =
AC
899
∗ bM βH SH ∗ NH 0
CR IP T
48
Optimal Control Strategies for Dengue Transmission in Pakistan F.B. Agusto, M.A. Khan PII: DOI: Reference:
S0025-5564(18)30453-X https://doi.org/10.1016/j.mbs.2018.09.007 MBS 8125
To appear in:
Mathematical Biosciences
Received date: Revised date: Accepted date:
29 July 2018 5 September 2018 6 September 2018
Please cite this article as: F.B. Agusto, M.A. Khan, Optimal Control Strategies for Dengue Transmission in Pakistan, Mathematical Biosciences (2018), doi: https://doi.org/10.1016/j.mbs.2018.09.007
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.
ACCEPTED MANUSCRIPT
Highlights • A mathematical model of dengue outbreak in Peshawar, Pakistan is considered. • Model formulations and its basic properties are presented.
CR IP T
• Sensitivity analysis and optimal control problem is formulated and discussed
AC
CE
PT
ED
M
AN US
• Control strategies are presented for disease elimination.
1
ACCEPTED MANUSCRIPT
CR IP T
Optimal Control Strategies for Dengue Transmission in Pakistan F. B. Agusto1 , M.A. Khan2∗ 2 1
Department of Ecology and Evolutionary Biology, University of Kansas, USA.
Department of Mathematics, City University of Science and Information Technology, 25000, KP, Pakistan.
6 7 8 9 10 11 12 13 14
15 16
17
18 19 20 21 22 23
M
5
ED
4
PT
3
This paper presents a deterministic model for dengue virus transmission. The model is parameterized using data from the 2017 dengue outbreak in Pakistan. We estimated the basic reproduction number (R0 ) without any interventions for the 2017 dengue outbreak in Peshawar district of Pakistan as R0 ≈ 2.64, the distribution of the reproduction number lies in the range R0 ∈ [1.21, 5.24] (with a mean R0 ≈ 2.64). Optimal control theory is then applied to investigate the optimal strategy for curtailing the spread of the disease using two time-dependent control variables determined from sensitivity analysis. These control variables are insecticide use and vaccination. The results show that the two controls avert the same number of infections in the district regardless of the weights on the costs this is due to the reciprocal relationship between the cost of insecticide use and vaccination. A strong reciprocal relationship exists between the use of insecticide and vaccination; as the cost of insecticide increases the use of vaccination increases. The use of insecticide on the other hand slightly increases when vaccination level decreases due to increase in cost.
CE
2
Abstract
Key Words: Dengue model, Vaccination, Basic reproduction number, Parameters estimation, Sensitivity analysis, optimal control strategies, numerical results.
AC
1
AN US
September 12, 2018
1
Introduction
Dengue fever is a mosquito-borne disease caused by the dengue virus. It has four different serotypes, i.e., DEN1, DEN2, DEN3, and DEN4, of the genus Flavivirus. The dengue is endemic in almost of 100 countries of the world, such as Africa, Eastern Mediterranean, the Americas, the subtropical regions, and especially the South-East Asia and Western Pacific regions are the most seriously affected [23, 55, 65]. After malaria, dengue has been recognized as the next deadliest mosquito-borne disease, with more than 100 millions ∗
Corresponding Author email: [email protected]
ACCEPTED MANUSCRIPT
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
CR IP T
30
The dengue virus transmission to the human is by the bites of mosquitoes (female mosquitoes, of type Aedes aegypti) [10], acquiring infection by taking the blood meal from an infected individual, then the infected mosquitoes pass the virus to healthy individuals (susceptible individuals). The recovery of an individual from one serotype leads to the permanently immune to it and temporarily or partially-immune to the others serotypes [63].
AN US
29
There is no specific effective treatment for the dengue fever until now but only the fluid replacement therapy at the beginning of infection diagnosis, and some traditional medicine (Brazilian traditional medical) in the dengue is treated with cat’s claw herb, but no such treatment observed [22, 24]. Besides treatment, there is no such effective and safe vaccine available for the dengue fever in the market. A report WHO shows some recent advancement about the vaccine development for the dengue fever in [64]. The first dengue vaccine, Dengvaxia (CYD-TDV) by Sanofi Pasteur, was first registered in Mexico in December, 2015. CYD-TDV is a live recombinant tetravalent dengue vaccine that has been evaluated as a 3-dose series on a 0/6/12 month schedule in Phase III clinical studies. It has been registered for use in individuals 9-45 years of age living in endemic areas [64]. Besides this, undergoing vaccine development includes, including live attenuated mono- and tetra-valent formulation, inactivated whole virus vaccines, and recombinant subunit vaccines[14, 13]. Any future dengue vaccine that perfectly protect against all serotype may or may not be possible, thus it is expected that any future dengue vaccine would be imperfect [10].
M
28
ED
27
PT
26
infections and some reports estimates about 390 million [12] and 20,000 deaths [55]. A study from 2012, shows that 3.9 billion people in 128 countries are at risk of infection with dengue viruses [17]. It is spread by many types of mosquitoes of the type Aedes, principally A. aegypti. Break bone fever or classical dengue fever relatively causes both mild morbidity and mortality, and the sufferers recover within a short span of time of about one or two weeks from the onset of fever [35, 36], some people develop hemorrhagic fever (DHF) or dengue shock syndrome (DSS) [37]. The disease severity drastically increased with mortality ranging from (5-15%) [35, 36, 38]. World Health Organization (WHO) report shows the occurrence of high number of cases of DHF annually worldwide [13, 14].
A number of mathematical models on the transmission dynamics of dengue fever have been proposed and analyzed, see [11, 23, 25, 26, 29, 30, 31, 42, 50, 49, 52, 57, 58, 60, 70] and the references therein. For instant, Rodrigues et al. [50] formulated and analyzed a vaccination models and optimal control strategies to dengue. First they studied the vaccination process by incorporating a new compartment in the model, and investigated the different ways of distributing the vaccines: pediatric and random mass vaccines, with distinct levels of efficacy and durability. They found that dengue eradication success depends on the type of vaccine as well as on the vaccination coverage. Simulating both epidemic and endemic scenarios, they showed using vaccine as a control variable, that vaccine when available can be a promising strategy to fight the disease. Rodrigues et al. [49] applied optimal control theory to study Dengue epidemics. Their proposed cost functional depends on the costs of medical treatment of the infected people, and the costs related to educational and sanitation campaigns. They showed efficient ways to decrease the number of infected mosquitoes and individuals in less time and with lower costs using current computational tools. Shim [57] studied the dynamics of dengue and cost-effectiveness of vaccine in the
CE
25
AC
24
3
ACCEPTED MANUSCRIPT
68 69 70
71 72 73 74 75 76
Philippines using an age-structured model. Their results indicate that the higher the cost of vaccination, the less cost-effective the dengue vaccination program. Furthermore, their results demonstrated that at a relatively low vaccine efficacy an age-targeted vaccination may be cost-effective provided the cost is sufficiently low vaccination. Motivated by the above studies, we develop a mathematical model to investigate the transmission dynamics of Dengue fever in Pakistan. We aim to study the impact of an imperfect vaccine in the bid to control dengue in Pakistan. Therefore, incorporated into the model is an imperfect vaccine. We parameterize the model using the 2017 dengue outbreak from Peshawar district of Khyber Pakhtunhwa Province in Pakistan. We then apply optimal control theory to study the best control strategies to curtail the disease in the community.
CR IP T
67
84
2
85 86 87 88 89 90 91 92
Model Formulation
M
82
The model features both human and mosquito populations. The human population is divided into five different classes, namely, susceptible individuals, SH (t), vaccinated, VH (t), exposed, EH (t), infected individuals, IH (t), and recovered individuals, RH (t), at any time t. Thus, the total human population denoted by NH (t) is given as NH (t) = SH (t) + VH (t) + EH (t) + IH (t) + RH (t). The population of mosquitoes is similarly divided into three different classes, namely, susceptible mosquitoes, SM (t), exposed mosquitoes, EM (t), and infected mosquitoes, IM (t), at any time t. So that, the total population of vector (mosquitoes) denoted by NM (t) is given by NM (t) = SM (t) + EM (t) + IM (t).
ED
81
PT
80
CE
79
AC
78
AN US
83
The paper is organized as follows: A model formulation and the basic mathematical properties of the model are presented in Section 2. The parameterization of the model is presented in Section 3. This is followed by a sensitivity analysis in Section 4. Information glean from the sensitivity analysis are used to identify the control variables applied in Section 5 on optimal control problem. Numerical results with detailed discussions are presented in Section 6. Finally, a detailed discussion and conclusion is presented in Section 7.
77
4
ACCEPTED MANUSCRIPT
The model that govern the system of differential equations is given by,
94
where λH =
= νH SH − (1 − ε)λH VH − (ωH + µH )VH , = λH [SH + (1 − ε)VH ] − (σH + µH )EH , = σH EH − (γH + µH + δH )IH , = γH IH − µH RH , = ΛM − λM SM − µM SM , = λM SM − (σM + µM )EM , = σM EM − µM IM ,
bM (αH EM + βH IM ) bM (αM EH + βM IH ) , and λM = . NH NH
AC
CE
PT
ED
M
95
= ΛH + ωH VH − λH SH − (νH + µH )SH ,
CR IP T
dSH dt dVH dt dEH dt dIH dt dRH dt dSM dt dEM dt dIM dt
AN US
93
Figure 1: Flow diagram of the dengue model (1). 96
The variables and their parameters are defined in Table 1.
5
(1)
ACCEPTED MANUSCRIPT
97
Analysis of the model
98
Basic qualitative properties
99
Positivity and boundedness of solutions
102
103 104 105
CR IP T
101
For the dengue model (1) to be epidemiologically meaningful, it is important to prove that all its state variables are non-negative for all time. In other words, solution of the system (1) with non-negative initial data will remain non-negative for all time t > 0. Lemma 1. Let the initial data F (0) ≥ 0 , where F (t) = (SH (t), VH (t), EH (t), IH (t), RH (t), SM (t), EM (t), IM (t)). Then the solutions F (t) of the dengue model (1) are non-negative for all t > 0. Furthermore lim sup NH (t) ≤ t→∞
ΛM ΛH , and lim sup NM (t) ≤ µH µM t→∞
with
AN US
100
NH (t) = SH (t) + VH (t) + EH (t) + IH (t) + RH (t), and NM = SM (t) + EM (t) + IM (t).
ED
Description Population of susceptible individuals Population of vaccinated individuals Population of exposed individuals Population of infected individuals Population of recovered individuals Population of susceptible mosquitoes Population of exposed mosquitoes Population of infected mosquitoes Description Vaccination rate Vaccine efficacy Vaccine waning rate Human recruitment rate Transmission probability from EM to SH . Transmission probability from IM to SH Natural mortality rate in humans Disease progression rate of infectious of exposed humans Mosquito biting rate Recovery rate of infected individuals in human compartments IH Disease related death rate Mosquito recruitment rate Transmission probability from EH to SM Transmission probability from IH to SM Natural death rate of mosquitoes Disease progression rate of infectious of exposed mosquitoes
CE
PT
Variable SH VH EH IH RH SM EM IM Parameter νH ε ωH ΛH αH βH µH σH bM γH δH ΛM αM βM µM σM
M
The proof of Lemma 1 is given in Appendix A.
AC
106
Table 1: Description of the parameters and variables for the dengue model (1). 6
ACCEPTED MANUSCRIPT
107
108 109
Invariant regions The dengue model (1) will be analyzed in a biologically-feasible region as follows. Consider the feasible region Ω = ΩH × ΩM ⊂ R5+ × R3+ with,
ΛH 4 , ΩH = (SH (t), VH (t), EH (t), IH (t), RH (t)) ∈ R+ : NH (t) ≤ µH
and
ΩM
112
113
ΛM 3 = (SM (t), EM (t), IM (t)) ∈ R+ : NM (t) ≤ . µM
Lemma 2. The region Ω ⊂ R8+ is positively-invariant for the model (1) with non-negative initial conditions in R8+ .
AN US
111
CR IP T
110
The prove of Lemma 2 is given in Appendix B.
116
Stability of disease-free equilibrium (DFE)
M
115
In the next section, the conditions for the existence and stability of the equilibria of the model (1) are stated.
114
ED
The dengue model (1) has a disease free equilibrium (DFE). The DFE is obtained by setting the right-hand sides of the equations in the model (1) to zero, which is given by
120 121
k2 ΛH Λ H νH ΛM , , 0, 0, 0, , 0, 0 . k1 k2 − νH ωH k1 k2 − νH ωH µM
CE
118 119
where k1 = νH + µH and k2 = ωH + µH . The stability of E0 can be established using the next generation operator method on system (1). Taking EH (t), IH (t), EM (t), and IM (t) as the infected compartments and then using the notation in [61], the Jacobian F and V matrices for new infectious terms and the remaining transfer terms, respectively, are defined as:
AC
117
PT
0 0 E0 = (SH , VH0 , 0, 0, 0, SM , 0, 0) =
0 0 F = 0 bM αM SM N0 H 0
0 0 0 bM βM SM NH0 0
0 bM αH [SH + (1 − ε)VH0 ] NH0 0
0 0
7
0 bM βH [SH + (1 − ε)VH0 ] NH0 0 , 0 0
ACCEPTED MANUSCRIPT
(σH + µH ) 0 0 0 −σH (γH + µH + δH ) 0 0 . V = 0 0 (σM + µM ) 0 0 0 −σM µM
129
130 131 132 133 134 135 136
137 138 139 140 141 142 143
CR IP T
AN US
128
Lemma 3. The disease-free equilibrium (DFE) of the dengue model (1) is locally asymptotically stable (LAS) if R0 < 1 and unstable if R0 > 1.
Backward Bifurcation Analysis
Disease transmission models typically undergo a simple transcritical bifurcation at R0 = 1 where there is an exchange of model’s stability as the model equilibria moves from the DFE to an endemic equilibrium. Some models such as vaccination models, on the otherhand exhibit backward bifurcation, where the stable DFE co-exists with a stable endemic equilibrium when the reproduction number is less than unity (see [7, 4, 8, 18, 27, 28, 34, 54, 53]). In a backward bifurcation setting, disease control is only feasible if R0 is reduced further to values below another sub-threshold less than unity.
M
127
ED
126
The expression R0 is the number of secondary infections in completely susceptible population due to infections from one introduced infectious individual with dengue. Further, using Theorem 2 in [61], the following result is established.
PT
124 125
where, ρ is the spectral radius and k3 = σH + µH , k4 = γH + µH + δH , k5 = σM + µM .
The important implication of this phenomenon on public health is that the classical requirement of having the reproduction number less than unity, although necessary, is no longer sufficient for disease control. This implies that effective disease control is dependent on the initial sizes of the sub-populations of the model. It is instructive, therefore, to explore whether or not the model (1) exhibits the phenomenon of backward bifurcation. In determining this possibility in the model (1), we use the Centre Manifold theory [19], as described in Theorem 4.1 by Castillo-Chavez and Song [21].
CE
123
Therefore, the reproduction number is s 0 0 b2M SM (αH µM + βH σM )(αM k2 + βM σH )[SH + (1 − ε)VH0 ] , R0 = ρ(F V −1 ) = k3 k4 k5 µM NH0
AC
122
145
Theorem 1. The dengue model (1) undergoes a backward bifurcation at R0 = 1 whenever the inequality a > 0 holds.
146
The prove of Theorem 1 is given in Appendix C.
144
8
ACCEPTED MANUSCRIPT
147
3
148
3.1
156 157 158 159 160 161
162 163 164 165 166 167 168
169 170 171 172
CR IP T
155
Between January - September 2017, the Ministry of National Health Services, Regulations and Coordination in Pakistan started reporting a high number of dengue fever cases in seven provinces of the country. From July - August 2017, a total of 2,199 cases had been reported by the ministry, Khyber Pakhtunkhwa (KPK) province had the highest number of reported laboratory confirmed cases (1,279 cases), this was followed by Sindh province with 793 cases, AJK province had the least number of reported cases (4 cases) [68]. However, by November 2017, a total of 24,807 laboratory confirmed cases with 69 associated deaths were reported in Khyber Pakhtunkhwa province [69].
AN US
154
A demographic break down of the 1,279 number of cases reported in July - August 2017 from the KPK province [68] show the 25-64 age group having the highest number of reported cases (668 cases, about 52.22% of the total), this was followed by the 15-24 age group with 371 reported cases (29.04%). A total of 199 cases were reported amongst the 0-14 age group (15.53%); the least affected age group was the 65+ group, with 41 reported cases (3.2%). 62.9% of the cases (804) are male while females represented 37.1% (475) of the reported cases.
M
153
ED
152
In order to parametrize our model, we obtained from the dengue response unit in Peshawar the 2017 dengue outbreak data from various hospitals in the Peshawar district, namely, Lady Reading Hospital (LRH), Hayatabad Medical Complex (HMC), and Khyber Teaching Hospital (KTH), see Figure 2.
PT
151
The dengue is a frequent epidemic in Pakistan. In 2011 and 2013 dengue outbreak occurred in Lahore, Punjab and in Swat Khyber Pakhtunkha district respectively. The 2011 dengue outbreak in Lahore, Punjab lead to 300 death while the 2013 outbreak in Swat Khyber Pakhtunkha district caused a total of 57 death cases (see, [39, 40] and the references therein).
CE
150
The 2017 dengue outbreak data
AC
149
Parameter Estimation
9
ACCEPTED MANUSCRIPT
5000
0
Jul 24 17
Sep 12 17 Time (days)
CR IP T
10000
Nov 01 17
AN US
Reported numberof cases
15000
Figure 2: Number of new cases in the 2017 dengue outbreak in Pakistan. The data were obtained from Peshawar district hospitals.
175
In order to parametrize the dengue model (1), first, we study the worst-case scenario where public health interventions such as vaccination are not implemented in the community. In the absence of such interventions, the model (1) reduces to the following basic model: dSH dt dEH dt dIH dt dRH dt dSM dt dEM dt dIM dt
bM (αH EM + βH IM )SH − µH SH , NH bM (αH EM + βH IM )SH = − (σH + µH )EH , NH = ΛH −
AC
CE
PT
176
Pre-intervention model
M
174
3.2
ED
173
177
= σH EH − (γH + µH + δH )IH , = γH IH − µH RH , bM (αM EH + βM IH )SM − µM SM , NH bM (αM EH + βM IH )SM = − (σM + µM )EM , NH = ΛM −
= σM EM − µM IM .
The dengue model (2) is bounded and positively invariant in the region Ω = ΩH × ΩM ⊂ R4+ × R3+ 10
(2)
ACCEPTED MANUSCRIPT
178
with, ΩH = and
180 181 182 183
184 185 186
ΛH : NH (t) ≤ , µH
ΛM 3 = (SM (t), EM (t), IM (t)) ∈ R+ : NM (t) ≤ . µM
In Appendix D, we state the conditions for the stability of the disease-free equilibrium ˜ 0 . Further analysis about backward of model (2) based on the reproduction number, R bifurcation of model (2) are given in Appendix E. Backward bifurcation is the phenomenon in which two stable equilibria (in this case disease-free and endemic equilibria) co-exist when the reproduction number is less than one.
AN US
179
(SH (t), EH (t), IH (t), RH (t)) ∈
R4+
CR IP T
ΩM
A number of vector-borne diseases models [4, 5, 15, 59] have shown that in the absence of disease related death, the phenomenon of backward bifurcation can be lost and the disease-free equilibrium will become globally stable.
192
3.3
194 195 196 197 198 199 200
201 202 203 204 205
ED
193
Data fitting
In order to parametrize dengue model (2), we obtained some of the parameter values from literature (see Table 2), others were estimated or fitted to the Peshawar dengue outbreak data obtained from Peshawar district hospital using the classical least-squares method (see Figure 3). For instance, the demographic parameter, µH , is estimated as µH = 1/67.7 per year, where 67.7 years is the average lifespan in Pakistan [62]. The other demographic parameter, ΛH , is then estimated as follows. Since the total urban population of Peshawar district as at 2017 was 1, 970, 042 [47], we assumed that ΛH /µH = 1, 970, 042 is the limiting total human population in the absence of the disease, so that ΛH = 29, 099.59 per year.
PT
190
CE
189
AC
188
M
191
There is no specific treatment for dengue or severe dengue [67]; however with early detection and access to proper medical care fatality rates can be lower to rate below 1% [67] or may vary between 2%-5% [56]. But, when left untreated, the mortality rate can be as high as 20% [56]. Hence, in Appendix F, we ignore the disease induce death rate and show that the disease-free equilibrium corresponding to model (2) is globally stable.
187
Thus, we obtained after fitting the model (2) to the 2017 dengue outbreak data the following estimates for bM = 0.8759, αH = 0.8728, βH = 0.8535, σH = 0.01298, γH = 0.01837, αM = 0.8593, βM = 0.7646, σM = 0.005728, δH = 0.000015857. Consequently, ˜ 0 (stated in using these parameter estimates, we obtained, in this study, that the value of R ¯ 0 ≈ 2.65. Appendix D) for the 2017 dengue outbreak in Peshawar district of Pakistan is R
11
15000
Model Data
CR IP T
10000
5000
0
AN US
Cumulative number of infected cases
ACCEPTED MANUSCRIPT
Jul 24 17
Sep 12 17 Time (days)
Nov 01 17
M
Figure 3: Data fitting of the number of cases using model (1), for the 2017 dengue outbreaks in Pakistan. The parameters fitted are given as bM = 0.8759, αH = 0.8728, βH = 0.8535, σH = 0.01298, γH = 0.01837, αM = 0.8593, βM = 0.7646, σM = 0.005728, δH = 0.000015857. Baseline value 29,099.59 0.8728 0.8535 1/67.7 0.01298 0.8759 0.01837 0.000015857 0.3 0.8593 0.7646 1/42 0.005728
AC
CE
PT
ED
Parameter ΛH αH βH µH σH bM γH δH ΛM αM βM µM σM
Reference Estimated Fitted Fitted [62] Fitted Fitted Fitted Fitted [43] Fitted Fitted [43] Fitted
Table 2: Numerical values for the parameters of the dengue model (1). 206
207 208
3.4
Quantify the expected disease burden
In order to quantify the expected burden of the disease in the country (under the worst˜ 0 is generated, using the parameter values and case scenario), the distribution profile of R 12
ACCEPTED MANUSCRIPT
211 212 213 214
ranges in Table 2. The results obtained, depicted in Figure 4, show the distribution of the ˜ 0 ∈ [1.21, 5.24] (with a mean R ˜ 0 ≈ 2.71, suggesting reproduction number in the range R the potential for a larger dengue outbreaks, where one infected case infects, on average, about 2.71 others). Furthermore, it is worth noting that for the worse-case scenario, the ˜ 0 values are concentrated in the interval (2.2036, 3.1518), central 50% of the generated R ˜ with the median value of R0 = 2.64 been close to the mean.
CR IP T
209 210
5 4.5
3.5 3 2.5 2 1.5 200
300
400
500
600
700
800
900
1000
M
100
AN US
R0
4
Runs
ED
˜ 0 , for the models (2). Parameter Figure 4: The Scatterplots of the reproduction number, R values (baseline) and ranges used are as given in Table 2.
218
219 220
221
222 223 224 225 226 227 228
PT
217
CE
216
˜ 0 for the 2017 Hence, to summarize, we have obtained in this study, that the value of R ˜ dengue outbreak in Peshawar district of Pakistan is R0 ≈ 2.65 and the distribution of the ˜ 0 ∈ [1.2127, 5.2426] (with a mean R ˜ 0 ≈ 2.7052) for reproduction number is in the range R the estimated worse case scenario. In the next section, we apply optimal control theory to determine the optimal control strategy that will effectively curtail the spread of dengue in the community.
AC
215
4
Sensitivity analysis
The outputs of deterministic models are governed by the model input parameters, which may exhibit some uncertainty in their determination or selection. We employed a global sensitivity analysis to assess the impact of uncertainty and the sensitivity of the outcomes of the numerical simulations to variations in each parameter of the model (1) using Latin Hypercube Sampling (LHS) and partial rank correlation coefficients (PRCC). LHS is a stratified sampling without replacement technique which allows for an efficient analysis of parameter variations across simultaneous uncertainty ranges in each parameter [16, 44, 45, 13
ACCEPTED MANUSCRIPT
237 238 239 240 241 242 243
244 245 246 247 248 249 250 251
252 253 254 255 256 257 258
CR IP T
236
Figure 5 depict the PRCC values for each parameter of the model using the reproduction ˜ 0 and R0 , as the response functions. Parameters with the highest PRCC values numbers, R have the largest impact on the reproduction numbers. Therefore, the key parameters ˜ 0 and R0 ) are separated into those that decrease influencing there production numbers (R ˜ 0 and R0 when increased (bars extending to the left for negative PRCC values) and those R ˜ 0 and R0 to increase when increased (bars extending to the right for positive that causes R PRCC values).
AN US
235
From Figures 5(a) and 5(b), it follows that the parameters that have the most negative influence on dengue transmission dynamics are human recruitment rate (ΛH ), disease progression rate in humans (σH ), human recovery rate (γH ), mosquito natural death rate ˜0 (µM ), and vaccine efficacy (ε). While the parameters with the most positive impact on R and R0 are the transmission probability per contact of infectious mosquitoes with susceptible humans (αH and βH ), mosquito recruitment rate (ΛM ), mosquito biting rate (bM ), transmission probability per contact of susceptible mosquitoes with infectious humans (αM and βM ), disease progression rate in mosquitoes (σM ), and vaccine waning rate (ωH ).
M
233 234
ED
232
Identification of these key parameters is important to the formulation of effective control strategies necessary for combating the spread of disease in the community. In particular, the results of this sensitivity analysis suggest that a strategy that reduces the parameters with positive PRCC values (i.e., βH , αH , ΛM , bM , βM , αM , σH , ωH ) will adequately reduce the spread of dengue in the community. Furthermore, a strategy that increases the parameters with negative PRCC values (i.e., ΛH , σH , γH , µM , ε) will be effective in curtailing the spread of dengue in the community.
PT
231
51]. PRCC measures the strength of the relationship between the model outcome and the parameters, stating the degree of the effect that each parameter has on the outcome [16, 44, 45, 51]. Thus, sensitivity analysis determines the parameters with the most significant impact on the outcome of the numerical simulations of the model [16, 44, 46]. To generate the LHS matrices, we assume that all the model parameters are uniformly distributed. Then a total of 1,000 simulations of the models per LHS run were carried out, using the baseline values tabulated in Table 2 and the ranges as 20% from the baseline values (in either direction).
CE
230
AC
229
14
ǫ ωH ν µMH σM β αMM b ΛMM δ µHH γ σHH β αHH ΛH
µM σM βM αM bM ΛM δH µH γH σH βH αH ΛH
-1
-0.5
0
0.5
1
-1
-0.5
0
0.5
1
(b)
AN US
(a)
CR IP T
ACCEPTED MANUSCRIPT
Figure 5: PRCC values for the dengue models (2) and (1), using as response functions: (a) The
˜ 0 . (b) The reproduction number R0 . Parameter values (ranges) used are reproduction number R as given in Table 2.
263 264 265 266 267 268
269 270 271 272 273 274 275
276 277 278 279 280
M
ED
262
The result of the sensitivity analysis suggest that a control strategy that reduces the transmission probability per contact of infectious mosquitoes with susceptible humans (αH and βH ), mosquito recruitment rate (ΛM ), mosquito biting rate (bM ), transmission probability per contact of susceptible mosquitoes with infectious humans (αM and βM ), disease progression rate in mosquitoes (σM ), and vaccine waning rate (ωH ) will adequately reduce the spread of dengue in the community. Furthermore, a strategy that increases the human recruitment rate (ΛH ), disease progression rate in humans (σH ), human recovery rate (γH ), mosquito natural death rate (µM ), and vaccine efficacy (ε) will be effective in curtailing the spread of dengue in the community.
PT
261
The optimal control problems
CE
260
5
Therefore, we introduce into the transmission model (1) a time dependent control variable uV (t), representing the vaccination rate which was previously taken as a constant parameter (νH ). Although, this parameter does not significantly impact the reproduction number (R0 ), we nevertheless consider this parameter as time dependent variable since the vaccine efficacy (ε) have a significant impact of the (R0 ). We also consider as time dependent control variable uM (t) that impact the mosquito population since mosquito death rate µM has a strong significant impact on R0 .
AC
259
The mosquito biting rate bM and transmission probabilities αH , βH , αM and βM , are another set of candidate parameters to use as a time dependent controls. Using them as control measures however, will require the use of personal protection which is hard to implement since the mosquito species causing dengue only bite during the day. So, our control efforts using control variables uV (t) and uM (t) will indirectly impact these 15
ACCEPTED MANUSCRIPT
parameters. Thus, model (1) with the time dependent control variables uV (t) and uM (t) becomes:
283
bM (αH EM + βH IM )SH − (uV (t) + µH )SH NH (1 − ε)bM (αH EM + βH IM )VH = uV (t)SH − − (ωH + µH )VH NH bM (αH EM + βH IM )[SH + (1 − ε)VH ] = − (σH + µH )EH NH = ΛH + ωH VH −
= σH EH − (γH + µH + δH )IH = γH IH − µH RH ,
bM (αM EH + βM IH )SM − µM SM − uM (t)SM NH bM (αM EH + βM IH )SM = − (σM + µM )EM − uM (t)EM NH = ΛM −
= σM EM − µM IM − uM (t)IM .
Our goal is to minimize the cost function defined as Z
0
289 290 291 292
293
+
A4 u2M (t)
dt
And we seek to find optimal controls u∗V , u∗M , such that J(u∗V , u∗M ) = min{J(uV , uM )} U
294
(4)
subject to the differential equations (3), where tf is the final time. This performance specification involves minimizing the numbers of exposed and infected individuals, along with the cost of applying the controls (uV (t), uM (t)). In this paper, a quadratic objective functional is implemented for measuring the control cost, since the costs of the intervention are nonlinear. This assumption is based on the fact that there are no linear relationship between the effects of intervention and the cost of intervention of the infective populations, such quadratic costs have been frequently used [1, 2, 9]. The coefficients, Ai , i = 1, · · · , 4, are balancing cost factors. The controls (uV (t), uM (t)) is a bounded Lebesgue integrable functions [1, 2, 9].
PT
288
A1 EH + A2 IH +
A3 u2V (t)
CE
287
AC
285 286
tf
ED
J(uV , uM ) = 284
(3)
CR IP T
dSH dt dVH dt dEH dt dIH dt dRH dt dSM dt dEM dt dIM dt
AN US
282
M
281
where the control set, U = {(uV , uM ) : [0, tf ] → [0, 1], (uV (t), uM (t)), is Lebesgue measurable}.
16
ACCEPTED MANUSCRIPT
295
296 297 298 299
The necessary conditions that an optimal control quintuple must satisfy come from the Pontryagin’s Maximum Principle [48]. This principle converts (3) and (4) into a problem of minimizing pointwise a Hamiltonian H, with respect to the controls (uV (t) and uM (t)). First we formulate the Hamiltonian from the cost functional (4) and the governing dynamics (3) to obtain the optimality conditions.
CR IP T
300
Characterization of optimal controls
305 306 307
308
PT
304
∗ ∗ ∗ Theorem 2. Given an optimal control quintuple (u∗V , u∗M ) and solutions SH , VH∗ , EH , IH , ∗ ∗ ∗ ∗ ∗ ∗ RH , SM , EM , IM of the corresponding state system (3) that minimizes J(uV , uM ) over U. Then there exists adjoint variables λSH , λVH , λEH , λIH , λRH , λSM , λEM , λIM , satisfying
CE
303
where λSH , λVH , λEH , λIH , λRH , λSM , λEM , λIM are the associated adjoints for the states SH , VH , EH , IH , RH , SM , EM , IM . The system of adjoint equations is found by taking the appropriate partial derivatives of the Hamiltonian (5) with respect to the associated state and control variables.
AC
301 302
ED
M
AN US
H = A1 EH + A2 IH + A3 u2V + A4 u2M i h bM (αH EM + βH IM )SH − (uV + µH )SH + λSH ΛH + ωH VH − NH n o (1 − ε)bM (αH EM + βH IM )VH + λVH uV SH − − (ωH + µH )VH NH o n b (α E + β I )[S + (1 − ε)V ] M H M H M H H − (σH + µH )EH + λEH NH h i h i + λIH σH EH − (γH + µH + δH )IH + λRH γH IH − µH RH h i bM (αM EH + βM IH )SM + λSM ΛM − − µM SM − uM (t)SM NH h b (α E + β I )S i M M H M H M + λEM − (σM + µM )EM − uM (t)EM NH h i + λIM σM EM − µM IM − uM (t)IM
−
∂H dλi = dt ∂i
and with transversality conditions λi (tf ) = 0, where i = SH , VH , EH , IH , RH , SM , EM , IM .
309
The optimality conditions is given as ∂H = 0, ∂uj 310
(5)
j = 1, 2.
Furthermore, the control (u∗V , u∗M ) is given as 17
(6)
ACCEPTED MANUSCRIPT
u∗M
313 314
−
dλSH dt
dλRH − dt −
dλSM dt
dλI − M dt
=
∂H , ∂RH
λRH (tf ) = 0,
=
∂H , ∂SM
λSM (tf ) = 0,
∂H , ∂IM
λIM (tf ) = 0,
···
··· =
ED
= (λSH − λVH )uV + µH λSH + λ∗H +λ∗M
dλVH dt
PT
dλSH dt
∗ SM (λEM − λSM ), NH∗
CE
316
λSH (tf ) = 0,
evaluated at the optimal controls and corresponding state variables, results in the stated adjoint system (5) and (6) given by
+(λVH − λEH )(1 − ε)λ∗H
dλEH dt
∗ ∗ ) (NH∗ − SH ∗ VH (λ − λ ) + (1 − ε)λ (λEH − λVH ) SH EH H NH∗ NH∗
= ωH (λVH − λSH ) + µH λVH + (λEH − λSH )λ∗H
AC
315
∂H , ∂SH
=
AN US
312
Proof. The existence of an optimal control is guaranteed using the result by Fleming and Rishel [33]. Thus, the differential equations governing the adjoint variables are obtained by the differentiation of the Hamiltonian function, evaluated at the optimal controls. Thus, the adjoint system can be written as:
M
311
(7)
CR IP T
u∗V
" # ∗ SH (λ∗SH − λ∗VH ) , = min 1, max 0, 2A3 " # ∗ ∗ ∗ λ∗SM SM + λ∗EM EM + λ∗IM IM = min 1, max 0, . 2A4
= (λEH − λSH )
∗ ∗ SH ∗ SM + (λ − λ )λ E S M M M NH∗ NH∗
(NH∗ − VH∗ ) , NH∗
∗ ∗ ∗ ∗ λ∗H SH ∗ VH ∗ (NH − EH ) + (λ − λ )(1 − ε)λ + (λ − λ )b α S EH VH SM EM M M M H NH∗ NH∗ NH∗ 2
+σH (λEH − λIH ) + µH λEH − A1 ,
dλIH dt
= (λEH − λSH )
∗ ∗ ∗ λ∗H SH λ∗H VH∗ ∗ (NH − IH ) + (λ − λ )(1 − ε) + (λ − λ )b β S EH VH SM EM M M M NH∗ NH∗ NH∗ 2
+γH (λIH − λRH ) + (µH + δH )λIH − A2 , dλRH dt
= (λEH − λSH )
∗ λ∗H ∗ λ∗H ∗ ∗ SM ) )λ S + (λ − λ V (1 − ε) + (λ − λ + µH λRH , EH VH EM SM M NH∗ H NH∗ H NH∗
18
ACCEPTED MANUSCRIPT
dλSM dt dλEM dt
= (λSM − λEM )λ∗M + (µM + uM )λSM , = (λSH − λEH )bM
αH αH ∗ S + (1 − ε)(λVH − λEH )bM ∗ VH∗ + σM (λEM − λIM ) NH∗ H NH
+(µM + uM )λEM ,
319
Furthermore, differentiating the Hamiltonian function with respect to the control variables in the interior of the control set and then solving for controls (u∗V , u∗M ) result in the optimality conditions given as ∂H ∂uV ∂H ∂uM
320
∗ = 2A3 uV − SH (λ∗SH − λ∗VH ) = 0,
∗ ∗ ∗ = 2A4 uM − λ∗SM SM + λ∗EM EM + λ∗IM IM = 0.
Solving for u∗V and u∗M , we have,
u∗M
326 327 328 329 330
331 332
PT
Remark 5.1. Due to the a priori boundedness of the state and adjoint functions and the resulting Lipschitz structure of the ODE’s, the uniqueness of the optimal control for small time (tf ) was obtained. The uniqueness of the optimal control quintuple follows from the uniqueness of the optimality system, which consists of (3) and (5), (6) with characterization (7). The restriction on the length of the time interval is to guarantee the uniqueness of the optimality system, the smallness in the length of time is due to the opposite time orientations of (3), (5), and (6); the state problem has initial values, and the adjoint problem has final values. This restriction is very common in control problems (see [1, 2, 9]).
CE
324 325
Using the bounds on the controls, the characterization (7) can be derived.
AC
322 323
∗ (λ∗SH − λ∗VH ) SH , 2A3 ∗ ∗ ∗ + λ∗IM IM + λ∗EM EM λ∗SM SM . = 2A4
=
ED
u∗V
321
AN US
318
bM β H ∗ bM β H SH + (λVH − λEH ) ∗ VH∗ (1 − ε) + (µM + uM )λIM . ∗ NH NH
M
317
= (λSH − λEH )
CR IP T
dλIM dt
Next, we discuss the numerical solutions of the optimality system, the corresponding optimal control and the interpretations from various cases.
19
ACCEPTED MANUSCRIPT
340 341
342 343 344 345 346 347 348 349 350 351 352 353 354 355 356
357 358 359 360 361 362 363
CR IP T
339
Thus, to illustrate the optimal control strategies, the following values were taken for the initial conditions to mimic the Peshawar District of Pakistan. According to the Pakistan Bureau of Statistics the population of Peshawar District in the 2017 census is estimated as 1,970,042 [47]. We take the total initial population as NH (0) = 1970 per 1000 of the total population. We assume that initial population of vaccinated individuals VH (0) = 10; the initial population of exposed is given as EH (0) = 10. The initial infected is as given in the data is I(0) = 1. We assume there are no recovered individuals initially, so R(0) = 0. Hence, the initial susceptible population is given as S(0) = NH (0) − V (0) − EH (0) − IH (0) − RH (0) = 1959. We assume the initial population for the mosquitoes are given as SM (0) = 1338, EM (0) = 10, IM (0) = 152. For the weight factors we choose A1 = 1.0, A2 = 1.0, A3 = 0.0001, and A4 = 0.01. It should be pointed out that the weights in the simulations here are only of theoretical sense to illustrate the control strategies proposed in this paper. The reproduction number in the absence of intervention is given ˜ 0 ≈ 2.65, this was obtained from using the estimated and fitted parameter values in as R Table 2. This value indicate that the disease is epidemic in the district.
AN US
338
M
337
ED
336
The following algorithm was used to compute the optimal controls and state values using a Runge-Kutta method of the fourth order. First, an initial estimate for the control quintuple is made. Then the state variables are solved forward in time using the dynamics (3). The results obtained for the state variables are plugged into the adjoint equations (5). These adjoint equations with given final conditions (7) are then solved backward in time, employing the backward fourth order Runge-Kutta method. Both the state and adjoint values are then used to update the control, and the process is repeated until the current state, adjoint, and controls values converge sufficiently [41].
The results of the optimal control simulations of the dengue model (3) are depicted in Figure 6. With optimal control, most of the susceptible populace are either protected or vaccinated from the virus (see figures 6(a) and 6(b)) thereby leading to fewer individuals at the risk of infection (see Figure 6(c)) and eventually resulting in fewer individuals recovering from the infection as shown in Figure 6(d). This is not the case in the absence of the control strategy. For the mosquitoes, there are fewer of them with the control than in the absence of control as observed in Figures 6(e) and 6(f).
PT
335
Numerical Illustration
CE
334
6
AC
333
20
1500
1000
500
0
20
40
60
80
1500
without control with control
1000
500
100
0
20
Time (days)
60
80
ED
1000
PT
800 600 400 200 0
20
40
60
80
100
80
100
400
200
0
20
40
60
Time (days)
(d)
1200
without control with control
1000 800 600 400 200 0
100
20
40
60
80
100
Time (days)
Time (days)
AC
80
600
0
without control with control
CE
Susceptible mosquitoes
(c)
1200
without control with control
M
Time (days)
100
Recovered individuals
500
800
AN US
1000
Infectious mosquitoes
Infectious individuals
1500
40
60
(b)
without control with control
20
40
Time (days)
(a)
0
CR IP T
without control with control
Vaccinated individuals
Susceptible individuals
ACCEPTED MANUSCRIPT
(e)
(f)
Figure 6: Simulation results of model (3) with and without controls. (a) Susceptible individuals; (b) Vaccinated individuals; (c) Infectious individuals (exposed and Infected); and (d) Recovered individuals; (e) Susceptible mosquitoes; and (f) Infectious mosquitoes (exposed and Infected). Using A1 = 1.0, A2 = 1.0, A3 = 0.0001, A4 = 0.01.
364 365 366
The corresponding time-dependent controls (uV (t) and uM (t)) are depicted in Figure 7. The time-dependent controls uV (t) is at the upper bound of unity for about 10 days before gradually decreasing to the lower bound at the end of the simulation period; while the 21
ACCEPTED MANUSCRIPT
368 369 370
control uM (t) starts at the upper bound for over 60 days and slowly decrease to the lower bound at the end of the simulation period. These results suggest that to prevent an outbreak, vaccination control should start at a relatively high levels and while the vector control is maintain at high levels for a long time.
CR IP T
367
uV
1
uM
0.6
AN US
control
0.8
0.4
0
0
20
M
0.2
40
60
80
100
ED
Time (days)
373 374
375
CE
372
Varying cost weight on controls uV (t) and uM (t) In this section, we investigate the effect of varying the weights A3 , and A4 on the controls uV (t) and uM (t). These coefficients, Ai , i = 3, 4, are balancing cost factors on the cost of implementing the control strategies.
AC
371
PT
Figure 7: Simulation results of the control profile of model (3). Using A1 = 1.0, A2 = 1.0, A3 = 0.0001, A4 = 0.01.
Varying the weight A3 on control uV (t)
381
First we vary the weight A3 on the control uV (t) representing the use of vaccination. We use the weights A3 = 0.001, A3 = 0.01, A3 = 5.0 and A3 = 10.0 and keep constant A1 = A2 = 10.0 and A4 = 1.00. These costs correspond to the case with very cheap, cheap, expensive and very expensive vaccination strategies. The solution profile are depicted in Figure 8. The solution profile for the other variables are similar to what we observed in Figures 6 and are therefore not show here.
382
Figure 8 show that as the weight A3 increases from low to high (i.e., from cheap to very
376 377 378 379 380
22
ACCEPTED MANUSCRIPT
384 385 386 387
expensive cost), control uV (t) decreases from the higher levels to very low levels. The insecticide control uM (t) remain high, starting at the upper bound for about 10 days. The control uM slightly increases as weights A3 increases. This slight adjustment in control uM (t) help to keep the community protected and the infection low as observed in the sub-populations depicted in Figure 6.
uV
1
CR IP T
383
uV
1
uM
uM
0.8
control
control
0.8 0.6 0.4 0.2
20
40
60
80
100
0
uV
M
0.8
ED
0.6
20
PT
0
40
40
60
80
100
60
80
uV
1
uM
0.2
20
(b)
1
0.4
0
Time (days)
(a)
control
AN US
0
Time (days)
0
0.4 0.2
uM
0.8
control
0
0.6
0.6 0.4 0.2 0
100
Time (days)
0
20
40
60
80
100
Time (days) (d)
CE
(c)
Figure 8: Simulation results of the control profile of model (3) using A1 = A2 = 10.0, A4 = 1.00
388
AC
with cost weights (a) A3 = 0.001; (b) A3 = 0.01; (c) A3 = 5.0; and (d) A3 = 10.0.
Varying the weight A4 on control uM (t)
394
Next, we vary the weight A4 on control uM (t) which corresponds to insecticide use. We keep constant the weights A1 = A2 = 10.0, A3 = 0.001 and use the weights A4 = 0.01, A4 = 0.10, A4 = 5.0, and A4 = 10.0. These costs as with vaccination correspond to the case with very cheap, cheap, expensive and very expensive insecticide use. The control solution profile are depicted in Figure 9. The solution profile for the other variables are similar to what we observed in Figures 6 and are also not show here.
395
At low A4 values (i.e., A4 = 0.01) which corresponds to really cheap vaccine, the control
389 390 391 392 393
23
ACCEPTED MANUSCRIPT
397 398 399
400 401 402 403 404
uM (t) is usually at the upper bound and remain at the upper bound longer than control uV (t), the vaccination rate. However, as A4 gets bigger, the level of the control uM (t) falls below the control uV (t) as the level of control uV starts to rise. Indicating the existence of a reciprocal type relation between the weight A4 and control uV (t) [32]. This reciprocal relation between the weight A4 and control uV (t) is strongerthan the reciprocal relationship between A3 and control uM (t) (see Figure 8). Note, this relationship help to keep the community protected and the infection low as observed in Figure 6. It is possible to get a reciprocal relationship with a different strength depending on the parameter set used.
1
control
uM
0.4 0.2
0.6 0.4 0.2
0
20
40
60
80
100
0
ED
(a)
1 0.8
0.2
CE
0.4
AC
40
60
80
60
80
100
80
100
(b)
0.6 0.4 0.2
uM
20
40
Time (days)
0.8
uV
0
20
1
PT
0.6
0
M
Time (days)
control
AN US
0.8
control
control
uM
uV
0.6
0
uV
1
0.8
0
CR IP T
396
0
100
Time (days)
uV uM
0
20
40
60
Time (days)
(c)
(d)
Figure 9: Simulation results of the control profile of model (3) using A3 = 0.0001 with cost weights (a) A4 = 0.001; (b) A4 = 0.01; (c) A4 = 10.0; and (d) A4 = 13.0.
405 406 407
Next, we compare in Figure 10 the impart of these weights on the infection averted ratio (IAR), average cost-effectiveness ratio (ACER), infection averted, and total cost discussed in details in the next section below.
24
ACCEPTED MANUSCRIPT
408
Benefits gained from varying controls uV (t) and uM (t)
413
Infection averted ratio
414
The infection averted ratio (IAR) is stated as:
410 411
IAR =
CR IP T
412
Furthermore, we will investigate the associated benefits gained from implementing these controls while varying Ai , i = 3, 4, using cost-effectiveness analysis [2, 6]. To achieve this, we will consider three approaches, the infection averted ratio (IAR), and the total cost obtain from the objective functional in (4).
409
Number of infection averted . Number of recovered
(8)
418
Average cost-effectiveness ratio (ACER)
416
419 420
Next, we consider the average cost-effectiveness ratio (ACER) which deals with a single intervention, evaluating it against the no intervention baseline option. ACER is calculated as Total cost produced by the intervention . (9) ACER = Total number of infection averted
425 426
427 428 429 430 431
PT
424
The results depicted in Figures 10(c) and 10(d) are expected, lower weights are more cost-effective than higher weights, surprisingly, insecticide led to lower cost compare to vaccination at higher weights. This is due to the strong reciprocal effect we observed in Figures 9(c) and 9(d) when A4 was been varied and the control uV adjusted to compensate for the low levels of uM use.
CE
423
The two weights averted approximately the same number of infection regardless of the cost (see Figure 10(a), and 10(b)). This is due to the reciprocal effect we observed in Figures 8 and 9, as the two controls uV and uM adjust to compensate for the low levels in one when the other is high. Thus, the two controls exhibit this sea-saw type behavior, thus helping to keep the level of infection the same regardless of the cost.
AC
422
ED
M
421
AN US
417
The number of infection averted is given as the difference between the total infectious individuals without control and the total infectious individuals with control. The strategy with the highest ratio is the most cost-effective.
415
25
ACCEPTED MANUSCRIPT
5
A3
Infection Averted
A4
10
5
0
0.001
0.01
5.0
A4
4 3 2 1 0
10.0
A3
0.001
Weights
10.0
0.6
0.2
0.01
5.0
Weights
M
0.4
A3 A4
10
Total Cost
0.8
×105
AN US
A4
12
10.0
ED
Average Cost-Effectiveness Ratio (ACER)
A3
(c)
5.0
(b)
1
0.001
0.01
Weights
(a)
0
CR IP T
×105
Infection Averted Ratio (IAR)
15
8 6 4 2 0
0.001
0.01
5.0
10.0
Weights
(d)
PT
Figure 10: Simulation results of the control profile of model (3) using cost weights A3 = A4 =
433 434 435 436 437 438 439 440
441 442 443
Control within interval of uncertainty for bM and µM In carrying out the sensitivity analysis test we set the values of the parameters to be within ± 20% range of their baseline values and found the parameters bM , and µM among ˜ 0 and R0 . However, the time others to have a high impact on the reproduction numbers R dependent controls uV (t) and uM (t) incorporated into system (3) did not directly act on these parameters. Hence, we vary first, the parameter bM using a range of ± 50% from the baseline value and determine its impact on the system and the control variables. We have used the ± 50% range for emphasis purpose, the result is expected to be the same in the ± 20% range.
AC
432
CE
0.0001, A3 = A4 = 0.001, A3 = A4 = 10.0, and A3 = A4 = 13 (a) Infection averted; (b) Infection averted ratio (IAR); (c) Average cost-effective ratio (ACER); (d) Total cost.
The results of the optimal control simulations of the dengue model (3) in the uncertainty interval for the parameter bM are depicted in Figures 11(a) - 11(b). The parameter bM is the mosquitoes biting rate, since Aedes egypyti bite during the day, knowledge of the 26
ACCEPTED MANUSCRIPT
impact of the control in the uncertainty range of this parameter is vital. 1400
Lower bound Baseline values Upper bound
1500
1000
500
0
0
20
40
60
80
Lower bound Baseline values Upper bound
1200 1000 800 600 400 200 0
100
0
20
Time (days)
CR IP T
Infectious individuals
2000
Infectious mosquitoes
444
40
60
80
100
Time (days) (b)
AN US
(a)
Figure 11: Simulation of the dengue model (3) as a function of time without control (solid lines) and with optimal control (dash lines) for: (a). Total number of infectious humans; (b). Total number of infectious mosquitoes. The parameter bM is set within ± 50% interval of its baseline values given in Table 2 and A1 = 1.0, A2 = 1.0, A3 = 0.0001, A4 = 0.01.
448
1
PT
0.6
0.8
CE
control uV
0.8
0.4
AC
0.2 0
1
Lower bound Baseline values Upper bound
control uM
447
M
446
We observed that the upper bound led to the most number of infectious in humans (EH and IH ), and mosquitoes. However, due to the actions of the control variables depicted in Figures 12(a)-12(b), the number of the infectious humans was drastically reduced while the infectious mosquitoes remain the same within the period of simulation.
ED
445
0
20
40
0.6 0.4 Lower bound Baseline values Upper bound
0.2
60
80
0
100
Time (days)
0
20
40
60
80
100
Time (days)
(a)
(b)
Figure 12: The optimal controls of the dengue model (3) in the ± 50% uncertainty interval of parameter bM for: (a). control uV ; (b). control uM .
449 450
Hence, regardless of the parameter values chosen within the uncertainty interval of bM the control profile for the mosquitoes remained the same. Although, the duration of the 27
ACCEPTED MANUSCRIPT
455
456 457 458 459 460 461 462
Next, we investigate the optimal control of system (3) for different values of the mosquito life span; we set µM = 1/14, 1/21, 1/42 [43]. The results are depicted in Figures 13(a)13(b). We observed the same profile for infectious humans regardless of the mosquito life span either with or without control. The control led to fewer infectious humans. For the mosquitoes, there were more infected mosquitoes as expected when the mosquitoes live longer (µM = 1/42) and fewer infected ones when they live for a shorter period of time (µM = 1/14).
CR IP T
453 454
controls at the maximum level are longer when parameter bM is at the upper bound and shorter at the lower bound. The duration of the control uV is marginally different at the upper bound. Note, that we have used very low weights, the profile is the same even for higher weights and tighter intervals as the ± 20% uncertainty interval used for the sensitivity analysis.
µ = 1/42 M
µ = 1/21 M
µM = 1/14
1500
1000
0
M
500
1400
0
20
40
60
80
100
ED
Infectious individuals
2000
AN US
452
µ = 1/42 M
1200
µ = 1/21 M
µ = 1/14
1000
M
800 600 400 200 0
0
20
40
60
80
100
Time (days)
Time (days) (a)
Infectious mosquitoes
451
(b)
PT
Figure 13: The optimal controls of the dengue model (3) when µM = 1/14, 1/21, 1/42 for : (a).
464 465 466
For the optimal controls uV (t) and uM (t), there are marginal differences in the control profiles in Figures 14(a)-14(b). The implication of these results are that, regardless of the mosquito life span the control efforts are the same and they will adequately reduce the infectious populations both in humans and mosquitoes.
AC
463
CE
Total number of infectious humans; (b). Total number of infectious mosquitoes. The parameter values used are given in Table 2 and A1 = 1.0, A2 = 1.0, A3 = 0.0001, A4 = 0.01.
28
ACCEPTED MANUSCRIPT
1
1 µM = 1/42 µ = 1/21
0.8
0.6 0.4 0.2 0
µ = 1/42 M µ = 1/21
0.6
M
µ = 1/14 M
0.4 0.2
0
20
40
60
80
0
100
Time (days)
0
20
CR IP T
M
µM = 1/14
control uM
control uV
0.8
40
60
80
100
Time (days)
(a)
(b)
471 472 473 474
475 476 477 478 479 480 481
482 483 484 485 486 487
488 489
M
ED
470
Dengue is a viral infection transmitted via infected female Aedes mosquito bites resulting in either a mild dengue fever or haemorrhagic fever [65]. There are four distinct dengue virus serotypes (namely DEN-1, DEN-2, DEN-3 and DEN-4) [65]. Since 2010, three provinces in Pakistan namely Khyber Pakhtunkhwa, Punjab and Sindh has been experiencing dengue fever epidemic [66]. The first confirmed dengue fever outbreak in the country was in reported in 1994; however, the first annual epidemic trend and sudden rise in cases occurred in Karachi in November 2005 [66].
PT
469
Discussion and Conclusion
The government of Pakistan solicited technical assistance from the World health Organization (WHO) to combat the epidemic by emphasizing on vector control, case management and community awareness [66]. Unfortunately, these strategies were not effective, as at 2010, there was another sudden rise in cases with more reported cases than the preceding five years combined [66]. The country has since faced a number of repeated outbreaks in different provinces [69]. In 2017, the country experienced another outbreak starting from the province of Khyber Pakhtunkhwa spreading to Sindh province, and AJK province [69].
CE
468
7
AC
467
AN US
Figure 14: The optimal controls of the dengue model (3) when µM = 1/14, 1/21, 1/42 for: (a). Control uV ; (b). Control uM .
In this paper, we applied optimal control theory and investigate the impact of using an imperfect vaccine and mosquito reduction strategy in the bid to control dengue in Pakistan. For our investigation we use the dengue model (1) proposed in Section 2. The analysis shows that the disease-free equilibrium of the model is locally asymptotically stable whenever the associated reproduction number (R0 ), is less than unity and unstable otherwise. ˜0 To identify the parameters with the strongest impact on the reproduction numbers (R and R0 ), the model outcome, we used Latin Hypercube Sampling (LHS) and partial rank 29
ACCEPTED MANUSCRIPT
498 499 500 501 502 503 504 505
506 507 508 509 510 511 512 513 514
515 516 517 518 519 520 521 522 523 524 525 526 527 528
529 530
CR IP T
496 497
In particular, the results of this sensitivity analysis suggest that a strategy that reduces the transmission probability per contact of infectious mosquitoes with susceptible humans (αH and βH ), mosquito recruitment rate (ΛM ), mosquito biting rate (bM ), transmission probability per contact of susceptible mosquitoes with infectious humans (αM and βM ), disease progression rate in mosquitoes (σM ), and vaccine waning rate (ωH ). Furthermore, a strategy that increases the human recruitment rate (ΛH ), disease progression rate in humans (σH ), human recovery rate (γH ), mosquito natural death rate (µM ), and vaccine efficacy (ε) will be effective in curtailing the spread of dengue in the community.
AN US
495
Hence, we introduced into the model two time-dependent control variables representing vaccination and the use of insecticide. We found that the total number infectious (exposed and infected) individuals can be reduced in the community by the application of time dependent controls. Our results are in line with the results obtained in [34,35,39] which respectively showed that vaccine when available can be a promising tool to fight the disease in a community. Our study used a simple model and strategy to reach the same conclusion unlike Rodrigues et al. [35] which investigated using optimal control strategies the different ways of distributing vaccines: pediatric and random mass vaccines, with distinct levels of efficacy and durability.
M
494
ED
493
PT
492
correlation coefficients (PRCC) methods. The results of the sensitivity analysis of the model show that the parameters that have the most influence the reproduction number, R0 , are the human recruitment rate (ΛH ), disease progression rate in humans (σH ), human recovery rate (γH ), mosquito natural death rate (µM ), and vaccine efficacy (ε), the transmission probability per contact of infectious mosquitoes with susceptible humans (αH and βH ), mosquito recruitment rate (ΛM ), mosquito biting rate (bM ), transmission probability per contact of susceptible mosquitoes with infectious humans (αM and βM ), disease progression rate in mosquitoes (σM ), and vaccine waning rate (ωH ).
The application of optimal control theory requires an objective/cost functional. The optimal control objective or cost functional we utilized, minimized the numbers of exposed and infected individuals along with the cost of applying the controls. We choose a quadratic objective functionalin this study, since the costs of the intervention are nonlinear. This assumption is based on the fact that there are no linear relationship between the effects of intervention and the cost of intervention of the infective populations. Rodrigues et al. [34] used a cost functional which depends on the costs of medical treatment of the infected people, and the costs related to educational and sanitation campaigns. They then went on to show efficient ways to decrease the number of infected mosquitoes and individuals in less time and with lower costs using current computational tools. Although, that was not our approach in this study, we nevertheless investigated the associated benefits of the control strategies using cost-effectiveness analysis. This was achieved, using three approaches, the infection averted ratio (IAR), the average cost-effectiveness ratio (ACER) and the objective functional.
CE
491
AC
490
We observed while investigating the effect of the costs on the cost weights a reciprocal relation exist between the cost of insecticde use and vaccination; as the cost of insecticide 30
ACCEPTED MANUSCRIPT
538 539 540 541 542 543
544 545 546 547
548 549 550 551 552
553
554 555 556 557
558 559 560 561 562 563 564 565
566 567 568
569
CR IP T
537
AN US
536
Furthermore, our study showed that lower weights on the control are more cost-effective than higher weights. That is controls with lower weights are cheaper strategies to use compare to controls with higher weights. Insecticide led to lower cost compare to vaccination at higher weights; therefore, insecticide is more cost-effect that vaccination at lower weight. Shim [39] in their study of the dynamics of dengue and cost-effectiveness of vaccine in the Philippines using an age-structured model showed that the higher the cost of vaccination, the less cost-effective the dengue vaccination program. Furthermore, their results demonstrated that at a relatively low vaccine efficacy an age-targeted vaccination may be cost-effective provided the cost is sufficiently low. Although, we did not investigate the effect of vaccine efficacy on cost-effectiveness and overall performance; it is easy to see that due to the reciprocal relationship between insecticide use and vaccination the reduce performance in vaccination due to reduced efficacy will be off-set by increase in insecticide use. To conclude, in this study, we have presented a deterministic model of ordinary differential equations for dengue infection transmission dynamics. And we have studied using optimal control theory the use of insecticide, and vaccination as effective control measures against the epidemics using data obtained from hospital from the Peshawar District in Pakistan. The following results were observed from our analysis and numerical simulations:
M
535
ED
534
(i) The model has a DFE that is locally asymptotically stable if R0 < 1; (ii) Hence, we estimated in this study the basic reproduction number (R0 ) without any interventions for the 2017 dengue outbreak in Peshawar district of Pakistan as R0 ≈ 2.64, the distribution of the reproduction number lies in the range R0 ∈ [1.21, 5.24] (with a mean R0 ≈ 2.64);
PT
533
increases the use of vaccination increases. Only a slight increase in insecticide use is observed when vaccination level decreased due to increase in cost. The two controls averts about the same number of infections in the district regardless of the weights on the costs this is due to the reciprocal relation existing between the cost of insecticde use and vaccination.
CE
532
(iii) The sensitivity analysis of the model shows that the dominant parameters influencing ˜ 0 and R0 ) are the human recruitment rate (ΛH ), disthe reproduction numbers s (R ease progression rate in humans (σH ), human recovery rate (γH ), mosquito natural death rate (µM ), and vaccine efficacy (ε), the transmission probability per contact of infectious mosquitoes with susceptible humans (αH and βH ), mosquito recruitment rate (ΛM ), mosquito biting rate (bM ), transmission probability per contact of susceptible mosquitoes with infectious humans (αM and βM ), disease progression rate in mosquitoes (σM ), and vaccine waning rate (ωH ).
AC
531
(iv) The application of two time-dependent controls uV (t) and uM (t) obtained from the sensitivity analysis can reduce the total number of infected humans and mosquitoes in the district; (v) Investigating the effect of the costs on the weights A3 , and A4 we observed that 31
ACCEPTED MANUSCRIPT
571 572 573 574 575 576 577 578
579
(a) A reciprocal relation exist between the cost of insecticde use and vaccination; as the cost of insecticide increases the use of vaccination increases. Only a slight increase in insecticide use is observed when vaccination level decreased due to increase in cost. (b) The two controls averts about the same number of infections in the district regardless of the weights on the costs this is due to the reciprocal relation exist between the cost of insecticde use and vaccination.
CR IP T
570
(c) Lower weights are more cost-effective than higher weights with insecticide leading to cheaper cost than vaccination.
Acknowledgments
582
References
588
589 590 591
592 593 594
595 596
597 598
599 600
601 602 603
M
ED
587
[3] F.B. Agusto, S. Bewick, and W.F. Fagan. Mathematical model of zika virus with vertical transmission. Infectious Disease Modelling, 2(2):244–267, 2017. [4] F.B. Agusto, S.Y. Del Valle, K.W. Blayneh, C.N. Ngonghala, M.J. Goncalves, N. Li, R. Zhao, and H. Gong. The impact of bed-net use on malaria prevalence. Journal of theoretical biology, 320:58–65, 2013.
PT
586
[2] F.B. Agusto. Optimal isolation control strategies and cost-effectiveness analysis of a two-strain avian influenza model. BioSystems, 113(3):155–164, 2013.
[5] F.B. Agusto, S. Easley, K. Freeman, and M. Thomas. Mathematical model of three age-structured transmission dynamics of chikungunya virus. Computational and mathematical methods in medicine, 2016, 2016.
CE
585
[1] F. Agusto and S. Lenhart. Optimal control of the spread of malaria superinfectivity. Journal of Biological Systems, 21(04):1340002, 2013.
AC
583 584
AN US
581
One of the authors, FBA, was supported by the Strategic Environmental Research and Development Program under grant RC-2639.
580
[6] F.B. Agusto and I.M. ELmojtaba. Optimal control and cost-effective analysis of malaria/visceral leishmaniasis co-infection. PloS one, 12(2):e0171102, 2017.
[7] F.B. Agusto and A.B. Gumel. Theoretical assessment of avian influenza vaccine. DCDS Series B, 13(1):1–25, 2010. [8] F.B. Agusto and A.B. Gumel. Qualitative dynamics of lowly-and highly-pathogenic avian influenza strains. Mathematical biosciences, 243(2):147–162, 2013.
[9] F.B. Agusto, N. Marcus, and K.O. Okosun. Application of optimal control to the epidemiology of malaria. Electronic Journal of Differential Equations, 2012(81):1–22, 2012. 32
ACCEPTED MANUSCRIPT
605
[10] T.L Bancroft. On the aetiology of dengue fever. Australian Medical Gazette, 25:17–18, 1906.
606
[11] G. B´eraud. Mathematical models and vaccination strategies. Vaccine, 2017.
611 612 613
614 615 616 617
618 619 620
621 622 623
624 625 626 627
CR IP T
610
[13] J.E. Blaney, J.M. Matro, B.R. Murphy, and S.S. Whitehead. Recombinant, liveattenuated tetravalent dengue virus vaccine formulations induce a balanced, broad, and protective neutralizing antibody response against each of the four serotypes in rhesus monkeys. Journal of virology, 79(9):5516–5528, 2005. [14] J.E. Blaney, N.S. Sathe, C.T. Hanson, C.Y. Firestone, B.R. Murphy, and S.S. Whitehead. Vaccine candidates for dengue virus type 1 (den1) generated by replacement of the structural genes of rden4 and rden4δ30 with those of den1. Virology journal, 4(1):23, 2007.
AN US
609
[12] S. Bhatt, P.W. Gething, O.J. Brady, J.P. Messina, A.W. Farlow, C.L. Moyes, J.M. Drake, J.S. Brownstein, A.G. Hoen, O. Sankoh, et al. The global distribution and burden of dengue. Nature, 496(7446):504, 2013.
[15] K.W. Blayneh, A.B. Gumel, S. Lenhart, and T. Clayton. Backward bifurcation and optimal control in transmission dynamics of west nile virus. Bulletin of mathematical biology, 72(4):1006–1028, 2010. [16] S.M. Blower and H.I. Dowlatabadi. Sensitivity and uncertainty analysis of complex models of disease transmission: an hiv model, as an example. International Statistical Review/Revue Internationale de Statistique, pages 229–243, 1994.
M
608
ED
607
[17] O.J. Brady, P.W. Gething, S. Bhatt, J.P. Messina, J.S. Brownstein, A.G. Hoen, C.L. Moyes, A.W. Farlow, T.W. Scott, and S.I. Hay. Refining the global spatial limits of dengue virus transmission by evidence-based consensus. PLoS neglected tropical diseases, 6(8):e1760, 2012.
PT
604
630
[19] J Carr. Applications of Centre Manifold Theory. Springer-Verlag, New York, 1981.
632 633
634 635
636 637
638 639 640
[20] C. Castillo-Chavez, S. Blower, P. van den Driessche, D. Kirschner, and A.A. Yakubu. Mathematical approaches for emerging and reemerging infectious diseases: an introduction, volume 1. Springer Science & Business Media, 2002.
AC
631
CE
629
[18] F. Brauer. Backward bifurcations in simple vaccination models. Journal of Mathematical Analysis and Applications, 298(2):418–431, 2004.
628
[21] C. Castillo-Chavez and B. Song. Dynamical models of tuberculosis and their applications. Mathematical biosciences and engineering, 1(2):361–404, 2004.
[22] Centres for Disease Control and Prevention. Dengue. https: // www. cdc. gov/ dengue/ index. html , 2016. [23] G. Chowell, P. Diaz-Duenas, J.C. Miller, A. Alcazar-Velazco, J.M. Hyman, P.W. Fenimore, and C. Castillo-Chavez. Estimation of the reproduction number of dengue fever from spatial epidemic data. Mathematical biosciences, 208(2):571–589, 2007. 33
ACCEPTED MANUSCRIPT
649
650 651 652
653 654
655 656
657 658
659 660
661 662
663 664
665 666
667 668 669
670 671 672
673 674 675 676
CR IP T
648
[27] J. Dushoff, W. Huang, and C. Castillo-Chavez. Backwards bifurcations and catastrophe in simple models of fatal diseases. Journal of mathematical biology, 36(3):227–248, 1998. [28] E.H. Elbasha and A.B. Gumel. Theoretical assessment of public health impact of imperfect prophylactic hiv-1 vaccines with therapeutic benefits. Bulletin of mathematical biology, 68(3):577, 2006.
AN US
647
[29] C. Esteva, L.and Vargas. Influence of vertical and mechanical transmission on the dynamics of dengue disease. Mathematical biosciences, 167(1):51–64, 2000. [30] L. Esteva and C. Vargas. Analysis of a dengue disease transmission model. Mathematical biosciences, 150(2):131–151, 1998. [31] L. Esteva and C. Vargas. A model for dengue disease with variable human population. Journal of mathematical biology, 38(3):220–240, 1999.
M
646
[26] I. Dorigatti, C. McCormack, G. Nedjati-Gilani, and N.M. Ferguson. Using wolbachia for dengue control: Insights from modelling. Trends in parasitology, 2017.
[32] F.B. Agusto. Optimal control and cost-effective analysis of the 2017 meningitis outbreak in nigeria. Submitted manuscript, 2018.
ED
645
[33] W.H. Fleming and R.W. Rishel. Deterministic and stochastic optimal control, volume 1. Springer Science & Business Media, 2012.
PT
644
[25] M. Derouich and A. Boutayeb. Dengue fever: Mathematical modelling and computer simulation. Applied Mathematics and Computation, 177(2):528–544, 2006.
[34] S.M. Garba, A.B. Gumel, and M.R.A. Bakar. Backward bifurcations in dengue transmission dynamics. Mathematical biosciences, 215(1):11–25, 2008.
CE
643
[24] Dengue Virus Ne. Treatment of dengue. http: // www. denguevirusnet. com/ treatment. html , 2018.
[35] D.J. Gubler. Dengue and dengue hemorrhagic fever. Clinical microbiology reviews, 11(3):480–496, 1998. [36] D.J. Gubler and G. Kuno. Dengue and dengue hemorrhagic fever: its history and resurgence as a global public health problem. Dengue and dengue hemorrhagic fever. CAB international, London, United Kingdom, pages 1–22, 1997.
AC
641 642
[37] S.B. Halstead, S. Nimmannitya, and S.N. Cohen. Observations related to pathogenesis of dengue hemorrhagic fever. iv. relation of disease severity to antibody response and virus recovered. The Yale journal of biology and medicine, 42(5):311, 1970. [38] I. Kawaguchi, A. Sasaki, and M. Boots. Why are dengue virus serotypes so distantly related? enhancement and limiting serotype similarity between dengue virus strains. Proceedings of the Royal Society of London B: Biological Sciences, 270(1530):2241– 2247, 2003. 34
ACCEPTED MANUSCRIPT
683
684 685 686
687 688 689
690 691 692
693 694 695
696 697 698
[41] S. Lenhart and J.T. Workman. Optimal control applied to biological models. CRC Press, 2007.
CR IP T
682
[42] S.B. Maier, X. Huang, E. Massad, M. Amaku, M.N. Burattini, and D. Greenhalgh. Analysis of the optimal vaccination age for dengue in brazil with a tetravalent dengue vaccine. Mathematical biosciences, 294:15–32, 2017. [43] C.A. Manore, K.S. Hickmann, S. Xu, H.J. Wearing, and J.M. Hyman. Comparing dengue and chikungunya emergence and endemic transmission in a. aegypti and a. albopictus. Journal of theoretical biology, 356:174–191, 2014.
AN US
681
[44] S. Marino, I.B. Hogue, C.J. Ray, and D. E. Kirschner. A methodology for performing global uncertainty and sensitivity analysis in systems biology. Journal of theoretical biology, 254(1):178–196, 2008. [45] M.D. McKay, R.J. Beckman, and W.J. Conover. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics, 42(1):55–61, 2000.
M
680
[40] M.R. Khanani, A. Arif, and R. Shaikh. Dengue in pakistan: Journey from a disease free to a hyper endemic nation. Journal of the Dow University of Health Sciences Karachi, 5(3):81–84, 2011.
[46] R.G. McLeod, J.F. Brewster, A.B. Gumel, and A. Slonowsky. Sensitivity and uncertainty analyses for a sars model with time-varying inputs and outputs. Mathematical Biosciences and Engineering, 3(3):527, 2006.
ED
679
[39] J. Khan, A. Ghaffar, and S.A. Khan. The changing epidemiological pattern of dengue in swat, khyber pakhtunkhwa. PloS one, 13(4):e0195706, 2018.
PT
677 678
703
[48] L.S. Pontryagin. Mathematical theory of optimal processes. CRC Press, 1987.
701
704 705 706
707 708
709 710 711
AC
700
CE
702
[47] Pakistan Bureau of Statistics. Pakistan’s 6th census: Population of major cities census. http: // www. pbscensus. gov. pk/ sites/ default/ files/ population_ of_ major_ cities_ census_ 2017% 20_ 0. pdf , 2017.
699
[49] H.S. Rodrigues, M.T.T. Monteiro, and D.F.M. Torres. Dynamics of dengue epidemics when using optimal control. Mathematical and Computer Modelling, 52(9-10):1667– 1673, 2010.
[50] H.S. Rodrigues, M.T.T. Monteiro, and D.F.M. Torres. Vaccination models and optimal control strategies to dengue. Mathematical biosciences, 247:1–12, 2014. [51] M.A. Sanchez and S.M. Blower. Uncertainty and sensitivity analysis of the basic reproductive rate: tuberculosis as an example. American Journal of Epidemiology, 145(12):1127–1137, 1997.
35
ACCEPTED MANUSCRIPT
714
715 716 717
718 719 720
[52] T. Sardar, S. Rana, and J. Chattopadhyay. A mathematical model of dengue transmission with memory. Communications in Nonlinear Science and Numerical Simulation, 22(1-3):511–525, 2015. [53] O. Sharomi, C. Podder, A. Gumel, and B. Song. Mathematical analysis of the transmission dynamics of hiv/tb coinfection in the presence of treatment. Mathematical Biosciences and Engineering, 5(1):145, 2008.
CR IP T
712 713
[54] O. Sharomi, C.N. Podder, A.B. Gumel, E.H. Elbasha, and J. Watmough. Role of incidence function in vaccine-induced backward bifurcation in some hiv models. Mathematical Biosciences, 210(2):436–463, 2007. [55] C. Shekhar. Deadly dengue: new vaccines promise to tackle this escalating global menace. Chemistry & biology, 14(8):871–872, 2007.
723
[56] S.M. Shepherd. Dengue. Medscape Logo, 2017.
731
732 733
734 735 736
737 738 739
740 741
742 743 744 745
746 747
M
730
[59] M.I. Teboh-Ewungkem, C.N. Podder, and A.B. Gumel. Mathematical study of the role of gametocytes and an imperfect vaccine on malaria transmission dynamics. Bulletin of mathematical biology, 72(1):63–93, 2010.
ED
729
[60] J.J. Tewa, J.L. Dimi, and S. Bowong. Lyapunov functions for a dengue disease transmission model. Chaos, Solitons & Fractals, 39(2):936–941, 2009.
PT
727 728
[58] C.J. Struchiner, P.M. Luz, C.T. Codeco, F.C. Coelho, and E. Massad. Current research issues in mosquito-borne diseases modeling. Contemporary Mathematics, 410:349–366, 2006.
[61] P. Van den Driessche and J. Watmough. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical biosciences, 180(1):29–48, 2002.
CE
726
[57] E. Shim. Dengue dynamics and vaccine cost-effectiveness analysis in the philippines. The American journal of tropical medicine and hygiene, 95(5):1137–1147, 2016.
[62] WHO Health Organization (WHO). Who country cooperation strategic. http: // apps. who. int/ iris/ bitstream/ 10665/ 136607/ 1/ ccsbrief_ pak_ en. pdf , 2016.
AC
724 725
AN US
722
721
[63] World Health Organization (WHO). Dengue and severe dengue. http: // www. who. int/ mediacentre/ factsheets/ fs117/ en/ , 2017.
[64] World Health Organization (WHO). Dengue vaccine research. Immunization, Vaccines and Biologicals http: // www. who. int/ immunization/ research/ development/ dengue_ vaccines/ en/ , 2017. [65] World Health Organization (WHO). Epidemiology. http: // www. who. int/ denguecontrol/ epidemiology/ en/ , 2018. 36
ACCEPTED MANUSCRIPT
757
758 759 760 761
762 763
CR IP T
756
[68] World Health Organization (WHO), Eastern Mediterranean Regional Office (EMRO). Dengue fever in pakistan. Weekly Epidemiological Monitor https: // reliefweb. int/ sites/ reliefweb. int/ files/ resources/ Epi_ Monitor_ 2017_ 10_ 36. pdf , 10(36), 03 September 2017. Accessed March 16, 2018. [69] World Health Organization (WHO), Eastern Mediterranean Regional Office (EMRO). Deaths from dengue in pakistan. Weekly Epidemiological Monitor. https: // reliefweb. int/ sites/ reliefweb. int/ files/ resources/ Epi_ Monitor_ 2017_ 10_ 46. pdf , 10(46), 12 November 2017. Accessed March 16, 2018.
AN US
754 755
[70] H.M. Yang and C.P. Ferreira. Assessing the effects of vector control on dengue transmission. Applied Mathematics and Computation, 198(1):401–413, 2008.
M
753
[67] World Health Organization (WHO). Dengue and severe dengue. http: // www. who. int/ news-room/ fact-sheets/ detail/ dengue-and-severe-dengue , 2018. Accessed September 3, 2018.
ED
752
PT
751
CE
750
[66] World Health Organization (WHO). Pakistan. http: // www. emro. who. int/ pak/ programmes/ dengue-fever. html , 2018. Accessed July 24, 2018.
AC
748 749
37
ACCEPTED MANUSCRIPT
765 766 767
A
Appendix A: Proof of Lemma 1
Lemma 1: Let the initial data F (0) ≥ 0 , where F (t) = (SH (t), VH (t), EH (t), IH (t), RH (t), SM (t), EM (t), IM (t)). Then the solutions F (t) of the model (1) are non-negative for all t > 0. Furthermore lim sup NH (t) ≤ t→∞
ΛH ΛM , and lim sup NM (t) ≤ . µH µM t→∞
CR IP T
764
where
NH (t) = SH (t) + VH (t) + EH (t) + IH (t) + RH (t), and NM (t) = SM (t) + EM (t) + IM (t). 768
772
773
AN US
M
771
dSH bM (αH EM + βH IM )SH = ΛH + ωH VH − − (νH + µH )SH . (A-1) dt NH bM (αH EM + βH IM ) Let λH (t) = , then (A-1) becomes NH dSH = ΛH + ωH VH − λH (t)SH − (νH + µH )SH dt which can be re-written as ( Z t1 ) Z t1 d SH (t) exp λH (ζ)dζ + (νH + µH )t = (ΛH + ωH VH ) exp λH (ζ)dζ(νH + µH )t . dt 0 0
ED
770
Proof. Let t1 = sup{t > 0 : F (t) > 0 ∈ [0, t]}. Thus, t1 > 0. It follows from the first equation of the system (1), that
Hence, t1
λH (ζ)dζ + (νH + µH )t1
CE
SH (t1 ) exp
Z
PT
769
0
so that,
"
SH (t1 ) = SH (0) exp − "
+ exp − ×
exp
> 0.
− SH (0) = ×
AC
774
Z
" Z
Z
Z
t1
(ΛH + ωH VH (p)) Z p exp λH (ζ)dζ + (νH + µH )p dp 0
0
t1
λH (ζ)dζ + (νH + µH )t1
0
t1
λH (ζ)dζ + (νH + µH )t1
0
p
λH (ζ)dζ + (νH + µH )p
0
38
#
# Z
0
dp
#
t1
(ΛH + ωH VH (p))
ACCEPTED MANUSCRIPT
776 777 778
779
Similarly, it can be shown that F > 0 for all t > 0. For the second part of the proof, note that 0 < SH (0) ≤ NH (t), 0 ≤ VH (0) ≤ NH (t), 0 ≤ EH (0) ≤ NH (t), 0 ≤ (IH 0) ≤ NH (t), 0 ≤ RH (0) ≤ NH (t), 0 < SM (0) ≤ NM (t), 0 ≤ EM (0) ≤ NM (t), 0 ≤ IM (0) ≤ NM (t). Adding the human and mosquito components of model (1) gives
CR IP T
775
dNH (t) = ΛH − µH NH − δH IH dt ≤ ΛH − µH NH and
AN US
780
(A-2)
dNM (t) = ΛM − µM NM dt Hence, lim sup NH (t) ≤ t→∞
784
Lemma 2: The region Ω = ΩH × ΩM ⊂ R5+ × R3+ is positively-invariant for the model (1) with non-negative initial conditions in R8+ . Proof. It follows from summing up the human and mosquito equations of model (1) that
786
787
788
789
AC
CE
785
Appendix B: Proof of Lemma 2
ED
783
B
PT
782
M
as required.
781
ΛH ΛM , and lim sup NM (t) ≤ µH µM t→∞
dNH (t) = ΛH − µH NH − δH IH dt ≤ ΛH − µH NH
and
dNM (t) = ΛM − µM NM dt
dNH (t) ΛH dNM (t) ΛM ≤ 0, if NH (0) ≥ and ≤ 0, if NM (0) ≥ . dt µH dt µM ΛH ΛM Thus, NH (t) ≤ NH (0)e−µH t + (1 − e−µH t ) and NM (t) ≤ NM (0)e−µM t + (1 − e−µM t ). µH µM ΛH ΛM In particular, NH (t) ≤ and NM (t) ≤ . µH µM
Hence,
39
ACCEPTED MANUSCRIPT
794 795 796 797 798 799
CR IP T
793
C
Appendix C: Proof of Theorem 1
The phenomenon of backward bifurcation is shown using the concept of centre manifold theory on system (1). First, we make the following changes in variables by setting SH = x1 , VH = x2 , EH = x3 , IH = x4 , RH = x5 , SM = x6 , EM = x7 and − IM = x8 Furthermore, employing the vector notation → x = (x1 , x2 , x3 , x4 , x5 , x6 , x7 , x8 ), − the dengue model (1) with vaccination can be reformulated in the form dx/dt = F → x , with F = (f1 , f2 , f3 , f4 , f5 , f6 , f7 , f8 )T as given below = ΛH + ωH x2 − λH x1 − (νH + µH )x1 ,
(C-1)
= νH x1 − (1 − ε)λH x2 − (ωH + µH )x2 , = λH [x1 + (1 − ε)x2 ] − (σH + µH )x3 , = σH x3 − (γH + µH + δH )x4 , = γH x4 − µH x5 , = ΛM − λM x6 − µM x6 ,
where,
= λM x6 − (σM + µM )x7 , = σ M x7 − µ M x8 ,
AC
800
CE
PT
dx1 dt dx2 dt dx3 dt dx4 dt dx5 dt dx6 dt dx7 dt dx8 dt
AN US
792
M
791
ΛM ΛH and NM (0) > , µH µM ΛH then either the solutions enters Ω in finite time, or NH (t) approaches and NM (t) µH ΛM asymptotically. Hence, the region Ω attracts all solutions in R8+ . approaches µM Thus, the region Ω is positively-invariant. Furthermore, if NH (0) >
ED
790
λH =
801
bM (x7 αH + βH x8 ) bM (x3 αM + x4 βM ) , λM = , NH NH
and NH = x1 + x2 + x3 + x4 + x5 .
The Jacobian of system (C-1) evaluated at the disease-free equilibrium, E0 , is given by
40
ACCEPTED MANUSCRIPT
νH 0 0 J1 = 0 0 0 0
807
∗ −bM βH 0 β∗ (1−ε)bM VH H 0 SH 0 +(1−ε)V 0 β ∗ bM (SH H) H 0 SH
0
0
0
0
0
−k3
0
0
0
0 +(1−ε)V 0 α bM (SH H) H 0 SH
0
σH
−k4
0
0
0
0
γH
−µH
0
0
−µM
0
0
0
0
0
−
0 α bM SM M 0 SH
−
0 β bM SM M 0 SH
0
0 α bM SM M 0 SH
0 β bM SM M 0 SH
0
0
0
0
0
0
0
0
0
−k5
0
σM
−µM
Consider the case R0 = 1. Suppose that βH is chosen as a bifurcation parameter. Setting R0 = 1 and solving for βH using the expression for R0 gives o n 0 0 02 + (1 − ε)VH0 ] (σH βM + k4 αM ) [SH − b2M αH SM µM k3 k4 k5 SH 0 0 b2M σM SM (SH + (1 − ε)VH0 ) (σH βM + k4 αM )
.
PT
∗ The Jacobian matrix J1 evaluated at βH = βH has a simple zero eigenvalue with the other eigenvalues having negative real parts. Hence the centre manifold theorem can be applied ∗ . to analyze the dynamics of (C-1) near βH = βH
− Computing the right and the left eigenvector of J1 , denoted respectively by → w = [w1 , w2 , w3 , w4 , → − T w5 , w6 , w7 , w8 ] and v = [v1 , v2 , v3 , v4 , v5 , v6 , v7 , v8 ], we have 0 k3 k4 k5 w7 SH (k2 SH + (1 − ε)VH0 ωH ) , 0 0 bM S M (k1 k2 − νH ωH ) (SH + (1 − ε)VH0 ) (σH βM + k4 αM ) 0 0 k3 k4 k5 w7 SH (k1 (1 − ε)VH0 + νH SH ) , = 0 0 0 bM SM (k1 k2 − νH ωH ) (SH + (1 − ε)VH ) (σH βM + k4 αM ) 0 0 0 k4 k5 w7 SH k5 w7 σH SH k5 w7 γH σH SH , w = , w = , = 4 5 0 0 0 bM SM (σH βM + k4 αM ) bM SM (σH βM + k4 αM ) bM µ H S M (σH βM + k4 αM ) k5 w 7 w7 σM = − , w7 = w7 > 0, w8 = . µM µM
w1 = −
w2 w3 w6 809
−bM αH
−k2
AC
808
0
CE
805 806
0
0α (1−ε)bM VH H − 0 SH
∗ βH =
804
0
M
803
0
ED
802
ωH
CR IP T
−k1
AN US
and 41
ACCEPTED MANUSCRIPT
0 k3 k4 v3 SH k3 v3 βM , v7 = , v1 = v2 = v5 = v6 = 0, v3 = v3 > 0, v4 = 0 σH βM + k4 αM bM S M (σH βM + k4 αM ) 0 0 bM αH (SH − εVH + VH ) k3 k4 k5 SH . v8 = −v3 − 0 0 SH σM bM σM SM (σH βM + k4 αM )
811
Using the result from Theorem 4.1 by Castillo-Chavez and Song [21], we calculate the coefficients a and b as
CR IP T
810
n 2bM 0 0 v7 (w3 αM + w4 βM ) (w1 + w2 + w3 + w4 + w5 ) SM − w6 SH + VH0 0 0 2 (SH + VH ) o 0 (w2 ε + w3 + w4 + w5 ) . −v3 (w7 αH + w8 βH ) VH0 [w1 ε + (w3 + w4 + w5 ) (ε − 1)] − SH
The coefficient b is given as follows:
812
b=
AN US
a = −
v3 w7 bM σM ((1 − ε)νH + k2 ) . µM (νH + k2 )
815
D
M
814
The coefficient b is clearly positive; the presence of backward bifurcation in model (1) is determined by the sign of coefficient a.
813
Appendix D: Stability of model (2) DFE
818 819
The stability of E0 can be established using the next generation operator method on system (2). Taking EH (t), IH (t), EM (t), and IM (t) as the infected compartments and then using the notation in [61], the Jacobian F and V matrices for new infectious terms and the remaining transfer terms, respectively, are defined as: 0 0 bM α H S H bM βH SH 0 0 NH0 NH0 0 0 0 0 , F = 0 0 bM α M S M bM βM SM 0 0 0 S0 S H H 0 0 0 0 (µH + σH ) 0 0 0 −σH (γH + µH + δH ) 0 0 . V = 0 0 (µM + σM ) 0 0 0 −σM µM
CE
817
AC
816
PT
ED
The dengue model (2) has a disease free equilibrium (DFE). The DFE is obtained by setting the right-hand sides of the equations in the model (2) to zero, which is given by Λ ΛM H 0 0 , 0, 0, 0, , 0, 0 . E0 = (SH , 0, 0, 0, SM , 0, 0) = µH µM
42
ACCEPTED MANUSCRIPT
822 823 824
825 826
827
where, g1 = σH + µH , g2 = γH + µH + δH , g3 = σM + µM . ˜ 0 is the number of secondary infections in completely susceptible popuThe expression R lation due to infections from one introduced infectious individual with dengue. Further, using Theorem 2 in [61], the following result is established.
CR IP T
821
Therefore, the reproduction number is s 2 0 ˜ 0 = bM (αH µM + βH σM )(αM g2 + βM σH )SM , R 0 g1 g2 g3 µM SH
Lemma 4. The disease-free equilibrium (DFE) of the dengue model (2) is locally asymp˜ 0 > 1. totically stable (LAS) if R˜0 < 1 and unstable if R
E
828
Appendix E: Backward bifurcation analysis of the pre-intervention model (2)
AN US
820
Let,
∗ ∗ ∗ ∗ E1 = (SH , EH , IH , RH , SV∗ , EV∗ , IV∗ )
denote the endemic equilibrium of the dengue pre-intervention model (2) where λ∗H ΛH λ∗H ΛH σH λ∗H γH ΛH σH ΛH ∗ ∗ ∗ , E = , I = , R = , H H H µH + λ∗1 g1 (µH + λ∗1 ) g1 g2 (µH + λ1 ) g1 g2 µH (µH + λH ) ΛM λ∗M ΛM λ∗M ΛM σM ∗ ∗ = , E = , I = ; (E-1) M λ∗M + µM g3 (λ∗M + µM ) M g3 µM (λ∗2 + µM )
M
829
∗ SM
with λ∗1 =
833
Substituting the expression in (E-1) into (E-2), we obtain after some simplifications the quadratic polynomial
CE
832
(E-2)
z1 λ∗1 2 + z2 λ∗1 + z3 = 0,
(E-3)
AC
831
∗ ∗ ∗ ∗ ) βM ) αH + βH IM αM + IH bM (EM bM (EH ∗ , λ = . 2 ∗ ∗ NH NH
PT
830
ED
∗ SH =
where
z1 = g3 ΛH µM [g2 µH + σH (γH + µH )] [g2 µH (bM αM + µM ) + σH (bM µH βM + µM (γH + µH ))] , z2 = g1 g2 bM µ2H (g2 αM + σH βM ) [g3 ΛH µM − bM ΛM (αH µM + βH σM )] +2g1 g2 g3 ΛH µH µ2M [g2 µH + σH (γH + µH )] ,
834 835
e 2 ). z3 = g12 g22 g3 ΛH µ2H µ2M (1 − R 0
(E-4)
e 0 < 1 and negative when R e 0 > 1. Here, z1 is obviously positive and z3 positive whenever R Thus, the following result is established: 43
ACCEPTED MANUSCRIPT
836
837
838
Theorem 3. The dengue pre-intervention model (2) has: e 0 > 1; (i) a unique endemic equilibrium exists if and only if z3 < 0 and R
(ii) a unique endemic equilibrium exists if z2 < 0, and z3 = 0 or z22 − 4z1 z3 = 0; (iii) two endemic equilibria exists if z3 > 0, z2 < 0 and z22 − 4z1 z3 > 0;
840
(iv) otherwise no endemic equilibrium exists.
842 843 844 845
In Theorem 3, case (i) obviously shows that the model (2) has a unique endemic equilibrium e 0 > 1. Furthermore, case (iii) of the aforementioned Theorem 3, shows the whenever R e 0 < 1. To show this possibility, we set the possibility of backward bifurcation whenever R 2 e 0 , denoted by R e c and is polynomial z2 − 4z1 z3 = 0, and solve for the critical values of R given as ec = R
M
0.009 0.008
ED
0.01
PT
849
force of infection λ *H
848
z22 . 4z1 g12 g22 g3 ΛH µ2H µ2M
0.007 0.006
CE
847
1−
(E-5)
e 0 such that R ec < R e 0 < 1. Using So, backward bifurcation can occur for the values of R some of the parameter values given in Table 2 except for δH = 0.4097, µH = 0.0288, ΛH = 39, µM = 0.15 and ΛM = 30, we illustrate the possibility of backward bifurcation in Figure 15 for the dengue pre-intervention model (2).
AC
846
s
AN US
841
CR IP T
839
0.005 0.004 0.003 0.002 0.001 0
0.6
0.8
Rc
1
1.2
1.4
1.6
R0
Figure 15: Backward bifurcation curve for the dengue model (2).
44
ACCEPTED MANUSCRIPT
F
851
852 853 854 855
Global Asymptotic Stability of the DFE of the preintervention model (2) when δH = 0
Next, we consider the case when the disease-induced mortality rate is ignored. With proper medical care dengue fatality rates is below 1% [67]. In this case, we can either ignore the disease-induced death rate(δH = 0). Thus, we consider the case when δH = 0 in order to study the global stability of the disease-free equilibrium.
CR IP T
850
For dengue elimination to be independent of the initial sizes of the sub-populations of the model, the global asymptotic stability of the DFE must be established. This is considered next. Consider the feasible region ∗ ∗ Γ1 = {X ∈ Γ : SH ≤ SH , SM ≤ SM },
where , X = SH , EH , IH , RH , SM , EM , IM .
857
Lemma 5. The region Γ1 is positively invariant for model (2).
858
∗ Proof. It follows from the first equation of dengue model (2), (where SH =
AN US
856
ΛH ), so that µH
M
bM (αH EM (t) + βH IM (t))SH (t) dSH (t) = ΛH − − µH SH SB (t) dt NH (t) ≤ ΛH − µH SH (t)
ΛH − SH (t) µH
(F-1)
ED
≤ µH
PT
Thus,
CE
859
∗ = µH [SH − SH (t)].
∗ ∗ SH (t) ≤ SH − [SH − SH (0)]e−µH t .
∗ ∗ Thus, if SH (0) ≤ SH for all t ≥ 0, then SH (t) ≤ SH for all t ≥ 0.
861
∗ Similarly, it follows from the fifth equation of the dengue model (2) (where SM =
862
AC
860
that
ΛM ), µM
dSM (t) bM (αM EH (t) + βM IH (t))SM (t) = ΛM − − µM SM (t), dt NH (t) ≤ ΛM − µM SM (t), ≤ µM
ΛM − SM (t) , µM
∗ = µM [SM − SM (t)].
45
(F-2)
ACCEPTED MANUSCRIPT
867 868 869 870 871
872 873
874 875
876 877
∗ ∗ and SM (0) ≤ SM for all t ≥ 0, then SM (t) ≤ SM for all t ≥ 0.
Thus, in summary, it has been shown that the region Γ1 is positively invariant and attracts all solutions in R7+ for the dengue model (2).
CR IP T
866
ΛM µM
One of the sufficient conditions for the global asymptotic stability (GAS) of the diseasefree equilibrium (DFE) is a constant population level. This is typically the case when the disease-induced death rate is ignored (δH = 0). In this case the GAS result could be established using the second generation approach given in Castillo-Chavez et al. (2002) [20].
AN US
865
∗ Hence, if SM =
∗ ∗ SM (t) ≤ SM − [SM − SM (0)]e−µM t .
Theorem 4. The DFE, E0 , of model (2), is globally asymptotically stable (GAS) in Γ1 ˜ 0 ≤ 1. whenever R The above result shows that dengue will be eliminated from the community if the threshold ˜ 0 can be brought to a value less than unity. quantity R
M
864
Thus,
Theorem 4: The DFE, E0 , of model (2), is globally asymptotically stable (GAS) in Γ1 ˜ 0 ≤ 1 and δH = 0. whenever R
ED
863
880
Let X = (SH , RH , SM ) and Z = (EH , IH , EM , IM ) and group the system (2) into
CE
PT
879
Proof. To prove the global stability of the disease-free equilibrium, we will following the approach in [4, 3].
878
882
883
= F (X, 0)
dZ dt
= G(X, Z),
(F-3)
where F (X, 0) is the right hand side of S˙H , R˙H , and S˙M with EH = IB = IH = EM = ˙ . IM = 0 and G(X, Z) is the right hand side of E˙H , I˙H , E˙M and IM
AC
881
dX dt
Next, consider the reduced system: dSH dt dRH dt dSM dt
dX = F (X, 0) given as: dt
= ΛH − µH SH (t) = −µH RH (t)
(F-4)
= ΛM − µM SM (t). 46
ACCEPTED MANUSCRIPT
Let ∗
X =
∗ ∗ ∗ (SH , RH , SM )
=
ΛH ΛM , 0, . µH µM
886
To do this, solve the first and second equations of (F-4), this gives
SB (t) =
Z
t µH u
e
0
RH (t) =
Z
t µH u
e
887
Integrating SH (t) in (F-5) gives:
888
du + RH (0) e−µH t .
ΛH [1 − eµH t ] + SH (0)e−µH t . µH
ΛH and lim RH (t) = 0. t→∞ µH
ED
t=→∞
Similarly, solving for SM (t) in (F-4) gives SM (t) =
890
converges to
894
895
896
897 898
ΛM ΛM + e−µM t [SM (0) − ] which µM µM
CE
ΛM , as t → ∞. µM
These asymptotic dynamics are independent of initial conditions in Γ. Hence, the convergence of solutions of (F-4) is global in Γ1 . Next, following [20], we require G(X, Z) to satisfy the two stated conditions:
AC
893
PT
889
892
(F-6)
Taking the limit of SH (t) in (F-6) and RH (t) in (F-5) as t → ∞, we have lim SH (t) =
891
(F-5)
M
SH (t) =
ΛH du + SH (0) e−µH t
AN US
0
CR IP T
885
be an equilibrium of the reduced system (F-4), we show that X ∗ is a globally stable equilibrium in Γ1 .
884
(i). G(X, 0) = 0 and
b b (ii). G(X, Z) = DZ G(X ∗ , 0)Z − G(X, Z), G(X, Z) ≥ 0, ΛH ΛM ∗ where (X , 0) = , 0, , 0, 0, 0, 0 and DZ G(X ∗ , 0) is the Jacobian of G(X, Z) taken µH µM with respect to (EH , IH , EM , IM ) and evaluated at (X ∗ , 0), which is an M-matrix (the diagonal elements are nonnegative). Thus,
47
ACCEPTED MANUSCRIPT
−k1
0
∗ bM αH SH ∗ NH 0
σH −k2 DZ G(X ∗ , 0) = ∗ ∗ bM αM SM bM βM SM −k3 ∗ ∗ NH NH 0 0 σM ∗ SH (bM αH EM + bM βH IM ) 1− ∗ NH 0 b G(X, Z) = ∗ SM (b α E + b β M M H M M IH ) 1− ∗ NH 0
904
SH NH ∗ NH SH
∗ SM NH ∗ NH SM
AN US
Therefore, the disease-free equilibrium is globally asymptotically stable by the theorem in [20] (page 246).
M
903
ED
902
PT
901
−µM ∗
,
NH∗ SH Hence, if the human population is at equilibrium level, we have that 1 − ∗ > 0, SH NH ∗ SM NM b and 1 − > 0; thus, G(X, Z) ≥ 0. ∗ NM SM
CE
900
0
ΛH ΛH ΛM ∗ ∗ ∗ , NH∗ = and SM = . We have in Γ1 that, SH ≤ SH , and SM ≤ SM . µH µH µM Therefore, it follows that SH ≤ NH . ∗ where, SH =
AC
899
∗ bM βH SH ∗ NH 0
CR IP T
48