Sub-Optimal Control in the Zika Virus Epidemic Model Using Differential Evolution

A dynamical model of Zika virus (ZIKV) epidemic with direct transmission, sexual transmission, and vertical transmission is developed. A sub-optimal control problem to counter against the disease is proposed including three controls: vector elimination, vector-to-human contact reduction, and sexual contact reduction. Each control variable is discretized into piece-wise constant intervals. The problem is solved by Differential Evolution (DE), which is one of the evolutionary algorithm developed for optimization. Two scenarios, namely four time horizons and eight time horizons, are compared and discussed. The simulations show that models with controls lead to decreasing the number of patients as well as epidemic period length. From the optimal solution, vector elimination is the prioritized strategy for disease control.


Introduction
Despite being discovered in Africa and named after a forest in Uganda in 1947, there are only several reported cases of Zika virus (ZIKV) infection before 2007.The virus itself is a mosquito-borne virus from the genus flavivirus that is transmitted primarily through the bite of Aedes mosquitoes.At first, the ZIKV infection was not considered a threat due to its mild symptoms such as rashes, headache, malaise, and fever [1].Besides, the mortality rate of ZIKV infection is considered negligible, unlike other Flavivirus infections such as dengue and yellow fever.However, the announcement [2] on the possible association between the ZIKV infection and the increasing number of newborn babies with microcephaly in 2015 raised the alarm and public awareness of this outbreak.There is some correlative evidence reported from Brazil, Colombia, and other Latin America countries [3,4].
Some previous studies suggest that ZIKV could pass from pregnant mothers to their baby through the intrauterine infection, and microcephaly also happened during this process.Beside vector-mediated transmission and intrauterine (or vertical) transmission, ZIKV could also transmit through sexual transmission [5].Apart from being detected in blood and semen [6], researchers found that the virus persists in semen and urine after it disappeared from the serum [7].Moreover, it could pass on even before the patient develops the symptoms [8].To control the outbreak, it can be seen that only relying on reducing the number of mosquito population or preventing the contact between host and vectors might be insufficient.The control strategy should include reducing the sexual transmission and vertical transmission as well.
A dynamical model was a popular tool for describing and analyzing various diseases.It was used to describe the mechanism functions and causes abnormality in some non-infectious diseases such as diabetes [9,10].In prior publications [11][12][13], models with systems of ordinary differential equations were utilized for describing the behavior of infectious virus such as Zika.Optimal control is one of the popular methods for finding the solution that will minimize the objective function value.According to some publications [14,15], optimal control can be used to find the best strategy associated with the implemented cost and the provided situations.The optimal solutions for some flavivirus control such as chikungunya [16] and dengue [17] have been studied.
In this paper, we develop the ZIKV virus infection model using some previous works [12,13,18].The model consists of three populations: adult humans, newborn babies human, and mosquitoes.Three controls are selected from the suggestions of the World Health Organization(WHO) [19] and health authorities [20]: using larvicide and adulticide to reduce the mosquito population, using insect repellents to prevent mosquito biting, and abstaining from sexual activity or using protection to prevent pregnancy.
Regarding optimal control and sub-optimal control problem, there are direct and indirect methods for solving the problems.The indirect approaches, such as Pontryagin's Maximum Principle [21] and shooting method, are popular to find optimal solutions of some flavivirus control problems, for example, chikungunya [16] and dengue [17].Previously, the authors also studied the optimal control problem of the Zika virus infection [18] by the indirect approach.However, the implementations for health authorities by these studies [16,17] are difficult to practically apply since the policies have to be varied continuously.Considering this drawback, we propose a sub-optimal control problem and find the solution that is both close to optimal and practical for real situation.
A sub-optimal control problem can be proposed through Model Predictive Control (MPC).MPC is the mathematical tool for controlling the optimization process in various fields such as industrial and medical technology [22].It can be used from the simple to complex dynamic process, and the result control law is easy to implement.In our work, the control framework consists of a nonlinear system without constraints on spending cost or the final result of the system variables.This framework is related to the nonlinear unconstrained MPC schemes, which are those without terminal constraints and cost [23].Furthermore, the study of Herty [24] estimated the suboptimality condition and suggested the guideline for the system containing a quadratic cost function.
In sub-optimal control problem, the previously continuous time horizon is discretized into N time horizons.Finding the solution of the controls in each time horizon helps archive the strategy which only applies the control at the constant level in each period.Caetano [25] and Yan [26] used the direct approach in finding the solutions to the sub-optimal disease control problem, and Yan [26] also used the Genetic Algorithm (GA).
Differential Evolution (DE) algorithm is an improved GA by Storn and Price [27].DE algorithm is simple and fast, requires only a few control parameters and has high convergence attribute [28].Thus, we find the solution of the sub-optimal control problem using the DE algorithm.The dynamical model is developed in Section 2 and the Differential Evolution is introduced in Section 3. The sub-optimal control problem is proposed through MPC and the numerical solutions are solved in Section 4. Lastly, the discussion and conclusion are presented in Section 5.

Dynamical Control Model for ZIKV Infection
To describe the behavior of ZIKV infection, the dynamical compartment model with control is developed from the models of Gao et al. [12] and Agusto et al. [13].The developed model consists of two human populations, adults and babies, and one mosquito (vector) population.The total newborn baby population consists of six classes: susceptible S b (t), exposed E b (t), asymptomatically infected A b (t), symptomatically infected without microcephaly, referred to as infected I b (t), infected with microcephaly I bm (t), and recovered R b (t).The total adult population consists of six classes: susceptible S w (t), exposed E w (t), asymptomatically infected A w (t), symptomatically infected, referred to as infected I w (t), infected with microcephaly I wm (t), and recovered R w (t).The total vector population consists of three classes: susceptible S v (t), exposed E v (t), and infected I v (t).
Hence, the total newborn baby population is N b (t) = S b (t) + E b (t) + A b (t) + I b (t) + I bm (t) + R b (t), the total adult population is N w (t) = S w (t) + E w (t) + A w (t) + I w (t) + I wm (t) + R w (t), the total vector population is N v (t) = S v (t) + E v (t) + I v (t), and the total human population is N h (t) = N b (t) + N w (t).The flow diagram of the compartment model is shown in Figure 1.Therefore, the ZIKV disease model with sexual transmission and vertical transmission is written by a system of differential equations under defined assumptions. where are the force of infection rates.
In the system, a susceptible adult can be infected either through the bite of an infected vector or sexual contact with an exposed or symptomatically infected adult.The infected class is more contagious due to the higher load of viremia [29], hence the infection modification parameters ρ w , ρ s are introduced for vector transmission and sexual transmission of an adult.
The newborn baby could be infected directly by mosquitoes or intrauterine via pregnant mother [30,31].We assumed that the exposure rate of babies is different from adults, represented by the modification parameter η > 0. The infection modification ρ b is introduced as well.In the vertical transmission in babies, it occurred in some cases that newborn baby is infected during pregnancy [31,32].We assumed that babies are born with infection due to vertical transmission by the parameter 0 ≤ q E , q I , q R ≤ 1.Let the newborn birth rate be π b , the fraction of healthy newborn babies are , and the infected babies are the remaining π b q E E w , π b q I Iw, and π b q R R w fractions.A recovered mother might give birth to the baby with microcephaly [30], and the baby born with microcephaly is assumed by the fraction (1 − r)π b q R R w in this paper.
Individuals with microcephaly experience a delay in development [33] and severe neurological disorders, so they are likely to have a short lifespan.Hence, we assumed that the adults with microcephaly do not reproduce or are involved in sexual transmission.The asymptomatically infected babies and adults are assumed to be noncontagious for both human and vector.Furthermore, the mortality rate from the disease is considered negligible.All parameter descriptions are summarized in Table 1.
Based on WHO [19] and Centers for Disease Control and Prevention (CDC) [20] suggestions, we considered vector, vertical, and sexual transmission the three control parameters selected for the model.The control variable u 1 (t) is the use of mosquito repellents to prevent or reduce the biting from mosquitoes.The control variable u 2 (t) is the procedure such as using protection for preventing pregnancy or unsafe sexual activity.The control variable u 3 (t) is the use of insecticide to control or eliminate the mosquito population.By applying control strategy, forces of infection in the adult and baby populations decrease by the fraction of (1 − u 1 (t)) and (1 − u 2 (t)).Correspondingly, the force of infection in vector population decreases by a fraction of (1 − u 1 (t)).Lastly, the recruitment rate of mosquito decreases by a fraction of (1 − u 3 (t)).It was also assumed that the mosquito death rate increases by r 0 u 3 (t), where r 0 > 0. Remark: If the value of each control parameters u 1 , u 2 , and u 3 is equal to zero for all given period, the model in the system in Equation ( 1) would be the Zika virus model without controls.
The controls u 1 , u 2 , and u 3 are determined to minimize the given objective function: subject to the state of system in Equation (1).In Section 4, we find a set of controls that minimize the number of exposed humans, symptomatically infected humans, all mosquitoes and the costs related to the controlling implementation.Firstly, we introduce constants A 1 , A 2 , A 3 , A 4 , and A 5 as the weighted constants related to the exposed newborn babies, symptomatically infected newborn babies, exposed adults, symptomatically infected adults, and all mosquitoes, respectively.The constants B 1 , B 2 , and B 3 are the weighted constants of the control variables u 1 , u 2 , and u 3 , respectively.The terms 1  2 B 1 u 2 1 , 1 2 B 2 u 2 2 , and 1 2 B 3 u 2 3 are the implementation costs of the three controls.The cost included in the first control might be the prices of using insect repellent, mosquito net, and herbal spray for instance.The cost included in the second control might be the prices of providing protections or warning leaflets about the safe sexual activity.The cost included in the last control could come from the expenses of using mosquito pesticides and the process implementation.

Differential Evolution
The differential evolution algorithm is a population-based optimization algorithm, introduced by Storn and Price [27,28].It is also one of the Evolutionary algorithms developed to find the parameter values and optimize real parameters or functions.DE is suitable for various practical problems which are nonlinear, non-differentiable, and multi-dimensional.The process consists of four parts, namely initialization, mutation, crossover, and selection.Suppose optimizing a problem consisting of D predicted parameters with N p number of population.Then, a D−dimensional vector is solved by the following main steps until the best solution vector reaches the termination criteria.
Step 1. Initialization The DE algorithm randomly selected a population of the parameter vectors {x G 1,i , x G 2,i , . . ., x G D,i }, i = 1, 2, . . ., N p } where G is the generation number.Define the initial population as , where x L j is the lower bound and x U j is the upper bound of the j−th parameter and φ i is a uniformly distributed random number between 0 and 1.Then, each of the N p parameter vectors will undergo mutation, crossover, and selection.
Step 2. Mutation This process helps expand the search space.For each target vector x G i , randomly select three vectors x G r1 , x G r2 , and x G r3 such that each index i, r1, r2, and r3 is different.Then, a mutant vector v G i is generated by where F ∈ [0, 2] is the mutation vector and v G i is called the donor vector.Step 3. Crossover Crossover or recombination process involves successful solutions from the previous generation.First, a trial vectors if rand j ≤ CR, or j = I rand x G i,j if rand j > CR and j = I rand where rand j ∼ U[0, 1] for a comparison to the crossover rate CR ∈ [0, 1], and I rand is a random index in {1, . . ., D} that u G i gets at least one component.

Step 4. Selection
The next generation vectors are selected by comparing the target vector x G i and trial vector u G i and selecting the one with the lowest objective function value.The objective function is given by comparing the value of E b , I b , E h , I h , and N v in the system in Equation ( 1) by applying vector u G i and x G i .Thus, the next generation vector x G+1 i is replaced by u G i or x G i under the following condition, where J is the objective function (see Section 4).
Step 5. Repeating After obtaining a new target vector, the mutation, crossover, and selection steps are repeated until the termination criteria are met.

Sub-Optimal Control Problem
Given the objective function as defined in Equation ( 2), we assume there exists a solution (u * 1 , u * 2 , u * 3 ) of the optimal control problem characterized by the Pontryagins Maximum Principle [21].
Since the implemented result of the optimal control problem is complicated for the real situation, a sub-optimal control problem is proposed.The idea of sub-optimal control problem is based on the works of Rodrigues et.al. [16] and Moulay et al. [17], and the Model Predictive Control [22].Regarding MPC, its strategy is based on the model of the system, a cost function which penalizes the undesirable behaviors, and constraints which represent the system limits.Furthermore, the stability and suboptimality of MPC can be estimated from the control horizon and the cost function [23,24].
For our nonlinear unconstrained model with a finite horizon, the continuous time horizon is discretized into N time horizon as Then, each of the controls u 1 (t), u 2 (t), and u 3 (t), in each time horizon, is approximated by a piece-wise constant control such that where u 1i , u 2i , and u 3i are constants and are referred to as the control parameters.Thus, for our model with now N control parameters for each u 1 , u 2 , and u 3 , we have total 3N parameters to be minimized.Let Then, the sub-optimal control problem is to subject to the system in Equation (1) and J is as defined in Equation (2).

Numerical Simulations
The simulations in this work are constructed using the parameter values in Table 2. Since the Zika virus is still under investigation and discovery, some parameter values are obtained from the previous study of other mosquito-borne viruses.The initial conditions are also selected for the theoretical sense and illustration purpose.The initial conditions for the system in Equation ( 1) are given by S b (0) = 10, 000, , and E 3 = 20.The total time span is 120 days, which is approximately four months.
From [24], the performance bound, which depends on the time horizon N and the cost function, could be computed to estimate the distance between MPC cost and the optimal cost.It also suggests that the time horizon N should not be too small.Hence, in our simulation with the total time of 120 days, we consider two scenarios with the different time horizon N as follows: (i) with N = 4; and (ii) with N = 8.All controls must be non-negative and let the maximum value be 0.9 since it is almost impossible to fully implement each control strategy.The numerical simulation from each scenario is generated with the DE algorithm with the maximum iteration of 150, F = 0.850, and CR = 1.000.Number of population N p is 120 and 240 for Scenarios N = 4 and N = 8, respectively.σ b , σ w 1 7.5 [36] γ b , γ w 1 8.5 [36] µ b 1 18.60×365 [3] µ w 1 70×365 [35] π v 500 [36] b 0.5 [35,37,38] σ v 1 10 [37] [38] r 0 0.1 [18] (i) Scenario N = 4 In this scenario, the period is 30 days or one month for each time horizon.Figure 2 represents the simulation of three controls over 120 days.From the graph, the levels of u 1 and u 2 are at the highest around 0.36 and 0.37 in the first month, respectively.After that, the levels of both controls u 1 and u 2 discretely drop in the second month and become approximately zero in the third month.Differently, the level of control u 3 stays at the maximum value 0.9 during the entire time.
The objective function value is 2.293810 × 10 5 , where the control values are:  (ii) Scenario N = 8 In this scenario, the period is 15 days or approximately two weeks for each time horizon.Figure 3 represents the simulation of three controls over 120 days.From the graph, u 1 is at the highest level around 0.6 in the first two weeks and discretely drops to approximately zero within one month.For u 2 , it reaches around 0.45 for the first two weeks, then drops to nearly zero after the second month.Differently, the level of control u 3 stays at the maximum value 0.9 during the entire time.
The objective function value is 2.293640 × 10 5 , where the control values are: The simulations in Figures 4-9 represent the number of exposed babies, infected babies, exposed adults, infected adults, exposed vectors and infected vectors, respectively.In each graph, the dotted line represents the model without control, the solid line represents the control model with Scenario N = 4, and the dashed line represents the control model with Scenario N = 8.In Figures 4-9, the numbers of exposed and infected individuals from the model with controls are less than those of the model without control.This implies that control strategy ceases the serve of the disease.The models from both scenarios yield very close results where N = 8 is slightly better.Notice in Figures 4-7 that the disease spreads heavily during the first month.Later, those patients decrease and become stable after two months, corresponding to the drop in control u 1 and u 2 after two months.Thus, both u 1 and u 2 contribute to the control of the disease.The control u 3 contributes the most and becomes the top priority in the strategy.

Conclusions
In this work, we developed the dynamical model of ZIKV disease with sexual transmission and vertical transmission.We proposed sub-optimal control problem with three control parameters: vector elimination, human-vector and human-human contact reduction.To find the control solution which is easy to implement, we partitioned the controls into discrete intervals.By using Differential Evolution, we solved and presented the numerical scenarios with four time horizons and eight time horizons.The simulations of both controls efficiently reduce the number of exposed and infected patients where the model with eight intervals yields a slightly better result.Controls u 1 , using repellents to prevent mosquito biting, and u 2 , avoiding pregnancy and unsafe sexual activity, should be applied mostly during the peak of the disease.After the disease dies out or becomes stable, the level of both controls will decrease.Mosquito elimination should be the first focus in controlling the disease, corresponding to the primary procedure announced by general health authorities.However, the cost associated with each control in this simulation is set to be equal.The consideration of adjusting the relative cost according to alternative scenarios might provide different results.Additionally, if more information is available in the future, the simulations should be conducted again using more precise data.

Figure 1 .
Figure 1.Flow diagram of the Zika virus model.

Figure
Figure Simulation results of infected vectors from: (i) the model without controls; (ii) the control model with N = 4; and (iii) the control model with N = 8.

Table 1 .
Definition of variables and parameters.Remaining fraction of symptomatic infection q E , q I , q R Fraction of newborn babies who are infected from pregnant adult of each class 1 − r Fraction of newborn babies who are infected with microcephalyβ b ,β w Transmission probability per contact from infected vectors to susceptible newborn babies and adults β v Transmission probability per contact from infected humans to susceptible vectors β s Transmission probability per sexual contact from infected humans to susceptible humans η Exposure modification parameter in babies ρ b , ρ w Infectivity modification parameters in exposed babies and adults ρ s Sexual infectivity modification parameters in exposed adults σ b , σ w Progression rate of exposed newborn babies and adults γ b , γ w Recovery rate of newborn babies and adults µ b , µ w