Multi-Strain Age-Structured Dengue Transmission Model: Analysis and Optimal Control

TL;DR

This paper develops a multi-strain, age-structured dengue transmission model to analyze the impact of vaccination and control strategies, revealing that traditional interventions outperform Dengvaxia in reducing cases.

Contribution

It introduces a novel multi-strain, age-structured mathematical model for dengue transmission and applies optimal control theory to evaluate intervention strategies.

Findings

Traditional control strategies outperform Dengvaxia in case reduction.

The model highlights the importance of age-specific interventions.

Optimal control simulations provide insights into effective dengue management.

Abstract

Dengue is a serious health problem in the Philippines. In 2016, the Department of Health of the country launched a dengue vaccination campaign using Dengvaxia. However, the campaign was mired with controversy and the use of Dengvaxia was banned in the country. This study proposes a mathematical model that represents the dynamics of the transmission of dengue with its four strains. Considering that the Dengvaxia vaccine was intended to be given only to people aging from 9 to 45 years old, the human population is divided into two age groups: from 9-45 years old and the rest of the population. Using this model and optimal control theory we simulate what could have been the effect of Dengvaxia in the number of dengue cases and then compare this with the other usual dengue intervention strategies. Results show that the best implementation of the usual strategies is better than that of Dengvaxia.

Equations

$$\frac{d{S}_{4a}}{dt}$$

$$\frac{d{S}_{ia}}{dt}$$

$$\frac{d{I}_{ia}}{dt}$$

$$\frac{d{R}_a}{dt}$$

$$\frac{d{S}_{4b}}{dt}$$

$$\frac{d{S}_{ib}}{dt}$$

$$\frac{d{I}_{ib}}{dt}$$

$$\frac{d{R}_b}{dt}$$

$$S_{iA}=\frac{S_{ia}}{N},I_{iA}=\frac{I_{ia}}{N},S_{iB}=\frac{S_{ib}}{N},I_{iB}=\frac{I_{ib}}{N},R_A=\frac{R_a}{N},\text{ and }{R}_B=\frac{R_b}{N},$$

$$\sum _1^4(S_{iA}+S_{iB}+I_{iA}+I_{iB})+R_A+R_B=1$$

$$S_{4A}=1-\left(\sum _{i=1}^3(S_{iA})+\sum _{i=1}^4(S_{iB}+I_{iA}+I_{iB})+R_A+R_B\right).$$

$$\frac{d{I}_{1A}}{dt}$$

$$\frac{d{S}_{iA}}{dt}$$

$$\frac{d{I}_{iA}}{dt}$$

$$\frac{d{R}_A}{dt}$$

$$\frac{d{S}_{4B}}{dt}$$

$$\frac{d{I}_{iB}}{dt}$$

$$\frac{d{S}_{iB}}{dt}$$

$$\frac{d{R}_B}{dt}$$

$$D$$

$$\displaystyle\hskip 28.45274ptI_{4B},R_{B})\in\mathbb{R}_{+}^{17}\,\,\bigg{|}\,\,S_{iA}\geq 0,I_{ia}\geq 0,\,\,\,R_{A}\geq 0,\,\,\,S_{iA}\geq 0,\,\,\,I_{iA}\geq 0,\,\,\,$$

$$R_B\ge 0,\text{where}i=1,2,3,4\text{and}$$

$$\displaystyle\hskip 56.9055pt\left(\displaystyle\sum_{i=1}^{3}(S_{iA}+S_{iB})+S_{4B}+\sum_{i=1}^{4}(I_{iA}+I_{iB})+R_{a}+R_{b}\right)\leq 1\bigg{\}}.$$

$$\sum _{i=1}^3\frac{d{S}_{iA}}{dt}+\sum _{i=1}^4\frac{d{I}_{iA}}{dt}+\sum _{i=1}^4\frac{d{S}_{iB}}{dt}+\sum _{i=1}^4\frac{d{I}_{iB}}{dt}+\frac{d{R}_A}{dt}+\frac{d{R}_B}{dt}<0.$$

$$\left(0,0,0,0,0,0,0,0,\frac{{ϵ}_a}{{ϵ}_a+{ϵ}_b+\mu },0,0,0,0,0,0,0,0\right).$$

$$R_0$$

$$+{\alpha }_{1b}\left(\frac{{ϵ}_1}{{ϵ}_1+{ϵ}_2+\mu }\right)\left(\frac{({ϵ}_1+\mu +{\gamma }_{1a}+{\beta }_{1a}+{ϵ}_2)}{({ϵ}_1+\mu +{\gamma }_{1a}+{\beta }_{1a})({ϵ}_2+\mu +{\gamma }_{1b}+{\beta }_{1b})-{ϵ}_1{ϵ}_2}\right),$$

$$F(I_{1A},I_{2A},\dots ,I_{4B})=$$

$$\left[\begin{array}{cccccccc}{\alpha }_{1a}x& {\alpha }_{1a}x& {\alpha }_{1a}x& {\alpha }_{1a}x& {\alpha }_{1a}x& {\alpha }_{1a}x& {\alpha }_{1a}x& {\alpha }_{1a}x\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ {\alpha }_{1b}{S}_{4B}^{\ast }& {\alpha }_{1b}{S}_{4B}^{\ast }& {\alpha }_{1b}{S}_{4B}^{\ast }& {\alpha }_{1b}{S}_{4B}^{\ast }& {\alpha }_{1b}{S}_{4B}^{\ast }& {\alpha }_{1b}{S}_{4B}^{\ast }& {\alpha }_{1b}{S}_{4B}^{\ast }& {\alpha }_{1b}{S}_{4B}^{\ast }\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]$$

$$V(I_{1A},I_{2A},\dots ,I_{4B})=\left[\begin{array}{cccccccc}{v}_{11}^{\prime }& 0& 0& 0& -{ϵ}_b& 0& 0& 0\\ 0& v_{22}^{\prime }& 0& 0& 0& -{ϵ}_b& 0& 0\\ 0& 0& v_{33}^{\prime }& 0& 0& 0& -{ϵ}_b& 0\\ 0& 0& 0& v_{44}^{\prime }& 0& 0& 0& -{ϵ}_b\\ -{ϵ}_a& 0& 0& 0& v_{55}^{\prime }& 0& 0& 0\\ 0& -{ϵ}_a& 0& 0& 0& v_{66}^{\prime }& 0& 0\\ 0& 0& -{ϵ}_a& 0& 0& 0& v_{77}^{\prime }& 0\\ 0& 0& 0& -{ϵ}_a& 0& 0& 0& v_{88}^{\prime }\end{array}\right]$$

$$v_{11}^{\prime }$$

$$v_{22}^{\prime }$$

$$v_{33}^{\prime }$$

$$v_{44}^{\prime }$$

$$V^{-1}(I_{1A},I_{2A},\dots ,I_{4B})=\left[\begin{array}{cccccccc}{d}_1& 0& 0& 0& b_1& 0& 0& 0\\ 0& d_2& 0& 0& 0& b_2& 0& 0\\ 0& 0& d_3& 0& 0& 0& b_3& 0\\ 0& 0& 0& d_4& 0& 0& 0& b_4\\ c_1& 0& 0& 0& d_5& 0& 0& 0\\ 0& c_2& 0& 0& 0& d_6& 0& 0\\ 0& 0& c_3& 0& 0& 0& d_7& 0\\ 0& 0& 0& c_4& 0& 0& 0& d_8\end{array}\right]$$

$$K=F(I_{1A},I_{2A},\dots ,I_{4B})V^{-1}(I_{1A},I_{2A},\dots ,I_{4B}).$$

$$\left|K-\lambda I_8\right|$$

$$=({\lambda }^6)\left[-\lambda {\alpha }_{1a}(1-S_{4B}^{\ast })(d_1+c_1)-\lambda {\alpha }_{1b}{S}_{4B}^{\ast }(d_5+b_1)+{\lambda }^2\right]$$

$$=({\lambda }^7)\left[-{\alpha }_{1a}x(d_1+c_1)-{\alpha }_{1b}{S}_{4b}^{\ast }(d_5+b_1)+\lambda \right]$$

$$R_0$$

$$+{\alpha }_{1b}\left(\frac{{ϵ}_a}{{ϵ}_a+{ϵ}_b+\mu }\right)\left(\frac{({ϵ}_a+\mu +{\gamma }_{1a}+{\beta }_{1a}+{ϵ}_b)}{({ϵ}_a+\mu +{\gamma }_{1a}+{\beta }_{1a})({ϵ}_b+\mu +{\gamma }_{1b}+{\beta }_{1b})-{ϵ}_a{ϵ}_b}\right).$$

$$\frac{d{S}_{4a}}{dt}$$

$$\frac{d{S}_{ia}}{dt}$$

Full text

AMS Subject Classification (2010): 92D30, 37N25, 49N90

Keywords: dengue, dengvaxia, SIR, optimal control

1 Introduction

Dengue is a a mosquito-borne viral infection for which there is no specific anti-viral treatment. Dengue virus is transmitted by female mosquitoes mainly of the species Aedes aegypti and, to a lesser extent, Aedes albopictus. There are four dengue virus serotypes, namely DEN-1, DEN-2, DEN-3, and DEN-4. Infection by one serotype confers life-long immunity to that serotype, but there is no cross-protective immunity to the other serotypes.

An effective technique used to detect dengue virus and the specific serotype has been developed. This laboratory test involves taking clinical samples and analyzing it through the use of Polymerase Chain Reaction (PCR) [9]. Though the method is very efficient, the laboratory test costs from 4,500 to 7,000 pesos which is expensive for an average Filipino [2]. Data shows that there were a total of 69,088 of dengue cases reported nationwide from January 1 to July 28, 2018 but only 301 cases were confirmed via PCR [8].

According to the World health Organization, 40 percent of the global population is estimated to be at risk of dengue fever [29]. Hence, there is a growing public health need for effective preventive interventions against dengue. Dengvaxia is the first dengue vaccine to be licensed. The vaccine was intended to prevent all four dengue types in individuals from 9 to 45 years of age living in endemic areas. In 2017, the manufacturer recommended that the vaccine only be used in people who have previously had a dengue infection, as outcomes may be worsened in those who have not been previously infected.This has caused a scandal in the Philippines where more than 733,000 children were vaccinated regardless of serostatus [30]. News have exposed deaths due to Dengvaxia and because of this the use of the vaccine is now permanently banned in the country.

In this paper, we propose a compartmental model describing the dynamics of dengue transmission considering its four strains. The model is age-structured taking into account the age requirement for the implementation of the Dengvaxia vaccine. Structuring our model this way helps us easily incorporate the control variable for Dengvaxia. There are existing models considering the four strains of dengue but our proposed model have the advantage of focusing only on the number of strains a person have been infected to and do not need any information on what specific strain(s) the person have been infected to. This is mathematically advantageous because the model now needs fewer equations than when one needs to be specific about the strains. This study aims to know how effective Dengvaxia could have been if it was implemented succesfully. This could be done by using optimal control theory on our model. Since Dengvaxia is now banned in the country, we also investigate on the other usual alternatives, like vector control and seeking early medical help and compare these with Dengvaxia.

2 Model Formulation

The proposed model considers eight susceptible classes depending on how many times a person was infected with dengue. The researcher considers primary, secondary, tertiary and quaternary infection considering the four serotypes of dengue. But unlike the model in [1], the specific kind of strain that infects the person will not matter anymore. The model divides the human population into two compartments: 9-45 years old human population and the rest of the human population taking into account the age implemetantation of the Dengvaxia vaccine. Pandey et al. [25] compared SIR model and vector-host model in dengue transmission and finds it that explicitly incorporating the mosquito population may not be necessary in modeling dengue transmission for some populations.

The total population at time $t$, denoted by $N(t)$, is subdivided into eighteen compartments. Age group $a$ consist of individuals of age less than 9 years or greater than 45 years. Age group $b$ consist of individuals of age 9 up to 45 years. Let $i=1,2,3,4$ and $j=a,b$. We denote by $S_{ij}(t)$ the number of individuals of age group $j$ who are susceptible to $i$ strain(s) at time $t$; $I_{ij}(t)$ the number of individuals of age group $j$ who were infected with any $i-1$ strain(s) in the past and is currently (at time $t$) infected with a different strain; and $R_{j}$ is the number of individuals of age group $j$ who had been infected by all the strains and has now (at time $t$) completely recovered. It is assumed that all individuals are born susceptible to any of the dengue strains and thus enter the class $S_{4a}$ and a recovered individual from a strain of dengue will have complete immunity to that strain. Humans leave the population through natural death rate and through per-capita death rate due to infection. Humans of age less than 9 years in the susceptible class $(S_{ia})$ leave the population and move to infected class $(I_{ia})$ or grow older to susceptible class $S_{ib}$. Furthermore, humans of age greater than 45 years in the susceptible class $(S_{ib})$ leave the population and move to infected class $(I_{ib})$ or grow older to susceptible class $S_{ia}$. In general, any member of the human population remains susceptible with the four strains of dengue and stays in the class $S_{4j}$ for a certain period unless bitten and infected with one strain. Once in the infective class $I_{1j}$, an individual may die or recover from that strain. If recovered from the first infection, susceptibility remains to the other 3 strains and then transfered to $S_{3j}$. If the same individual is infected with another strain, then he will be moved to the class $I_{2j}$. Once recovered from that strain, he is still susceptible to the other two strains and now belongs to the class $S_{2j}$. The cycle repeats until the individual gets infected four times and identified into the class $I_{4j}$. If ever the person survived the fourth infection, he then completely recovers and remains in the class $R_{j}$.

The parameters used are transmission coefficient ($\alpha_{ij}$), recovery rate ($\beta_{ij}$), birth rate ($\mu$), natural death rate ($\delta$), per-capita death rate through infection ($\gamma_{ij}$), rate of progression from age group $a$ to $b$ ($\epsilon_{a}$), and rate of progression from age group $b$ to $a$ ($\epsilon_{b}$). The following figure gives the flow diagram of the dengue transmission model.

The model can be mathematically described by the following system of 18 differential equations:

[TABLE]

where $I=\displaystyle\sum_{i=1}^{4}(I_{ia}+I_{ib})$ and $N=\displaystyle\sum_{i=1}^{4}(S_{ia}+S_{ib}+I_{ia}+I_{ib})+R_{a}+R_{b}$. Note that $\dfrac{dN}{dt}=\mu N-\delta N-\displaystyle\sum_{i=1}^{4}(\gamma_{ia}I_{ia}+\gamma_{ib}I_{ib}).$ To analyze the model, fractional quantities will be used and this could be done by scaling the population of each class with the total population. The new variables will be as follows:

[TABLE]

for $i=1,2,3,4$. Note that capital letters are used as subscripts to denote the scaled quantities. Furthermore, observe that

[TABLE]

or

[TABLE]

Using this and differentiating each fractional quantities with respect to time, the system (2) is reduced to the following system of 17 differential equations:

[TABLE]

where $I_{N}=\displaystyle\sum_{i=1}^{4}\left(I_{iA}+I_{iB}\right)$ and $J=\sum_{i=1}^{4}(\gamma_{ia}I_{iA}+\gamma_{ib}I_{iB})$.

This system of equations is epidemiologically and mathematically well-posed on the domain

[TABLE]

The space $\mathbb{R}_{+}^{17}$ denotes the positive orthant in $\mathbb{R}^{17}.$

2.1 Model Analysis

In this section, we obtain the disease-free equilibrium point and the reproductive number of the model. A disease-free equilibrium (DFE) is a steady state solution of an epidemic model with all infected variables equal to zero [6]. Over time, we want to achieve a disease-free state. A threshold that determines whether a disease-free state is achievable or not is the reproductive number $R_{0}$. The reproductive number is the expected number of individuals infected by a single infected individual over the duration of the infectious period in a population, which is entirely susceptible [11]. The reproduction number also gives an idea of which parameters in our model may be significant in our dynamics. The next generation operator approach defined by Diekmann et al. [12] and Driessche and Watmough [11] will be used to compute $R_{0}$.

Theorem 2.1

Assuming that the initial conditions lie in $D$, the system of equations (4) has a unique solution that exists and remains in $D$ for all time $t\geq 0$.

Proof: Note that the right-hand side of the system of equations (2) is continuous with continuous partial derivatives in $D$, by Cauchy-Lipschitz Theorem, the system of equations (2) has a unique solution. To show that $D$ is forward invariant, note that if $I_{iA}=0$ then $\dfrac{dI_{iA}}{dt}\geq 0$ for $i=1,2,3,4$. If $S_{iA}=0$ then $\dfrac{dS_{iA}}{dt}\geq 0$ for $i=1,2,3$. If $S_{iB}=0$ then $\dfrac{dS_{iB}}{dt}\geq 0$ for $i=1,2,3,4$. If $I_{iB}=0$ then $\dfrac{dI_{iB}}{dt}\geq 0$ for $i=1,2,3,4$. Note also that if $\displaystyle\sum_{i=1}^{3}(S_{iA}+S_{iB})+S_{4B}+\sum_{i=1}^{4}(I_{iA}+I_{iB})+R_{A}+R_{B}$ = 1, it follows that

[TABLE]

Hence, none of the orbits can leave $D$ and a unique solution exists for all time. $\blacksquare$

Theorem 2.2

The disease-free equilibrium point of the epidemic model (2) is

[TABLE]

Proof: Set the right-hand side of system (2) to zero and let $I_{iA}=I_{iB}=0$ for $i=1,2,3,4$. Since $\mu>0$ and the DFE point is in $D$, we have $S_{3B}=S_{3A}=S_{2B}=S_{2A}=S_{1B}=S_{1A}=R_{A}=R_{B}=0$. With these values, we are left with $\epsilon_{1}(1-S_{4B})-(\epsilon_{2}+\mu)S_{4B}=0$. Thus, we have $S_{4B}=\dfrac{\epsilon_{1}}{\epsilon_{1}+\epsilon_{2}+\mu}$. $\blacksquare$

Theorem 2.3

The reproductive number for system (2) is

[TABLE]

and the system is asymptotically stable if $R_{0}<1$ and unstable if $R_{0}>1$.

Proof: We form the matrices F and V where F is the matrix of the rates of appearance of new infections and V is the matrix of the rates of transfer of individuals out of the compartments. Let $X$ be the vector of infected classes and $Y$ be the vector of susceptible and recoverd classes.

Let $\mathcal{F}(X,Y)$ be the vector of new infection rates (flows from $Y$ to $X$). Let $\mathcal{V}(X,Y)$ be the vector of all other rates (not new infection). These rates include flows from $X$ to $Y$ (for instance, recovery rates), flows within $X$ and flows leaving from the system (for instance, death rates). The next generation operator formed is $K=FV^{-1}$ where $F=\left[\dfrac{\partial\mathcal{F}}{\partial X}\right],V=\left[\dfrac{\partial\mathcal{V}}{\partial X}\right]$. Evaluated at the disease-free equilibrium $(0,0,0,0,0,0,0,0,S_{4B}^{*},0,0,0,0,0,0,0,0)$ where $S_{4B}^{*}=\dfrac{\epsilon_{a}}{\epsilon_{a}+\epsilon_{b}+\mu}$, this becomes

[TABLE]

where $x=1-\dfrac{\epsilon_{a}}{\epsilon_{a}+\epsilon_{b}+\mu}=\dfrac{\epsilon_{b}+\mu}{\epsilon_{a}+\epsilon_{b}+\mu}$. We also have,

[TABLE]

where

[TABLE]

Moreover,

[TABLE]

where $d_{1}=\dfrac{v_{55}^{\prime}}{v_{11}^{\prime}v_{55}^{\prime}-\epsilon_{a}\epsilon_{b}},\,d_{2}=\dfrac{v_{66}^{\prime}}{v_{22}^{\prime}v_{66}^{\prime}-\epsilon_{a}\epsilon_{b}}$,   $d_{3}=\dfrac{v_{77}^{\prime}}{v_{33}^{\prime}v_{77}^{\prime}-\epsilon_{a}\epsilon_{b}},\,d_{4}=\dfrac{v_{88}^{\prime}}{v_{44}^{\prime}v_{88}^{\prime}-\epsilon_{a}\epsilon_{b}}$,

$d_{5}=\dfrac{v_{11}^{\prime}}{v_{11}^{\prime}v_{55}^{\prime}-\epsilon_{a}\epsilon_{b}},\,d_{6}=\dfrac{-v_{22}^{\prime}}{v_{22}^{\prime}v_{66}^{\prime}-\epsilon_{a}\epsilon_{b}}$,   $d_{7}=\dfrac{v_{33}^{\prime}}{v_{33}^{\prime}v_{77}^{\prime}-\epsilon_{a}\epsilon_{b}},\,d_{8}=\dfrac{v_{44}^{\prime}}{v_{44}^{\prime}v_{88}^{\prime}-\epsilon_{a}\epsilon_{b}}$,

$c_{1}=\dfrac{\epsilon_{a}}{v_{11}^{\prime}v_{55}^{\prime}-\epsilon_{a}\epsilon_{b}},\,c_{2}=\dfrac{\epsilon_{a}}{v_{22}^{\prime}v_{66}^{\prime}-\epsilon_{a}\epsilon_{b}}$,   $c_{3}=\dfrac{\epsilon_{a}}{v_{33}^{\prime}v_{77}^{\prime}-\epsilon_{a}\epsilon_{b}},\,c_{4}=\dfrac{\epsilon_{a}}{v_{44}^{\prime}v_{88}^{\prime}-\epsilon_{a}\epsilon_{b}}$,

$b_{1}=\dfrac{\epsilon_{b}}{v_{11}^{\prime}v_{55}^{\prime}-\epsilon_{a}\epsilon_{b}},\,b_{2}=\dfrac{\epsilon_{b}}{v_{22}^{\prime}v_{66}^{\prime}-\epsilon_{a}\epsilon_{b}},\,b_{3}=\dfrac{\epsilon_{b}}{v_{33}^{\prime}v_{77}^{\prime}-\epsilon_{a}\epsilon_{b}},\,b_{4}=\dfrac{\epsilon_{b}}{v_{44}^{\prime}v_{88}^{\prime}-\epsilon_{a}\epsilon_{b}}$.

Now, the next generation matrix $K$ is formed as

[TABLE]

Using cofactor expansion to get the eigenvalue of $K$, we have

[TABLE]

Hence we have, $\lambda=\alpha_{1a}x(d_{1}+c_{1})+\alpha_{1b}S_{4b}^{*}(d_{5}+b_{1})$ and $\lambda=0$ of multiplicity 7. Taking the dominant eigenvalue of matrix $K$, we get $R_{0}=\alpha_{1a}(1-S_{4B}^{*})(d_{1}+c_{1})+\alpha_{1b}S_{4B}^{*}(d_{5}+b_{1}).$ Substituting the values

$S_{4B}^{*}=\dfrac{\epsilon_{a}}{\epsilon_{a}+\epsilon_{b}+\mu}$, $x=\dfrac{\epsilon_{b}+\mu}{\epsilon_{a}+\epsilon_{b}+\mu}$, $d_{1}=\dfrac{v_{55}^{\prime}}{v_{11}^{\prime}v_{55}^{\prime}-\epsilon_{a}\epsilon_{b}}$, $d_{5}=\dfrac{v_{11}^{\prime}}{v_{11}^{\prime}v_{55}^{\prime}-\epsilon_{a}\epsilon_{b}}$, $b_{1}=\dfrac{\epsilon_{b}}{v_{11}^{\prime}v_{55}^{\prime}-\epsilon_{a}\epsilon_{b}}$, $c_{1}=\dfrac{\epsilon_{a}}{v_{11}^{\prime}v_{55}^{\prime}-\epsilon_{a}\epsilon_{b}}$, $v_{11}^{\prime}=\epsilon_{a}+\mu+\gamma_{1a}+\beta_{1a}$,

and $v_{55}^{\prime}=\epsilon_{b}+\mu+\gamma_{1b}+\beta_{1b}$, we finally have

[TABLE]

$\blacksquare$

3 Optimal Control Strategies

To mitigate the spread of dengue in the Philippines, three control intervention strategies $u_{i}(t)$, $i=1,2,3$ were identified to be incorporated in our model to see its best effects in reducing the number of infected individuals at its best implementation with the minimum cost. First is the transmission reduction conrol $u_{1}(t)$. This consists of holistic and effective methods in fighting mosquito at every stage of their life cycle, such as searching and eliminating breeding grounds, larvicide treatment, using adult mosquito trap, insecticides and thermal fogging during outbreaks. This also encompasses environmental management and sources reduction, self-protection measures which is just strengthening the capability of a person to avoid dengue and health education to provide awareness for their protection against dengue. Next strategy is proper medical care control $u_{2}(t)$. This means an effort for patients to seek early medication and immediate reporting to the nearest health care facility as symptoms are observed. With this control, chances of recovery increases specially from severe dengue. Last control is the Dengvaxia vaccine $u_{3}(t)$. This controversial Dengvaxia vaccine could have been a major breakthrough for dengue prevention in the country, but, the Philippines holds a lot of pending issues against its implementation. Moreover, it has been banned and not allowed for marketing in the Philippines. In this paper, we would like to know what could be the best effect of Dengvaxia (number of infected individuals are minimized, while the intervention costs are kept low) if it was properly implemented, that is, if it had been given only to susceptible population from 9-45 years old and to those who have had previous infection. We also investigate on the best possible impact of the other intervention strategies to compare it with Dengvaxia. This study is influenced by a similar study done for tuberculosis in the Philippines in [17].

The dengue state dynamics incorporating the different controls can now be described as follows:

[TABLE]

Transmission reduction control $u_{1}(t)$ is applied on all compartments having the transmission rates and proper medical control $u_{2}(t)$ is incorporated on all compartments involving recovery rates. The Dengvaxia control $u_{3}(t)$ is only applied on age group $b$ since it has restrictions on its implementation to people of age 9-45. Moreover, it is not incorporated to the compartment $S_{4b}$ since it must only be applied to people who have had previous infection.

The controls range is the interval $(0,1)$. The case when $u_{i}(t)\equiv 0$ means no control effort is employed while the case $u_{i}(t)\equiv 1$ indicates maximum administration of the control intervention. In this study, the optimal control problem aims to minimize the number of infected individuals with the minimum implementation cost of the control measure(s). It is assumed that the controls are quadratic functions to incorporate nonlinear societal cost associated with the implementation of control measures. This quadratic control is a common form of an objective functional in an optimal control problem. If three controls are considered $u_{1}(t),u_{2}(t)$ and $u_{3}(t)$, the objective functional to be minimized is represented by

[TABLE]

where $I(t)=\displaystyle\sum_{i=1}^{4}(I_{ia}(t)+I_{ib}(t))$. The bounds $T_{0}$ and $T_{1}$ are taken as 2020 and 2040, respectively, to consider a 20-year dengue intervention program and $B_{i}$’s, $i=1,2,3$, represent the weight constants associated to the relative costs of implementating the respective control strategy. These constants balanced the size and importance of each term in the integrand. In this work, we seek to find optimal controls $u_{1}^{*},u_{2}^{*}$ and $u_{3}^{*}$ satisfying

[TABLE]

where \Omega=\left\{(u_{1},u_{2},u_{3})\big{|}\,\,a\leq u_{i}(t)\leq b,\,\,\,u_{i}\in\mathcal{L}^{2}(2020,2040),\,\,i=1,2,3\right\}. Parameters $a$ and $b$ are the upper and lower bounds of the controls and are assumed to be 0.05 and 0.95, respectively. The weight parameters are set as $B_{1}=B_{2}=B_{3}=10^{6}$.

3.1 Main Theorem

The goal is to minimize the number of infected individuals and corresponding costs. From Pontryagin’s Maximum principle, optimal controls should satisfy the necessary conditions. Pontryagin’s Maximum Principle changes into a problem that minimize pointwise a Hamiltonian $H$, with respect to the control. We have $H=I(t)+\dfrac{B_{1}}{2}u_{1}^{2}(t)+\dfrac{B_{2}}{2}u_{2}^{2}(t)+\dfrac{B_{3}}{2}u_{3}^{2}(t)+\displaystyle\sum_{i=1}^{19}\lambda_{i}g_{i},$ where $g_{i}$ is the right hand side of the differential equation of the $i$th state variable. Applying Pontryagin’s Maximum Principle, we obtain the following theorem.

Theorem 3.1

There exist optimal controls $u_{1}^{*}(t),u_{2}^{*}(t)$ and $u_{3}^{*}(t)$ minimizing the objective functional

[TABLE]

Given these optimal solutions, there exist adjoint variables, $\lambda_{1}(t),\lambda_{2}(t),\ldots,\lambda_{19}(t)$, which satisfy

[TABLE]

[TABLE]

[TABLE]

with transversality conditions $\lambda_{i}(t_{f})=0$, for $i=1,2,\ldots,19.$ Furthermore,

[TABLE]

where

[TABLE]

and

[TABLE]

Proof: The existence of optimal controls $u_{1}^{*}(t),u_{2}^{*}(t),$ and $u_{3}^{*}(t)$ such that

$J(u_{1}^{*}(t),u_{2}^{*}(t),u_{3}^{*}(t))=\displaystyle\min_{\Omega}(u_{1},u_{2},u_{3})$, with state system (3) is given by the convexity of the objective functional integrand. By Pontryagin’s Maximum Principle, the adjoint equations and transversality conditions can be obtained. Differentiation of the Hamiltonian $H$ with respect to the state variable gives the following system:

$\dfrac{d\lambda_{1}}{dt}=-\dfrac{\partial H}{\partial S_{4a}},\,\,\,\,\dfrac{d\lambda_{2}}{dt}=-\dfrac{\partial H}{\partial I_{1a}},\,\,\,\,\dfrac{d\lambda_{3}}{dt}=-\dfrac{\partial H}{\partial S_{3a}},\,\,\dfrac{d\lambda_{4}}{dt}=-\dfrac{\partial H}{\partial I_{2a}},\,\,\,\,\\ \dfrac{d\lambda_{5}}{dt}=-\dfrac{\partial H}{\partial S_{2a}},\,\,\,\,\dfrac{d\lambda_{6}}{dt}=-\dfrac{\partial H}{\partial I_{3a}},\,\,\,\,\dfrac{d\lambda_{7}}{dt}=-\dfrac{\partial H}{\partial S_{1a}},\,\,\,\,\dfrac{d\lambda_{8}}{dt}=-\dfrac{\partial H}{\partial I_{4a}},\,\,\,\,\\ \dfrac{d\lambda_{9}}{dt}=-\dfrac{\partial H}{\partial R_{a}},\,\,\,\,\dfrac{d\lambda_{10}}{dt}=-\dfrac{\partial H}{\partial S_{4b}},\,\,\,\,\dfrac{d\lambda_{11}}{dt}=-\dfrac{\partial H}{\partial I_{1b}},\,\,\,\,\dfrac{d\lambda_{12}}{dt}=-\dfrac{\partial H}{\partial S_{3b}},\,\,\,\,\\ \dfrac{d\lambda_{13}}{dt}=-\dfrac{\partial H}{\partial I_{2b}},\,\,\,\,\dfrac{d\lambda_{14}}{dt}=-\dfrac{\partial H}{\partial S_{2b}},\,\,\,\,\dfrac{d\lambda_{15}}{dt}=-\dfrac{\partial H}{\partial I_{3b}},\,\,\,\,\dfrac{d\lambda_{16}}{dt}=-\dfrac{\partial H}{\partial S_{1b}},\,\,\,\,\\ \dfrac{d\lambda_{17}}{dt}=-\dfrac{\partial H}{\partial I_{4b}},\,\,\,\,\dfrac{d\lambda_{18}}{dt}=-\dfrac{\partial H}{\partial R_{b}},\,\,\,\,\dfrac{d\lambda_{19}}{dt}=-\dfrac{\partial H}{\partial N},$

with $\lambda_{i}(t_{f})=0$, for $i=1,2,\ldots,19.$

Optimal controls $u_{1}^{*}(t),u_{2}^{*}(t),$ and $u_{3}^{*}(t)$ are derived by the following optimality conditions:

[TABLE]

at $u_{1}^{*}(t),u_{2}^{*}(t),$ and $u_{3}^{*}(t)$ on the set $\Omega$. On this set

[TABLE]

where

[TABLE]

and

[TABLE]

Taking into account the bounds on controls, we obtain the characterization of $u_{1}^{*}(t),u_{2}^{*}(t),$ and $u_{3}^{*}(t)$. We have

[TABLE]

$\blacksquare$

3.2 Estimating Parameters

Pandey et al. [25] compared the SIR models and the vector-host models for dengue transmission and found that explicitly incorporating the mosquito population may not be necessary in modeling dengue transmission for some populations. In their paper, by comparing the equilibria of the vector host model and the SIR model they obtained the transmission coefficient $\beta$ in terms of the parameters of the vector host model: $\beta\approx\dfrac{mc^{2}\beta_{H}\beta_{V}}{\mu_{V}}$ where $m$ is the number of mosquitos per person, $c$ is the biting rate, $\beta_{H}$ is the mosquito-to-human transmission probability, $\beta_{V}$ is the human-to-mosquito transmission probability, $\mu_{V}$ is the mosquito mortality rate and $\beta$ is the composite human to human transmission rate. Note that in our paper, we use $\alpha$ as our transmission coefficient. From [7], we can get the following values: biting rate $c=1$, mosquito-to-human transmission probability $\beta_{H}=0.375$, human-to-mosquito transmission probability $\beta_{V}=0.75$, mosquito mortality rate $\mu_{V}=0.1$. From [21], we also get the value of the number of mosquito per person $m$ =1. Substituting these parameter values to $\alpha$, we have

[TABLE]

This now becomes the value for $\alpha_{1a}$ and $\alpha_{1b}$.

Data from [8] reveals that 54 % of the dengue cases is caused by DENV 3, 25 % is caused by DENV 1, 18 % is caused by DENV 2 and 3 % of the dengue cases is caused by DENV 4. Assuming that this distribution directly affect the transmission coefficients for the second, third and fourth infections, we can arrive at the following estimates for the other transmission coefficients: $\alpha_{2a}=0.6125999\times\alpha_{1a}$, $\alpha_{3a}=0.195696\times\alpha_{1a}$, $\alpha_{4a}=0.017496\times\alpha_{1a}$, $\alpha_{2b}=0.6125999\times\alpha_{1b}$, $\alpha_{3b}=0.195696\times\alpha_{1b}$, and $\alpha_{4b}=0.017496\times\alpha_{1b}$. Thus, we have $\alpha_{2a}=1.7229$, $\alpha_{3a}=0.550395$, $\alpha_{4a}=0.0478125$, $\alpha_{2b}=1.7229$, $\alpha_{3b}=0.550395$, and $\alpha_{4b}=0.0478125$.

From Syafruddin and Noorani model of dengue fever in [23], the recovery rate of infected humans is 0.32833. But data from [8] shows that the recovery of people of age less than nine years old, where most dengue cases occur, is very low . Thus, $\beta_{1a}$ should be lower compared to $\beta_{1b}$. We estimate $\beta_{1b}=0.32833$ and $\beta_{1a}=0.30$. However, dengue can become severe in the next infections thus we may have $\beta_{2b}=0.164165$, $\beta_{3b}=0.164165$, and $\beta_{4b}=0.164165$. We also have $\beta_{2a}=0.15$, $\beta_{3a}=0.15$, and $\beta_{4a}=0.15$. From the same literature, we also obtain the value for the death rate through infection given by $\gamma=0.0000002$. Moreover, data from the Department of Health shows that the dengue deaths is more prevalent in the age group 0-9 years old. This would imply that, $\gamma_{1a}$ should be higher compared to $\gamma_{1b}$. We let $\gamma_{1b}=2.0\times 10^{-7}$ and $\gamma_{1b}=3.0\times 10^{-7}$. In the next infections, dengue can become severe and may lead to more number of deaths . Thus, we let $\gamma_{2b}=\gamma_{3b}=\gamma_{4b}=4.0\times 10^{-7}$ and $\gamma_{2a}=\gamma_{3a}=\gamma_{4a}=6.0\times 10^{-7}$.

From [7], we have the birth rate $\mu$ equal to $8.5\times 10^{-4}$ and death rate $\delta$ equal to $4.5\times 10^{-4}$. Furthermore, we estimate the growth rate from age group $a$ to age group $b$ to be $\epsilon_{1}=6.92\times 10^{-5}$ and the growth rate from age group $b$ to age group $a$ to be $\epsilon_{2}=6.92\times 10^{-5}$.

3.3 Numerical Simulations

For our numerical simulations, we use the following initial values:

$S_{4a}=2.1\times 10^{7}$,   $I_{1a}=7.0\times 10^{6}$,   $S_{3a}=8.0\times 10^{6}$,   $I_{2a}=4.0\times 10^{6}$,   $S_{2a}=4.5\times 10^{6}$,   $I_{3a}=2.5\times 10^{6}$,   $S_{1a}=1.0\times 10^{6}$,   $I_{14a}=3.5\times 10^{5}$,   $R_{a}=1.5\times 10^{6}$,   $S_{4b}=2.1\times 10^{7}$,   $I_{1b}=1.4\times 10^{7}$,   $S_{3b}=1.0\times 10^{7}$,   $I_{2b}=6.0\times 10^{6}$,   $S_{2b}=4.5\times 10^{6}$,   $I_{3b}=2.5\times 10^{6}$,   $S_{1b}=1.5\times 10^{6}$,   $I_{4b}=6.5\times 10^{5}$,   and $R_{b}=3.5\times 10^{5}$.

With the estimated parameters above and the set of initial conditions, and being assured of the existence of the optimal controls $u_{1}^{*}(t),u_{2}^{*}(t),$ and $u_{3}^{*}(t)$ by Theorem 4.1.1, we performed the numerical simulations. The results are given by the following figures. On each of the figures, the graph on the left shows how the control should be implemented to have the minimum value of the objective functional and the graph on the right shows the corresponding effect of the implementation on the total number of infected individuals. The black graph represents the total number of infected individuals over time when no control is used.

3.4 Discussion

As expected, we can see a significant decrease in the total number of infected individuals if Dengvaxia is implemented with the optimal strategy. However, one can clearly see that the optimal strategies for the other alternative controls have better impacts than that of Dengvaxia. The best impact if only one control is to be used can be seen in the optimal strategy for the transmission reduction control. These results is telling us that even without Dengvaxia we can still reduce the number of infected individuals and the reduction could be significantly better. It is interesting to note that the transmission reduction control is better than the proper medical care control, showing that prevention is still better than cure.

However, there is a very important point that should be made here. Note that the effects can only be seen if the optimal strategies are used. For the transmission reduction and proper medical care control, the graphs on the left shows us that the optimal strategy is almost maximum administration of the control for the duration of the program. This means that almost all of the affected population should do their part in implementing these strategies, if possible, at all times. This clearly calls for discipline, and this is where the Dengvaxia control has the advantage. The vaccination does not need discipline from the population to be implemented. Therefore, a concerted effort from the people and the government is needed to effectively control dengue without the aid of vaccination like Dengvaxia, and doing so seems rewarding because the expected effect seems to be way better than the effect of Dengvaxia and the implementation will not cause any Dengvaxia-like public fear.