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

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.


Introduction
Malaria is one of the most devastating infectious diseases in the world and is caused by Plasmodium 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 malariaendemic 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][4][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][8][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 2 International Scholarly Research Notices 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][14][15][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][19][20][21][22][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 vectorhost 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.

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 ℎ is divided into mutually exclusive subpopulations of susceptible naive , susceptible clinically immune , infected naive , and infected clinically immune , so that ℎ = + + + . The total mosquitoes population denoted by is divided into compartments of susceptible and infected mosquitoes, so that = V + V . All recruitment is assumed to be into the susceptible naive human population generated via birth and/or immigration at a rate ℎ ℎ . The population of naive susceptible individuals ( ) is increased by naive infected individuals that recovered without immunity at a rate , treated naive individuals that recovered without immunity at a rate , and clinically immune individuals that lost immunity at a rate . The population of naive susceptible individuals ( ) is decreased by natural death rate ℎ = ( ℎ2 ℎ + ℎ ) and force of infection (Λ ℎ ), following effective contacts with infected mosquito. Here ℎ2 and ℎ 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 Λ ℎ = ( ( )) V / ℎ .
Here is the probability of infection of susceptible human per bite by an infected mosquito and ≡ ( ( )) is the biting rate of mosquitoes on susceptible human; ≡ ( ) represents the proportion of ITN usage and depends on environmental temperature . See (7) for the functional form of . We assumed that temperature is a parameter that is time independent to make the analysis easier. Thus, 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 , infected clinically immune individuals that recovered with clinical immunity at a rate , and treated clinically immune infected individuals that recover at a rate . It is decreased by susceptible clinically immune individuals that lose immunity at rate , natural death at rate ℎ , 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 Λ ℎ . We assume that clinically infected individuals recover into the clinically susceptible compartment only. Thus, 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 , infected naive individuals that recover without immunity at a rate , infected naive individuals  that recover with immunity at a rate , the natural death ℎ , and disease induced death rate . Thus, 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 , population of clinically immune individuals that recover with immunity at a rate , the natural death ℎ , and disease induced death at a rate . Thus, The population of susceptible mosquitoes is generated by birth at a rate V . It is reduced by natural death V = ( V2 + V ), contact with ITNs at the rate V3 , and infection when in contact with infected human at rate Λ V = ( + )/ ℎ . The parameters V and V2 represent the density independent and density dependent parts of the mosquitoes death rate, respectively. Here = 1 , = 2 , and 1 and 2 are the probabilities that susceptible mosquitoes become infectious after biting an infected naive or clinically immune human, respectively. Thus, 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, 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 where = −ℎ( /( − 0 )) 2 represents the proportion of ITN usage and the parameters 0 and ℎ 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). It 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 Figure 2, and the state variables and parameters of the model are described in Tables 1 and 2, respectively:

Basic Properties of the Model
Lemma 1. Let * ℎ , * be the equilibrium solutions of the total human and mosquito populations, respectively. The closed set Proof. Adding the first four equations and the last two equations of model (8) we obtained International Scholarly Research Notices 5 Human deaths

Mosquitoes deaths
Mosquitoes population Human population The mosquito population is modeled by logistic growth with carrying capacity If ℎ (0) ≤ ℎ , then ℎ ( ) ≤ ℎ and if (0) ≤ V , then ( ) ≤ V . Thus, the region D is positively invariant for the model. Moreover, if ℎ (0) ≥ ℎ , (0) ≥ V , then either the solution enters the region D in finite time or ℎ ( ) → ℎ , ( ) → V , as → ∞. Thus the region attracts all solutions in 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.

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 We arbitrarily scale the time variables by V by introducing = V so that and so on for the rest of the variables. In the absence of the disease, the human population follows logistic growth with carrying capacity ℎ = ( ℎ − ℎ )/ ℎ2 . Following [22] we scaled the human and vector populations in the first 6 equations of model (8) using their respective carrying capacities as ℎ = ℎ * ℎ , = V * . So we have a new 6-dimensional system of equations with two additional dimensions for the two total population variables ℎ and as , and by substituting these in (12), two equations are eliminated leaving us with the reduced model: where = 1 V / ℎ , = * / * ℎ , = ( + 1 + + + ), = ( + ), = ( + + + ), = + , = − + 1 .

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) The spectral radius is We define the reproduction number as Proof. The Jacobian matrix evaluated at dfe is ] .
The characteristic equation is The eigenvalues are Clearly, , = 1, 2, 3, are negatives and that 4 is negative provided < 1.

Endemic Equilibrium Point and Backward Bifurcation.
Further analysis of the model is done by finding the endemic equilibrium solutions * , V * , * , * of the reduced model (13). First we state the following. Proof. To prove this, we first expressed the right hand side of the reduced model (13) in terms of the equilibrium solutions * , V * , * , * . We then eliminated the rest of the variables leaving only one equation in terms of * as * 2 ( 2 ( * 2 ( + + + ) This equation can be written in terms of as International Scholarly Research Notices 7 The nonzero solutions of (21) can be written as It can be seen from (22) that, for > 1, * can only have one positive value. We expressed the other variables in terms of * as * = * 2 (1 − * ) ( * + * + ) , It is clear from (23) that * , V * , * are positive whenever * is. Hence, the endemic equilibrium point exists when > 1 and is uniquely given by = ( * , V * , * , * ), = 1 or 2.
This establishes the second part of Theorem 3. It can easily be seen from (22) that the other possibilities 1 > 1, or 1 = 1, 1 < 3 , or 1 < 1, 1 < 3 will not result in any real positive value of * .
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 < 1. We state the result in the following lemma.

Lemma 4. The reduced model (13) undergoes backward bifurcation when
To prove this, we first solve for the value of such that Using condition (2) of Theorem 3 and (25), backward bifurcation will occur if < 1 and In analogous fashion we show that 3 / 1 > 1 . Hence we get the result.

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 1 = ( * 1 , V * 1 , * 1 , * 1 ) and we dropped the asterisks. The local stability of 1 is now investigated by linearizing the right hand side of model (13) about the equilibrium solution. This gives the Jacobian matrix ] .
The characteristic polynomial is where 3 = V 1 + 1 + 2 1 + + + + 1, Notice that 3 , 2 , 1 are all positives and 0 can be positive or negative depending on the sign of ( − 1). Hence, by Hurwitz criterion, for 1 to be locally stable, it is necessary that > 1.

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 = ( / )( / ) and = ( / )( / ), 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 are given, respectively, by (30) and (31).
where = ℎ /( − 0 ) 2 . 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 < 0 and will decrease with temperature whenever > 0 . 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.

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, , 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, . 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  > 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 < 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 . We conducted sensitivity analysis on the bifurcation parameter using PRCC with the baseline values. To achieve this, we set = 0.5532 and allow 1 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 < 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 = 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 ≤ 1 or we must have < . 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 1 , , 2 , 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 , ℎ2 , , ℎ, and . Increasing 1 and 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 and using these samples, we carry out sensitivity analysis on 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 , max , ℎ2 , 2 , and from Figure 6, while the following have least importance and are listed in decreasing order of importance: V , V , , and V3 .

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 , , . 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, , 0 , and ℎ were not sampled because is constant for a given run of simulation. We used MATLAB ode45 suite to integrate the     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 is the largest and positive, which indicates that rate of recovery without immunity can increase the number of infected humans. The PRCC for 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, ℎ2 and V2 are the most significant parameters. While increase in V2 will result in decreasing the number of infected humans, increase in ℎ2 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.

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 International Scholarly Research Notices Value of correlation coefficient at   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 make > . This can be achieved by considering strategies that will decrease or increase 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 . 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 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.

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.