MODELING THE EFFECT OF TEMPERATURE VARIABILITY ON MALARIA CONTROL STRATEGIES

In this study, a non-autonomous (temperature dependent) and autonomous (temperature independent) models for the transmission dynamics of malaria in a population are designed and rigorously analysed. The models are used to assess the impact of temperature changes on various control strategies. The autonomous model is shown to exhibit the phenomenon of backward bifurcation, where an asymptotically-stable disease-free equilibrium (DFE) co-exists with an asymptotically-stable endemic equilibrium when the associated reproduction number is less than unity. This phenomenon is shown to arise due to the presence of imperfect vaccines and disease-induced mortality rate. Threshold quantities (such as the basic offspring number, vaccination and host type reproduction numbers) and their interpretations for the models are presented. Conditions for local asymptotic stability of the disease-free solutions are computed. Sensitivity analysis using temperature data obtained from Kwazulu Natal Province of South Africa [K. Okuneye and A.B. Gumel. Mathematical Biosciences 287 (2017) 72–92] is used to assess the parameters that have the most influence on malaria transmission. The effect of various control strategies (bed nets, adulticides and vaccination) were assessed via numerical simulations. Mathematics Subject Classification. 37C75, 92D30. Received November 16, 2019. Accepted November 1, 2020.

surface wetness limits the population size of aquatic mosquitoes [37]. In addition, temperature is a key determinant of environmental suitability for transmission of malaria in human. The spatial limits of the distribution and seasonal activities of malaria are sensitive to climate factors, as well as the local capacity to control the disease [8]. The dynamics and distribution of malaria strongly depend on the interplay between the parasite, the mosquitoes and the environment [28,34], it has recently been shown that mosquito and parasite biology are influenced not only by average temperature, but also by the extent of the daily temperature variation [5]. At the extremes, temperature regimes constrain the geographical extent of the disease and, within this envelope, contribute to determining its intensity [22]. These constraints are temporally dynamic, with fluctuations in transmission suitability and intensity driven by seasonal and inter-annual temperature cycles. The importance of temperature as an environmental determinant of malaria endemicity arises from a series of effects on the life cycles of the plasmodium parasite and anopheles mosquitoes [22].
Although, there have been tremendous success in the reduction of malaria cases especially in Africa, nevertheless, mortality due to the disease incidences still remains high in comparison to other infections. Over the years, there have been several initiatives aimed at ending malaria cases, some of which include WHO's Roll Back Malaria program, the Multilateral Initiative in Malaria, the Medicines for Malaria Venture, the Malaria Vaccine Initiative, and the Global Fund to Fight AIDS, TB and Malaria, which supports the implementation of prevention and treatment programs [38]. In fact, there is no single way of preventing malaria, however, there are a number of ways to decrease the transmission of the disease which include the use of treated nets, as well as the use of larvicides and adulticides to clear mosquito breeding sites and kill adult mosquitoes, respectively. Although there is no specific effective vaccine for malaria at the moment, a number of candidate vaccines targeting different stages of the malaria parasite life-cycle have been developed or are currently under development, in particular, RTS,S/AS01 is a strong candidate for the prevention of Plasmodium falciparum infection, in fact phase 3 trials of the vaccine has been completed [41]. Predicting the public health impact of a candidate malaria vaccine requires using clinical trial data to estimate the vaccine's efficacy profile, initial efficacy after vaccination and the pattern through which the vaccine efficacy wanes over time. With an estimated vaccine efficacy profile, the effects of vaccination on malaria transmission can be simulated with the aid of mathematical models [41].
A number of mathematical models have been developed in the literature to gain insights into the effects of temperature, climate change (or seasonality) in the transmission dynamics of mosquito borne diseases in a community, see for instance [1-3, 7, 16, 17, 19, 21, 22, 25-28, 31, 32, 34, 36, 43]. Eikenberry and Gumel [17] extensively review the idea of mathematical modeling of climate change in the transmission dynamics of malaria. Greenhalgh and Moreim [21] study SIRS epidemic model with general seasonal variation in the contact rate. Moreim and Greenhalgh [26] applies generalised periodic vaccination strategy to model for control of the spread and transmission of an infectious disease with seasonal varying contact rate. Dushoff [16] incorporate seasonality of birth rate into the standard deterministic SIR and SEIR epidemic models and identify parameter regions in which birth seasonality can be expected to have observable epidemiological effects. Bury et al. [7] develop a simple socio-climate model by coupling an Earth system model to a social dynamics model, they concluded that socio-climate models should be included in the ensemble of models used to project climate change. In particular, malaria has received lots of attention. Mordecai et al. [27] considered a non-linear response of mosquito and malaria parasite to temperature which are closely consistent with field data, the work which changed predictions on how temperature change affects malaria predicts optimal malaria transmission at 25 • C. A malaria transmission model with periodic birth rate and age structure for the vector population was rigorously analysed by Loy and Zhao [25]. The examination of the process via which parasite development within the mosquito (extrinsic incubation period) is expected to vary over time and space, depending on the diurnal temperature range and baseline mean temperature in Kenya and across Africa was presented by Blanford et al. [5]. Agusto et al. [2] consider a temperature-dependent deterministic model that gave some qualitative insights into the effects of temperature variability on malaria transmission dynamics, the model incorporates gradual increase in infection-acquired immunity via repeated exposure to malaria infection. The impact of variability in temperature and rainfall on the transmission dynamics of malaria in age-structured population, with the dynamics of immature and mature mosquitoes was also considered by Okuneye and Gumel [31]. A malaria model that qualitatively studied the effect of seasonal variations (wet and dry seasons) on the spread of malaria was introduced and analysed by Dembele et al. [13], the authors present an explicit formulation for the basic reproduction ratio and stability analysis of the disease free equilibrium. Olaniyi et al. [33] present malaria model that takes into account time-dependent control functions such as personal protection, treatments and mosquito reduction intervention strategies.
In this study, we consider both autonomous and non-autonomous deterministic models for the transmission dynamics of malaria with control in the presence of temperature variability. Notable feature of the model is that the host population is basically divided according to their vaccination status and the use of multiple control and prevention strategies. The model further assumed that recovered individuals do not become wholly susceptible. The paper is organized as follows. A non-autonomous malaria model is developed and some of its dynamical properties are discussed in Section 2. The autonomous model is analysed for its dynamical properties in Section 3. Endemic equilibrium and the backward bifurcation analysis of the autonomous model is presented in Section 4. Section 5 presents the analysis of the non-autonomous system. The effect of various control strategies is provided in Section 6. Sensitivity analysis and numerical simulations are presented in Section 7.

Model formulation
The total human population at time t, denoted by N H (t), is divided into populations of vaccinated and non-vaccinated individuals. The sub-population of non-vaccinated individuals are further sub-divided into 5 mutually exclusive sub-population of wholly susceptible (without ever been infected) (S U (t)), susceptible after recovery (W U (t)), exposed (E U (t)), infected (I U (t)) and recovered (R U (t)) humans. Similarly, the sub-population of vaccinated individuals are sub-divided into wholly susceptible (S V (t)), susceptible after recovery (W V (t)), exposed (E V (t)), infected (I V (t)) and recovered (R V (t)) humans, so that the total human population at time t is given by In order to assess the potential effect of temperature dependent oviposition of mosquitoes, the population of mosquitoes are divided into aquatic and non-aquatic stages. The aquatic stage (which involves egg, larva and pupa) is represented by a single equation (A M (t)). The non-aquatic (adult) stage is further divided into susceptible (M U (t)), exposed (M E (t)) and infected (M I (t)) mosquitoes. so that the total adult mosquito population (in the non-aquatic stage) is given by Note that, in this study only female mosquitoes are considered as male mosquitoes are non-infectious. The model incorporates the use of larvicides (to clear aquatic mosquitoes) and adulticides (to kill matured mosquitoes). In either case, the death rates due to the use of larvicides and adulticides are proportional to successful rates of applications of larvicides and adulticides.

Dynamics of humans
The population of wholly susceptible humans is generated by birth (or immigration) at a constant rate Π H . This population increases through the loss of vaccination-acquired immunity by wholly vaccinated individuals (at a waning rate ω V ). It is decreased by vaccination (at a rate ξ V which move to the class of wholly vaccinated humans). Proportion of this individuals (in S U class) acquire malaria infection following effective contact with infectious mosquitoes in M I class at a temperature dependent rate λ H (T ), given by where the parameter β V H is a transmission probability from infectious mosquitoes to susceptible humans. The parameter a M (T ) is the temperature dependent biting rate of mosquitoes, 0 < B < 1 represents efficacy rate of bed nets and α B measure compliance rate in the use of bed nets (low values of α B imply limited use of bed nets by the members of the public, whereas values of α B close to one imply widespread use of bed nets). Therefore B α B represents the use of insect repellents to minimize contacts with mosquitoes. Natural mortality occurs in all human classes at a rate µ H , so that The population of a wholly vaccinated individuals S V is generated by vaccination of susceptible individuals at the rate ξ V . This population decreases due to waning of vaccine (at the rate ω V ), infection at the rate λ H (T )(1 − V ) (where 0 < V < 1 is a vaccine efficacy) and by natural death (at the rate µ H ). This gives The populations of non-vaccinated (W U ) and vaccinated (W V ) susceptible individuals (who are partially immune due to prior infection) are generated following loss of partial immunity by recovered individuals that are non-vaccinated and vaccinated at the rates τ U and τ V , respectively. These populations decrease by secondary infection at a reduced rate λ H (T )(1 − W ) (where W is a protection rate due to prior malaria infection) and by natural death. So that The populations of non-vaccinated exposed (E U ) and vaccinated exposed (E V ) individuals are generated by the infection of non-vaccinated (S U , W U ) and vaccinated (S V , W V ) individuals at the rates λ H , λ H (T )(1 − W ) and λ H (1 − V ), λ H (T )(1 − W ), respectively. These populations are reduced by progressing to non-vaccinated and vaccinated infectious classes at the rates σ U and σ V , respectively, and by natural death. Thus The populations of non-vaccinated infectious (I U ) and vaccinated infectious (I V ) individuals are generated by progression of non-vaccinated and vaccinated exposed individuals to the infectious classes at the rates σ U and σ V , respectively. These populations decreases by recovery at the rates γ U and γ V , respectively, natural death and disease-induced death at the rates δ U and δ V , respectively. This gives The populations of non-vaccinated recovered (R U ) and vaccinated recovered (R V ) individuals are generated by the recovery of non-vaccinated and vaccinated infectious individuals at the rates γ U and γ V , respectively. The populations are reduces by loss of partial immunity at the rates τ U and τ V , respectively, and natural death. So that

Dynamics of mosquitoes
In the presence of intervention (using larvicides), the population of aquatic mosquitoes (eggs, larvae and pupae) increases through oviposition by reproductive mosquitoes at a temperature dependent rate φ A (T ) 1 − α L L , where α L is a rate of applying larvicides and L is an efficacy of larvicides (so that α L L = c L accounts for effectiveness of larvicides). The growth of aquatic mosquitoes is moderated by a constant environmental carrying capacity K. This population decreases by maturation at a temperature dependent rate σ A (T ), die naturally and due to the use of larvicides at a rate µ A (T ) 1 + α L L , where µ A (T ) is a temperature dependent death rate. Thus The population of susceptible adult female mosquitoes (M U (t)) is generated by maturation of aquatic mosquitoes at the temperature dependent rate σ A (T ). It decreases by acquiring infection following a substantial contact with an infectious human at a temperature dependent infection rate λ V (T ), given by where, β HV is the probability of infection from infectious humans to susceptible mosquitoes. The parameters η I , η U and η V are modification parameters which account for the reduction in the infectivity of individuals in I V , R U and R V classes in comparison to those in I U class, respectively. Similarly, the population is further decreases at a rate µ V (T ) 1 + α A A , where µ V (T ) is a temperature dependent death in the absence of intervention, α A is a rate of applying adulticides and A is an efficacy of adulticides (so that α A A = c A accounts for effectiveness of indoor residual spraying). Thus The population of exposed mosquitoes in the M E (t) class is generated by the infection of adult mosquitoes in the M U (t) class at the rate λ V (T ). This population decreases by progression to infectious class at a temperature dependent rate σ M (T ) and die naturally and due to the use of adulticides at the rate µ V (T ) 1 + α A A . Hence Finally, the population of infectious mosquitoes in the M I (t) class is generated by progression of mosquitoes in the M E (t) class to M I (t) class at the temperature dependent rate σ M (T ). It decreases due to natural death and the use of adulticides at the rate µ V (T ) 1 + α A A . This gives Thus, the above formulations for human and mosquito dynamics are represented by the following nonautonomous deterministic system of non-linear differential equations (a flow diagram of the model is depicted in Fig. 1 and the state variables and parameters of the model are described in Tabs. 1 and 2): (2.1) For simulation purpose, a generalized temperature function given by will be used, where T 0 is the mean annual temperature, T 1 represents the variation about the mean, ω measures the periodicity of the function and φ is the phase shift of the function. Therefore the time dependent temperature T = T (t), the temperature dependent parameters φ A (T ), σ A (T ), µ A (T ), µ V (T ), a M (T ) and σ M (T ) are continuous, bounded, positive and ω−periodic functions. That is they belong to L ∞ + (0, ω, R + ). Using similar argument to those in [10-12, 18, 25, 29] (and some of the references therein) and the basic fact that for mosquito-borne diseases (such as malaria), the total number of bites made by mosquitoes must equal the total number of bites received by humans. Thus, for the number of bites to be conserved, the following equation must hold Thus, the force of infection in human populations is now given by The model extends (in certain ways) some malaria transmissions models in the literature, such as those in [1-3, 17, 22, 25, 28, 31, 34, 43], by inter alia: (I) Assuming that recovered individuals do not revert to wholly-susceptible class because they enjoy reduced susceptibility to new malaria infection [2,3]; (II) Incorporating vaccination and the use of treated bed nets in humans (this was not considered in [2,3,31]); (III) Including both the use of larvicides and adulticides in mosquito populations (this was not considered in [1-3, 17, 22, 25, 28]); (IV) Dividing human population into compartments based on malaria infection in line with their vaccination status (this was not considered in [2,3,31]); (V) Considering a reduced disease induced death rate, faster recovery rate and slower waning of immunity for vaccinated humans, i.e δ V ≤ δ U , γ V ≥ γ U and τ V ≤ τ U , this was also not considered in [1-3, 17, 22, 25, 28, 34]; (VI) Incorporating the effect of endemicity of malaria by differentiating wholly susceptible from susceptible with prior infection, this was not considered in [2,3,22,25,28,33,34,43].

Temperature dependent parameters
Temperature is known to directly affects vector borne diseases in host; insects are poikilothermic and hence their internal temperature is greatly influenced by environmental temperature, which affect their physiology, as well as exposing the pathogen they carry to environmental temperature [34]. Using the formulations in [6,27,31,36], the temperature dependent parameters of malaria model (2.1) are defined either as Briere or Quadratic functions as follows. The temperature dependent: 1. Oviposition rate (φ A (T )) of mosquitoes is given by

Maturation rate (σ A (T )) of aquatic mosquitoes is given by
3. Death rate of aquatic mosquitoes is obtained from the formulation in [36] as 6. Finally the temperature dependent vector competence defined as the product of the proportion of the bites by infective mosquitoes that infect susceptible humans and the bites by susceptible mosquitoes on  infectious humans that infect susceptible mosquitoes [27] is given by The graphical representations of the temperature dependent parameters are presented in Figures 2, 3, 4 and 5.

Basic properties of the model (2.1)
The basic dynamical properties of the non-autonomous system given by (2.1) is explored. Adding the first ten equations of the model (2.1) gives  Since Thus, a standard comparison theorem can be used to show that , µ V (T )}, and c M = min{α L L , α A A }, the total mosquito population (in both aquatic and non-aquatic stages) satisfies Using similar approach to that of [25] and Gronwall's lemma, the mosquito population has a globallyasymptotically periodic solution satisfying Also from [25], it is assumed that the mosquito population stabilizes at a periodic state, thus for the continuous, bounded, positive and ω-periodic Lemma 2.1. Consider the non-autonomous model (2.1) with non-negative initial condition for all t ≥ 0. Then for any x ∈ C([0], R 14 + ), the model has a unique non-negative solution through x that is ultimately bounded and uniformly bounded.
It is clear that for all x ∈ C([0], R 14 + ), G(t, x) is continuous, and Lipschitzian in x, in addition, > 0 for all i = 1, 2, . . . 14 whenever x i ≥ 0 and x j = 0, thus existence, uniqueness and positivity of solution through (0, X) is guaranteed [25]. Following the approach of [14,15,31] and the fact that B(X) is Metzler and Z > 0, the system is positively invariant in C([0], R 14 + ). Moreover it follows from (2.5) and (2.7), that lim sup Hence all solutions for the system given by (2.1) are ultimately and uniformly bounded [25].

Analysis of the autonomous model
In this section, we analyse the dynamics of the autonomous form of the model (2.1). That is the case when the model parameters are temperature independent. Thus We first of all analysed mosquito-only model in the absence of interaction with humans for its basic dynamical features.

Mosquito-only population model
In this section, we carry out analysis of mosquito-only population model in the absence of interaction with humans and no intervention (adulticide or larvicide) is apply. In the absence of humans, the model (2.1), reduces to the following mosquito-only system: (3.1) The system (3.1) has a threshold quantity called the basic offspring number denoted by N 0 , given by The threshold quantity (N 0 ) can be interpreted as follows: The average time spent by mosquito in the aquatic stage is given by 1/(σ A + µ A ), where σ A is the rate at which aquatic mosquitoes develop into an adult mosquito, so that the probability that an aquatic mosquito develop into an adult female mosquito is given by In the absence of disease, the average life expectancy of an adult female mosquito is given by 1 µ V , so that the average eggs laid by an adult female mosquito throughout her life span is given by Thus, the product of (3.3) and (3.4) gives the basic offspring number of the mosquito-only population model. The mosquito-only model (3.1) has two equilibria depending on N 0 . If N 0 ≤ 1, then, the system (3.1) has only the trivial equilibrium called an extinction equilibrium (E 0 ) given by If N 0 > 1, then, the system (3.1) has a non-trivial equilibrium given by It is worth mentioning that the trivial equilibrium (E 0 ) is biologically less attractive since mosquitoes goes extinct in the population. Rewriting the mosquito-only model (3.1) in the formẋ = f (x), where Ω * ⊆ R 2 + and f : Ω * → R 2 + is continuous. Then we have the following results (The proof is given in Appendix A). Theorem 3.1. The extinction equilibrium (E 0 ) is globally-asymptotically stable (GAS) when N 0 ≤ 1 and unstable otherwise. The non-trivial equilibrium (E 1 ) exists and is locally-asymptotically stable (LAS) when N 0 > 1.
The epidemiological implication of Theorem (3.1) is that if the basic offspring number can be brought to a value less than unity, then mosquito population will goes to extinction and the diseases dies out in time (since no horizontal transmission).

Disease-free equilibrium
The disease-free equilibrium (DFE) is the steady-state solution of the autonomous system (form of model (2.1)) obtained in the absence of disease. The autonomous form of (2.1) has two disease-free equilibria depending on the magnitude of N 0 . Suppose N 0 ≤ 1 and the diseased compartments are zero, then the model has a mosquito extinction DFE, E 2 , given by , , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 . (3.5) If N 0 > 1, then the model has a non-mosquito extinction DFE, E 3 given by , (3.6) Notice that E 2 is biologically less attractive (due to absence of mosquitoes), thus we concentrate on E 3 . The local stability of E 3 can be established using the next generation method [39]. Let It follows then that the vaccinated reproduction number, denoted by, R 0V , is given by Following Theorem 2 of [39], the following result is established.
The threshold quantity R 0V , is the vaccinated reproduction number of the disease. It represents the average number of secondary malaria cases that one infected case can generate if introduced into a population where fraction are vaccinated and the aforementioned control strategies (bed nets, adulticides and larvicides) are used. It can be interpreted as follows.
Infection in humans occurs either in the class of non-vaccinated susceptible or vaccinated susceptible. Susceptible humans acquire infections following effective contact with infectious mosquito in M I class. The number of new human infections in S U class generated by an infectious mosquito in the M I class is the product of the infection rate of infectious mosquito ( β HV a M (1− B α B ) N * H ), probability that a mosquito survives the exposed class (M E ) and move to infected class ( σ M K11 ) and the average duration in the infectious stage ( 1 K10 ). Thus the average number of new human infection generated in S U class by an infectious mosquito (near DFE) is given by (noting . (3.10) Similarly, the number of infections generated by an infected mosquito in the S V class is the product of the infection rate of infected mosquito , the probability that a mosquito survives the exposed class and move to infected class ( σ M K11 ) and the average duration in the infectious stage ( 1 K10 ). So that the average number of new human infection generated in S V class by an infectious mosquito (near DFE) is given by .
, the probability that an individual survives the exposed class E U and move to infectious ) and the average time spent in I U class ( 1 K5 ). Furthermore, the number of new mosquito cases generated by non-vaccinated recovered humans in R U class (near the DFE) is the product of the infection rate , the probability that an individual survives E U and move to I U class ( σ U K3 ), the probability that an individual survives I U and move to R U class ( γ U K5 ) and the average duration in R U compartment ( 1 K7 ). Thus, the average number of new mosquitoes infection generated by non-vaccinated infectious humans (infected or recovered) is given by (noting that M * (3.12) Using similar argument as above, with η I and η V accounting for reduction in infectivity of vaccinated infectious and vaccinated recovered humans in comparison to non-vaccinated infectious humans, the average number of new mosquitoes infections caused by infectious humans in I V or R V compartments is given by Thus, the total number of infections generated by infectious mosquito in the non-vaccinated humans is given by the product of (3.10) and (3.12) denoted by (3.14) Similarly, the total number of infections generated by infectious mosquito in the vaccinated humans is given by the product of (3.11) and (3.13) represented by Thus, the square root of the sum of (3.14) and (3.15) represented by (3.8) gives the vaccinated reproduction number.

Type reproduction number
For a homogeneous system, the vaccinated reproduction number can be seen as the control threshold required to eliminate the disease from a community. The case is different in the case of multiple host types. The typereproduction number (T) is defined as where I is an identity matrix, P is a projection matrix, e is a unit vector with all elements equal to zero except the ith term and K = F V −1 is the next generation matrix with large domain. The type reproduction number correctly determines the critical control effort for heterogeneous populations [23]. We should note from the next generation matrix (K = F V −1 ) that new infections occur only in compartments E U , E V and M E , and therefore it can not be used to compute the type reproduction number for other infected/infectious compartments without new infection. Let the type reproduction numbers of compartments E U , E V and M E be respectively denoted by T 1 , T 2 and T 3 . From (3.16), it can be shown that Similarly, and,

Endemic equilibrium
Let the endemic equilibrium (the case when λ H = 0 and λ V = 0) of the model (2.1) be denoted by E * * = so that, and,

Backward bifurcation
Here we apply the method described in [9,39] which is based on the use of Centre Manifold Theory to prove the existence of backward bifurcation for the autonomous version of model (2.1). We claim the following result (The proof is given in Appendix B): The epidemiological implication of the phenomenon of backward bifurcation is that the classical requirement of R 0V < 1 is, although necessary, no longer sufficient for disease elimination.

Non-existence of backward bifurcation
It is well known that disease induced death rate and imperfect vaccination are some of the major causes of backward bifurcation in vector borne disease models. Let consider the case where disease induced mortality is negligible (δ U = δ V = 0) and the vaccine is perfect (ω V = 0), then the bifurcation coefficient given by a in equation (B.6) reduces to which is less than zero provided W ≥ max

Analysis of the non-autonomous model
Consider the non-autonomous model given by (2.1) with T (t) = T 0 1 + T 1 cos 2π 365 (ωt + φ) as defined in (2.2) and the time dependent basic offspring number given by is strictly greater than 1. To find the disease-free state of the system given by (2.1), we let E U (t) = E V (t) = I U (t) = I V (t) = R U (t) = R V (t) = M E (t) = M I (t) = 0 and obtained a non-trivial disease-free state given by , which is obtained when N 0 (t) > 1. On the other hand, a unique positive trivial non-periodic solution is obtained when N 0 (t) ≤ 1.

Basic reproduction ratio
The local-asymptotic stability of the positive periodic disease-free state (E 4 ) can be established using a threshold parameter called the basic reproduction ratio [40]. We computed the basic reproduction ratio for the model (2.1) using the theory developed in [40] and used for many periodic models such as [1,2,25,28,31,43] and some of the references therein.
Consider the disease compartments of the model (2.1) given by (where T (t) is as defined in (2.2)) where λ H (T ) and λ V (T ) are as defined earlier. The matrix of new infection terms F (t) and that of transfer in and out of infectious compartments V (t) are respectively given by where, T , then the linearization of (5.2) can be re-written in the form Following the approach of [25,40], let Y (t, s) and Φ T = Y (t, 0) respectively be the evolution operator and monodromy matrix of the linear ω-periodic system dy dt = −V y(t), t ≥ s, that is for each s ∈ R, the 8 × 8 matrix Y (t, s) satisfies where I is the identity matrix of order 8. Let Z T be the Banach space of all ω-periodic functions equipped with the maximum norm and an ω-periodic function of s denoted by α(s) be the initial distribution of infectious individuals in the community, then the rate at which new infections are produced by an infected individual in the community who were introduced at time s is given by F (s)α(s) [2,25,40]. Likewise the distribution of new infected individuals from infections at time s and remain in the infected compartments at a later time t is Y (t, s)F (s)α(s). Therefore gives the cumulative distribution of new infections at time t that are produced by all infected individuals (α(s)) introduced at sometimes before t. Define the linear operator L : Suppose ρ(L) is the spectral radius of L, then the basic reproduction ratio (R 0T ) is given by ρ(L) [40]. It is easy to show that, in addition to Assumptions A1 to A5 of [40] satisfied by the autonomous system, the non-autonomous model (2.1) can be shown to satisfy the additional Assumptions A6 and A7 of [40]. Thus the following stability result is obtained.

Effect of control strategies
In this study we consider five main control strategies, namely:  Table 2 (unless otherwise stated), the model will be simulated to assess the effectiveness of these strategies (implemented singly or in combination). Since we are interested in exploring the feasibility of disease elimination, these simulations are carried out for the special case of the model where backward bifurcation does not occur. In order to analyse the effect of the aforementioned control strategies in the presence of temperature changes, three different sets of numerical simulations are carried out, that is the cases where temperature (T ) is 20 • C, 25 • C and 30 • C.

Bed nets-only strategy
Here, the effect of using bed nets-only is assessed, by setting all other parameters related to vaccination and mosquito control to zero, that is the case when ξ V = V = α A = A = α L = L = 0. The model is simulated using the following levels of bed nets effectiveness (with bed nets efficacy of 0.5): (i) Low bed nets effectiveness: α B = 0.1 (i.e., only 10% of individuals uses bed nets effectively); (ii) Moderate bed nets effectiveness: α B = 0.3 (i.e., only 30% of individuals uses bed nets effectively); (iii) High bed nets effectiveness: α B = 0.5 (i.e., only 50% of individuals uses bed nets effectively).
As expected, an increase in the rate of bed net use leads to a decrease in the number of malaria cases. For instance, the resulting number of infected individuals corresponding to the low, moderate and high bed nets use (for the case when the temperature is taken to be 20 • C) is 18,000, 12,200 and 6,100, respectively (Tab 3). Figure 8 shows the simulations of the model using different temperature levels and bed nets use, from which it is evident that the use of bed nets is more effective when temperature is 20 • C (as in Fig. 8A), where the total number of infected humans is lower and approach the DFE faster. The DFE is reached by the total infected humans when T = 30 • C (as shown in Fig. 9C) at almost half the time taken to reach the DFE when T = 25 • C (where the total infected humans are at their peak when T = 25 • C; Fig. 8B). Figure 9D depicts the cumulative number of new cases in humans with low, medium and high rates of applying bed nets when T = 20 • C.
A contour plot of the reproduction number R 0V (with δ U = δ V = ω V = 0, so that backward bifurcation does not exist) as a function of rate of bed nets applications (α B ) and the vaccine efficacy ( B ) is depicted in Figure 23. As expected, the plot show decrease in R 0V values with increasing values of α B and B . Furthermore, based on the parameter values used in the simulation as given in Table 2, a high use of bed nets and highly effective bed nets are required to reduce the associated reproduction number to a value below unity (at least 60% each).

Vaccination-only strategy
The   Similar to the use of bed nets, vaccination is more effective in reducing the total number of infected humans and time taken to reach the DFE when the temperature is 20 • C (Fig. 10A), where as similar effect for both temperatures of 25 • C (Fig. 10B) and 30 • C (Fig. 11C) are obtained. A cumulative number of new cases in humans with low, medium and high rates of vaccine application at T = 20 • C is depicted in Figure 11D.    vaccinated, the disease will still persist in the population, unless the vaccine efficacy is at least 50% (to reduce the vaccination reproduction number to a value below unity).

Mosquito control-only strategy (adulticides strategy)
Here, the effect of the use of adulticides only with efficacy of 0.5 is simulated, that is when α B = B = α L = L = ξ V = V = 0. Notice that the use of larvicides has less or no effect in the dynamics of infected mosquitoes (it mainly affects the basic offspring number).
This control strategy shows the biggest positive effect in both reducing the total number of infected individuals and the effect of using different levels of controls for different temperatures compared to other forms of single controls. For T = 20 • C (as shown in Fig. 12A), mosquito control measure at the rate α A = 0.5 pushes the total  number of infected humans to reach the DFE at the shortest period in comparison to when T = 25 • C (as in Fig. 12B), and when T = 30 • C (Fig. 13C). The cumulative number of new cases in humans with low, medium and high rates of applying adulticides at T = 20 • C is depicted in Figure 13D.
A contour plot of the reproduction number R 0V (with δ U = δ V = ω V = 0) as a function of rate of adulticides applications (α A ) and the efficacy of adulticides ( A ) is depicted in Figure 25. The figure show that, based on the parameter values used, a high rate of application of adulticides with high efficacy rate (at least 85% each) is required to reduce the associated reproduction number to a value below unity, so that malaria could be control from the community.

Bed nets and vaccination strategy
The combined effect of the use of bed nets and imperfect vaccination are explored here. Vaccine efficacy of 0.75 and bed net efficacy of 0.5 are used. For this simulation, parameters related to other control strategies are set to zero (i.e., A = α A = L = α L = 0).
As expected, the combined use of bed nets and vaccination is more effective than the singular use of bed nets or vaccine. Figures 14 and 15 depict the simulations of the model (2.1) with the use of bed nets and vaccination, from which it is evident that the total infected humans reach the DFE (for low, medium and high application of controls) faster than the separate use of the controls for both T = 20 • (Fig. 14A), T = 25 • C (Fig. 14B), and T = 30 • C (Fig. 15C). The cumulative number of new cases in humans with low, medium and high rates of applying bed nets and vaccines at T = 20 • C is depicted in Figure 15D.

Hybrid strategy
The potential impact of using all controls are examined, using the following effectiveness levels:  Table 2 with δ U = δ V = ω V = 0.  Table 2 with δ U = δ V = ω V = 0. For high levels of intervention, the total number of infected humans is lower and reach the DFE at a very short time for the different temperature levels used. For instance, the total infected humans reach the DFE in less than 1,000 days when T = 20 • C as depicted in (Fig. 16A) and it reaches the DFE in about 3,000 days when T = 25 • C (Fig. 16B) and T = 30 • C (Fig. 17C). The cumulative number of new cases in humans with  Table 2 low, medium and high rates of applying bed nets, vaccines and adulticides at T = 20 • C is depicted in Figure  17D.

Sensitivity analysis and numerical simulation
In this section, the partial rank correlation coefficient (PRCC) of the parameters of the vaccinated reproduction number for different temperature ranges are given. In addition, numerical simulations of the model are also presented.

Sensitivity analysis
The partial rank correlation coefficient is a sampling based sensitivity index that measures the strength of the linear associations between a dependent variable (in this case the vaccinated reproduction number), and independent variables (its parameters) after removing the linear effect of other parameter values. We consider the cases for constant and various temperatures. Using the vaccinated reproduction number as the output for temperature values for ranges 15 • C−20 • C, 20 • C−25 • C, 25 • C−30 • C, 30 • C−35 • C and that of constant parameter values with a confidence interval of 95% and 1000 number of boots, the PRCC are obtained. Tables 4 -8 show the PRCC values, bias, standard error, minimum and maximum confidence interval for each of the 27 parameters of the basic reproduction ratio with the aforementioned temperature intervals and constant temperature, respectively. Similarly, Figures 18, 19 and 20 show the bar plot of the PRCC of the parameters of the basic reproduction ratio as temperature varies and constant temperature. For the different temperature ranges, the rate of vaccination ξ V , efficacy of vaccination V and rate of successful use of bed nets C B show little variation as temperature varies. The rate of vaccination is positively correlated to R 0V due to the imperfect vaccine (not 100% effective), thus vaccinated individuals can still acquire infections.
The use of larvicides is positively correlated to R 0V for temperature ranging from 15 • C to 20 • C with PRCC value of +0.4099, but when the temperature range is from 20 • C to 25 • C, the PRCC becomes negatively correlated to R 0V with value of −0.0179, it remains negative with PRCC value of −0.1008 when the temperature is between 25 • C-30 • C and returns to positive for temperature of 30 • C-35 • C having PRCC of +0.1099, where as Disease induced death rates, which have been shown to be main causes of backward bifurcations in mosquito borne diseases such as [10,12,18,20,31]

Numerical simulations
We fix temperature values for the non-autonomous model (2.1), so that each of the temperature dependent parameter becomes constant. For simulation purposes, two different sets of parameter values that give R 0V = 0.297 < 1, N 0 = 14.2 and R 0V = 6.457 > 1, N 0 = 104.4 were chosen within the ranges given in Table 2.   In order to simulate the effect of seasonal variation, we use the generalized temperature function given by (2.2). Figures 6 and 7 show the solution profile of the model (2.1) for the total number of infectious humans (I U + I V ) when R 0V < 1 and when R 0V > 1, respectively. The solution profile of total infectious humans approach the DFE when R 0V < 1 (Fig. 6) and approach an endemic equilibrium when R 0V > 1 (Fig. 7). Similar results hold for the non-autonomous system for the case when R 0T < 1 and R 0T > 1 as presented in Figures 21A and 21B, respectively. Figure 22 shows effect of vaccination on disease incidence in humans with seasonal variation. The figure show that the oscillation pattern differs between the vaccinated and the non-vaccinated infectious (Fig. 22A) and exposed (Fig. 22B) humans, both in their subharmonic periods and the relative phase of cycles.

Conclusion
A new mathematical model to assess the effect of temperature on control strategies of malaria, is constructed and analysed. Some of the main findings of this study are summarized below: (I) The autonomous mosquito-only model has a threshold quantity called the basic offspring number with the property that, if the threshold quantity (N 0 ) is less than or equal to unity, the mosquito population goes to extinction, and it establishes if (N 0 ) > 1. (II) The autonomous model version of system (2.1) has two disease-free equilibria, the mosquito-extinction equilibrium (E 2 ) which is globally-asymptotically stable (GAS) when the basic offspring number (N 0 ) is less than unity and the non-mosquito-extinction equilibrium (E 3 ) which is locally asymptotically stable when R 0V ≤ 1. (III) The autonomous model undergoes the phenomenon of backward bifurcation, which could be removed for a special case when malaria induced death rates (δ U = 0 and δ V = 0) and the vaccine waning rate are negligible (ω V = 0). (IV) Relationship between the vaccinated reproduction number and the type reproduction numbers is established, where it is shown that T i < 1 (i = 1, 2, 3), provided R 0V < 1 (and T i ⇔ R 0V ). This result suggest that, malaria can be control by targeting certain groups in the population. (V) The non-autonomous model (2.1) has a disease-free equilibrium (E 4 ), which is shown to be locallyasymptotically stable whenever the associated reproduction ratio is less than unity and unstable otherwise.  (VIII) For non-control parameters, the most positively correlated parameters within all ranges are the mosquito carrying capacity (K), probability of disease transmission (β HV ), reduction parameter in the transmission of infected vaccinated humans and rate of vaccination (ξ V ). On the other hand, recovery rates γ U , γ V and human recruitment rate (Π H ) are the most negatively correlated in all the temperature ranges. (IX) Numerical simulation of the model, using appropriate demographic and epidemiological data for Kwazulu Natal province of South Africa, show that (for the case where backward bifurcation does not occur), the hybrid strategy which combines all the strategies (that is combined use bed nets, vaccination and adulticides) is more effective than singular use of the aforementioned control strategies, that can lead to malaria elimination in the province. (X) Seasonal variation may cause fluctuation in the number of people becoming infected with malaria in both the vaccinated and non-vaccinated populations. (XI) Simulations of the model show that high vaccine efficacy is required to reduce the vaccinated reproduction number to a value below unity. Further, a singular effective use of bednets can result in effective control of malaria in a community provided the bed net coverage and the bed net efficacy are high enough (at least 60 % each). so that f (b) < 0, provided N 0 ≤ 1.
Therefore f (b) ≤ 0 ≤ f (0), whenever N 0 < 1. Thus the mosquito component of model (2.1) defines a positive dynamical system on [0, b] and E 0 is GAS on [0, b]. But since q is arbitrary, b can be extended and the result holds on R 2 + . The second part of the proof follows straightforward by linearization.
Appendix B. Proof of theorem 4.1 Proof. The proof is based on using centre manifold theory [9,39]. To apply the theorem, we carry out the following changes of variables. Let The transformed malaria model (2.1) is represented by with the associated forces of infection given by By letting R 0V = 1 we have, suppose, further that β HV = β * HV is chosen to be the bifurcation parameter. The Jacobian matrix (J * ) at the DFE with β HV = β * HV is given by