Mathematical study for Zika virus transmission with general incidence rate

1.   Introduction

Zika fever is caused by an arbovirus (virus transmitted by insects), belonging to the Flaviviridae family, of the flavivirus genus, such as dengue viruses or yellow fever ^[1,2,3,4]. The species currently capable of transmitting Zika virus is Aedes aegypti, which originates from Africa ^[5,6]. The Aedes albopictus (tiger mosquito, native to Asia) could also be found to be a vector of the Zika virus, as it is already for dengue and Chikungunya ^[7,8]. The mosquito becomes infected with the virus during a blood meal, when it bites an infected person. The virus multiplies within the mosquito without affecting the insect. Then, the next time it is bitten, the mosquito releases the virus into the blood of a new person. Symptoms appear 3 to 12 days after the bite, but during this time a person can infect other mosquitoes if they get bitten again. Therefore, patients with Zika should avoid being bitten in order to interrupt the cycle of viral transmission. The majority of people infected with the virus (an estimated 70 to 80% of cases) do not develop any symptoms. In the rest of the population, the symptoms caused by the Zika virus are flu-like: fatigue, fever (not necessarily high), headaches, muscle and joint pain in the limbs. In addition to these symptoms, there are different types of rashes. Conjunctivitis, pain behind the eyes, digestive problems or even edema of the hands or feet may appear. In most cases, the symptoms are moderate and do not require hospitalization. These symptoms are not very specific, and the Zika virus being found in the same regions as those of dengue and Chikungunya, make an exact diagnosis difficult.There is no specific antiviral treatment, nor an active vaccine against the Zika virus. Only symptomatic treatment can be prescribed (analgesics against pain and fever of the paracetamol type).

Zika virus was first detected in a monkey in Uganda in 1947. A year later, it was isolated in the same region from an Aedes mosquito. The first human cases appeared in the 1970s in other African countries (Uganda, Tanzania, Egypt, Central African Republic, Sierra Leone, Gabon and Senegal), then, in certain Asian countries (India, Malaysia, Philippines, Thailand, Vietnam and Indonesia) ^{[1,2,5,6,9,10,11]}. In 2007, a real epidemic broke out in Micronesia (Yap Island in the Pacific), causing 5,000 infections ^[12,13,14]. In 2013 and 2014, in French Polynesia, 55,000 cases of Zika were reported. The epidemic then spread to other islands in the Pacific, including New Caledonia, the Cook Islands and Easter Island ^[3,4]. The Zika virus was first detected in north-western Brazil in May 2015 and its presence is spreading very quickly to other regions of the country. Brazil reports the highest number of Zika cases ever described: between 440,000 to 1,500,000 reported suspected cases ^[15,16]. The virus has been present since October 2015 in Colombia, Salvador, Guatemala, Mexico, Panama, Paraguay, Surinam, Venezuela and Honduras ^[17]. In 2015, some cases were detected in French Guyana and in Martinique ^[17]. In 2016, Martinique had 16,650 suggestive cases (in the process of biological confirmation). Guyana has 3,620 cases and Guadeloupe 1,090 cases. In France, 176 cases have been biologically confirmed in people returning from the Zika virus circulation area, including 7 pregnant women and 1 case of neurological complications. A person has been infected with the Zika virus through sexual intercourse. The spread of Zika fever could take place in areas where the Aedes mosquito is already implanted and where a person already infected with Zika is staying. In mainland France, the Aedes albopictus mosquito (tiger mosquito) is present in 30 departments. The vectorial expansion period usually occurs in May and its period of activity (and therefore the risk of transmission of the virus) is between May and November.

The worst epidemics and pandemics that have ravaged humanity brought out the role of mathematical models in political and health decision-making. In fact, mathematical models are essential tools investigating the dynamics of the spread of an epidemic that help us to determine a threshold called reproduction number which usually provides information on how infection will be sustained. The potential factors for such a virus spread can be well estimated via mathematical models. Infectious disease control is key to the World Health Organization and optimal time control is able to provide useful theoretical information about both the prevention and treatment of diseases. The information obtained from optimal control modeling can help decision makers to plan well and provide the best services for the management and treatment of diseases. The optimal time control has been used to study several mosquito-related diseases ^{[18,19,20,21,22,23,24,25]}. To our knowledge, there is no optimal control model for ZIKA virus. Several models have been proposed to describe the transmission dynamics vector-borne diseases see ^{[26,27,28,29,30,31,32,33,34]} for example. The case of Dengue, has is the subject of several studies ^{[12,13,35,36]} for example. In ^{[12,13,14,35,36]}, the authors propose the study of SI-SIR models to describe the transmission dynamics of Dengue fever in which the total human population is assumed to be constant. In ^[13], the author does not use a standard infection rate in transmitting the infection towards the human population, which this time depends on the size of the vector. In all these models, the growth dynamics of the vector, and especially immature stages, is not taken into account. In ^[37], describing the transmission of Dengue, the dynamics of the immature stage is taken into account, and the egg stage is regulated by the effect of a reception capacity of the medium. This model also considers an infection rate dependent on the vector population. This article only deals with the global stability of the endemic equilibrium. In ^[29], the authors proposed a model describing this time, the transmission of Chikungunya with a classic infection rate, or the growth dynamics of an immature stage, grouping together the eggs, larvae and nymphs but in which only local stability is studied.

The bilinear incidence rate $\beta SI$ and the standard incidence rate $\beta (S/N)I$ for the classical epidemiological models are often considered ^[29,30,36]. However, for collected data from the real phenomenon, the disease dynamics is not always follows these type of rates and many of the epidemiological mechanisms are more appropriate with nonlinear transmission rates, in particular the saturated rates of the form incidence rate $g(I)S$ ^[23,24] or of more general form $g(S, I)$ ^[38]. When a large number of infective involves in the population then exposure to the disease agent is virtually certain and the transmission rate may slow down. This happened because of the number of effective contact between the susceptible population and infected population may saturate at high infection level due to overcrowding of infective or due to protective measures by susceptible. Since a large number of infective involves in vector borne disease, saturated incidence rate has been considered as more suitable for the vector borne diseases like dengue than linear incidence rates. Hence, the goal of this paper is to generalize the results of ^[26,27,28] to nonlinear incidence rates. In ^[23,24], the authors have used similar saturated incidence rate in order to describe vector borne diseases.

Treatment is significant in every infectious disease for the infected population to become recovered. Usually, the treatment rate is considered to be proportional to the number of infective individuals (linear dependence on the number of infected individuals). In this paper we propose a more general mathematical model for ZIKA virus transmission. We analysed the model and we studied the sensitivity analysis with respect to the input parameters. Then we used optimal control to examine the effect of treatment. The purpose of this control is to minimize the number of infected individuals with optimal treatment cost by applying Pontryagin's maximum principle.

This paper is organised as follows. The mathematical model, its general properties, the calculation of the basic reproduction number using next generation method, the existence and the number of equilibrium points are presented in section 2. In section 3, we discussed the local stability of the steady states using the Jacobian linearisation method. Later in section 4, we discussed the global stability of the equilibrium points using the basic Lyapunov theory. In section 5, the sensitivity analysis of the basic reproduction number with respect to each of the model parameters was carried out. In section 6, we presented the optimal strategy to minimize the number of infected individuals and we formulated the optimal control problem then we expressed the control in terms of the adjoint variables. Numerical simulations have been given in section 7. Finally, some conclusion remarks have been presented in section 8.

2.   Nonlinear Zika virus dynamics model

The virus is spread as follows. Mosquitoes contract the virus by biting animals or men infected with it. They then spread the virus by biting uninfected people. It has been observed that mosquitoes, once infected, remain so until they die. Conversely, men become immune after contracting the disease. We will therefore use an SI model for describe the transmission of the virus in the vector population and a SIR-type model to describe it within the human population ^[39]. We then obtain the generalization of the models given in ^[26,27,28]).

$S_h, I_h$ and $R_h$ describe the susceptible, infected and recovered in the human population, respectively. $S_v$ and $I_v$ describe, respectively, the susceptible and infected in the mosquito population.

In the table hereafter, we give a description of the model parameters.

In the human population, the uninfected become infected at a rate $f_h(I_h)S_h$, where $f_h$ is the infected to uninfected incidence rates.

Assumption 1. $f_h$ is non-negative $C^1(\mathbb{R}_+)$, increasing concave function such that $f_h(0) = 0$.

Lemma 1. $If_h'(I) \leq f_h(I) \leq If_h'(0), \; \forall I \in \mathbb{R}_+$.

Proof. For $I, I_1 \in \mathbb{R}_+$, let $g_1(I) = f_h(I)-If_h'(I)$. Since $f_h$ is an increasing concave function then $f_h'(I) \geq 0$ and $f_h''(I) \leq 0$. Therefore $g_1'(I) = -If_h''(I) > 0$ and $g_1(I) \geq g_1(0) = 0$ or also $f_h(I)\geq If_h'(I)$. Samelly, let $g_2(I) = f_h(I)-If_h'(0)$ then $g_2'(I) = f_h'(I)-f_h'(0) < 0$ since $f_h$ is concave. Then $g_2(I) \leq g_2(0) = 0$ then $f_h(I)\leq If_h'(0)$.

2.1.   Basic properties

Lemma 2.

is a positively invariant attractor of all solutions of system (2.1) where $N_h = S_h+I_h+R_h$ and $N_v = S_v+I_v$ are the total population sizes of human and mosquitoes, respectively.

Proof. We have

Thus one can deduce that the closed non-negative cone $\mathbb{R}_+^5$ is positively invariant by system (2.1). From Eqs (2.1) we get $\dot N_h = \Lambda_h -\mu_h N_h$. Hence $N_h = \dfrac{\Lambda_h}{\mu_h}$ if $N_h(0) = \dfrac{\Lambda_h}{\mu_h}$. Similarly, $\dot N_v = \Lambda_v -\mu_v N_v$. Hence $N_v = \dfrac{\Lambda_v}{\mu_v}$ if $N_v(0) = \dfrac{\Lambda_v}{\mu_v}$.

2.2.   Steady states: Existence and uniqueness

For any disease, a major health issue is whether it is spreading in the population and at what speed (doubling time). This amounts to calculating the average number of people an infectious person infects while they are contagious. This rate is called the basic reproduction number, and is denoted $\mathcal{R}_{0}$ (ratio). This rate is intuitively easy to understand, but if it is linked to the pathogen, its calculation is complex. $\mathcal{R}_{0}$ should be used with caution, as it can lead to misinterpretations, both on the real role that $\mathcal{R}_{0}$ has on the spread of an infectious disease and on the ability to control an epidemic. The calculation of the $\mathcal{R}_{0}$ presupposes a population where all the individuals are healthy, except the infectious individual introduced (patient zero). If $\mathcal{R}_{0} < 1$, then the infected individual infects less than one other individual on average, which means that the disease is disappearing from the population. If $\mathcal{R}_{0} > 1$, then the disease spreads in the population and becomes epidemic. Determining $\mathcal{R}_{0}$ according to the parameters of the model thus makes it possible to calculate the conditions under which the disease spreads. As van den Driessche and Watmough note ^[40], "in the case of a single infected compartment, $\mathcal{R}_{0}$ is simply the product of the infection rate and its mean duration". When the model is simple, it is often possible to find an exact expression for $\mathcal{R}_{0}$. If there are several compartments representing infectious individuals, other methods are necessary, such as the next-generation matrix method introduced by Diekmann ^[41]. In the approach of van den Driessche and Watmough for this method ^[40], the population is divided into $n$ compartments whose $m$ first represent the infected individuals. In other words, instead of focusing on a single compartment of infected individuals as previously, the method considers that the infected individuals are distributed over $m$ compartments. The goal is therefore to see the rate of change in the population established in each of these compartments.

The matrix approach to calculate the basic reproduction number in complex models shows the relationship between compartmental models and population matrix models. In our case, $F = \left(\begin{array}{c} 0&f_h'(0)\dfrac{\Lambda_h}{\mu_h}\\ f_v'(0) \dfrac{\Lambda_v}{\mu_v} & 0 \end{array} \right)$ and $V = \left(\begin{array}{c} r_h+u+\mu_h & 0 \\ 0 & \mu_v \end{array} \right)$.

The determinent of $V$ is given by det$(V) = \mu_v(r_h+u+\mu_h) > 0$ and therefore $V^{-1} = \left(\begin{array}{c} \dfrac{1}{r_h+u+\mu_h} & 0 \\ 0 & \dfrac{1}{\mu_v} \end{array} \right)$ and the next generation matrix is given by

Then the basic reproduction number for model (2.1) is given by the spectral radius of the matrix $FV^{-1}$:

Lemma 3.  ● If $\mathcal{R}_0\leq 1,$ then model (2.1) has a trivial steady state$E_{0}$.

● If $\mathcal{R}_0 > 1,$ then model (2.1) has two steady states $E_{0}$ and $E^*$.

Proof. Let $E(S_h, I_h, R_h, S_v, I_v)$ be any steady state satisfying

which is equivalent to

From the second equation of system (2.2), we have

If $I_h = 0$, then we get a steady state given by the ZIKV-free steady state

$E_{0} = \left(\dfrac{\Lambda_h}{\mu_h}, 0, 0, \dfrac{\Lambda_v}{\mu_v}, 0\right)$. If $I_h\not = 0$, define the function

Then we obtain

Now, we have

The derivative of the function $g$ on $(0, \dfrac{\Lambda_h}{\mu_h})$ is given by

Thus $g$ is an decreasing function. Then the equation $g(I_h) = 0$ admits a unique solution $I_h^*\in(0, \dfrac{\Lambda_h}{\mu_h})$. Therefore,

and the infected steady state $E^* = (S_h^*, I_h^*, R_h^*, S_v^*, I_v^*)$ exists only if $\mathcal{R}_0 > 1$.

3.   Local stability of steady states

The $\mathcal{R}_0$ rate is important to determine the risk of a new pathogen causing an epidemic (impossible if $\mathcal{R}_0 < 1$, possible if $\mathcal{R}_0 > 1$) and to estimate the final size of the epidemic, with or without control measures.

Theorem 1. If $\mathcal{R}_0 < 1$, then the trivial equilibrium point $E_0$ is locally asymptotically stable.

Proof. The Jacobian matrix of system (2.1) at a point $(S_h, I_h, R_h, S_v, I_v)$ is given by:

The Jacobian matrix evaluated at $E_0$ is then given by:

$J_0$ admits five eigenvalues. The first three eigenvalues are given by $\lambda_1 = \lambda_2 = -\mu_h < 0$ and $\lambda_3 = -\mu_v < 0$. The other two eigenvalues are those of the following sub-matrix

where the trace is given by

and the determinant is given by

Then all eigenvalues must have negative real parts once $\mathcal{R}_0 < 1$ and the trivial equilibrium point $E_0$ is then locally asymptotically stable once $\mathcal{R}_0 < 1$.

Theorem 2. If $\mathcal{R}_0 > 1$, then the endemic equilibrium point $E^* = (S_h^*, I_h^*, R_h^*, S_v^*, I_v^*)$ is locally asymptotically stable.

Proof. The Jacobian matrix at the endemic equilibrium $E^* = (S_h^*, I_h^*, R_h^*, S_v^*, I_v^*)$ is given by:

The characteristic polynomial is then given by:

Then $\lambda_1 = -\mu_v < 0$ and $\lambda_2 = -\mu_h < 0$ are two eigenvalues. The other three eigenvalues are the roots of

where

Using the fact that $f_h'(I_v^*) \leq \dfrac{f_h(I_v^*)}{I_v^*}, f_v'(I_h^*) \leq \dfrac{f_v(I_h^*)}{I_h^*}, (r_h+u+\mu_h) = \dfrac{f_h(I_v^*)S_h^*}{I_h^*}$, and $\mu_v = \dfrac{f_v(I_h^*) S_v^*}{I_v^*}$, we obtains

Then the eigenvalues have negative real parts. Therefore, the equilibrium $E^*$ (exists if $\mathcal{R}_0 > 1$) is locally asymptotically stable. The proof is completed.

4.   Global properties of steady states

Theorem 3. By considering the dynamics (2.1), if $\mathcal{R}_0 \leq1,$ then $E_{0}$ is globally asymptotically stable.

Proof. Assume that $\mathcal{R}_0\leq1$ and define the following Lyapunov function $U_{0}(S_h, I_h, R_h, S_v, I_v)$:

Clearly, $U_{0}(S_h, I_h, R_h, S_v, I_v) > 0$ for all $S_h, I_h, R_h, S_v, I_v > 0$ and $U_{0}\left(\dfrac{\Lambda_h}{\mu_h}, 0, 0, \dfrac{\Lambda_v}{\mu_v}\right) = 0$. The derivative of $U_{0}$ with respect to time along system (2.1)is given by:

If $\mathcal{R}_{0}\leq1,$ then $\dfrac{dU_{0}}{dt}\leq0$ for all $S_h, I_h, R_h, S_v, I_v > 0.$ Let $W_{0} = \left\{(S_h, I_h, R_h, S_v, I_v):\dfrac{dU_{0}}{dt} = 0\right\}$. It can be easily shown that $W_{0} = \left\{E_{0}\right\}$. Applying LaSalle's invariance principle ^[42] (see ^[43,44,45] for other applications), we deduce that $E_{0}$ is GAS when $\mathcal{R}_0\leq 1.$

Define the set

Theorem 4. By considering the dynamics (2.1), if $\mathcal{R}_0 > 1,$ then the endemic equilibrium $E^* = (S_h^*, I_h^*, R_h^*, S_v^*, I_v^*)$ is globally asymptotically stable in $\Omega$.

Proof. Define a function $G(x) = x-1-\ln x$. Let a function $U^*(S_h, I_h, R_h, S_v, I_v)$ be defined as:

Clearly, $U^*(S_h, I_h, R_h, S_v, I_v) > 0$ for all $S_h, I_h, R_h, S_v, I_v > 0$ and $U^*(S_h^*, I_h^*, R_h^*, S_v^*, I_v^*) = 0.$ Calculating $\dfrac{dU^*}{dt}$ along the trajectories of (2.1), we obtain

Applying the steady state conditions for $E^*$

we get

Therefore, $\dfrac{dU^*}{dt}\leq 0$ for all $S_h, I_h, R_h, S_v, I_v \in \Omega$ and $\dfrac{dU^*}{dt} = 0$ if and only if $(S_h, I_h, R_h, S_v, I_v) = (S_h^*, I_h^*, R_h^*, S_v^*, I_v^*) = 0$. By applying the LaSalle's invariance principle ^[42], one can easily deduce that $E^*$ is globally stable (see ^[23,24,25] for other applications).

5.   Sensitivity analysis

Sensitivity analysis consists of identifying and quantifying the contribution of the input parameters of the model to illness transmission ^[21,46]. This type of analysis is particularly considered in the context of a reliability study, the overall methodology for dealing with uncertainties or even in a robust model. For the engineer responsible for reducing the variability of a quantity of interest, the sensitivity analysis makes it possible to distinguish the modeling parameters that need to be better controlled. For the numericalist, it allows a reduction in the numerical cost by a reduction in the number of parameters, and therefore of the size of the model, by eliminating those that do not contribute to the variability of the response. It is used regularly to determine the robustness of model predictions for parameter values since the data collection can contain noises^[21]. The variable sensitivity with respect to parameters is generally measured by a sensitivity index which allow us to determine the relative change in a variable when once the parameter value changes. Thus, the importance of each parameter to the disease transmission is fixed by the sensitivity analysis. In particular, we can discover parameters that have a highest impact on the basic reproduction number $\mathcal{R}_0$. Since $\mathcal{R}_0$ is a differentiable function to some parameters, the sensitivity index is calculated using partial derivatives. If the index is negative, then the relationship is inversely proportional. However, if the sensitivity index is positive then an increase in the value of a parameter imply an increase in the value of $\mathcal{R}_0$ ^[46].

Definition 1. (see ^[21,47]). The sensitivity index of $\mathcal{R}_0$ with respect to a parameter $\rho$ is defined by

Proposition 1. The explicit expression of $\mathcal{R}_0$ is given by

The sensitivity indices of $\mathcal{R}_0$ with respect to the model parameters are given by Table 1:

Only the sensitivity indices of $\mathcal{R}_0$ with respect to $\Lambda_h$ and $\Lambda_v$ are positive. The other sensitivity indices are negative. Note that the indices $Y^{\mathcal{R}_0}_{\mu_h}, Y^{\mathcal{R}_0}_{r_h}, Y^{\mathcal{R}_0}_u$ and $Y^{\mathcal{R}_0}_{}$ are not constant but depend on other parameters. Note also that $\mathcal{R}_0$ is most sensitive to the recruitment rate of human population $\Lambda_h$ and the recruitment rate of mosquito population $\Lambda_v$. This means that an increase in one of the recruitment rates increases the disease transmission in the human population.

From the Table 1, it is to be noted that the sensitivity indices of $\mathcal{R}_0$ with respect to the parameters $\Lambda_h$ and $\Lambda_v$ are positive and hence an increase of the values of parameters $\Lambda_h$ and $\Lambda_v$ leads to an increase of the value of $\mathcal{R}_0$. The sensitivity indices of $\mathcal{R}_0$ with respect of the remaining parameters $\mu_v, \mu_h, r_h$ and $u$ are negative and hence an increase of the values of parameters $\mu_v, \mu_h, r_h$ and $u$ leads to a decrease of the value of $\mathcal{R}_0$.

The sensitivity analysis can be confirmed by the Figure 1 where we presented the behaviour of the basic reproduction number, $\mathcal{R}_0$, with respect to the model parameters. $\mathcal{R}_0$ increases only with the parameters $\Lambda_h$ and $\Lambda_v$, however it decreases with the remaining parameters $\mu_v, \mu_h, r_h$ and $u$.

6.   Optimal strategy via treatment rate

Intuitively, one wonders how to act on transmission to control the epidemic. The theory of optimal control (or optimal control) goes further by allowing this intensity of control to vary over time. We therefore seek the best control value to use. In our case, we consider a time-varying treatment function ($u(t)$).

Assume moreover that $f_h$ and $f_v$ are globally Lipschitz with an upper bound $\bar f_h = \sup_{I > 0} f_h(I)$ and $\bar f_v = \sup_{I > 0} f_v(I)$, respectively and with Lipschitz constants $L_h$ and $L_v$, respectively. The control set $\bf{P_{ad}}$ is

The aim is to find the best control $u(t)$ and associated state variables $S_h(t$), $I_h(t)$, $R_h(t)$, $S_v(t)$ and $I_v(t)$ to apply on the epidemic to minimize an objective functional:

where $\alpha_1, \alpha_2$ and $\alpha_3$ are positive balancing constants. The goal is to optimise the number of infected population, and maximize the number of recovered population while minimizing the cost of the treatment control. If the strategy is to minimize the infected population and not be concerned with the cost, one can take $\alpha_3 = 0$.

Existence and uniqueness

Since the system (2.1) is linear in the control $u(t)$ with bounded states, we can use standard results ^[48] to show the existence of the optimal control.

For $\varphi = (S_h, I_h, R_h, S_v, I_v)^t$, the model (2.1) can be modified in the simple form

where $A = \left(\begin{array}{cccccc} -\mu_h & 0 & 0 & 0 & 0\\ 0 & -(r_h+u+\mu_h+) & 0 & 0 & 0\\ 0 & (r_h +u) & -\mu_h & 0 & 0 \\ 0 & 0 & 0 & -\mu_v & 0\\ 0 & 0 & 0 & 0 & -\mu_v \end{array} \right)$ and $F(\varphi) = \left(\begin{array}{ccc} \Lambda_h -f_h(I_v)S_h\\ f_h(I_v)S_h\\ 0 \\ \Lambda_v -f_v(I_h) S_v\\ f_v(I_h) S_v \end{array} \right)$.

Proposition 2. The function $K$ is continuous and uniformly Lipschitz.

Proof. We can easily prove that the function $F$ is continuous and uniformly Lipschitz as follows

where $M = 2 \max(\bar f_h, \; \dfrac{L_h \Lambda_h}{\mu_h}, \; \bar f_v, \; \dfrac{L_v \Lambda_v}{\mu_v})$. Since,

where $\big\|A\big\|_1: = \sup_{X\neq 0}\frac{\big\|AX\big\|_1}{\big\|X\big\|_1}$ is the norm of the matrix $A$ subordinate to the vector norm $\big\|\cdot\big\|_1$. Therefore

Here $K = \max(M, \big\|A\big\|)$. Then the continuous functions $F$ and $G$ are uniformly Lipschitz. Therefore the model (6.1) admits a unique solution.

Now, we will derive some necessary conditions that permits to calculate optimal control and corresponding states using Pontryagin's Maximum Principle ^[48,49,50]. The Hamiltonian is

For a given optimal control $u^*$, there exist adjoint functions $\lambda_1, \lambda_2, \lambda_3, \lambda_4$ and $\lambda_5$ corresponding to the states $S_h, I_h, R_h, S_v$ and $I_v$ such that:

where the final conditions are given by $\lambda_1(T) = 0$, $\lambda_2(T) = 0$, $\lambda_3(T) = 0$, $\lambda_4(T) = 0$ and $\lambda_5(T) = 0$.

Note that the partial derivative of the Hamiltonian with respect to the control variable at $u$ is given by

$\dfrac{\partial H}{\partial u}$ is linear with respect to the control $u$ and $\dfrac{\partial H}{\partial u} = 0$ for

if

Thus, the strategy is to choose the control $u$ as following

7.   Numerical simulations

In this section, we adopted the numerical simulations validating analytical findings. We consider incidence rates of Michaelis-Menten type: $f_h(I) = \dfrac{a_h I}{b_h+I}$ and $f_v(I) = \dfrac{a_vI}{b_v+I}$ where $a_h, b_h, a_v$ and $b_v$ are constants. $f_h$ and $f_v$ are continuous functions globally Lipschitz on $\mathbb{R}_+$ with Lipschitz constants $\dfrac{a_h}{b_h}$ and $\dfrac{a_v}{b_v}$, respectively. Their upper bounds are given by $a_h$ and $a_v$, respectively. For all numerical simulations, we consider the same set of parametric values given by $\mu_h = 20$, $r_h = 4$, $u = 5$, $\mu_v = 10$, $a_h = 300$, $b_h = 220$, $a_v = 60$, $b_v = 40$.

7.1.   Direct problem

For $\Lambda_h = 300$ and $\Lambda_v = 150$ then the reproduction number $\mathcal{R}_0 = 2.94 > 1$, and the solution of the model (2.1) converges to the equilibrium $E^*$ (Figure 2). This validate the global stability of $E^* = (S_h^*, I_h^*, R_h^*, S_v^*, I_v^*)$ when $\mathcal{R}_0 > 1$.

For $\Lambda_h = 120$ and $\Lambda_v = 40$ then the reproduction number $\mathcal{R}_0 = 0.41 < 1$, and the solution of the model (2.1) converges to the equilibrium $E_0 = (\dfrac{\Lambda_h}{\mu_h}, 0, 0, \dfrac{\Lambda_h}{\mu_h}, 0) = (6, 0, 0, 4, 0)$ (Figure 3). This confirms the global stability of $E_0 = (\dfrac{\Lambda_h}{\mu_h}, 0, 0, \dfrac{\Lambda_h}{\mu_h}, 0)$ when $\mathcal{R}_0 \leq 1$.

7.2.   Control problem

Concerning the resolution of the control problem, we used a numerical scheme that improves the Gauss-Seidel-like implicit finite-difference method where we applied a first-order backward-difference for the adjoint variables (see Appendix 8 for more details and ^{[18,23,24,25]} for other applications). We used the same parameters as in the direct problem (2.1): $\mu_h = 20$, $r_h = 4$, $\mu_v = 10$, $a_h = 300$, $b_h = 220$, $a_v = 60$, $b_v = 40$. Here $u$ is a variable such that the initial condition $u(0) = 1.5$ and with bounds $u^{\min} = 0$ and $u^{\max} = 3$.

In Figures 4, 5 and 6, the behaviours of $S_h(t), I_h(t), R_h(t), S_v(t)$ and $I_v(t)$ (right) and the adjoint states are plotted with respect to time for different values of $\alpha_1, \alpha_2$ and $\alpha_3$.

As we can see in Figures 4, 5 and 6, for the same parameters and initial conditions used for the resolution of the direct problem, the number of uninfected individuals increases and the number of infected individuals decreases which is the goal of the optimal strategy. This confirm that this strategy is very efficient and effective in the control of the number of infected individuals.

8.   Conclusions

We considered a 5D dynamical system modelling Zika virus dynamics in a human population with two general nonlinear saturated incidence rates. We take into account a treatment for the human population. We analysed the mathematical model in terms of local and global stability. We deduced that the endemic equilibrium point is globally asymptotically stable once $\mathcal{R}_0 > 1$ and the disease-free equilibrium point is globally asymptotically stable once $\mathcal{R}_0 \leq 1$. We studied the sensitivity of the model to all the parameters. An optimal strategy relative to this epidemic model based on the given treatment was considered. Finally, some numerical simulations were proposed to validate the obtained results.

The simplicity of the compartmental models makes it easy to compute, but also likely oversimplifies complex disease processes. This explain the limitations of compartmental models. In fact, these models assume that the population is homogeneously mixed, however, it is more realistic to consider heterogeneous number of contacts for each individual, i.e., network epidemic models ^[51]. Other open questions related to the limitations of compartmental models are also of great interest, those who deal with the role of time-delay and the role of intrinsic fluctuations.

Appendix

Adapted numerical scheme for solving the control problem

Let subdivide the interval of time $[0, T]$ such that $[0, T] = \cup_{n = 0}^{N-1} [t_n, t_{n+1}], \quad t_n = n \; dt, \quad dt = T/N$. Let $S_h^n, I_h^n, R_h^n, S_v^n, I_v^n, \lambda^n_1, \lambda^n_2, \lambda^n_3, \lambda^n_4, \lambda^n_5$ and $u^n$ be an approximation of $S_h(t), I_h(t), R_h(t), S_v(t), I_v(t)$, $\lambda_1(t)$, $\lambda_2(t)$, $\lambda_3(t)$, $\lambda_4(t)$, $\lambda_5(t)$ and the control $u(t)$ at the time $t_n$. $S_h^0, I_h^0, R_h^0, S_v^0, I_v^0$, $\lambda^0_1, \lambda^0_2, \lambda^0_3, \lambda^0_4, \lambda^0_5$ and $u^0$ as the initial values of the state variables, the adjoint variables and the control. $S_h^N, I_h^N, R_h^N, S_v^N, I_v^N$, $\lambda^N_1, \lambda^N_2, \lambda^N_3, \lambda^N_4, \lambda^N_5$ and $u^N$ as the final values of the state variables, the adjoint variables and the control. Concerning the resolution of the control problem, we used a numerical scheme that improves the Gauss-Seidel-like implicit finite-difference method where we applied a first-order backward-difference for the adjoint variables

Then the explicit expressions will take the following form

Conflict of interest

The authors declare no conflict of interests.

Acknowledgements

This Project was funded by the Deanship of Scientific Research (DSR) at King Abdulaziz University, Jeddah, Saudi Arabia under grant no. (G: 119-130-1442). The authors, therefore, acknowledge with thanks DSR for technical and financial support.

The authors are also grateful to the unknown referees for many constructive suggestions, which helped to improve the presentation of the paper.