# Modeling the Impact of Bed-Net Use and Treatment on Malaria Transmission Dynamics

Authors: Bello Gimba; Saminu Iliyasu Bala
Journal: International Scholarly Research Notices
 (2017)

---

## Abstract

We modeled the impact of bed-net use and insecticide treated nets (ITNs), temperature, and treatment on malaria transmission dynamics using ordinary differential equations. To achieve this we formulated a simple model of mosquito biting rate that depends on temperature and usage of insecticides treated bed nets. We conducted global uncertainty and sensitivity analysis using Latin Hypercube Sampling (LHC) and Partial Rank Correlation Coefficient (PRCC) in order to find the most effective parameters that affect malaria transmission dynamics. We established the existence of the region where the model is epidemiologically feasible. We conducted the stability analysis of the disease-free equilibrium by the threshold parameter. We found the condition for the existence of the endemic equilibrium and provided necessary condition for its stability. Our results show that the peak of mosquitoes biting rate occurs at a range of temperature values not on a single value as previously reported in literature. The results also show that the combination of treatment and ITNs usage is the most effective intervention strategy towards control and eradication of malaria transmissions. Sensitivity analysis results indicate that the biting rate and the mosquitoes death rates are the most important parameters in the dynamics of malaria transmission.


---

## Body

## 1. Introduction

Malaria is one of the most devastating infectious diseases in the world and is caused byPlasmodium parasite, which is transmitted via the bites of infected mosquitoes. Among the high risk groups are pregnant women, nonimmune travelers, and children [1]. In pregnant women, malaria has adverse effect on birth outcome which includes low birth weight, abortion, and still born [1]. Apart from health related problems, it also imposes huge socioeconomic burden in malaria-endemic nations. As discussed in Forouzannia and Gumel [2] the annual economic burden of malaria in Africa alone was estimated to be around US $8 billion. These necessitated the formation of several intervention strategies in many countries to mitigate the impact of malaria disease. This includes the use of insecticide treated bed nets (ITNs), intermittent preventive treatment (IPT) especially, for pregnant women during antenatal period, reducing mosquitoes population through the destruction of breeding sites or killing of the larva stage at breading sites that cannot be destroyed [3–5]. Other interventions strategies are the use of indoor residual spraying (IRS) in killing infected mosquitoes resting indoors after blood meal and the use of sterile insect technique [6]. Despite broad efforts for eradication, malaria remains a significant problem resulting in the death of millions of people [2, 7–9]. Most malaria cases and deaths occur in sub-Saharan Africa with Nigeria and Democratic Republic of Congo accounting for about 40% of malaria mortality worldwide [9]. There are a number of characteristics of malaria disease that complicates control efforts. Typical among them is clinical immunity, which is a situation where protection against the clinical symptoms of the disease is developed despite the presence of the parasites [7, 10, 11]. Others include seasonality [12, 13] and treatment failure that might occur due to wrong dosage of medication; see [2] and the references therein. In the context of malaria transmission, seasonality encapsulates complex phenomenon whose definition varies in many studies. Temperature variations have been reported by many to play significant role in the dynamics of malaria transmissions. For example, the report of Roll Back Malaria 2015 indicates that a rise in temperature by 2-3°C will increase the number of people at climatic risk of malaria by 3–5%. Furthermore, the abundance of mosquitoes and the transmission risk have been reported to be influenced by temperature [13–16]. At high temperature, studies have indicated that people are unlikely to use ITNs much [17].Mathematical models of malaria transmission have been developed by several researchers to gain insight into the dynamics of the disease transmission so as to contribute towards its eradication. Some of these models can be found in [18–23]. These models are in varying degree of complexity. For example, the model of [2] is an age-structured model with several compartments. The model of [7] is made up of four compartments comprising of only the human population. In the model of [5, 17] the authors introduced explicit equation for the proportion of ITNs use as a function of mosquito biting rates.Some of the problems that complicate malaria control include (1) the presence of individuals who are clinically immune to the disease but can transmit it through bite from susceptible mosquitoes and (2) hot weather which can lead to reduction in the use of ITN. For these reasons, it is important to study the qualitative impact of treatment, immunity, and seasonality on the dynamics of malaria transmission. In this work, we present a vector-host model of malaria transmission dynamics of immune and nonimmune human populations that accounts for the impact of ITN usage and seasonality on the disease. We propose a model of mosquito biting rate as a nonlinear function of temperature and ITN usage to mimic seasonality. This will help in devising optimal intervention strategies that will offer more realistic predictions to control malaria spread. To the best of our knowledge, this is the first vector-host mathematical model for malaria transmission, which explores the impact of daily temperature variations and ITNs usage on control of malaria transmissions. The current study extends the work of [7] by designing a vector-host model for malaria transmission dynamics. The study also extends the work of [17] by modeling the mosquito biting rate as a function of temperature to mimic seasonality. The paper is organized as follows. We formulated the model in Section 2 and analyzed it qualitatively in Section 3, in Section 4 we conducted global uncertainty and sensitivity analysis, Section 5 is the discussion part, and in Section 6 we present our conclusions.


## 2. Model Formulation

In this section, we modify an existing mathematical model for malaria transmission dynamics developed by [7]. The process of the modification is presented below. Following [7] we defined naive individuals as those who have never been infected with malaria, or those who have been infected but have not developed clinical immunity, or those who have lost all immunity. Similarly, clinically immune individuals are those with immunity to clinical symptoms. The total human population denoted by Nh is divided into mutually exclusive subpopulations of susceptible naive Sn, susceptible clinically immune Sc, infected naive In, and infected clinically immune Ic, so that Nh=Sn+Sc+In+Ic.The total mosquitoes population denoted byNm is divided into compartments of susceptible and infected mosquitoes, so that Nm=Sv+Iv.All recruitment is assumed to be into the susceptible naive human population generated via birth and/or immigration at a rateλhNh. The population of naive susceptible individuals (Sn) is increased by naive infected individuals that recovered without immunity at a rate γNN, treated naive individuals that recovered without immunity at a rate ϵnn, and clinically immune individuals that lost immunity at a rate αc. The population of naive susceptible individuals (Sn) is decreased by natural death rate ϕh=μh2Nh+μh and force of infection (Λh), following effective contacts with infected mosquito. Here μh2 and μh represent the density dependent and density independent part of human death rate and emigration, respectively. We model the force of infection from mosquitoes to human as Λh=a(β(T))bIv/Nh.Hereb is the probability of infection of susceptible human per bite by an infected mosquito and a≡a(β(T)) is the biting rate of mosquitoes on susceptible human; β≡β(T) represents the proportion of ITN usage and depends on environmental temperature T. See (7) for the functional form of β. We assumed that temperature is a parameter that is time independent to make the analysis easier. Thus,(1)dSndt=λhNh-μh2Nh+μhSn+ϵnn+γNNIn+αcSc-ΛhSn.The clinically susceptible population is generated by the treated naive individuals that become clinically immune, infected naive individuals that recovered with clinical immunity at a rate γNC, infected clinically immune individuals that recovered with clinical immunity at a rate γC, and treated clinically immune infected individuals that recover at a rate Tc. It is decreased by susceptible clinically immune individuals that lose immunity at rate αc, natural death at rate ϕh, and the force of infection that pushed out susceptible clinically immune human into infected clinically immune population as a result of contact with infected mosquitoes at rate Λh. We assume that clinically infected individuals recover into the clinically susceptible compartment only. Thus,(2)dScdt=Tn-ϵnn+γNCIn-ΛhSc+Tc+γCIc-αc+μh2Nh+μhSc.The population of infected naive humans is generated by the population of the infectious susceptible naive human that become infected. It is decreased by the treated infected naive individuals at a rate Tn, infected naive individuals that recover without immunity at a rate γNN, infected naive individuals that recover with immunity at a rate γNC, the natural death ϕh, and disease induced death rate δN. Thus,(3)dIndt=ΛhSn-Tn+γNN+γNCIn-μh2Nh+μhIn-δNIn.The population of infected clinically immune human is generated by the population of susceptible clinically immune humans that become infected. It is decreased by the treated infected clinically immune individuals at a rate Tc, population of clinically immune individuals that recover with immunity at a rate γC, the natural death ϕh, and disease induced death at a rate δC. Thus,(4)dIcdt=ΛhSc-γCIc-μh2Nh+μhIc-TcIc-δCIc.The population of susceptible mosquitoes is generated by birth at a rate λv. It is reduced by natural death ϕv=μv2Nm+μv, contact with ITNs at the rate βμv3, and infection when in contact with infected human at rate Λv=Icτc+τnIn/Nh. The parameters μv and μv2 represent the density independent and density dependent parts of the mosquitoes death rate, respectively. Here τn=ac1, τc=ac2, and c1 and c2 are the probabilities that susceptible mosquitoes become infectious after biting an infected naive or clinically immune human, respectively. Thus,(5)dSvdt=λvNm-Icτc+τnInSvNh-μv2Nm+μv+βμv3Sv.The population of infected mosquitoes is increased by infected susceptible mosquitoes at rate Λv. It is decreased by the natural death ϕv and death when they come into contact with ITN βμv3. Thus,(6)dIvdt=Icτc+τnInSvNh-μv2Nm+μv+βμv3Iv.In the report of [17], the authors model the mosquitoes biting rate as a linear function of ITN usage while [24] considers a more general form. None of these authors consider the impact of temperature on bed-net use despite its significance. In this work we model the biting rate as(7)a=βmax-ββmax-βmin,where β=e-hT/(T-T0)2 represents the proportion of ITN usage and the parameters T0 and h are location and scale parameters measured in °C, respectively. The choice of temperature as a parameter in the biting rate is to mimic seasonality. The justification of this novel approach is due to many reports in literature on the relative importance of temperature in malaria transmission dynamics as outlined in the introduction. In Figure 1, we study impact of shape and scale parameters on the biting rate. From Figure 1, we observe that the optimum temperature for the biting rate is not a single temperature value as discussed in [25] but a range of values. In [15], the authors review some calibrated models of temperature variations in terms of mosquitoes biting rates. One of the findings is that biting rates are optimal at certain temperature values. Typical values reported are 24.4°C, 25.0°C, 26.3°C, and 27.5°C. In this work we are reporting a range of values that encapsulates individual values from several reports. From Figure 1(a), it can be seen that before the maximum biting rate is attained, high values of location parameter will predict relatively lower biting rates. This finding is in contrast to the result of increasing the scale parameter as depicted on Figure 1(b).Figure 1
(a) Effects of location parameter on the biting rate forT0=27°C. (b) Effects of scale parameter on the biting rate for h=0.1°C.
(a)
Effect of location parameter on the biting rate
(b)
Effect of scale parameter on the biting rateIt follows, based on the above derivations and assumptions, that the model for the transmission dynamics of malaria is given by the following deterministic system of nonlinear differential equations. Flow diagram of the model is depicted in Figure2, and the state variables and parameters of the model are described in Tables 1 and 2, respectively:(8)dSndt=λhNh-μh2Nh+μhSn+ϵnn+γNNIn+αcSc-abIvSnNh,dIndt=abIvSnNh-Tn+γNN+γNCIn-μh2Nh+μhIn-δNIn,dScdt=Tn-ϵnn+γNCIn-abIvScNh+Tc+γCIc-αc+μh2Nh+μhSc,dIcdt=abIvScNh-γCIc-μh2Nh+μhIc-TcIc-δCIc,dSvdt=λvNm-Icτc+τnInSvNh-μv2Nm+μv+βμv1Sv,dIvdt=Icτc+τnInSvNh-μv2Nm+μv+βμv1Iv,dNhdt=λhNh-μh2Nh+μhNh-δNIn-δCIc,dNmdt=λvNm-μv2Nm+μv+βμv1Nm.Table 1
Description of state variable of the model.
Variable
Description
Sn
Population of susceptible naive human
In
Population of infected naive human
Sc
Population of susceptible clinically immune human
Ic
Population of infected clinically immune human
Nh
Total human population
Sv
Population of susceptible mosquitoes
Iv
Population of infected mosquitoes
Nm
Total population of mosquitoesTable 2
Parameters and their descriptions.
Parameters
Description and dimension
λh
Recruitment rate into human population(humans×day-1)
λv
Recruitment rate into mosquitoes population (mosquitoes×day-1)
Tn
Treatment rate of infected naive humans (day-1)
Tc
Treatment rate of clinically infected immune humans (day-1)
c1
Probability of transmission of infection from infected naive human to susceptible mosquito
c2
Probability of transmission of infection from infected immune human to susceptible mosquito
ρ
Number of mosquitoes per human host
b
Probability of transmission of infection from an infected mosquito to susceptible human
T
Temperature (°C)
T0
Location parameter (°C)
βmax
Maximum biting rate per mosquito (day-1)
βmin
Minimum biting rate per mosquito (day-1)
h
Scale parameter (°C)
ϵnn
Rate at which treated naive individuals recovered without immunity (day-1)
γNN
Rate at which naive individuals recovered without immunity (day-1)
γNC
Rate at which naive individuals recovered with clinical immunity (day-1)
αc
Rate at which clinically susceptible individuals lose immunity (day-1)
δN
Disease induced death rate for a naive individual (day-1)
δC
Disease induced death rate for a clinically immune individual (day-1)
γC
Recovery rate for clinically immune individuals (day-1)
μh2
Density dependent part of the death and emigration rate for humans (humanday-1)
μh
Density independent part of the death rate for humans (humanday-1)
μv
Density independent part of the death rate for mosquitoes (mosquitoesday-1)
μv2
Density dependent part of the death rate for mosquitoes (day-1)
μv3
ITN induced death rate for mosquitoes (day-1)Figure 2
Susceptible humans,Sn, Sc. can be infected by infectious mosquitoes. They then pass to the respective infectious compartments, In, Ic, before reentering the susceptible classes again or die due to the disease. The susceptible mosquitoes, Sv, can become infected when they bite infectious humans. The infected mosquitoes then move to the infectious class Iv.


## 3. Model Analysis

### 3.1. Basic Properties of the Model

Lemma 1.
LetNh∗,Nm∗ be the equilibrium solutions of the total human and mosquito populations, respectively. The closed set D={Sn,In,Sc,Ic,Sv,Iv∈R+6:Nh≤Nh∗,Nm≤Nm∗} is positively invariant and attracting.Proof.
Adding the first four equations and the last two equations of model (8) we obtained(9)dNhdt=λhNh-μh2Nh+μhNh-δNIn-δCIc,dNmdt=λvNm-μv2Nm+μv+βμv1Nm.The mosquito population is modeled by logistic growth with carrying capacity kv=(λv-μv-μv1)/μv2. It is easy to see that dNh/dt≤λh-μhNh1-Nh/kh and dNm/dt≤λv-μvNm1-Nm/kv. This shows that dNh/dt≤0 if 1≥Nh/kh and it approaches kh. Similarly, dNm/dt≤0 if 1≥Nm/kv and it approaches kv. Hence, using comparison theory [2] (10)Nht≤khNh0Nh01-e-λh-μht+khe-λh-μht,Nmt≤kvNm0Nm01-e-λv-μvt+kve-λv-μvt.If Nh(0)≤kh, then Nh(t)≤kh and if Nm(0)≤kv, then Nm(t)≤kv. Thus, the region D is positively invariant for the model. Moreover, if Nh(0)≥kh, Nm(0)≥kv, then either the solution enters the region D in finite time or Nh(t)→kh, Nm(t)→kv, as t→∞. Thus the region attracts all solutions in R+6. Now that we have shown that D is positively invariant, the requirement for existence and uniqueness of solutions holds for the system [2]. The dynamics of the system in the region D will henceforth be investigated.


### 3.2. Scaling

To analyze the malaria model (8), we think it is easier to work with fractional population instead of actual populations by scaling the population of each class by the total species population. We let Sn=uNh, In=vNh, Sc=wNh, Ic=xNh, Iv=zNm, Sv=yNm. We arbitrarily scale the time variables by λv by introducing τ=tλv so that(11)dudτ=1λvNhdSn/dt-SndNh/dtNh2,dvdτ=1λvNhdIn/dt-IndNh/dtNh2and so on for the rest of the variables. In the absence of the disease, the human population follows logistic growth with carrying capacity kh=(λh-μh)/μh2. Following [22] we scaled the human and vector populations in the first 6 equations of model (8) using their respective carrying capacities as Nh=khNh∗, Nm=kvNm∗. So we have a new 6-dimensional system of equations with two additional dimensions for the two total population variables Nh and Nm as(12)dudτ=-a1bkvNm∗khNh∗z-λu+αw+vγnn+λ,dvdτ=-Tj+γnc1+γnn+λ+δnv+a1bkvNm∗khNh∗uz,dwdτ=-wa1bkvNm∗khNh∗z+α+λ+vTj-ϵn+γnc1+xTi+γc,dxdτ=a1bkvNm∗khNh∗wz-xλ+Ti+γc+δc,dydτ=-vytn-xytc-y+1,dzdτ=vytn+xytc-z,dNhdτ=1-Nh1-ω1Nh-Nhvδn+xδc,dNmdτ=1-Nm1-ω2Nm-Nmω3,where λ=λh/λv, a1=a/λv, b1=b/λv, α=αc/λv, γnn=γNN/λv, Ti=Tc/λv, Tj=Tn/λv, ϵn=ϵnn/λv, δc=δC/λv, δn=δN/λv, γnc1=γNC/λv, γc=γC/λv, tn=τn/λv, tc=τc/λv, ω1=μh/λv, ω2=μv/λv, and ω3=μv3/λv.The analytic solution of model (12) will be more tractable if the number of equations is reduced. To this end, we let u=1-v-w-x, y=1-z, and by substituting these in (12), two equations are eliminated leaving us with the reduced model:(13)dvdτ=1-v-w-xξz-Av,dwdτ=-ξz+Fw+mx+vγnc,dxdτ=-Bx+ξwz,dzdτ=-vtn+xtc+1z+vtn+xtc,where ξ=a1bkvρ/kh, ρ=Nm∗/Nh∗, A=Tj+γnc1+γnn+λ+δn, m=Ti+γc, B=Ti+γc+λ+δc, F=α+λ, γnc=Tj-ϵn+γnc1.


### 3.3. Reproduction Number and Disease-Free Equilibrium

The disease-free equilibrium point (DFE) is the solution of the reduced model (13) when the disease classes are identically zero. By setting the right hand sides of (13) and the disease classes to zero, then solving the resulting equations simultaneously, we obtained the DFE of the reduced model (13) as Edfe=w,v,x,z=0,0,0,0. The nonnegative matrix for the new infection (ϝ) and the matrix for the transition terms (V) are given by(14)ϝ=00ξ000tntc0,V=A000B0001,ϝV-1=00ξ000tnAtcB0. The spectral radius is(15)ρϝV-1=ξtnA.We define the reproduction number as(16)Rτ=ξtnA.Lemma 2.
The equilibrium pointEdfe of the reduced model (13) is locally asymptotically stable (LAS) if Rτ<1 and unstable if Rτ>1.Proof.
The Jacobian matrix evaluated atEdfe is(17)JEdfe=-FγNCm00-A0ξ00-B00tntc-1.The characteristic equation is(18)F+ηB+ηηA+1+η2-tnξ+A=0.The eigenvalues are (19)η1=-F,η2=-B,η3=-A+1-A+12-4A1-Rτ2,η4=-A+1+A+12-4A1-Rτ2.Clearly, ηi, i=1,2,3, are negatives and that η4 is negative provided Rτ<1.


### 3.4. Endemic Equilibrium Point and Backward Bifurcation

Further analysis of the model is done by finding the endemic equilibrium solutionsx∗,v∗,w∗,z∗ of the reduced model (13). First we state the following.Theorem 3.
The reduced model has the following:(1)
A single endemic equilibrium solution ifRτ≥1, or R1=1, and P1Rτ>P3.(2)
Two endemic equilibrium solutions ifP1Rτ>P3 and R1=4P2P41-Rτ/P1Rτ-P32<1.(3)
No endemic equilibrium solution otherwise.Here P1=Aλtn+γNCtc, P2=Aλtn+γNCtc+λ+γNC, P4=BFtn2, and P3=tnBFtn+λA+BF+BγNC.Proof.
To prove this, we first expressed the right hand side of the reduced model (13) in terms of the equilibrium solutions x∗, v∗, w∗, z∗. We then eliminated the rest of the variables leaving only one equation in terms of z∗ as(20)z∗2ξ2z∗2λtn+γNCtc+λ+γNC+z∗-λtn-γNCtc+BFtn+λA+BF+γNCzi∗-BFtnξ+ABF=0.This equation can be written in terms of Rτ as(21)P2Rτ2z∗2+-P1Rτ2+P3Rτz∗-P4Rτ+P4=0.The nonzero solutions of (21) can be written as(22)zi∗=P1Rτ-P32P2Rτ1±1-4P2P41-RτP1Rτ-P32.It can be seen from (22) that, for Rτ>1, zi∗ can only have one positive value. We expressed the other variables in terms of zi∗ as(23)xi∗=zi∗2ξγNC1-zi∗λξzi∗tn+ξzi∗γNCtc+BFtn,vi∗=λξzi∗+BFzi∗1-zi∗λξzi∗tn+ξzi∗γNCtc+BFtn,wi∗=Bzi∗γNC1-zi∗λξzi∗tn+ξzi∗γNCtc+BFtn.It is clear from (23) that xi∗, vi∗, wi∗ are positive whenever zi∗ is. Hence, the endemic equilibrium point exists when Rτ>1 and is uniquely given by Ei=(xi∗,vi∗,wi∗,zi∗), i=1or2.
Now supposeP1Rτ-P3>0 and R1<1. It is clear from (22) that both values of zi∗ are positive. In this case, the equilibria can be written as Ee=(xi∗,vi∗,wi∗,zi∗), i=1,2.
This establishes the second part of Theorem3. It can easily be seen from (22) that the other possibilities R1>1, or R1=1, P1Rτ<P3, or R1<1, P1Rτ<P3 will not result in any real positive value of zi∗.Backward bifurcation, which is a situation where the locally asymptotically stable DFE coexists with a locally asymptotically stable endemic equilibrium point, has been observed in many epidemic models. See, for instance, [2, 7, 28, 29]. In this scenario, the requirement for the basic reproduction number to be less than one for the disease to be eradicated no longer holds. Condition (2) of Theorem 3 provides possibility of backward bifurcation in our model. Here we provide parameter range within which backward bifurcation is likely to happen. To do this, it is enough to show that there is a positive endemic equilibrium when the basic reproduction number Rτ<1. We state the result in the following lemma.Lemma 4.
The reduced model (13) undergoes backward bifurcation when P3/P1<Rτb<Rτ≤1, where(24)Rτb=P1P3-2P2P4+2P12P2P41-P3/P1+P22P42P12.To prove this, we first solve for the value of Rτ such that P1Rτ-P32+4P2P4Rτ-1=0, where we obtained two possibilities as(25)Rτb=P1P3-2P2P4+2P12P2P41-P3/P1+P22P42P12,Rτ1=P1P3-2P2P4-2P12P2P41-P3/P1+P22P42P12.Using condition (2) of Theorem 3 and (25), backward bifurcation will occur if Rτ<1 and Rτ>max{P3/P1,Rτb,Rτ1}. Now it is clear that Rτb>Rτ1. Now consider the difference(26)Rτb-P3P1=2-P2P4+P12P2P41-P3/P1+P22P42P12>0.In analogous fashion we show that P3/P1>Rτ1. Hence we get the result.


### 3.5. Stability of the Unique Endemic Equilibrium

Under the hypothesis of condition(1) of the theorem, the endemic equilibrium point of the reduced model (13) exists and is unique. Some of its local stability properties can now be investigated. Without loss of generality, we assume that the endemic equilibrium point is E1=(x1∗,v1∗,w1∗,z1∗) and we dropped the asterisks. The local stability of E1 is now investigated by linearizing the right hand side of model (13) about the equilibrium solution. This gives the Jacobian matrix(27)JE1=-ξz1-FγNCm-ξw1-ξz1-ξz1-A-ξz1-ξv1-ξw1-ξx1+ξξz10-Bξw10-tnz1+tn-tcz1+tc-tcx1-tnv1-1.The characteristic polynomial is(28)η4+A3η3+A2η2+A1η+A0=0,where(29)A3=v1tn+x1tc+2ξz1+A+B+F+1,A2=z1-11-v1-x1-w1tn+z1+2v1z1tnξ+A+B+Fv1tn+2x1ξz1+A+B+Fx1tc+z1-1w1tc+z1ξ+ξ2z12+A+B+F+λ+γNCξz1+AB+AF+BF+A+B+F,A1=-z1+1ξ2z1+ξK1w1+x1tn+Fξz1-1w1tc+tn+z1+ξBz1+z1-1tn+z1ξ2z1-1w1tc+tn+z1+Aξz1-1w1tc+z1+ξ2z12+B+K2+λ+γNCξz1+K1A+BFx1tc+ξ2z1+A+λ+γNCz1+K1ξ+K1A+BFv1tn+λ+γNCξ2z12+Aλ+BF+BγNC+λ+γNCξz1+BF+K1A+BF,A0=Rτz1P1Rτ+P4ART-11-z1tn2,K1=B+F,K2=A+F,RT=z1RτQ1+P42Rτz1P1Rτ+P42,Q1=z1P3Rτ+P4P1+z1-2z1+3P2P4Rτ+-z1+2z12P1P2Rτ2+P3P4. Notice that A3, A2, A1 are all positives and A0 can be positive or negative depending on the sign of RT-1. Hence, by Hurwitz criterion, for E1 to be locally stable, it is necessary that RT>1.


### 3.6. Effects of ITNs Use and Temperature on the Basic Reproduction Number

We evaluate the effects of ITNs use and temperature on the dynamics of malaria on the reduced model (13) through their elasticity indices Rτβ=(β/Rτ)(∂Rτ/∂T) and RτT=(T/Rτ)(∂Rτ/∂β), respectively; see [24, 26] for more explanation on this procedure. The analytical formulations of the elasticity indices of the proportion of ITNs use and temperature with respect to Rτ are given, respectively, by (30) and (31). (30)Rτβ=2βαβmax-βminαβmax-βmin-βmax<0,(31)RτT=2ThT+T0e-ζβmax-βminT-T03e-ζβmax-βmin-βmax,T≠T0,where ζ=hT/T-T02. The implication of (30) is that an increase in the number of ITNs used will bring about decrease in the reproduction number and vice versa. This finding corroborates the report of [24]. From (31), the reproduction number will increase with temperature whenever T<T0 and will decrease with temperature whenever T>T0. In other words, the reproduction number is not a monotonic function of temperature. It should be understood that this type of sensitivity analysis only focuses on specific parameters. It does not indicate the effect of concurrent, large perturbations in all model parameters which is almost always the case in epidemiology. In order to gain more insight into the sensitivities of the conglomerations of the parameters, we carry out more investigations in the next section.


## 3.1. Basic Properties of the Model

Lemma 1.
LetNh∗,Nm∗ be the equilibrium solutions of the total human and mosquito populations, respectively. The closed set D={Sn,In,Sc,Ic,Sv,Iv∈R+6:Nh≤Nh∗,Nm≤Nm∗} is positively invariant and attracting.Proof.
Adding the first four equations and the last two equations of model (8) we obtained(9)dNhdt=λhNh-μh2Nh+μhNh-δNIn-δCIc,dNmdt=λvNm-μv2Nm+μv+βμv1Nm.The mosquito population is modeled by logistic growth with carrying capacity kv=(λv-μv-μv1)/μv2. It is easy to see that dNh/dt≤λh-μhNh1-Nh/kh and dNm/dt≤λv-μvNm1-Nm/kv. This shows that dNh/dt≤0 if 1≥Nh/kh and it approaches kh. Similarly, dNm/dt≤0 if 1≥Nm/kv and it approaches kv. Hence, using comparison theory [2] (10)Nht≤khNh0Nh01-e-λh-μht+khe-λh-μht,Nmt≤kvNm0Nm01-e-λv-μvt+kve-λv-μvt.If Nh(0)≤kh, then Nh(t)≤kh and if Nm(0)≤kv, then Nm(t)≤kv. Thus, the region D is positively invariant for the model. Moreover, if Nh(0)≥kh, Nm(0)≥kv, then either the solution enters the region D in finite time or Nh(t)→kh, Nm(t)→kv, as t→∞. Thus the region attracts all solutions in R+6. Now that we have shown that D is positively invariant, the requirement for existence and uniqueness of solutions holds for the system [2]. The dynamics of the system in the region D will henceforth be investigated.


## 3.2. Scaling

To analyze the malaria model (8), we think it is easier to work with fractional population instead of actual populations by scaling the population of each class by the total species population. We let Sn=uNh, In=vNh, Sc=wNh, Ic=xNh, Iv=zNm, Sv=yNm. We arbitrarily scale the time variables by λv by introducing τ=tλv so that(11)dudτ=1λvNhdSn/dt-SndNh/dtNh2,dvdτ=1λvNhdIn/dt-IndNh/dtNh2and so on for the rest of the variables. In the absence of the disease, the human population follows logistic growth with carrying capacity kh=(λh-μh)/μh2. Following [22] we scaled the human and vector populations in the first 6 equations of model (8) using their respective carrying capacities as Nh=khNh∗, Nm=kvNm∗. So we have a new 6-dimensional system of equations with two additional dimensions for the two total population variables Nh and Nm as(12)dudτ=-a1bkvNm∗khNh∗z-λu+αw+vγnn+λ,dvdτ=-Tj+γnc1+γnn+λ+δnv+a1bkvNm∗khNh∗uz,dwdτ=-wa1bkvNm∗khNh∗z+α+λ+vTj-ϵn+γnc1+xTi+γc,dxdτ=a1bkvNm∗khNh∗wz-xλ+Ti+γc+δc,dydτ=-vytn-xytc-y+1,dzdτ=vytn+xytc-z,dNhdτ=1-Nh1-ω1Nh-Nhvδn+xδc,dNmdτ=1-Nm1-ω2Nm-Nmω3,where λ=λh/λv, a1=a/λv, b1=b/λv, α=αc/λv, γnn=γNN/λv, Ti=Tc/λv, Tj=Tn/λv, ϵn=ϵnn/λv, δc=δC/λv, δn=δN/λv, γnc1=γNC/λv, γc=γC/λv, tn=τn/λv, tc=τc/λv, ω1=μh/λv, ω2=μv/λv, and ω3=μv3/λv.The analytic solution of model (12) will be more tractable if the number of equations is reduced. To this end, we let u=1-v-w-x, y=1-z, and by substituting these in (12), two equations are eliminated leaving us with the reduced model:(13)dvdτ=1-v-w-xξz-Av,dwdτ=-ξz+Fw+mx+vγnc,dxdτ=-Bx+ξwz,dzdτ=-vtn+xtc+1z+vtn+xtc,where ξ=a1bkvρ/kh, ρ=Nm∗/Nh∗, A=Tj+γnc1+γnn+λ+δn, m=Ti+γc, B=Ti+γc+λ+δc, F=α+λ, γnc=Tj-ϵn+γnc1.


## 3.3. Reproduction Number and Disease-Free Equilibrium

The disease-free equilibrium point (DFE) is the solution of the reduced model (13) when the disease classes are identically zero. By setting the right hand sides of (13) and the disease classes to zero, then solving the resulting equations simultaneously, we obtained the DFE of the reduced model (13) as Edfe=w,v,x,z=0,0,0,0. The nonnegative matrix for the new infection (ϝ) and the matrix for the transition terms (V) are given by(14)ϝ=00ξ000tntc0,V=A000B0001,ϝV-1=00ξ000tnAtcB0. The spectral radius is(15)ρϝV-1=ξtnA.We define the reproduction number as(16)Rτ=ξtnA.Lemma 2.
The equilibrium pointEdfe of the reduced model (13) is locally asymptotically stable (LAS) if Rτ<1 and unstable if Rτ>1.Proof.
The Jacobian matrix evaluated atEdfe is(17)JEdfe=-FγNCm00-A0ξ00-B00tntc-1.The characteristic equation is(18)F+ηB+ηηA+1+η2-tnξ+A=0.The eigenvalues are (19)η1=-F,η2=-B,η3=-A+1-A+12-4A1-Rτ2,η4=-A+1+A+12-4A1-Rτ2.Clearly, ηi, i=1,2,3, are negatives and that η4 is negative provided Rτ<1.


## 3.4. Endemic Equilibrium Point and Backward Bifurcation

Further analysis of the model is done by finding the endemic equilibrium solutionsx∗,v∗,w∗,z∗ of the reduced model (13). First we state the following.Theorem 3.
The reduced model has the following:(1)
A single endemic equilibrium solution ifRτ≥1, or R1=1, and P1Rτ>P3.(2)
Two endemic equilibrium solutions ifP1Rτ>P3 and R1=4P2P41-Rτ/P1Rτ-P32<1.(3)
No endemic equilibrium solution otherwise.Here P1=Aλtn+γNCtc, P2=Aλtn+γNCtc+λ+γNC, P4=BFtn2, and P3=tnBFtn+λA+BF+BγNC.Proof.
To prove this, we first expressed the right hand side of the reduced model (13) in terms of the equilibrium solutions x∗, v∗, w∗, z∗. We then eliminated the rest of the variables leaving only one equation in terms of z∗ as(20)z∗2ξ2z∗2λtn+γNCtc+λ+γNC+z∗-λtn-γNCtc+BFtn+λA+BF+γNCzi∗-BFtnξ+ABF=0.This equation can be written in terms of Rτ as(21)P2Rτ2z∗2+-P1Rτ2+P3Rτz∗-P4Rτ+P4=0.The nonzero solutions of (21) can be written as(22)zi∗=P1Rτ-P32P2Rτ1±1-4P2P41-RτP1Rτ-P32.It can be seen from (22) that, for Rτ>1, zi∗ can only have one positive value. We expressed the other variables in terms of zi∗ as(23)xi∗=zi∗2ξγNC1-zi∗λξzi∗tn+ξzi∗γNCtc+BFtn,vi∗=λξzi∗+BFzi∗1-zi∗λξzi∗tn+ξzi∗γNCtc+BFtn,wi∗=Bzi∗γNC1-zi∗λξzi∗tn+ξzi∗γNCtc+BFtn.It is clear from (23) that xi∗, vi∗, wi∗ are positive whenever zi∗ is. Hence, the endemic equilibrium point exists when Rτ>1 and is uniquely given by Ei=(xi∗,vi∗,wi∗,zi∗), i=1or2.
Now supposeP1Rτ-P3>0 and R1<1. It is clear from (22) that both values of zi∗ are positive. In this case, the equilibria can be written as Ee=(xi∗,vi∗,wi∗,zi∗), i=1,2.
This establishes the second part of Theorem3. It can easily be seen from (22) that the other possibilities R1>1, or R1=1, P1Rτ<P3, or R1<1, P1Rτ<P3 will not result in any real positive value of zi∗.Backward bifurcation, which is a situation where the locally asymptotically stable DFE coexists with a locally asymptotically stable endemic equilibrium point, has been observed in many epidemic models. See, for instance, [2, 7, 28, 29]. In this scenario, the requirement for the basic reproduction number to be less than one for the disease to be eradicated no longer holds. Condition (2) of Theorem 3 provides possibility of backward bifurcation in our model. Here we provide parameter range within which backward bifurcation is likely to happen. To do this, it is enough to show that there is a positive endemic equilibrium when the basic reproduction number Rτ<1. We state the result in the following lemma.Lemma 4.
The reduced model (13) undergoes backward bifurcation when P3/P1<Rτb<Rτ≤1, where(24)Rτb=P1P3-2P2P4+2P12P2P41-P3/P1+P22P42P12.To prove this, we first solve for the value of Rτ such that P1Rτ-P32+4P2P4Rτ-1=0, where we obtained two possibilities as(25)Rτb=P1P3-2P2P4+2P12P2P41-P3/P1+P22P42P12,Rτ1=P1P3-2P2P4-2P12P2P41-P3/P1+P22P42P12.Using condition (2) of Theorem 3 and (25), backward bifurcation will occur if Rτ<1 and Rτ>max{P3/P1,Rτb,Rτ1}. Now it is clear that Rτb>Rτ1. Now consider the difference(26)Rτb-P3P1=2-P2P4+P12P2P41-P3/P1+P22P42P12>0.In analogous fashion we show that P3/P1>Rτ1. Hence we get the result.


## 3.5. Stability of the Unique Endemic Equilibrium

Under the hypothesis of condition(1) of the theorem, the endemic equilibrium point of the reduced model (13) exists and is unique. Some of its local stability properties can now be investigated. Without loss of generality, we assume that the endemic equilibrium point is E1=(x1∗,v1∗,w1∗,z1∗) and we dropped the asterisks. The local stability of E1 is now investigated by linearizing the right hand side of model (13) about the equilibrium solution. This gives the Jacobian matrix(27)JE1=-ξz1-FγNCm-ξw1-ξz1-ξz1-A-ξz1-ξv1-ξw1-ξx1+ξξz10-Bξw10-tnz1+tn-tcz1+tc-tcx1-tnv1-1.The characteristic polynomial is(28)η4+A3η3+A2η2+A1η+A0=0,where(29)A3=v1tn+x1tc+2ξz1+A+B+F+1,A2=z1-11-v1-x1-w1tn+z1+2v1z1tnξ+A+B+Fv1tn+2x1ξz1+A+B+Fx1tc+z1-1w1tc+z1ξ+ξ2z12+A+B+F+λ+γNCξz1+AB+AF+BF+A+B+F,A1=-z1+1ξ2z1+ξK1w1+x1tn+Fξz1-1w1tc+tn+z1+ξBz1+z1-1tn+z1ξ2z1-1w1tc+tn+z1+Aξz1-1w1tc+z1+ξ2z12+B+K2+λ+γNCξz1+K1A+BFx1tc+ξ2z1+A+λ+γNCz1+K1ξ+K1A+BFv1tn+λ+γNCξ2z12+Aλ+BF+BγNC+λ+γNCξz1+BF+K1A+BF,A0=Rτz1P1Rτ+P4ART-11-z1tn2,K1=B+F,K2=A+F,RT=z1RτQ1+P42Rτz1P1Rτ+P42,Q1=z1P3Rτ+P4P1+z1-2z1+3P2P4Rτ+-z1+2z12P1P2Rτ2+P3P4. Notice that A3, A2, A1 are all positives and A0 can be positive or negative depending on the sign of RT-1. Hence, by Hurwitz criterion, for E1 to be locally stable, it is necessary that RT>1.


## 3.6. Effects of ITNs Use and Temperature on the Basic Reproduction Number

We evaluate the effects of ITNs use and temperature on the dynamics of malaria on the reduced model (13) through their elasticity indices Rτβ=(β/Rτ)(∂Rτ/∂T) and RτT=(T/Rτ)(∂Rτ/∂β), respectively; see [24, 26] for more explanation on this procedure. The analytical formulations of the elasticity indices of the proportion of ITNs use and temperature with respect to Rτ are given, respectively, by (30) and (31). (30)Rτβ=2βαβmax-βminαβmax-βmin-βmax<0,(31)RτT=2ThT+T0e-ζβmax-βminT-T03e-ζβmax-βmin-βmax,T≠T0,where ζ=hT/T-T02. The implication of (30) is that an increase in the number of ITNs used will bring about decrease in the reproduction number and vice versa. This finding corroborates the report of [24]. From (31), the reproduction number will increase with temperature whenever T<T0 and will decrease with temperature whenever T>T0. In other words, the reproduction number is not a monotonic function of temperature. It should be understood that this type of sensitivity analysis only focuses on specific parameters. It does not indicate the effect of concurrent, large perturbations in all model parameters which is almost always the case in epidemiology. In order to gain more insight into the sensitivities of the conglomerations of the parameters, we carry out more investigations in the next section.


## 4. Global Uncertainty and Sensitivity Analysis

One of the important components of epidemic modeling is parameter estimation. This is because many factors combine to form variability in inputs into the model output. These factors may include erroneous parameter estimation and uncertainty in the exact parameter values. For this reason, it is important to determine parameters that have substantial influence on the results. This can be achieved through sensitivity and uncertainty measurement that can be done more easily due to rapid growth in computing technology. For example, Sampling and Sensitivity Analysis Tools (SaSAT) is software developed for sensitivity analysis [30]. In our model, the proportion of infectious human population, v, x, are regulated by various malaria-related epidemiological parameters that are shrouded with uncertainty. Following [30], we use the Latin Hypercube Sampling (LHS) and Partial Rank Correlation Coefficient (PRCC) techniques to perform a global uncertainty and sensitivity analysis of the reduced model (13). We sample all the twenty-five parameters of the model and then carry out simulations to measure their statistical influence on the proportion of infectious humans v, x. Using published results in literature, we assigned upper and lower bounds for each parameter as shown on Table 3. As in [24, 30], we assumed that each parameter follows a uniform distribution and then draw 1000 samples from the distribution. This gives a 1000×25 matrix whose rows consist of unique collections of parameters. For each row of the matrix, the reduced model (13) is integrated using MATLAB ode45 and we record the proportion of infected humans at each time step. Typical results of such calculations are depicted on Figure 3. The legend and the axes labels on A4 are the same for all the subfigures. From Panel A1 of the figure Rτ>1 and there is no bifurcation. This implies that the disease will invade and will persist. In this case, it is possible to eradicate the disease by using any intervention that will decrease the reproduction number to below unity. Panels A2 and A4 depict the scenario Rτ<1 and there is no bifurcation. This means that the disease will not persist and will not invade. In Panel A3 we have typical scenario where the endemic equilibrium point and the DFE coexist. In this case, the disease will persist but will not invade. The results suggest that it is the presence of the immune individuals that makes the disease persist. The strategy that will be adopted here will be examined through sensitivity analysis. From the 1000 collections of unique parameters, 469 satisfy the condition for backward bifurcation. We selected one of these collections as our baseline values as shown on Table 3 except for Ti. We conducted sensitivity analysis on the bifurcation parameter using PRCC with the baseline values. To achieve this, we set Ti=0.5532 and allow c1 to vary in the range [0.2299,0.35]. The results are depicted on Figure 4. The existence of backward bifurcation implies that the classical requirement for Rτ<1 for an epidemic to be eradicated no longer holds. The region to the left of the bifurcation point represents the area where malaria will not invade and will not persist (comfort zone), whereas the region enclosed by the solid red and the dash-dotted curves and Rτ=1 represent the area where malaria will persist but will not invade. If malaria is to be eradicated, the comfort zone should be maximized subject to Rτ≤1 or we must have Rτ<Rτb. The natural question is what are those parameters that are the most important to the bifurcation parameter? To investigate this, we conducted sensitivity analysis on the bifurcation parameter using the 469 samples mentioned earlier. The results of this analysis are depicted on Figure 5. From the figure, the parameters c1, Ti,  c2, Tj are the most sensitive to the bifurcation parameter in that order, while in decreasing order of importance the following have the least significance in terms of bifurcation; βmin, μh2, b, h, and ϵnn. Increasing c1 and Ti will bring about increase in the bifurcation parameter. The public health implication of this result is that increase in treatment of the clinically immune individuals will increase the comfort zone. However, parameter sensitivities on bifurcation parameter should not be taken in isolation. Thus, we draw new 1000 samples for all the 19 parameters affecting Rτ and using these samples, we carry out sensitivity analysis on Rτ using LHS and PRCC. The results of these analyses are depicted on Figure 6. We list the following parameters in order of their importance on Rτ, βmax, μh2, μ2, and T from Figure 6, while the following have least importance and are listed in decreasing order of importance: μv, λv, δN, and μv3.Table 3
Parameters, their baseline values, and the ranges for sensitivity analysis.
Parameter
Baseline value
Range
Reference
λh
9.3614×10-5
[
2.7,14
]
×
1
0
-
5
[26]
λv
0.4478
[
0.27,0.7
]
[8, 26]
Tn
0.1201
[
0.03704,0.4
]
[21]
Tc
0.0476
[
0.003704,0.2
]
[21]
c1
0.2299
[
0.1,0.7
]
[24]
c2
0.5342
[
0.072,0.64
]
[24]
ρ
7
[
2,8
]
[27]
b
0.8
[
0.1,0.8
]
[8, 26]
T
22
[
20,35
]
Estimated
T0
30.6
[
27,31
]
Estimated
βmax
0.6334
[
0.1,1
]
[24]
βmin
0.0696
[
0,0.1
]
[24]
h
0.3165
[
0.1,0.5
]
Estimated
ϵnn
0.2201
[
0.1,0.5
]
Estimated
γNN
0.0805
[
0.056,0.2
]
[24]
γNC
0.0106
[
0.0014
,
0.017
]
[7, 24]
αc
0.0022
55
×
10
-
6
,
11
×
10
(
-
3
)
[7, 24]
δN
5.7341
×
1
0
-
5
[
0,4.1
]
×
1
0
-
4
[26]
δC
3.2084
×
1
0
-
4
[
0,4.1
]
×
1
0
-
4
[26]
γC
0.0101
[
0.0014,0.017
]
[7, 24]
μh2
6.0146
×
1
0
-
7
[
1,100
]
×
1
0
-
8
[26]
μh
1.6728
×
1
0
-
5
[
1,20
]
×
1
0
-
6
[26]
μv
0.0668
[
1,100
]
×
1
0
-
3
[26]
μv2
6.8754
×
1
0
-
4
[
1,1000
]
×
1
0
-
6
[26]
μv3
0.0995
[
0.0476,0.1
]
[24]Figure 3
Time course of the reduced model with initial condition(w,v,x,z)=(0.1,0.2,0.005,0.3). Panel A1 corresponds to case Rτ=4.0573, Rbτ=0.3920. Panel A2 corresponds to case Rτ=0.5341, Rbτ=0.9952. Panel A3 corresponds to case Rτ=0.3397, Rbτ=0.0901. Panel A4 corresponds to case Rτ=0.092, Rbτ=0.4143. The legend and labels on A4 are the same for all the subfigures.Figure 4
(a) Proportion of infected naive population. (b) Proportion of clinically infected population.
(a)
Bifurcation diagram for the proportion of naive population
(b)
Bifurcation diagram for the proportion of clinically immune populationFigure 5
Tornado plot of sensitivity analysis of the bifurcation parameterRτb.Figure 6
PRCC results for all the parameters affecting the threshold parameterRτ.


### 4.1. Intervention Strategies

In our model, there are many biological parameters that influence disease dynamics. Also, there are intervention parameters that are crucial for eliminating an epidemic. These areβ, Ti, Tj. We consider a range of different intervention strategies by taking 1000 samples using LHS for 19 parameters of the model. The intervention parameters were not sampled but are given specific values as shown on Table 4. Similarly, T, T0, and h were not sampled because β is constant for a given run of simulation. We used MATLAB ode45 suite to integrate the reduced system (13) using the initial condition: (w,v,x,z)=(0.1,0.2,0.005,0.3). For each run of simulation, we recorded the sum of proportions of infected human over the entire time steps. The results of these intervention strategies are presented in the form: NNN≡Tj=0, Ti=0, β=0; NNY≡Tj=0, Ti=0, β=0.9; NYN≡Tj=0, Ti=0.2, β=0, and so forth. The first means no treatment of both the naive and clinically infected populations and no ITN coverage. The second means coverage, no treatment of both the naive and clinically infected populations with 90% ITN coverage. Altogether there are 8 different interventions which correspond to 8000 simulations. The box plot of the results is depicted on Figure 7. The strategy YYY appears to be the best. It skews to the right and is characterized by low median and small interquartile range and the presence of large upper outliers. The strategies NNY and NYY are the better options after YYY. These further demonstrate the effectiveness of ITNs usage and that any intervention with ∗∗N (any intervention that does not constitute ITNs coverage) is relatively ineffective. The public health interpretation of this result is that reducing contacts between humans and mosquitoes is important in controlling the size of the infectious human population. We conducted further sensitivity analysis by integrating the time dependent dynamics of (13) using 1000 parameter sets used in plotting Figure 3 and the same initial conditions. To do that, the nondimensional time [0,30] was divided into 50 times steps and we recorded the number of infected humans at each time step. The results of these analyses are shown on Figure 8 for the nondimensional time τ=20. The PRCC for γNN is the largest and positive, which indicates that rate of recovery without immunity can increase the number of infected humans. The PRCC for δC is the second largest in magnitude and has negative sign. This indicates that death rate of clinically immune individuals will decrease the number of infectious humans. The parameter βmin has the least importance. The PRCCs results depicted on Figure 8 correspond to a single time step; however, the dynamics of the PRCCs is likely to vary over time. Hence, to fully characterize how sensitive the infectious human population is to the parameters of the reduced system (13), we look into the evolution of the PRCCs over time and the results are depicted on Figures 9 and 10 for 0≤τ≤30. We grouped the parameters in the figures for the purpose of clarity of presentation only. The results on these figures indicate that, initially, μh2 and μv2 are the most significant parameters. While increase in μv2 will result in decreasing the number of infected humans, increase in μh2 will lead to increase in the number of infected human. See panel 1 on Figure 9. The maximum biting rate βmax and the density dependent part of mosquitoes death rate appear to be the most important parameters in the long run. The public health implication of this is that increased ITNs coverage and death rate of mosquitoes can lead to the control of malaria transmission.Table 4
Parameter values used for interventions.
Parameter/symbolic representation
N
Y
Tj
0
0.5
Ti
0
0.2
β
0
0.9Figure 7
Box plots showing the sum of proportion of infected human populations subjected to different intervention strategies.Figure 8
PRCC results for all the parameters of the reduced model, simulated at timeτ=20.Figure 9
Time dependent PRCC results for 16 parameters of the model. We can infer from the figure that the maximum biting rate and the mosquitoes death rate have the most significant impact in the long run.Figure 10
Time dependent PRCC results for 9 parameters of the model. We can infer from the figure that treatment rate for naive infected individuals can have both positive and negative values of PRCCs.


## 4.1. Intervention Strategies

In our model, there are many biological parameters that influence disease dynamics. Also, there are intervention parameters that are crucial for eliminating an epidemic. These areβ, Ti, Tj. We consider a range of different intervention strategies by taking 1000 samples using LHS for 19 parameters of the model. The intervention parameters were not sampled but are given specific values as shown on Table 4. Similarly, T, T0, and h were not sampled because β is constant for a given run of simulation. We used MATLAB ode45 suite to integrate the reduced system (13) using the initial condition: (w,v,x,z)=(0.1,0.2,0.005,0.3). For each run of simulation, we recorded the sum of proportions of infected human over the entire time steps. The results of these intervention strategies are presented in the form: NNN≡Tj=0, Ti=0, β=0; NNY≡Tj=0, Ti=0, β=0.9; NYN≡Tj=0, Ti=0.2, β=0, and so forth. The first means no treatment of both the naive and clinically infected populations and no ITN coverage. The second means coverage, no treatment of both the naive and clinically infected populations with 90% ITN coverage. Altogether there are 8 different interventions which correspond to 8000 simulations. The box plot of the results is depicted on Figure 7. The strategy YYY appears to be the best. It skews to the right and is characterized by low median and small interquartile range and the presence of large upper outliers. The strategies NNY and NYY are the better options after YYY. These further demonstrate the effectiveness of ITNs usage and that any intervention with ∗∗N (any intervention that does not constitute ITNs coverage) is relatively ineffective. The public health interpretation of this result is that reducing contacts between humans and mosquitoes is important in controlling the size of the infectious human population. We conducted further sensitivity analysis by integrating the time dependent dynamics of (13) using 1000 parameter sets used in plotting Figure 3 and the same initial conditions. To do that, the nondimensional time [0,30] was divided into 50 times steps and we recorded the number of infected humans at each time step. The results of these analyses are shown on Figure 8 for the nondimensional time τ=20. The PRCC for γNN is the largest and positive, which indicates that rate of recovery without immunity can increase the number of infected humans. The PRCC for δC is the second largest in magnitude and has negative sign. This indicates that death rate of clinically immune individuals will decrease the number of infectious humans. The parameter βmin has the least importance. The PRCCs results depicted on Figure 8 correspond to a single time step; however, the dynamics of the PRCCs is likely to vary over time. Hence, to fully characterize how sensitive the infectious human population is to the parameters of the reduced system (13), we look into the evolution of the PRCCs over time and the results are depicted on Figures 9 and 10 for 0≤τ≤30. We grouped the parameters in the figures for the purpose of clarity of presentation only. The results on these figures indicate that, initially, μh2 and μv2 are the most significant parameters. While increase in μv2 will result in decreasing the number of infected humans, increase in μh2 will lead to increase in the number of infected human. See panel 1 on Figure 9. The maximum biting rate βmax and the density dependent part of mosquitoes death rate appear to be the most important parameters in the long run. The public health implication of this is that increased ITNs coverage and death rate of mosquitoes can lead to the control of malaria transmission.Table 4
Parameter values used for interventions.
Parameter/symbolic representation
N
Y
Tj
0
0.5
Ti
0
0.2
β
0
0.9Figure 7
Box plots showing the sum of proportion of infected human populations subjected to different intervention strategies.Figure 8
PRCC results for all the parameters of the reduced model, simulated at timeτ=20.Figure 9
Time dependent PRCC results for 16 parameters of the model. We can infer from the figure that the maximum biting rate and the mosquitoes death rate have the most significant impact in the long run.Figure 10
Time dependent PRCC results for 9 parameters of the model. We can infer from the figure that treatment rate for naive infected individuals can have both positive and negative values of PRCCs.


## 5. Discussion

Mosquitoes are very efficient vectors of human diseases and are responsible for transmitting some of the many devastating diseases today. For many of these diseases, climate change is key in determining the ability of a mosquito to transmit the disease through the population effectively. One of such diseases is malaria, which is widely considered as the most devastating and the most prevalent human vector borne disease, with one-half of the world population living in areas where there is risk of infection [14]. Two of the factors that complicate malaria control are as follows.(1)
Seasonal changes: This has strong influence on malaria transmission which makes it difficult to predict future malaria intensity accurately. There is no single definition of seasonality in relation to malaria in the literature. Malaria metrics such as the mosquito biting rates have been investigated in many studies in relation to temperature; see for instance [15]. Numerous values of mosquitoes biting rates have been reported in literature. In many studies, the biting rate has constant values and these values are mostly different. Given the importance of biting rate as a driving force for malaria transmission, there is the need to conduct further studies on this so as to discern the factors that bring about changes in the biting rate. Typical of these factors is seasonality. This is because transmission rate of malaria has been reported to peak in certain period of the year. It is quite reasonable to model the biting rates as a function of temperature in order to mimic seasonality.(2)
Clinical immunity: Another factor that complicates malaria control is that individuals living in regions where malaria is endemic can develop immunity to malaria which enables them to remain asymptomatic while still carrying the parasites [7, 31]. The development of acquired clinical immunity by individuals will result in such individuals not seeking treatment for a long time. Thus, they will harbor the disease and can transmit it when bitten by mosquitoes.From Panels A1 and A3 of Figure 3, the infectivity of the clinically infected individuals does not wane with time. This might be partly attributed to the reason outlined in bullet (2) above. In view of this finding, one may pause the question; how do we control malaria without treatment? We believe that treatment is inevitable once malaria cases are established. In this case one can consider a number of other intervention strategies that can be combined together with treatment towards control and eradication of malaria transmission. Modeling the dynamics of malaria by separating human population into immune and nonimmune classes, incorporating the temperature and ITNs use as part of mosquitoes biting rates and sensitivity analysis is the main focus of this paper.Our simple model of temperature dependent biting rate generalizes calibrated model results reviewed by [15]. The local sensitivity analysis indices results are in line with literature reports. For instance, it is common knowledge that the use of fan in our houses could reduce contact between human and mosquitoes and this can lead to reduction in transmission of malaria. The local sensitivity analysis of the reproduction number in relation to temperature results also supports this observation and other similar findings in many reports. See [15] and the references therein.To eradicate malaria in bistable regions, we need to makeRτb>Rτ. This can be achieved by considering strategies that will decrease Rτ or increase Rτb sufficiently. The sensitivity analysis results reveal that we need to embark on combination of strategies that will decrease the possibilities of immune individuals becoming infective and treating them in order to effectively reduce Rτ. This is quite difficult because clinically immune population are not likely to seek treatment for a long time. Moreover, failing to detect the parasite in them may not necessarily mean the absence of the disease. In view of these difficulties, the best thing to do is to embark on maximum ITNs usage. Is there any way in which the rate at which immune individuals lose their clinical immunity can be accelerated? if so, then it will help towards increasing Rτb and hence malaria eradication in a bistable region. Other metrics that will assist significantly in this direction are the increase in mosquitoes death rates and temperature. In some African countries such as Nigeria where malaria is endemic and electricity supply is erratic, improved supply of electricity will enable people to use fans and air conditioning system to provide relatively colder environment in their houses. It will also help to provide enabling environment for effective usage of ITNs and also reduces human-mosquitoes contact.In general, the sensitivities of the model parameters are time dependent. In other words, the most sensitive parameter of the model at the onset of the epidemic may not necessarily maintain the same status at later times. Of particular interest is the fact that the maximum biting and death rates of mosquitoes are the most important parameters in the long run. The public health implication of this is that combined effort should be put in place in reducing human-mosquitoes contact and reduction of mosquito population by any means possible.


## 6. Conclusion

In this paper, we formulated a deterministic model of malaria transmission by dividing the human population broadly according to whether an individual is immune to clinical symptoms of malaria or not. We proposed a simple model with mosquitoes biting rate that is temperature dependent so as to mimic seasonality. We investigated stability conditions for the disease-free equilibrium point and we have found the necessary condition for the unique endemic equilibrium point to be locally stable. We have also shown the range of parameter values in which backward bifurcation can occur. We conducted global sensitivity analysis using Latin Hypercube Sampling and Partial Rank Correlation Coefficient in order to identify the most important parameters that govern the dynamics of malaria transmissions.We find that the combination of treatment and usage of insecticide treated net is the most effective strategy for malaria control. Provision of relatively colder environment will positively impact malaria control. Mosquitoes biting and death rates are the most important parameters of malaria transmissions. The peak of mosquitoes biting rates occurs not at a single temperature value but as a range of values. Temperature is positively correlated to the reproduction number. The disease-free equilibrium point of the model is locally asymptotically stable when the reproduction number is less than one and unstable when it is greater than one.



---

*Source: 6182492-2017-08-01.xml*