Optimal control analysis of a malaria transmission model with applications to Democratic Republic of Congo

. In this paper, a dynamical model of malaria transmission with vector-bias and time-dependent controls is investigated. The controls include the RTS,S malaria vaccine, using insecticide-treated mosquito net, treatment of infectious human, and indoor spraying. For constant controls, the existence and stability of equilibrium, as well as the existence of backward bifurcation, are obtained. The sensitivity analysis quantiﬁes the impact of parameters and controls on the basic reproduction number. For time-dependent controls, by using the Pontryagin’s maximum principle the existence and expression of optimal controls are established. As an application of the model and control strategies, the malaria transmission and controls in Democratic Republic of Congo are discussed. To be speciﬁc, we simulate the reported cases of Democratic Republic of Congo by World Health Organization and predict the trends. Cost-effectiveness analysis and numerical simulations show that combining all controls can minimize the number of infected humans to the full extent, using insecticide-treated mosquito net is the most cost-effectiveness strategy, combining RTS,S malaria vaccine with using insecticide-treated mosquito net and treatment of infectious human is also cost-effective among all the strategies with good effect.


Introduction
Malaria, caused by protozoan parasites belonging to the genus Plasmodium, is a mosquitoborne disease, which affects health systems and economies significantly. According to the World Health Organization (WHO) 2021 reports, 241 million cases of malaria and 627 thousand malaria deaths occurred globally in 2020, 14 million cases more and 69,000 more deaths compared to 2019 [23]. Humans get infected from the bite of infected female Anopheles mosquitoes, especially, from these infected Plasmodium falciparum and Plasmodium vivax. In the WHO African Region, Plasmodium falciparum accounts for 99.7% of estimated malaria cases, while Plasmodium vivax is responsible for 74.1% of malaria cases in the WHO Region of Americas.
Take Democratic Republic of Congo (Congo, DR), an African country, for example. Based on data in Table 1, we can get Fig. 1, which shows that the death rate keeps decreasing from 1% to 0.1% during 2010-2020 with the advancement of society. However, we can see that in Congo, DR, there are still about 5% people infected with malaria in the year 2010, and that the number of population infected with malaria in Congo, DR keeps increasing surprisingly to about 25% in 2020. Therefore, controlling the transmission of malaria in Congo, DR is desperately in need.
Mathematical models have become vital tools for understanding the dynamics of infectious diseases long time ago. The first mathematical model depicting the transmission process of malaria was introduced by Ross [17] and refined by MacDonald [14] later. From then on, the malaria transmission model was developed extensively [4,15,16]. Particularly, the vector-bias effect, namely, the greater attractiveness of infectious humans to mosquitoes than susceptible ones [5,12], was firstly introduced to a malaria transmission model by Kingsolver [11] in 1987. Subsequently, feeding bias by vectors toward infected hosts and incubation time in mosquitoes included, Hosack et al. [7] found Table 1. Reported cases (from [23]) and the population (from [20]) of Congo, DR.

Year Reported cases Death cases Population
Year https://www.journals.vu.lt/nonlinear-analysis parasite modified behavior by a refined malaria model. Further, Chamchod and Britton [3] improved the previous models by defining the attractiveness in a different way. Based on these works, a lot of researchers investigated the malaria models with vector bias [1,10]. All these results show that the vector bias has an important impact on the epidemiology of malaria.
Furthermore, based on the optimal control theory, malaria transmission models are also used in the decision-making of prevention and control of malaria [8,9,15,16,18]. Kim et al. [10] applied two optimal controls (treatment and media awareness) and indicated that the combination of the two controls was the most effective strategy to monitor the disease. It should be pointed that few researcher focus on controlling the malaria with vaccine since the uncertainty of malaria vaccination. However, recently, the World Health Organization has recommended to use the RTS,S malaria vaccine (a vaccination mainly for children under 5) broadly since if the RTS,S vaccine introduced widely and urgently, tens of thousands of children's lives could be saved every year [23].
In this paper, we developed a malaria transmission model with vector-bias effect and time-dependent controls: RTS,S vaccine, using insecticide-treated mosquito net, treatment of infectious human, indoor spraying simultaneously. We will use our model to simulate the reported cases of Congo, DR and propose several strategies to control the spread of malaria in Congo, DR. The research topic in this paper is novel and characteristic. In fact, few research analyze the control strategies with area-specific parameters.
The paper is arranged as follows. In Section 2, we formulate the nonautonomous malaria model with vector bias. Then we analyze the dynamics of the autonomous version in Section 3. Specifically speaking, the existence and stability of the disease-free and endemic equilibria, the bifurcation analysis, as well as the sensitivity analysis, are given. Furthermore, optimal control analysis of the nonautonomous model is performed in Section 4. In Section 5, we apply the model to the malaria transmission in Congo, DR. To be specific, we simulate the reported cases of Congo, DR by WHO and predict the trends in five years. Numerical simulations are used to observe the outcomes, characterizations as well as the cost of these strategies. Furthermore, cost-effectiveness analysis and numerical simulations show that combining all controls can minimize the number of infected humans to the full extent, using insecticide-treated mosquito net is the most cost-effectiveness strategy, combining RTS,S malaria vaccine with using insecticidetreated mosquito net and treatment of infectious human is also cost-effective among all the strategies with good effect. Finally, conclusions are summarized in Section 6.

Model formulation
The human population N h (t) is divided into three compartments, that is, susceptible S h (t), infected I h (t), and recovery R h (t). The susceptible will be vaccinated RTS,S (a malaria vaccine for children under 5) by rate of u 1 (t) ∈ [0, 0.4], δ is the efficiency of vaccination. Assume that Λ 1 is the birth or immigration rate of human population. They either die naturally or diminished following infection with malaria at a rate (infection rate) λ h = σβ 1 (1 − u 2 (t))I m /(νI h + S h + R h ), where σ is the biting rate, β 1 is the transmission probability of per bite, ν is the bias parameter, and u 2 (t) ∈ [0, 0.89] refers to a time-dependent control function representing the decrease of transmission rate by the use of insecticide-treated mosquito nets. Infectious individuals are assumed to recover at a rate γ + τ u 3 (t), where γ is the spontaneous rate, u 3 (t) ∈ [0, 0.88] also is a timedependent control function representing the treatment of infectious human with malaria symptoms, and τ is the efficacy to treatment. Infectious individuals who does not recover die naturally at rate d 1 , and die of the disease at rate α.
The female anopheles mosquito population N m (t) is divided into two compartments, susceptible S m (t) and infectious I m (t). For mosquito population, the recruitment rate is Λ 2 , the infection rate is where β 2 is the probability for a vector to get infected by an infectious human. Infectious mosquito die at a rate d 2 + cu 4 (t), where d 2 is the natural death rate, and u 4 (t) ∈ [0, 1] is the control function on mosquito population by spraying insecticide.
Based on the above details, the malaria mathematical model is given by the following nonautonomous and nonlinear equations: where all parameters are positive and described in Table 2.  Firstly, as a basic property of solutions for model (1), we have the following lemma.
for all t > 0 also are positive. Furthermore, The proof of this lemma is omitted for simplicity.

Analysis of model with constant controls
In this section, we assume that all the controls u i (i = 1, 2, 3, 4) in model (1) are constants. We will investigate the dynamical behaviors of model (1), including the calculation of basic reproduction number, the existence of equilibria and their stability, the backward bifurcation analysis, and the sensitivity analysis.

Equilibrium and stability
It is easy to see that model (1) has a unique disease-free equilibrium Using the next generation matrix method (see [19]), we can easily obtain the basic reproduction number of model (1) .
Based on the result in [19], the following theorem is established immediately.
Theorem 1. The disease-free equilibrium E 0 is locally asymptotically stable if R 0 < 1, and is unstable if R 0 > 1.
be the endemic equilibrium of model (1). Denote Then, after some simple computations, we get Therefore, we can get λ * h satisfies the quadratic equation Thus, the existence of endemic equilibrium of model (1) is equivalent to the existence of positive roots of equation (2). Based on Vieta's theorem, when either of the following conditions hold, then equation (2) has positive roots It is easily to see that b 2 > 0 (=, <) is equivalent to R 0 < 1 (=, >). For the convenience, define the constants: https://www.journals.vu.lt/nonlinear-analysis Then ∆ > 0 (=, <) is equivalent to R 0 > R 1 (=, <), and b 1 < 0 (=, >) is equivalent to R 0 > R 2 (=, <). Therefore, we can establish the following conclusion on the existence of endemic equilibrium for model (1).
Remark 1. Theorem 2 shows that model (1) can generate a backward bifurcation under the condition R c < R 0 < 1. The analysis will be given in Section 3.2 in detail.
Moreover, when R 0 > 1, the local stability of endemic equilibrium E * can be obtained by the linearization method. The characteristic equation of linearized system at equilibrium E * and the corresponding coefficients c i (i = 1, 2, 3, 4, 5) in characteristic equation are presented in Appendix. We have the following conclusion.
Theorem 3. When R 0 > 1, then endemic equilibrium E * is locally asymptotically stable if the following conditions corresponding to Hurwitz criterion hold: Table 2

Bifurcation analysis
In this subsection, the forward and backward bifurcations are investigated for model (1). The main method on the bifurcation analysis is presented in [2] based on the centre manifold theory in [2]. We can establish the following result.
Proof. We choose the transmission rate β 1 as a bifurcation parameter. Clearly, R 0 = 1 is equivalent to According to Theorem 1, equilibrium E 0 is locally stable if β 1 < β * 1 , and unstable if It is not difficult to obtain that the Jacobian matrix of system (4) and the corresponding characteristic equation is It is clear that equation (5) admits a simple zero eigenvalue and all other eigenvalues are negative.
In the following, we calculate the right eigenvector w = (w 1 , w 2 , w 3 , w 4 , w 5 ) T and the left eigenvector Calculate the constants a and b as follows: From system (4) we can obtain , , Therefore, we have Clearly, b > 0 and a > 0 (< 0) if and only if A > 0 (< 0). Thus, from the results in [2] we can obtain that the conclusions in Theorem 6 hold. This completes the proof.

Sensitivity of basic reproduction number to parameters and controls
Sensitivity analysis is a technique that permits exploration of complex models by evaluating how the quantities of interest (QOI) change with the variation of parameters of interest (POI). Generally, the growth of an epidemic is partly characterized by the basic reproduction number R 0 . For this, we analytically calculate the sensitivity indices of R 0 to parameter p as follows: [4]: In addition, the sign of the sensitivity index suggests whether the quantities of interest, R 0 increases (> 0) or decreases (< 0) with the parameter of interest p.
The sensitivity index of R 0 with respect to parameters Λ 1 , Λ 2 , σ, β 1 , β 2 , and ν are constant as shown in Table 3. For others, we get where A 2 = γ+d 1 +α+τ u 3 , which are related to controls. The ranges of these indices are summarized in Table 3, which shows the most sensitive parameter is the average number of times one mosquito would bite a human per year, σ, while the least sensitive parameter is the disease induced death rate of human, α. In general, R 0 decreases with the increase of parameters γ, d 2 , and Λ 1 . Besides, R 0 increases with the increase of parameters d 1 , Λ 2 , β 1 , β 2 , and ν.   The sensitivity indices of R 0 with respect to the constant controls are given by Applying parameter values from Table 2, we get Fig. 3, which shows that the effect of the four controls are closely related to the control levels. To be specific, using of bednets u 2 affects transmission mostly, followed by u 3 , u 4 . Besides, from Fig. 4 we know that using of insecticide-treated bednets and treatment of infectious can reduce the value of basic reproduction number under 1.

Remark 4.
It should be emphasized that to eliminate the disease, we must adjust some constant controls u i large enough to make sure R 0 < R c (but R 0 < 1) since the backward bifurcation may happen when R c < R 0 < 1. Furthermore, in the reality, it is not easy to keep the controls constant all the time. We know the fact that constant controls can be regarded as an approximation of time-dependent ones. Therefore, in the next section, we consider the model with time dependent controls and investigate the optimal controls according to optimal control theory.

Optimal control
The purpose of this section is to achieve the optimal control of model (1) with time dependent controls. We introduce the feasible control set as follows: The optimal control for model (1) aims at minimizing the number of malaria-infected humans I h and mosquito I m and the cost needed during the intervention period. Therefore, the objective functional is defined by is the solution of model (1), A 1 , A 2 are the weight constants with respect to the number of infected humans and infected mosquitoes; ξ i 0 (i = 1, 2, 3, 4) are the weight constants on the benefit and cost; ξ i u 2 i /2 (i = 1, 2, 3, 4) is the cost of corresponding control.
In the following, we need to solve the optimal control problem: find a control u * ∈ ∆ such that min J(u) = J(u * ). Firstly, the existence of an optimal control is settled by the following theorem.
Theorem 5. There exists an optimal control u * ∈ ∆ for the objective functional J(u) subject to model (1) with positive initial conditions such that J(u * ) is minimal.
Theorem 5 can be proved by using the similar arguments as in [6,16]. Hence, we here omit the proof for simplicity.

Applications to malaria transmission in Congo, DR
As shown in Fig. 1 and Table 1, the situation of malaria spread in Congo, DR is quite serious. To make suitable control strategies, we will first simulate the reported cases with model (1). Then we can get that the basic reproduction number of Congo, DR is estimated to be R 0 . = 1.6038 with parameter values in Table 2. According to Theorem 2, the disease in Congo, DR will be endemic, and Fig. 5(b) further verifies it. Therefore, it is necessary to carry out research on optimal control strategies, so that reasonable responses can be made from prevention and treatment to eradicate the disease as soon as possible and minimize the loss as much as possible.

Data and parameters estimation
In this subsection, model (1) is used to simulate the reported annual malaria human cases in WHO [23]. Here we use parameter values in Table 2. Some parameters values are chosen based on references, and some are to match the data. We explain part of them in the following.
1. The Birth rate of Congo, DR is 4.103% [21], and the total population of Congo, DR in 2020 is 8.9561 · 10 7 , so that the recruitment is 3.67 · 10 6 . 2. The number of mosquito biting per day is 0.25-0.4 [4]. To fit the model better, we take biting 0.3 · 365 per year. 3. The life span of human in Congo, DR is 62 [22]. So the corresponding death rate is 1/62. 4. The average life expectancy of adult mosquito is about 15 to 20 days. Here we take 1/d 2 to be 15/365 year. 5. The recovery rate γ is 3 months to 50 years assumed by [4]. Here we take the immunity period to malaria of humans 1/γ to be 2 years. 6. The average disease induced death rate is 3.3 · 10 −3 by data in Table 1.
Based on these parameter values, we carry out the numerical simulations of our model and obtain a reasonable match between the infected human of model (1) and the malaria data of Congo, DR from 2010 to 2020 in Fig. 5. It indicates that the transmission of malaria in Congo, DR has not arrived at a stable period yet and the disease will become more serious without further control measures.
Moreover, on the basis of parameter values in Table 2, the basic reproduction number of Congo, DR are estimated to be 1.6038. Therefore, according to Theorem 2, the disease in Congo, DR will be endemic and Fig. 5(b) further verifies it. Thus, it is necessary to carry out research on optimal control strategies, so that reasonable responses can be made from prevention and treatment to eradicate the disease as soon as possible and minimize the loss as much as possible.

Effect of some control strategies
To eradicate the spread of malaria in Congo, DR as soon as possible and to minimize the corresponding loss as much as possible, we carry out research on optimal control, so the reasonable responses can be made by the government of Congo, DR. In this section, we propose the following strategies to control the spread of malaria in Congo, DR. In the following, we obtain the optimal controls numerically by solving model (1), adjoint system (6), (7) and using the characterization of optimal controls (8) by the forward-backward sweep method [13]. In detail, the forward fourth-order Runge-Kutta method is used to solve the state system, and backward fourth-order Runge-Kutta method is used for solving the adjoint system. The adjoint system is solved under the initial assumption of zero controls and obtained solutions of the state system. The controls are updated by taking average of the previous result and the characterizations in (8). This condition continues repeatedly up to the consecutive iteration are negligibly close.
Remark 5. It needs to point out that the parameter values used above are in Table 2, which are closely related to the situation in Congo, DR, and the initial values of human population above are taken based on real data in Table 1.
For all the scenarios, we can see from Figs. 6, 7, 9, 10 that the magnitudes of infected humans reduce to lower levels in a way. The control profiles of u * 1 , u * 2 , u * 3 , and u * 4 are also depicted in Fig. 6 and Figs. 8-10.
Moreover, from Fig. 11 and Table 4 we can get the following result: using RTS,S only is the cheapest but the effect is poor; combining all the controls could reduce the infected population to the full extent but is the most expensive; the longer these controls go on the more effective they are. Remark 6. Actually, for strategy A, that is, using the RTS,S vaccine u 1 only, the effect is poor; see Figs. 6(a), 6(b). This implies that using RTS,S vaccine only is not good enough, people need to change their way of life and apply other interventions such as using insecticide-treated mosquito net in the daily life. Although the interventions raised above all have good outcomes to some extent, due to limited resources in a country or a city, we need to evaluate which intervention is more economical. Cost-effectiveness analysis is widely used to do this evaluation [15,16,18].

Cost-effectiveness analysis
In this subsection, two approaches are applied to evaluate which intervention is more economical, that is, the average cost-effectiveness ratio (ACER) and incremental costeffectiveness ratio (ICER). Cost-effectiveness analysis is widely used to do this evaluation [15,16,18].

Average cost-effectiveness ratio (ACER)
The average cost-effectiveness ratio (ACER) deals with a single intervention and evaluates that intervention against its baseline option, namely, no intervention.

ACER =
Total cost produced of scenarios i Total number of infection averted of scenarios i .
From Table 4 we can see the ACER of each scenarios, and Fig. 11(c) shows more clearly that the most cost-effective strategy is strategy B, and strategy A is the least costeffective one. To further investigate the cost-effectiveness of the various control strategies, we evaluated the incremental cost-effectiveness ratio (ICER).

Incremental cost-effectiveness ratio (ICER)
The incremental cost effectiveness ratio (ICER)

ICER =
Difference in infection averted costs in scenarios i and j Difference in total number of infection averted in scenarios i and j is used to measure up the changes between costs and benefits of two alternative control strategies and is generally described as the additional cost per additional health outcome. The numerator describe the difference of costs 4 of interventions i and j, while the denominator is the difference of health outcomes. To calculate ICER, we rank the control strategies in increasing order effectiveness according to the total infection averted (see Table 4). Firstly, compare the cost effectiveness of strategy A and strategy G. We have From ICER(A) and ICER(G) we can see that strategy A is strongly dominated 5 and strategy G saves 1.2387 · 10 −4 than strategy A. Therefore, it is better to exclude strategy A from the set of control strategies and alternative interventions to implement to keep the limited resources. Then we compute the ICER between strategy G and strategy D and get 4 T 0 (1/2)Σ 4 i=1 ξ i u 2 i dt is the total cost. ICER(G) = 1.5708 · 10 −4 , ICER(D) = 0.0094, so we exclude strategy D. Repeating this process up to the final strategy, we obtain that strategy B is the most cost-effectiveness strategy, which is consistent with the result by ACER. However, from Fig. 9(b) we can see directly that strategy B reduce 32.73% less infected population than strategy H, K, N, O. In other words, from Fig. 11 we can see that strategy H, K, N, O cause reduction of infected population to the most extent, which is what Congo, DR needs at present. Besides, Fig. 9(b) and Table 4 indicate that strategy O cause reduction of infected population mostly. So we compare the incremental costeffectiveness ratio between these four strategies and found that strategy K is the most costeffectiveness. Therefore, if the government see reducing the number infected population as the primary goal, then strategy K is a good choice. Now the policy maker must decide which strategy to use.

Conclusions and discussions
As we all know, various control measures are applied to wipe out or control malaria since the great loss caused by malaria for the world. Recently, there are some good news on malaria vaccine, that is, WHO recommend to use the RTS,S malaria vaccine broadly [23]. In this paper, controls of RTS,S vaccine, using insecticide-treated mosquito net, treatment of infectious humans, indoor spraying are incorporated to a malaria transmission model with vector-bias effect.
Our model is used to simulate the reported human malaria cases of Congo, DR. Fifteen strategies are proposed to control the spread of malaria in Congo, DR. In addition, costeffectiveness analysis suggests that the use of insecticide-treated mosquito net minimizes malaria infection and costs needed for implementation. In summary, based on our analysis, we have some suggestions on the control for malaria in Congo, DR. c 3 = (q + ρ + d 1 ) A 1 (A 6 + d 2 + cu 4 ) − A 2 A 7 + d 1 m(A 3 + A 6 + q + d 2 + cu 4 ) with m = ρ + d 1 + δu 1 . Thus, by the Hurwitz principle we have Theorem 3 on the local stability of the endemic equilibrium.