Mathematical analysis and optimal control of Dengue fever epidemic model

In this article, we are working on an SEIR-SI type model for dengue disease in order to better observe the dynamics of infection in human beings. We calculate the basic reproduction number R 0 and determine the equilibrium points. We then show the existence of global stability in each of the diﬀerent states depending on the value of R 0 . Moreover, to support the theoretical work, we present numerical simulations obtained using Python. We also study the sensitivity of the parameters included in the expression of R 0 with the aim of identifying the most inﬂuential parameters in the dynamics of dengue disease spread. Finally, we introduce two functions u and v, respectively indicating the treatment of the infected people and any prevention system minimizing contact between humans and the disease causing vectors. We present the curves of the controlled system after calculating the optimal pair of controls capable of reducing the dynamics of the disease spread, still using Python.


Introduction
Dengue is a viral and vector-borne disease.It is transmitted to humans by the bite of infected mosquitoes.Dengue is a fairly regular disease in tropical and subtropical regions, especially in urban and semi-urban areas.The vector causing the dengue disease are mosquitoes mainly of Aedes aegypti species.These mosquitoes are also vectors of chikungunya.
According to the World Health Organization, dengue fever is manifested by a wide range of symptoms that vary from basic to severe flu-like symptoms.Although rare, some people develop a severe form of dengue fever disease, which is characterized by various complications such as severe nose bleeding.The severe form of dengue, associated with high risk of death, was first discovered in 1950 during epidemics in the Philippines and Thailand.Nowadays, the severe form of dengue affects most countries in Africa, Asia, and Latin America, where it has become a major cause of hospitalization and death.There are four distinct but closely related stereotypes that cause dengue fever: DENV-1, DENV-2, DENV-3, and DENV-4.After recovery, lifelong immunity is obtained against the type of virus that caused the infection.A dengue patient who subsequently developed the disease is at high risk of developing a severe form of dengue.Dengue fever has existed in Burkina Faso since 1925.At the time of the last epidemic in 2016, the frequency of severe forms was 33% [20] in the city of Ouagadougou, the capital of Burkina Faso.
The mathematical modeling of phenomena is essential in the field of applied sciences.In physics, chemistry, biology, and many other fields, mathematical models are the pillars of a rigorous scientific study in order to make predictions about observed phenomena.Infectious diseases are increasingly becoming a major focus of mathematical modeling.malaria, dengue, Ebola virus disease, schistosomiasis, COVID-19 [6,26] are among many other infectious diseases that are regularly the subject of mathematical studies [7,8,17].Guiro et al. [16] analyzed a model of dengue transmission with general incidence.Ivorra et al. [18] worked on an optimal control model to reduce the spread of Ebola virus disease.In the modeling of schistosomiasis, Traore et al. [5] presented a study of a discrete class of schistosomiasis models with general delay and incidence.
There are numerous mathematical models of dengue fever [2,3,[14][15][16]21], most consider the group of infected people as one simple homogeneous group.However, the infected group is almost always made up of people with simple dengue disease and some people with a more severe case of the disease.Also, several optimal control studies on epidemiological models exist [1,7,8,17].Guiro et al. [1] developed a problem of optimal control of an SIR epidemic model with a general incidence function and time delays, while Kumar et al. proposed an optimal control problem of age-structured SEIRV model with imperfect vaccination [8].
As a contribution, we work on an SEIR-SI dengue model with an infected population made up of two groups: one group with simple infected cases and another one with hemorrhagic cases.The split of the infected people into two groups is necessary to better observe the dynamics of different types of dengue patients in the population.In addition to the stability study, we perform an optimal control study on an SEIR-SI dengue model.
In this work, we subdivide the human population into five compartments: susceptible people (S H ), exposed people (E H ), people infected with simple cases (I R H ), people infected with severe cases (I D H ), cured people (R H ), and deceased people (D).We present the dengue vectors in two classes, one of which designates the susceptible (S v ) and the other the infected (I v ) mosquitoes.
Our paper is organized as follows.In Sect.2, we present the system of ordinary differential equations after presenting the transfer diagram.In Sect.3, we make sure that our mathematical model is well defined.Indeed, we exhibit the mathematical properties of the model, calculate the equilibrium points and the basic reproduction number.We then study the stability of the equilibrium points in Sect. 4. We numerically simulate the model with some real data and some estimated data using Python 3.7 in Sect. 5. We also study the sensitivity of the parameters that appear in the expression of the basic reproduction number R 0 in Sect.6.An optimal control work is done in Sect.7, and we conclude in Sect.8.

Mathematical model
At a time t in each compartment i (i ∈ N), there are input movements E i and output movements S i of individuals.The input movements are counted positively and the outputs negatively.The dynamics in the compartment i at time t denoted by ˙ i is mathematically translated by the equation: With the aim to well observe the dynamics of humans infectious classes, we propose the following diagram.
In the compartment S H , for example, we have Then, according to Fig. 1, we obtain the following system of eight differential equations: with the initial conditions: • H is the recruitment of humans, which is assumed to be constant; • V is the recruitment of mosquito population, which is assumed to be constant; • β H is the average number of humans infected by the bites of infected mosquitoes per day and per Aedes mosquito (human.day - .vector - ).The flow of people from the class of susceptible humans to the exposed humans class is therefore β H S H N H I V ; • β V is the average number of Aedes mosquitoes infected by biting infected humans per day and per humans (vector.day - .human - ).The flow of mosquitoes from the class of susceptible mosquitoes to the infected mosquitoes class is therefore • μ H is the natural mortality rate of humans (persons.day - ).This means that mortality is not due to the disease.Thus a number μ H N H of people leaves the dynamics of disease propagation at each time t; • μ V is the mortality rate of vectors (mosquitoes.day - ).Thus a number μ V N V of mosquitoes leaves the dynamics of disease propagation at each time t.The infectious period of mosquito individual ends with their death; • η H denotes the transition rate (day -1 ) from state E H to I R H .We therefore consider a • θ is simple cases infected individuals proportion (day -1 ) who will be recovered (R H ); • 1θ is simple cases infected individuals proportion (day -1 ) who will be in severe cases • ρ is the proportion of severe cases infected individuals who will be recovered (R H ); • 1ρ is the proportion of severe cases individuals who will die (D H ); • γ H 1 is the transition rate (day 3 Properties of the mathematical model
We have Then it follows that according to Lemma 3.1, R 7 + is an invariant set for model (3).
Proposition 3.2 System (3) solutions are bounded in the region Proof We observe that Since t → +∞, we have then 0 solutions are bounded in .

Equilibrium points
In this subsection, we determine the equilibrium points of (3). Let By resolving the equations of (5), we get Let X 0 and X * be respectively the disease-free equilibrium (DFE) point and the endemic equilibrium point of model (3).At the disease-free equilibrium point, there are no infectious persons ( At the endemic equilibrium point X * , we have where

Basic reproduction number
Proof We use the next-generation matrix method [23] to calculate the reproduction number R 0 of model (2).Let F and V, the transmission and flow matrix between the infectious compartments E, I R H , I D H , and I V : On the disease-free equilibrium X 0 = ( H μ H , 0, 0, 0, 0, V μ V , 0), we obtain , and Then we get , and where , , .
The basic reproduction number R 0 [23] is defined as the dominant eigenvalue of the matrix -FV -1 .Therefore, where Any solution of the equation ϒ(I R H , I V ) = 0 in (0, S 0 H ) × (0, S 0 V ) corresponds to an equilibrium point with S H , I R H , S V , I V > 0. We remark that ϒ(0, 0) = 0 and ϒ(S 0 H , S 0 V ) < 0. The sufficient condition for the equation ϒ(I R H , I V ) = 0 to have a solution in (0, S 0 H ) × (0, S 0 V ) is the increasing of ϒ at the point (0, 0).This condition is equivalent to For all (I R H , I V ) ∈ R 2 + , we have: At the point (0, 0) we get: It follows that: Eq. ( 19) ⇐⇒ R 2 0 > 1.
Since R 0 > 1 by hypothesis, the conditions of Eq. ( 19) are therefore verified.These conditions implies the increasing of the function ϒ at the point (0, 0).The system (3) has then an unique endemic equilibrium point X * when R 0 > 1.

Global stability of equilibrium points
We study here the global stability of the disease-free equilibrium X 0 and the endemic equilibrium X * , when R 0 < 1 and R 0 > 1 respectively.

Global stability of disease-free equilibrium point X 0
Theorem 4.1 The disease-free equilibrium X 0 of system (3) is globally asymptotically stable when R 0 < 1.

Proof Let us consider the infected classes E H , I R H , I D H , and I V .
By the equations corresponding to these states, we have at X 0 the following system: System ( 20) can be rewritten as follows: We consider Ȳ = ( ĒH , ĪR H , ĪD H , ĪV ) and the following system: where We remark that F ≥ 0 and V is an asymptotic stable Metzler invertible matrix.
By applying Lakshmikantham standard comparison theorem in [9], we get We use ) is therefore globally asymptotically stable for R 0 < 1.

Global stability of endemic equilibrium point X *
Theorem 4. 2 The endemic equilibrium X * of system (3) is globally asymptotically stable when R 0 > 1.

Lemma 4.1 Let us consider the following function g defined by
We have that g(I R H ) ≥ 0. Indeed, by limited development we get Proof Consider the Lyapunov function candidate: Differentiating with respect to time yields ) is then a Lyapunov function.
In addition, we get V (X(t)) = 0 for According to the LaSalle invariance theorem [10,11], the endemic equilibrium point X * is globally asymptotically stable.

Numerical simulations and comments
In this section we have performed some numerical simulations to corroborate the theoretical work in the disease-free case and in the endemic case.In our simulations, we assume that the population is annually affected by dengue fever.This implies that the number of hemorrhagic cases in the event of dengue is considerable.
Depending on whether the number of basic reproduction is less than or greater than unity, the dynamics of the different states are observed using the Python version 3.7 programming software.
Comments In each of the subfigures (a)-(d) of Fig. 2 and (a)-(c) of Fig. 3, we have presented the evolution of the population of susceptible humans, exposed humans, simple infected humans, severely infected humans, recovered humans, susceptible vectors, and infected vectors respectively.Also, we present the evolution of susceptible persons, simple infected persons, and severe infected persons, all in the same subfigure Fig. 3(d).The evolution of infected vectors and infected persons, both severe and simple cases, are all shown in subfigure Fig. 3(e).
In a single subfigure, namely subfigure Fig. 3(g), we show the evolution of all human compartments.Similarly, we propose in the single curve Fig. 3(f ) the evolution of all compartments of vectors.
Theses curves point out that the disease tends to disappear in the population.Indeed, the human and vector infectious classes I R H , I D H , and I V are getting closer and closer to zero with the evolution of time t.This implies the disappearance of the disease in the population.The next curves are obtained with the estimated values grouped in Table 1.We keep the value of H , V , μ H , μ V , γ H 1 , γ H 2 and make the following changes: Indeed, with the estimated value rho = 0.3, we assume that at the endemic equilibrium point for ten (10) haemorrhagically infected people, three (03) will recover and seven (07) will die of dengue.Also, with the estimated value theta = 0.1, we assume that for every ten (10) people infected with simple dengue, one (01) will recover and nine (09) will experience a complication.Indeed, with the estimated value rho = 0.3, we assume that at the endemic equilibrium point for ten (10) haemorrhagically infected people, three (03) will recover and seven (07) will die of dengue.Also, with the estimated value theta = 0.1, we assume that for every ten (10) people infected with simple dengue, one (01) will recover and nine (09) will experience a complication.
Comments By the curves (a)-(e) of Fig. 4 and (a)-(f ) of Fig. 5, we show respectively the dynamics as a function of time (days) of the population of susceptible humans, exposed humans, simple infected humans, severely infected humans, recovered humans, susceptible vectors, and infected vectors respectively.
We show in the same figure, subfigure 5(c), the curve of the evolution of the simple infected and the severe infected.Also we show in the single Fig. 5(d) the three curves of the three infectious classes (simple infected humans, severe infected humans, infected mosquitoes).In 5(c) and 5(d), the persistence of the three infectious classes during the dynamics of the disease can be seen.
5(e) and 5(f ) show the dynamics of the different classes of vectors and the different classes of humans, respectively.
Overall, the different curves show a persistence of the disease in the population when the basic reproduction number is greater than 1.

Global sensitivity of R 0 parameters
The normalized forward sensitivity index of R 0 , which depends differentiably on the parameter ζ [13,16], is defined by Recall ( 18) that the basic reproduction number R 0 is given by For each parameter ζ of R 0 , thanks to (23), we evaluate the impact on R 0 of the parameters variation.By calculation, the sensitivity indices of R 0 with respect to β H , β V , μ H , μ V , η H , γ H 1 , and γ H 2 are as follows: ; Figure 5 The dynamics of humans and vectors, R 0 > 1 .
For the different parameters of R 0 , we give a value (Table 2) and obtain a numerical value of the sensitivity index (Table 2) corresponding to the parameter value.
The diagram (Fig. 6) is a graphic illustration of the sensitivity indices we have calculated.Indeed, by Fig. 6, we have in picture positive or negative impact of the variation of R 0 parameters.The fact that R 0 β H = +0.5 means that 1% increase in β H , keeping all other parameters fixed, will increase the value of R 0 by 0.5%.Also, when each of the parameters β V , η H , and γ H 1 increases by 1% while the other parameters remain constant, the value of R 0 increases by 0.5%, 0.285714285%, and 0.631883655%, respectively.On the other hand, when each of the parameters μ H , μ V , θ , and γ H 2 increases by 1% while the other parameters remain constant, the value of R 0 decreases by 0.532921918%, 0.5%, 0.00020842%, and 0.089233609%, respectively.

Statement of the optimal control problem
In this section, we compute the optimal function of the control (u(t), v(t)) to determine the best measures in terms of treatment and any prevention method to minimize the population of infected individuals.Indeed: • u consists of early supportive care with rehydration and symptomatic treatment.These methods improve the patient's survival; • v represents any method that can reduce vector-human contact (destruction of egg nests, mosquito nets, awareness-raising, etc.).This optimal couple minimizes at the same time the cost of implementing the treatment and the prevention strategies.So we consider the following optimal control problem: where subject to the equation The two functions u(t) and v(t) represent respectively the treatment of infected persons and the different ways to prevent dengue fever.The use of insecticides, destruction of egg nests, and any other measure that reduces mosquito bites are ranked in the order of prevention methods.The functions u and v are elements of the set U defined by The two constants A 1 > 0, A 2 > 0 are measures of the relative cost of the interventions associated with the controls u and v respectively.
Theorem 7.1 (Existence of optimal control) Consider the optimal control problem (24) subject to (26).Then there exists an optimal control ( ū, v) in and a corresponding solution that minimize J(u, v) over a set of admissible controls .
Proof The existence of an optimal control can be obtained by using a result of Flemning and Rishel [4].To use an existence result, Theorem III.4.1 from Lukes [12], we must check if the following properties are satisfied: (1) The set of controls and the corresponding solutions are not empty.
(2) The set of admissible controls U is convex and closed in L 2 (0, T).
(3) The vectors field of the state system is borned by a linear control function.

Proposition 7.1 (Hamiltonian characterization of minimization problem)
The minimization problem (24) induces to a problem of minimization of Hamiltonian H defined by where: • f i is the right-hand side of the differential equation of ith state variable, • p(•) is absolutely continuous application defined to [0, Proof Let where and Then the Hamiltonian of the optimal problem is defined by and verify Proof According to Theorem (7.1) the couple of controls ( ū, v) associated with the solution X minimize J(u, v) over U.According to the maximum principle of Pontryagin, there exists an absolutely continuous application Then Then we have By the same method, we have The condition of transversality at final time t f is p(t f ) = 0. Then Finally, the characteristics of the vector p(•) : t − → λ 1 (t), λ 2 (t), λ 3 (t), λ 4 (t), λ 5 (t), λ 6 (t), λ 7 (t), λ 8 (t) are (31) Theorem 7.2 (Characterization of optimal control) The optimal control ( ū, v) is defined by where Proof To prove the characterizations of optimal control, we define the Lagrangian associated with the problem.It corresponds to Hamiltonian increased by coefficients of penalty.
where w ij (t) ≥ 0 are penalization coefficients that verify The partial derivative of H in relation to u is given by The partial derivative of H in relation to v is given by We obtain and at ū and v, we have: Let be the set {t : 0 < ū < 1}.We have Let be the set {t : ū = 0}.We have w 12 (1 -ū) = 0 ⇒ w 12 = 0, therefore Since w 11 ≥ 0, then -M( X(t)) Thus, on the set {t : 0 ≤ ū < 1}, ū is defined as follows: Let be the set {t : ū = 1}.We have Since w 12 ≥ 0, then -M( X(t)) On the set {t : 0 ≤ ū ≤ 1}, ū is defined by By the same method, we get the expression of v: Finally, on the set U, the optimal control ( ū, v) is given by where

Numerical simulations
In this section, the optimal system and the results of simulations obtained by Python 3.7 (see the Annex) are presented.The optimal system is given by With the values of Table 1, we obtain the following curves:  Figures 7(c), 7(d), and 7(h) represent the different dynamics of infected vectors and infected people.The orange color curves represent the dynamics when there are treatment and prevention (u = 0 and v = 0) and the blue curve represents the population when there are no vaccination and treatment (u = 0 and v = 0).This shows that without prevention and treatment infected individuals will be more numerous.
Again using the data in Table 1, we present the curves expressing the dynamics in the different infectious classes.
Curves (a), (b), and (c) in Fig. 8 show the population dynamics of infected humans in single cases, infected humans in severe cases, and infected vectors (mosquitoes) as a function of time, respectively.
In each of the three subfigures, it can be seen that when prevention alone is applied, the trend is the same as when both prevention and treatment are applied.In fact, when prevention alone is applied, the infectious classes tend towards zero after a certain time.When treatment is added, the number of individuals in the infectious classes continues to tend towards zero, but after a shorter time.

Conclusion
In our paper, we worked on an SEIR-SI model of dengue fever in a population regularly affected by this disease.There are considerable hemorrhagic cases of the dengue disease because of the regularity of the disease.We took full account of the hemorrhagic cases of dengue denoted by I D H .We have shown the global stability of the equilibrium points depending on whether R 0 is greater or less than 1.Moreover, to corroborate the theoretical study, we presented numerical simulations using Python.We then calculated the sensitivity of each parameter affecting the basic reproduction number R 0 .We found that a variation of +1% in each of the parameters β H , β V , η H , γ H 1 separately increases the value of R 0 and a variation of +1% in each of the parameters μ H , μ V , θ , γ H 2 decreases the value of R 0 .Also, γ H 1 and μ H when varied by +1% are respectively the parameters that have the most positive and negative influence on the value of R 0 .Finally, optimal control work was carried out for finding how the optimal pair (u(t); v(t)) is able to drastically reduce the spread of dengue in a population regularly affected by this disease.

Annex
Numerical method with python 3.7 In first, we implement (2), the model without control by using the function odeint of PYTHON.We obtain for example def F0(Y , t) : where Y0 is the initial condition and T is the time.
Secondly, we implement model (32) by using the method of shoot [22].Let The solution of (34) depends on T and y 0 , and is written y(T, y 0 ).At final time T, y T, y 0 = y(T), (35) and this means y T, y 0y(T) = 0. (36)

Figure 1
Figure 1 Transfer diagram

Figure 6
Figure 6 Sensitivity indices diagram

Figure 7
Figure 7 Humans and vectors dynamics with and without control (u(t);v(t))

Figure 8
Figure 8 Dynamics of infectious classes with v(t) and with (u(t);v(t)) = odeint(F0, Y 0, T), S H , E H , I R H , I D H , R H , S V , I V , λ 1 , . . ., λ 8 ) y = (y 1 , y 2 , y 3 , y 4 , y 5 , y 6 , y 7 , y 8 , y 9 , y 10 , . . ., y 15 , y 16 ).(33)By re-writing model (32), we get the two point boundary value problem ) = F(t, y(t)),y(0) = (S 0 H , E 0 H , I 0 R H , I 0 D H , R 0 H D 0 H , S 0 V , I 0 V , λ 1 (0), .. ., λ 8 (0))) = y 0 , y(T) = (S H (T), E H (T), . . ., D H (T), S V (T), I V (T), 0, . . ., 0) = y T .(34) where • S H , E H , I R H , I D H , I D H , and R H denote the number of susceptible (S H ), latent (E H ), infectious (I R H , I D H ), and recovered (R H ) persons at time t.I R H and I D H designate respectively people with classical dengue fever and those with dengue hemorrhagic fever; • S V and I V denote the number of vectors (mosquitoes) in susceptible and infectious states at time t, respectively; • N H = S H + E H + I R H + I D H + R H and N V = S V + I Vare the human and vector population size at time t, respectively; -1 ) from state I R H to R H with proportion θ and to I D H with proportion 1θ .So θγ H 1 I R H leaves I R H for R H and (1θ )γ H 1 I R H leaves I R H for I D H ; • γ H 2 is the transition rate (day -1 ) from state I D H to R H with proportion ρ and to D H with proportion 1ρ.So ργ H 2 I D H leaves I D H for R H and (1ρ)γ H 2 I D H leaves I D H for D H .