Optimal control of a malaria model with long-lasting insecticide-treated nets

1.   Introduction

Malaria is a vector-borne disease of global public health concern with high level of morbidity and mortality in the tropical regions of the world. The disease is caused by several species of parasites of the Plasmodium genus type and transmitted to humans by the bite of a female anopheles mosquito when taking the blood meal necessary for egg production. In 2018, there was an estimated 228 million cases of malaria worldwide with about 405,000 deaths attributable to malaria^[39]. Various mathematical models of the transmission dynamics of malaria and its control have been proposed in the literature ^{[2,3,5,8,10,15,16,20,34,38]}. The very first model is that of Ross-MacDonald who laid the foundations for modeling malaria ^[26,35]. Models that include therapeutic (treatment and vaccination) and non-therapeutic (insecticide-treated bed net) measures have also flourished in the literature ^{[12,17,19,33]}. Because insecticide-treated nets (ITNs) reduce human/mosquito contacts, distribution campaigns have been organized in affected countries, including Cameroon. However, the use of these mosquitoes treated bednets have not always been satisfactory as several people let holes in the bednets, do not use them every night (poor adherence), or use these bednets for other activities such as fishing ^[36]. Mosquitoes insecticide-treated bednets could influence the force of infection, the rate of recruitment of new females mosquito, or the death rate of mosquitoes ^{[6,10,11,17,19]}. Moreover, the use of these treated bednets could influence the rate of loss of immunity.

We formulate a mathematical model for the transmission dynamics of malaria in human populations, which takes the (good or poor) use of bednets as a control measure. First, we formulate the autonomous model with a constant proportion of bednets usage as control strategy. Next, we compute the basic reproduction number $\mathcal{T}_0$, and investigate the existence and stability of the equilibria. Analytical results show that both model equilibria; the disease-free and the endemic states are locally asymptotically stable when $\mathcal{T}_0 <$ and when $\mathcal{T}_0 > 1$, respectively. However, the model could exhibit the phenomenon of backward bifurcation when $\mathcal{T}_0 < 1$, an epidemiological situation where, although necessary, having the basic reproduction number less than unity is not sufficient for malaria elimination ^[38].

We then extend our autonomous model by considering a time-dependent control of the proportion of bednets usage. Optimal control theory is used to establish conditions under which the spread of malaria can be mitigated. The characterization of the optimal control is obtained by the application of Pontryagin's maximum principle. We use the Forward-backward fourth-order Runge-Kutta method for numerical simulations to determine an optimal control strategy. In addition, we focus on a bednet control strategy since the other controls measures are expensive. By other vector controls, we mean outdoor application of larvicides (chemical or biological), breeding habitat reduction (e.g., draining standing water), outdoor vector control (mosquito fogging, attractive toxic sugar bait (ATSB)), indoor residual spraying (IRS), repellents, including topical repellents, mosquito coils, etc, rapid diagnosis and treatment (RDT), preventative drugs like seasonal malaria chemo-prevention (SMC), intermittent preventative treatment (IPT) ^[19]. Generally, the bednet control in the literature concerns the bednets usage, including insecticide-treated bed nets (ITNs), long-lasting insecticide-treated nets (LLINs), and untreated bednets (UBNs) ^[19]. Some of these other vector controls do not respectful the ecological population environment. But, the use of these other vector control require periodic actions on a short time (1 day, 1 week, 1 month, 1 year) that is less effective and practicable than the three years' use of the LLNs. We observe that the disadvantage (loss of immunity, even at its maximum) in the use of bednets is compensated by the benefit of the number of susceptible/infected individuals excluded from the malaria disease dynamics, the only danger being the poor use of the long-lasting bednets. In fact, the model suggests that if one family fails to use accurately the LLNs, then the number of infectious will increase (sometimes exponentially) since the protection of the bednets would drop down. However, it is possible to come close to the optimal results with a realistic strategy.

The rest of the paper is organized as follows. In Section 2, we present the mathematical model for malaria transmission dynamics with a parameter $w$ that represents the proportion of persons having and using the treated mosquito bednets correctly. In Section 4, we propose an optimal control problem for the minimization of the number of infected humans while controlling the cost of control interventions with bednets. Finally, in Section 4.4, some numerical simulations provided to support the analytical results are interpreted from the epidemiological point of view. Section 5 is the conclusion.

2.   Modeling of the LLITN controlled dynamics of malaria

2.1.   Description of involved phenomena

In this section, we consider two populations, namely human hosts and female mosquitoes that are assumed homogeneously distributed. We suppose that female mosquitoes only feed with human blood. The human population is subdivided into three classes, namely the susceptible $S_h$, the infectious $I_h$ and the immune $R_h$ as shown in Figure 1. According to ^[20], we do not consider an exposed compartment because it does not significantly influence qualitatively the evolution of infections in the human population. The recruitment is done only in the susceptible class at a rate $\Lambda_h$. In all compartments, there is an output of $\mu_h$ due to natural death, with an additional death rate $\mu_d$ in the infectious compartment. When in contact with an infectious mosquito, a susceptible humans become infected at a rate $\alpha_{h}$ representing the force of infection. Infectious humans recover and gain immunity at a rate $\delta_h$, while the rate of loss of immunity is $\gamma_h$.

Following the work in ^[17,19,37], we distinguish the questing anopheles (looking for the blood meal) from the resting one. As depicted in Figure 2, the new mosquitoes arrive through the compartment of questing susceptible at rate $\Lambda_v$. The natural death rate is $\mu_v$. When in contact with an infectious human, a susceptible mosquito can become infected at a rate of $\alpha_{v}$ corresponding to the strength of infection of the mosquitoes. Once infected, a mosquito will alternate (at most $12$ times) during its period of latency between the resting status and the questing status. Hence, we consider $6$ resting infected compartments and $6$ questing infected compartments. The transition rate from resting status to questing status is $\chi$, and the transition rate from questing to resting is $\beta$. Assuming that the number of bites on a human by female mosquitoes per day is $a$, the fraction of persons having a long-lasting insecticide-treated net is $b\in\left[0, 1 \right]$, and the sub-fraction of persons using it effectively is $u\in\left[0, 1 \right]$. Hence, the proportion of people who own a mosquito net and use it adequately is $w: = b\times u\in\left[0, 1 \right]$. We suppose that the dependence of some parameters with respect to $w$ can be express as

Indeed, the adequate use of mosquito bednet reduces transmission forces $\alpha_{h}$ and $\alpha_{v}$ through the factor $a$. It also limits indirectly the horizontal immigration of mosquitoes in the human environment via the rate $\Lambda_v$. Finally, $w$ indirectly increases the rate of death in the population of mosquitoes since it reduces mosquito-human contacts, and then there is a reduction of mosquito's blood meal leading to additional death rate.

The variables and parameters of the model are presented in detail in Table 1, Table 2 and Table 3.

2.2.   Mathematical model and preliminary properties

According to the above description and assumptions, we proposed the following system of non-linear ordinary differential equations (5).

By setting $\mathbf{x}_S = \left(S_h, S_q\right)\in \mathbb{R}^{2}$ and $\mathbf{x}_I \equiv \left((E_r^i, E_q^i)_{1 \leq i \leq 6}, E_r^7, I_r, I_q, I_h, R_h \right)\in \mathbb{R}^{17}$, the system (5) can be rewritten as

More precisely, we have

$\mathbf{A}_{SI}$ is the $2\times 17$ matrix with all its coefficients null, except $\mathbf{A}_{SI}\left(1, 15 \right) = -\dfrac{amS_h}{N_h}$, $\mathbf{A}_{SI}\left(1, 17 \right) = \gamma_h$, $\mathbf{A}_{SI}\left(2, 16 \right) = -\dfrac{acS_q}{N_h}$ and $\mathbf{A}_{SI}\left(2, 17 \right) = -\dfrac{\tilde{ac}S_q}{N_h}$.

$\mathbf{A}_{11}$ is a $13 \times 13$ matrix satisfying

$\mathbf{A}_{12}$ is the $13 \times 4$ matrix defined by $\mathbf{A}_{12}\left(1, 3 \right) = \dfrac{acS_q}{N_h}$, $\mathbf{A}_{12}\left(1, 4 \right) = \dfrac{a\tilde{c}S_v}{N_h}$ and $\mathbf{A}_{12}\left(i, j \right) = 0$ for other cases. The matrix $\mathbf{A}_{21}$ is the $4 \times 13$ matrix defined by $\mathbf{A}_{21}\left(2, 13 \right) = \chi$ and $\mathbf{A}_{21}\left(i, j \right) = 0$ for the other cases. The matrix $\mathbf{A}_{22}$ is the $4 \times 4$ square matrix defined by

Proposition 2.1. The non-negative cone $\left(\mathbb{R}_+\right) ^{19}$ is positively invariant for system $(7)$.

Proof 2.0.1. The result comes from the fact that $\mathbf{A(x)}$ is a Metzler matrix for all$\mathbf{x} \in \left(\mathbb{R}_+\right)^{19}$, i.e., $a_{i j} (\mathbf{x}) = \mathbf{A(x)}\left(i, j \right) \geq 0$ for $i \neq j$. Indeed, since the solution of the system $(7)$ is continuous, it suffices to check that every trajectory starting at the boundary of $\left(\mathbb{R}_+\right)^{19}$ remains in $\left(\mathbb{R}_+\right)^{19}$. That boundary is defined as $\cup _{i = 1}^{19} \mathcal{H}_i$ where $\mathcal{H}_i\equiv \left\lbrace \mathbf{x} \in \left(\mathbb{R}_{+}\right)^{19}, \mid x_i = 0\right\rbrace$. For $x \in \mathcal{H}_i$, $x'_i = \sum_{j = 1}^{19} a_{i j} (\mathbf{x})x_i +b_i = \sum_{j = 1, i \neq j }^{19} a_{i j} (\mathbf{x})x_i +b_i \geq 0$. This means that $x_i$ moves from $0$ to non-negative values.

Proposition 2.2. The set $\Omega = \bigg\lbrace \left(S_h, S_q, (E_r^i, E_q^i)_{1 \leq i \leq 6}, E_r^7, I_r, I_q, I_h, R_h \right)\in \mathbb{R}^{19}_{+} / 0 \leq N_h$ $\leq \dfrac{\Lambda_h}\mu_h, 0 \leq N_{v} \leq \dfrac{\Lambda_v}{\mu_v}\bigg\rbrace$ is a compact forward-invariant and absorbing set for the model system $(5)$.

Proof 2.0.2. Adding respectively the human sub-populations equations and then the mosquitoes sub-populations equations, we have

From Proposition 2.1, we have $N_h, N_v\geq 0$. Applying Gronwall's inequality, we have

Assuming that $N_h\left(0\right) \leq \dfrac{\Lambda_h}{\mu_h}$ and $N_v\left(0\right) \leq \dfrac{\Lambda_v}{\mu_v}$. it follows that $N_h \leq \dfrac{\Lambda_h}{\mu_h}$ and $N_v \leq \dfrac{\Lambda_v}{\mu_v}$.

Proposition 2.3. At any time $t$, the model system $(5)$ is well-posed, i.e., it admits a non-negative global unique $C^{1}$-solution on the set $\mathbb{R}_{+}$.

Proof 2.0.3. Applying the Cauchy-Lipschitz theorem, Proposition 2.1 and Proposition 2.2, the result is immediate.

3.   Asymptotic and bifurcation analyses

In this section we study the existence of equilibria and their stability using the theory of bifurcation.

3.1.   Disease-free equilibrium and threshold condition

This section is devoted to local and global stability of the disease-free equilibrium which unconditionally exists.

Proposition 3.1. System $(5)$ admits a disease-free equilibrium (DFE) given by $\mathbf{x}^{*} = (\mathbf{x}^{*}_S, \mathbf{x}^{*}_I)$ with

Proof 3.0.1. An equilibrium is obtained by solving $\mathbf{x}' = \mathbf{A}(\mathbf{x})\mathbf{x}+\mathbf{b} = 0$. A DFE corresponds to any solution satisfying $\mathbf{x}_I = 0_{\mathbb{R}^{17}}$.

Proposition 3.2. The system $\mathbf{x}' = \mathbf{A}_S(\mathbf{x}^{*}).\left(\mathbf{x}-\mathbf{x}_S^{*}\right)$ is globally asymptotically stable (GAS) at $\mathbf{x}^{*}_S$ on $\mathbb{R}^{2}_+$.

Proof 3.0.2. The proof is immediate since the matrix

has all its eigenvalues $- \mu_h$ and $- \tilde{\mu_v}$ negative.

In the following, we determine a stability threshold condition using a technique well described and used in ^[17,37]. In our case, this threshold can be biologically interpreted as the basic reproduction number $\mathcal{T}_0$ ^[19].

Let us define the threshold

Theorem 3.1. $\mathcal{T}_0$ is equivalent to the basic reproduction number $\mathcal{T}_0$ of the model system $(5)$. Moreover, the DFE is locally asymptotically stable (LAS) if $\mathcal{T}_0\leq 1$.

Proof 3.1.1. According to Gautier-Sallet's algorithm in ^[18] and Proposition 3.2, it is sufficient to obtain condition ensuring that $\mathbf{A}_{I}(\mathbf{x}^{\ast})$ is stable (all its eigenvalues are negative). Looking carefully at $\mathbf{A}_{11}$ and $\mathbf{A}_{22}$, they are Metzler stable. Thus, let$\mathbf{N} = \mathbf{A}_{22}(\mathbf{x}^{\ast}) -\mathbf{A}_{21}(\mathbf{x}^{\ast}) \times \mathbf{A}_{11}^{-1}(\mathbf{x}^{\ast}) \times \mathbf{A}_{12}(\mathbf{x}^{\ast})$. Using formal calculus (under Sagemath software for instance), we have

Let us consider the $2\times 2$ blocks in $\mathbf{N}$. Since $\mathbf{N}_{11}$ is Metzler stable, let $\mathbf{L} = \mathbf{N}_{22}-\mathbf{N}_{21}\times \mathbf{N}_{11}^{-1} \times \mathbf{N}_{12}$.

$\mathbf{A}_{I}(\mathbf{x}^{\ast})$ is stable if

Since $N^{\ast}_h = S^{\ast}_h$, this is equivalent to

Theorem 3.2. The DFE is GAS in $\Omega$ when $\mathcal{T}_0 < \dfrac{\mu_h}{\mu_h+\mu_d}\equiv \zeta$.

Proof 3.2.1. Our proof relies on Theorem 4.3 in ^[18], which establishes global asymptotic stability (GAS) for epidemiological systems that can be expressed in matrix form $(7)$. The demonstration is completely similar to that made in ^[37].

3.1.1.   Endemic equilibrium

Theorem 3.3. There exists $\mathcal{R}_{-}, \mathcal{R}_c, \mathcal{R}_{+}\in\mathbb{R}$ such that the model system $(5)$ has

(a) a unique endemic equilibrium if $\mathcal{R}_0 > 1$,

(b) a unique endemic equilibrium if $\mathcal{R}_0 = 1$ and $\mathcal{R}_c < 1$,

(c) two endemic equilibria if $\mathcal{R}_c < \mathcal{R}_0 < min(1, \mathcal{R}_{-})$ or $max(\mathcal{R}_{c}, \mathcal{R}_{+}) < \mathcal{R}_{0} < 1$,

(d) No endemic equilibrium elsewhere.

Proof 3.3.1. An endemic equilibrium is any non-zero and positive solution of the following system

The equations $(15)$ and $(16)$ allow us to write $S_h^{\star}$ and $R_h^{\star}$ in function of $I_h^{\star}$ as follows: $S_h^{\star} = \dfrac{\delta_h+ \mu_h + \mu_d}{\alpha_h^{\star}}I_h^{\star}$ and $R_h^{\star} = \dfrac{\delta_h}{\gamma_h + \mu_h}I_h^{\star}.$

To simplify the expressions, let $D = \delta_h + \mu_h + \mu_d$, $C = \dfrac{\delta_h}{\gamma_h + \mu_h}$ and $F = \gamma_h C$.

By subsequently replacing $S_h^\star$ and $R_h^\star$ by their values in $(14)$, we then obtain the expression of $I_h^\star$ with respect to $\alpha_h^\star$.

Also, $\alpha_v^{\star} = \dfrac{a(cI_h^{\star}+\tilde{c}R_h^{\star})}{N_h^{\star}} = \dfrac{a(C\tilde{c}+c)\alpha_h^{\star}}{(C+1)\alpha_h^{\star}+D}.$ Using equations $(17)$, $(18)$, $(19)$, $(20)$, $(21)$ and $(22)$, we have

$S_q^{\star} = \dfrac{\Lambda_v(w)}{\alpha_v^{\star}+(\mu_v+w\Delta_{\mu_v})}$, $E_r^{1 \star} = \dfrac{\alpha_v^{\star}S_q^{\star}}{\chi + \mu_v}$, $E_q^{i \star} = \dfrac{\chi E_r^{i \star}}{\beta +(\mu_v+w\Delta_{\mu_v})}\; \; for\; \; 1\leq i \leq 6$, $E_r^{j \star} = \dfrac{\beta E_q^{(j-1) \star}}{\chi +\mu_v}\; \; for\; \; 2\leq j \leq 7$; $I_q^{\star} = \dfrac{\chi f_q E_r^{7 \star}}{\beta(1-f_q f_r)}$ et $I_r^{\star} = \dfrac{\beta I_q^{\star}}{\chi + \mu_v}.$

So all our unknowns are expressed in terms of $\alpha_h^\star$, and it only remains to determine the value of $\alpha_h^\star$.

By definition, we have $\alpha_h\star = \dfrac{amI_q^\star}{N_h^\star}$and by replacing $I_q^\star$ and $N_h^\star$ by their values, and after simplification and re-arrangement, we obtain

where

Equation $(23)$ has solution $\alpha_h^{\star} = 0$ and solutions of the equation $(E): P_2 (\alpha_h^{\star})^2 + P_1 \alpha_h^{\star} + P_0 = 0$.

The case $\alpha_h^{\star} = 0$ leads us to the equilibrium without disease. We are interested in the equation $(E)$, of which we are going to analyze the number of positive solutions as a function of the value of $\mathcal{T}_0$.

1. If $\mathcal{T}_0 > 1$ then, $P_0 > 0$ and since $P_2 < 0$, the discriminant $\Delta = P_1^2 -4 P_2 P_0$ of the equation $(E)$ is positive, hence the equation $(E)$ has two different real solutions. In addition, the product of the solutions is $p = \dfrac{P_0}{P_2} < 0$. Hence, equation $(E)$ has a unique positive solution.

2. if $\mathcal{T}_0 = 1$, then, equation $(E)$ has two different real solutions, which are zero and $-\dfrac{P_1}{P_2}$. But $P_2 < 0$, so this solution is positive if $P_1 > 0$, that is to say if $\mathcal{T}_0 > \mathcal{T}_c$.

3. if $\mathcal{T}_0 < 1$ and $\Delta = P_{1}^2-4P_{2}P_{0} > 0$ and $\mathcal{T}_0 > \mathcal{T}_c$, then, equation $(E)$ admits two different positive solutions.

Let $P_2 = -b_2$ $P_1 = b_1(\mathcal{T}_0-\mathcal{T}_c)$ and $P_0 = b_0(\mathcal{T}_0-1)$; $b_2, \; \; b_1\; \; and\; \; b_0$ are all positive coefficients. We have, $\Delta = P_1^2-4P_2P_0 = b_1^2\mathcal{T}_0^{2}-(2b_1^{2}\mathcal{T}_c-4b_0b_2)\mathcal{T}_0 - 4b_0b_2+b_1\mathcal{T}_c^{2}$.

The last condition can be re-written as follows $\left\lbrace \begin{array}{l}\mathcal{T}_c < \mathcal{T}_0 < 1, \\\Delta = b_1^2\mathcal{T}_0^{2}-(2b_1^{2}\mathcal{T}_c-4b_0b_2)\mathcal{T}_0 - 4b_0b_2+b_1\mathcal{T}_c^{2} > 0.\end{array}\right.$

Let us study the sign of $\Delta$ in relation to the values of $\mathcal{T}_0$. Consider the equation

$(E_{\mathcal{T}_0})$ has as discriminant $\Delta_r = (2b_1^{2}\mathcal{T}_c-4b_0b_2)^2-4b_1^2(-4b_0b_2+b_1^{2}\mathcal{T}_c^{2}) = 16b_2b_0\left[ b_2b_0+b_1^{2}(1-\mathcal{T}_c)\right]$ which is positive for $\mathcal{T}_c < 1$, and the equation $(E_{\mathcal{T}_0})$ has two solutions $\mathcal{T}_{-}$ and $\mathcal{T}_{+}$.

We then have $\left\lbrace \begin{array}{l}\mathcal{T}_c < \mathcal{T}_0 < 1, \\\mathcal{T}_0 \in \left]-\infty, \mathcal{T}_{-} \right[ \cup \left]\mathcal{T}_{+}, +\infty \right[, \end{array}\right.$

which yields $\mathcal{T}_c < \mathcal{T}_0 < min(1, \mathcal{T}_{-})$, where $max(\mathcal{T}_{c}, \mathcal{T}_{+}) < \mathcal{T}_{0} < 1$.

Remarks 3.1. For the (global) stability of the endemic equilibrium, one could follow the approach in ^[19] by using a suitable Lyapunov-type functional along the positive flow of the model 5 on a "two domains" subdivision of the phase state $\mathbb{R}^{19}$, under appropriate conditions.

Note that the disease-free equilibrium is only globally asymptotically stable when $\mathcal{T}_0 < \zeta < 1$, so it is possible that if this condition is violated, bistability could occur. That is, for $\zeta < \mathcal{T}_0 < 1$, a stable DFE could co-exist with a stable endemic equilibrium, a phenomenon known as backward bifurcation ^{[7,9,14,21,38]}. In this case, the condition $\mathcal{T}_0 < 1$, although necessary is insufficient to mitigate the spread of malaria in a community.

The following figure depicts the backward bifurcation for model system (5), representing the plotting of $I_h^{\star}$ as a function of $\mathcal{T}_0$ for values of the bifurcation parameter $m$ ranging from $0.001$ to $0.01$, $\mu_d = 4.54 \times 10^{-3}$ and $\mathcal{T}_c = 0.388975$.

4.   Optimal control model

To find the solution to the model (5), following four steps are followed. (i) Description the optimal control (ii) Proof of the existence of optimal control problem (iii) Proof of the uniqueness of the optimal control (iv) Numerically solve the optimal control and show these graphically.

4.1.   Description of optimal control

There are several methods to mitigate the prevalence of malaria in a community by reducing the mosquito density, contact, longevity and competence. Among all these methods, the possession and correct use of insecticide-treated mosquito bednets is the strategy that considers three of the biological elements mentioned above ^[27]. Therefore, we consider as control the function $w$ representing the fraction of population that owns and properly use a mosquito net.

Consider the following objective functional

The terms $A_1 I_h$ and $A_2 \left(\sum\limits_{i = 1}^{6}E^{i}_{q}+ I_q\right)$ are the cost of infection while $Bw^{2}(t)$ is a quadratic cost related to the effort of using bednets. Our main goal is to find an optimal control function $w^{*}$ such that $J (w^{*}) = min\left\lbrace J(w) \mid w \in \Gamma(T)\right\rbrace$, with $\Gamma(T)$ the set of admissible controls defined as

4.2.   Existence and characterization of an optimal control

The aim of this section is to prove the existence of an optimal control for the model system (5) and then derive the optimality system. The existence of optimal control of the system (5) will be considered by applying the following theorem ^{[4,22,23,25,29]}.

Theorem 4.1. Consider the objective functional $J$ given by equation $(25)$, with $w \in \Gamma$ subject to the constraint state system $(5)$. There exists $w^{*} \in \Gamma(T)$ such that$J (w^{*}) = min\left\lbrace J(w) \mid w \in \Gamma(T)\right\rbrace$ subject to thecontrol system $(5)$ with initial conditions at $t = 0$.

Proof 4.1.1. The state and control variables of the system $(5)$ are positive and thecontrol set $\Gamma(T)$ is closed and convex. Therefore, the integrand of the objectivefunctional $J$ in which it was expressed in the system $(5)$ is a convex function of $\omega$ on the control set $\Gamma(T)$. Since the state solutions are bounded, the Lipschitz property of the state system with respect to the state variables is satisfied. It can also be seen that there exist positive numbers $\eta_1$ and $\eta_2$, and a constant $\varepsilon > 1$ such that $J (\omega) \geq \eta_1 \mid \omega \mid ^{\varepsilon}-\eta_2.$ Therefore, the state variables are bounded and the existence of optimal control of the model system $(5)$ is concluded.

4.3.   The uniqueness of optimal control

To derive the necessary conditions that the three controls and corresponding state variables must satisfy, we use Pontryagin's maximum principle ^[31]. To this end, we define the Hamiltonian function for the system, where $\lambda_i, \; \; i = 1, \ldots, 19$ are the adjoint variables or co-state variables

The following result presents the adjoint system and control characterization.

Theorem 4.2. Given an optimal control $w^*$, and the corresponding state solutions

of the corresponding state system $(5)$, there exists adjoint variables, $\lambda_i, \; \; i = 1, \ldots, 19$, satisfying

The control $w^{*}$ satisfies the optimality condition.

where $\alpha_v^{*} = \dfrac{mI_q^{*}}{N_h^{*}}$ and $\alpha_h^{*} = \dfrac{cI_h^{*}+\tilde{c}R_h^{*}}{N_h^{*}}.$

Proof 4.2.1. The differential equations governing the adjoint variables are obtained by differentiation of the Hamiltonian function, evaluated at the optimal control. Then, the adjoint system can be written as

with zero final time conditions (transversality) $\lambda_i(T) = 0$. The characterization of the optimal control given by $(28)$ is obtained by solving the equations on the interior of the control set, where $0 < w < 1$. That is,

with $\dfrac{\partial \mathbb{H}}{\partial \omega} = 0$, where $\alpha_v^{*} = \dfrac{mI_q^{*}}{N_h^{*}}$ and $\alpha_h^{*} = \dfrac{cI_h^{*}+\tilde{c}R_h^{*}}{N_h^{*}}.$Hence, we obtain

4.4.   Numerical simulations: the biological admissibility and approximate controls

We numerically solve the optimal transmission parameter control for the malaria model. The optimal control is obtained by solving the optimality system, consisting of $19$ non-linear ordinary differential equations from the state and adjoint equations. An iterative scheme is used for solving the optimality system ^[24]. For the simulations, we consider the initial (and arbitrary) number of individuals at time $t = 0$: $S_h(0) = 100000, I_h(0) = 100, R_h(0) = 1000, Sq(0) = 100000$, $E^{1}_r(0) = 10$, $E^{2}_r(0) = 9, E^{3}_r(0) = 8$, $E^{4}_r(0) = 7, E^{5}_r(0) = 6, E^{6}_r(0) = 5, E^{7}_r(0) = 4$, $E^{1}_q(0) = 3, E^{2}_q(0) = 3, E^{3}_q(0) = 3, E^{4}_q(0) = 3$, $E^{5}_q(0) = 3$, $E^{6}_q(0) = 2, I_r(0) = 35, I_q(0) = 800$.

For the cost weight in the objective functional $J$, we take $B = ＄4.5\; \; USD$ (for three years) which represents what the state of Cameroon spends on the purchase of an insecticide-treated mosquito net for two individuals ^[32]. This is somehow comparable to $B = ＄3.95\; \; USD$ for the average cost for a household (with about 5.5 individuals) per month (for the first two largest cities of Cameroon in terms of population - Douala and Yaounde ^[28]). The main practical problem is the difficulty to provide bednets to everybody in the household as well as individuals complain of feeling excessive heat when sleeping under a bednet ^[28], the latter being a potential reason why some individuals use other vector control measures. We could consider the cost per $household^2$ for $T = 3 years$, and finally discuss the impact of the optimal controls associated to $B_1 = 4.5*\frac{5.5}{2}$ and $B_2 = 3.95*36$. Clearly, from the economic stand point, the other vector control strategies are more expensive than the bednet control. Thus, we focus on the bednet control strategies. The simulations are carried out with $T$ = $3$ years, the duration of one LLIN (Long-lasting insecticide-treated bednet) efficacy ^[32]. We define the uniform control $u_{unif}(t) = k$ and the multi-intervals ("stage") one $u_{stages}(t) = \left\lbrace \begin{array}{ccc} u_1 & \text{ffl} & t\in\left[0;365.25\textbf{ days } \right], \\ u_2 & \text{ffl} & t\in\left]365.25;730.5\textbf{ months } \right], \\ u_3 & \text{ffl} & t\in\left]730.5;1080.75\textbf{ months } \right], \end{array} \right.$ over three years with $k, u_1, u_2, u_3\in\mathbb{R_+}$ and $\frac{1}{1080.75}\int^{1080.75}_0 u_{stages}(t)dt = k$ as the mean value. $u_{optimal}$ is an optimal control for our optimal problem in Theorem 4.1. $u_{forced}$ is an administrative control of distribution of the bednets over three years (it is imposed or "forced"); it is either $u_{unif}$ or $u_{stages}$. We also define the following in percentage:

$1.$ $t^{u_{optimal}}_s(u_{forced}) = \frac{100\times\textbf{Total of susceptible humans on [ 0;T] under }u_{forced} }{\textbf{ Total of susceptible humans on [ 0;T] under }u_{optimal}}$;

$2.$ $t^{u_{optimal}}_{I_h}(u_{forced}) = \frac{100\times\textbf{Total of infectious humans on [ 0;T] under }u_{optimal}}{\textbf{ Total of susceptible humans on [ 0;T] under }u_{forced}}$;

$3.$ $t^{u_{optimal}}_{R_h}(u_{forced}) = \frac{100\times\textbf{Total of recovered humans on [ 0;T] under }u_{optimal} }{\textbf{Total of recovered humans on [ 0;T] under }u_{forced}}$.

Definition 4.1. Let $\Gamma(T)$ be the set of admissible controls relative to a dynamical system $D_{(u(.))}$, $u(.)\in\Gamma(T)$. An optimal control, mathematically admissible, is biologically admissible if $t^{u_{optimal}}_s(u_{forced})\leq 100$, $t^{u_{optimal}}_{I_h}(u_{forced})\leq 100$ and $t^{u_{optimal}}_{R_h}(u_{forced})\leq 100$.

It is easy (even numerically) to study the biological admissibility to an (mathematically) admissible control $u_{forced}$. But (for all $u_{forced}$), the biological admissibility is a challenge related to the choice of the objective function.

Numerically, for $T = 1080.75$ days, $u_{unif}(t) = 0.65$ and $u_{stages}(t) = \left\lbrace \begin{array}{lll} u_1 = 0.9 & & t\in\left[0;365.25\textbf{ days } \right], \\ u_2 = 0.6 & & t\in\left]365.25;730.5\textbf{ days } \right], \\ u_3 = 0.45 & & t\in\left]730.5;1080.75\textbf{ days } \right]. \end{array} \right.$ For all our numerical simulations, graphs related to optimal control are in black solid lines, while those linked to the "forced" control are in solid green lines. The effects of the "uniform" control are graphically represented in Figures 4-9 while the "stage" control effects are shown in Figures 10-15.

Practically, the common strategies $u_{forced}$ in malaria affected countries are decreasing functions of time (as trends), due to the difficulty to maintain a constant or high ($> 90\%$) level of possession and use of bednets throughout a 3-year campaign of LLINs distribution.

Table 4 suggests that, even if the "stage" and uniform controls have the same mean, it is better to use a "stage" control with emphasis on the first few months of the 3 years. Clearly, the effort should be done to cover the gap (between the results of the system state following $u_{forced}$ and $u_{optimal}$) in Table 4. This approach, based on the reality of the malaria programs in each country, could support health policies and decision-makers in order to obtain an accurate threshold in the percentage $\epsilon_{u_{forced}}$ of the "administrative/public planners controls" $u_{forced}$ ($u_{stage}$ or $u_{unif}$) applications compared to optimal effects, such that $100-t^{u_{optimal}}_s(u_{forced})\leq \epsilon_{u_{forced}}$, $100-t^{u_{optimal}}_{I_h}(u_{forced})\leq \epsilon_{u_{forced}}$ and $100-t^{u_{optimal}}_{R_h}\leq \epsilon_{u_{forced}}$. This allows us to introduce the definitions of the $\epsilon_{u_{forced}}$-approximate weak or strong "sub-optimal" controls.

Definition 4.2. (Approximate controlability) Let $\Gamma(T)$ be the set of admissible controls relative to a dynamical system $D_{(u(.))}$, $u(.)\in\Gamma(T)$, for $T > 0$. For $V^{u_{optimal}}_{u_{forced}} = \left(t^{u_{optimal}}_s(u_{forced}), t^{u_{optimal}}_{I_h}(u_{forced}), t^{u_{optimal}}_{R_h}(u_{forced})\right)$, let define

also written as

and

That is,

A biologically admissible control $u_{forced}$ is $\epsilon_{u_{forced}}$-approximate weak "sub-optimal" if

A biologically admissible control $u_{forced}$ is $\epsilon_{u_{forced}}$-approximate strong "sub-optimal" if

Remarks 4.1. These definitions in 4.2 improve on the efficiency index ^[1]. It is possible to consider the reduction of noise $N_{mosq}$ (similar to $Norm_{.}$ for mosquitoes) produced by mosquitoes as the percentage of mosquitoes with optimal control compared to the states with forced control: then the new index would be $Norm^{\alpha, \beta}_{.}: = \alpha Norm_{.}+\beta N_{mosq}$ such that $\alpha+\beta = 1$. The coefficients $\alpha$ and $\beta$ relate respectively the importance of the humans' group and mosquitoes' group. Herein, we focus on the optimal impact on humans and consider $\alpha = 1$.

Straightforward computations lead to the following proposition 4.1.

Proposition 4.1. There is an equivalence between $Norm_{weak}$ and $Norm_{weak}$:

In Table 4,

and

We see that for $\epsilon_{u_{forced}} = 35\%$, the $V^{u_{optimal}}_{u_{unif}}$ is 35$\%$-approximate weak "sub-optimal" like $V^{u_{optimal}}_{u_{stage}}$. But only $V^{u_{optimal}}_{u_{stage}}$ is 35$\%$-approximate strong "sub-optimal" and not $V^{u_{optimal}}_{u_{unif}}$. Another interesting point is the fact that $Norm_{weak}(V^{u_{optimal}}_{u_{stage}}) = 18.186947$, and this comes from the fact that the "weak" deviation from the optimal strategy is only about $18.187\%$ in total (with the collective effort/contribution of $S_h, I_h$ and $R_h$ to reach the optimal strategy). By the way, $Norm_{strong}(V^{u_{optimal}}_{u_{stage}}) = 31.410$, and this corroborates the fact that the "weak" deviation from the optimal strategy is only about $31.410\%$ following the individual efforts/contributions of $S_h, I_h$ and $R_h$ to reach the optimal objective. We observe that the disadvantage (loss of immunity) in the use of bednets is compensated by the benefit of the number of susceptible/infected individuals excluded from the malaria disease dynamics. Then, it is possible to get closer to the optimal results with a realistic strategy.

5.   Conclusion

We formulated and rigorously analyzed a vector multi-stage malaria model with the use of mosquito treated bednets as preventive measure. The proposed model is biologically meaningful and mathematical well-posed. We investigated the local and global stability of equilibria. The analytical results reveal the possibility of bistability when $\mathcal{T}_0 < \zeta < 1$ (see subsection 3.1.1 with additional mortality $\delta_h$ less than $10^{-5}$ see discussion in ^[6]). That is, the model could exhibit the phenomenon of backward bifurcation, an epidemiological situation which although necessary, having the basic reproduction number less than unity is no longer sufficient to mitigate the malaria transmission dynamics ^[38]. Thus, a low level of additional (disease-induced) mortality could lead to the existence of an endemic equilibrium even if the basic reproduction number is less than one.

Next, an optimal control strategy is investigated with the proper usage of LLINs (during three years compared to a "forced" control) as the control parameter. We observe that the disadvantage (loss of immunity) in the use of bednets is compensated by the benefit of the number of susceptible/infected individuals excluded from the malaria disease dynamics. Moreover, it is possible to get close to the optimal results with a realistic strategy.

Results from this study could help inform health policy and decision-makers on the potential optimum strategies to mitigate malaria transmission dynamics in affected communities by designing reachable malaria program implementation objectives "close" to the optimal strategies $\epsilon_{u_{.}}\%$ by the "weak" collective contribution or the "strong" individual effort to achieve the optimal objective. The notions of $\epsilon_{u_{forced}}$-approximate strong/weak "sub-optimal" control are more practical than the theoretical optimal control which remains a daunting task to health officials. The upper bound $\epsilon_{u_{forced}}$ of the gap, from the "sub-optimal" results to the optimal ones, is of great interest practically since it delineates the acceptable error one could essentially make if we apply $u_{forced}$ instead of $u_{optimal}$.

The proposed model is not exhaustive. One could consider splitting the human population into adults and juveniles since malaria disproportionately affects children under five years of age. The model could then be fitted with real data from an affected region/country with most model parameter values estimated for the purpose. Also, because of uncertainties in most parameter values, a detailed sensitivity and uncertainty analysis could be carried out to understand the dependence of the basic reproduction number and model state variables on their components ^[30].

Conflict of interest

The authors declare that they have no conflicts of interest to this work.