Periodic oscillations and backward bifurcation in a model for the dynamics of malaria transmission

Mathematical Biosciences 240 (2012) 45–62

Contents lists available at SciVerse ScienceDirect

Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs

Periodic oscillations and backward bifurcation in a model for the dynamics of malaria transmission Calistus N. Ngonghala a, Gideon A. Ngwa b, Miranda I. Teboh-Ewungkem c,⇑ a

National Institute for Mathematical and Biological Synthesis (NIMBioS), University of Tennessee, Knoxville, TN 37996, USA Department of Mathematics, University of Buea, Cameroon c Department of Mathematics, Lafayette College, Easton, PA 18042, USA b

a r t i c l e

i n f o

Article history: Received 8 January 2012 Received in revised form 23 May 2012 Accepted 10 June 2012 Available online 23 June 2012 Keywords: Vector demography and life style Malaria control Oscillations Hopf bifurcation Backward bifurcation

a b s t r a c t A deterministic ordinary differential equation model for the dynamics of malaria transmission that explicitly integrates the demography and life style of the malaria vector and its interaction with the human population is developed and analyzed. The model is different from standard malaria transmission models in that the vectors involved in disease transmission are those that are questing for human blood. Model results indicate the existence of nontrivial disease free and endemic steady states, which can be driven to instability via a Hopf bifurcation as a parameter is varied in parameter space. Our model therefore captures oscillations that are known to exist in the dynamics of malaria transmission without recourse to external seasonal forcing. Additionally, our model exhibits the phenomenon of backward bifurcation. Two threshold parameters that can be used for purposes of control are identified and studied, and possible reasons why it has been difficult to eradicate malaria are advanced. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction Malaria is one of the leading causes of morbidity and death of humans in many parts of the world and its impact has been felt for decades. Scientists recently determined that malaria played a role in the death of the boy Pharaoh Tutankhamen in about 1324 BC (See [1]). Despite its longevity, the malaria disease still constitutes a major health menace in many developing countries, especially within sub-Saharan African countries. According to the 2010 WHO report [2], an estimated 225 million cases, leading to about 781,000 malaria related deaths (90% in Africa), occurred in 2009. Most of the deaths were of children under five [3]. In addition to being the paramount source of morbidity and mortality in malaria endemic regions, malaria also weakens the active and potential workforce, thereby impacting economic growth negatively in the affected areas. Hence continuous research to find ways to effectively control the disease will save lives for future generations and also improve the economic conditions of the nations at risk.

⇑ Corresponding author. Tel.: +1 610 330 5328; fax: +1 610 330 5721. E-mail address: [email protected] (M.I. Teboh-Ewungkem). 0025-5564/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.mbs.2012.06.003

The disease, malaria, is caused by a micro-parasitic organism of the genius Plasmodium. Plasmodium falciparum is the most pernicious species of the malaria parasite and is transmitted indirectly1 from human to human by the female Anopheles sp. mosquito as it persistently quests for blood within the human population by bitting humans to obtain the blood that she needs for the development of her eggs. During each successful blood meal sortie, a female Anopheles mosquito may pick up the form of the parasite that is transmissible from humans to mosquitoes (called gametocytes) from an infectious human or deposit the form of the parasite that is transmissible from mosquitoes to humans (called sporozoites) to the human (who may be susceptible or already infected). Thus, a bite from an infectious female Anopheles mosquito can initiate the human phase of the parasite’s life cycle that will culminate with the production of gametocytes in the human’s blood system, or initiate the vector phase which will culminate with the release of sporozoites into the mosquito’s salivary glands. So, the Plasmodium parasite has an intricate life process, well adapted, with one part invested within the human host and the other part within the vector host

1 Unlike directly transmissible diseases such as HIV, communicable diseases of humans such as malaria are transmitted from one human to another by an intermediate transmission agent, the disease vector. This leads to a vicious cycle in which humans are infected by vectors and vectors are infected by humans in turn. However, except in the rare case of blood transfusion, humans do not infect other humans directly while vectors do not infect other vectors directly.

46

C.N. Ngonghala et al. / Mathematical Biosciences 240 (2012) 45–62

[4]. Therefore, any mathematical model for malaria must take this fact into consideration. Many mathematical models and paradigms have been proposed to study the dynamics of malaria transmission, starting with the pioneering work of Sir Ronald Ross [5], through the introduction of the notions of superinfection [6] and acquired immunity by continuing exposure and other variations of the central theme that the mosquito has a human biting habit [7–15]. Some recent models on malaria have instead focused on the malaria parasite as the basic unit of study, [16–20]. A lot has therefore been said about the different factors that can affect the dynamics and prevalence of malaria. These include factors that affect the developmental stages of the Plasmodium sp. parasite [17,18,21], conflicting arguments on the role of global warming, [22–26], availability of breeding sites for the mosquito [12,27–29]. Despite the large number of mathematical models geared towards understanding the biology, dynamics and transmission of malaria, little focus has been placed on models that combine disease dynamics together with the population dynamics of the malaria vector, and its interaction with the human population. The interaction between the mosquitoes and humans would invariably introduce a high variability in the mosquito population density and thus it is no longer realistic to trade information about infected mosquitoes for information about all mosquitoes and infected humans for the sake of mathematical tractability, as is often the case in most of the studies cited above. In fact, the most widely available control measure that is employed and that applies to all mosquito transmitted diseases of humans is the sanitation based control measure, which attacks the breeding habits of the mosquito and includes the use of insecticides both by spraying and on insecticide impregnated nets. These sanitation based control measures invariably point to the fact that over relatively reasonable time scales, the actor in the malaria transmission process that would show great variability in population densities would be the mosquito and not the human. We thus posit that research in malaria that integrates the disease transmission dynamics, the demographic properties of the vector and the different developmental stages of the parasite may provide novel insights towards better malaria control and eradication and consequently, assist in planning and implementing new intervention and vector control strategies. We therefore introduce here an approach to the development of models for malaria transmission wherein the mosquito vector is rightly placed at the center of the entire transmission process; an idea hitherto not fully exploited but partially ignored by simply viewing the mosquito population density as a parameter in models where the pseudo steady state approximation theory is evoked [30]. Our main objective, is to develop a mathematical model for the dynamics of malaria transmission which takes into consideration the population dynamics of the malaria vector and how these vectors interact with the human population. The Anopheles mosquito, undergoes complete metamorphosis and the average life span of the flying adult mosquito ranges from 2–3 weeks. Like any other disease vectors, the mosquito is capable of locating humans [21,31] and attempting as many times as possible to feed on the blood of the targeted humans. In the course of interaction with a human during a typical attempt, the mosquito may emerge victorious by acquiring the much needed blood meal, fail and go in for another attempt or may be killed. Nourishment is usually followed by resting and a choice of breeding site where the mosquito eventually lays its eggs. Hence, the feeding and reproductive habit of the malaria vector (taking a blood meal, resting, laying eggs and then going in for another blood meal, . . .) constitute a major factor to be considered when designing malaria transmission models and consequently, malaria intervention strategies. Here, we will consider these aspects and update the model originally developed by Ngwa [32] and further studied in [27,33] to

Table 1 Model variables and their definitions. Variable

Definition

Sh ðtÞ Ih ðtÞ N h ðtÞ Su ðtÞ Sv ðtÞ Sw ðtÞ Iu ðtÞ Iv ðtÞ Iw ðtÞ N mv ðtÞ

Susceptible humans at time t. Infected/infectious humans at time t. Total human population at time t. Susceptible vectors of class U at time t. Susceptible vectors of class V at time t. Susceptible vectors of class W at time t. Infected/infectious vectors of class U at time t. Infected/infectious vectors of class V at time t. Infected/infectious vectors of class W at time t. Total vector population at time t.

include the malaria disease stages. It is well known that the most general form of a compartmental model for the dynamics and transmission of malaria should be an SEIRS type model [12]. However, to illustrate our proposed new route to periodic oscillations in the dynamics of malaria transmission, we shall develop and analyze an SIS type model, which we believe will introduce an alternate framework for studying malaria control. The rest of the paper is organized as follows: In Section 2, we define the variables and parameters to be used in our model, derive the general SIS non-linear differential equation model and show that it is well posed with bounded solutions in the appropriate domain. In Section 3, we examine issues relating to the existence and stability of steady state solutions. In Section 4, we present results from computer simulations and then round up the paper with a discussion, biological implications and conclusions in Section 5.1.

2. The model parameters and derivation 2.1. Variables and parameters The mathematical model will take the form of a non-linear system of ordinary differential equations involving both the human and vector populations. The equations are derived based on the additional fact that, in the presence the malaria disease in the populations, both mosquitoes and humans can infect each other upon contact. While infected humans can recover from the malaria infection, it is assumed that once a mosquito is infected, it remains infected until death. To capture the life style of the mosquito vector in the model, the vector population is divided into three biologically realistic compartmental classes representing physiological status. These classes are: The class of fed and reproducing vectors returning from human habitats to vector breeding sites represented by the variable U; the class of unfed and resting vectors present at vector breeding sites represented by the variable V; and the class of unfed vectors questing (or foraging) for food (blood meal) in human habitats represented by the variable W. Note that type V vectors comprise previously fertilized female vectors at the breeding site that have just laid their eggs and all unfertilized, unfed and non-questing female vectors that are swarming at the breeding site (some of these are the newly emerging vectors). This is a new feature that is being introduced into malaria and other mosquitoborne disease modeling. The human population and the vector population in each compartmental class are again classified into two epidemiologically realistic compartmental classes2 representing disease status, namely: the susceptible class and the infected/ infectious class. Therefore, there are altogether eight dependent variables that continuously depend on time t. All variables and parameters in the human system will carry the subscript h, while those in 2 The basis for this approximating SIS (two classes) instead of the usual SEIRS (four classes) model is as explained above.

47

C.N. Ngonghala et al. / Mathematical Biosciences 240 (2012) 45–62

the vector subsystem will carry one of the subscripts v ; u or w depending on the physiological status of the vector. Thus, at any time t we have humans of type Sh ðtÞ and Ih ðtÞ, and vectors of type Sv ðtÞ; Sw ðtÞ; Su ðtÞ; Iv ðtÞ, Iw ðtÞ, and Iu ðtÞ as explained in Table 1. The total human and vector populations, N h ðtÞ and N mv ðtÞ are then given as

Nh ðtÞ ¼ Sh ðtÞ þ Ih ðtÞ and N mv ðtÞ ¼ Sv ðtÞ þ Iv ðtÞ þ Sw ðtÞ þ Iw ðtÞ þ Su ðtÞ þ Iu ðtÞ: Recruitment into and deaths from each of the above classes together with transitions between classes occur at rates given by the parameters as defined on Table 2. 2.2. Model derivation The model derivation uses a restricted form of homogeneous mixing based on the idea that the mosquito has a human biting habit. We assume that all new born humans and newly emerged mosquitoes are susceptible (no vertical transmission). A conceptual framework for the model is shown in Fig. 1. Flow from susceptible humans to infectious humans is as a result of contact between the susceptible humans, Sh , and an infectious questing vector, Iw . Recruitment of new susceptible humans is from births by both susceptible and infectious humans at rate kh . Deaths can occur amongst susceptible humans naturally at rate lh . In the presence of the malaria disease in the human population, infectious humans can either die naturally at rate lh or due to malaria infection at rate ch or recover and join the susceptible class at rate r h . The transmission rate of the parasite from mosquitoes to humans depends on the interaction between susceptible humans Sh and infectious questing mosquitoes Iw . When an infectious questing mosquito successfully takes a blood meal from a susceptible human with contact rate bhis , the human becomes infectious. Hence, the equations governing the human populations are

S_ h ¼ kh ðNh ÞN h þ r h Ih  bhis ðNh ÞSh Iw  lh ðNh ÞSh ; I_h ¼ b ðNh ÞSh Iw  ðl ðNh Þ þ r h þ c ÞIh his

h

h

ð2:1Þ ð2:2Þ

and the equation for the total human population is:

Fig. 1. The conceptual framework for the new SIS model for malaria transmission. Mosquitoes interact with humans and can transfer the malaria parasite from one human to another. A successful blood meal without being killed means a successful interaction and thus lead to mosquito population growth. Broken lines represent interactions between humans and questing mosquitoes. The different flow rates are explained in Table 2 and in the text.

  N_ h ¼ kh ðNh Þ  lh ðNh Þ Nh  ch Ih :

Next, we derive the equations that describe the population dynamics in each of the vector classes. Resting susceptible vectors are recruited when type Su or type Iu vectors reproduce at rate ~ kv per fed and reproducing vector. We account only for the vector populations involved in disease transmission, the female vectors, by setting kv ¼ g~ kv , where g 2 ð0; 1Þ is a constant. Reproduction by a mosquito is possible only when the mosquito visits a human being at the human habitat, takes a blood meal, rests, and then lays eggs at a chosen breeding site. Eggs hatch into larvae that further metamorphose into pupae and eventually develop into adult vectors of class Sv . Following the procedure in [32], it can be shown that the rate of emergence of new adult mosquito vectors at the breeding site is an expression of the form

av kv ðSu ðtÞÞSu ðtÞ |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl}

Table 2 System parameters representing the transition rates from one class to another together with recruitment and removal rates for both the human and mosquito populations. Parameter

Description

av

The rate at which successfully fed vectors return to breeding site. The rate at which vectors are attracted to human habitats. Natural human death rate. Natural mosquito death rate from each mosquito population class. Human birth rate. Mosquito birth rate.

av lh lu;v ;w kh ~v k rh bv bh p q p1 q1

g

Recovery rate amongst humans. Flow rate from susceptible questing mosquitoes to susceptible or infectious humans. Flow rate from infectious questing mosquitoes to susceptible or infectious humans. Probability that a type Sw vector successfully takes a blood meal from a susceptible human. Probability that a type Sw vector successfully takes a blood meal from an infectious human. Probability that a type Iw vector successfully takes a blood meal from an infectious human. Probability that a type Iw vector successfully takes a blood meal from a susceptible human. Proportion of emerging adult mosquitoes at the breeding site that are female.

ð2:3Þ

offspring from susceptible vectors

þ

av kv ðIu ðtÞÞIu ðtÞ |fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl}

:

ð2:4Þ

offspring from infected vectors

Additional recruitment of susceptible resting vectors is obtained via the susceptible fed vectors (Su ) that return at rate av from human habitats to the breeding sites to lay eggs. On the other hand, if the fed vectors that returned to the breeding site are infected with the malaria parasite (i.e. they fed on type Ih humans and thus are considered type Iu vectors), upon arrival at the breeding site and laying their eggs, they become vectors of type Iv , and thus serve as recruiters of type Iv vectors. After acquiring some rest, and after swarming, newly emerged adult female vectors (see expression (2.4)), and all previously fertilized adult susceptible (Su ) and infectious (Iu ) vectors can then return to human habitats at rate av ðN h Þ to seek blood meals. Additional exit out of each of the susceptible and infectious resting classes occur via natural death at rate lv . Taking into consideration the descriptions above and all assumptions made, the time rates of change of the Sv and Iv classes are, respectively,

S_ v ¼ av ½kv ðSu ÞSu þ kv ðIu ÞIu  þ av Su  ðlv þ av ðNh ÞÞSv I_v ¼ av Iu  ðl þ av ðNh ÞÞIv : v

ð2:5Þ ð2:6Þ

In Eq. (2.5), the birth rate function kv : ½0; 1 ! R, which is assumed to be continuously differentiable and strictly monotonic decreasing, can be modeled using any of the common birth

48

C.N. Ngonghala et al. / Mathematical Biosciences 240 (2012) 45–62

functions in the biological and ecological literature, such as, the Ricker function [34], and the Beverton–Holt function [35], etc. See, for example, [36,37] for more birth functions. The specific form of the birth rate function used in our paper will be described below. Additionally, in Eqs. (2.5) and (2.6), the term av ðN h Þ models the rate at which type V vectors from the breeding site (susceptible or infectious) are attracted to humans at human habitats. To derive an appropriate functional form for av ðN h Þ, let Bh and Ba be the blood preference indices per vector. That is, preference for either human or animal blood, respectively. Then Bh Sv and Ba Sv susceptible unfed vectors feed on humans and animals, respectively. The fractions of susceptible unfed vectors Sv that feed on humans and animals are, respectively, given by

on was susceptible or infectious, at the rates described above. However, when infectious questing vectors (Iw ) successfully feed on humans, they become infectious fed vectors of type Iu . Exits from the questing vector classes occur via natural deaths at constant rate lu and through vectors returning to the breeding site at rate av . Thus, the equations for the time rate of change of Su and Iu are, respectively,

Bh Sv Bh Nh ¼ Nh Bh Sv þ Ba Sv Bh þ Ba

N_ mv ¼ av kv ðSu ÞSu þ av kv ðIu ÞIu þ p1 bhis Sh Iw þ q1 bhii Ih Iw  ð1    pÞbv ss Sh Sw  ð1  qÞbv si Ih Sw  bhis Sh þ bhii Ih Iw    lv Sv þ lw Sw þ lu Su þ lv Iv þ lw Iw þ lu Iu ;

and

Ba Sv Ba Na ¼ Na ; Bh Sv þ Ba Sv Bh þ Ba

where N h is the total human population and Na is the total animal population. Therefore, susceptible unfed vectors at the breeding site are attracted to humans at rate

av ðNh Þ ¼ bNh =ðNh þ kÞ;

ð2:7Þ

where k ¼ Ba N a =Bh is a positive constant accounting for the existence of an alternative blood source for the vector, N h =ðN h þ kÞ is the proportion of susceptible unfed vectors that prefer human blood meals to animal blood meals and b is a positive constant to measure the flow rate of vectors from breeding sites to human habitats. We are assuming that the mosquito’s preference for human or other blood is not altered when the mosquito is infected. Once a susceptible (Sv ) or infectious (Iv ) vector from the breeding site arrives at a human habitat, it becomes a questing vector of type W (respectively, susceptible (Sw ) or infected (Iw )). A questing susceptible or infectious mosquito can either feed on a susceptible or infectious human. When a susceptible questing mosquito successfully takes a blood meal from an infectious human with contact rate bv si and success probability q, this mosquito becomes an infectious fed vector of type U, represented as (Iu ), else it fails with probability 1  q and is assumed killed. These contacts may lead to the production of new infectious vectors. If a questing susceptible vector instead successfully feeds on a susceptible human with contact rate bv ss and with success probability p, it becomes a susceptible fed vector of type U, represented as (Su Þ, else it fails with probability 1  p and is assumed killed. These contacts do not lead to any new infections. When an infectious questing mosquito successfully takes a blood meal from an infectious human with contact rate bhii and success probability p1 , this mosquito becomes an infectious fed vector of type U, (Iu ), else it fails with probability 1  p1 and is assumed killed. No new infections ensue from these contacts. If an infectious questing vector instead successfully feeds on a susceptible human, with contact rate bhis and probability q1 , it becomes an infectious fed vector of type U, (Iu Þ, else it fails with probability 1  q1 and is assumed killed. These contacts may lead to new infections in the human population but not in the vector population. Natural deaths in the Sw and Iw classes occur at rate lw which is considered to be constant at all times. From the above description, the respective governing equations for the population densities of the susceptible and infectious questing vectors assume the forms:

S_ w ¼ av ðNh ÞSv  ðlw þ bv ss Sh þ bv si Ih ÞSw ; I_w ¼ av ðNh ÞIv  ðl þ b Sh þ b Ih ÞIw : w

his

hii

ð2:8Þ

S_ u ¼ pbv ss Sh Sw  ðav þ lu ÞSu ; I_u ¼ p1 bhis Sh Iw þ q1 bhii Ih Iw þ qbv si Ih Sw  ðav þ lu ÞIu :

ð2:10Þ ð2:11Þ

Summing up Eqs. (2.5) and (2.6) and (2.8)–(2.11), we obtain the following equation that governs the total vector population:

ð2:12Þ

where av is given by (2.7). Now, in a general setting, it is natural to have the different vectors in different classes to have different physiological needs and properties and population parameters and study the general SIS model as derived above. In the present analysis it assumed that p1 ¼ q1 , and that deaths in each of the mosquito classes occur at a uniform constant rate lv . That is lw ¼ lu ¼ lv ¼ constant. We take into consideration the fact that the entire malaria transmission process is essentially driven by the human biting habit of the mosquito and that in classical epidemic models, contact rates can be assumed to be directly proportional to the population numbers as in [38] or constant as in [39] or some intermediate form as in [12,16,29,40] and assume mass action contact rates in all cases so that bv ss ¼ bv si ¼ bv ðN h Þ, bhis ¼ bhii ¼ bh ðN h Þ, which we shall simply refer to as bv and bh respectively. Next, we choose a specific functional form for the birth rate function kv . To this effect, we consider the logistic function,

  # ; kv ð#Þ ¼ k0 1  L

ð2:13Þ

where k0 is a positive growth constant, L > 0 is the environmental carrying capacity3 and # ¼ Su or # ¼ Iu . Notice that this is a positive real-valued continuously differentiable and strictly monotonic decreasing function of its argument provided # < L. The motivation for using the logistic birth rate function for the mosquito lies in the fact that it is linear, and in fact, is a general form of the linear approximation to any monotone decreasing continuously differentiable function. More nonlinear forms such as the Maynard-SmithSlatkin function, have been used to study the effect of nonlinear birth in the dynamics of mosquito populations [33]. Furthermore, we assume that human recruitment and deaths both occur at an equal constant rate lh (that is, kh ðNh Þ ¼ lh ðNh Þ ¼ lh , a constant) and that there are no disease-induced human deaths, that is, ch ¼ 0. These assumptions permit us to safely consider a constant human population (especially if there is always a human for a mosquito to attempt to feed on, in which case the total human population can be regarded essentially as constant (a buffer)). Combining Eqs. (2.1), (2.2), (2.5), (2.6) and (2.8)–(2.11) taking into consideration the above simplifying assumptions, the model we study then takes the form:

ð2:9Þ

When susceptible questing vectors of type Sw successfully feed on humans, they either become fed and reproducing vectors of type Su or of type Iu , depending on whether the human they fed

3 By environmental carrying capacity, we are referring to the constant maximum population size of a specific species that a designated habitat or breeding site can accommodate over an indefinite period of time without any deleterious repercussions.

49

C.N. Ngonghala et al. / Mathematical Biosciences 240 (2012) 45–62

S_ h ¼ lh Nh þ rh Ih  bh Sh Iw  lh Sh ;

Theorem 2.1. The closed set

I_h ¼ bh Sh Iw  ðlh þ rh ÞIh ;

W ¼ ðSu ; Sv ; Sw ; Iu ; Iv ; Iw Þ 2 R6þ : Nmv ¼ Su þ Sv þ Sw þ Iu þ Iv þ Iw



S_ u ¼ pbv Sh Sw  ðav þ lv ÞSu ; S_ v ¼ av kv ðSu ÞSu þ av kv ðIu ÞIu þ av Su  ðlv þ av ðNh ÞÞSv ; S_ w ¼ av ðNh ÞSv  ðlv þ bv Nh ÞSw ;

6 ð2:14Þ

I_u ¼ p1 bh Nh Iw þ qbv Ih Sw  ðav þ lv ÞIu ; I_v ¼ av Iu  ðl þ av ðNh ÞÞIv ;

2av k0 L

ð2:18Þ

lv

is positively-invariant and attracting with respect to the system (2.14). Proof. From Eq. (2.16) we see that

v

I_w ¼ av ðNh ÞIv  ðlv þ bh Nh ÞIw ;

N_ mv 6 av kv ðSu ÞSu þ av kv ðIu ÞIu  lv Nmv ;

together with the following equations for the total populations:

N_ h ¼ 0; _ v ¼ av kv ðSu ÞSu þ av kv ðIu ÞIu  l Nmv  ð1  pÞb Sh Sw Nm v v  ð1  qÞbv Ih Sw  ð1  p1 Þbh Nh Iw ;

ð2:15Þ ð2:16Þ

where av ðNh Þ; kv ðSu Þ and kv ðIu Þ are given by (2.7) and (2.13), respectively. To complete the formulation, system (2.14) must be solved with appropriate initial conditions. The simplest type of initial conditions would take the form

ðSh ð0Þ; Ih ð0Þ; Sv ð0Þ; Sw ð0Þ; Su ð0Þ; Iv ð0Þ; Iw ð0Þ; Iu ð0ÞÞ   ¼ S0h ; I0h ; S0v ; S0w ; S0u ; I0v ; I0w ; I0u :

since all human and mosquito population variables are non-negative and bv P 0; bh P 0 and 0 6 p; q; p1 6 1. From equation (2.13), if # 6 L then kv ð#Þ 6 k0 . Thus if Su 6 L, then Iu 6 L and

N_ mv 6 2av k0 L  lv Nmv ; ) Nmv ðtÞ   2av k0 L 2av k0 L lv t e þ Nmv ð0Þ  ; 6

lv

lv

Eq. (2.15) immediately shows that N h , the total human population, is constant and thus may appear as a parameter in the model. We note, in this case, that the mass action contact term is actually bh ðN h ÞSh Iw , where, to bring out the fact that it is really the proportion of susceptible humans that get infected, we can have bh ðN h Þ ¼ b=N h , and the constant b may be viewed as a mass action coefficient. Since N h is constant, we write bh ðN h Þ ¼ bh , a constant, and the above simplification becomes reasonable. The same applies to bv ¼ bv ðN h Þ. Now, for a constant human population, if we write Sh ¼ N h  Ih , we can compute Sh once we know Ih , which can be computed from the I_h equation in (2.14). 2.3. Well posedness, positivity and boundedness Since the system (2.14) models human and mosquito populations, it is assumed that all variables in the system are non-negative. This assumption yields the epidemiologically reasonable domain

D ¼ ðSh ; Ih ; Su ; Sv ; Sw ; Iu ; Iv ; Iw ;Þ 2 R8 : Sh P 0; Ih P 0; 0 6 Sh þ Ih 6 Nh ; Su P 0; Sv P 0; Sw P 0; Iu P 0; Iv P 0; Iw P 0; 0 6 Su þ S v þ S w þ I u þ I v þ I w 6 N m v Clearly, since the right-hand sides of (2.14), (2.15), (2.16) and their partial derivatives are continuous in D, it can be shown using standard techniques described in [41] that if initial conditions are specified within D for each of the state variables at time t ¼ 0 with Sv ð0Þ þ Sw ð0Þ þ Su ð0Þ þ Iv ð0Þ þ Iw ð0Þ þ Iu ð0Þ ¼ N mv ð0Þ, then there exist a unique solution satisfying these initial conditions with Sv ðtÞ þ Sw ðtÞ þ Su ðtÞ þ Iv ðtÞ þ Iw ðtÞ þ Iu ðtÞ ¼ N mv ðtÞ, for all t P 0. It can also be verified that if Nmv ð0Þ > 0, then N mv ðtÞ > 0 for all t, whereas if Nmv ð0Þ ¼ 0 then N mv ðtÞ ¼ 0 for all t. Of course similar arguments apply to the human equations with corresponding expressions. Thus, the system (2.14) is well posed from a mathematical and physical standpoint. Since the total human population is constant, we will focus on the total mosquito population (given by Eq. (2.16)) and present the following theorem:

8t P 0;

ð2:20Þ

where Nmv ð0Þ is the initial total mosquito population. Hence, if Nmv ð0Þ 6 2alv k0 L, Nmv ðtÞ 6 2alv k0 L implying that W is a positively-invariv

ð2:17Þ

ð2:19Þ

v

ant set. Notice that positivity is preserved if both Su P L and Iu P L since in this case N_ mv ðtÞ 6 0 for t P 0 and so the mosquito population is decreasing, in which case, Nmv ðtÞ 6 Nmv ð0Þ forall t P 0. Next, if N mv ðtÞ > 2av k0 L then from Eq. (2.20) we see that N_ mv ðtÞ < 0 and so lv

the mosquito population is again decreasing. From these analyses, solutions enter W in finite time and Nmv ðtÞ is bounded as t ! 1. Hence W is an attracting set. Therefore the solutions of model are considered epidemiologically and mathematically well posed in W. h 2.4. Nondimensionalization and reparametrization It is convenient to scale the model by introducing dimensionless variables. Starting with the fact that Sh ¼ N h  Ih , we introduce the new variables:

s ¼ ðav þ lv Þt; I ¼

Ih ; Nh

vs ¼

pbv Nh av ðNh Þ Sv ; Lðav þ lv Þðlv þ bv Nh Þ

ws ¼

pbv Nh Sw ; Lðav þ lv Þ

us ¼

Su ; L

wi ¼

av Lav ðNh Þ Iw ; ðlv þ av ðNh ÞÞðlv þ bh N h Þ

vi ¼

ð2:21Þ av L

lv þ av ðNh Þ

Iv ; ui ¼

Iu L

ð2:22Þ

and the dimensionless parameter groupings

b¼

av av ðNh Þbh L ; ðav þ lv Þðlv þ av ðNh ÞÞðlv þ bh Nh Þ

d¼

av av ðNh Þbh Nh p1 ; ðav þ lv Þðlv þ av ðNh ÞÞðlv þ bh Nh Þ

¼

lv þ bh Nh av av ðNh Þbh Nh p ; a¼ ; av þ lv ðav þ lv Þ2 ðlv þ bv Nh Þ

c¼

lv þ bv Nh l þ av ðNh Þ l þ rh ; q¼ v ; l¼ h ; a v þ lv av þ lv av þ lv

r ¼ q=p;

50

C.N. Ngonghala et al. / Mathematical Biosciences 240 (2012) 45–62

Q ¼ c þ q þ 1 > 0;

as the new parameters of the model. System (2.14) then takes the form:

I_ ¼ bð1  IÞwi  lI; u_ s ¼ ð1  IÞws  us ;

v_ s ¼ ako ðus ð1  us Þ þ ui ð1  ui ÞÞ þ aus  qv s ; _ s ¼ cðv s  ws Þ; w

ð2:23Þ

u_ i ¼ dwi þ rIws  ui ; v_ i ¼ qðui  v i Þ; _ i ¼ ðv i  wi Þ; w where the dots on the new variables represent differentiation with respect to s and the subscripts s and i, respectively, denote new variables in the susceptible and infectious vector compartments. It is clear from the scaling that 0 6 I < 1. Since 0 < p; p1 ; q < 1, it can easily be verified that 0 6 d < 1; 0 < a < q; b < L=Nh and ar < q. Assuming that b ¼ av (i.e. the flow rates to and from human habitats are the same) implies av ðNh Þ 6 av . Thus, 0 < q < 1. Hence, a feasible region for the parameter space is

X ¼ f0 6 b < L=Nh ; l > 0; 0 < a < q < 1; 0 6 d < 1;

 > 0;

r > 0; c > 0; k0 > 0g:

ð2:24Þ

Remark 2.2. Note that  and c can be bounded as follows: lv l bh N h and av þvl 6 c 6 1 þ blv Nh . av þl 6  6 1 þ l v

v

v

v

3. Existence and linear stability of steady states 3.1. The disease free system In the absence of the malaria disease, I ¼ 0; ui ¼ 0; wi ¼ 0, and system (2.23) reduces to

u_ s ¼ ws  us ; v_ s ¼ ak0 us ð1  us Þ þ aus  qv s ;

v i ¼ 0;

R ¼ c þ q þ cq > 0;

P ¼ cðq  aÞ > 0 since q > a:

ð3:4Þ

Proof. The steady state solutions of system (3.1) are obtained by setting the right hand side of that system to zero and solving the ensuing system of algebraic equations. It is easy to verify that the existence of realistic steady state solutions is determined by the parameter R , as given by (3.2) and that if R > 1, then E1 , as given in (3.3), exists. If R ¼ 1; E1 reduces to the trivial steady state E0 and when R < 1; E1 does not exist, as a realistic steady state. Stability of steady state solutions of system (3.1) is determined by the signs of the eigenvalues of the linearized system at the steady state, namely

0

1 0 10 1 1 0 1 us u_ s B _ C B CB C @ v s A ¼ @ aðk0 þ 1Þ  2ak0 us q 0 A @ v s A; _s 0 c c ws w |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Jðus ;v s ;ws Þ

where Jðus ; v s ; ws Þ is the Jacobian matrix of the system evaluated at the steady state. Thus if k is an eigenvalue of Jðus ; v s ; ws Þ, then we have

     Jðu ; v ; w Þ  kI ¼ 0 () k3 þ Q k2 þ Rk þ qc þ 2cak0 u s s s s  caðk0 þ 1Þ ¼ 0;

ð3:5Þ

where Q and R are as stated. It is then straightforward calculation to show that at E0 ; us ¼ 0 and all solutions of Eq. (3.5) have negative real parts when R 6 1 and there is a growing solution when R > 1, and that at E1 ; us ¼ 1  1=R , and solutions of Eq. (3.5) have negative real parts when QR  PðR  1Þ > 0, and that a Hopf Bifurcation occurs where PðR  1Þ ¼ QR. See [32] for details. h

ð3:1Þ

Remark 3.2. The condition QR > PðR  1Þ in Theorem 3.1 is equivalent to the condition 0 < k0 ðcÞ < ðcþqþ1ÞðcþqacþcqÞþcðqaÞ ¼ kH0 in the ðk0 ; cÞ space.

A similar system has been studied in [32]. Here, for completeness, we highlight the following important result, which we state and briefly proof in the form appropriate to the reparametrization given above:

Theorem 3.3. The trivial steady state for (3.1) is globally and asymptotically stable whenever R 6 1.

_ s ¼ cðv s  ws Þ: w

Theorem 3.1. In the absence of the malaria disease there is a threshold parameter R given by

R ¼ ak0 =ðq  aÞ

ð3:2Þ

with the following properties: (i) If R 6 1, there exist a unique steady state, the trivial steady state E0 ¼ ðus ; v s ; ws Þ ¼ ð0; 0; 0Þ which is linearly stable to small perturbations. (ii) If R > 1, there exists the trivial steady state E0 , and a non-trivial steady state E1 given by

E1 ¼ ðus ; v s ; ws Þ ¼

 1

 1 1 1 : ; 1  ; 1  R R R

ð3:3Þ

More so, for R > 1, we have that (a) The trivial steady state, which always exist, is linearly unstable to small perturbations. (b) The non-trivial steady state E1 is linearly stable to small perturbations whenever QR  PðR  1Þ > 0 and can be driven to instability via a Hopf Bifurcation at the point in parameter space where QR  PðR  1Þ ¼ 0, where

Proof. See Ngwa et. al. [33]. h Instability of E0 , for R > 1, is reasonable from a physical stand point and rules out any possibility of easily eradicating the mosquito population permanently whenever R > 1. The interesting feature here is that the non-trivial steady state, E1 , can lose stability at k0 ¼ kH0 and the instability is accompanied by the emergence of periodic solutions. The approximate period and amplitude of these oscillations can be computed as in [32,27]. We see here that the disease free model displays richer properties than is normally the case with the disease free system in previously studied mathematical models for malaria. 3.2. Analysis in the presence of the disease 3.2.1. A discussion on the basic reproduction number, R0 A threshold parameter that is of essential importance to infectious disease transmission is the basic reproduction number denoted by R0 , which measures the average number of secondary clinical cases of infection generated in an absolutely susceptible population by a single infectious individual through out the period within which the individual is infectious [42,43,12,16]. Generally, it is assumed that the disease will eventually disappear from the

51

C.N. Ngonghala et al. / Mathematical Biosciences 240 (2012) 45–62

community if R0 < 1 and may possibly establish itself within the community if R0 > 1. However, in some situations, there is the possibility of the occurrence of backward bifurcation which then complicates disease control even when R0 < 1. The critical case R0 ¼ 1 represents the situation in which the disease reproduces itself thereby leaving the community with a similar number of infection cases at any time. The definition of R0 specifically requires that initially everybody but the infectious individual in the population be susceptible. Thus, this definition breaks down within a population in which some of the individuals are already infected or immune to the disease under consideration. In such a case, the notion of reproduction number R becomes useful. Unlike R0 which is fixed, R may vary considerably with disease progression. However, R is bounded from above by R0 and it is computed at different points depending on the number of infected or immune cases in the population. One way of calculating R0 is to determine a threshold condition for which endemic steady state solutions to the system under study exists or for which the disease free steady state is unstable. Another method is the next generation approach where R0 is the spectral radius of the next generation matrix [43]. Using the next generation approach, we compute the basic reproduction number, e 0 as in [16,42,43] to be the spectral radius of the here denoted by R next generation Matrix M ¼ FV 1 where

0

0 B  B qbv Sw B F ¼B B 0 B @ 0 0 B B B V¼B B B @

0 0 bh Nh 0 0 0 0 0 0

lh þ r h

0

0

av þ lv

0

av

0

0

1

where R is the disease free threshold parameter defined in (3.2); and zero, one, or two endemic steady state solutions defined via

Ee ¼ ðI ; us ; v s ; ws ; ui ; v i ; wi Þ ¼ ðI ðwi Þ; us ; ws ðwi Þ; ws ðwi Þ; wi ; wi ; wi Þ;

whose existence and value is determined by the size of the threshold parameter R0 as well as the sign of the two parameter groupings A1 and A2 , respectively, defined as

A1 ¼ 1 

0

ð3:10Þ

that are coefficients of a uniquely determined nonlinear equation. Additionally, when Ee exists, it can be written in terms of wi , where wi , when strictly positive, is the scaled endemic steady state value of the infected questing mosquitoes. Proof. The steady states of the full malaria model are obtained by equating the right hand side of (2.23) to zero and solving the resulting algebraic equations. Thus if I ; us ; v s ; ws ; ui ; v i , and wi are the constant solutions of (2.23), then some algebra reveals that these can be written in terms of wi as

I ðwi Þ ¼

bwi ; bwi þ l bus

l

us ¼

wi ;

lð1  dÞ ; v s ðwi Þ ¼ ws ðwi Þ br

ui ¼ v i ¼ wi ;

 w2 i  A1 wi  A2 ¼ 0;

0 C p1 bh Nh C C C: lv þ av ðNh Þ 0 C av ðN h Þ lv þ bh Nh C A

ð3:11Þ

wi1;2 ¼

From these, we have

That is,

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u b b N S a a ðN Þ h w e0 ¼ u  v v h :q: t v h : R lv ðrh þ lh Þ av av ðNh Þ 1 þ ð1p1 Þbh Nh þ ðl þ b Nh Þðav þ av ðNh Þ þ l Þ

wi1;2

v

h

v

ð3:6Þ

which is positive since R > 1 and 0 6 d < 1. In this case, since e 0 6 1, implies 0 < R e 2 6 1 and R e 0 > 1 implies R e 2 > 1, we use 0
E0 ¼ ðI0 ; us0 ; v s0 ; ws0 ; ui0 ; v i0 ; wi0 Þ ¼ ð0; 0; 0; 0; 0; 0; 0Þ;

¼

A1 

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi A21 þ 4A2 2

ð3:13Þ

:

bðak0 r  qð1  dÞÞ 

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðbðak0 r  qð1  dÞÞ2 þ ð2ak0 lð1  dÞÞ2 ðR0  1Þ 2abk0 r

;

ð3:14Þ

e 0 takes the form In dimensionless parameters, R

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi brws brðR  1Þ e0 ¼ ¼ ; R lð1  dÞ lð1  dÞR

ð3:12Þ

with the coefficients A1 and A2 given by (3.10). We thus identify the coefficients of (3.12) with the parameter groupings A1 and A2 . Solving (3.12) for the positive value of wi , we obtain:

0

ð3:7Þ

a disease  free steady state solution,

 Edfe ¼ Idfe ; usdfe ; v sdfe ; wsdfe ; uidfe ; v idfe ; widfe   1 1 1 ¼ 0; 1   ; 1   ; 1   ; 0; 0; 0 ; R R R

  1 2 A2 ¼ us 1    us ¼ us ðR0  1Þ; R

where wi is the positive solution of the quadratic equation

1

lv

bqus ; alk0

¼ us þ

C C C C and C C 0 A 0 0

ð3:9Þ

ð3:8Þ

whose existence as a real and positive solution will be determined by the size, sign and different combinations of A1 ; A2 and R0 when referring to Eq. (3.13) or k0 and R0 when referring to Eq. (3.14), leading to a possibility of zero, one or two real solutions. h Remark 3.5. The endemic steady state postulated by Theorem 3.4 will exist only when wi > 0. Thus, if wi ¼ 0, that is, there is no disease in the population, then by (3.11), ui ; v i and I also vanish and we cannot have an endemic steady state. Moreover, us ¼ v s ¼ ws ¼ lð1  dÞ=br, from which we recover the trivial equilibrium solution if d ¼ 1. However, for 0 6 d < 1, since wi ¼ 0 leads to A2 ¼ 0 (or R0 ¼ 1, comparing with the non-trivial disease free equilibrium, (3.3), we have us ¼ lð1dÞ ¼ 1  R1 which leads to the br condition

lð1  dÞ q  a þ ¼ 1: ak0 br

ð3:15Þ

Note, however, that if A2 ¼ 0 wi can be different from zero since, by Eq. (3.12), if A2 ¼ 0; wi ¼ 0 or wi ¼ A1 . We can completely characterize the various realistic or biologically feasible endemic steady states in closed form by studying the sizes of R0 and k0 or the signs of A1 ; A2 and D, where

52

C.N. Ngonghala et al. / Mathematical Biosciences 240 (2012) 45–62

D ¼ A21 þ 4A2 ¼

ðbðak0 r  qð1  dÞÞ2 þ ð2ak0 lð1  dÞÞ2 ðR0  1Þ ðabk0 rÞ2

:

ð3:16Þ

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1  dÞ R0  1 lð1  dÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; us ¼ I ¼ ; v s ¼ ws br r þ ð1  dÞ R0  1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi   lð1  dÞ ð1  dÞð R0  1Þ ; ui ¼ v i ¼ wi ¼ 1þ br r pffiffiffiffiffiffiffiffiffiffiffiffiffiffi lð1  dÞð R0  1Þ ¼ : br 

We note the following possibilities: (a) A2 > 0 when R0 > 1; A2 < 0 when R0 < 1 and A2 ¼ 0 when R0 ¼ 1. Thus, the sign of A2 is dependent on the size of R0 ; A1 (b) Let kA01 ¼ qð1dÞ ar . Then A1 > 0 when k0 > k0 ; A1 < 0 when

k0 < kA01 and A1 ¼ 0 when k0 ¼ kA01 , where r ¼ q=p. Thus, the sign of A1 is dependent on the size of k0 . Various possibilities on conditions for the existence of endemic steady states and the number of steady states are illustrated in Fig. 2. Note that it is possible to have realistic solutions even when A1 < 0 or k0 < kA01 . Additionally, there are bifurcation points at which the the number of endemic steady states switch between zero, one or two steady states. The fact that for A2 < 0; i.e. R0 < 1 there are two non trivial endemic steady states, is a prerequisite for the occurrence of a backward bifurcation. We shall study the backward bifurcation in Subsection 3.2.3, but first, we analyze Eq. (3.14) further to obtain the following non-negative non trivial endemic steady states in closed form: 1. When A2 > 0 or R0 > 1 and A1 > 0 or k0 > unique endemic steady state is

since A1 is positive in (3.17) while here, A1 < 0 (i.e. k0 < kA01 ), the two steady state values will differ numerically. 3. When A2 > 0; A1 ¼ 0, and D > 0, the only endemic steady state is

kA01

and D > 0, the

4. When A2 ¼ 0; A1 > 0, and D > 0, the only endemic steady state is

I ¼

2. When A2 > 0; A1 < 0, and D > 0, the only endemic steady state has the same representation as in equation (3.17). However,

bðak0 r  qð1  dÞÞ ; bðak0 r  qð1  dÞÞ þ ak0 lr

us

lð1  dÞ  ; v s ¼ ws br ak0 lr þ bðak0 r  qð1  dÞÞ ¼ ð1  dÞ ; ui ¼ v i abk0 r2 ak0 r  qð1  dÞ ¼ wi ¼ : ak0 r ¼

ð3:19Þ

5. When A2 < 0; A1 > 0, and D ¼ 0, we have the single endemic steady state

I ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi bðak0 r  qð1  dÞÞ þ ðabk0 rÞ2 D lð1  dÞ  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; us ¼ I ¼ ; br 2 2ak0 lr þ bðak0 r  qð1  dÞÞ þ ðabk0 rÞ D qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 a k lr þ bð a k r  q ð1  dÞÞ þ ðabk0 rÞ2 D 0 0 v s ¼ ws ¼ ð1  dÞ ; 2 2abk0 r ð3:17Þ

ð3:18Þ

bðak0 r  qð1  dÞÞ ; 2ak0 lr þ bðak0 r  qð1  dÞÞ

us

lð1  dÞ  ; v s ¼ ws br 2ak0 lr þ bðak0 r  qð1  dÞÞ ¼ ð1  dÞ ; ui ¼ v i 2abk0 r2 bðak0 r  qð1  dÞÞ ¼ wi ¼ : ð3:20Þ 2abk0 r ¼

6. When A2 < 0; A1 > 0, and D > 0, we have the following two realistic endemic steady states:

Fig. 2. Regions in the A2 versus A1 parameter space where zero, one or two biologically feasible endemic steady states exist. Lines on which biologically feasible endemic steady states exist are solid, while lines on which no feasible endemic steady states exist are broken. Note that A2 > 0 when R0 > 1; A2 < 0 when R0 < 1 and A2 ¼ 0 when A1 A1 A1 R0 ¼ 1; A1 > 0 when k0 > qð1dÞ ar ¼ k0 ; A1 < 0 when k0 < k0 and A1 ¼ 0 when k0 ¼ k0 .

53

C.N. Ngonghala et al. / Mathematical Biosciences 240 (2012) 45–62

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 bðak0 r  qð1  dÞÞ  ðabk0 rÞ D lð1  dÞ  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; us1;2 ¼ I1;2 ¼ ; br 2 2ak0 lr þ bðak0 r  qð1  dÞÞ  ðabk0 rÞ D

v s1;2 ¼ ws1;2 ¼ ð1  dÞ ui1;2 ¼ v i1;2 ¼ wi1;2 ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2ak0 lr þ bðak0 r  qð1  dÞÞ  ðabk0 rÞ D

2abk0 r2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi bðak0 r  qð1  dÞÞ  ðabk0 rÞ2 D 2abk0 r

;

0

l B  B us B B 0 B B A¼B B 0 B  B rus B B @ 0

0

0

0

0

0

1

0

0

akc0 ð1  2us Þ þ a q 0 akc0 c c 0 0

0

:

0 1

0 0

0

0

0

1

0

0

0

0

q

q

0

0

0

0



1

b C 0 C C 0 C C C 0 C C C d C C C 0 A  ð3:22Þ

ð3:21Þ

and the characteristic equation of A is given by The full malaria model in the presence of the disease thus displays richer possibilities in terms of steady state existence that is not the case with the classical models for malaria that have been studied thus far. 3.2.3. Characterization of backward bifurcation Epidemiological models can exhibit some interesting dynamical behavior near R0 ¼ 1. In most cases, there is a supercritical or forward bifurcation at R0 ¼ 1. However, situations do arise in epidemiological modeling where we run into a double bifurcation at R0 ¼ 1. In this case, we have both a supercritical or forward bifurcation and a subcritical or backward bifurcation. The occurrence of a forward bifurcation is accompanied by the emergence of an asymptotically stable endemic equilibrium and loss of stability of the disease free equilibrium. A pre-requisite for the occurrence of a backward bifurcation is the existence of two endemic equilibria when R0 < 1. The other requirement is that one of the endemic equilibria together with the disease free equilibrium have to be stable while the second endemic equilibrium has to be unstable. The phenomenon of backward bifurcation was first established by Huang et al. [44] and further exploited by [45–48] to name just a few. The existence of a backward bifurcation tells us that simply reducing R0 below unity might not be enough to eradicate certain diseases. Here, we apply the approach described in [48] to explore the existence of a backward bifurcation in model (2.23) for R0 < 1.

3  k þ Q k2 þ Rk þ PðR  1Þ k4 þ ðQ 1 þ lÞk3 þ ðlQ 1 þ R1 Þk2 ð3:23Þ þðlR1 þ P1 Þk þ P1 lð1  R0 Þ ¼ 0; where P; Q and R are as defined in (3.4) and

Q 1 ¼  þ q þ 1 > 0;

Proof. The matrix of the linearized system about the disease free equilibrium Edfe defined in (3.8) is

n X

a¼

0.3

v k wi wj

k;i;j¼1

n X

¼

v k wi

I

0.0 0.0

when

R0 ¼ 1,

and

@ 2 fk ðx ; kc0 Þ @xi @xj

@ 2 fk ðx ; kc0 Þ ¼ av 3 w2 ð1  2us Þ þ av 3 w5 ; @xi @k0

w ¼ ðw1 ; w2 ; w3 ; w4 ; w5 ; w6 ; w7 ÞT  T b A1 bus A1 bus A1 ¼ ; ; þ ; þ  ; 1; 1; 1 : l us l us l us

Unstable Ee Stable Ee Stable Edfe Unstable Edfe

0.4

bb

0.1

R > 1

with the right eigenvector given by the solution of the system Aw ¼ 0 and defined as

0.8

0.2

Since

¼ 2bv 1 w1 w7  2v 2 w1 w4  2akc0 v 3 ðw22 þ w25 Þ þ 2rv 5 w1 w4 ; b

(b) 1.2

Stable Ee Unstable Ee Stable Edfe Unstable Edfe

qa . a½1lð1dÞ br 

ðQ 1 þ lÞðlQ 1 þ R1 Þ  ðlR1 þ P1 Þ ¼ lQ 1 ðl þ Q 1 Þ þ ð þ qÞð1 þ R1 Þ > 0, the Routh Hurwitz conditions assure us that all other eigenvalues of A have negative real parts when 0 < k0 ðcÞ < ððc þ q þ 1Þðc þ q þ cqÞþ cðq  aÞÞ=ðacÞ ¼ kH0 . Next, by the methods of Theorem 4.1 in [48], we have to determine a and b, as well as a nonnegative right eigenvector w ¼ ðw1 ; w2 ; w3 ; w4 ; w5 ; w6 ; w7 ÞT and a left eigenvector v ¼ ðv 1 ; v 2 ; v 3 ; v 4 ; v 5 ; v 6 ; v 7 ÞT that correspond to the simple eigenvalue of the matrix A. It can be shown that

ui = vi = wi

(a) 0.4

ð3:24Þ

Note that when R0 ¼ 1, there is a simple zero eigenvalue and aÞ kc0 ¼ a½bbrrðqlð1dÞ ¼

k;i¼1

Theorem 3.6. Let k0 ðcÞ < ððc þ q þ 1Þðc þ q þ cqÞ þ cðq  aÞÞ= ðacÞ ¼ kH 0 . Then model (2.23) undergoes a backward (subcritical) bifurcation at R0 ¼ 1.

R1 ¼  þ q þ q > 0;

P1 ¼ qð1  dÞ > 0; since 0 6 d < 1:

R0

bb

R0 (a) 0.5

R0

1.0

1.5

0.0 0.0

(b) 0.5

R0

1.0

1.5

Fig. 3. A plot depicting both a backward (subcritical) bifurcation and a forward (supercritical) bifurcation for L ¼ 1  104 ; N h ¼ 2  106 and the other parameters in Table 3 bb and Eq. (4.1). When R0 < Rbb 0 ¼ 0:45, there exists a unique stable disease free steady state. For R0 < R0 < 1, there exists a stable disease free steady state, a stable endemic bb steady state (solid line above the R0 axis between R0 ¼ Rbb 0 and R0 ¼ 1), and an unstable endemic steady state (broken line between R0 ¼ R0 and R0 ¼ 1). Finally, for R0 > 1 there exists a stable endemic steady state (solid line segment to the right of R0 ¼ 1) and an unstable disease free steady state (broken line segment to the right of R0 ¼ 1). For the parameter regime used here, the values of a and b are, respectively, 0:5175 and 0:0776.

54

C.N. Ngonghala et al. / Mathematical Biosciences 240 (2012) 45–62

v through the equation v :w ¼ 1 to obtain brw3 brw3 ; 2ð2brw23 þ akc0 lðw22 þ 1ÞÞ 2ð2brw23 þ akc0 lðw22 þ 1ÞÞ  brw23 þ akc0 lðw22 þ 1Þ ; 0; 0 ; 2brw23 þ akc0 lðw22 þ 1Þ

We then compute

v

 ¼ 0; 0;

Note that kc0 ¼

qa að1lð1dÞ br Þ

qa qa c c  ¼ að1u  Þ ) ak0 ¼ 1u and ak0 ð1  2us Þþ s

s

Proposition 3.8. The trivial Equilibrium E0 ¼ ð0; 0; 0; 0; 0; 0; 0Þ, which always exist, is linearly stable to small perturbations whenever R 6 1 and unstable when R > 1. Proof. Linear stability of E0 is determined by the eigenvalues of the Jacobian matrix, J evaluated at the trivial equilibrium, E0 ¼ ð0; 0; 0; 0; 0; 0; 0Þ. Here,

Moreover, 0 < us < 1 and 0 < 1  us < 1 since

0

when R0 ¼ 1; us ¼ 1  R1 < 1 and 1  us ¼ R1 < 1. Also note that since we are in the right half plane of the bifurcation diagram (Fig. 2), A1 > 0. With the vectors w and v , we can simplify a and b as

B B B B B B JðE0 Þ ¼ B B B B B @

a ¼ q  ðq  a

u Þ 1us  . s

A21 þ1 u2 s

a ¼ 2akc0

!

r brðbu2 s þ l A1 Þ v 3 þ 2 w3 v 5 ¼ > 0; b 2l2 us l

  A1 A1 ¼ aðð1  þ w5 Þv 3 ¼ a   2  us þ 1 v 3 us us    A1 a ¼ a   A1  A1 þ 1 v 3 ¼  A1 ð1  us Þ þ ð1  A1 Þus v 3 : us us 2us Þw2

Since 0 < us < 1 and 1  A1 > 0; b > 0. Hence, the model exhibits a backward bifurcation. h Numerical results confirming the occurrence of a backward bifurcation are presented in Fig. 3. Remark 3.7. Note that there exists a threshold value of R0 , namely, bb Rbb 0 such that for R0 < R0 , we have a unique stable disease free steady state and for Rbb < R0 < 1, a stable disease free steady state 0 coexists with two endemic steady states, one of which is stable. The value of Rbb 0 can be computed by setting D ¼ 0, where D is the discriminant of Eq. (3.12) and solving for R0 . This yields

Rbb 0 ¼ 1



A1 2us

2

¼1

 2 bðak0 r  qð1  dÞÞ ; 2ak0 lð1  dÞ

from which we deduce that

k0 < kbb 0

Rbb 0

l

0

0

0

0

0

1

0

1

0

0 0

aðk0 þ 1Þ q 0 ak0 0 c c 0

0

0

0

0

1

0

0

0

0

q

0

0

0

0

0

b

1

0 C C C 0 0 C C C 0 0 C C 0 d C C C q 0 A   0

If k is an eigenvalue of JðE0 Þ, then we easily establish that k is the solution of the equation

  ðk þ lÞ k3 þ Q k2 þ Rk þ Pð1  R Þ k3 þ Q 1 k2 þ R1 k þ P 1 ¼ 0; ð3:26Þ where P; Q and R are as defined in (3.4) and P 1 ; Q 1 and R1 are as defined in (3.24). Thus the characteristic polynomial in this case decouples into three polynomials. We observe that the coefficients of the second cubic polynomial on the left hand side of Eq. (3.26) together with Q 1 R1  P1 are all positive and hence this part of the equation can only contribute eigenvalues with negative real parts to the spectrum of JðE0 Þ. Thus eigenvalues with positive real parts, in view of growing solutions in the linear regime, can only be generated from the first cubic polynomial of (3.26). But this is precisely the polynomial (3.5) with us ¼ 0 and the result in the relevant part of Theorem 3.1 carries over. h

ð3:25Þ

< 1 when

bqð1  dÞ qð1  dÞ   ¼ ¼ aðbr  2lð1  dÞÞ ar 1  2 lð1dÞ br   lð1dÞ bq br  : ¼ al 1  2 lð1dÞ br

What this indicates is that when the disease is endemic in the community, a control strategy that focuses on reducing this value to a value below one may not be successful, unless R0 < Rbb 0 < 1. 3.2.4. Stability of Steady States To investigate the linear stability of the steady state solutions for system (2.23), we begin by linearizing the system about the steady state ðI ; us ; v s ; ws ; ui ; v i ; wi Þ. We obtain the following linearized system:

Proposition 3.9. The disease free equilibrium, which exists only when R > 1, is linearly stable to small perturbations whenever QR PðR  1Þ > 0 and R0 < 1, but can be driven to instability via a Hopf Bifurcation, even when R0 < 1, at the point in parameter space where QR  PðR  1Þ ¼ 0. In addition, the disease free equilibrium is always linearly unstable to small perturbations whenever R0 > 1. Proof. The disease-free equilibrium Edfe is given in (3.8) and the Jacobian matrix and the characteristic equation are given by (3.22) and (3.23), respectively. That is,

k3 þ Q k2 þ Rk þ PðR  1Þ ¼ 0;

ð3:27Þ

or

k4 þ ðQ 1 þ lÞk3 þ ðlQ 1 þ R1 Þk2 þ ðlR1 þ P1 Þk þ P 1 lð1  R0 Þ ¼ 0;

1 0 1 0 1 I ðbwi þ lÞ 0 0 0 0 0 bð1  I Þ I_ C B _ C B Bu C   ws 1 0 1I 0 0 0 C B us C B B sC C B C B B C   C B v_ s C B B vs C 0 a k ð1  2u Þ þ a  q 0 a k ð1  2u Þ 0 0 0 0 s i C B C B B C C       B C B C B _ C¼B 0 0 c c 0 0 0 CJðI ;us ;v s ;ws ;ui ;v i ;wi Þ B ws C: Bw C B sC B B C C B u_ C B B ui C rws 0 0 rI  1 0 d C B iC B B C C B C B B C q q 0 0 0 0 0 A @ v_ i A @ @ vi A _ wi 0 0 0 0 0    wi |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} 0

0

ð3:28Þ

55

C.N. Ngonghala et al. / Mathematical Biosciences 240 (2012) 45–62 Table 3 Some physically meaningful ranges of parameter values for our model. Parameter

Range of values

Value used

Reference

lh , (per day)

½1=ð72  365Þ; 1=ð45  365Þ ½1=120; 1=80 ½1=21; 1=14 0:5 6 av 6 1 0:5 6 av ðN h Þ < av Variable Variable Variable

1=ð60  365Þ 1=120 1=21 1.0 0.8 variable variable variable

[49] [12] [12,17] [12] [32]

rh , (per day) lv , (per day) av av ðNh Þ k0 L Nh

where P; Q and R are as defined in (3.4), and P1 ; Q 1 and R1 are as defined in (3.24). Now, when R0 > 1, lP1 ð1  R0 Þ < 0 and hence there is one sign change in the sequence of coefficients of Eq. (3.28). Thus there exist at least one positive real value of k and so the disease free equilibrium Edfe is linearly unstable to small perturbations whenever R0 > 1. When R0 6 1, all coefficients of Eq. (3.28) are positive and thus by Descartes rule of signs, there are no positive real solutions k. We use the Routh-Hurwitz criteria to ascertain that when R0 6 1 all the roots of the Eq. (3.28) have negative real parts. To do this, we first notice that Q 1 þ l > 0, lR1 þ P1 > 0 and lP1 ð1  R0 Þ > 0 when R0 < 1 and deduce that all we need to show is that

ðQ 1 þ lÞðlQ 1 þ R1 ÞðlR1 þ P1 Þ > ðlR1 þ P1 Þ2 þ ðQ 1 þ lÞ2 lP1 ð1  R0 Þ: Indeed, subtracting the expression on the right-hand side of the inequality from the expression on the left-hand side and simplifying yields

ðQ 1 þ lÞðlQ 1 þ R1 ÞðlR1 þ P1 Þ h i  ðlR1 þ P1 Þ2 þ ðQ 1 þ lÞ2 lP1 ð1  R0 Þ ¼ ðlR1 þ P1 ÞðQ 1 R1  P1 Þ þ ðQ 1 þ lÞ2 lP1 R0 2

þ ðQ 1 þ lÞl ðQ 1 R1  P1 Þ > 0; since

Q 1 R1  P1 ¼ ð1 þ  þ qÞðq þ  þ qÞ  qð1  dÞ ¼ ð þ qÞðq þ  þ qÞ þ  þ q þ qd P 0: Hence, all solutions of Eq. (3.28) have negative real parts when R0 < 1 and hence any eigenvalue with a positive real part, in view of growing solutions in the linear regime, can only come from

polynomial (3.27). But this is precisely Eq. (3.5) with us given in (3.8), taking into consideration the provisions of Remark 3.5. Hence the result in the relevant part of Theorem 3.1 carries over. Furthermore, at R0 ¼ 1, except for one eigenvalue which is zero, all other eigenvalues of the Jacobian matrix at the disease free equilibrium will have negative real parts. h Whenever R > 1, Theorem 3.4 provides conditions for the existence of endemic steady states. We expect the dynamics of solutions of the full nonlinear system near the endemic steady states to be influenced by the behaviour of the disease free system when R > 1, and the size of R0 . We recall from Proposition 3.9 that the disease free steady state can be driven to instability at the point in parameter space where k0 > kH0 . As the next result shows, whenever the endemic steady state is linearly unstable to small perturbations, there exists a parameter regime where it can only admit oscillatory instabilities. Proposition 3.10. Let the parameters be such that at least one endemic equilibrium solution, as postulated by Theorem 3.4, exists. Then any instabilities to small perturbations near such an endemic steady will be oscillatory instabilities whenever

/ P 0;

wi P

1 2

and R0 P 2

ð3:29Þ

where

/ ¼ P þ ac½ðI  k0 ð1  I ÞÞ þ 2k0 us ð1  I Þ; wi is the steady state solution of the questing infected mosquitoes e2 whose value is determined by the solution of Eq. (3.12), and R0 ¼ R 0 e 0 given by (3.6) is the unique disease threshold parameter. with R Proof. Linear stability of the endemic equilibrium Ee is determined by the eigenvalues of the Jacobian matrix, J evaluated at Ee ¼ ðI ; us ; v s ; ws ; ui ; v i ; wi Þ. Here, the characteristic equation is

 k þ ðbwi þ lÞ ½ðk þ 1Þðk þ Þðk þ qÞ þ dq½ðk þ 1Þðk þ cÞðk þ qÞ   cðak0 ð1  2us Þ þ aÞð1  I Þ þ ak0 ð1  2wi ÞcrI ðk þ ðbwi þ lÞÞðk þ 1Þðk þ Þðk þ qÞ

 þ bqrð1  I Þws ðk þ 1Þðk þ cÞðk þ qÞ  cðak0 ð1  2us Þ þ aÞ ¼ 0:

ð3:30Þ This can be rewritten in terms of the positive parameters P; P 1 ; Q ; Q 1 ; R and R1 and expanded to obtain

A7 k7 þ A6 k6 þ A5 k5 þ A4 k4 þ A3 k3 þ A3 k3 þ A2 k2 þ A1 k þ A0 ¼ 0; ð3:31Þ

Fig. 4. (a) Relationship between k0 and c. Stable solutions arise within the region labeled SS, while oscillatory solutions occur within region OS. (b) Relationship between the threshold parameters R and R0 and k0 . Notice that R0 attains a maximum of 1=us or ðbrÞ=ðlð1  dÞÞ. Since R depends on k0 linearly, a plot of R0 versus R is similar to that of R0 versus k0 .

56

C.N. Ngonghala et al. / Mathematical Biosciences 240 (2012) 45–62

Fig. 5. Plot of solutions to system (2.23) in the absence of endemic steady states depicting a stable disease free equilibrium for k0 ¼ 10; L ¼ 5  103 and N h ¼ 4  106 , with all other parameters as given on Table 3 and Eq. (4.1). For this parameter regime, R ¼ 30; R0 ¼ 0:7050, and usdfe ¼ v sdfe ¼ wsdfe ¼ 0:9667. In Figure (a) the initial data includes some infectious humans and infectious mosquitoes and in (b), the initial data does not include infectious individuals in any of the populations.

1.2

0.0005

1.0

I

0.0010

0.0000

0

5000 Time

0.8 9950

10000

1.6

0.010

1.0

0.005

0.4 9950

9975 Time

0.000

10000

0

9975 Time

10000

5000 Time

10000

Fig. 6. Time series plot of solutions to (2.23) for k0 ¼ 13; L ¼ 5  103 ; N h ¼ 4  106 , and all other parameters as shown on Table 3 and Eq. (4.1), showing oscillatory behavior for us ; v s and ws . For these parameters, R ¼ 39 and R0 ¼ 0:7106. Notice that each of the other variables relaxes to zero and that even though R0 < 1, we still have oscillatory solutions since k0 > kH0  12:8. The approximate period of oscillation is 5 days.

where A7 ¼ 1; A6 ¼ Q 1 þ Q þ bwi þ l; A5 ¼ ðQ 1 þ QÞðbwi þ lÞ þ QQ 1 þ R þ R1 ; A4 ¼ ðQ 1 Q þ R þ R1 Þðbwi þ lÞ þ P 1 þ Q 1 R þ R1 Q þ / þ ak0 crð2wi  1ÞI ; A3 ¼ ðQ 1 R þ R1 Q þ /Þðbwi þ lÞ þ ack0 rð2wi  1ÞI ððbwi þ lÞ þ Q 1 Þ þ P 1 ðbwi þ Q Þ þ /Q 1 þ R1 R; A2 ¼ ð/Q 1 þ R1 RÞðbwi þ lÞ þ ack0 rð2wi  1ÞI ðQ 1 ðbwi þ lÞ þ R1 Þ þ P 1 ðbQwi þ RÞ þ /R1 ; A1 ¼ /R1 ðbwi þ lÞ þ ack0 rð2wi  1ÞI ðP 1 ðbwi þ lÞ þ qÞ þ bP 1 Rwi þ /P1 ; /P1 ðbwi 

k0 ð2wi 

 1Þwi

A0 ¼ þ lÞ þ bqrðac / ¼ P þ caI  cak0 ð1  2us Þð1  I Þ ¼ P þ ac

þ Pðus Þ2 R ðR0  2ÞÞ; ½ðI  k0 ð1  I ÞÞ þ 2k0 us ð1  I Þ:

signs assures us that there are no positive real values for k that are zeros of (3.31). This indicates that any instabilities that can arise as a result of perturbations in parameters are oscillatory instabilities. At this stage, it is then a straightforward, but tedious exercise to apply the Routh-Hurwitz criteria on the characteristic polynomial (3.31) to see that there are indeed eigenvalues with positive real parts which will indicate growing oscillatory solutions in the linear regime. h Illustrative Example. Considering the endemic steady state presented in Eq. (3.18), we have

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1  dÞ R0  1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; I ¼ r þ ð1  dÞ R0  1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi lð1  dÞð R0  1Þ : wi ¼ br 

We note that ðbwi þ lÞI ¼ bwi and ð1  I Þws ¼ us and that bqrus ¼ lðqð1  dÞÞ ¼ lP 1 ; bqrQus ¼ lP 1 Q and bqrRus ¼ lP1 R, and then easily verify that all coefficients of the characteristic polynomial (3.31) are positive if / P 0; wi P 1=2 and R0 P 2. But I then / > 0 if I  k0 ð1  I Þ P 0 or if k0 6 1I  . Thus Descartes rule of

At this steady state,

us ¼

lð1  dÞ and br

57

C.N. Ngonghala et al. / Mathematical Biosciences 240 (2012) 45–62

0.4

1.0

0.2

0.5

0.0

0

750

1500

T ime

0.0

I us vs ws ui vi wi

0

750

1500

T ime

Fig. 7. Characterization of the endemic steady states in regions I and II of Fig. 2 for d ¼ 0:4; L ¼ 5  103 ; N h ¼ 105 ; 0:35 6 k0 6 0:6857, and the other parameters in Table 3 and Eq. (4.1). (a) Region (I) plot: For k0 ¼ 0:5 < kA01 ¼ 0:6857; R0 ¼ 3:2415; A1 ¼ 0:3714; A2 ¼ 0:0237; and D ¼ 0:2328. (b) Region (II) plot: For k0 ¼ 0:6857 ¼ kA01 ; R0 ¼ 5; A1 ¼ 0; A2 ¼ 0:0423; and D ¼ 0:1691.

2

2

1.5

1.5

1

1

0.5

0.5

0 0

200

400 600 Time

800

1000

0 0

10

20

Time

30

40

50

Fig. 8. Plot of solutions to system (2.23) for k0 ¼ 12; L ¼ 5  103 ; N h ¼ 105 . These and the rest of the parameters in Table 3 and Eq. (4.1) yield R ¼ 36; R0 ¼ 28:3631; A1 ¼ 0:9810; A2 ¼ 0:0322; and D ¼ 1:0909. Note that as predicted by our analysis, v s ¼ ws and ui ¼ v i ¼ wi .

0.83554

0.045

0.83547

0.035

0.83540 9950

9975

10000

0.025 9950

Time 1.070

0.2

1.015

9975

Time

10000

Time

0.4

0.0 9950

9975

10000

0.960 9950

9975

10000

Time

Fig. 9. Time series plots of solutions to system (2.23) showing oscillatory phenomena for k0 ¼ 14:8; L ¼ 5  103 ; N h ¼ 105 and the rest of the parameters are as given on Table 3 and Eq. (4.1). For these parameters, R ¼ 44:4; R0 ¼ 28:5164; A1 ¼ 0:9846; A2 ¼ 0:0323; and D ¼ 1:0987. The approximate period of oscillation is 5 days.

58

C.N. Ngonghala et al. / Mathematical Biosciences 240 (2012) 45–62

Fig. 10. Backward/forward bifurcation diagram illustrating the behavior of the full disease model for L ¼ 5  103 ; N h ¼ 105 ; 0:25 < k0 6 18:35, and the other parameters in Table 3 and Eq. (4.1) as the right half plane of Fig. 2 is traversed from region III through region V. Bounded oscillatory solutions that emerge via a Hopf bifurcation (HB) are c observed for k0 P 14:75. As k0 is reduced from kH1 0 ¼ 14:75, stable steady state solutions, which persist until when k0 ¼ k0 ¼ 0:35 emerge (see region LS on Fig. 10(b)). At k0 ¼ 0:35, a second endemic steady state, which is unstable emerges. Further decrease in k0 from 0:35 to 0:27 leads to bistability, whereby a stable endemic steady state coexists with a stable disease-free steady state (BB). Both endemic steady state branches disappear as k0 is reduced past 0:27, leaving only a disease-free steady state (see Fig. 10(b)). Note that R0 P 28:5142 when k0 P 14:75; R0 < 1 when k0 < 0:35, and that Fig. 10(b) is a continuation of Fig. 10(a).

Fig. 11. Characterization of the endemic steady states in regions I; II and III of Fig. 2 for d ¼ 0:4; 0:4 6 k0 6 18:4, and the other parameters in Table 3 and Eq. (4.1). (a) Plot of maximum and minimum values of wi against k0 . For 0:4 6 k0 < 0:6857; A1 < 0 and a stable endemic equilibrium solution is observed (Region I of Fig. 2). At k0 ¼ 0:6857; A1 ¼ 0 and again a stable endemic equilibrium solution is observed (Region II of Fig. 2). For k0 > 0:6857; A1 > 0 ((Region III of Fig. 2)). Within this region, a stable endemic equilibrium solution is observed as k0 is increased. A Hopf bifurcation occurs at kH1 0 ¼ 12:8 and as k0 is increased beyond 12:8, bounded periodic oscillations emerge. This plot represents a vertical cut through Fig. 4(a), where c is fixed at 0:7849, while k0 is varied from 0:4 to 18:4. (b) This plot represents a horizontal cut through Fig. 4(a), where k0 is fixed at 13:5 and c is varied from 0 through 1:6. Stable solutions are observed for 0 6 c < 0:5 (see region i of Fig. 4(a)). A Hopf bifurcation occurs at c ¼ 0:5 and further increase in c leads to oscillatory solutions (see region OS of Fig. 4(a)). Stable solutions reappear as c is increased beyond 1.2 (see region SS of Fig. 4(a)). Note that the point c ¼ 0:5 corresponds to the point where a horizontal line through Fig. 4(a) cuts the curve the first time while c ¼ 1:2 corresponds to when the same line cuts the curve the second time.

wi P 1=2 () P1þ



pffiffiffiffiffiffiffiffiffiffiffiffiffiffi

lð1  dÞð R0  1Þ 1 P () R0 br 2  2 2 br 1

2lð1  dÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1  dÞ R0  1 ¼ ;

r

¼1þ

q2

a

¼1þ

and if k0 6

2us  2 k0 r R0 P 1 þ 1d

4. Numerical simulations

I 1  I

in which case / > 0:

It is evident that there is richer dynamics in the region in parameter space where the restrictions in Proposition 3.10 are lifted for the steady state in (3.18) or even in the region where 0 6 R0 6 2. This can equivalently be shown for the other endemic steady states. However, since the goal of this present paper is to present a new route to modeling malaria dynamics, we do not investigate these regions further in this write up. To proceed then we use physically reasonable parameter 4 values as shown in Table 3 to illustrate Proposition 3.10.

4

The meaning of all the parameters in the model have been explained in Table 2.

Given the size of the system and the complex nature of the coefficients of Eq. (3.31), for the purposes of illustrating the potential of our model, we used MATLAB and the symbolic manipulation package MATHEMATICA to numerically evaluate and search the parameter space in view of exploring the nature of the solutions of the characteristic polynomial, and also having an insight into the possible solutions that exist for the full model. Our numerical evaluations indeed show that there exists parameter regimes whereby the steady state solutions can be driven to instability via a Hopf bifurcation. We investigated the relationship between the two threshold parameters R and R0 as the parameters c and k0 are varied. It is worth noting that in terms of the original system parameters, c depends on bv , the contact rate between a susceptible questing mosquito and a susceptible or infectious human and k0 is the growth rate of newly emerging mosquitoes. By varying c and computing the corresponding values of k0 ; R and R0 , the relationship between these threshold parameters and c can be determined. On the other hand, by varying k0 and computing the corresponding values of R and R0 , we can also determine a relationship between

C.N. Ngonghala et al. / Mathematical Biosciences 240 (2012) 45–62

these threshold parameters and k0 . See Fig. 4 for graphical illustrations of these. Fig. 4(b) illustrates that there are potentially four regions of interest for the threshold parameters that could be of great importance in potentially informing policy makers on control strategies. Using a MATLAB code, the full system (2.23) is simulated using the physically reasonable parameter values presented on Table 3 to illustrate conditions when oscillations occur in the dynamics of malaria transmission, as well as conditions when there is a backward bifurcation. Based on the parameter values shown on Table 3, and the dimensionless parameter groupings in Section 2.4, the values for q and l can be computed to be 0:809 and 0:008, respectively. In addition, from Remark 2.2, we see that  is between 0:0455 and 1 þ 21bh N h , while c is between 0:0455 and 1 þ 21bv N h . Then using these relationships and those for the scaled parameter values as given in Eq. (2.24)), as well as the fact that b < L=N h and ar < q, we can fix the values of the remaining scaled parameters as follows:

a ¼ 0:4636; d ¼ 0:8; c ¼ 0:7849;  ¼ 0:9; r ¼ 1:1667;

b ¼ 0:8L=Nh :

ð4:1Þ

These values will be used in the discussions below, unless otherwise stated (as will be in the case of d). When there is no disease in the community, Proposition 3.9 gives conditions for the stability of the disease free steady state. If we set L ¼ 5  103 and N h ¼ 4  106 , where N h is specifically selected such that R0 < 1, we obtain kH0  12:8, which is the value at which QR  PðR  1Þ ¼ 0 (see Remark 3.2). By choosing k0 appropriately and maintaining the other parameters as provided in Table 3 and Eq. (4.1), the non-trivial disease free equilibrium solution is either linearly stable or oscillatory. For instance, when k0 ¼ 10 (which is less than kH0  12:8), the above parameter regime yields R ¼ 30 and R0 ¼ 0:7050 and system (2.23) relaxes to a stable disease free steady state within the long time limit, as illustrated in Fig. 5. For k0 ¼ 13 > kH0 ; R ¼ 39 and R0 ¼ 0:7106 and we obtain oscillatory behavior around the non-trivial disease free equilibrium. Hence, kH0 is a critical value below which the disease free equilibrium is linearly stable and above which it is unstable with the instability giving rise to oscillatory phenomena in the susceptible mosquito populations. See Fig. 6 for a graphical demonstration of the oscillatory behavior. Note that, here, the infectious human and mosquito populations die out to their stable disease free values of zero, since R0 < 1. Thus, when the disease is not yet present in the community, even small introductions of infections in the community do not lead to successful transmissions. Next, we study the dynamics of the solutions to Eq. (2.23), when there exist at least one endemic disease steady state for the different regions as illustrated in Fig. 2, where either R0 > 1 ðA2 > 0Þ or R0 < 1 ðA2 < 0Þ (i.e respectively the top and bottom halves of Fig. 2) or with either A1 > 0 or A1 < 0 (equivalent to k0 < kA01 ð¼ qð1  dÞ=ðarÞÞ or k0 > kA01 , which are, respectively, the left and right halves of Fig. 2), using the physically reasonable parameter values presented on Table 3 and Eq. (4.1). Our focus will be on the regions (I to V) where there is at least one endemic disease steady state. First, we select parameters such that A1 6 0 and R0 > 1, so that we are in regions I (the second quadrant) and II (the line A1 ¼ 0 or the A2 axis) of Fig. 2. In particular, for d ¼ 0:4; L ¼ 5  103 ; N h ¼ 105 and the remaining parameters as given in Table 3 and Eq. (4.1), if k0 ¼ 0:5, then k0 < kA01 ¼ 0:6857; R0 ¼ 3:2415; A1 ¼ 0:3714; A2 ¼ 0:0237, and Fig. 7(a) shows that there exist a stable endemic disease steady state in region I. On the other hand, if k0 ¼ kA01 ¼ 0:6857, then R0 ¼ 5; A1 ¼ 0; A2 ¼ 0:0423, and

59

Fig. 7(b) shows that there exist a stable endemic disease steady state in region II. Note that the malaria disease establishes itself within the community in both cases. Next, we select parameters such that A1 > 0 and R0 > 1, so that we are in region III (first quadrant) of Fig. 2. In this region, there is critical value kH1 of k0 above which oscillatory dynamics are ob0 served and below which we have linearly stable endemic equilibrium solutions. Note that kH1 is not necessarily the same as kH0 , 0 the critical value of k0 for the case in which there is no biologically feasible endemic steady state. In particular, with the parameters as 3 H1 in Table 3 and Eq. (4.1), kH1 0 ¼ 14:75. For k0 ¼ 12 < k0 ; L ¼ 5  10 5 and N h reduced to N h ¼ 10 so that R0 > 1, we can compute A1 ¼ 0:9810; A2 ¼ 0:0322; R ¼ 36 and R0 ¼ 28:2239 and Fig. 8 shows a stable endemic equilibrium solution representing the situation in which malaria establishes itself within the community.  On the other hand, for k0 ¼ 14:8 > kH1 0 ; R ¼ 44:4; R0 ¼ 28:5164; A1 ¼ 0:9846; A2 ¼ 0:0323, and unstable endemic equilibrium solutions arise that are accompanied by the emergence of oscillatory solutions as shown in Fig. 9. The sustained bounded oscillations are due to a Hopf bifurcation that occurs at k0 ¼ kH1 0 . From observation, natural oscillations in the mosquito population introduced the oscillations in the malaria dynamics. Most of the oscillations in the infectious human population and the mosquito populations are in phase but with different amplitudes and thus the total mosquito population is also oscillating with a larger e 2 , where R e 0 is given by (3.6), amplitude. We showed that R0 ¼ R 0 so although the values of R0 computed here for the particular parameter regime may seem large ð< 29Þ, the corresponding e 0 are less than 6. values in terms of R Until now, the way in which oscillations in malaria dynamics has mostly been captured was by the inclusion of a seasonal forcing term or via delay. However, here, we have demonstrated that it is possible to capture oscillations in the dynamics of malaria transmission by appropriately interpreting the life style and behavior of the transmitting vector and incorporating the relevant aspects into the malaria model without the consideration of delay nor the inclusion of a seasonal forcing term.

5. Discussion, biological implications and conclusions 5.1. Discussion In this paper, we have developed a new SIS model for endemic malaria which explicitly integrates the demography and life style of the vector that transmits the disease together with its interaction with the human population, and incorporates the disease dynamics. Our model is an improvement over SIS type models for malaria transmission based on the Ross–Macdonald framework as follows:  In the Ross type models, Su þ Iu is just called N v ; Sw is just Sv . In our model, the demography and life style of the mosquito is important and so it is well emphasized. Only those female mosquitoes that need to feed and are questing are involved in disease transmission and they are appropriately accounted for. Additionally, only breeding vectors populate the susceptible class. This is true for fed vectors and breeding site vectors. The novelty is especially visible when control is taken into consideration because with our model we can say what happens when the questing or breeding site mosquitoes are reduced or even the fed mosquitoes, how that will impact the model dynamics. Questing mosquitoes are controlled through residual indoor spraying with insecticides or using insecticide impregnated bed and window nets.

60

C.N. Ngonghala et al. / Mathematical Biosciences 240 (2012) 45–62

 The basic reproduction number provides another comparison between a Ross type SIS model and our model. Using the notation adopted here, the basic reproduction number for a Ross type SIS model for malaria transmission takes the form

Rclassical ¼ 0

bh bv Nh Nv  ;

ð5:1Þ

lv rh þ lh

while the basic reproduction number for our model is R0 ¼

bv bh Nh Sw a a ðN Þ  v v h :q: : lv ðrh þ lh Þ av av ðNh Þ 1 þ ð1p1 Þbh Nh þ ðlv þ bh Nh Þðav þ av ðNh Þ þ lv Þ l v

ð5:2Þ It is evident that R0 < Rclassical since Sw 6 Nv ; 0 < q < 1, and 0

0<

a a ðN Þ   v v h ð1p1 Þbh N h þ ðlv þ bh Nh Þðav þ av ðNh Þ þ lv Þ av av ðNh Þ 1 þ l v

< 1: Hence, our model brings out an important scaling factor for the basic reproduction number R0 obtained from the classical SIS model so that R0 ¼ .Rclassical , where 0 6 . < 1 is a scaling factor easily ob0 tained from Eq. (5.2). This shows that the R0 for our new model is smaller than that computed for the Ross Macdonald-type SIS models since not all mosquitoes take part in the transmission process, a fact that has hitherto not been quantified in previous models.  Additionally, two threshold parameters have been identified in our model: one directly linked to the demography and life style of the vector (R ) and the other linked to the disease dynamics (R0 ), which as highlighted in Fig. 4 (b) illustrate four regions of interest for the threshold parameters that could be of great importance in potentially informing policy makers on control strategies.  Furthermore, our model captures the natural fluctuations known to occur in mosquito populations and malaria dynamics, as well as presents a new route to periodic oscillations in the dynamics of malaria transmission, a feature that has been shown [50] to be impossible with the Ross type models without recourse to external forcing (by means of periodic contact rates, delays, etc). These oscillations provide a plausible framework for mosquito control and hence, malaria control.  Our model also exhibits a backward bifurcation (an important phenomenon to consider when designing control schemes), highlighting the reason why it is more difficult to eradicate malaria. This indicates that when the disease is endemic, control measures that focus on reducing R0 so that R0 < 1 may not be successful unless R0 falls below a certain threshold value; i.e. R0 < Rbb 0 < 1. The difficulty in eradicating malaria may therefore lie in the complicated life style of the vector responsible for the transmission process, the complicated nature of the life style of the malaria parasite in its ability to share its life cycle so that part of it is in the mosquito and the other in the human. For example when on the right side of Fig. 2 (A1 > 0 which is equivalent to k0 > kA01 ), there are parameter regimes where the system goes from no endemic disease steady state through a situation in which a linearly stable and an unstable endemic disease steady state coexist with a stable non-trivial disease steady state, to the case where there exist only one linearly stable endemic steady state to the case where the endemic steady state is oscillating, as we traverse the right half plane from bottom to top or as we move through regions VII ! VI ! V ! IV ! III in that order. This is equivalent to k0 increasing bb from kA01 to kbb 0 (where R0 ¼ R0 < 1), to some critical value of k0 , say c k0 (where R0 ¼ 1), to the value kH1 0 , where a Hopf bifurcation exists

and oscillatory solutions emerge. See, for example, Fig. 10 that is generated using the parameter regime given on Table 3 and Eq. (4.1) and for L ¼ 5  103 ; N h ¼ 105 , with k0 increasing in the order c 0:25 which is less than kA01 ¼ 0:2286 ! kbb 0  0:27 ! k0  0:35 ! kH1  14:75. Note that Fig. 10 (b) is a continuation of 0 Fig. 10(a). If d ¼ 0:4 and R0 > 1 with the remainder of the parameters as given on Table 3 and Eq. (4.1) such that we are in the region where A2 > 0 (the top half of Fig. 2), Fig. 11 (a) shows the dynamics of the solutions as we traverse Fig. 2 in the order regions I ! II ! III. Notice that we move from a linearly stable endemic steady state H to oscillatory dynamics arising when k0 ¼ kH1 0 ¼ k0 ¼ 12:8. Note that for the reduced value of d here, the numerical values of the threshold parameters kH0 and kH1 0 at which Hopf bifurcations occur in the disease free and endemic disease systems coincide. Fig. 11 shows the solution curves against c as we traverse Fig. 2 in the same order of regions I ! II ! III, showing the system moving from linearly stable endemic solutions through limit cycles and back to linearly stable endemic constant solutions. This is equivalent to c moving across a fixed k0 in Fig. 4 (a) from region SS to OS and back to SS. An important question to consider is one involving the origin of the oscillatory dynamics as well as the backward bifurcation. There are many potential factors but, as indicated by our analysis and numerical simulations, we believe that both the oscillations and backward bifurcation are generated due to changes in the mosquito growth rate k0 , as well as the feedback mechanisms captured in the requirement that successfully fed and fertilized mosquitoes return to the breeding site to lay eggs at rate k0 and then start the process again, while unsuccessful ones are assumed killed. Now, there are several factors that affect k0 at different stages of the mosquito breeding cycle. In particular, the cyclic behavior of mosquitoes reproducing, in which mosquitoes must lay eggs, emerge, quest, feed and return to the breeding sites to lay more eggs affect k0 and the dynamics of the malaria disease and likely constitutes a basis for oscillations. If there is a healthy blood meal with many mosquitoes returning to the breeding site to breed and lay healthier eggs, k0 gets large and the mosquito population at the breeding site gets large. However, when these mosquitoes and the newly emerged mosquitoes leave the breeding site and go to the human habitat site to quest and feed on blood, the population of mosquitoes at the breeding site reduces and hence the cyclic and oscillatory patterns in the mosquito populations that drives the oscillations into the malaria dynamics. However, when k0 gets really small, the steady state at the breeding site is likely sustained by the batch of mosquitoes that are returning to the breeding site. These are probably the mosquitoes that drive the malaria infection in the system until a critical value of k0 ¼ kbb 0 below which the endemic steady state ceases to exist and the nontrivial disease free steady state is stable. It is worth noting that other parameters beside k0 detect the different regions in space. Note that the following inequality is always true: H kA01 < kbb 0 < k0 . But, there is a minimal value of d required to ensure A1 bb that k0 < k0 < kc0 < kH0 (which for our analysis was true when d ¼ 0:8 but not so when d ¼ 0:4). 5.2. Biological implications and conclusions Our analysis points to the fact that to fully understand the dynamics of mosquito populations, and consequently the transmission dynamics of the diseases they carry, it would be meaningful to characterize the demography and life style of the adult mosquito vector in a higher dimensional framework including adult mosquitoes in the different compartments shown in the model, as opposed to the standard two-stage (aquatic juveniles and flying adult) approach. The fact that the mosquito demography and life style has a multidimensional character confirms the

C.N. Ngonghala et al. / Mathematical Biosciences 240 (2012) 45–62

challenging nature of the malaria control problem and at the same time presents a wonderful opportunity for research and other tools in order to understand what this might mean for malaria control. In particular, control schemes need to be well thought out to capture all three phases of the mosquito population. A desirable control scheme should be one that target each of the three stages of the mosquito population. It should seek to limit contact between questing mosquitoes and humans and disrupt breeding sites or move them further away from human habitats. For example, when k0 > kA01 where we are on the right half of (Fig. 2), Fig. 10 indicates that as the growth rate of the mosquito population (k0 ), is reduced from a large value 18.35 to about 0.35 (indicated by regions III to IV on Fig. 2), the infectious human and mosquito populations (which both depend on wi ), move from stable oscillating solutions through limit cycles to linearly stable constant solutions. A further reduction of k0 from 0:35 to about 0:27 (regions IV to V to VI on Fig. 2), leads to the backward bifurcation scenario in which the linearly stable constant solution coexists with stable disease free steady state populations, in which the mosquito population is non-zero and oscillating. Thus, at high growth rates, control measures should take into consideration the oscillatory nature of the solutions and be implemented to capture the bottom swing of the infectious human and mosquito populations. Some common practical measures for controlling the malaria transmitting vector (the Anopheles sp.) are: (1) by spraying the adult mosquitoes and their larval forms with non-resistant insecticides and other chemicals, as well as the usage of insecticide treated bednets. This form of control though desirable, may pose problems to the public because some of the chemicals may be highly toxic to humans; (2) the use of biological control such as the introduction of a possible predator at the mosquito breeding sites so that the predators can feed on the returning and newly emerging mosquito vectors at the vector breeding sites, as well as their larval forms; (3) the use of preventive and sanitation schemes designed to limit or eliminate stagnant man-made water reservoirs such as those in empty bottles, cans, old tires and potholes from unpaved roads, as well as some natural ones from around human dwellings. We note that d is also an important parameter in itself and in that it affects k0 . It depend on N h ; bh , and p1 . Increasing (or decreasing) any of these parameters increases (decreases) d as well. However, for control purposes, a realistic parameter to vary in other to introduce a change in d is bh (the flow rate from infectious questing mosquito to susceptible humans), which can be reduced by reducing contacts between infectious questing mosquitoes and humans through the use of insecticide treated bednets, mosquito repellents and by indoor spraying, or by moving breeding sites further away from human habitats. The parameter c is also important in that it depends on bv (the flow rate from infectious humans to susceptible questing mosquitoes). An increase (decrease) in bv increases (decreases) c. Disease control through the reduction of bv is similar to that through the reduction of bh discussed above. For an effective and successful control campaign, our results suggest that these control measures should be continuous and sustained over longer periods until a small mosquito population is attained with small growth rates, because even a slight increase in the mosquito population might lead to the backward bifurcation situation in which we have a bistable endemic steady state and a disease free state coexisting. Note that within this region, a slight increase in the mosquito reproduction rate can push the mosquito population to the stable endemic steady state branch of the curve. Thus, there should be continuous monitoring and evaluation of the control measures and continuous education of the public to ensure that the disease burden is continuously reducing before eradication may be attained.

61

Our model suggests that to fully understand the dynamics of indirectly transmitted diseases such as malaria fever, dengue, yellow fever, for which the transmitting vector resides in different habitats from the humans, mathematical models should not be developed in isolation but should incorporate the demography and life style of the arthropods that transmit the agent that cause these diseases. From our results, we have only identified regions in parameter space where different interesting solutions can be found. The complete characterization of the solutions in these regions and inclusion of more realistic factors such as the creation of breeding sites in view of spatial dynamics, the extension of the model to the full SEIRS malaria model that also includes a disease related death rate in the human population and others are under investigation. Additionally, the impact of climate related factors to the oscillations is under investigation as well. For now, we have presented (as well as analytically and numerically analyzed) here, a framework that researchers can build upon and apply to other indirectly transmitted diseases and hope it provides information that could potentially lead to eradicating malaria or significantly reducing the disease burden. Acknowledgments CNN acknowledges the support of the National Institute for Mathematical and Biological Synthesis (NIMBioS), an Institute sponsored by the National Science Foundation, the U.S. Department of Homeland Security, and the U.S. Department of Agriculture through NSF Award #EF-0832858, with additional support from The University of Tennessee, Knoxville, a WOC grant from the Society of Mathematical Biology (SMB), and travel grants from West Virginia University (WVU). GAN acknowledges the support of the University of Buea Grant for the 2008/2009 academic year. MIT-E acknowledges support of NSF Award # OISE-0855380. The authors thank Kenneth Showalter of WVU and Suzanne Lenhart of NIMBioS and University of Tennessee, Knoxville for their valuable comments. References [1] Z. Hawass, Y.Z. Gad, S. Ismail, et al., Ancestry and pathology in King Tutankhamun’s family, JAMA 303 (7) (2010) 638. [2] World Health Organisation, The World Malaria Report 2009, WHO Press (2010,) Available at http://www.who.int/malaria/world_malaria_report_2010/ en/index.html. Assessed August 27, 2011. [3] K.E. Mace, M.F. Lynch, J.R. MacArthur, S.P. Kachur, L. Slutsker, R.W. Steketee, T. Popovic, The opportunity for and challenges to malaria eradication, morbidity & mortality weekly report, Centers for Dis. Control Prevention (CDC) 60 (15) (2011) 476. [4] I.W. Sherman, Malaria: Parasite biology, Pathogenesis, and Protection, ASM Press, 1998. [5] R. Ross, The Prevention of Malaria, John Murray, London, 1911. [6] G. Macdonald, The analysis of infection rates in diseases in which superinfection occurs, Trop. Dis. Bull. 47 (1950) 907. [7] K. Dietz, L. Molineaux, A. Thomas, A malaria model tested in the Africa Savanna, Bull WHO 50 (1974) 347. [8] T.J.N. Bailey, The Mathematical Theory of Infectious Diseases and Its Application, Griffin, London, 1975. second ed.. [9] J.L. Aron, Mathematical modeling of immunity to malaria, Math. Biosci. 90 (1988) 385. [10] J.L. Aron, R.M. May, The population dynamics of malaria, in: R.M. May (Ed.), Population Dynamics of Infectious Diseases, Chapman and Hall, London, 1982. [11] J.C. Koella, On the use of mathematical models of malaria transmission, Acta Trop. 49 (1991) 1. [12] G.A. Ngwa, W.S. Shu, A mathematical model for endemic malaria with variable human and mosquito populations, Math. Comput. Model. 32 (7) (2000) 747. [13] G.A. Ngwa, C.N. Ngonghala, N.B.S. Wilson, A model for endemic malaria with delay and variable populations, J. Cameroon Acad. Sci. 1 (3) (2001) 169. [14] N. Chitnis, J.M. Cushing, J.M. Hyman, Bifurcation analysis of a mathematical model for malaria transmission, SIAM J. Appl. Math. 67 (1) (2006) 24. [15] N. Chitnis, J.M. Hyman, J.M. Cushing, Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model, Bull. Math. Biol. 70 (2008) 1272. [16] M.I. Teboh-Ewungkem, C.N. Podder, A.B. Gumel, Mathematical study of the role of gametocytes and an imperfect vaccine on malaria transmission dynamics, Bull. Math. Biol. 72 (1) (2010) 63.

62

C.N. Ngonghala et al. / Mathematical Biosciences 240 (2012) 45–62

[17] M.I. Teboh-Ewungkem, T. Yuster, A within-vector mathematical model of plasmodium falciparum and implications of incomplete fertilization on optimal gametocyte sex ratio, J. Theory Biol. 264 (2010) 273. [18] M.I. Teboh-Ewungkem, M. Wang, Male fecundity and optimal gametocyte sex ratios for P. falciparum during incomplete fertilization, J. Theory Biol. 307 (2012) 183. [19] J. Labadin, C.M.L. Kon, S.F.S. Juan, Deterministic malaria transmission model with acquired immunity, Proc. World Congress Eng. Comput. Sci. 2 (2009) 779. [20] Ana Paula P. Wyse, L. Bevilacqua, M. Rafikovd, Simulating malaria model for different treatment intensities in a variable environment, Ecol. Model. 206 (2007) 322. [21] H.M. Giles, D.A. Warrel, Bruce–Chwatt’s Essential Malariology, fourth ed., Hodder Arnold Publication, London, 2002. [22] P. Martens, R.S. Kovats, S. Nijhof, P. de Vries, M.J.T. Levermore, D.J. Bradley, J. Cox, A.J. McMichael, Climate change and future populations at risk of malaria, Global Environ. Change 9 (1999) S89. [23] G. Zhou, N. Minakawa, A.K. Githeko, G. Yan, Association between climate variability and malaria epidemics in the East African highlands, Proc. Natl. Acad. Sci. 101 (2004) 2375. [24] J.A. Patz, S.H. Olson, Malaria risk and temperature: influences from global climate change and local land use practices, Proc. Natl. Acad. Sci. 103 (15) (2006) 5635. [25] P. Reiter, Global warming and malaria: knowing the horse before hitching the cart, Malaria J. 7 (Suppl 1: S3) (2008). [26] J.L. Gallup, J.D. Sachs, The economic burden of malaria, Am. J. Trop. Med. Hyg. 64 (2001) 85. [27] S. Nourridine, M.I. Teboh-Ewungkem, G.A. Ngwa, A mathematical model of the population dynamics of disease transmitting vectors with spatial consideration, J. Biol. Dynam. (2011) 1751. [28] S.G. Staedke, E.W. Nottingham, J. Cox, M.R. Kamya, P.J. Rosenthal, G. Dorsey, Proximity to mosquito breeding sites as a risk factor for clinical malaria episodes in an urban cohort of Ugandan children, Am. J. Trop. Med. Hyg. 69 (3) (2003) 244. [29] M.I. Teboh-Ewungkem, Malaria control: The role of local communities as seen through a mathematical model in a changing population-Cameroon, Advances in Disease Epidemiology, Nova Science Publishers, 2009 (Chapter 4, pp. 103– 140). [30] J. Nedelman, Introductory review some new thoughts about some old malaria models, Math. Biosci. 73 (2) (1985) 159. [31] P.J. McCall, D.W. Kelly, Learning and memory in disease vectors, Res. Update. Trends Parasito 18 (10) (2002) 429.

[32] G.A. Ngwa, On the population dynamics of the malaria vector, Bull. Math. Biol. 68 (8) (2006) 2161. [33] G.A. Ngwa, A.M. Niger, A.B. Gumel, Mathematical assessment of the role of non-linear birth and maturation delay in the population dynamics of the malaria vector, Appl. Math. Comput. 217 (2010) 3286. [34] W.E. Ricker, Stock and recruitment, Can. J. Fisheries Aquatic Sci. 11 (1954) 559. [35] R.J.H. Beverton, S.J. Holt, On the dynamics of exploited fish populations, vol. 19 of 2, Fisheries investment, London, Her Majesty’s Stationary Office, 1981. [36] K. Cook, P. van den Driessche, X. Zou, Interaction of maturation delay and nonlinear birth in population and epidemic models, J. Math. Biol. 39 (1999) 332. [37] Å. Brännström, D.J.T. Sumpter, The role of competition and clustering in population dynamics, Proc. R. Soc. B 272 (2005) 2065. [38] R.M. Anderson, R.M. May, Population biology of infectious diseases: Part I, Nature 280 (1979) 361. [39] R.M. Anderson, Mathematical and statistical study of the epidemiology of HIV, AIDS 3 (6) (1989) 333. [40] H.R. Thieme, Epidemic and demorgraphic interaction in the spread of a potentially fatal disease in a growing population, Math. Biosci. 111 (1992) 99. [41] H.K. Hale, Ordinary Differential Equations, John Wiley & Sons, New York, 1969. [42] P. van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci. 180 (2002) 29. [43] O. Diekmann, J.A.P. Heesterbeek, J.A.J. Metz, On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations, J. Math. Biol. 28 (1990) 365. [44] W. Huang, K.L. Cook, C. Castillo-Chavez, Stability and bifurcation for a multiple-group model for the dynamics of hiv/aids transmission, SIAM J. Appl. Math. 52 (3) (1992) 835. [45] J. Dushoff, W. Huang, C. Castillo-Chavez, Backward bifurcations and catastrophe in simple models of fatal diseases, J. Math. Biol. 36 (1998) 36:227. [46] C.M. Krisb-Zeleta, J.X. Velasco-Hernandez, A simple vaccination model with multiple endemic states, Math. Biosci. 164 (2) (2001) 183. [47] P. van den Driessche, J. Watmough, Reproductive numbers and sub-threshold endemic equilibria for compartment models of disease trsmission, Math. Biosci. 180 (2001) 29. [48] C. Castillo-Chavez, B. Song, Dynamical models of tuberculosis and their applications, Math. Biosci. Eng. 1 (2) (2004) 361. [49] Central Intelligence Agency, Country Comparison: Life expectancy at birth, The world fact book (Assessed September 2010). [50] J.C.B.S. Ruan, D. Xiao, On the delayed Ross–Macdonald model for malaria transmission, Bull. Math. Biol. 70 (4) (2008) 1098.