Role of Differential Susceptibility and Infectiousness on the Dynamics of an SIRS Model for Malaria Transmission

: A deterministic model for the transmission dynamics of SIRS-type malaria in hosts and SI in mosquito populations is proposed. The host population is differentiated between naive, primary, and secondary susceptible individuals. Primary and secondary infected individuals (and also recovered) are differentiated from each other according to their degree of infectiousness. The impact of changing the relative susceptibilities of primary and secondary (with respect to naive) susceptible individuals on the dynamics is investigated. Also, the impact of changing the relative infectiousness of secondary infected, primary, and secondary recovered individuals (with respect to primary infected) on the transmission dynamics of malaria is studied.


Introduction
Malaria, a mosquito-borne infectious disease caused by the protozoan parasites of the genus Plasmodium, is transmitted to humans via bites from infected female Anopheles mosquitoes, which introduce the organisms from its saliva into the host's circulatory system, which is potentially fatal.In humans, malaria is caused by Plasmodium falciparum, Plasmodium malariae, Plasmodium ovale, and Plasmodium vivax, with deaths caused primarily by P. falciparum and P. vivax, while P. ovale and P. malariae are typically responsible for milder forms of malaria.Malaria is prevalent around the equator, in tropical and subtropical regions that include much of Sub-Saharan Africa, Asia, and the Americas.The prevalence of malaria is tied in to circumstances connected to rainfall and temperature patterns, altitude, and the availability of appropriate larvae breeding sites (stagnant waters).
Malaria, an ancient disease [1,2], has instigated sustained efforts aimed at ameliorating its impact on populations using what has been learned from scientific research and epidemiological studies, which are frequently complemented with the results of analyses and simulations of mathematical models over the past few decades [3][4][5][6].
Malaria mathematical modeling at the population level essentially started with the work of Sir Ronald Ross [7] and was expanded and applied by Macdonald [8] and many others, including Dietz [9], Aron and May [10], Torress-Sorando and Rodriguez [11], Koella and Anita [12], and Filipe et al. [13].Sir Ronald Ross shared the 1902 Nobel award in malaria through research that connected the life-cycle of the malaria parasite in hosts and vectors [14].His 1911 paper introducing a nonlinear system of differential equations to study the transmission dynamics of malaria was groundbreaking and arguably the paper that established the field that is now known as mathematical epidemiology, which is tightly connected to the study of host-parasite dynamics, disease evolution, and vector control [7].Macdonald extended Ross's model by adding biological realism, with the aim of connecting questions to models and models to data [8,[15][16][17][18].In an early attempt to study the theory of eradication of malaria, Macdonald and his collaborators [17] analyzed the factors affecting the basic reproduction number, R 0 .Moreover, they discussed the epidemiology of malaria [18], with special reference to the pattern observed in equatorial Africa.In [16], Macdonald and Goeckel studied the effect of considering a complete interruption of transmission on the regular progressive decrease in the falciparum parasite over a specific period of time.Macdonald and Goeckel also estimated the slowest acceptable vivax fall-off rate over a 12-month period [16].
Certainly, the research that connected the life cycle of vectors, hosts, and parasites naturally brought exceeding interest to the study of the role of the immune system within hosts and across vectors and human malaria interactions.How would such interactions within and between individuals influence the levels of infectiousness, susceptibility, or re-infection after recovery?The challenges raised in connecting multiple levels of malaria dynamics and infections were probably addressed first by Dietz et al. [9], when these researchers incorporated the role of immunity at the population level.The use of malaria models to study specific questions exploded after the seminal work of Sir Ronald Ross and the innovative extensions carried out by Macdonald [15][16][17][18] and Dietz and his collaborators [9] (see also [19][20][21][22] and the references therein).Some reviews of the mathematical modeling and analysis of malaria are shown in [3,4,10,[23][24][25]. Recent population-level efforts have included the role of genetics, with specific research carried out on the role of sickle cell anemia and the presence of large populations of asymptomatics [26,27].More studies on the impact of climatic change and weather on malaria dynamics are shown in [20,28,29].
Ngwa and Shu [22] considered an SEIRS model for the host and SEI for the vector under the assumption that the population is closed and of varying sizes.It was assumed that recovered individuals are partially immune and yet still capable of transmitting the infection.They derived conditions for the existence of endemic equilibria and studied the global stability of the malaria-free equilibrium for R0 < 1, where R0 is their basic reproduction number.
Chitnis et al. [30] considered an SEIRS model for the host and an SEI model for the vector.The authors considered an open population with recruitment through births and immigrating (susceptible) individuals.Their model supports the existence of an infectionfree equilibrium that is locally asymptotically stable whenever R 0 < 1.In their paper, it was shown that more than one endemic equilibrium may exist if R 0 > 1.In a special case, where the malaria-induced mortality rate vanishes, the bifurcation at R 0 = 1 is forward.
Niger and Gumel [31] considered a model for malaria transmission dynamics with repeated exposure.The model has been analyzed and shown to exhibit the existence of multiple endemic equilibria when R 0 < 1 for certain parameter values, in case of standard incidence.However, the authors showed that if the repeated exposure is ignored and if the incidence function is of a mass-action type, the model is reduced to a simple SIRS model in the host and SI in the vector.The reduced model is shown to not have multiple endemic equilibria.
This paper uses mathematical models to expand on the role of the human immune system's variability in the transmission dynamics of malaria and control.Specifically, we focus on the interplay between differences in susceptibility, ability by humans to increase or reduce the 'health' of the parasite population, and differential levels of infectiousness; that is, the ability of humans to harbor effectively limited or large populations of parasites.The incorporation of heterogenous and variable levels of infectiousness and susceptibility bring unexpected, possibly unwanted, consequences to the forefront.Models that consider variation in susceptibility and infectiousness will invariably support multiple endemic states; that is, they naturally support disease dynamics that depend on initial conditions and that cannot be characterized by the basic reproduction number, a dimensionless quantity that has brought simplicity and clarity to those interested in finding ways of characterizing effectively simple and direct modes of malaria elimination or prevalence reductions.The fact is that variation in host susceptibility and infectiousness, the norm in biology, limits the preponderant role given to the basic reproduction number; the result of the theoretical work carried out by models that omit variation (a most important pre-requisite for selection).
More specifically, we developed a mathematical model where we differentiate between individuals who have never experienced the infection (called naive), those who experienced it only once (called primary), and those who experienced it at least twice (called secondary).The aim of this paper is to study the impact of changing the relative susceptibility and infectiousness of primary and secondary individuals on the disease dynamics in the human population.
The paper is organized as follows.In Section 2, we formulate the model and study the existence and stability of the malaria-free equilibrium.We further compute the basic reproduction number and find an equation that defines the endemic equilibrium human force of infection λ h .In Section 3, we study the direction of bifurcation at R 0 = 1 and investigate the possible bifurcation diagrams in the plane (R 0 , λ h ).In Section 4, we study the impact of changing the relative susceptibilities and infectiousness in the population dynamics of the host and parasites.In Section 5, we study the role of differential mortality.We end the paper with a summary and conclusion in Section 6.

Null Model
In this section, we introduce a null model, which allows us to state basic malaria model results in the absence of variation in infectiousness and susceptibility.It is an SIRS model for human populations and an SI model for mosquito populations.Here, the total human population (with size N h (t)) is subdivided into susceptible, infected, and recovered, whose size at time t is given by S h (t), I h (t), and R h (t), respectively.The mosquito population (of size N v (t)) is subdivided into susceptible S v (t) and infected I v (t), see Table 1 for more detailed description of the physical meaning.Both populations of interaction are assumed to be homogeneously and symmetrically fully mixed.The dynamics of interaction between both (human and mosquito) populations is governed by the model where Λ h , Λ v , µ h , µ v , p h , p v and β hv are as defined in Table 2.The per-capita rate γ h denotes the transition from I h to R h state, while the transition from R h to S h state is denoted by the per-capita rate ν h .The malaria-induced mortality rate (in humans) is denoted by δ h .Recovered individuals are assumed to be weakly infectious.Here, q denotes the infectiousness reductions for Rwith respect to S-individuals.Equation ( 1) is defined on the closed set D , where In the simple case, where deaths due to malaria are neglected (i.e., δ h = 0), Equation (1) supports a malaria-free equilibrium that is stable if and only if the basic reproduction number for Equation ( 1) is less than one.Moreover, Equation (1) has a unique endemic equilibrium that exists only if the basic reproduction number is bigger than one; an equilibrium that is stable whenever it exists.On the other hand, whenever human deaths due to malaria are considered, we see that Equation (1) exhibits backward bifurcation for certain parameter values.The derivation of these results is deferred to Appendix A.
If we extend the above model by further subdividing the epidemiological classes S h , I h and R h into naive susceptible representing those who have never experienced the infection, primary susceptible representing those who experienced it only once, and secondary susceptible to denote those who experienced malaria at least twice, then the dynamics supported by the model are enriched; an expected outcome when individuals are differentiated from each other by their susceptibilities.Similarly, infected humans are subdivided into primary infected, those who experience malaria for the first time, and secondary infected, those who are infected at least for the second time, with individuals in both classes supporting different degrees of infectiousness.Also, individuals in the class of recovered are categorized into primary recovered and secondary recovered.The extended model is shown in the next section.The aim of this paper is to study the impact of variation in susceptibility and infectiousness on the dynamics of malaria.(1) Do these extensions affect endemic behavior?(2) What is the impact of changing the relative susceptibilities of primary and secondary (with respect to naive) susceptible individuals on the dynamics and, therefore, on the elimination effort?(3) What is the impact of changing the relative infectiousness of secondary (with respect to primary) infected individuals?(4) What makes it possible to support multiple endemic equilibria?

Extended Model
To study the impact of differential susceptibility, infectiousness, and mortality on the dynamics of malaria (at the population level), we consider a model of type SIRS for humans and SI for mosquitoes.In such a model, individuals are epidemiologically asymmetric, but they are symmetric in their mixing behavior in the sense that they mix homogeneously so that each mosquito bite has an equal chance of transmitting the virus to susceptible humans in the population.The total host (human) population at time t, whose size is denoted by N h (t), is subdivided into seven mutually exclusive sub-populations: naive susceptible S 0 (t), primary infected I 1 (t), primary recovered R 1 (t), primary susceptible after experiencing the infection once before S 1 (t), secondary infected I 2 (t), secondary recovered R 2 (t), and secondary susceptible S 2 (t), so that N h (t) = S 0 (t) + I 1 (t) + R 1 (t) + S 1 (t) + I 2 (t) + R 2 (t) + S 2 (t).Naive susceptible human individuals are those who never experienced a malaria infection, while primary infected human individuals are hosts who acquired a malaria infection for the first time and are capable of transmitting malaria to susceptible mosquitoes while the mosquitoes are taking their meals from human individuals.Primary (secondary) recovered human individuals are human hosts who experienced malaria infection only once (at least twice) and became immunized, but are still slightly infectious [32].Primary (secondary) susceptible individuals are susceptible human hosts who experienced malaria only once (at least twice) and recovered with partial immunity, in the sense that they are less susceptible to acquire malaria than naive susceptible ones.Secondary infected individuals are human individuals who acquired a malaria infection at least twice and are capable of transmitting it to susceptible mosquitoes, with differential transmissibility of the infection as opposed to primary infected ones.Table 1 shows briefly the physical meaning of the model state variables.
Similarly, the total vector (mosquitoes) population at time t, whose size is denoted by N v (t), is split into two sub-populations: susceptible mosquitoes S v (t) and infectious mosquitoes I v (t), so that N v (t) = S v (t) + I v (t).
It is supposed that new recruits of the host enter the naive susceptible human subpopulation class with recruitment rate Λ h and that deaths occur at the natural death per-capita rate µ h in all human sub-population classes except those, for primary and secondary infected sub-populations, where excess deaths due to malaria are included.The malaria-induced mortality for primary and secondary infected individuals occurs at rates δ 1 and δ 2 , respectively.Naive susceptible humans either die of natural death (µ h ) or acquire a malaria infection with force of infection λ h (t), to become primary infected individuals who either die or recover to be primary recovered individuals with rate γ 1 .Therefore, that rate of change in the number of naive susceptible and primary infected humans is described by two ordinary differential equations Human individuals are supposed to stay primary recovered for a period of time 1/ν 1 before losing their acquired immunity to become primary susceptible, and then die naturally or become infected with rate r 1 λ h (t) to be secondary infected, where r 1 is the relative susceptibility of primary with respect to naive susceptible individuals.Therefore, we have the following subsystem of equations It is worth explaining that the parameter r 1 (r 2 ) is a rescaling dimensionless parameter that accounts for the reduction in the susceptibility of individuals due to enhancing their immunity after experiencing a malaria infection only once (at least twice).
Secondary infected individuals either die or recover with rate γ 2 to become secondary recovered, who are supposed to either die naturally or lose their immunity to become secondary susceptible.The rate at which secondary recovered individuals lose their immunity is ν 2 .Finally, secondary susceptible individuals either die or become infected again with force of infection r 2 λ h (t), where r 2 is the relative susceptibility of S 2 compartment individuals with respect to naive susceptible individuals.Therefore, we have the following sub-dynamical system On the other hand, susceptible mosquito populations recruit with rate Λ v and either die naturally with rate µ v , or become infected with force of infection λ v (t), while infected mosquitoes are assumed to die naturally with the same rate µ v .A schematic flow diagram for the model states is depicted in Figure 1, while a full description for the model variable states and parameters is given in Tables 1 and 2, respectively.Following the same approaches shown in [31,33], the incidence rates λ h and λ v are given by where the dimensionless parameters q 1 , q 2 and q 3 account for the reduction in the infectiousness of human hosts who, respectively, acquired a malaria infection at least twice, recovered primarily from their malaria infection (and are assumed to be asymptomatic carriers), and recovered after experiencing the infection at least twice (assumed to be an asymptomatic carrier too), in comparison to individuals who were infected only once.
S 0 (t)    The relative susceptibility of secondary susceptible with respect to naive susceptible individuals Dimensionless The relative infectivity of secondary infected with respect to primary infected individuals Dimensionless The relative infectivity of primary recovered with respect to primary infected individuals Dimensionless The relative infectivity of secondary recovered with respect to primary infected individuals Dimensionless Collecting the above formulated equations together, we arrive at the following mathematical model of ordinary differential equations where λ h and λ v are given by Equations ( 3) and ( 4), respectively.Moreover, all parameters are assumed to be strictly positive, with the exception of the parameters δ 1 and δ 2 , representing the malaria-induced human mortality rates, which are assumed to be nonnegative.
Following the same approach shown in [33], it can be shown, for Equation ( 5), that the closed set where is positively invariant, and the Equation ( 5) is mathematically and epidemiologically wellposed, in the sense that solutions starting in D remain in it all the time.

Existence and Stability of the Infection-Free Equilibrium and the Basic Reproduction Number
To find the equilibria, we put the derivatives in the left-hand side in Equation ( 5) equal to zero.The malaria-free equilibrium (MFE) is obtained by setting and solving the resulting system.Therefore, the MFE is where denotes vector transpose.

The Basic Reproduction Number
The basic reproduction number is obtained using the next generation operator method [34].Accordingly, we evaluate the matrices F (for the new infection terms) and V (for the transition terms) as It follows then that the basic reproduction number is given by where ρ is the spectral radius (dominant eigenvalue in magnitude) of the matrix FV −1 .It is clear that the basic reproduction number R 0 depends on the relative infectiousness of secondary infected (with respect to primary infected) individuals and depends neither on the relative susceptibilities r 1 , r 2 nor on the relative infectiousness q 2 , q 3 of recovered individuals.

Stability of the Malaria-Free Equilibrium:
To establish the stability of the malaria-free equilibrium, we linearize the model around it to obtain the Jacobian matrix (evaluated at the infection-free equilibrium) as Using the columns number one, four, seven, and eight, one may compute four eigenvalues −µ h , −µ h , −µ h , −µ v for the Jacobian matrix J 0 .Then, we use row number five to obtain an eigenvalue of −D 21 = −(δ 2 + γ 2 + µ h ).Thereafter, we use row number six to obtain an eigenvalue of −D 22 = −(µ h + ν 2 ).All these six eigenvalues are negative.Thus, the remaining eigenvalues of J 0 are those of the sub-matrix Hence, the stability of the Malaria-free equilibrium E 0 is determined by studying the roots of the characteristic equation of J sub that is given by where ξ is the eigenvalue of J sub .On simplifying and collecting terms together, we arrive at a cubic equation in the form where To this end, we apply the Routh-Hurwitz criterion [35] .We notice that A 3 > 0 and Thus, with the use of Equation ( 9), we deduce that A 1 > 0 and A 0 > 0, for R 0 < 1.We further notice that However, Hence, with the use of Equation ( 9), we deduce that both T 2 and T 3 are positive.Hence, based on the Routh-Hurwitz criterion [35], all roots of Equation ( 8) have a negative real part for R 0 < 1.
We end this subsection by stating the following proposition on the local stability of the malaria-free equilibrium.
Proposition 1.The malaria-free equilibrium is locally stable if and only if R 0 < 1, while it is unstable otherwise.

Endemic Equilibria
Persistent solutions (endemic equilibria) are solutions where the infection is assumed to persist, in the sense that the prevalence of infection at equilibrium is positive.So, let us assume that the equilibrium human and vector force of infections λ h = 0 and λ v = 0. Hence, on setting the derivatives in the left hand side of Equation ( 5) equal zero and solving the resulting system with respect to the state variables, one can obtain at equilibrium and On using Equation (18) into Equation (3), we get at equilibrium Now, we use Equation (20) into Equation (4) to get where, with the use of Equations ( 11), ( 12), (14), and ( 15) Thus, on using Equations ( 19) and (22) in Equation ( 21), we get On simplifying Equation ( 23), we get a nonlinear algebraic polynomial equation of degree six in the endemic host's force of infection λ h .Its term-free of λ h is while the coefficients of the other six terms with λ h have very complicated mathematical expressions.

Center Manifold Analysis near the Malaria-Free Equilibrium
The model we study has an endemic equilibrium equation with a polynomial of degree six in the endemic force of infection λ h .This polynomial has six solutions.However, of interest are only feasible solutions, in the sense that these solutions are real and positive.Also, it is more interesting if these feasible solutions exist for values of 0 ≤ R 0 < 1.
On setting λ h = 0 in the equilibrium Equation ( 23), we get β hv = β 0 hv , where This value corresponds to R 0 = 1.Thus, in the plane (β hv , λ h ), there is a bifurcation point given by (β 0 hv , 0), at which we compute the direction of bifurcation.Using the theorem 4.1 of Castillo-Chavez and Song [36], the conditions ensuring the existence of backward bifurcation have been derived and are shown in Appendix B. The computations show that Equation (5) undergoes backward bifurcation at R 0 = 1 if and only if K > 0, where On performing a sensitivity analysis of the expression K with respect to the parameters r 1 , q 1 , q 2 , q 3 , δ 1 , and δ 2 , we find that q 2 , r 1 and q 3 are the most sensitive parameters.Namely speaking, q 2 is more sensitive than r 1 , which in turn is more sensitive than q 3 .Therefore, the inequality a 1 > 0 required for the existence of backward bifurcation is expressed as Figure 2 shows the area in the planes (q 2 , r 1 ) and (q 3 , r 1 ) for which the bifurcation is subcritical (backward) as well as that in which it is supercritical (forward).It shows that the possibility for the existence of backward bifurcation increases with the increase in q 3 , while it decreases as q 2 increases.Based on the above analysis, we summarize our results in the following proposition.Proposition 2. Equation (5) exhibits backward bifurcation if and only if condition (26) holds.

Possible Bifurcation Diagrams
Equation ( 23) could be seen as a function in the contact rate β hv .Since the expressions F 1 (λ h ), F 2 (λ h ), and F 3 (λ h ) do not depend on β hv , we can rewrite Equation (23) in terms of where The area for backward as well as forward bifurcation.Part (a) shows these regions in the plane (q 2 , r 1 ) for several values of q 3 .It shows that the backward bifurcation region extends with the increase in the q 3 value.However, part (b) shows that, in the plane (q 3 , r 1 ), the backward bifurcation region shrinks with the increase in the q 2 value.Simulations have been done for the parameter values given in Table 3. 1.2 assumed q 2 0.4 assumed q 3 0.2 assumed It is easy to check that D 21 D 22 − γ 2 ν 2 > 0. Therefore, F 3 (λ h ) > 0 and F 2 (λ h ) > 0. Since N h > 0, then F 1 (λ h ) > 0. Thus, all coefficients C0 (λ h ), C1 (λ h ) and C2 (λ h ) are positive.Hence, Equation (27) has the unique positive solution Assume now that the equilibrium human force of infection λ h is known.In order to study the possible bifurcation diagrams, we investigate the number of feasible turning points.By a feasible turning point we mean a turning point that occurs at a positive value of the force of infection λ h .To do so, we study and count the number of possibly feasible solutions of the equation If Equation ( 29) has no feasible solution, then the model shows only forward bifurcation, Figure 3a.However, if it has a unique feasible solution, then it shows a backward bifurcation, Figure 3d.If it has two feasible solutions, then the model shows the existence of hysteresis, in the sense that it shows multiple super-and subcritical endemic states where the bifurcation curve looks like those shown in Figure 3b,c.3, except for r 1 = 0.0194 and q 2 , q 3 as stated above.
Since the functions F 1 , F 2 are polynomials of degree 3, while F 3 is a polynomial of degree 2 in the variable λ h , then according to formula (28) the expression under the square root in (28) would be a polynomial of degree 8.Moreover, this expression is added and multiplied with other polynomials and the whole thing is divided by some other polynomials.All of these, along with the Equation ( 29), lead to a complicated nonlinear algebraic (with square root function) equation from which we compute the turning points.It is difficult, if not impossible, to find exact formulae for the turning points and also for the constraints on model parameters which assure the existence of such feasible points.Therefore, numerical simulations have been performed.
For our model, extensive simulations have been made and we found out that the model may show forward, multiple supercritical, multiple subcritical, and backward bifurcations, as we shall see in the next section, where we study the impact of the relative susceptibilities r 1 , r 2 , relative infectiousness q 1 , q 2 , q 3 , and the differential mortalities δ 1 and δ 2 on developing all these kinds of bifurcations.

Role of Differential Susceptibility and Infectiousness
Differential susceptibility and differential transmissibility/infectiousness play a key role not only in the type of bifurcation but also in the shape of the bifurcation curve.In the following we try to, numerically, figure out the role of changing the relative susceptibilities r 1 , r 2 as well as the relative infectiousness q 1 , q 2 and q 3 in the overall bifurcation dynamical behavior.This is assessed by letting the parameter of interest have changing values (in various ways) while fixing the remaining model parameters and noticing the type of the resulting bifurcation diagram.It is worth mentioning that the bifurcation diagram gives information on the number of endemic equilibria and the minimum value of the basic reproduction number under which the malaria infection dies out.In other words, it gives information about the minimum elimination effort required to eliminate the malaria infection [39].We first consider the case that recovered individuals are not capable of transmitting malaria which, mathematically, means that q 2 = q 3 = 0.
4.1.Non-Infectious Recovered Individuals (q 2 = q 3 = 0) In this case, the basic reproduction number R 0 reduces to while the condition (26), on the relative susceptibility of individuals in compartment S 1 , with respect to naive susceptible individuals, for backward bifurcation reduces to Relation (31) shows that the critical level for backward bifurcation r1 c does not depend on the relative susceptibility (of secondary with respect to naive susceptible human individuals) r 2 , but depends on the relative infectiousness (of secondary with respect to primary infected human individuals) q 1 and other model parameters.Thus, for fixed r 2 , the backward bifurcation condition ( 31) is a decreasing function in q 1 , see Figure 4, while for fixed q 1 , it is constant.On the other hand, Equation ( 23), from which we determine the equilibrium human force of infection, remains at degree six.Numerical simulations have been made to detect the progression of the shape of the bifurcation curve with the increase in one of the relative susceptibilities r 1 , r 2 and the relative infectiousness q 1 while fixing the other two parameters, with the remaining model parameter values as shown in Table 3. Bifurcation diagram in the plane (q 1 , r 1 ), for q 2 = q 3 = 0 and with parameter values as in Table 3. High enough values of q 1 and r 1 exhibits the existence of backward bifurcation.

Impact of the Relative Infectivity Ratio q 1
In order to assess the impact of changing the relative infectiousness q 1 of secondary (with respect to primary) infected individuals on the model dynamics during a long-time run, we let q 1 change and fix the other model parameters.In this case, the resulting bifurcation diagram looks like that shown in Figures 3 and 5.For fixed high-enough values of the relative susceptibility of primary susceptible with respect to naive susceptible human individuals r 1 < r1 c , while increasing q 1 , the bifurcation curve progresses from the case where there is exactly one supercritical endemic equilibrium, Figure 5a, to the case where there are two subcritical endemic equilibria but only one supercritical endemic equilibrium (i.e., backward bifurcation), Figure 5b.However, if we choose r 1 to be small enough (e.g., r 1 = 0.0194) and letting q 1 increase (Figure 3), the shape of the bifurcation curve changes from a simple forward bifurcation (Figure 3a) to a backward bifurcation (Figure 3d), passing through the formulation of hysteresis (Figure 3b,c).Bifurcation diagram for q 2 = q 3 = 0.It shows that with the increase in q 1 , while fixing r 1 , r 2 , and all other parameters, the model develops forward (a) to backward (b) bifurcation.Simulations have been done for the parameter values given in Table 3, except for r 1 = 4.0194 and q 2 , q 3 as stated above.

Impact of the Relative Susceptibility Ratio r 1
On the other hand, if we fix the relative infectiousness ratio (of secondary infected relative to primary infected human individuals) q 1 and let the relative susceptibility ratio r 1 increase, the model shows the following two behaviors concerning the bifurcation curve.The first is that for a fixed value of q 1 = 0.5 in a certain range (e.g., q 1 < q1 ), if we let r 1 be small enough ( r 1 = 0.00001 < r1 c = 2.0267), the bifurcation curve represents the situation where there is exactly one supercritical endemic equilibrium (i.e., simple forward bifurcation), Figure 6a.However, if we let r 1 increase (r 1 = 0.005 < r1 c ), then the model exhibits the existence of multiple supercritical endemic equilibria, Figure 6b.If we further let r 1 increase (r 1 = 2 < r1 c ), the model shows again the existence of a unique supercritical endemic equilibrium, Figure 6c.Also, if we let r 1 be bigger than r1 c (r 1 = 2.5), then a backward bifurcation appears, Figure 6d. .Bifurcation diagram for q 2 = q 3 = 0.It shows that for small q 1 = 0.5, while fixing r 2 and all other parameters, but increasing r 1 , the model develops forward (a-c) to backward (d) bifurcation.
Simulations have been done for the parameter values given in Table 3, except for r 2 = 0.7 and q 2 , q 3 as stated above.
Unlike the above case, if we choose the relative infectiousness q 1 to be large, then, with the increase in the relative susceptibility r 1 , the bifurcation curve develops from forward bifurcation (Figure 7a) to forward with multiple supercritical endemic equilibria (Figure 7b), forward with multiple subcritical endemic equilibria (Figure 7c), and then backward bifurcation (Figure 7d).Bifurcation diagram for q 2 = q 3 = 0.It shows that for q 1 = 1.6 and keeping r 2 and all other parameters fixed, except r 1 , which is allowed to increase, the model develops forward (a), multiple supercritical (b), multiple subcritical (c), and backward (d) bifurcation.Simulations have been done for the parameter values given in Table 3, except for r 2 = 0.7 and q 2 , q 3 as stated above.

Impact of the Relative Susceptibility r 2
On fixing the relative infectiousness (of secondary infected with respect to primary infected human individuals) q 1 , the condition (31) for backward bifurcation says that the critical level of the relative susceptibility (of primary with respect to naive susceptible human hosts) r1 c is constant.Thus, for any fixed value of the relative infectiousness q 1 (e.g., q 1 = 1.6), then the critical value of the relative susceptibility r 1 above which the bifurcation is backward would be constant ( r1 c = 0.6333).Here, we have three situations.
• The first one is to consider r 1 to be small enough (e.g., r 1 = 0.001).Here, if we let r 2 increase, then the bifurcation curve would progress as follows.For small values of r 2 (e.g., r 2 = 0.0001), the bifurcation diagram shows the existence of exactly one supercritical endemic equilibrium (similar to that shown in Figure 6a).If we increase r 2 (e.g., r 2 = 0.01), then we get a bifurcation curve for which a forward bifurcation as well as multiple supercritical endemic equilibria are exhibited (similar to that shown in Figure 6b).However, on letting r 2 increase further, the bifurcation curve becomes forward with a unique endemic equilibrium, similar to that shown in Figure 6c.Thus, in this case, multiple subcritical and supercritical endemic equilibria are impossible to exist.

•
The second situation is to consider slightly larger values of r 1 (e.g., r 1 = 0.0194), but still less than r1 c = 0.6333.In this case, for small values of the relative susceptibility (of secondary with respect to naive susceptible human individuals) r 2 (e.g., r 2 = 0.005), the model shows the existence of a unique endemic equilibrium, Figure 8a.However, for larger values of r 2 (r 2 = 0.5), the model exhibits forward as well as multiple supercritical endemic equilibria, but no subcritical endemic equilibria, Figure 8b.If r 2 is increased further (r 2 = 5), the model shows the existence of forward bifurcation as well as sub-and supercritical endemic equilibria, Figure 8c.If r 2 is assumed to further increase, then the critical value of R 0 (at the turning point) separating between nonexistence and existence of endemic equilibria decreases (Figure 8d), which implies an increase in the effort required to eliminate the malaria infection [39].

•
The last situation is to consider values of the relative susceptibility (of primary with respect to naive susceptible human hosts) r 1 > r1 c = 0.6333 (e.g., r 1 = 0.7).In this case, the model exhibits backward bifurcation, but what is important is that increasing the value of the relative susceptibility (of secondary with repect to naive susceptible human hosts) r 2 implies a decrease in the value of the critical R 0 separating between the nonexistence and existence of endemic equilibria, which in turn increases the minimum effort required to eliminate malaria, Figure 8e-h. .Bifurcation diagram for q 2 = q 3 = 0 and for fixed q 1 = 1.6.It shows that for a value of r 1 = 0.0194 < 0.6333 = r c 1 , increasing the relative susceptibility r 2 , while fixing all other parameters, the model develops forward (a), multiple supercritical (b), and multiple subcritical (c,d) bifurcation, but it cannot exhibit backward bifurcation.However, for r 1 = 0.7 > 0.6333 = r c 1 , the model exhibits backward bifurcation where the critical value of R 0 at which the turning point occurs decreases, (e-h).This induces an increase in the elimination effort.Simulations have been done for the parameter values given in Table 3, except for r 1 = 0.0194 and q 2 , q 3 as stated above.
We summarize our results as the following: if the infectivity of recovered individuals is ignored (i.e., q 2 = q 3 = 0), then, 1.
For small enough values of the relative infectiousness q 1 (i.e., for values of q 1 on the left of the dotted vertical line in Figure 9), the model undergoes backward bifurcation as the relative susceptibility r 1 increases.

2.
If the relative infectiousness q 1 of secondary infected individuals is allowed to be fairly increased (e.g., for values of q 1 between the two vertical lines in Figure 9), then (as r 1 increases) the progression of the bifurcation diagram is as follows: a forward bifurcation (with the existence of a unique endemic equilibrium for R 0 > 1) is exhibited.Then, a forward bifurcation with multiple supercritical endemic equilibria (hysteresis) is shown.Then, it exhibits again a forward bifurcation with unique endemic equilibrium for R 0 > 1.Finally, it undergoes backward bifurcation.

3.
If q 1 is further increased (e.g., for values corresponding to higher than that at the vertical dashed line in Figure 9), then the model shows transcritical bifurcation for small values of r 1 , which is followed by forward bifurcation with multiple supercritical endemic equilibria as r 1 increases.If r 1 is further increased, then the model shows the existence of forward with multiple super and subcritical endemic states (hysteresis) and finally, it undergoes backward bifurcation.A schematic bifurcation diagram for the case q 2 = q 3 = 0 in the plane (q 1 , r 1 ).The notation: 1 means a forward bifurcation with unique endemic equilibrium for R 0 > 1; 2 means backward bifurcation; 3 means forward bifurcation with multiple supercritical endemic states only (i.e., multiple equilibria exist only for R 0 > 1, but no endemic equilibrium exists for R 0 < 1); 4 means forward bifurcation with multiple super-and subcritical endemic states (i.e., multiple equilibria exist for R 0 < 1 and for R 0 > 1).
The interpretation of the appearance of backward bifurcation and hysteresis phenomena is that in some season the mosquitoes are feeding and breeding and therefore their population size is growing.During such a season (like summer) they live long enough to cause malaria.This explains the sudden increase in infection levels.However, in other seasons where the female mosquitoes don't live enough to cause malaria, we get a sudden decrease in the level of infection.
4.2.Infectious Recovered Individuals (q 2 = 0, q 3 = 0) Now, we assess the impact of changing the relative infectiousness parameters of primary and secondary recovered individuals on the endemic behavior of Equation (5).Therefore, we consider the general case where none of the model parameters are neglected.More precisely, we assume that recovered individuals are capable of transmitting the infection with relative infectiousness q 2 , q 3 > 0.

Impact of the Relative Infectiousness q 3
To figure out the role of the relative infectiousness of secondary recovered individuals (with respect to primary infected human individuals) q 3 , we assume that q 2 as well as the other model parameters are constant, except q 3 , which is assumed to be increased, see Figure 10.The figure shows how the shape of the bifurcation curve can change as the value of q 3 increases.We notice that for small values of q 3 , the model exhibits forward bifurcation with a unique endemic equilibrium, Figure 10a.On increasing the value of q 3 , the model shows the existence of hysteresis, where there are two turning points, Figure 10b,c.In part (b) of the figure, there is a forward bifurcation with multiple supercritical, but not subcritical, endemic equilibria.However, with the increase in q 3 we, additionally, get multiple subcritical endemic equilibria, Figure 10c.If we allow q 3 to further increase, the model exhibits backward bifurcation, where there exist two subcritical endemic equilibria but only one supercritical endemic equilibrium, Figure 10d.It is clear that the increase in q 3 value implies a decrease in the critical basic reproduction level below which the infection does not persist, and therefore, it implies an increase in the minimum elimination effort of the malaria infection [39].

Impact of the Relative Infectiousness q 2
Similar to above, if we assume the relative infectiousness (of primary recovered with respect to primary infected human individuals) q 2 changes, and we let all other parameters be fixed (see the legend of Figure 11 to get to know the parameter values), then the bifurcation curve progression has a different behavior compared with the case of q 3 .For small values of q 2 , the model shows the existence of backward bifurcation, Figure 11a.On letting q 2 increase, the model exhibits forward bifurcation with multiple supercritical with (and without) subcritical endemic equilibria.A further increase in q 2 values implies the existence of forward bifurcation, where no subcritical endemic equilibria exist but exactly one endemic equilibrium does exist.This implies that a reduction in q 2 values implies an increase in the minimum elimination effort of malaria. .Bifurcation diagram for q 2 = 0, q 3 = 0.It shows that for q 2 = 0.5, and keeping all other parameters fixed, except q 3 , which is allowed to increase, the model develops forward (a), multiple supercritical (b), multiple subcritical (c), and backward (d) bifurcation.Simulations have been done for the parameter values given in Table 3, except for q 1 = 1.6, r 1 = 0.3 and r 2 = 0.7.Bifurcation diagram for q 2 = 0, q 3 = 0.It shows that for q 3 = 0.7, and keeping all other parameters fixed, except q 2 , which is allowed to increase, the model develops backward (a), multiple subcritical (b), multiple supercritical (c), and forward (d) bifurcation.Simulations have been performed for the parameter values given in Table 3, except for q 1 = 1.6, r 1 = 0.3 and r 2 = 0.7.

Role of Differential Mortality
In this section, we figure out the dynamics of the model if the infection-induced mortality rates δ 1 and δ 2 are neglected.In this case, the equilibrium human population size is simply Λ h /µ h .Hence, F 1 (λ h ) = F 2 (λ h ) and thus the left-hand size of Equation ( 23) would be reduced to a polynomial of degree three.If we further assume that recovered individuals are not capable of transmitting malaria, then the basic reproduction number reduces to while the condition (26) for backward bifurcation reduces to On the other hand, the Equation ( 23) from which we determine the equilibrium human force of infection λ h reduces to where It is clear that the coefficient C 3 ≥ 0, while the other three coefficients, namely C 0 , C 1 and C 2 , can be positive or negative.Thus, Equation (34) could have up to three feasible solutions depending on the model parameter values, see Figure 12 .This could be verified by studying the number of possible turning points in the bifurcation curve in the plane (R 0 , λ h ).To this end, Equation ( 34) is rearranged in a way such that the parameter Λ v (considered as a bifurcation parameter) is written as a function of the equilibrium force of infection λ h , say Λ v = H 1 (λ h ).The number of feasible solutions for the equation dH 1 /dλ h = 0 determines the number of turning points.Thus, the equilibrium force of infection at a turning point satisfies the equation 0.92 0.96 1.00 1.04 10 The figure is depicted for parameter values q 1 = 1.6 and keeping all other parameters fixed as in Table 3, except r 1 , which is allowed to take values 0.9 and 0.6, as shown in the head of the sub-figures.The model develops backward bifurcation, (a).Moreover, it shows forward bifurcation with the existence of multiple supercritical and subcritical endemic equilibria, (b).
It is clear that the coefficients of the two terms λ 4 h and λ 3 h are non-negative, while the other coefficients can be positive or negative.Thus, the number of sign changes in the coefficients of Equation ( 35) determines the number of feasible solutions.On the other hand, it is easy to check that at R 0 = 1, the bifurcation is backwards if and only if the term-free of λ h in Equation ( 35) is negative.This means C 11 < C 01 C 12 /C 02 .Therefore, the coefficient of λ 2 h in Equation ( 35) is If we assume the coefficient of λ h in Equation ( 35) to be negative, then (irrespective of the sign of the λ 2 h -coefficient) the number of sign changes will be one, and therefore there is a unique feasible solution of Equation (35).However, if we assume that the λ h -coefficient is positive, then, according to Equation (36), the λ 2 h -coefficient will be positive, and thus a unique solution for Equation (35) exists.Hence, if the bifurcation is backward, then there is a unique feasible turning point in the bifurcation curve drawn in the (R 0 , λ h ) plane.
Similarly, the case of forward bifurcation is argued.In this case, the term-free of λ h in Equation ( 35) is assumed to be positive.Thus, two feasible solutions for Equation (35) may exist if either or both coefficients of the λ 2 h and λ h are negative, while otherwise there is no solution.We summarize our result in the following proposition.Proposition 3. If the malaria-induced mortality rates as well as the infectivity of recovered individuals are neglected, the model exhibits backward bifurcation if and only if the condition (33) holds.Moreover, it shows the existence of multiple sub-and supercritical (hysteresis) endemic equilibria, with two turning points at the most.
Simulations have been performed for this case and these show that the model can show backward bifurcation and hysteresis (the existence of supercritical and subcritical endemic states) phenomena.

Summary and Conclusions
A mathematical model for the transmission dynamics of malaria is proposed.The model is of type SI in the vector and SIRS in the host.For the host population, it is differentiated between individuals who have never experienced the infection (called naive), individuals who experienced the infection only once (called primary), and those who experienced it at least twice (called secondary).It is assumed that primary and secondary infected individuals have different infectiousness, where the relative infectivity of secondary (with respect to primary) infected individuals is q 1 .Also, the susceptibility of primary and secondary susceptible individuals is different from naive susceptible individuals (those who have never experienced the infection).The relative susceptibility of primary and secondary (with respect to naive) susceptible individuals are r 1 and r 2 , respectively.Moreover, it is assumed that recovered individuals have some immunity to malaria.They do not become clinically ill, but they do harbor low levels of parasites in their blood streams, and therefore they can pass the infection to mosquitoes.It is further assumed that the relative infectiousness of primary and secondary recovered (with respect to primary infected) individuals are q 2 and q 3 , respectively.
The analysis shows that the model has a malaria-free equilibrium, which is proven to be locally asymptotically stable for R 0 < 1, while it is unstable if R 0 > 1.Here, R 0 is the basic reproduction number and is given by Equation (7).Center manifold analysis shows that the model exhibits backward bifurcation if the expression in Equation ( 25) is positive.Local sensitivity analysis for this expression shows that the most sensitive parameters are q 2 , r 1 and q 3 .Therefore, the backward bifurcation condition is expressed in the form of condition (26).This backward bifurcation condition (26) does depend on the parameters r 1 , q 1 , q 2 , and q 3 , but it does not depend on the relative susceptibility r 2 .
The equilibrium human force of infection is shown to be determined from the roots of a polynomial of degree six, which supports the possible existence of more than two endemic equilibria for certain model parameter values.However, in the special case δ 1 = δ 2 = 0 and q 2 = q 3 = 0, the analytical results show that the bifurcation diagram in the plane (R 0 , λ h ) can have up to two turning points, which means that the model shows the existence of forward hysteresis (i.e., the existence of forward bifurcation with multiple super-and subcritical endemic steady states).On the other hand, in the general case, numerical simulations have been employed to investigate the impact of changing the relative susceptibilities r 1 , r 2 , the relative infectiousness q 1 , q 2 , q 3 , and the differential mortalities δ 1 , δ 2 on the development of the bifurcation diagram and therefore on the persistence of malaria.
The implication of the appearance of multiple super and subcritical endemic states is that the basic reproduction number R 0 is no longer an indicator in the malaria control process.This has an effect, especially on the value of the critical threshold value of R 0 , below which malaria infection does not persist, especially in the case of multiple subcritical endemic states.In this case, malaria does persist, even for values of R 0 < 1.Consequently, the minimum effort required to control (eliminate) malaria is increased.Throughout the analysis, we came up with the following concluding remarks:

•
The higher the value of either or both of the relative susceptibilities r 1 (of primary susceptible) and infectiousness q 1 (of secondary infected) is, the higher the effort needed to eliminate malaria is.

•
The higher the value of the relative infectiousness q 3 of secondary recovered individuals is, the higher the effort required to eliminate malaria is.

•
The higher the value of the relative infectiousness q 2 of primary recovered individuals is, the lower the effort required to eliminate malaria is.
It is worth mentioning that our model could be extended to include factors representing the control of malaria.Among the strategies that could be considered in a future extended model are the use of insecticide-treated bed nets by the human hosts and studying the optimization of who should use these treated bed nets.

Limitation of the Work and Probable Future Work
It is worth mentioning that our model neglects various factors that significantly affect malaria dynamical spread.Among these factors are environmental factors, like temperature, rainfall, and humidity.These factors impact mosquito and parasite vital rates, and thus affect the transmission intensity, seasonality, and geographical distribution of malaria.Also, the latency period of malaria has been neglected and may cause inaccurate estimations in the effort needed to control malaria.Other factors include biological and physiological differences between human males and females.Some or all of these factors will be considered in a forthcoming work.

•
The second scenario is λh = 0.In this case, the rearrangement of the terms results in the following quadratic equation in the equilibrium force of infection λh where Once we obtain a feasible solution of (A12), we substitute in (A4)-(A8) to obtain the corresponding equilibrium.At λh = 0, we may derive the basic reproduction number R0 for Equation (1) as It is noteworthy that Equation (A12) is quadratic in λh and the number of feasible solutions could be determined by employing the implicit function theorem by following the approach shown in [40].Moreover, the left hand side of Equation (A12) could be seen as a bifurcation equation in the equilibrium force of infection λh and the bifurcation parameter β hv .Moreover, based on the use of the implicit function theorem, we may have Hence Thus, Equation (1) exhibits backward bifurcation if and only if The backward bifurcation condition (A15) implies that two malria-endemic equilibria do exist for values of R0 < 1.However, in the simple case δ h = 0, the condition (A15) does not hold and, consequently, the model has a unique malaria-endemic equilibrium if and only if R0 > 1.
The stability analysis of the malaria-free equilibrium for Equation (1) could be established by following the same approach shown in Section 2.3.However, the stability analysis of the malaria endemic equilibrium of Equation ( 1) for values of R0 > 1 could be established by following the same approach shown in [33].(5) To compute the direction of bifurcation for Equation ( 5) at R 0 = 1 , we make use of theorem 4.1 of Castillo-Chavez and Song [36].It can be checked that the Jacobian J(E 0 ) evaluated at β hv = β 0 hv (i.e., R 0 = 1) has: • a simple eigenvalue 0, while all other eigenvalues have negative real parts.Therefore, the center manifold theory can be used to analyze the dynamics of Equation ( 5), • a left eigenvector v = (v 1 , v 2 , ...., v 9 ) and a right eigenvector w = (w 1 , w 2 , ...., w 9 ) , where

Figure 1 .
Figure 1.A schematic diagram for the interaction and transition between the various states of Equation (5).

Figure 2 .
Figure2.The area for backward as well as forward bifurcation.Part (a) shows these regions in the plane (q 2 , r 1 ) for several values of q 3 .It shows that the backward bifurcation region extends with the increase in the q 3 value.However, part (b) shows that, in the plane (q 3 , r 1 ), the backward bifurcation region shrinks with the increase in the q 2 value.Simulations have been done for the parameter values given in Table3.

Figure 3 .
Figure3.Bifurcation diagram for q 2 = q 3 = 0.It shows that with the increase in q 1 , while fixing r 1 , r 2 and all other parameters, the model develops forward (a), multiple supercritical (b), multiple subcritical (c), and backward (d) bifurcations.Simulations have been done for the parameter values given in Table3, except for r 1 = 0.0194 and q 2 , q 3 as stated above.
Figure 4. Bifurcation diagram in the plane (q 1 , r 1 ), for q 2 = q 3 = 0 and with parameter values as in Table3.High enough values of q 1 and r 1 exhibits the existence of backward bifurcation.

Figure 5 .
Figure5.Bifurcation diagram for q 2 = q 3 = 0.It shows that with the increase in q 1 , while fixing r 1 , r 2 , and all other parameters, the model develops forward (a) to backward (b) bifurcation.Simulations have been done for the parameter values given in Table3, except for r 1 = 4.0194 and q 2 , q 3 as stated above.

Figure 6
Figure 6.Bifurcation diagram for q 2 = q 3 = 0.It shows that for small q 1 = 0.5, while fixing r 2 and all other parameters, but increasing r 1 , the model develops forward (a-c) to backward (d) bifurcation.Simulations have been done for the parameter values given in Table3, except for r 2 = 0.7 and q 2 , q 3 as stated above.

Figure 7 .
Figure7.Bifurcation diagram for q 2 = q 3 = 0.It shows that for q 1 = 1.6 and keeping r 2 and all other parameters fixed, except r 1 , which is allowed to increase, the model develops forward (a), multiple supercritical (b), multiple subcritical (c), and backward (d) bifurcation.Simulations have been done for the parameter values given in Table3, except for r 2 = 0.7 and q 2 , q 3 as stated above.

Figure 8
Figure 8. Bifurcation diagram for q 2 = q 3 = 0 and for fixed q 1 = 1.6.It shows that for a value of r 1 = 0.0194 < 0.6333 = r c 1 , increasing the relative susceptibility r 2 , while fixing all other parameters, the model develops forward (a), multiple supercritical (b), and multiple subcritical (c,d) bifurcation, but it cannot exhibit backward bifurcation.However, for r 1 = 0.7 > 0.6333 = r c 1 , the model exhibits backward bifurcation where the critical value of R 0 at which the turning point occurs decreases, (e-h).This induces an increase in the elimination effort.Simulations have been done for the parameter values given in Table3, except for r 1 = 0.0194 and q 2 , q 3 as stated above.

Figure 9 .
Figure 9.A schematic bifurcation diagram for the case q 2 = q 3 = 0 in the plane (q 1 , r 1 ).The notation: 1 means a forward bifurcation with unique endemic equilibrium for R 0 > 1; 2 means backward bifurcation; 3 means forward bifurcation with multiple supercritical endemic states only (i.e., multiple equilibria exist only for R 0 > 1, but no endemic equilibrium exists for R 0 < 1); 4 means forward bifurcation with multiple super-and subcritical endemic states (i.e., multiple equilibria exist for R 0 < 1 and for R 0 > 1).

Figure 10
Figure10.Bifurcation diagram for q 2 = 0, q 3 = 0.It shows that for q 2 = 0.5, and keeping all other parameters fixed, except q 3 , which is allowed to increase, the model develops forward (a), multiple supercritical (b), multiple subcritical (c), and backward (d) bifurcation.Simulations have been done for the parameter values given in Table3, except for q 1 = 1.6, r 1 = 0.3 and r 2 = 0.7.

Figure 11 .
Figure11.Bifurcation diagram for q 2 = 0, q 3 = 0.It shows that for q 3 = 0.7, and keeping all other parameters fixed, except q 2 , which is allowed to increase, the model develops backward (a), multiple subcritical (b), multiple supercritical (c), and forward (d) bifurcation.Simulations have been performed for the parameter values given in Table3, except for q 1 = 1.6, r 1 = 0.3 and r 2 = 0.7.

Figure 12 .
Figure12.Bifurcation diagram for δ 1 = δ 2 = q 2 = q 3 = 0.The figure is depicted for parameter values q 1 = 1.6 and keeping all other parameters fixed as in Table3, except r 1 , which is allowed to take values 0.9 and 0.6, as shown in the head of the sub-figures.The model develops backward bifurcation, (a).Moreover, it shows forward bifurcation with the existence of multiple supercritical and subcritical endemic equilibria, (b).

Table 1 .
Description of the state variables of Equation (5).
Number of susceptible humans who never experienced the infection at time t.I 1 (t) Number of primary infected humans (i.e., individuals who experience malaria for the first time) at time t.R 1 (t) Number of primary recovered (assumed carrier) humans from primary infection at time t, but are still slightly infectious.S 1 (t) Number of humans who recovered totally from the first infection and become susceptible again (after experiencing the infection only once) at time t.I 2 (t) Number of secondary infected humans (i.e., human individuals who acquired malaria infection at least for the second time) at time t.R 2 (t) Number of secondary recovered humans (i.e., recovered from secondary infection) at time t.S 2 (t) Number of secondary susceptible humans (i.e., susceptible individuals who have experienced the infection at least twice before) at time t.S v (t) Number of susceptible mosquitoes at time t.I v (t) Number of infectious mosquitoes at time t.N h (t) Total human population size at time t.N v (t) Total mosquito population size at time t.

Table 2 .
Description of the model parameters and their dimensions.
h Natural death rate of humans (i.e., deaths due to causes other than malaria) 1 Malaria-induced mortality rate for humans due to primary infections Time −1 δ 2 Malaria-induced mortality rate for humans due to secondary infections Time −1 ν 1 Rate of loss of immunity acquired by primary infection Time−1

Table 3 .
Table showing numerical values of parameters used in the numerical manipulations and the references (when applicable).