The Dynamics of Vector-Host Feeding Contact Rate with Saturation: A Case of Malaria in Western Kenya ()
1. Introduction
Malaria is an infectious disease caused by a parasite of the genus, Plasmodium. It is transmitted between human hosts by female anopheles mosquitoes as they seek blood meal for their eggs development. When a mosquito bites an infected person, a small amount of blood is taken in, which contains microscopic malaria parasites. When the mosquito takes its next blood meal, these parasites mix with the mosquito saliva and are injected into the person being bitten and the transmission process is perpetuated.
Malaria constitutes a big health problem especially within sub-saharan Africa and Asia. World Malaria Report 2011 estimates that it causes between 250 - 260 million infections and more than a million deaths (mostly among children in Africa), annually. In Kenya, reports show that despite the many control strategies to eliminate malaria, it has re-emerged and increased in incidence. The disease continues to wreck havoc on millions especially from the poor countries [1].
Vector abundance in Western Kenya is driven by temperature variation, ecosystem characteristics and human activities. The population varies depending on the site, the season and the species of the vector. Some sites in Western Kenya have 12.7 fold indoor resting densities during the long rainy season (March-June) and 23.3 fold during the dry season (January-March) [2]. This implies that the vector population is never constant as assumed in many models. On the other hand host population changes due to seasons and economic activities, natural deaths and death due to diseases like malaria and migration to urban centers and other regions for greener pastures. For our model to capture the reality of the epidemics in Western Kenya highlands, we assume that the host and mosquito populations change with time.
Most malaria models assume a constant human biting rate in their models, which means that hosts are freely available whenever a mosquito wants to bite, but in practice, this is more of a simplifying assumption. Research shows that for small host population, this rate is proportional to the host population size, and for large host population, it is constant [3,4]. The feeding cycle of a mosquito involves, host-seeking, feeding, resting, site-seeking, oviposition and host seeking resumes [5,6]. The probability of finding a host and successfully obtaining a blood meal depends on many factors. Among them is human avoidance and defensive behaviour [7,8].
If the mosquito survives this process, and has a successful blood meal, it rests, finds a larval habitat, oviposits and continues host seeking. Since this process drives malaria transmission, it is necessary to address the specific form of the mosquito-human contact process. Arditi [9], argues that the rates of (successful predation) contact between a predator and a prey is most properly a function of the ratios of their proportions. This would fit into malaria mosquitoes which inhabit homesteads and other areas where human hosts are available, like farms and urban areas [10]. It is also clear that this contact rate does not increase without bound, as the predator-prey ratio increases, this is because once a mosquito is fed, it rests before ovipositing, to resume host seeking and biting again [11].
When the predator-prey ratio value is low, the contact rate will be limited by the predators ability to find the prey, on the other hand if the ratio is high, the contact rate is limited by the predators satiation (desired predation rate) [4,9]. For malaria, the contact rate takes a similar course where the mosquito bites will increase as a function of host-vector ratio until the ratio reaches a critical level [12].
Saturation models are also not lacking in literature. A cholera model with saturation in the incidence was proposed by Capasso [13]. They argue that when there is a real threat to infected people becoming cautious and taking preventive measures which control further infection. Heesterbeek [14] formulated a saturated individual contact rate in relationships such as courting and marriage, where they assume that the population mixes randomly. Zu and Ma [15] analyzed a SEIR epidemic model whose latent period is described by delay and included a saturated incidence rate. Zhang and Ma [16] studied a SEIR model with saturation in contact rates and did a thorough analysis of its global dynamics.
In 2010, Ming and Li [17] formulated vector borne disease model, where they argue that increasing the density of the susceptible hosts with respect to infected ones leads to Holling type II saturation on the force of infection of host. In this model the biting rate of vectors and both populations are assumed to be constant. Further the model neglects disease related deaths a very crucial factor in malaria infection.
A model by [18] on dengue with variable human population is formulated and analyzed for both local and global dynamics. Ngwa et al. [19] analyzed the stability of a malaria model, with disease deaths, recovery and variable host and vector populations. However they assumed that the biting rate of vectors is constant hence their infection term is the one described in [20]. Realising then the need to predict the dynamics and transmission of malaria with great precision, we are motivated to engage in this study, as we pay particular attention to saturation in mosquito feeding habits and the varying host and vector populations.
The rest of the paper is subdivided as follows. Section 2 covers vector-host contact with saturation. In Section 3, the saturated contact process model is formulated. Section 4 is dedicated to the existence of equilibiria, while Section 5 studies the stability of the Disease Free and the Endemic equilibrium. Finally in Section 6 we give some results on numerical simulation.
2. Vector-Host Contact with Saturation
Vector-host contact results from the need for mosquitoes to obtain a blood meal for their eggs development. A given vector’s biting rate is limited by both host population density and its own feeding frequency [21].
Therefore the per vector biting rate should increase as a function of the host-vector ratio until the ratio reaches a critical threshold, which we denote
, above which, biting rate saturates and the average vector can feed at its preferred rate 
(contacts per vector per time). Below this threshold, the relative scarcity of hosts constrains the rate at which a vector can feed on the given type of hosts (it must seek other sources).
We assume that an average host can receive bites at a maximum rate 
beyond which it successfully defends itself against the vector (including leaving the place altogether) [12]. Then this threshold density ratio is given by
There are two ways of modeling the biting rate which increases for small host population and then approaches a maximum for large populations. The first is using a smooth verhulst-type function of the form
here 
is the preferred biting rate, and A, the population density ratio at which there is 
saturation, and it measures how soon the saturation occurs. Dietz [22] and Ming [17] have used an equivalent form of saturation [23].
The other way of modeling the vector biting rate saturation is by using a continuous and piecewise function with a “switch point”, 
deleanating the boundary between the two ranges as used by kribs in [23,24]. This function takes the form
(1)
where 
is the maximum preferred biting rate and 
is the threshold density ratio. There are other Hollingtype responses suggested in [14]. The difference between the two saturation models is the saturation sharpness, that is, how quickly the per vector biting rate levels off as the hosts become plentiful, that a vector can feed at its preferred rate,
. The saturation is gradual for the function 
model and sudden for the model described by 
see Figure 1.
The sharpness is crucial for our malaria model because mosquito populations can change drastically in a short span of time [23]. We also choose the predator-prey type or response since the biting in malaria transmission is vector initiated. The fact that the second model is also easy and captures more dynamics adds to the reason we apply it here in this malaria model.
From the reasons above, we shall assume a saturated contact process as used in [23,24] with the so-called Holling Type 1 form. Under this assumption, the per-vector contact rate can be described as a function of the hostvector density ratio 
as
When 
(many host per vector) the rate completely saturates at the maximum desired biting rate
, while for 
(i.e. few hosts per vector) 
, and the rate rises linearly with the hostvector ratio. The later is our interest in this study since the saturated contact process 
has been used in the classical Ross model for malaria [20]. Substituting 
in the function 
and multiply by
, we obtain the total biting rate which is
Figure 1. A graphical comparison between f1 and f2.
which can be rewritten in the form
We note that for the current host density the maximum number of vectors that can effectively bite hosts at one given time is
. Therefore the parameter 
is an important determinant in our model. It determines which of the two population densities is driving the biting contact rate.
To examine the rate of appearance of new malaria infections from the rate of mosquito feeding contacts, we have to take into account the probability of infection resulting from an effective contact where one party (host or vector) is infected with malaria parasite and the other is not. Let 
be the probability that such a contact between an infected vector and an uninfected host results in infecting the host, 
as the proportion of blood meal contacts between infected hosts and uninfected vectors which result in an infected vector. Further if 
are the susceptible and infectious hosts respectively and 
are the susceptible and infectious vectors respectively, the new infections will be given as defined in [12].
• For the Hosts
• For the Vectors
since 
In many malaria models, saturation has been assumed to be constant [20,25,26]. Here we consider the density dependent biting rate, and the populations ratio plays a vital role in the transmission. If the vectors-host ratio is low, the bites are few and hence the probability of transmission reduces too. As the ratio increases the infection will rise as a function of the ratio until it reaches the threshold and it becomes a constant. We wish to model this change in biting rate as the vector-host ratio changes and study its effect on the basic reproduction ratio.
3. The Model Equations
The model we derive here is mathematically equivalent to the classical Ross model [20]. The saturation in contact processes will address how the infection rates depend on this ratio. We assume that mosquito has a variable population growth such that birth
. For the human population, we assume a density dependent mortality rate, such that the total population vary with time and is modified by a logistic equation that include disease induced deaths.
A description of the variables and parameters used in the model follows in Tables 1 and 2 respectively.
The dynamics of our model will be governed by the following set of equations:
(2)
The term 
in the susceptible hosts compartment corresponds to a constant recruitment of susceptible hosts by natural birth. The transmission term 
corresponds to frequency dependent infection of susceptible hosts by infectious mosquitoes, on infection they move to the infectious compartment. The infected hosts who recover 
become susceptible again as malaria has no permanent immunity. The last terms 
represents per capita deaths of the susceptible, infected
Table 1. Variables used in the model related to infection contact process.
Table 2. Parameters used in the model related to infection contact process.
hosts respectively. In the susceptible mosquito vectors, 
represent the recruitment of susceptible mosquitoes by birth. The term 
corresponds to the transmission of malaria to an susceptible mosquito by and infected host. Both the susceptible and infectious mosquitoes are subject to natural deaths as defined in the terms 
respectively. Infective period of mosquitoes ends with their death due to their relatively short life-cycle so we do not have recovery or immune term in the vector equations [27,28].
All the parameters in the model are non negative and the model equations are well posed. For initial values 
in
, the solutions exist and remains in the region for all
.
In the absence of disease the host population dynamics is given by
. In this kind of demographic structure, the total human and mosquito population size
approaches a carrying capacity 
for any non zero initial population size.
The mosquito population 
also approaches a carrying capacity
. For ease of studying the system, we let setting 
and
, and the equation now takes the form
(3)
which is defined in feasible region (i.e. where the model makes biological sense)
where 
denotes the non-negative cone of 
including its lower dimensional faces. It is clear that 
is positively invariant with respect to (3). We denote the boundary and the interior of 
by 
and 
respectively.
4. Global Stability of the Disease Free Equilibrium
We recall the equations of the model and use the relation 
and
, to study the system
(4)
4.1. A Compact Positively Invariant Set
Using Barrier theorems (e.g. [29,30]) we prove that the following set
is a positively invariant compact set for system (4). Moreover 
is a global attractor on the nonnegative orthant
.
Since the ODE is Lipschitz, it is sufficient to check that the vector field induced by the system is either tangent or entering 
on the boundary
.
Clearly we have the following implications:
1) 
and
;
2)
;
3) Since 
we have
;
4) 
and
;
5) When 
and 
we have
;
6) When 
and 
we have
.
The preceding relations prove that all trajectories tends to
, which ends the proof of our claim. This also implies that all the trajectories are forward bounded.
We denote the demographic equilibria by 
and
.
4.2. Reduction of the Model
We remark that in the last equation only the state 
is appearing. Hence our system is a triangular system. Using Vidyasagar theorem on 
([31], Theorem A.4 given in annex) we can reduce the stability study to the stability of the equivalent system
(5)
This system is considered now on
A similar argument, as in the preceding section, shows that 
is a global attractor on the nonnegative orthant 
for system (5), as shown in Figure 2.
4.3. Basic Reproduction Ratio
For system (5) there exists a disease free equilibrium (DFE), which is 
in
. According to the technique by van den Driessche-Watmough [32], we define the transmission vector and the vector of other transfers defined on the “infected components”
The Jacobian of these vectors computed at the DFE gives
hence
The spectral radius of 
gives finally
4.4. Equilibria
Let 
denote an equilibria. We have the following system of equations
(6a)
(6b)
(6c)
The following relations are satisfied :
• From (6a)
(7)
• From (6c)
(8)
Replacing these values in (6a) gives
We see that if
, then 
and therefore
. In other words, we obtain the DFE. Then we suppose
. Then, from the last expression we obtain
Therefore
(9)
If 
we have 
and, from relation (7),
. From (6a) we obtain
Now from relation (8) we get 
and finally from (7), we get
.
In conclusion, if
, we have an unique endemic equilibrium in the interior of
.
4.5. Global Stability of the DFE
The DFE is on the boundary of
.
Theorem 4.1
The disease-free equilibrium 
of system (5) is globally asymptotically stable in the nonnegative orthant if and only if 
and is unstable if
.
Proof.
When 
the unstability of the DFE is a consequence of the theorem in [32,33].
To prove the global stability, we consider a LasalleLyapunov function 
[34-36] on the positively invariant compact set. By this definition we mean that on
, 
is continuous and nonnegative. The function 
is not positive definite. We define
We take advantage on the fact that the ODE (5 ) can be written
That we can write, in short, 
If we define
, then, the derivative along the trajectories 
is given by
But
, then 
For the inequality in the last row we assume 
on
.
Finally, with
, since
Now the largest invariant set contained in the set 
defined by
is certainly contained in the set of points for which 
or
. We have two situations:
• To be in an invariant set with 
has for consequences that 
and consequently
, hence to be at the DFE;
• To be in an invariant set with 
implies either
, thus to be at the DFE, or
, which is cannot be contained in an invariant set by the last equation system (5).
We have proved that the largest invariant set contained in 
and in the positively invariant set 
is the singleton constituted by the DFE. By Lasalle’s theorem [35,36] this proves that the DFE is globally asymptotically stable in
, hence in the nonnegative orthant. 
4.6. Stability of the Endemic Equilibrium
In this section we will prove that, when 
then the unique endemic equilibrium which is asymptotically stable.
Theorem 4.2
The endemic equilibrium 
of system (5), is given by Equations (7), (8) and (9), and is asymptotically stable in
.
Proof. We consider the Jacobian of system (5)
We can see that our system is not monotone and cannot be monotone for any orthant of
, since we have the 
entry of J negative and the 
entry positive.
The Jacobian computed at the EE, using the relation of (6), can be expressed as
We will examine the characteristic polynomial of
, denoted by 
.
It is well known that
The coefficient 
is equal to the sum of the principal minors
.
The Routh-Hurwitz criteria ensures that 
is Hurwitz iff 
for 
and
.
The upper 
-block of 
is
Hence the determinant is given by
The determinant of 
is
We consider now the trace of 
which is negative
We have now to compute 
the sum of the 
principal minors.
The first three requirements of Routh-Hurwitz’s criteria are satisfied. We have now to conduct a long calculation to prove 
Let compute 
This proves the asymptotic stability of the endemic equilibrium. 
4.7. Global Stability of the Endemic Equilibrium, when There Is No Disease Induced Mortality
We suppose in this section that
. In this case we can reduce the system to a two dimensional system, thanks to Vidyasagar’s theorem. The original system is then, for the stability point of view, equivalent to
(10)
In this case this system is a Ross system, and the global stability of the endemic equilibrium is well known, see for example [25,35,37].
4.8. Estimating 
The saturation threshold ratio 
is an important parameter that determines which of the two populations, that is the host and the mosquito, is driving malaria infection. We will assume that the indoor resting density for mosquitoes is 2.5 for each household that has approximately 6.5 persons (Githeko, personal communication), then
. We also assume that
, and 
, this gives
. Recalling our function for the feeding rate with saturation, 
implies that
, and the vector feeds at its preferred rate,
. When
, 
must be greater than
, so the biting rate will depend on the vector and not the host population. This would result to high vector host contact and high malaria transmission. In this case 
and malaria will persist in the population as shown in Figures 3(b) and (d). If however the host irritation from vector bites is high, then we assume a
, then
. This implies from our model that there are few host for each vector and the relative scarcity of hosts constrains the rate at which an average vector can feed. The biting rate therefore will be driven by the host and not the vector population. For its survival it may be forced to seek other sources for its blood meal. Reduced host biting translates to lower malaria infection and a reduced basic reproduction number. Then 
, the infection dies out, and the infectious population goes to zero as shown in Figures 3(a) and (c). Using the values in Table 3 we simulate system 2, and show the dynamics of various populations in the figures below.
Table 3. Parameters used in the model related to infection contact process.
5. Conclusion and Discussion
In this study, we developed a vector feeding habits saturation model for the spread of malaria with disease induced deaths and varying human and host populations. We have shown that the two populations drive the entire infection process through the threshold population density ratio
, which plays a vital role in the basic reproduction number. Our model captures the natural fluctuations known to occur in mosquito and host populations in malaria dynamics, and the effect of varying contacts between the vector and the host. The inclusion of diseaseinduced deaths was also of importance. Having in mind that majority of malaria deaths occur in children in Kenya, Mathematical analysis was done to establish that in the absence of the disease, a disease-free equilibrium will always exist if
. In the presence of the disease, that is when
, an endemic equilibrium is established with the infectious populations greater than zero.
We observe from Figure 3 that a decrease in 
increases 
and vice versa. When the human population is very low, mosquitoes will turn to other bloodmeal source, and malaria transmission goes down. Then 
controls the magnitude of malaria transmission. This implies that the best methods of controlling malaria in the highlands should target the adult mosquito, its biting habits and alternative sources of blood meals. Our results are consistent with results in literature that, 
is a threshold that completely determines the global dynamics of disease transmission [38].
6. Acknowledgements
We wish to acknowledge guidance from Calistus Ngonghala in the initial model formulation. Also Inria, France, the French Embassy in Nairobi and the University of Nairobi, Kenya, for their financial, logistic and moral, support during the writing of this article.