Optimal Control of Dengue Transmission with Vaccination

: Dengue disease is caused by four serotypes of the dengue virus: DEN-1, DEN-2, DEN-3, and DEN-4. The chimeric yellow fever dengue tetravalent dengue vaccine (CYD-TDV) is a vaccine currently used in Thailand. This research investigates what the optimal control is when only individuals having documented past dengue infection history are vaccinated. This is the present practice in Thailand and is the latest recommendation of the WHO. The model used is the Susceptible-Infected-Recovered (SIR) model in series conﬁguration for the human population and the Susceptible-Infected (SI) model for the vector population. Both dynamical models for the two populations were recast as optimal control problems with two optimal control parameters. The analysis showed that the equilibrium states were locally asymptotically stable. The numerical solution of the control systems and conclusions are presented.


Introduction
The dengue epidemic first occurred in the Philippines in 1954. It reached Thailand in 1958. The disease is caused by an infection by any one of the four serotypes of dengue virus, which are labeled as DEN-1, DEN-2, DEN-3, and DEN-4. The dengue viruses are transmitted by two species of the Aedes mosquitoes, the Aedes aegypti and the Aedes albopictus. All four serotypes have a common antigen, resulting in cross-reaction and cross-protection of the four serotypes. The cross protection is not permanent. A person infected by one of the serotypes will have permanent immunity to that serotype, but only partial immunity to the other three. Some of the immunity will last for a short period, approximately 6-12 months. Those people might be re-infected if they happen to meet a different serotype of dengue virus. This second infection is different from the initial infection and is labeled as a secondary dengue infection [1][2][3][4] since the symptoms and outcomes of the primary and secondary dengue infections can be quite different. In some individuals (infants or young children), infection by the dengue virus may lead to undifferentiated fever (uf). The individuals are said to have the viral syndrome of dengue fever, which can only be detected through laboratory tests. In older children and adults, infection by the dengue viruses leads to what is usually labeled as dengue fever (DF). People with DF exhibit symptoms such a mild fever, headaches, pain around the eyes, muscular pain, and pain in the bones. If an individual experiences the clinical symptoms of high fever accompanied by bleeding, enlarged liver, and severe shock, he is said to have dengue hemorrhagic fever (DHF). During the fever, there will also be a low platelet count and plasma leakage. If large amounts of plasma leak out, the patient will have a shock condition called dengue shock syndrome (DSS) [5][6][7]. The last two (DHF and DSS) are the symptoms that the individuals in the secondary dengue infection group experience. These symptoms can be viewed as an allergic reaction.
Since millions of people can be infected by the dengue virus, there is an economic cost to these people becoming sick, and vaccines have been developed. Chimeric yellow fever dengue tetravalent dengue vaccine (CYD-TDV) is one of these vaccines. It was first registered as a dengue vaccine in Mexico. It is now registered in 13 countries around the world, including Thailand, which had a role in its development. This dengue vaccine is the first and only vaccine in the world at the moment to cover all four strains of the dengue virus and is called Dengvaxia ® , which is developed by Sanofi Pasteur to help protect against the dengue disease. The vaccine has an overall effectiveness of 56.5%, which is 74% more effective in older children, 12-14 years old, and 75% effective for DEN-3 and DEN-4 [8,9]. It is reported that the efficacy of the vaccine is higher in children who have previously been infected with dengue fever. The dengue vaccine reduces the severity of the disease by 88.5% and hospitalization by 67.2% [10][11][12]. In December 2017, the WHO [13] issued a new recommendation that states that the WHO recommends vaccination (with Dengvaxia) only in individuals with a documented past dengue infection. This should be taken into consideration in any models used.
Esteva and Vargas [4,14] were among the first to study the transmission of dengue disease. They developed a mathematical model in which there were compartmental models for both the human and mosquito populations. The human population was described by a Susceptible-Infected-Recovered (SIR) model while the mosquito population was described by an SI (no recovered) model. Pongsumpun et al. [15][16][17] have also studied the transmission of dengue virus. Most of Pongsumpun's work has been centered on the situation in Thailand since the dengue fever is of major concern to Thailand. She has included an exposed class (E) to the model, making the Susceptible-Infected-Exposed-Recovered (SEIR) model, to describe the dynamics of the human population. In Ref. [16], the author used the SIR model to simulate the possible outcomes of vertical transmission of the virus among mosquitos. She and her coworkers [17] included vertical transmission in a SEIR model. Syafruddin and Noorani [18] studied the mathematical model for dengue transmission and applied it to the situations in Indonesia and Malaysia. Yaacob [19] studied the mathematical model of the dengue disease in people who have no immunity.
Singh et al. [20] and Tasman et al. [21] considered the effects of vaccination on a model in which the human population is divided into children and adults. They also considered that there were two types of infections, primary and secondary dengue infection. It was assumed that individuals experiencing a secondary infection were at a higher risk. In these studies, the adults were further divided in two groups, so the human population consisted to three groups: less than 9 years, between 9 and 45 years, and between 45 and 65 years [22,23]. Using a similar model to study the transmission of another disease, melioidosis, Tavaen and Viriyapong [24] studied the local and global stability analyses and optimal control for this disease. There are many studies about the effects of the dengue vaccination on the spread of the dengue disease. [25][26][27][28]. They have introduced various models to simulate the dynamic of the programs when there is complete vaccination, random mass vaccination, imperfect random mass vaccination, and random mass vaccination with waning immunity levels. They have used optimal control strategies to simulate the results of the programs.
The number of cases and deaths by month and the number of cases and deaths each year in Thailand from 2003-2020 data from the Bureau of Epidemiology at the Ministry of Health is shown in Figure 1. It can be seen in the figure that the dengue fever is prevalent in the rainy season from June to September. The incidence is the highest in July. The number of cases, which fluctuated month to month, tended to increase yearly from 2003 to 2020. The same is true for the number of deaths. When the number of cases is high, there will be more deaths. The percentage of deaths is very small. Using the sources that gave these results, we were interested in the outcome of a vaccination program in which only individuals with a documented past dengue infection (i.e., an individual who would have a secondary dengue infection if bitten by a mosquito infected with a different serotype of the virus) are vaccinated. We used the double Susceptible-Infected-Recovered (SIR) model for the human population and the Susceptible-Infected (SI) model for the vector population. The analysis of the stability of the model was carried out by using dynamic analysis. The Routh-Hurwitz criteria were applied to analyze the system model for stability. The reproductive number was calculated. The optimal control theory was applied in the transmission model in order to minimize the number of infected humans with primary and secondary infections. Numerical simulation was performed. Results and conclusion are presented in this paper.

Mathematical Model
The basic mathematical model was the SEIR model presented in ref. [15]. The equations in that model describe the dynamics of the spread of dengue fever when there is only one serotype of the virus present. In this work, the basic SEIR model of [15] was extended to include secondary infection of a different serotype, whereby the members of the recovered population become the susceptible members in the second SIR, effectively providing a framework for describing a vaccination program in which only people who have been infected are considered. In Thailand, the medical status of each Thai citizen is kept at the District Office in each province in the country. It is easy to determine from the past medical histories anyone who was infected with the dengue virus. Dengue fever is one of the five diseases that must be reported to the Thai Ministry of Health. The susceptible human population in the second SIR model used here are the not sick humans who have been infected by the serotype A virus, since they will be the only ones given the vaccine. A person who has no prior history of any dengue infection is not considered to be a candidate for the vaccination. The vector population was divided into two compartments: susceptible and infected (SI). The infected mosquito was the subset of infected mosquitoes transmitting virus B. The human population was subdivided into six population groups. It should be remembered that all of the recovered individuals have the antibodies to a particular serotype of the dengue virus at the end of primary infection. Susceptible people of this kind are not born into this group; they emerge after several months of being infected by a serotype virus. The vector population was classified into two subclasses. The variables are defined in Table 1. Table 1. The definition of the variables used in the differential system of Equations (1)- (10). The dynamic transmission of dengue disease with the vaccination model is shown in Figure 2. The dynamics of human and vector populations and the system of differential equations are given by:

S HP
with the conditions: S HP + I HP + R HP + S HS + I HS + R HS = N H (9) where the parameters of Equations (1)-(10) are defined in Table 2. Table 2. The definition of the parameters of the differential system of Equations (1)- (10).

Parameters Definition α
The biting rate of the vector population ψ The vaccine efficiency θ The recurrent infection rate N H The total number of humans in the study population N V The total number of vectors in the study population β H The transmission rate of dengue virus from vector to human β V The transmission rate of dengue virus from human to vector b The birth and natural mortality rate of the human population d V The natural mortality rate of the vector population d d The mortality rate from infection of the human population d k The mortality rate from infection of the vector population γ P The recovery rate of those with a primary infection γ S The recovery rate of those with a secondary infection The rate of change of both the total population of humans and vectors is zero and given by: with conditions, we get: Normalizing the equations by introducing the following normalized variables: with the additional condition: S HP + I HP + R HP + S HS + I HS + R HS = 1 (17) Mathematics 2021, 9, 1833 6 of 33 The mathematical model of Equations (1)- (8) is now reduced to the following equation:

The Equilibrium Points
Definition 1 ([17]). The point X ∈ R n is an equilibrium point for the differential equation Since epidemiological models are inherently dynamical systems, the knowledge of the equilibrium points is vital for determining the behavior of long-term dynamics. Towards this goal, the most important parameter for determining whether an outbreak will occur or not is the basic reproductive number R 0 . The equilibrium points are obtained by setting the right-hand side of Equations (19)- (24) to zero. This system model now admits two equilibrium points, namely the disease-free point and an endemic equilibrium point. The disease-free equilibrium point E 1 is: The endemic equilibrium point E 2 = S * HP , I * HP , R * HP , S * HS , I * HS , I * V is: where:

The Basic Reproductive Number
Definition 2 ([24]). Basic reproductive number (R 0 ) is defined as the average number of secondary infections when a single infective enters an entirely susceptible population.
The basic reproductive number (R 0 ) is obtained using the next-generation matrix method [30][31][32]. We selected I HP , I HS , and I V to be the classes to construct the F and V matrices, which are important to this method. For our system, the matrices F and V contain new infection terms and transition terms. We evaluated the Jacobian matrices F and V at the disease-free equilibrium point E 1 = (1, 0, 0, 0, 0, 0), where F is non-negative and V is non-singular. The F (gains) and V (losses) matrices are: The required matrices are: R 0 is the dominant eigenvalue of the matrix R; then, the value of the basic reproductive number is given by:

Local Stability of Equilibrium Points
has all its eigenvalues with negative real parts. The equilibrium point E 0 is not stable if at least one of the eigenvalues of the matrix J has a positive real part.
The local stability of each equilibrium point states of this model is determined from the Jacobian matrix at that equilibrium point of the system of Equations (19)- (24). The Jacobian matrix is: Theorem 1. At the equilibrium point E 1 , the disease-free state is locally asymptotically stable when R 0 < 1.
Proof. See the proof of Theorem A1 in Appendix A.
Theorem 2. The equilibrium point E 2 is locally asymptotically stable when R 0 > 1.
Proof. See the proof of Theorem A2 in Appendix A.

Numerical Simulation
In this section, the numerical analysis of the transmission of dengue disease with the vaccination will only consider individuals who have a documented history of past dengue infection or have died from the infection. The parameter values within this model are listed in Table 3. Note that though most of the parameters are taken from the literature, some of the values need to be assumed for the purpose of this investigation. The use of N H and N V values of 10,000 for both cases reflects a small rural town, for example, one like Mae Hong Son, with a lower socio-economic status, which is very useful for investigating the efficiency of the vaccination program. Note that we also assume a constant vector population, which is appropriate for this case since the town is small, and far away from larger cities such as Bangkok. The value of d d will be assumed to be 1/180 for both the disease-free and endemic equilibrium cases. This value is close to the reported values of 500 deaths per population of 100,000 [33]. The value of d k is assumed to be the same as the reported natural mortality rate. The model is then simulated with the use of the Runge-Kutta differential equation solver in MATLAB. The numerical results are shown in Figures 3-6. Table 3. The parameters used in the numerical simulation.

Parameters
Disease-Free Endemic References α 1/7 1/7 [1,9], [15][16][17], [29,33] [1,9], [15][16][17], [29,33] [1,9], [15][16][17], [29,33] For the case of larger cities such as Chiang Mai or Khon Kaen, etc., we assume that: Note that Case 2 simulates the situation where the size of the vector population is not constant but is rather a function of time. The exponential power of −1/14 is also used to reflect the natural death rate of the vector population. It can be seen from Figures 7 and 8 that there is a small change in trajectory from Case 1 to Case 2, which is to be expected since the total vector population slowly decays.

Sensitivity Analysis of Parameters
The index tells us which parameters have a high impact on the basic reproductive number, R 0 , and should be targeted by intervention strategies. In addition, the sensitivity index allows us to measure relative changes in variables as parameters change. The normalized forward sensitivity index of a variable related to a parameter is the ratio of the relative change in the variable to the relative change in the parameter [34].
Definition 4 (Chitnis, Hyman, and Cushing [34]). The normalized forward sensitivity index of R 0 , which depends differentiably on a parameter p, is defined by: Given an expression for R 0 in Equation (27), the sensitivity index of Equation (29) can now be used to evaluate the index of each parameter. Note that the sensitivity index could be a function of some parameters, or a number, signifying its independence of any parameter. Table 4 shows the sensitivity index of R 0 to all the parameters of Table 3. Table 4. Sensitivity indices R 0 p evaluated at the baseline parameter values of Table 1. Figure 9 plots a comparison of the R 0 values against the changes in d d and γ P . It can be seen that as the d d and γ P values change, the increment of R 0 in d d is more pronounced than that of γ P . This is reflected in the

Parameters Sensitivity
Note that a similar definition to Equation (29) with regards to the sensitivity of I * V with respect to the parameters could also be given. However, due to the algebraic complexities of I * V in Equation (26), a numerical plot is instead given. It is apparent in Figure 10 that as β V changes, the endemic equilibrium I * V also changes somewhat linearly.

The Optimal Control Problem
In this section, the solution of the optimal control problem is discussed. The Pontryagin Minimum Principle (PMP) theory is used to solve this problem. Equations (19)-(24) will be recast as a control problem. The purpose of this is to recast the problem as one of minimizing the number of infected humans to achieve an optimal outcome. Since, the system consists of two dynamics, one for the humans and the other for the vectors, two control parameters will be needed, u 1 for the human population and u 2 for the vector population. u 1 is the vaccination rate and u 2 is the rate at which the breeding of the Aegypti mosquitoes is annihilated. This model can be written as the system of the equation as follows: In terms of the normalized compartments, the differential equations become: All parameters retain same definitions as before. The optimal control problems of Equations (38)-(43), require a definition of the objective function given as: with initial condition S HP ≥ 0, I HP ≥ 0, R HP ≥ 0, S HS ≥ 0, I HS ≥ 0, R HS ≥ 0, S V ≥ 0, and I V ≥ 0. The constants X 1 , X 2 , X 3 , and X 4 are weight constants and the terms X 3 u 2 1 (t) and X 4 u 2 2 (t) represent the costs associated with the control variables u 1 and u 2 , respectively. We can assign an optimal solution of this model optimal control problem by using the Lagrangian and the Hamiltonian of the problems. The Lagrangian of the optimal control problem is given by: Theorem 3. We consider the objective function J given by Equation (44) with (u 1 , u 2 ) ∈ U subjecting to the control system of Equations (38)-(43) with initial condition. There exists u Proof. See the proof of Theorem A3 in Appendix A.

Theorem 4.
There exists the adjoint variables λ 1 , λ 2 , λ 3 , λ 4 , λ 5 , and λ 6 under the control that satisfies the following: with the boundary conditions: In addition, the optimal control variables are given by: Proof. See the proof of Theorem A4 in Appendix A.
The simulation results for the optimal states are presented in Figures 11-15 and the optimal controls are shown in Figure 14 using parameter values according to Table 1. Figures 11-15 show the simulation results of the system of Equations (38)-(43) with and without controls of S HP , I HP , R HP , S HS , I HS , and I V . The plot of Figure 11 is obtained by setting the weight X 1 to be equal to X 2 . Figure 12 presents the case in which the weight X 1 less than X 2 . Figure 13 presents the case where X 1 is at least 10 times greater than X 2 . Figure 14 presents the case in which the weight X 1 is greater than X 2 . Figure 15 presents the case where X 1 is at least 10 times greater than X 2 . Note that the main goal of the control is to minimize the number of infected humans with primary infection I HP , as well as the secondary infection I HS . The main emphasis, however, is on the secondary infectious individuals. For each of these scenarios, a comparison is also made with the case where no control is applied. If X 1 = X 2 , the convergence time to equilibrium is significantly faster than if X 1 is less than X 2 , X 1 is much less than X 2 , X 1 is greater than X 2 , and X 1 much greater than X 2 . Moreover, as the weight of X 1 increases, the steady state of infected individuals of the secondary population I * HS gradually decreases to zero. The trajectory of the infected vectors gradually deviates from the no-control case. No significant change seems to occur to the trajectories as the weight X 2 is increased, except for the infected vectors, which appear to approach zero for a large weight of 100. With no control measures, it will take longer for equilibrium to be reached.    Figure 16a shows that to maintain the optimal control of the infected population, u 1 (t) would have to be at 50%, 70%, and 90% for the first 25 days, after which the control required will show an exponential decline to zero. Figure 16b shows the values required to maintain optimal control when u 2 (t) is used as the controlling factor. The application of 50%, 70%, and 90% of u 2 (t) for achieving optimal control should be maintained for the first 47 days, after which the amount of control steadily decreases to zero.

Discussion and Conclusions
In this paper, we analyzed the effects of different vaccination strategies to prevent secondary dengue infection in order to reduce the severity of the disease. The dynamic of the dengue fever transmission model assumes that the human and vector population are constant. The analysis was based on the use of the Routh-Hurwitz criteria to establish local asymptotic stability. The equilibrium points that we found were disease-free convergence to E 1 = (1, 0, 0, 0, 0, 0) and the endemic equilibrium point. The basic reproductive number was defined as R 0 . If R 0 is less than one, then the disease-free state exits and is locally asymptotically stable; it is unstable when R 0 is greater than one. We simulated the numerical solution of the parameters with different values, as shown in Figure 5. We can see that if the transmission rate of dengue virus from vector to human β H is large, then the convergence to an equilibrium point becomes slower for susceptible humans to primary and secondary infection. The number of humans recovered from primary infection, the number of humans with primary or secondary infection, and the number of infected vectors will converge to an equilibrium point more rapidly. Likewise, in Figure 6, we see that if the transmission rate of the dengue virus from human to vector β V is large, then there is a slower convergence to an equilibrium point of the susceptible humans. Similarly, humans who recovered from a primary infection, humans infected with a primary or secondary infection, and infected humans will converge to an equilibrium point more quickly, i.e., E 2 (S HP = 0.00003, I HP = 0.00022, R HP = 0.00810, S HS = 0.00002, I HS = 0.00031). However, the infected vector will converge to a different equilibrium point. In addition, this will also make the basic reproductive number R 0 greater than one. Changing β V to a different value affects the convergence time for reaching the equilibrium point, as shown in Figure 8. To investigate whether there is a limitation on the model when there is a non-constant total vector population N V , two further cases were also considered. The first case had the largest total human population N H at 500,000, and a large N V of 100,000. The second case kept N H the same, whilst the N V was treated as an exponential function of time. Results showed that although having a non-constant N V does have some effect on the trajectories, such a change is quite small compared to that of the constant vector population case.
We adopted an optimal control approach using the vaccination rate and a rate for destroying the breeding of the Aegypti mosquito in order to minimize the number of humans with primary or secondary infections. To do this, we used the Pontryagin Minimum Principle (PMP) method to solve the optimal control problem with conditions X 1 equal X 2 , X 1 less than X 2 , X 1 much less than X 2 , X 1 greater than X 2 , and X 1 much greater than X 2 . We can see that, if there are no controls, the number of humans with primary or secondary infections increases. This will cause the solution to converge to the equilibrium point over a longer time. With the controls in place, the number of mosquitoes in the infected vector population and the number of people who need to be vaccinated will decrease over time, as shown in Figures 11-18.  Table 3. Using Routh-Hurwitz criteria [35,36] for n = 4, the endemic equilibrium point is stable if conditions (i) − (iv) below are satisfied. Since symbolic computation may be difficult, to illustrate whether conditions (i) − (iv) are indeed satisfied due to the algebraic numerical complexities, numerical simulations of conditions (i) − (iv) are instead given. It is seen from Figure 1A that conditions (i) − (iv) are indeed satisfied. This means that the polynomial of Equation (A3) is Hurwitz.
Proof. We apply the existence of an optimal control problem from [37,38].
The control set U is closed and convex by its definition above and the integrand of the function Equation (38) is also convex in U. It is obvious that these states and the control variable are nonnegative. Since the solution to the systems given by Equations (38)-(43) are bounded, the control function will be convex in U. Let q 1 and q 2 be two positive constants and ζ > 1. If we now set q 2 = min(I HP (t), I HS (t)), q 1 = min(X 3 , X 4 ), and ζ = 2, the Lagrangian function L can be rewritten as: L(I HP , I HS , u 1 , u 2 ) = X 1 I HP + X 2 I HS + 1 2 X 3 u 2 1 (t) + X 4 u 2 2 (t) ≥ q 2 (I HP + I HS ) + q 1 |u 1 | 2 + |u 2 | 2 = q 2 + q 1 |u 1 | 2 + |u 2 | 2 (A4) The optimal control of this model is obtained by applying Pontryagin's Minimum Principle [38].