# Modeling and Optimal Control Analysis for Malaria Transmission with Role of Climate Variability

Authors: Temesgen Duressa Keno; Lemessa Bedjisa Dano; Oluwole Daniel Makinde
Journal: Computational and Mathematical Methods
 (2022)

---

## Abstract

In this paper, we present a nonlinear deterministic mathematical model for malaria transmission dynamics incorporating climatic variability as a factor. First, we showed the limited region and nonnegativity of the solution, which demonstrate that the model is biologically relevant and mathematically well-posed. Furthermore, the fundamental reproduction number was determined using the next-generation matrix approach, and the sensitivity of model parameters was investigated to determine the most affecting parameter. The Jacobian matrix and the Lyapunov function are used to illustrate the local and global stability of the equilibrium locations. If the fundamental reproduction number is smaller than one, a disease-free equilibrium point is both locally and globally asymptotically stable, but endemic equilibrium occurs otherwise. The model exhibits forward and backward bifurcation. Moreover, we applied the optimal control theory to describe the optimal control model that incorporates three controls, namely, using treated bed net, treatment of infected with antimalaria drugs, and indoor residual spraying strategy. The Pontryagin’s maximum principle is introduced to obtain the necessary condition for the optimal control problem. Finally, the numerical simulation of optimality system and cost-effectiveness analysis reveals that the combination of treated bed net and treatment is the most optimal and least-cost strategy to minimize the malaria.


---

## Body

## 1. Introduction

Malaria is a vector-borne disease caused by a parasite calledPlasmodium and transmitted between humans and mosquitos via bites of infected female Anopheles mosquitoes. It is a major public health problem, particularly on the African continent [1]. This disease can also be transmitted by blood transfusion or congenital. According to the world malaria report published in December 2021, an estimated 241 million cases and 627000 deaths occurred worldwide, with the African Region accounting for 95 percent of all cases in 2020 [2]. Climate variability is recognized to have a significant impact on the malaria vector’s life cycle, and the life of the mosquito is governed by temperature and rainfall [3]. The most effective way to avoid the malaria dynamics are treated bed net, treatment with antimalaria drugs, and indoor residual spraying [4].The mathematical model for malaria dynamics is the most useful in comprehending the disease’s presence in the human population. Ross devised the first mathematical model for malaria dynamics transmission [5]. Various mathematical models based on Ross’s basic malaria model were described by different scholars with different factors such as including the exposed class in mosquitoes and humans [6–8] and the impact of climate variations on malaria epidemics in terms of death rate, birth rate, and prevalence of mosquito population [9–11].Many investigations of malaria dynamics models were presented using optimum control problems employing optimal control theory. For example, Agusto et al. [12] investigated optimum solutions for decreasing the dynamics of malaria using a malaria transmission model with an optimal control problem. The authors proposed that the combination of the treated bed net, therapy, and indoor residual spraying is the optimum technique for illness reduction. Okosun et al. [13] created a malaria model that includes optimum control techniques of malaria disease transmission with two measures that regulate treatment infectious and immunization. The authors determined that the combination of therapy and immunization is the most effective and least expensive way to control a malaria condition. Leiton et al. [14] presented an SEIRS-SI optimal control model for malaria transmission in Colombia considering three optimal control strategies. The authors concluded that integrated of control measures treated bed net, intermittent prophylactic treatment in pregnancy and effective case management is the best strategy to prevent the malaria dynamics. Olaniyi et al. [15] formulated a malaria transmission dynamic model with an optimal control analysis using four time dependent continuous controls. The authors concluded that the combination of all these controls is the best measure to minimize the spread of malaria. Keno et al. [16] formulated an optimal control and cost effectiveness strategies of malaria transmission with temperature variability factor by considering three optimal control strategies. The model analysis provides that the combination of treatment and indoor residual spraying is the most efficient and less costly strategy to minimize the malaria. Olaniyi et al. [17] proposed SEIRS for the human population and SEI malaria model with optimal control and cost-effectiveness analysis in the presence of reinfection and relapse in malaria dynamics. The authors suggested that the strategy with a combination of treated bed net and indoor residual spraying is the most effective and least cost to eradicate the malaria.To the best of our knowledge, the impact of climate variability on malaria epidemics with optimal control and cost effectiveness analysis with a logistic growth of climate variation with respect to mosquitoes breeding and malaria infection was not considered in their models. In the malaria transmission model [11], we considered impact of temperature variability on malaria epidemics. In this paper, we extended model [11] to the impact of climate variability (temperature and rainfall) with respect to mosquitoes breeding rate and malaria infection with optimal control and cost effectiveness analysis.The following is the format of this paper: in Section2, we develop a malaria transmission model that incorporates the impact of climate variability on malaria epidemics. The model’s qualitative analysis is described in Section 3. We perform sensitivity analysis in Section 4. The optimal control analysis of the malaria transmission model is described in Section 5. In Section 6, we compared numerical simulation results to analytical results. The cost-effectiveness analysis is depicted in Section 7. Finally, Section 8 discusses the study’s conclusion.


## 2. Model Description and Formulation

In this section, we proposed and developed a malaria transmission model that takes into account the human and vector populations, with the total human population denoted byNht and the vector population denoted by Nmt. The SEIRS model describes the human population and divides it at time t into the following subpopulations: susceptible human Sht, exposed human Eht, infected human Iht, and recovered human Rht. The SEI model describes the mosquito population, and the total mosquito population at time t, denoted by Nmt, is subdivided into susceptible mosquito Smt, exposed mosquito Emt, and infected mosquito Imt. The total human and mosquito populations are then given by Nht=Sht+Eht+Iht+Rht and Nmt=Smt+Emt+Imt, respectively. Assume that all of the parameters in the model are positive. The new recruits (assumed to be susceptible) enter the human population by birth (migration) at the rate Ψ. An infectious mosquito transmits malaria to the susceptible human, and the susceptible human moves to the exposed human at a climate dependent rate of βhT,R, β0h is the contact rate of humans with mosquitos when there is no variation in climate, and β1h is the incremental contact rate of humans with mosquitos due to climate variation. The progression of exposed humans to infected humans is αh, the natural death rate of all humans is μh, and induced death is δ. The infected human recovered due to the use of antimalarial drugs through a treatment rate of γh. The recovered populations of humans have short period immunity that can be lost and become susceptible to reinfection rate of ωh. The recruitment of mosquito populations at a rate of ΦT,R is climate-dependent; Φ0 is the birth rate of mosquitos when there is no variation in climate, and ϕ1m is the incremental birth rate of mosquitos when there is variation in climate. Mosquito population contact with an infected human with malaria and moved to exposed class at a climate dependent rate of βmT,R, β0m is mosquito contact rate with human when there is no variation in climate, and β1m is incremental mosquito contact rate with a human due to climate variation. The exposed mosquito moved at a rate of αm to an infected mosquito. Every mosquito dies at the same rate μm. Mosquitoes do not recover from malaria because an infected mosquito is infectious until it dies. The temperature growth rate r follows a logistic function, m is the temperature dependent rate of precipitation, Tmax is the maximum temperature for the mosquito to be most active, and T0 is the minimum temperature for the mosquito to be less active. All the description of parameters are given in Table 1, and the diagram of malaria disease transmission is shown in Figure 1.Table 1
Parameters description used for the model (1).
ParametersParameter descriptionμhHuman death rate due to natural causeδInduced death rate of humanγhRecover rate of infected humansΨHuman recruitment rateΦ0Mosquito recruitment rateμmMosquito natural death rateωhHuman immune deficiency rateβ1mMosquito contact rate (incremental)β1hHuman contact rate (incremental)Φ1mIncremental birth rate of mosquitoβ0hContact rate of human with mosquitoβ0mContact rate of mosquito with humanαhExposed human progression rateαmExposed mosquito progression ratemTemperature dependent rate of precipitationz∗Percentage of an antibody produced by a humanm∗Percentage of an antibody produced by a mosquitorTemperature increase rateT0Minimum temperature for the mosquito to be less activeTmaxMaximum temperature for the mosquito to be most activeFigure 1
Flow diagram of malaria transmission.The transmission dynamics of malaria is described by the following system of nonlinear differential equations, based on the flow diagram depicted in Figure1:
(1)dShdt=Ψ−βhT,RShIm1+z∗Im−μhSh+ωhRh,dEhdt=βhT,RShIm1+z∗Im−μh+αhEh,dIhdt=αhEh−μh+δ+γhIh,dRhdt=γhIh−μh+ωhRh,dSmdt=ΦT,R−βmT,RSmIh1+m∗Ih−μmSm,dEmdt=βmT,RSmIh1+m∗Ih−μm+αmEm,dImdt=αmEm−μmIm,dTdt=r1−TTmaxT−T0,dRdt=mdTdt⇒Rt=mTt−T0+ε,where βhT,R=β0h+β1hk,βmT,R=β0m+β1mk,ΦT,R=Φ0+Φ1mk and k=1+mT−T0/Tmax.
(2)WithinitialconditionsSh0=Sh0,Eh0=Eh0,Ih0=Ih0,Rh0=Rh0,Sm0=Sm0,Em0=Em0,Im0=Im0,T0=Tm0,R0=Rm0.Base on the fact that as average temperatures at the earth’s surface rise more, evaporation occurs, which in turn increases overall precipitation (rainfall). Therefore, a warming climate is expected to increase precipitation (rainfall) in many areas. This implies that rainfall pattern is an increasing function of temperature. Note thatm is the temperature dependent rate of precipitation, and ε is the amount of precipitation at T=T0. In essence,
(3)R0=εatT=T0andRmax=mTmax−T0+R0.


## 3. Qualitative Analysis of the Model

### 3.1. Invariant Region

The invariant region is used to determine where the model’s solution is bounded. The model (1) is divided into two parts: the host population and the vector population. The total human population is represented by Nht=Sht+Eht+Iht+Rht. Then, by differentiating Nht both sides with respect to time and combining the first four equations from the model (1), we get
(4)ddtSh+Eh+Ih+Rh=Ψ−μhNh−δIh.Then, equation (4) becomes
(5)ddtSh+Eh+Ih+Rh≤Ψ−μhNh.By solving the equation (5), we can see that 0≤Nh≤Ψ/μh. Thus, for the human population, the bounded region of the system (1) is given by
(6)Ωh=Sh,Eh,Ih,Rh∈ℝ+4:0≤Sh+Eh+Ih+Rh≤Ψμh.Similarly, the equation gives the total mosquito population in the system (1) is
(7)Nmt=Smt+Emt+Imt.Obviously, the result we obtain after differentiatingNmt with respect to time t is
(8)ddtSm+Em+Im=ΦT,R−μmNm.By solving the equation (8), we get 0<Nm≤ΦT,R/μm. As just a result, for the mosquito population, the invariant region of the system (1) is given by
(9)Ωm=Sm,Em,Im∈ℝ+3:0≤Sm+Em+Im≤ΦT,Rμm,is a positive invariant. As a result, the system’s invariant region (1) is given by
(10)Ω=Ωh×Ωm=Sh,Eh,Ih,Rh,Sm,Em,Im∈ℝ+7:Nh≤Ψμh,Nm≤ΦT,Rμm,is a positive invariant set, and all of the solution set of system (1) is bounded in Ω within the region.


### 3.2. Positivity of the Solution

The purpose of this subsection is to demonstrate that all solutions of the model (1) will remain nonnegative in the future if their initial data is nonnegative.Theorem 1.
The model solution (1) given by Sht,Eht,Iht,Rht,Smt,Emt and Imt with nonnegative initial conditions Sh0,Eh0,Ih0,Rh0,Sm0,Em0 and Im0 remain nonnegative for all time t≥0.Proof.
First, take the begin equation from system (1) as given by
(11)dShdt=Ψ−βhT,RShIm1+z∗Im−μhSh+ωhRh.
Consequently, this indicates that(12)dShdt≥−βhT,RIm1+z∗Im+μhSh.
We obtain by integrating equation (12) with respect to time and solving it using the technique of variable separation with initial condition.
(13)Sht≥Sh0exp−μht+∫0tβhT,RImx1+z∗Imxdx≥0.
The other state variablesEht,Iht,Rht,Smt,Emt,andImt are nonnegative for all time t≥0 by the same procedure. As a result, the malaria transmission model stated in system (1) is both epidemiologically significant and mathematically well posed in Ω.


### 3.3. Basic Reproduction NumberR0

To find the steady-state solution of the model (1), we equated the left-hand side of the system (1) to zero with the value of Eh=0,Ih=0,Em=0,Im=0, and obtained the disease free-equilibrium of system (1) denoted by E1 or E2, where
(14)E1=Ψμh,0,0,0,ΦT0μm,0,0,T0,R0orE2=Ψμh,0,0,0,ΦTmaxμm,0,0,Tmax,Rmax.The basic reproduction number is defined as the average amount of secondary infectious caused by a primary infectious over a given time period [18], and it has been calculated using the next-generation matrix method [19]. Then, to get R01 and R02, at E1 and E2, respectively, we rewrite the model (1) beginning with newly infective classes of human and mosquito population as
(15)dEhdt=βhT,RShIm1+z∗Im−μh+αhEh,dIhdt=αhEh−μh+δ+γhIh,dEmdt=βmT,RSmIh1+m∗Ih−μm+αmEm,dImdt=αmEm−μmIm.The right hand side of system (15) can then be written as f−v, where
(16)f=βhT,RShIm1+m∗Im0βmT,RSmIh1+m∗Ih0andv=μh+αhEhμh+δ+γhIh−αhEhμm+αmEmμmIm−αmEm.The Jacobian matrices off and v at the disease-free equilibrium points give F and V, respectively, where
(17)F=000βhT∗,R∗Ψμh00000βmT∗,R∗ΦT∗,R∗μm000000andV=αh+μh000−αhμh+δ+γh0000αm+μm000−αmμm.The next-generation matrix from the product of equation (17) calculated by FV−1 is obtained as
(18)FV−1=00βhT∗,R∗Ψαmμhμmαm+μmβhT∗,R∗Ψμhμm0000βmT∗,R∗ΦT∗,R∗αhμmαh+μhαh+γh+μhβmT∗,R∗ΦT∗,R∗μmμh+δ+γh000000.Thus, the fundamental reproduction numberR0=ρFV−1 where ρ is the dominant eigenvalue of the product FV−1 and the R01 and R02 at disease free equilibrium E1 and E2 are obtained, respectively, in equations (19) and (20) as given by
(19)R01=β0hΨβ0mΦ0αhαmμhμm2μh+δ+γhαh+μhαm+μm,(20)R02=β0h+β2hΨβ0m+β2mΦ0+Φ2mαhαmμhμm2μh+δ+γhαh+μhαm+μm,where β2h=β1hk,β2m=β1mk,Φ2m=Φ1mk, and k=1+mTmax−T0/Tmax.As a result, the basic reproduction number (R02) at maximum temperature can be written in terms of R01 in the equation (21), as described below:
(21)R02=R012+Φ0β0mβ2mΨαhαm+Ψβ0m+β2mβ2hΦ0+Φ2m+Φ2mβ0hαhαmμhμm2μh+δ+γhαh+μhαm+μm,where R01 denotes the basic reproduction number at T0 for the mosquito to be less active in breeding.


### 3.4. Local Stability of Disease-Free Equilibrium

Theorem 2.
IfR01<R02<1; then, the disease free equilibrium point(s) of the system (1) is locally asymptotically stable in Ω.Proof.
The Jacobian matrix of equation (1) at the disease-free equilibrium point is given as
(22)zE∗=z1100ωh00−βhT∗,R∗Ψμh0z220000βhT∗,R∗Ψμh0αhz33000000γhz4400000−βmT∗,R∗ΦT∗,R∗μm0z550000βmT∗,R∗ΦT∗,R∗μm00z66000000αmz77,where z11=−μh,z22=−αh+μh,z33=−μh+δ+γh,z44=−μh+ωh,z55=−μm,z66=−αm+μm,z77=−μm.
Settingc1=αh+μh,c2=μh+δ+γh,c3=αm+μm,c4=μm and from the equation (22) the Jacobian matrix obtained as the polynomial function given by
(23)−λ−μh−λ−μm−λ−ωh+μhλ4+b1λ3+b2λ2+b3λ+b4=0,where
(24)b1=c1+c2+c3+c4,(25)b2=c1c2+c2c3+2c1c4+c3c4+c2c4,(26)b3=2c1c2c3+c1c2c4+c2c3c4,(27)b4=c1c2c3c4−βhT∗,R∗ΨβmT∗,R∗ΦT∗,R∗αhαmμmμh.
From equation (23), we obtain
(28)λ1=−μh<0,λ2=−μm<0,λ3=−ωh+μh<0,and we get from the final polynomial equation,
(29)λ4+b1λ3+b2λ2+b3λ+b4=0.
Using the Routh-Hurwitz criteria [20, 21], we can see that all of the equation’s eigenvalues (29) have negative roots or imaginary roots with a negative real part if
(30)b1>0,b2>0,b3>0,b4>0,b1b2−b3>0,b1b2b3−b32>0,b1b2b3−b32−b4b12>0.
Clearly, we have seen thatb1>0,b2>0, and b3>0, because they are a sum of positive parameters and at T∗=Tmax the value of b4 is described by
(31)b4=μmγh+δ+μhαh+μhαm+μm−β0h+β2hΨβ0m+β2mΦ0+Φ2mαhαmμmμh=1−R022.
However, forb4 to be positive, 1−R022 must also be positive, resulting in R02<1. Furthermore, when R02>1, then, b4<0, implying that disease-free equilibrium (DFE) is unstable. As a result, since R01<R02, the DFE is locally asymptotically stable in Ω if R01<R02<1.


### 3.5. Global Stability of Disease-Free Equilibrium

Theorem 3.
IfR01<R02<1; then, the disease free-equilibrium point(s) of the system (1) is globally asymptotically stable in Ω.Proof.
We used the technique implemented by Lyapunov theorem [19]; first, we developed the following Lyapunov function defined as
(32)M=αhαh+μhμh+δ+γhEh+Ihμh+δ+γh+μmR02β3mΦ3mEm+μmαm+μmR02β3mΦ3mαmIm.
Then differentiating the Lyapunov function with respect to timet, the obtained result is given by
(33)dMdt=αhαh+μhμh+δ+γhdEhdt+1μh+δ+γhdIhdt+μmR02β3mΦ3mdEmdt+μmαm+μmR02β3mΦ3mαmdImdt,=αhαh+μhμh+δ+γhβ3hShIm−μh+αhEh+1μh+δ+γhαhEh−μh+δ+γhIh+μmR02β3mΦ3mβ3mSmIh−μm+αmEm+μmαm+μmR02β3mΦ3mαmαmEm−μmIm,=β3hαhShImαh+μhμh+δ+γh−Ih+μmSmIhR02Φ3m−μm2αm+μmR02Imβ3mαmΦ3m,≤β3hαhΨαh+μhμh+δ+γhμh−μm2αm+μmR02β3mαmΦ3mIm+R02−1Ih,=β3hαhΨμm2αm+μmαh+μhμh+δ+γhμhβ3mαmΦ3mIm+IhR02−1,where β3h=β0h+β2h,β3m=β0m+β2m, and Φ3m=Φ0+Φ2m.
Obviously,dM/dt<0 if R02<1 and dM/dt=0 iff Eh=0,Ih=0,Em=0,Im=0. This demonstrates that the dominant compact invariant set in Sh,Eh,Ih,Rh,Sm,Em,Im∈Ω:dM/dt=0 represents the singleton set DFE in Ω. As a result of LaSalle’s invariant principle [22], every solution which begins in the region approaches DFE as t (time) tends to infinity. Since R01<R02, the DFE is globally asymptotically stable in Ω if R01<R02<1.


### 3.6. Malaria Present Equilibrium

The endemic equilibrium pointE∗=Sh∗,Eh∗,Ih∗,Rh∗,Sm∗,Em∗,Im∗,T∗,R∗ is the condition in which malaria is found in the human population, which can be obtained by equating all the model equations of system (1) to zero. At this point, let λh∗=βhT∗,R∗Im∗ and λm∗=βmT∗,R∗Ih∗ be force of infection of human and mosquito, respectively, and from the system (1), the endemic equilibrium point at T∗,R∗=T0,R0 is provided by
(34)Sh∗=Ψαh+μhμh+δ+γhμh+ωhαh+μhμh+δ+γhμh+ωhλh∗+μh−ωhαhγhλh∗,(35)Eh∗=Ψμh+δ+γhμh+ωhλh∗αh+μhμh+δ+γhμh+ωhλh∗+μh−ωhαhγhλh∗,(36)Ih∗=Ψαhμh+ωhλh∗αh+μhμh+δ+γhμh+ωhλh∗+μh−ωhαhγhλh∗,(37)Rh∗=Ψαhγhλh∗αh+μhμh+δ+γhμh+ωhλh∗+μh−ωhαhγhλh∗,(38)Sm∗=Φ0λm∗+μm,(39)Em∗=λm∗Φ0λm∗+μmαm+μm,(40)Im∗=λm∗αmΦ0λm∗+μmαm+μmμm.By substitutingIh∗ and Im∗ from equation (34) into λm∗ and λh∗, respectively, and λh∗ is obtained by solving the equation
(41)A1λh∗2+B1λh∗=0,where
(42)A1=μmαm+μmβ0mαhΨμh+ωh+μmμh+δ+γhαh+μhμh+ωh−μmωhγhαh,(43)B1=ωh+μhμhμm2αh+μhαm+μmμh+δ+γh1−R012.Therefore,A1>0 and B1>0 whenever R01<1, so that λh∗=−B1/A1<0. As a result, positive endemic equilibrium does not occur when R01<1. This shows that endemic equilibrium is found if B is less than zero which means that there is endemic equilibrium for the model. Similarly, applying the same procedure, the endemic equilibrium point the model (1) with T∗,R∗=Tmax,Rmax is described in the following as obtain:
(44)Sh∗=Ψαh+μhμh+δ+γhμh+ωhαh+μhμh+δ+γhμh+ωhλh∗+μh−ωhαhγhλh∗,(45)Eh∗=Ψμh+δ+γhμh+ωhλh∗αh+μhμh+δ+γhμh+ωhλh∗+μh−ωhαhγhλh∗,(46)Ih∗=Ψαhμh+ωhλh∗αh+μhμh+δ+γhμh+ωhλh∗+μh−ωhαhγhλh∗,(47)Rh∗=Ψαhγhλh∗αh+μhμh+δ+γhμh+ωhλh∗+μh−ωhαhγhλh∗,(48)Sm∗=Φ0+Φ2mλm∗+μm,(49)Em∗=λm∗Φ0+Φ2mλm∗+μmαm+μm,(50)Im∗=λm∗αmΦ0+Φ2mλm∗+μmαm+μmμm.Similarly, substitutingIh∗ and Im∗ from equation (44) into λm∗ and λh∗, respectively, and λh∗ is obtained by solving the equation
(51)Cλh∗2+Dλh∗=0,where
(52)C=μmαm+μmβ0m+β2mαhΨμh+ωh+μmμh+δ+γhαh+μhμh+ωh−μmωhγhαh,(53)D=ωh+μhμhμm2αh+μhαm+μmμh+δ+γh1−R022.As a result of equation (21), if R02<1, it follows that R01<1 exists, and DFE exists for both R01 and R02. However, R01<1 does not automatically imply R02<1, because the value of R02 can be greater than 1, implying that R02 may present backward bifurcation while R01 can only present forward bifurcation.


## 3.1. Invariant Region

The invariant region is used to determine where the model’s solution is bounded. The model (1) is divided into two parts: the host population and the vector population. The total human population is represented by Nht=Sht+Eht+Iht+Rht. Then, by differentiating Nht both sides with respect to time and combining the first four equations from the model (1), we get
(4)ddtSh+Eh+Ih+Rh=Ψ−μhNh−δIh.Then, equation (4) becomes
(5)ddtSh+Eh+Ih+Rh≤Ψ−μhNh.By solving the equation (5), we can see that 0≤Nh≤Ψ/μh. Thus, for the human population, the bounded region of the system (1) is given by
(6)Ωh=Sh,Eh,Ih,Rh∈ℝ+4:0≤Sh+Eh+Ih+Rh≤Ψμh.Similarly, the equation gives the total mosquito population in the system (1) is
(7)Nmt=Smt+Emt+Imt.Obviously, the result we obtain after differentiatingNmt with respect to time t is
(8)ddtSm+Em+Im=ΦT,R−μmNm.By solving the equation (8), we get 0<Nm≤ΦT,R/μm. As just a result, for the mosquito population, the invariant region of the system (1) is given by
(9)Ωm=Sm,Em,Im∈ℝ+3:0≤Sm+Em+Im≤ΦT,Rμm,is a positive invariant. As a result, the system’s invariant region (1) is given by
(10)Ω=Ωh×Ωm=Sh,Eh,Ih,Rh,Sm,Em,Im∈ℝ+7:Nh≤Ψμh,Nm≤ΦT,Rμm,is a positive invariant set, and all of the solution set of system (1) is bounded in Ω within the region.


## 3.2. Positivity of the Solution

The purpose of this subsection is to demonstrate that all solutions of the model (1) will remain nonnegative in the future if their initial data is nonnegative.Theorem 1.
The model solution (1) given by Sht,Eht,Iht,Rht,Smt,Emt and Imt with nonnegative initial conditions Sh0,Eh0,Ih0,Rh0,Sm0,Em0 and Im0 remain nonnegative for all time t≥0.Proof.
First, take the begin equation from system (1) as given by
(11)dShdt=Ψ−βhT,RShIm1+z∗Im−μhSh+ωhRh.
Consequently, this indicates that(12)dShdt≥−βhT,RIm1+z∗Im+μhSh.
We obtain by integrating equation (12) with respect to time and solving it using the technique of variable separation with initial condition.
(13)Sht≥Sh0exp−μht+∫0tβhT,RImx1+z∗Imxdx≥0.
The other state variablesEht,Iht,Rht,Smt,Emt,andImt are nonnegative for all time t≥0 by the same procedure. As a result, the malaria transmission model stated in system (1) is both epidemiologically significant and mathematically well posed in Ω.


## 3.3. Basic Reproduction NumberR0

To find the steady-state solution of the model (1), we equated the left-hand side of the system (1) to zero with the value of Eh=0,Ih=0,Em=0,Im=0, and obtained the disease free-equilibrium of system (1) denoted by E1 or E2, where
(14)E1=Ψμh,0,0,0,ΦT0μm,0,0,T0,R0orE2=Ψμh,0,0,0,ΦTmaxμm,0,0,Tmax,Rmax.The basic reproduction number is defined as the average amount of secondary infectious caused by a primary infectious over a given time period [18], and it has been calculated using the next-generation matrix method [19]. Then, to get R01 and R02, at E1 and E2, respectively, we rewrite the model (1) beginning with newly infective classes of human and mosquito population as
(15)dEhdt=βhT,RShIm1+z∗Im−μh+αhEh,dIhdt=αhEh−μh+δ+γhIh,dEmdt=βmT,RSmIh1+m∗Ih−μm+αmEm,dImdt=αmEm−μmIm.The right hand side of system (15) can then be written as f−v, where
(16)f=βhT,RShIm1+m∗Im0βmT,RSmIh1+m∗Ih0andv=μh+αhEhμh+δ+γhIh−αhEhμm+αmEmμmIm−αmEm.The Jacobian matrices off and v at the disease-free equilibrium points give F and V, respectively, where
(17)F=000βhT∗,R∗Ψμh00000βmT∗,R∗ΦT∗,R∗μm000000andV=αh+μh000−αhμh+δ+γh0000αm+μm000−αmμm.The next-generation matrix from the product of equation (17) calculated by FV−1 is obtained as
(18)FV−1=00βhT∗,R∗Ψαmμhμmαm+μmβhT∗,R∗Ψμhμm0000βmT∗,R∗ΦT∗,R∗αhμmαh+μhαh+γh+μhβmT∗,R∗ΦT∗,R∗μmμh+δ+γh000000.Thus, the fundamental reproduction numberR0=ρFV−1 where ρ is the dominant eigenvalue of the product FV−1 and the R01 and R02 at disease free equilibrium E1 and E2 are obtained, respectively, in equations (19) and (20) as given by
(19)R01=β0hΨβ0mΦ0αhαmμhμm2μh+δ+γhαh+μhαm+μm,(20)R02=β0h+β2hΨβ0m+β2mΦ0+Φ2mαhαmμhμm2μh+δ+γhαh+μhαm+μm,where β2h=β1hk,β2m=β1mk,Φ2m=Φ1mk, and k=1+mTmax−T0/Tmax.As a result, the basic reproduction number (R02) at maximum temperature can be written in terms of R01 in the equation (21), as described below:
(21)R02=R012+Φ0β0mβ2mΨαhαm+Ψβ0m+β2mβ2hΦ0+Φ2m+Φ2mβ0hαhαmμhμm2μh+δ+γhαh+μhαm+μm,where R01 denotes the basic reproduction number at T0 for the mosquito to be less active in breeding.


## 3.4. Local Stability of Disease-Free Equilibrium

Theorem 2.
IfR01<R02<1; then, the disease free equilibrium point(s) of the system (1) is locally asymptotically stable in Ω.Proof.
The Jacobian matrix of equation (1) at the disease-free equilibrium point is given as
(22)zE∗=z1100ωh00−βhT∗,R∗Ψμh0z220000βhT∗,R∗Ψμh0αhz33000000γhz4400000−βmT∗,R∗ΦT∗,R∗μm0z550000βmT∗,R∗ΦT∗,R∗μm00z66000000αmz77,where z11=−μh,z22=−αh+μh,z33=−μh+δ+γh,z44=−μh+ωh,z55=−μm,z66=−αm+μm,z77=−μm.
Settingc1=αh+μh,c2=μh+δ+γh,c3=αm+μm,c4=μm and from the equation (22) the Jacobian matrix obtained as the polynomial function given by
(23)−λ−μh−λ−μm−λ−ωh+μhλ4+b1λ3+b2λ2+b3λ+b4=0,where
(24)b1=c1+c2+c3+c4,(25)b2=c1c2+c2c3+2c1c4+c3c4+c2c4,(26)b3=2c1c2c3+c1c2c4+c2c3c4,(27)b4=c1c2c3c4−βhT∗,R∗ΨβmT∗,R∗ΦT∗,R∗αhαmμmμh.
From equation (23), we obtain
(28)λ1=−μh<0,λ2=−μm<0,λ3=−ωh+μh<0,and we get from the final polynomial equation,
(29)λ4+b1λ3+b2λ2+b3λ+b4=0.
Using the Routh-Hurwitz criteria [20, 21], we can see that all of the equation’s eigenvalues (29) have negative roots or imaginary roots with a negative real part if
(30)b1>0,b2>0,b3>0,b4>0,b1b2−b3>0,b1b2b3−b32>0,b1b2b3−b32−b4b12>0.
Clearly, we have seen thatb1>0,b2>0, and b3>0, because they are a sum of positive parameters and at T∗=Tmax the value of b4 is described by
(31)b4=μmγh+δ+μhαh+μhαm+μm−β0h+β2hΨβ0m+β2mΦ0+Φ2mαhαmμmμh=1−R022.
However, forb4 to be positive, 1−R022 must also be positive, resulting in R02<1. Furthermore, when R02>1, then, b4<0, implying that disease-free equilibrium (DFE) is unstable. As a result, since R01<R02, the DFE is locally asymptotically stable in Ω if R01<R02<1.


## 3.5. Global Stability of Disease-Free Equilibrium

Theorem 3.
IfR01<R02<1; then, the disease free-equilibrium point(s) of the system (1) is globally asymptotically stable in Ω.Proof.
We used the technique implemented by Lyapunov theorem [19]; first, we developed the following Lyapunov function defined as
(32)M=αhαh+μhμh+δ+γhEh+Ihμh+δ+γh+μmR02β3mΦ3mEm+μmαm+μmR02β3mΦ3mαmIm.
Then differentiating the Lyapunov function with respect to timet, the obtained result is given by
(33)dMdt=αhαh+μhμh+δ+γhdEhdt+1μh+δ+γhdIhdt+μmR02β3mΦ3mdEmdt+μmαm+μmR02β3mΦ3mαmdImdt,=αhαh+μhμh+δ+γhβ3hShIm−μh+αhEh+1μh+δ+γhαhEh−μh+δ+γhIh+μmR02β3mΦ3mβ3mSmIh−μm+αmEm+μmαm+μmR02β3mΦ3mαmαmEm−μmIm,=β3hαhShImαh+μhμh+δ+γh−Ih+μmSmIhR02Φ3m−μm2αm+μmR02Imβ3mαmΦ3m,≤β3hαhΨαh+μhμh+δ+γhμh−μm2αm+μmR02β3mαmΦ3mIm+R02−1Ih,=β3hαhΨμm2αm+μmαh+μhμh+δ+γhμhβ3mαmΦ3mIm+IhR02−1,where β3h=β0h+β2h,β3m=β0m+β2m, and Φ3m=Φ0+Φ2m.
Obviously,dM/dt<0 if R02<1 and dM/dt=0 iff Eh=0,Ih=0,Em=0,Im=0. This demonstrates that the dominant compact invariant set in Sh,Eh,Ih,Rh,Sm,Em,Im∈Ω:dM/dt=0 represents the singleton set DFE in Ω. As a result of LaSalle’s invariant principle [22], every solution which begins in the region approaches DFE as t (time) tends to infinity. Since R01<R02, the DFE is globally asymptotically stable in Ω if R01<R02<1.


## 3.6. Malaria Present Equilibrium

The endemic equilibrium pointE∗=Sh∗,Eh∗,Ih∗,Rh∗,Sm∗,Em∗,Im∗,T∗,R∗ is the condition in which malaria is found in the human population, which can be obtained by equating all the model equations of system (1) to zero. At this point, let λh∗=βhT∗,R∗Im∗ and λm∗=βmT∗,R∗Ih∗ be force of infection of human and mosquito, respectively, and from the system (1), the endemic equilibrium point at T∗,R∗=T0,R0 is provided by
(34)Sh∗=Ψαh+μhμh+δ+γhμh+ωhαh+μhμh+δ+γhμh+ωhλh∗+μh−ωhαhγhλh∗,(35)Eh∗=Ψμh+δ+γhμh+ωhλh∗αh+μhμh+δ+γhμh+ωhλh∗+μh−ωhαhγhλh∗,(36)Ih∗=Ψαhμh+ωhλh∗αh+μhμh+δ+γhμh+ωhλh∗+μh−ωhαhγhλh∗,(37)Rh∗=Ψαhγhλh∗αh+μhμh+δ+γhμh+ωhλh∗+μh−ωhαhγhλh∗,(38)Sm∗=Φ0λm∗+μm,(39)Em∗=λm∗Φ0λm∗+μmαm+μm,(40)Im∗=λm∗αmΦ0λm∗+μmαm+μmμm.By substitutingIh∗ and Im∗ from equation (34) into λm∗ and λh∗, respectively, and λh∗ is obtained by solving the equation
(41)A1λh∗2+B1λh∗=0,where
(42)A1=μmαm+μmβ0mαhΨμh+ωh+μmμh+δ+γhαh+μhμh+ωh−μmωhγhαh,(43)B1=ωh+μhμhμm2αh+μhαm+μmμh+δ+γh1−R012.Therefore,A1>0 and B1>0 whenever R01<1, so that λh∗=−B1/A1<0. As a result, positive endemic equilibrium does not occur when R01<1. This shows that endemic equilibrium is found if B is less than zero which means that there is endemic equilibrium for the model. Similarly, applying the same procedure, the endemic equilibrium point the model (1) with T∗,R∗=Tmax,Rmax is described in the following as obtain:
(44)Sh∗=Ψαh+μhμh+δ+γhμh+ωhαh+μhμh+δ+γhμh+ωhλh∗+μh−ωhαhγhλh∗,(45)Eh∗=Ψμh+δ+γhμh+ωhλh∗αh+μhμh+δ+γhμh+ωhλh∗+μh−ωhαhγhλh∗,(46)Ih∗=Ψαhμh+ωhλh∗αh+μhμh+δ+γhμh+ωhλh∗+μh−ωhαhγhλh∗,(47)Rh∗=Ψαhγhλh∗αh+μhμh+δ+γhμh+ωhλh∗+μh−ωhαhγhλh∗,(48)Sm∗=Φ0+Φ2mλm∗+μm,(49)Em∗=λm∗Φ0+Φ2mλm∗+μmαm+μm,(50)Im∗=λm∗αmΦ0+Φ2mλm∗+μmαm+μmμm.Similarly, substitutingIh∗ and Im∗ from equation (44) into λm∗ and λh∗, respectively, and λh∗ is obtained by solving the equation
(51)Cλh∗2+Dλh∗=0,where
(52)C=μmαm+μmβ0m+β2mαhΨμh+ωh+μmμh+δ+γhαh+μhμh+ωh−μmωhγhαh,(53)D=ωh+μhμhμm2αh+μhαm+μmμh+δ+γh1−R022.As a result of equation (21), if R02<1, it follows that R01<1 exists, and DFE exists for both R01 and R02. However, R01<1 does not automatically imply R02<1, because the value of R02 can be greater than 1, implying that R02 may present backward bifurcation while R01 can only present forward bifurcation.


## 4. Model Parameter Sensitivity Analysis

Malaria eradication strategies should focus on key parameters that have a significant impact on the basic reproduction number. The parameters found in the model (1) influence the basic reproduction number. Because those parameters can increase or decrease a basic reproduction number R0 as their values change, and vice versa. The sensitivity analysis is used to identify the parameters that have a large influence on the basic reproduction number R0. To perform the sensitivity analysis, the method described in [23–25] is used.Definition 4 (see [23, 24]).
The normalized forward sensitivity index ofR0 that differentiable with respect to a given parameter Q is defined as
(54)τQR0=∂R0∂Q×QR0,for Q describes the basic parameters.For example, the sensitivity index ofR01 with respect to basic parameter β0m is obtained as
(55)τβ0mR01=∂R01∂β0m×β0mR01,=12β0hΨβ0mΦ0αhαm/μhμm2μh+δ+γhαh+μhαm+μm×β0hΦ0Ψαhαmμhμm2μh+δ+γhαh+μhαm+μm×β0mR01=12>0.Using a similar procedure for the remaining parameters,τβ0mR01,τΨR01,τΦ0R01,τμhR01,τμmR01, τδR01,τγhR01ταhR01,ταmR01 are obtained, and their sensitivity indices are shown in Table 2 as follows.Table 2
Sensitivity indices of parameters.
Parameter symbolSensitivity indexαh0.002346αm0.248098Ψ0.5Φ00.5β0h0.5β0m0.5μm-0.023650μh-0.540314δ-0.485258γh-0.026464The sensitivity index ofR02 with respect to parameters is obtained using a similar method, and their sensitivity indices are written in Table 3.Table 3
Sensitivity indices of parameters.
Parameter symbolSensitivity indexαh0.020346αm0.238096Ψ0.5Φ00.082438β0h0.321659β0m0.385000Φ1m0.068169β1h0.208341β1m0.214286μm-0.013650μh-0.500314δ-0.476258γh-0.026464


### 4.1. Expression of the Sensitivity Indices

The sensitivity analysis description of a basic reproductive numberR01 with respect to ten basic parameters was stated in the Table 2. We can conclude that the parameters αm,αh,Ψ,Φ0,β0h, and β0m with positive sensitivity indices will increase the value of R01 if their values increase while the other parameters remain constant. This demonstrates that the population’s infection rate is increasing as a result of secondary cases. In contrast, the basic parameters μm,μh,δ, and γh with negative indices decrease the value of R01 if their value increases while retaining the constant of rest parameters. For example, τβ0mR01==0.5 shows that decreasing (increasing) the mosquito contact rate by 10% reduces (increases) the R01 by 5%. Similarly, τμmR01=−0.023650 and indicates that decreasing (increasing) the mosquito death rate by 10% increases (decreases) the basic reproduction number R01 by 0.23650%. Also, the sensitivity analysis expression of the basic reproductive number R02 with respect to thirteen basic parameters was shown in Table 3. As a result, the parameters αm,αh,Ψ,Φ0,Φ1m,β0hβ1h,β0m, and β1m have positive sensitivity indices. They will have a large influence on the spread of the malaria disease in the population. Because the transmission of malaria disease increases as their value increases, so does the number of infections as secondary cases in the population. The sensitivity indices for the parameters μm,μh,δ, and γh are all negative. As their values increase, those parameters have the potential to reduce malaria disease transmission while remaining constant. For example, τβ0mR02=0.385 indicates that decreasing (increasing) the mosquito contact rate by 10% decreases (increases) the R02 by 3.85%. In contrast, τμmR02=−0.023650 demonstrates that decreasing (increasing) the mosquito death rate by 10% increases (decreases) the basic reproduction number R02 by 0.23650%.


## 4.1. Expression of the Sensitivity Indices

The sensitivity analysis description of a basic reproductive numberR01 with respect to ten basic parameters was stated in the Table 2. We can conclude that the parameters αm,αh,Ψ,Φ0,β0h, and β0m with positive sensitivity indices will increase the value of R01 if their values increase while the other parameters remain constant. This demonstrates that the population’s infection rate is increasing as a result of secondary cases. In contrast, the basic parameters μm,μh,δ, and γh with negative indices decrease the value of R01 if their value increases while retaining the constant of rest parameters. For example, τβ0mR01==0.5 shows that decreasing (increasing) the mosquito contact rate by 10% reduces (increases) the R01 by 5%. Similarly, τμmR01=−0.023650 and indicates that decreasing (increasing) the mosquito death rate by 10% increases (decreases) the basic reproduction number R01 by 0.23650%. Also, the sensitivity analysis expression of the basic reproductive number R02 with respect to thirteen basic parameters was shown in Table 3. As a result, the parameters αm,αh,Ψ,Φ0,Φ1m,β0hβ1h,β0m, and β1m have positive sensitivity indices. They will have a large influence on the spread of the malaria disease in the population. Because the transmission of malaria disease increases as their value increases, so does the number of infections as secondary cases in the population. The sensitivity indices for the parameters μm,μh,δ, and γh are all negative. As their values increase, those parameters have the potential to reduce malaria disease transmission while remaining constant. For example, τβ0mR02=0.385 indicates that decreasing (increasing) the mosquito contact rate by 10% decreases (increases) the R02 by 3.85%. In contrast, τμmR02=−0.023650 demonstrates that decreasing (increasing) the mosquito death rate by 10% increases (decreases) the basic reproduction number R02 by 0.23650%.


## 5. Optimal Control Model

In this section, we provide a thorough qualitative examination of the time-dependent malaria dynamic model (1). The Pontryagin’s maximum principle [28] is used to describe this analysis, with the goal of lowering the exposed human Eht, infected human Iht, infected mosquito Imt, and control costs uit. In the case of time-dependent control, we employ Pontryagin’s maximum principle to derive the necessary conditions for disease control. After incorporating the controls into the malaria transmission model (1), the optimal control problem is as follows:
(56)dShdt=Ψ−1−u1βhT,R1+z∗ImShIm−μhSh+ωhRh,dEhdt=1−u1βhT,R1+z∗ImShIm−μh+αhEh,dIhdt=αhEh−μh+δ+γh+u2Ih,dRhdt=γh+u2Ih−μh+ωhRh,dSmdt=ΦT,R−1−u1βmT,R1+m∗IhSmIh−μm+u3Sm,dEmdt=1−u1βmT,R1+m∗IhSmIh−μm+αm+u3Em,dImdt=αmEm−μm+u3Im,dTdt=r1−TTmaxT−T0,dRdt=mdTdt⇒Rt=mTt−T0+ε,where u1t represents protective control using treated bed net, u2t represents antimalaria drug treatment control, and u3t denotes mosquito control via indoor residual spraying. The objective functional (57) of the optimal control problem (56) is given by
(57)Ju1,u2,u3=min⏟u1,u2,u3∫0tfD1Eh+D2Ih+D3Im+12B1u12+B2u22+B3u32dt,where tf represents the final time of control implementation, and quantities D1,D2, and D3 are weights constants of the exposed human population, infected human population, and infected mosquito population, respectively, while B1,B2, and B3 are weight constants used for treated bed net, treatment control using antimalaria drugs, and indoor residual spraying. In this paper, we make the expression 1/2Biui2 for the cost-control functions to be quadratic in order to obtain a unique optimal expression for each of the control variables from the optimality condition as in other studies [29–34]. The main goal is to achieve the best control triple u1∗,u2∗, and u3∗ so that
(58)Ju1∗,u2∗,u3∗=minJu1,u2,u3:u1,u2,u3∈U,where U=u1,u2,u3:uit such that u1,u2, and u3 are Lebesgue measurable on t∈0,tf with 0≤uit≤1, for i=1,2,3 is the control set. We express the Hamiltonian (H), which consists of the state equations (56) and integrand of the objective functional (57), by using the Pontryagin’s maximum principle [28] to obtain the necessary condition.
(59)H=D1Eh+D2Ih+D3Im+12∑i=13Biui2+λidShdt+dEhdt+dIhdt+dRhdt+dSmdt+dEmdt+dImdt+dTdt+dRdt,i=1,⋯,9.By plugging the given equations (56) and (57) into a minimizing problem of a Hamiltonian function (H) with respect to controls. The Hamiltonian (59) becomes
(60)H=B1Eh+B2Ih+B3Im+12D1u12+D2u22+D3u32+λ1Ψ−1−u1βhT,R1+z∗ImShIm−μhSh+ωhRh+λ21−u1βhT,R1+z∗ImShIm−μh+αhEh+λ3αhEh−μh+δ+γh+u2Ih+λ4γh+u2Ih−μh+ωhRh+λ5ΦT,R−1−u1βmT,R1+m∗IhSmIh−μm+u3Sm+λ61−u1βmT,R1+m∗IhSmIh−μm+αm+u3Em+λ7αmEm−μm+u3+λ8r1−TTmaxT−T0+λ9mr1−TTmaxT−T0,where λi for i=1,⋯,9 is adjoint variables. Next, we present the adjoint system and control characterizations using Pontryagin’s maximum principle [30], in conjunction with the existence of the optimal control problem [35], the following result can be obtained.Theorem 5.
Suppose we have an optimal control setu∗=u1∗,u2∗,u3∗ and Sh∗,Eh∗,Ih∗,Rh∗,Sm∗,Em∗,Im∗,T∗,R∗ solutions of the respective state system (56) that minimize Ju1,u2,u3 over U; then, there exist costate variables λ1,λ2,λ3,λ4,λ5,λ6,λ7,λ8, and λ9 such that
(61)dλ1dt=−1−u1βhT,R1+z∗ImImλ2−λ1+μhλ1,dλ2dt=−D1+λ2μh+αh−αhλ3,dλ3dt=−D2−1−u1βmT,R1+m∗IhSmλ6−λ5+λ3μh+δ+γh+u2−λ4u2+γh,dλ4dt=−ωhλ1+λ4μh+ωh,dλ5dt=−1−u1βmT,R1+m∗IhIhλ6−λ5+λ5u3+μm,dλ6dt=λ6u3+μm+αm−λ7αm,dλ7dt=−D3−1−u1βhT,R1+z∗ImShλ2−λ1+λ7u3+μm,dλ8dt=−rλ8−rλ8T0Tmax−2TTmax,dλ9dt=−mrλ9−mrλ9T0Tmax−2TTmax.
with full time (transversality) conditions(62)λ1tf=λ2tf=λ3tf=λ4tf=λ5tf=λ6tf=λ7tf=λ8tf=λ9tf=0.
Moreover, the optimal controlsu1∗,u2∗,u3∗ with the optimality condition are represented by
(63)u1∗=max0,min1,λ2−λ1βhT,R/1+z∗ImSh∗Im∗+λ6−λ5βmT,R/1+m∗IhSm∗Ih∗B1,(64)u2∗=max0,min1,λ3−λ4Ih∗B2,(65)u3∗=max0,min1,λ5Sm∗+λ6Em∗+λ7Im∗B3.Proof.
Differentiate the Hamiltonian function H (??) with respect to the state variablesSh,Eh,Ih,Rh,Sm,Em,Im,T, and R, respectively, the adjoint equations stated at the system (61) are obtained. Thus, the adjoint equations are given by
(66)dλ1dt=−∂H∂Sh=−1−u1βhT,R1+z∗ImImλ2−λ1+μhλ1,dλ2dt=−∂H∂Eh=−D1+λ2μh+αh−αhλ3,dλ3dt=−∂H∂Ih=−D2−1−u1βmT,R1+m∗IhSmλ6−λ5+λ3μh+δ+γh+u2−λ4u2+γh,dλ4dt=−∂H∂Rh=−ωhλ1+λ4μh+ωh,dλ5dt=−∂H∂Sm=−1−u1βmT,R1+m∗IhIhλ6−λ5+λ5u3+μm,dλ6dt=−∂H∂Em=λ6u3+μm+αm−λ7αm,dλ7dt=−∂H∂Im=−D3−1−u1βhT,R1+z∗ImShλ2−λ1+λ7u3+μm,dλ8dt=−∂H∂T=−rλ8−rλ8T0Tmax−2TTmax,dλ9dt=−∂H∂R=−mrλ9−mrλ9T0Tmax−2TTmax,with transversality conditions
(67)λ1tf=λ2tf=λ3tf=λ4tf=λ5tf=λ6tf=λ7tf=λ8tf=λ9tf=0.
Further, the optimality conditions with respect to the controls are given by(68)∂H∂ui=0,fori=1,2,3.
Using optimality condition (68), the control variables are obtained as
(69)u1∗=λ2−λ1βhT,R/1+z∗ImSh∗Im∗+λ6−λ5βmT,R/1+m∗IhSm∗Ih∗B1,u2∗=λ3−λ4Ih∗B2,u3∗=λ5Sm∗+λ6Em∗+λ7Im∗B3.
Finally, the compact form with boundary condition, the controls in equation (69), becomes
(70)u1∗=max0,min1,λ2−λ1βhT,R/1+z∗ImSh∗Im∗+λ6−λ5βmT,R/1+m∗IhSm∗Ih∗B1,(71)u2∗=max0,min1,λ3−λ4Ih∗B2,(72)u3∗=max0,min1,λ5Sm∗+λ6Em∗+λ7Im∗B3.
In the next part, we will see the numerical simulation of optimality system to identify an optimum and least cost strategy for controlling the transmission dynamics of malaria.


## 6. Numerical Simulation

A numerical result is required in this section to validate the model’s analytical result. We use a numerical method to confirm the theoretical results obtained in our model’s optimal control model analysis. The optimality system is composed of two systems, the state system (56) and the adjoint system (61), each with its own initial and final-time conditions, with the control value state in equation (62). The iterative method known as the fourth forward-backward Runge-Kutta sweep is used to solve this optimality system. The state equation (56) is solved with the initial values of state variables using the fourth-order forward Runge-Kutta method. We used backward fourth order Runge Kutta to solve the adjoint equations once we had the solution of the state functions and the value of optimal controls. The controls are then updated by a convex combination of the previous controls and the value from the optimality condition (62). This process is repeated, and iterations are completed if the values from previous iterations are close to the values from current iterations [36]. The optimal control problem is simulated with the parameter values given in Table 4 when the basic reproduction numbers are R01=24.8 and R02=56.28. The initial conditions that we used for simulation of the optimal control are Sh0=120,Eh0=20,Ih0=10,Rh0=5,Sm0=400,Em0=60,Im0=30,T0=16.1oC, and R0=11 mm. The weight value constants of the state and controls that we used are given as D1=80,D2=60,D3=100,B1=60,B2=100, and B3=80. To determine the impact of each controls on the reduction of malaria, we used the following four strategies with different combination of two controls at a time and three controls at a time.Table 4
Parameter values from literature used for the model (1).
ParametersParameter descriptionEstimated valueReferenceμhHuman death rate due to natural cause0.00004[23]δInduced death rate of human0.068[23]γhRecover rate of infected humans0.0035[20]ΨHuman recruitment rate0.071[20]Φ0Mosquito recruitment rate0.041[24]μmMosquito natural death rate0.05[20]ωhHuman immune deficiency rate0.09[16]β1mMosquito contact rate (incremental)0.07[17]β1hHuman contact rate (incremental)0.05[16]Φ1mIncremental birth rate of mosquito0.09[17]β0hContact rate of human with mosquito0.03[23]β0mContact rate of mosquito with human0.04[14]αhExposed human progression rate0.058[24]αmExposed mosquito progression rate0.055[24]mTemperature dependent rate of precipitation0.08[26]z∗Percentage of an antibody produced by a human0.06[26]m∗Percentage of an antibody produced by a mosquito0.04[26]rTemperature increase rate0.007[11]T0The minimum temperature16°C[27]TmaxThe maximum temperature28°C[6]


### 6.1. Strategy A

Applying the treated bed netu1 and treatment of infected with drugs u2.In this strategy, we used a combination of two controls: protective by treated bed netu1 and treatment for infected humans with antimalaria drugs u2 to minimize the objective functional (57), while the control indoor residual spraying u3 is set to zero. According to Figures 2(a) and 2(b), if there are controls, the total population of exposed and infected humans is lower than if there are no controls, and it tends to be at its lowest value at the end of the intervention. Then, as shown in Figure 2(c), the number of infected mosquitos in the presence of controls is decreasing. In contrast, in the absence of controls, the number of infected mosquitos increases. In order to reduce the number of exposed humans, infected humans, infected mosquito populations, and the associated costs, corresponding with the controls u1 and u2, the control profiles shown in Figure 2(d) indicate that use of the treated bed net u1 maintained its maximum level bound 100 percent for 180 days entire time of strategy and treatment; the infected with antimalaria drugs u2 is intended upper bound 100 percent for 6 days before gradually declining to zero at the end time.Figure 2
Simulations of model with strategy of treated bed netu1 and treatment u2.
(a)(b)(c)(d)


### 6.2. Strategy B

Applying the treated bed netu1 and indoor residual spraying u3.To reduce the objective functional (57), a combination of two controls protective using treated bed net u1 and an indoor residual spraying u3 is implemented, while the control treatment for the infected with antimalaria drugs u2 is set to zero. Figures 3(a)–3(c) show that using controls reduces the total population of exposed humans, infected humans, and infected mosquitos faster than not using controls. At the end of the intervention, the total number of exposed humans Eh and infected mosquitoes Im is at its lowest. The profile of the control in Figure 3(d) indicates that the treated bed net u1 is maintained, it is upper bound at 100% for 178 days before gradually declining to its lower bound and an indoor residual spraying, and u3 maintains its high level at 100% until 18 days and gradually decreases to its lower bound in the final time.Figure 3
Simulations with strategy of treated bed netu1 and indoor spraying u3.
(a)(b)(c)(d)


### 6.3. Strategy C

Applying the treatment of infectedu2 and indoor spraying u3.The objective functional (57) is optimized in this strategy by combining two controls treatment of infected with anti-malaria drugs u2 and an indoor residual spraying u3, whereas the control treated bed net u1 is set to zero. Figures 4(a) and 4(b) show that in the absence of controls, the significant number of humans exposed and infected is higher than in the presence of controls. At the end of the intervention, the total population of exposed humans Eh and infected humans Ih has reached its lowest point. Figure 4(c) shows that infected mosquitos in the presence of controls are smaller than those in the absence of controls and fall to the minimum value at the end of the intervention. Similarly, in Figure 4(d), the control profiles show that both antimalaria drug treatment u2 and indoor residual spraying u3 maintained their upper bounds of 99 percent for 30 days and 60 days, respectively, before gradually decreasing to their minimum at the end of the intervention.Figure 4
Simulations model with strategy of treatmentu2 and indoor spraying u3.
(a)(b)(c)(d)


### 6.4. Strategy D

Using the treated bed netu1, treatment u2 and indoor spraying u3.To minimize the objective functional (57), we used all controls protective using treated bed net u1, treatment for the infected human with antimalaria drugs u2 and indoor residual spraying u3. It can be seen in Figures 5(a)–5(c) that the controls indicate that the total population of exposed human, infected human, and infected mosquito populations are increasing in the absence of controls while decreasing in the presence of controls. The total number of exposed humans Eh, infected humans Ih, and infected mosquitos Im is reduced to their lowest at times t=100,t=46, and t=10, respectively. The control profiles in Figure 5(d) show that with this strategy, the use of the treated bed net u1 nearly maintained its maximum level 100, and the treatment of an infected human with antimalaria drugs u2 maintains 100% coverage for 6 days before gradually decreasing to zero, whereas the control, an indoor residual spraying u3, maintains 100% coverage for 16 days before gradually decreasing to zero at the end of the final time.Figure 5
Simulations with strategy of treated bed netu1, treatment u2, and indoor spraying u3.
(a)(b)(c)(d)Figure6 illustrates the diagram of bifurcation for the case when T=T0 and T=Tmax for the malaria model problem. From this figure, we have seen that the model exhibits forward and backward bifurcation, respectively. The biological meaning of this indicates that in equation (21) if R02<1 then automatically implies that R01<1 and DFE exist for both R01 and R02. However, R01<1 does not automatically imply R02<1, since the value of R02 may be greater than 1, which implies that R02 may present backward bifurcation while R01 only present forward bifurcation.Figure 6
Bifurcation diagrams whenT=T0 and T=Tmax for the malaria model problem.


## 6.1. Strategy A

Applying the treated bed netu1 and treatment of infected with drugs u2.In this strategy, we used a combination of two controls: protective by treated bed netu1 and treatment for infected humans with antimalaria drugs u2 to minimize the objective functional (57), while the control indoor residual spraying u3 is set to zero. According to Figures 2(a) and 2(b), if there are controls, the total population of exposed and infected humans is lower than if there are no controls, and it tends to be at its lowest value at the end of the intervention. Then, as shown in Figure 2(c), the number of infected mosquitos in the presence of controls is decreasing. In contrast, in the absence of controls, the number of infected mosquitos increases. In order to reduce the number of exposed humans, infected humans, infected mosquito populations, and the associated costs, corresponding with the controls u1 and u2, the control profiles shown in Figure 2(d) indicate that use of the treated bed net u1 maintained its maximum level bound 100 percent for 180 days entire time of strategy and treatment; the infected with antimalaria drugs u2 is intended upper bound 100 percent for 6 days before gradually declining to zero at the end time.Figure 2
Simulations of model with strategy of treated bed netu1 and treatment u2.
(a)(b)(c)(d)


## 6.2. Strategy B

Applying the treated bed netu1 and indoor residual spraying u3.To reduce the objective functional (57), a combination of two controls protective using treated bed net u1 and an indoor residual spraying u3 is implemented, while the control treatment for the infected with antimalaria drugs u2 is set to zero. Figures 3(a)–3(c) show that using controls reduces the total population of exposed humans, infected humans, and infected mosquitos faster than not using controls. At the end of the intervention, the total number of exposed humans Eh and infected mosquitoes Im is at its lowest. The profile of the control in Figure 3(d) indicates that the treated bed net u1 is maintained, it is upper bound at 100% for 178 days before gradually declining to its lower bound and an indoor residual spraying, and u3 maintains its high level at 100% until 18 days and gradually decreases to its lower bound in the final time.Figure 3
Simulations with strategy of treated bed netu1 and indoor spraying u3.
(a)(b)(c)(d)


## 6.3. Strategy C

Applying the treatment of infectedu2 and indoor spraying u3.The objective functional (57) is optimized in this strategy by combining two controls treatment of infected with anti-malaria drugs u2 and an indoor residual spraying u3, whereas the control treated bed net u1 is set to zero. Figures 4(a) and 4(b) show that in the absence of controls, the significant number of humans exposed and infected is higher than in the presence of controls. At the end of the intervention, the total population of exposed humans Eh and infected humans Ih has reached its lowest point. Figure 4(c) shows that infected mosquitos in the presence of controls are smaller than those in the absence of controls and fall to the minimum value at the end of the intervention. Similarly, in Figure 4(d), the control profiles show that both antimalaria drug treatment u2 and indoor residual spraying u3 maintained their upper bounds of 99 percent for 30 days and 60 days, respectively, before gradually decreasing to their minimum at the end of the intervention.Figure 4
Simulations model with strategy of treatmentu2 and indoor spraying u3.
(a)(b)(c)(d)


## 6.4. Strategy D

Using the treated bed netu1, treatment u2 and indoor spraying u3.To minimize the objective functional (57), we used all controls protective using treated bed net u1, treatment for the infected human with antimalaria drugs u2 and indoor residual spraying u3. It can be seen in Figures 5(a)–5(c) that the controls indicate that the total population of exposed human, infected human, and infected mosquito populations are increasing in the absence of controls while decreasing in the presence of controls. The total number of exposed humans Eh, infected humans Ih, and infected mosquitos Im is reduced to their lowest at times t=100,t=46, and t=10, respectively. The control profiles in Figure 5(d) show that with this strategy, the use of the treated bed net u1 nearly maintained its maximum level 100, and the treatment of an infected human with antimalaria drugs u2 maintains 100% coverage for 6 days before gradually decreasing to zero, whereas the control, an indoor residual spraying u3, maintains 100% coverage for 16 days before gradually decreasing to zero at the end of the final time.Figure 5
Simulations with strategy of treated bed netu1, treatment u2, and indoor spraying u3.
(a)(b)(c)(d)Figure6 illustrates the diagram of bifurcation for the case when T=T0 and T=Tmax for the malaria model problem. From this figure, we have seen that the model exhibits forward and backward bifurcation, respectively. The biological meaning of this indicates that in equation (21) if R02<1 then automatically implies that R01<1 and DFE exist for both R01 and R02. However, R01<1 does not automatically imply R02<1, since the value of R02 may be greater than 1, which implies that R02 may present backward bifurcation while R01 only present forward bifurcation.Figure 6
Bifurcation diagrams whenT=T0 and T=Tmax for the malaria model problem.


## 7. Cost Effectiveness Analysis

It is critical to identify the most optimal and least expensive strategy, along with a cost-effectiveness analysis, in order to design a method of eradicating malaria transmission. A method known as the incremental cost-effectiveness ratio (ICER) is used to develop such a strategy. Basically, for two strategies say 1 and 2, the ICER will be computed as [37]
(73)Differenceinavertedcostsbetweentwostrategies1and2Differenceininfectionsavertedbewteentwostrategies1and2.In this section, the cost-effectiveness of two control strategies, including treated bed nets, antimalaria drug treatment of infected people, and indoor residual spraying, has been calculated. The total number of people averted is calculated by subtracting the total number of new cases of malaria infection with control from the total number of new cases of malaria without control, and the total cost of each strategy is calculated using their respective cost functions1/2B1u12,1/2B2u22, and 1/2B3u32, to calculate over the time of intervention [38]. We use the outcome of a numerical simulation from the malaria dynamics model (56) to compute the ICER. However, using only one strategy to prevent the spread of malaria disease in a community may not be effective. Then also from the numerical simulations of the model (56) with the parameters in Table 4, the computed total number of people saved and the total cost in increasing order were shown in Table 5.Table 5
Total amount of infection averted and total cost for all strategies.
StrategiesTotal infections avertedTotal cost ($)Strategy B3754.407012.84Strategy C4536.249098.19Strategy A4673.348306.89Strategy D4676.629154.62After obtaining the total amount of people averted and total cost of each strategy as given in Table5 to compare two intervention strategies, the incremental cost effectiveness ratio (ICER) for the each competing strategies is estimated as
(74)ICERB=7012.843754.40=1.86,(75)ICERC=9098.19−7012.844536.24−3754.40=266.72,(76)ICERA=8306.89−99098.194673.34−4536.24=−5.77,(77)ICERD=9154.62−8306.894676.62−4673.34=258.45.The number of people averted in strategy B, C, A, and D in an increasing rank is given in Table6.Table 6
Total amount of the infection averted and total cost with their ICER.
StrategiesTotal infections avertedTotal cost ($)ICERStrategy B3754.407012.841.86Strategy C4536.249098.19266.72Strategy A4673.348306.89-5.77Strategy D4676.629154.62258.45We can see that the ICER(B) is less than the ICER(C) from the strategies B and C in Table6. This implies that strategy B dominates strategy C. This implies that strategy C is more expensive than strategy C. As a result, we removed C from the strategies. Then, as shown in Table 7, recalculate the ICER for the remaining competing strategies B, A, and D.Table 7
Total amount of the infection averted and total cost with their ICER.
StrategiesTotal infections avertedTotal cost ($)ICERStrategy B3754.407012.841.86Strategy A4673.348306.891.40Strategy D4676.629154.62258.45The competition between the interventions B and A was depicted in Table7. As can be seen, ICER(A) is less than ICER(B). This demonstrates that strategy B was dominated by strategy A. As a result, strategy A is more efficient and less expensive than strategy B. As a result, we removed strategy B from the list of competing strategies and recalculated the ICER, as shown in Table 8.Table 8
Total number of infections averted and total cost with their ICER.
StrategiesTotal infections avertedTotal cost ($)ICERStrategy A4673.348306.891.78Strategy D4676.629154.62258.45Table8 compares strategies A and D. This predicts that ICER(D) is greater than ICER(A). This implies that strategy D is dominated by strategy A. As a result, strategy A is the least expensive and most optimal. As a result, we conclude that strategy A, which consists of a combination of treated bed nets and antimalarial drug treatment of infected humans, is the most optimal and least expensive strategy for preventing the spread of malaria transmission.


## 8. Conclusion

We proposed and developed a deterministic malaria transmission dynamic model with climate variability and nonlinear incidence in this paper. Malaria spread is influenced by climate change. The model analysis reveals that it is bounded, epidemiologically meaningful, and mathematically well-posed in a specific domain. Using the next-generation matrix method, we determined the basic reproduction number in relation to the disease-free equilibrium. The Routh-Hurwitz criterion is used to determine local stability of equilibria points, while the Lyapunov function is used to determine global stability. According to the model analysis, if the basic reproduction number is less than one, the disease-free equilibrium is locally and globally asymptotically stable, whereas if the basic reproduction number is greater than one, the unique endemic equilibrium exists. In addition, the sensitivity analysis of the basic reproduction number with respect to all parameters was obtained. The model has a forward and backward bifurcation. Based on the findings, we concluded that reducing human-mosquito contact, increasing mosquito mortality, and increasing the treated rate of infected humans can all help to reduce the malaria burden in the community. Furthermore, the malaria transmission model is extended to an optimal control problem by incorporating three continuous controls, namely, personal protection with treated bed nets, treatment of infected with antimalaria drugs, and indoor residual spraying for vector killing strategy. The maximum principle of Pontryagin is used to obtain the necessary condition of the optimal control problem. The cost-effectiveness of each control combination is then investigated. As a result, based on the optimality system simulation results and cost-effectiveness analysis, we proposed that the combination of treated bed nets, and treatment is the most optimal and least expensive strategy to reduce malaria transmission.



---

*Source: 9667396-2022-12-02.xml*