Modelling the effects of malaria infection on mosquito biting behaviour and attractiveness of humans

We develop and analyse a deterministic population-based ordinary differential equation of malaria transmission to consider the impact of three common assumptions of malaria models: (1) malaria infection does not change the attractiveness of humans to mosquitoes; (2) exposed mosquitoes (infected with malaria but not yet infectious to humans) have the same biting rate as susceptible mosquitoes; and (3) mosquitoes infectious to humans have the same biting rate as susceptible mosquitoes. We calculate the basic reproductive number, R0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0$$\end{document}, for this model and show the existence of a transcritical bifurcation at R0=1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0=1$$\end{document}, in common with most epidemiological models. We further show that for some sets of parameter values, this bifurcation can be backward (subcritical). We show with numerical simulations that increasing the relative attractiveness of infectious humans, increases R0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0$$\end{document} but reduces the equilibrium prevalence of infectious humans; decreasing the biting rate of exposed mosquitoes increases R0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0$$\end{document} and the equilibrium prevalence of infectious humans and mosquitoes; and increasing the biting rate of infectious mosquitoes has no impact on R0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0$$\end{document} or the equilibrium prevalence of infectious humans, but decreases the infectious prevalence of infectious mosquitoes. These analyses of a simple malaria model show that common assumptions around the relative attractiveness of infectious humans and the relative biting rates of exposed and infectious mosquitoes can have substantial and counter-intuitive effects on malaria transmission dynamics.


Introduction
Malaria is an infectious disease of humans, usually transmitted through the bites of Anopheline mosquitoes. Female mosquitoes bite humans (and other warm-blooded animals) for blood meals to provide protein for egg development. During each feed, the mosquito injects some saliva into the host to prevent clotting before sucking the blood into its stomach. After ingesting the blood, the mosquito rests for two to three days (depending on the ambient temperature) while it digests the blood. The mosquito then searches for a water body where it lays the eggs and then seeks another host to repeat its feeding cycle [23].
Mosquitoes experience different levels of mortality risk during the different phases of the feeding cycle. Host-seeking and attempting to feed on hosts is the most dangerous part of the cycle, where mosquitoes are most likely to die. The resting phase is usually the safest where mosquitoes are stationary. However, the mortality risks of the different stages can be altered by malaria control interventions. For example, insecticide-treated nets increase the mortality of host-seeking mosquitoes and indoor residual spraying with insecticides increases the mortality of resting mosquitoes [12].
Infectious mosquitoes contain malaria sporozoites in their salivary glands and infect humans when they inject saliva into them. In humans, the sporozoites pass through liver and asexual blood stages before transforming into infective sexual stages (called gametocytes). The incubation period in humans is relatively short (about 20 days) but the infectious period can last for multiple months. Humans infect mosquitoes when a mosquito sucks up a male and a female gametocyte during her blood meal. The gametocytes then fuse in the mosquito's stomach to eventually form an oocyst that releases sporozoites (that then travel to the mosquito's salivary glands to complete the malaria life cycle). The incubation period in the mosquito is about 10-12 days long, which is on the same order as the mosquito's life span so it has a substantial effect on malaria epidemiology and control [16].
Recent evidence has shown that malaria parasites can affect the mosquito's hostseeking behaviour to increase their probability of transmission [1,19,20]. Mosquitoes with oocysts (who are infected but not yet infectious) tend to be less mobile and spend more time resting. They therefore have longer feeding cycles and experience a lower mortality rate than uninfected mosquitoes. Mosquitoes with sporozoites tend to be more restless and are likely to take multiple blood meals in a single feeding cycle (that is, they bite multiple hosts before resting). Therefore, although they experience a higher mortality rate than uninfected mosquitoes, they also feed more frequently and thereby infect humans more often than they would if their feeding rate was unchanged.
Additionally, there is evidence that humans with a high density of gametocytes are more attractive to mosquitoes than humans without gametocytes [21], thereby making it easier for infectious humans to transmit to mosquitoes.
There is a long history of mathematical models of malaria starting with Ronald Ross' first ordinary differential equation (ODE) model for the proportion of infectious humans and mosquitoes [26]. Most models developed since then have been similar deterministic population-based models [2,3,10,24,25] although recently stochastic individual-based simulation models have increased in prominence [17,28]. However, these models have rarely considered the impact of malaria infection in mosquitoes on their feeding frequency or the impact of malaria infection in humans on their attractiveness to mosquitoes.
Here we develop a simple deterministic population-based model that extends existing models to include these effects. We first describe the model and perform a qualitative analysis to show the existence of disease-free and endemic equilibrium points, and the possibility of a backward transcritical bifurcation that is common to many such malaria models. We finally perform some numerical simulations that illustrate the effects of varying the infection-state dependent mosquito biting frequency and the relative attractiveness of infectious humans.

The model
We assume that the mosquito population is divided into three disjoint compartments: susceptible, exposed (mosquitoes carrying oocysts), and infectious (mosquitoes carrying malaria sporozoites). We denote the three populations at time t as S v (t), E v (t) and I v (t), respectively. In order to account for the effects of malaria parasites on the mosquito's host-seeking behaviour [1,19,20], we introduce two positive constants, ϕ E , and ϕ I , to represent the biting rate of exposed and infectious mosquitoes relative to susceptible mosquitoes, with ϕ E < 1 and ϕ I > 1.
Furthermore, by using the modelling approach employed first in [9] and later in [4][5][6]29], we divide humans into two compartments: susceptible, S h (t), who are not infected with malaria, and infectious, I h (t), who carry gametocytes. We denote the increased relative attractiveness of infectious humans as compared to susceptible humans (vector-bias parameter) by p. For simplicity, we assume that humans infected with malaria are always infectious; thus ignoring the incubation period and any variations in gametocyte density over time.
The balance equations lead to the following system of nonlinear ODES, We fix the human birth rate so that the stable human population size in the absence of disease is 1000 people. We fix the mosquito emergence rate so that the stable mosquito population size is 10,000 mosquitoes, leading to a density of 10 mosquitoes per person. Note that we only consider female mosquitoes in this model. All parameters are assumed to be positive except ϕ E , which is assumed to be non-negative where the upper dots denote the time derivative. The terms λ h (t) and λ v (t) denote the forces of infection on humans and vectors, respectively. The infection rate per susceptible human and per susceptible vector are given, respectively, by, The description and the baseline values of the parameters in (1) and (2) are given in Table 1. Denoting by N h the total human population and by N v the total vector population, i.e., whereR n + denotes the non-negative orthant of R n .

Local stability analysis of the disease-free equilibrium point
It is easy to check that model (1) admits the disease-free equilibrium (DFE) given by We use the next generation matrix method [13,14] to derive the basic reproduction number R 0 . The 'infected' compartments of model (1) are I h , E v and I v . The 'new infection' matrix F and the 'transition' matrix V defined in [14] are given, respectively, by The the basic reproduction number is the dominant eigenvalue of the next generation matrix F V −1 ,

Theorem 1
The disease-free equilibrium E 0 , given by (3), is locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1, where R 0 is given by (4).
This theorem follows from a straightforward application of Theorem 2 in [14] and showing that the model equations satisfy the five assumptions (A1)-(A5) in [14]. We note that the basic reproduction number may be written as The quantity K h is the number of mosquitoes infected by one human during his/her infectious life time. It is the product of the biting rate of a susceptible mosquito, the ratio of mosquitoes to humans at the disease-free equilibrium, the vector-bias parameter, the probability of disease transmission from infectious humans to susceptible mosquito per bite, and the average infectious period in humans. The quantity K v is the number of humans infected by one mosquito during her infectious life time. It is the product of the biting rate of infectious mosquitoes, the probability of disease transmission from infectious mosquito to susceptible human per bite, the average duration of the infectious period in mosquitoes, and the probability that a mosquito survives the exposed period. The basic reproduction number R 0 is the geometric mean of K h and K v since the infection from human to human goes through one generation of mosquitoes. We note that the increased frequency of biting of infectious, ϕ I , cancels out in K v so it does not influence R 0 .

Endemic equilibrium points
By endemic, we mean an equilibrium of system (1) where all components are positive. The components are solutions of the following system, represent any generic endemic equilibrium point of (1) and let be the forces of infection of humans and vectors at steady state, respectively. Solving equations (5) for the state variables provides where we have used the notation, Substituting (7) in (6) we get and By substituting (10) in (9), it can be shown that equilibria components satisfy the following equation, where the coefficients are given by, where R 0 is the basic reproduction number, given by (4). Equation (11) admits a trivial solution λ * h = 0 which corresponds to the disease-free equilibrium E 0 (3). Now, we assume λ * h = 0. The existence of endemic equilibria is regulated by the quadratic Note that the coefficient a 0 is always positive, and c 0 is positive (resp. negative) if R 0 is less than (resp. greater than) unity, respectively. The coefficient b 0 may be rearranged as: It follows that the number of endemic equilibria of model (1) depends on the coefficients a 0 , b 0 and c 0 as follows: (1) There is a unique endemic equilibrium if c 0 < 0 (i.e., R 0 > 1); (2) There is a unique endemic equilibrium if (4) There are no endemic equilibria otherwise.
The results of this section may be summarized in the following: A relevant aspect of this result is that conditions (15) indicate the occurrence of multiple endemic equilibria for R 0 < 1. Specifically, the disease-free equilibrium may co-exist with one or two endemic equilibria. From an epidemiological perspective, this means that the elimination of the disease in the population is no longer guaranteed by the classical threshold condition R 0 < 1. A new, smaller, threshold for R 0 must be determined. To this aim, we now express condition (15) in terms of the basic reproduction number R 0 as follows. Note that coefficient b 0 < 0 is equivalent to R 0 > R c , and c 0 > 0 is equivalent to R 0 < 1. Note also that b 2 0 − 4a 0 c 0 > 0 is equivalent to where, Equation (16) Choosing ψ < 0 ensures that μp − k 3 < 0, which implies > 0 and b 1 < 0. Setting it follows that condition (15) is equivalent to Remark 1 Condition (19) can be also expressed in terms of ω. Setting, and using the expression of R 0 , we obtain, where,

Bifurcation analysis
In this section we show that the sub-threshold occurrence of multiple endemic equilibria, as stated in Theorem 2, is the result of a transcritical backward (subcritical) bifurcation at R 0 = 1. This will also provide insight on the local stability properties of the endemic equilibria emerging from the bifurcation. We study the centre manifold near the criticality by using the approach developed in [7,14,15], based on general centre manifold theory [18]. In summary, this approach establishes that the normal form representing the dynamics of the system on the centre manifold is given byu = au 2 + bξ u, where, Here, the symbol ξ denotes a bifurcation parameter to be chosen, f i s denote the right hand side of system (1), x denotes the state vector, x 0 the disease-free equilibrium E 0 , D x denotes the differential operator with respect to x, D ξ denotes the differential operator with respect to ξ , and v and w denote the left and right eigenvectors, respectively, corresponding to the null eigenvalue of the Jacobian matrix of system (1) evaluated at x 0 for ξ = 0. We choose the mosquito biting rate, ω, as the bifurcation parameter. We observe that R 0 = 1 is equivalent to: so that the disease-free equilibrium E 0 is locally asymptotically stable when ω < ω * , and unstable when ω > ω * . Hence, ω = ω * is a bifurcation value. The Jacobian matrix of system (1), evaluated at E 0 for ω = ω * is given by where k 1 and k 2 are defined in (8). The eigenvalues are given by λ 1 = −μ, λ 2 = −η, and λ 3 , λ 4 , λ 5 eigenvalues of the submatrix, The characteristic equation ofJ is given by where and R 0 is given by (4). Since R 0 = 1 at the criticality, equation (26) becomes Thus, the Jacobian J * of the linearised system has a simple zero eigenvalue and the other eigenvalues have negative real part. Therefore the disease-free equilibrium E 0 is a nonhyperbolic equilibrium. In order to determine the coefficients (23) and (24), we look for the right and left eigenvectors corresponding to the zero eigenvalue. The components w i , for i = 1, . . . , 5, of the right eigenvector w are given by: Analogously, the components v i , for i = 1, . . . , 5, of the left eigenvector v are given by, Therefore for w 5 > 0, v 5 > 0, we have, Considering only the nonzero components of the left eigenvector, v, it follows that: where, where, The analysis performed in this sections may be summarized in the following theorem: Theorem 3 If A 0 < 0, then the malaria model (1) exhibits a backward bifurcation at R 0 = 1. If A 0 > 0, then the bifurcation is forward.

Remark 2
The change of stability of E 0 , from being stable to unstable, as R 0 crosses the critical value R 0 = 1, implies that b > 0 (see [14]). In our case, direct calculations shows that In the special case where there is no disease-induced death, i.e, α = 0, we have Since, by definition, p > 1 ⇒ 2 p − 1 > 0, it follows that in absence of disease induced death, the bifurcation is forward.
Finally, note also that when α = 0, the coefficients of equation (11) become, We see that a 0 is always positive. Additionally, since p > 1, it follows that b 0 and c 0 are always nonnegative whenever R 0 < 1. Therefore, there are no endemic equilibria bifurcating from the disease-free equilibrium when R 0 < 1.

Numerical simulations
We illustrate the backward bifurcation with numerical simulations in Figure 1 using the set of parameter values listed in Table 1 , are given by: E * = (20333, 1964, 99850354, 65808, 41919), which is locally asymptotically stable and E * * = (24079, 1650, 99884413, 50830, 32379), which is unstable. Figure 2 shows that increasing the relative attractiveness of infectious humans, p, increases R 0 . Figure 3 shows numerical simulations of the malaria model (1) for different values of p. The curve for p = 1 shows the result for the assumption used in most malaria models: that malaria infection makes no difference to the attractiveness  Table 1 except for ϕ I = 2, ϕ E = 1, and p varied as shown in the legend Fig. 4 One dimensional sensitivity analysis of R 0 to ϕ E with all other parameter values as specified in Table 1 except for ϕ I = 2 of humans. We see here that increasing p decreases the equilibrium prevalence of infectious humans, contrary to the effect of p on R 0 , but does not affect the equilibrium prevalence of infectious mosquitoes. Figure 4 shows that increasing the relative biting frequency of exposed mosquitoes, ϕ E , decreases R 0 . Figure 5 shows numerical simulations of (1) for different values of ϕ E . Consistent with the effect on R 0 , increasing ϕ E decreases the equilibrium prevalence of infectious humans and mosquitoes (with a stronger impact on the prevalence of mosquitoes).
The relative biting frequency of infectious mosquitoes does not appear in the expression for R 0 so changing this parameter has no effect on R 0 . Figure 6 shows numerical simulations of (1) for different values of ϕ I . Increasing ϕ I makes no difference to the equilibrium prevalence of infectious humans, but leads to a large decrease in the equilibrium prevalence of infectious mosquitoes. Figure 7 shows that the effects of simultaneously increasing p and ϕ E on R 0 are consistent with the one-dimensional sensitivity analyses shown in Figures 2 and 4. Additionally, as p increases, the dependence of R 0 on ϕ E increases.  (1) for parameter values as specified in Table 1 except for ϕ I = 1 and ϕ E varied as shown in the legend  Table 1 except for ϕ E = 1 and ϕ I varied as shown in the legend Fig. 7 Two dimensional sensitivity analysis of R 0 to p and ϕ E with all other parameter values as specified in Table 1 except for ϕ I = 2

Discussion and conclusions
We developed and analysed a simple malaria model to investigate the impact of relaxing three assumptions that mathematical models of malaria commonly make. We allowed the relative attractiveness of infectious humans (as compared to susceptible humans) to mosquitoes to increase; the feeding frequency of exposed mosquitoes (as compared to susceptible mosquitoes) to decrease; and the feeding frequency of infectious mosquitoes (as compared to susceptible mosquitoes) to increase.
We derived the basic reproductive number, R 0 , that provides a threshold condition for when the disease-free equilibrium loses stability. We showed that in common with many similar models of malaria, there is a transcritical bifurcation at R 0 = 1, that can be forward or backward depending on the parameter values, and is always forward if there is no disease-induced death rate. We provided threshold conditions for the direction of the bifurcation and the existence of zero, one or two endemic equilibrium points.
As may be expected, we showed that as the relative attractiveness of infectious humans, p, increases, R 0 increases. However, we also showed the surprising result that as p increases, the equilibrium proportion of infectious humans decreases, even as R 0 increases. This is because mosquitoes repeatedly bite infectious humans so they are less likely to pass the infection to susceptible humans, and therefore the same proportion of infectious mosquitoes leads to a lower proportion of infectious humans for higher values of p.
We showed that as the relative biting rate of exposed mosquitoes, ϕ E , increases, R 0 and the equilibrium proportion of infectious humans and mosquitoes decreases. This is reasonable because as the biting rate of exposed mosquitoes increases, their mortality increases so fewer exposed mosquitoes are likely to survive to become infectious.
We finally showed that varying the relative biting rate of infectious mosquitoes, ϕ I , has no impact on R 0 because the increased biting rate is cancelled by the shorter life span of the infectious mosquitoes. Correspondingly, varying ϕ I has no impact on the equilibrium proportion of infectious humans. However, increasing ϕ I leads to a substantial decrease in the equilibrium proportion of infectious mosquitoes because of their higher death rate. Importantly for malaria transmission, these fewer infectious mosquitoes are able to maintain transmission to humans at a higher level than if their biting rate was lower. These analyses of a simple malaria model show that common assumptions around the relative attractiveness of infectious humans and the relative biting rates of exposed and infectious mosquitoes can have substantial and counterintuitive effects on malaria transmission dynamics.