The transmission dynamics of a within-and between-hosts malaria model

Ecological Complexity 38 (2019) 31–55

Contents lists available at ScienceDirect

Ecological Complexity journal homepage: www.elsevier.com/locate/ecocom

Original Research Article

The transmission dynamics of a within-and between-hosts malaria model ⁎,a

b

F.B. Agusto , M.C.A. Leite , M.E. Orive a b

a

T

Department of Ecology and Evolutionary Biology, University of Kansas, Lawrence, KS, USA Mathematics and Statistics Program, University of South Florida St. Petersburg, St. Petersburg, FL, USA

ARTICLE INFO

ABSTRACT

Keywords: Malaria Within-host Between-host Feedback functions Backward bifurcation Oscillations

In this paper, we developed a novel deterministic coupled model tying together the effects of within-host and population level dynamics on malaria transmission dynamics. We develop within-host and within-vector dynamic models, population level between-hosts models, and a nested coupled model combining these levels. The unique feature of this work is the way the coupling and feedback for the model use the various life stages of the malaria parasite both in the human host and the mosquito vector. Analysis of the coupled and the within-human host models indicate the existence of locally asymptotically stable infection- and parasite-free equilibria when the associated reproduction numbers are less than one. The population-level model, on the other hand, exhibits backward bifurcation, where the stable disease-free equilibrium co-exists with a stable endemic equilibrium. A global sensitivity analysis was carried out to measure the effects of the sensitivity and uncertainty in the various model parameters estimates. The results indicate that the most important parameters driving the pathogen level within an infected human are the production rate of the red blood cells from the bone marrow, the infection rate, the immunogenicity of the infected red blood cells, merozoites and gametocytes, and the immunosensitivity of the merozoites and gametocytes. The key parameters identified at the population level are the human recovery rate, the death rate of the mosquitoes, the recruitment rate of susceptible humans into the population, the mosquito biting rate, the transmission probabilities per contact in mosquitoes and in humans, and the parasite production and clearance rates in the mosquitoes. Defining the feedback functions as a linear function of the mosquito biting rate, numerical exploration of the coupled model reveals oscillations in the parasite populations within a human host in the presence of the host immune response. These oscillations dampen as the mosquito biting rate increases. We also observed that the oscillation and damping effect seen in the within-human host dynamics fed back into the population level dynamics; this in turn amplifies the oscillations in the parasite population within the mosquito-host.

1. Introduction

Maier et al., 2009; White, 2017; Wiser, 2011). The injected sporozoites migrate to the liver invading the liver cells; they develop into schizonts which rupture releasing merozoites and thus marking the end of the exoerythrocytic phase. The erythrocytic phase starts with the released merozoites entering the blood stream and infecting the red blood cells (Baron, 1996; Centers for Disease Control and Prevention, 2015; Hoffman et al., 2012; Maier et al., 2009; Wiser, 2011). After about 48 hours, the infected cells rupture releasing merozoites which further invade the red blood cells to renew the cycle and maintain the infection. The merozoites are the asexual form of the parasite; however, some of the released merozoites differentiate into the sexual form known as gametocytes (Baron, 1996; Hoffman et al., 2012; Maier et al., 2009; Wiser, 2011). The gametocytes, when ingested by a susceptible mosquito, lead to the development of the parasite in the mosquito (Baron, 1996; Centers for Disease Control and Prevention, 2015; Hoffman et al.,

Malaria is a vector-borne disease caused by the Plasmodium parasite transmitted to humans through the bite of an infected mosquito. In 2015, an estimated 429,000 deaths occurred due to malaria infection from an estimated 212 million cases (World Health Organization, 2017). Most of the deaths occurred among children living in Africa where a child dies on average every minute from malaria (World Health Organization, 2017). Malaria in humans develops in two phases; the first is the exoerythrocytic phase (involving the liver cells) and the second is the erythrocytic phase (involving the red blood cells) (Baron, 1996; Centers for Disease Control and Prevention, 2015; Wiser, 2011). Malaria infection begins when an infected mosquito takes a blood meal from a susceptible human and thereby injects the malaria parasite (sporozoites) into the otherwise healthy human (Hoffman et al., 2012;

⁎

Corresponding author. E-mail address: [email protected] (F.B. Agusto).

https://doi.org/10.1016/j.ecocom.2019.02.002 Received 30 June 2017; Received in revised form 23 January 2019; Accepted 16 February 2019 1476-945X/ © 2019 Published by Elsevier B.V.

Ecological Complexity 38 (2019) 31–55

F.B. Agusto, et al.

2012; Maier et al., 2009; Wiser, 2011). In the study of infectious disease, recent work has focused on linking two vital levels of disease dynamics; this work has indicated that interactions between these levels may play an important role in shaping the host-pathogen relationship (Barfield et al., 2015; Cen et al., 2014; Day and Alizon, 2011; Feng et al., 2013; 2012; Gilchrist and Coombs, 2006; Gilchrist and Sasaki, 2002; Mideo et al., 2008; 2011; Park et al., 2013). These two levels are the between-host level (the host population level), which considers the epidemiological process of the disease transmission between hosts, and the within-host level, which focuses on pathogen dynamics within individual hosts and the immunological process between host cells and pathogens (Cen et al., 2014; Feng et al., 2013; 2012). In the case of malaria, these levels have the additional complexity of involving two hosts (humans and mosquitoes), and thus there are two dynamic systems to be considered at each level. A number of models (especially those for viral infections such as HIV) have studied the complex dynamics within hosts (see Beauchemin, 2006; Frost et al., 2001; Funk et al., 2005; Kelly et al., 2003; Murall et al., 2012; Orive et al., 2005 and references therein). Linking these types of models with classic epidemiological models at the host population level in the form of nested models requires careful consideration of how the dynamics of pathogen numbers within hosts may impact the transitions that take place at the higher level. The vast majority of the nested models developed have been used to determine the tradeoffs between epidemiological parameters, such as how a biological parameter important at both levels may experience conflicting selection pressure, and in evaluating the consequences of this conflicting selection pressure at different scales (Barfield et al., 2015; Gilchrist and Coombs, 2006; Gilchrist and Sasaki, 2002). For example, Feng et al. (2013, 2012) showed that new properties from the complex coupled system of the between-host epidemic process and within-host dynamics can emerge and unexpected dynamics such as backward bifurcations can occur (Cen et al., 2014). For malaria transmission, the two processes (within-host dynamics and population-level epidemic interactions) have been studied independently. The dynamics at the within-host level are considered in models in Chiyaka et al. (2008), Hellriegel (1992), Li and Ruan (2011), Niger and Gumel (2011) and Tumwiine et al. (2008), while the models in Agusto et al. (2013, 2012), Agusto and Lenhart (2013), Aneke (2002), Chitnis et al. (2006) and Niger and Gumel (2008) consider the epidemic process at the population level. The linkage of the dynamics of these two levels is a challenging task, and there is no standard way on how to model it. Legros and Bonhoeffer recently (Legros and Bonhoeffer, 2016) developed a coupled system considering the evolution of resistance to antimalarial drugs; this model coupled the withinhuman host dynamics to the population level dynamics using the gametocytes, the non-pathogenic form of the parasite responsible for the infection in mosquitoes and not in humans as suggested by the authors. In this paper, we develop a deterministic coupled model linking the dynamics of a within-host model with population-level dynamics, and consider the effects of the coupling for malaria transmission dynamics. The novel feature of this model lies in how the coupling and feedback for the model have been developed using the various life history stages (notably sporozoites, merozoites, and gametocytes) of the malaria parasite, and in consideration of both between- and within-host dynamics in both the human host and the mosquito vector. To the best of our knowledge, this is the first study to develop a coupled model involving two hosts, and in particular, it is the first model to explicitly couple within-host dynamics to an epidemic process at the population level for malaria using both human and mosquito hosts (see the model formulation section below). We observed from the theoretical and epidemiological results of the coupled, the within-human host, and population-level models that the models have locally asymptotically stable infection- and parasite-free equilibria when the associated reproduction numbers are less than one. In addition, we observe that the population-level model, but not the

within-host model, exhibits backward bifurcation. The results of the global sensitivity analysis indicate that the most important parameters driving the pathogen level within an infected human are the production rate of the red blood cells from the bone marrow, the infection rate, the immunogenicity (ability to provoke an immune response) of the infected red blood cells, merozoites and gametocytes, and the immunosensitivity (susceptibility to clearance by the immune system) of the gametocytes and merozoites. The key parameters identified at the population level are the human recovery rate, the death rate of the mosquitoes, the recruitment rate of susceptible humans into the population, the mosquito biting rate, the transmission probabilities per contact in mosquitoes and in humans, and the parasite production and clearance rates in the mosquitoes. Numerical exploration of the coupled model defining the feedback functions as a linear function of the mosquito biting rate, reveals oscillations which dampen with increasing biting rate amidst the parasite populations within a human host in the presence of the host immune response. 2. Model formulation 2.1. The within-human host model 2.1.1. The within-human host model without immune response The total population of uninfected red blood cells (RBCs) X(t) and infected red blood cells (IRBCs) Y(t), within human host at time t, denoted by NWH(t) is given as:

NWH (t ) = X (t ) + Y (t ). and the total population of the gametocytes G(t) and merozoites, M(t) within the human host at time t, is given by

NMG (t ) = M (t ) + G (t ). The population of the RBCs is increased through a constant production from the bone marrow at the rate πX. It is decreased through an infection from the merozoite at rate of β and decay at rate of μX. Hence the rate of change of the RBCs is given as:

X (t ) =

X

X (t ) M (t )

µX X (t ).

The population of the infected cells Y(t) is increased through a new infection from the merozoite at the rate β. The infected cells burst at the rate μY, releasing an average number of merozoites (usually 8–32 merozoites per bursting). Furthermore, some of the parasites differentiate and develop into the sexual parasite (some remain asexual parasites) at the rate c. Hence the equations for the rate of change of the infected cells are given as:

Y (t ) =

X (t ) M (t )

µ Y Y (t )

cY (t ).

The merozoites mature in the infected cell before they rupture (burst). Thus, the population of merozoites is increased through the average number of merozoites released (r) through the bursting of infected cells at the rate rμY. It is decreased through natural death at the rate μM and successful invasion of the RBCs at rate uβX(t)M(t). The parameter u is used to account for the invasion of the RBCs by the parasite, if u = 0, the invasion is unsuccessful and if u = 1, the invasion is successful (Iggidr et al., 2006; Niger and Gumel, 2011). In this study, we set u=1, indicating that invasion has taken place. In malaria endemic regions, individuals are continually been bitten by infected mosquitoes (Robert et al., 2003). Therefore, infected individuals who already have blood-stage parasitemia are continually receiving liver-tropic Plasmodium sporozoites, thereby increasing the population of merozoites (M(t)). This process is incorporated in the feedback function L(IM), where IM gives the number of infectious mosquitoes and IM > 0. Hence, the population of merozoites is increased by the release of the merozoites from the liver cells due to the feedback from the population level infection (i.e. infected mosquitoes). 32

Ecological Complexity 38 (2019) 31–55

F.B. Agusto, et al.

Thus, in the absence of the feedback (IM = 0 ), when the parasite load is independent of the population level infection, the parasite load within the human host reaches the steady state value (M*(0)) (Cen et al., 2014; Gilchrist and Coombs, 2006). Hence, the equation for the rate of change for the merozoites are written as follows:

M (t ) = L (IM ) + rµ Y Y (t )

µM M (t )

(i) Many different immune cells activated by the presence of the malaria parasite are assumed to be classified together in one population for simplicity, and this malaria-specific immune response is denoted by B(t). (ii) The parasites (gametocytes and merozoites) are eliminated from circulation by immune cells at rate kG, kM respectively and the merozoites are further reduced by infecting the RBCs at the rate of β. In particular, the human host acquired natural antibodies against gametocytes which help limit its development within the host circulation and interrupt the gamete development and fertilization within the mosquito after their ingestion (Gebru et al., 2017; McQueen et al., 2013; Naotunne et al., 1993; Sutherland, 2009). (iii) IRBCs are removed from the blood stream by immune cells at rate kY. (iv) Production of immune cells are stimulated by the presence of IRBCs, gametocytes and merozoites at rates ρY, ρG, ρM, respectively and deactivated at a constant rate (μB).

u X (t ) M (t ).

Following Cen et al. (2014), several functional forms can be considered L (IM ) = M IMz L (IM ) = M IM , for L(IM), such as and L (IM ) = M1 IM /(1 + M 2 IM ), where θM, θM1, θM2 > 0 and z < 1. We will consider the simplest case where L (IM ) = M IM . All the above functional forms have the following properties

L (0) = 0, L (IM )

0, L (IM ) > 0, L (IM )

0.

While we will here assume that the feedback function L(IM) is positive, so that the population of merozoites is increased with the number of infectious mosquitoes, it is possible to imagine situations where a negative relationship exists. For example, infection may trigger immune responses preventing other strains from establishing, leading to a type of priority effect (Portugal et al., 2011), and potentially leading to a negative relationship between IM(t) and M(t). Other possible negative feedbacks exists; for example, if parasite load increases mosquito mortality (such as was observed in a laboratory study by Dawes et al., 2009), there may be a negative feedback between the numbers of sporozoites within the mosquito host (SP(t), at the within mosquito host level) and the death rate for mosquitoes (μV, at the population level for mosquitoes). A more complex model for within-mosquito dynamics and the coupling of within- and between-host dynamics for the mosquito host will be a further extension of the model presented here. The population of gametocytes is increased when a proportion of the asexual parasite (found within the infected red blood cells) differentiate into the sexual form of the parasite (gametocytes) rather than merozoites; this occurs at rate c (Babiker et al., 2008; McQueen et al., 2013). The population of gametocytes is decreased through natural death at rate μG. Thus,

G (t ) = cY (t )

The above assumptions lead to the following system of differential equations, where (′) denotes derivative with respect to t.

X (t ) M (t )

Y (t ) =

X (t ) M (t )

M (t ) = L (IM ) + rµY Y (t ) G (t ) = cY (t )

µG G (t ),

X (t ) M (t )

µX X (t ),

µY Y (t )

µG G ( t ) YY

(t ) +

cY (t )

µM M (t )

kY B (t ) Y (t ), u X (t ) M (t )

kM B (t ) M (t ),

kG B (t ) G (t ), M M (t )

+

G G (t )]

µB B (t ). (2.2)

The current model extends the models in Chiyaka et al. (2008), Niger and Gumel (2011), Iggidr et al. (2006), Anderson et al. (1989) and Santos and Torres (2013) as follows: (i) the current model contains immune response but the model in Iggidr et al. (2006) does not. (ii) the current model incorporates the role of immune cells in suppressing the production of gametocytes and merozoites but models in Niger and Gumel (2011), Anderson et al. (1989) and Santos and Torres (2013) do not. (iii) the current model has a compartment for gametocytes while models in Chiyaka et al. (2008), Niger and Gumel (2011) and Anderson et al. (1989) do not.

cY (t ),

µM M (t )

X (t ) M (t )

B (t ) = B (t )[

µX X (t ),

µY Y (t )

Y (t ) =

G (t ) = cY (t )

Hence, the model describing the dynamics of within-human host infection is given by the system of ordinary differential equations. X

X

M (t ) = L (IM ) + rµY Y (t )

µG G (t ).

X (t ) =

X (t ) =

u X (t ) M (t ),

2.2. The within-vector model

(2.1)

where (′) denotes the derivative with respect to t. The current model extends the models in Chiyaka et al. (2008), Niger and Gumel (2011), Anderson et al. (1989) and Santos and Torres (2013) by incorporating a compartment for gametocytes that is not included in the models presented in Chiyaka et al. (2008), Niger and Gumel (2011) and Anderson et al. (1989). One of the important biological features of this model is that the within-human host dynamics occur on a faster time scale than the dynamics of the population level model; in other words, the betweenhuman host model (considered in the population level model section) occurs at a much slower time scale. These multiple time scales allow us to study the dynamical properties of the model by using the fast- and slow-systems framework, (see the time-scales of coupled model section for further discussions).

In this section, we consider a within-vector model that focuses on the sporozoites production within a mosquito host (see Baton and Ranford-Cartwright, 2005; Beier, 1998; Centers for Disease Control and Prevention, 2015; Kappe et al., 2003; Sreenivasamurthy et al., 2013; Wiser, 2011 and references within for details on sporozoites production). After the gametocytes (microgametocytes (male) and macrogametocytes (female)) are ingested by the mosquito, the microgametes penetrate the macrogametes forming zygotes. The zygotes, in turn, develop into ookinetes which invade the midgut wall of the mosquito where they develop into oocysts. The oocysts grow and rupture, releasing the sporozoites which make their way to the mosquitoes salivary glands (Baton and Ranford-Cartwright, 2005; Beier, 1998; Centers for Disease Control and Prevention, 2015; Kappe et al., 2003; Sreenivasamurthy et al., 2013; Wiser, 2011). In this within-mosquito model, we ignore the developmental stages of the parasite and consider only a simple birth and death model for the development of the sporozoites. However, we consider the feedback from the population level infection (i.e. infected humans) to describe further sporozoites formation and production following further uptake of infected blood meal by

2.1.2. The within-human host model with immune response We extend model system (2.1) by including malaria specific immune response. Following Chiyaka et al. (2008), we make the following assumptions on the model. 33

Ecological Complexity 38 (2019) 31–55

F.B. Agusto, et al.

an infected mosquitoes (Ferguson and Read, 2004; Teboh-Ewungkem, 2010). Thus, in the absence of feedback (IH = 0, the number of infectious humans), when the parasite load is independent of the population level infection, the parasite load within the mosquito host reaches the steady state value (SP* (0) ) (Cen et al., 2014; Gilchrist and Coombs, 2006). On the other hand, with feedback (IH > 0) there is a continuous release of the sporozoites. We let H(IH) take the simplest case, H (IH ) = H IH . As with L(IM), we can consider several other functional forms such as H (IH ) = H IH , H (IH ) = H IHz and H (IH ) = H 1 IH /(1 + H 2 IH ), where θH, θH1, θH2 > 0 and z < 1. These functional forms have the following properties

H (0) = 0, H (IH )

0, H (IH ) > 0, H (IH )

SH (t ) =

(µS +

S

t

µS +

(

H

[

H

µH SH (t ),

(M (t )) + µH +

H (M (t ))] IH (t ),

+ µH ) RH (t ), (2.4)

SH (0) = SH 0,

IH (0) = IH 0,

RH (0) = RH 0 .

The parameter πH in (2.4) is the recruitment rate into the population. In (2.4), the parameter μH is the natural death rate for all the human subpopulations, ωH is the rate at which recovered individuals lose their acquired immunity. The function βH(SP(t)) is the sporozoites-dependent transmission probability of an infected mosquitoes infecting susceptible individuals (Baton and Ranford-Cartwright, 2005; Centers for Disease Control and Prevention, 2015; Wiser, 2011), and this is given as H (SP (t ))

(2.3)

=

(2.5)

h SP (t ),

where βh is a constant. Note that in this paper the aim is to explore implications of coupling the dynamics of a within-host model to a population-level dynamics for malaria transmission dynamics. Thus, the functional form (2.5) was chosen to be linear to keep the model as simple as possible, so that it is amenable to theoretical analysis. Other functional forms such as the saturating function below can be considered, H (SP (t ))

=

h1 SP (t )

1+

h 2 SP (t )

,

where βh1 and βh2 are constants. The parameter bM is the mosquito biting rate. The corresponding parasite-dependent disease-induced death rate which is a function of the merozoites (since the pathology and clinical manifestations of malaria are due in part to the merozoites Baton and Ranford-Cartwright, 2005; Cunnington et al., 2013; Wiser, 2011) and this is given by

Lemma 1. Let the initial data SP(0) ≥ 0. Then the solution SP of the withinmosquito host model (2.3) is non-negative for all t > 0. Furthermore, S

(M (t )) IH (t )

H

SH (t )

SH (t )

with initial conditions given as

Note that (′) denotes the derivative with respect to t. The parameter πS is the recruitment rate of the sporozoites into the mosquito salivary glands (Baton and Ranford-Cartwright, 2005; Beier, 1998; Centers for Disease Control and Prevention, 2015; Kappe et al., 2003; Sreenivasamurthy et al., 2013; Wiser, 2011). The surviving sporozoites in the salivary glands undergo further modifications necessary to infect the human host (Sreenivasamurthy et al., 2013), thus, we assume that μS is the clearance (decay) rate of the sporozoites in the salivary glands and αS is the injection (depletion) rate of the sporozoites into a susceptible humans when the individual mosquito is taking a blood meal. In this section we will investigate the dynamic properties of model (2.3) in the absence of the feedback H (IH ) = 0 (i.e., the within-vector host parasite load is independent of the population level infection).

lim sup SP (t ) =

NBH (t )

NBH (t )

RH (t ) =

0.

S ) SP (t ).

H (SP (t )) bM IM (t )

H RH (t )

+

H (SP (t )) bM IM (t )

IH (t ) =

We again here assume that the feedback function is positive, so that there is an increase in the number of sporozoites within mosquitoes with an increase in the number of infected human hosts. Hence, the equation for the rate of change for the number of sporozoites (SP) are written as follows:

SP (t ) = H (IH ) +

H

H (M (t ))

.

(2.6)

= rH M (t ),

where δ is the parasite cost coefficient which is equal to the increase in host mortality rate caused by a single unit of parasite(merozoites) reproducing at rate rH = 1 (Gilchrist and Sasaki, 2002). The corresponding parasite-dependent recovery rate γH(M(t)) of the infectious individuals is given by

S

As in the within-human host model section, the within-vector dynamics occurs on a faster time scale than the dynamics of the population level model, i.e. the between-human host model occurs at a much slower time scale. The time-scales of coupled model section gives further details on the fast-slow framework of this model.

H (M (t ))

=

h

1 + M (t )

,

(2.7)

where γh is a scaling factor. The equations for the mosquitoes is given as follows

2.3. The population level between-hosts malaria model

SM (t ) =

For the population level between-hosts (human and mosquitoes) malaria model, we assume that individuals in the community are homogeneous both in their mixing and in their response to the infection. We further assume that the infected human classes are related to the within-human host dynamics of a particular individual. Thus, the total human population, NBH(t), at time t is split into mutually-exclusive sub-populations of individuals who are susceptible (SH(t)), infectious (IH(t)) and recovered (RH(t)), so that

SM (t )

M

IM (t ) = SM (t )

M (G (t )) bM

M (G (t )) bM

IH (t ) NBH (t )

IH (t ) NBH (t )

µV SM (t )

µV IM (t ).

(2.8)

with initial conditions given as

SM (0) = SM 0,

IM (0) = IM 0.

In (2.8), parameter πM is the mosquito recruitment rate while the natural death rate for mosquitoes is give as μV. The function βM(G(t)) is the parasite-dependent transmission probability of infected individuals infecting a susceptible mosquito, this is a function of the gametocytes (Baton and Ranford-Cartwright, 2005; Centers for Disease Control and Prevention, 2015; Wiser, 2011) and this is given as

NBH (t ) = SH (t ) + IH (t ) + RH (t ). Similarly, the total mosquito population size at time t, denoted by NM(t), is sub-divided into susceptible mosquitoes (SM(t)) and infectious mosquitoes (IM(t)), so that NM (t ) = SM (t ) + IM (t ) . The between-hosts malaria transmission dynamics is given by the following deterministic system of ordinary differential equations, where (′) denotes the derivative with respect to t.

M (G (t ))

=

m G (t ),

(2.9)

The function βM(G(t)) can also be represented by the saturating function such as: 34

Ecological Complexity 38 (2019) 31–55

F.B. Agusto, et al.

M (G (t ))

=

m1 G (t )

1+

m 2 G (t )

,

where βm, βm1 and βm2 are constants. Hence, this formulation via Eqs. (2.5), (2.6), (2.7) and (2.9) explicitly links the within-host and between host models, since the infection transmission probabilities from humans to mosquitoes and mosquitoes to humans both depend on the sporozoites and gametocytes respectively, while the death of the infected individuals and movement into the recovered class both depend on the number of merozoites in the infected individual. 2.4. The coupled within- and between-hosts malaria model Combining the Eqs. (2.2), (2.3) and (2.4)-(2.8) leads to the following system of ordinary differential equations for the coupled withinand between-hosts malaria model.

X (t ) =

X

X (t ) M (t )

Y (t ) =

X (t ) M (t )

µX X (t ),

µY Y (t )

M (t ) = L (IM ) + rµY Y (t ) G (t ) = cY (t ) B (t ) = B (t )[

SP (t ) = H (IH ) + SH (t ) = IH (t ) =

+

H

µM M (t )

µG G (t ) YY

(t ) +

M M (t )

NBH (t ) H

SM (t ) =

M

(M (t )) IH (t )

IM (t ) = SM (t )

SM (t )

kM B (t ) M (t )

+

u X (t ) M (t ),

NBH (t ) SH (t )

(

H

µB B (t ), Fig. 1. Flow diagram of the coupled model (2.10). The parameters indicated at the population level in the figure are the only ones used for linking the two levels of the model.

S ) Sp (t ),

[

H

SH (t )

µH SH (t ),

(M (t )) + µH +

Table 1 Description of the variables and parameters of the model (2.10).

H (M (t ))] IH (t ),

+ µH ) RH (t ),

M (G (t )) bM

M (G (t )) bM

G G (t )]

H (SP (t )) bM IM (t )

H RH (t )

H (SP (t )) bM IM (t )

RH (t ) =

kY B (t ) Y (t ),

kG B (t ) G (t ),

(µS +

S

cY (t )

IH (t ) NBH (t )

IH (t ) NBH (t )

µV SM (t )

µV IM (t ). (2.10)

The flow diagram of the model is depicted in Fig. 1, and the associated state variables and parameters are described in Table 1. 3. Analysis of the model To analyze the dynamic properties of model (2.10), we first investigate the properties of the isolated models (the within-human host in Eq. (2.2) and between-host models in Eq. (2.4)), followed by the coupled model in Eq. (2.10). 3.1. The within-human host model In this section we will investigate the dynamic properties of the within-human host model (2.2). First, we will consider model (2.2) in the absence of the feedback L (IM ) = 0 (i.e., the within-human host parasite load is independent of the population level infection) and then relax this assumption in subsequent analysis. 3.1.1. Basic properties of the within-human host model Positivity and boundedness of solutions For the within-human host malaria model (2.2) to be epidemiological meaningful in the absence of the feedback L (IM ) = 0, it can be shown (using the method in Appendix A of (Thieme, 2003)) that all its state variables are non-negative for all time. In other words, solutions of the model system (2.2) with non-negative initial data will remain nonnegative for all time t > 0. 35

Variable

Description of state variables

X(t) Y(t) G(t) M(t) B(t) SP(t) SH(t) IH(t) RH(t) SM(t) IM(t) Parameters πX μX, μY μG, μM μB β r kY, kG, kM Variable c ρY, ρG, ρM u πS μS αS πH βH βM bM ωH μH γH δH πM μV

Population of uninfected RBCs Population of infected RBCs Population of gametocytes Population of merozoites Population of malaria-specific immune cells Population of sporozoites Population of susceptible individuals in the community Population of infected individuals Population of recovered individuals in the community and hospital Population of susceptible mosquitoes Population of infectious mosquitoes Description Production rate of RBC from the bone marrow Natural death rate of uninfected and infected RBCs Natural death rate of gametocytes and merozoites Death rate of immune cells Rate of infection Average number of merozoites released per bursting IRBCs Immunosensitivity of parasites, gametocytes and merozoites Description of state variables Rate of differentiation into gametocytes from the asexual IRBCs Immunogenecity of IRBCs, gametocytes and merozoites Modification parameter that accounts for RBC invasion Sporozoites recruitment rate into the mosquito salivary glands Sporozoites clearance (decay) rate in the salivary glands Sporoziotes injection (depletion) rate into a susceptible humans Human recruitment rate Transmission probability of infection from humans to mosquitoes Transmission probability of infection from mosquitoes to humans Mosquito biting rate Loss of malaria immunity Natural death rate in humans Recovery rate of infection humans Malaria induced death rate Mosquitoes recruitment rate Natural death rate in mosquitoes

Ecological Complexity 38 (2019) 31–55

F.B. Agusto, et al.

Lemma 2. Let the initial data X(0) > 0, Y(0) ≥ 0, G(0) ≥ 0, M(0) ≥ 0 and B(0) ≥ 0. Then the solutions (X, Y, G, M, B) are all non-negative for all t > 0.

within a susceptible host (Anderson et al., 1989; Diekmann et al., 1990; Van den Driessche and Watmough, 2002; Hethcote, 2000). Thus, Lemma 4 implies that malaria can be eliminated from within a host (when WH <1) if the initial sizes of the sub-populations of the model (2.2) are in the basin of attraction of the DFE, WH .

The proof of Lemma 2 is given in Appendix A. Invariant regions Since the model (2.2) monitors the RBCs, IBRCs, parasite and the human host immune response against the parasite, all variables and parameters of the model are non-negative. The model (2.2) will be analyzed in a biologically-feasible region as follows. Consider the biological feasible region: WH

3.1.3. Within-human host model parasite-present equilibrium (PPE) The parasite-present equilibrium when B = 0 The within-human host model (2.2) parasite-present equilibrium, 1 WH , in the absence of the immune response (when B = 0 ) is given as 1 WH

5 +

=

where,

where k1 = µY + c and ˜ WH | B =0 =

=

(X (t ), Y (t ), M (t ), G (t ), B (t )) M

¯, G M

G¯ , B (t ) <

5 +:

NWH (t )

X

µmin

,

Lemma 5. The PPE, ˜ WH | B = 0 > 1.

,

NWH (t ) = X + Y ,

¯ = M

rµ Y X , µM µmin

G¯ =

c

X

µG µmin

, µmin = min{µX , µY }.

, 0, 0, 0, 0 .

The linear stability of WH can be established using the next generation operator method on the system (2.2). Using the notations in Van den Driessche and Watmough (2002), and taking Y, M, G, as the infected compartments, the matrices F and V, for the new infection terms and the remaining transfer terms, are, respectively, given by,

0 F= 0 0

X* 0 0 0 , 0 0

and,

V=

uk1)

.

is locally asymptotically stable (LAS) if

The proof of Lemma 6 is given in Appendix E. The dynamics of the within-human host model (2.2) in the presence of immune response is de2 picted in Fig. 2; in which case the PPE WH is locally asymptotically stable. Note that the simulations shown in Fig. 2 show values over 150 days to match the timescale in Fig. 3; longer timescales show values tending towards the equilibrium. The values used for most parameters are the same as the baseline values in Table 2 except X = 50, = 0.0005, c = 0.017, u = 0.013, µY = 1/40, µB = 1/25, Y = 0.001, G = 0.001, M = 0.001, and the initial conditions are taken as X (0) = 500, Y (0) = 10, M (0) = 1, G (0) = 1, B (0) = 1. However, the purpose of the immune response is to defend the infected human against the infection by preventing the proliferation of the parasite. Hence, we observe in Fig. 3 the elimination and clearance of the parasite due to the actions of the immune response and the eventual deactivation of the malaria-specific immune response (B(t)) once the parasites have all been cleared. The values used for most parameters are the same as the baseline = 0.0005, c = 0.017, u = 0.013, values in Table 2 except X = 500, µY = 1/40, µB = 3/50, Y = 0.001, G = 0.001, M = 0.001. Note, that we have chosen these values to show the elimination and clearance of the parasite in the infected individual which theoretically implies the instability 2 of the PPE WH . Thus, depending on the parameter set chosen, the PPE 2 may be stable or unstable. WH

3.1.2. Within-human host model parasite-free equilibrium (PFE) The within-human host model (2.2) has a PFE, obtained by setting the right-hand sides of the equations in the model to zero, given by X

k1 µM µX

2 , is locally asymptotically stable (LAS) if the Lemma 6. The PPE, WH 2 2 ) of model (2.2) evaluated at Jacobian matrix ( J WH WH have negative eigenvalues.

The proof of Lemma 3 is given in Appendix B.

µX

X (rµY

The parasite-present equilibrium when B ≠ 0 The within-human host model (2.2) has a parasite-present equili2 brium ( WH ) in the presence of the immune response provided B⁎⁎ ≠ 0 and M M ** + G G** < µB . The existence of the parasite-present equi2 librium WH is given in Appendix D.

5 Lemma 3. The region WH = + is positively-invariant for the model (2.2) with non-negative initial conditions in +5 .

= (X *, Y *, M *, G *, B *) =

1 WH ,

1) cµX µM ( ˜ WH 1) , ,0 . (rµY k1 u) µG

The proof of Lemma 5 is given in Appendix C.

with,

WH

= (X ** , Y **, M **, G **, B **) µ µ ( ˜ WH 1) µX ( ˜ WH X = , X M , (rµY k1 u) µX ˜ WH

µY + c + kY B * 0 0 rµY µM + kM B * + u X * 0 c 0 µG + k G B *

It follows that the basic reproduction number of the within-human host model (2.2), denoted by WH , is given by

The expression WH is the number of secondary infections within a completely susceptible individual due to the injection into a healthy individual of a single infected sporozoite. Using Theorem 2 in Van den Driessche and Watmough (2002), the following result is established.

3.1.4. Within-human host model with feedback L(IM) ≠ 0 In this section, we investigate the dynamic properties of model (2.2) in the presence of the feedback L (IM ) = M IM 0, (i.e., there is feedback from the population level infection into the within-human host parasite load). The within-human host model (2.2) has a parasite-present equilibrium ( 3WH (IM ) ) in the presence of the immune response which is dependent on the feedback from the population level infection. The existence of the parasite-present equilibrium 3WH (IM ), is given in Appendix F.

Lemma 4. The PFE of the within-human host model (2.2), given by WH , is locally asymptotically stable (LAS) if WH < 1, and unstable if WH > 1.

Lemma 7. The PPE 3WH (IM ) is locally asymptotically stable (LAS) provided B* ≠ 0 and Y Y ** (IM ) + M M ** (IM ) + G G ** (IM ) < µB .

The threshold quantity ( WH , i.e. the basic reproduction number) measures the average number of new infections generated by a parasite

The proof of Lemma 3.1.4 is given in Appendix G. In the next section we will investigate the dynamic properties of

WH

=

(FV 1) =

(µM

r µY X * + kM B * + u X *)(µY + c + kY B *)

36

Ecological Complexity 38 (2019) 31–55

F.B. Agusto, et al.

NBH (t ) = SH (t ) + IH (t ) + RH (t ) and NM (t ) = SM (t ) + IM (t ).

Invariant regions Since the model (2.4)-(2.8) monitors human and vector populations, all variables and parameters of the model are non-negative. The model (2.4)-(2.8) will be analyzed in a biologically-feasible region as follows. Consider the feasible region BH

=

H

3 +

M

×

2 +,

with, H

=

(SH (t ) + IH (t ) + RH (t ))

M

=

(SM (t ) + IM (t ))

2 +:

3 +:

H

NBH (t ) M

NM (t )

µV

µH

,

.

3 2 Lemma 9. The region BH = H M +× + is positively-invariant for the model (2.4)-(2.8) with non-negative initial conditions in +5 .

The proof of Lemma 9 is given in Appendix H. 3.2.2. The between-hosts model disease-free equilibrium (DFE) The disease-free equilibrium of model (2.4)-(2.8) is obtained by setting to zero the right-hand sides of the model equations in the absence of infection. This equilibrium is given as BH

Fig. 2. Simulation of within-human host model (2.2) in the presence of immune 2 response with stable PPE WH . (a) Uninfected red blood cell; (b) Infected red blood cell; (c) Merozoites; (d) Gametocytes; (e) Immune response; (f) Combined within-human host dynamics. Parameter values used are as given in Table 2. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

0 U=

BH

and lim sup NM (t ) t

M

µV

,0 .

* M (G **) bM SM

N BH *

0

,W=

H (M **)

+ µH + 0

H (M **)

0 . µV

=

(UW 1) =

2 * ** * M (G **) bM SM H (SP ) SH 2 * ) µV [ H (M **) + µH + H (M **)] (NBH

Further, using Theorem 2 in Van den Driessche and Watmough (2002), the following result is established. Lemma 10. The DFE of the between-host model (2.4)-(2.8), given by BH , is locally-asymptotically stable (LAS) if BH < 1, and unstable if BH > 1. The threshold quantity BH (i.e. the basic reproduction number) measures the average number of new infections generated by a single infected individual in a completely susceptible population (Anderson et al., 1989; Diekmann et al., 1990; Van den Driessche and Watmough, 2002; Hethcote, 2000). Thus, Lemma 10 implies that malaria can be eliminated from the population (when BH <1) if the initial sizes of the sub-populations of the model (2.4)-(2.8) are in the basin of attraction of the DFE, BH .

Lemma 8. Let the initial data SH(0) ≥ 0, IH(0) ≥ 0, RH(0) ≥ 0, SM(0) ≥ 0, IM(0) ≥ 0. Then the solutions (SH, IH, RH, SM, IH) of the between-host model (2.4)-(2.8) are non-negative for all t > 0. Furthermore, H

M

µV

It follows that the basic reproduction number of the between-host model (2.4)-(2.8), denoted by BH , is given by

3.2.1. Basic properties of the between-hosts model Positivity and boundedness of solutions For the between-hosts malaria model (2.4)-(2.8) to be epidemiological meaningful, it can be shown, following the approach in the section above, that all its state variables are non-negative for all time. In other words, solutions of the model system (2.4)-(2.8) with non-negative initial data will remain non-negative for all time t > 0.

µH

, 0, 0,

* H (SP**) bM SH

N BH *

3.2. The between-hosts model

t

H

µH

The linear stability of BH can be established using the next generation operator method. Using the notations in Van den Driessche and Watmough (2002), the matrices U and W, for the new infection terms and the remaining transfer terms, are, respectively, given by,

model (2.4)-(2.8). It should be noted that population level betweenhosts dynamics occur on a much slower time scale compare to both the within-vector and within-human hosts dynamics; following (Cen et al., 2014; Gilchrist and Coombs, 2006) we will ignore the transient dynamics of the parasite-cell interactions and focus only on the steady state dynamics of the within-vector and within-human hosts models (2.3) and (2.4)-(2.8). As a result the variables SP(t), G(t) and M(t) are assumed to be at the equilibrium (i.e. SP**, G ** and M⁎⁎) and the transmission probabilities and the disease induced death rate become ** H (SP (t )) = H (SP ), M (G (t )) = M (G **) and H (M (t )) = H (M **) .

lim sup NBH (t )

* , IM *) = = (SH* , IH* , RH* , SM

3.2.3. The between-hosts model endemic equilibrium points (EEP) Next, we consider the dynamics of the endemic equilibrium points of model (2.4)-(2.8). We take advantage of the different time scales of the within- and between-hosts dynamics; as in Cen et al. (2014) and Feng et al. (2013, 2012), we suppose that the dynamics of the withinhuman and within-mosquito host (the fast dynamics, see within-host

,

with, 37

16 0.5 0.06 0.06

Rate of differentiation into gametocytes from the asexual IRBCs Natural death rate of uninfected RBCs Natural death rate of infected RBCs Natural death rate of merozoites

Immunosensitivity of merozoites

Average number of merozoites released per bursting IRBCs Modification parameter that accounts for RBC invasion Immunogenecity of IRBCs

Immunogenecity of gametocytes

kM

r u ρY

38

Recovery rate of infection humans Mosquitoes recruitment rate Natural death rate in mosquitoes Humans to mosquitoes transmission probability scaling factor Mosquitoes to humans transmission probability scaling factor

γH πM μV βh βm

μH δH

Natural death rate in humans Malaria induced death rate

Within host feedback function (L(IM)) scaling factor Within vector feedback function (H(IH)) scaling factor Sporozoites recruitment rate into the mosquito salivary glands Sporozoites clearance (decay) rate in the salivary glands Sporoziotes injection (depletion) rate into a susceptible humans Human recruitment rate Mosquito biting rate Loss of malaria immunity

θm θh πS μS αS πH bM ωH

ρG

Immunogenecity of merozoites

Immunosensitivity of gametocytes

Natural death rate of gametocytes Death rate of immune cells Immunosensitivity of parasites

ρM

kG

μG μB kY

5

0.0092 0.525 1/18 0.356 0.32135

1/(65 × 365) 0.0003454

0.02 0.5 400 1/5 1/10 400 0.2980 0.0005275

0.60

0.030

0.06

3.4722 × 10 1/16 1/20 0.030

0.015 1/120 1/2

0.09

12020750

c μX μY μM

Rate of infection

Production rate of RBC from the bone marrow

β

πX

Baseline value

Parameter description

Parameter

Table 2 Values and ranges of the parameters of the models (2.2), [(2.4)-(2.8)], and (2.10).

2.5 × 10 8) cells /day

8

8

0.12) cells /day

0.06) cells /day

8

8

1.2) /day

0.12) /day

(2.7 × 10

3

0.64) /day

(1 × 10 15 4.1 × 10 4) /day (0.0014 0.017) /day (0.05 1.0) /day (1/21 1/3) /day (0.072 0.64) /day

(5.5 × 10 5 1.1 × 10 2) /day (1/(80 × 356) 1/(58 × 365)) /day

(2 × 10 8 0.12) /day (0 0.04) /day (0 1) /day (10 800) /day (0.05 1.0) /day (2/25 3/25) /day (10 800) /day (0.1 0.8) /day

(3 × 10

(2 × 10

(10 8 0.06) cells /day (12.0 20.0) /day (0 1) /day

(10

(10

(2.7778 × 10 5 4. 1666 5) /day (1/20 3/40) /day (1/25 3/50) /day

(2.5 × 10 10 0.12) /day (0.017 0.022) /day (0.00667 0.01) /day (0.0025 0.50) /day

(1

Range

Chiyaka et al. (2008), Niger and Gumel (2008) Laperriere and Rubel (2011), Lou and Zhao (2010), Rubel et al. (2008) Chiyaka et al. (2008), Laperriere and Rubel (2011), Niger and Gumel (2008), Rubel et al. (2008) Chiyaka et al. (2008), Niger and Gumel (2008) Chitnis et al. (2006)

Chiyaka et al. (2008), Niger and Gumel (2008) Chiyaka et al. (2008), Niger and Gumel (2008)

Estimated Estimated Estimated Estimated Estimated Variable Chiyaka et al. (2008), Niger and Gumel (2008) Chiyaka et al. (2008), Niger and Gumel (2008)

Hellriegel (1992), Li and Ruan (2011), Tumwiine et al. (2008)

Hellriegel (1992), Li and Ruan (2011), Tumwiine et al. (2008)

Hellriegel (1992), Li and Ruan (2011), Tumwiine et al. (2008) Niger and Gumel (2011) Hellriegel (1992), Li and Ruan (2011), Tumwiine et al. (2008)

Hellriegel (1992), Li and Ruan (2011), Tumwiine et al. (2008)

Hellriegel (1992); Li and Ruan (2011); Tumwiine et al. (2008)

Hellriegel (1992) Hellriegel (1992), Li and Ruan (2011), Tumwiine et al. (2008) Hellriegel (1992), Li and Ruan (2011), Tumwiine et al. (2008)

(1992) (1992), Li and Ruan (2011), Tumwiine et al. (2008) (1992), Tumwiine et al. (2008) (1992)

Hellriegel (1992), Li and Ruan (2011), Tumwiine et al. (2008)

Hellriegel Hellriegel Hellriegel Hellriegel

Hellriegel (1992), Li and Ruan (2011), Tumwiine et al. (2008)

Reference

F.B. Agusto, et al.

Ecological Complexity 38 (2019) 31–55

Ecological Complexity 38 (2019) 31–55

F.B. Agusto, et al.

model analysis above) are at the endemic equilibrium when S and limSP (t ) = . t (µS + S ) Let 1 BH

WH

>1

**, IM **), = (SH**, IH**, RH**, SM

be an arbitrary endemic equilibrium point of the between-hosts model (2.4)-(2.8) and let

** H

=

**) bM IM ** , NH**

H (SP

** M

**

M (G **) bM IH

=

(3.1)

NH**

be the force infections for susceptible humans and susceptible mosquitoes at steady-state, respectively. Solving the equations at steadystate as functions of the force infections gives: SH** = IH** = RH** =

** = SM

(

H

+ µH )[

( H + µH )[ H (M **) + µH + H (M **)] ** H (M **) + µH + H (M **)]( H + µH ) H

(

H

+ µH )[

H (M **)

+ µH +

(

H

+ µH )[

H (M **)

+ µH +

** ( H

H

H

** H (M **) H H

+ µH )

** H (M **)]( H

+ µH )

** H (M **) H H

+ µH )

** H (M **) H H

** H H (M **) H

M

** = IM

** + µV ) ( M

** H (M **)]( H

** M M . ** + µV ) µV ( M

(3.2) Substituting (3.2) into (3.1) and simplifying, leads to the following quadratic equation (3.3)

a0 ( H**)2 + b0 ( H**) + c0 = 0, where a0 = µV H [k2 + H (M **)](k2 µV + k2 bM M (G **) + H (M **) µV ), b0 = k2 bM µV H k1 k2 1 +

bM H (SP**) M H (M **) H k2 µV H k1

bM M H (SP**) µV H

Fig. 3. Simulation of within-human host model (2.2) in the presence of immune 2 response with unstable PPE WH . (a) Uninfected red blood cell; (b) Infected red blood cell; (c) Merozoites; (d) Gametocytes; (e) Immune response; (f) Combined within-human host dynamics. Parameter values used are as given in Table 2 with µB = 1/38. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

M (G **)

+ 2k1 µV2 H [k2 + H (M **)]}, c0 = k1 µV2 H (1

2 ) BH

with k1 = H (M **) + µH + H (M **), k2 = H + µH . The coefficient a0 of the quadratic (3.3) is always positive, however the coefficient c0 is positive or negative depending on whether BH is less than or greater than one. If it is greater than one, then c0 is positive otherwise c0 is negative. Thus, the following results have been established:

Remark 1. Case (iii) of Theorem 1 does not hold for mild diseaseinduced death rate (i.e. 0 ≤ δH(M⁎⁎) < μH). This case is relevant when the malaria induced death rate is greater than natural death rate. To reiterate, case (iii) is possible for (δH(M⁎⁎) > μH such that c BH < BH < 1).

Theorem 1. The model (2.4)-(2.8) has (i) Precisely one unique endemic equilibrium if b0 = 0, c0 < 0 (or equivalently BH > 1); (ii) Precisely one unique endemic equilibrium if b0 < 0, and either c0 = 0 or b02 4a 0 c0 = 0; (iii) Precisely two endemic equilibria if c0 > 0, b0 < 0 and b02 4a 0 c0 > 0; (iv) No endemic equilibrium otherwise.

3.3. Analysis of the coupled malaria model In this section we investigate the dynamics of the coupled model (2.10) using the linear functions L (IM ) = M IM , H (IH ) = H IH , H (SP (t )) = h SP (t ), M (G (t )) = m G (t ) and H (M (t )) = rH M (t ) .

Note that Case (iii) indicates the possibility of backward bifurcation in the model (2.4)-(2.8) when BH < 1. To find the backward bifurcation when BH < 1, set the discriminant b02 4a 0 c0 to zero and solve for the critical value of BH denoted by c BH

=

1

k 22 {bM

M (G **)[k2 µV H k1

4[µV

+ bM

H (k2

+

M H (G **)( H (M **) H

H (M **))(k2 bM M (G **)

Thus, backward bifurcation will occur for values of BH < 1. The above result is summarized below.

k1 k2)] + 2k1 µV2

+ k2 µV +

c BH

3.3.1. Time-scales of coupled model (2.10) The important biological features of the coupled model (2.10) are H [k2

H (M **) µV )]

+

H (M **)]}

2

.

the various time-scales of the infection at the within-human host level, within-vector host level and at the population level between humans and mosquitoes. The infection dynamics within-human host occur on a faster time scale than the dynamics at the population level since the death rate of the infected red blood cells are greater than the recovery rate in humans (i.e., μY > γH). In other words, the infection turnover in humans occurs at a faster rate than the turnover at the population level.

such that

Theorem 2. The model (2.4)-(2.8) undergoes a backward bifurcation when Case (iii) of Theorem 1 holds and cBH < BH < 1. The phenomenon of backward bifurcation is illustrated numerically in Fig. 4. 39

Ecological Complexity 38 (2019) 31–55

F.B. Agusto, et al.

Fig. 4. Backward bifurcation plot for the population level dynamics. (a) Infected humans; (b) infected mosquitoes. Parameter values used are as given in Table 2.

Thus, at the within-human host level, the time scale is faster compare to the much slower time-scale at the population level leading to the fastand slow-systems framework. In a similar manner, the time-scale at the within-vector dynamics occurs on a much faster time scale than the dynamics at the population level model since the sporozoites clearance rate in the mosquito salivary glands is greater than the mosquitoes natural death rate (i.e., μS > μV). The formulation of the coupled model (2.10) involves the feedback from infected vector to within-human level and from the infected humans into the within-vector level via the functions L(IM) and H(IH). Given, this unique feature of the model, implementing the fast and slow framework requires that the death rate of the infected red blood cells be greater than the natural death rate in mosquito (i.e., μY > μV) for the feedback from infected mosquitoes into within-human dynamics. For the feedback from infected humans into the within-vector level, the fast and slow framework require the sporozoites clearance rate in the mosquito salivary glands be greater than the recovery rate in humans (i.e., μS > γH). These fast-slow assumptions are satisfied by the baseline parameter values utilized in the numerical analyzes given below, listed in Table 2

(b) The disease free equilibrium of the population level model (2.4)(2.8) is locally asymptotically stable when the associated reproduction number is less than one; (c) The population-level model (2.4)-(2.8) exhibits backward bifurcation, where the stable disease-free co-exists with a stable endemic equilibrium, when the associated reproduction number is less than one (i.e., cBH < BH < 1 when δH(M⁎⁎) > μH); (d) The coupled model (2.10) has a locally asymptotically stable infection-free equilibrium (IFE); 4. Sensitivity analysis Deterministic model outputs are governed by its input parameters, which exhibit some uncertainty in the process of their selection. Exploration of the complex exploratory models developed above requires determination of which within- and between-host parameters most strongly affect key model outputs (such as the number of infected hosts or pathogen numbers within hosts). A global sensitivity analysis was carried out using Latin Hypercube Sampling (LHS) and partial rank correlation coefficients (PRCC) in order to assess the impact of parameter uncertainty and the sensitivity of these key model outputs of the numerical simulations to variations in each parameter. Global sensitivity analyses were carried out separately for the within-human host model (2.2), the population level model (2.4)-(2.8), and the coupled model (2.10). LHS is a stratified sampling without replacement technique which allows for an efficient analysis of parameter variations across simultaneous uncertainty ranges in each parameter (Blower and Dowlatabadi, 1994; Marino et al., 2008; McKay et al., 2000; Sanchez and Blower, 1997). PRCC measures the strength of the relationship between the model outcome and the parameters, stating the degree of the effect that each parameter has on the outcome (Blower and Dowlatabadi, 1994; Marino et al., 2008; McKay et al., 2000; Sanchez and Blower, 1997). The LHS-PRCC method is described in Appendix J. To generate the LHS matrices, we assume that all the model parameters are uniformly distributed. A total of 1000 simulations (runs) of the models for the LHS matrix were then carried out, using the parameter values given in Table 2 (with the stated ranges) and as response functions, the numbers of parasites (merozoites (M) and gametocytes (G)) both for the within-human host and coupled models. The total number of infected humans (IH) and mosquitoes (IM) were also used as response functions for population level and coupled models. These simulation runs were followed by the parameter ranking using PRCC. Parameters with the largest partial rank correlation coefficient (PRCC) values have the largest impact on the model outcome examined under each analysis. Fig. 5(a) and (b) depict the PRCC values for each

3.3.2. Stability of the infection-free equilibrium (IFE) Let 0 C

* , IM *) = (X *, Y *, M *, G*, B *, SP*, SH* , IH* , RH* , SM =

X

µX

, 0, 0, 0, 0, 0,

H

µH

, 0, 0,

M

µM

,0 ,

be the infection-free equilibrium of the coupled model (2.10); we establish in the Lemma 11 below the stability of the infection-free equilibrium. Lemma 11. The infection-free equilibrium (IFE) of the coupled model (2.10), given by C0 , is locally asymptotically stable (LAS) provided C

=

( X X * + µM g1)2 + 4rµY X X * X X * + µM + g1

< 1,

where g1 = µY + c . The proof of Lemma 11 is given in Appendix I. We summarized below some of the theoretical and epidemiological findings of the within-human host model (2.2), the population level model (2.4)-(2.8) and the coupled model (2.10): (a) The equilibria (parasite-free and-present) of the within-human host model (2.2) are locally asymptotically stable (LAS) when their associated reproduction number is less than one and unstable otherwise; 40

Ecological Complexity 38 (2019) 31–55

F.B. Agusto, et al.

Fig. 5. PRCC values of the within-human host model (2.2) (panels (a) and (c)) and coupled model (2.10) (panels (b) and (d)). (a) Total number of merozoites using the within-human host model; (b) Total number of merozoites using the coupled model; (c) Total number of gametocytes using within-human host model; (d) Total number of gametocytes using coupled model. Parameter values and ranges used are as given in Table 2. Vertical lines give PRCC values significantly different from zero with p-value = 0.01; negative values smaller than these values and positive values greater than these values are significant.

parameter of the within-human host model (2.2) and the coupled model (2.10) using the total number of merozoites (M) as the response function. We set the PRCC cutoff points at ± 0.0821 and ± 0.0827 respectively for the two models, corresponding to a p-value of 0.01 (as in Marino et al., 2008; cut-off values are indicated with vertical lines in Fig. 5(a) and (b)). Bars extending to the left (for negative values) or to the right (for positive values) of these vertical lines indicate model parameters with PRCC values significantly different from zero at a pvalue of 0.01. The critical values for the nonzero PRCC with p-value of 0.01 was determined using Eq. (J-1) (see also Eq. 7 in Marino et al., 2008 and p. 143 in Anderson, 2003) with 982 and 967 degrees of freedom; however, for our computation, we used a conservative value of 900 as the degrees of freedom to generate these cut-offs. Also, we set N = 1000, the number of simulations. The key parameters influencing the number of merozoites can be separated into two groups - those with significantly negative PRCC values (that thus decrease the number of merozoites when increased) and those with significantly positive PRCC values (that increase the number of merozoites when increased). Key parameters from both models with negative PRCC values include the immunogenicity of merozoites (ρM) and infected red blood cells (ρY), the infection rate (β), and the production of RBC from the bone marrow (πX). Additionally, we see a significant effect of the immunogenicity of gametocytes (ρG) for the coupled model only. Only one parameter show significantly positive PRCC values - the immunosensitivity of the gametocytes (kG). It is interesting to note, from Figs. 5(a) and 5(b), that the magnitude of the PRCC of the within-human host model (2.2) is amplified when the model is coupled with the population-level model (2.4)-(2.8) (hence leading to the coupled model (2.10)). The sensitivity analysis was also carried out using as the response function the total number of gametocytes (G) for the within-human host model (2.2) and coupled model (2.10) and the results obtained are given in Fig. 5(c) and (d). The PRCC cutoff point in this case are set to ± 0.0821. and ± 0.0827 respectively (again, indicated with vertical lines in the figures). Parameters with significantly negative PRCC values for both models are again the immunogenicity of merozoites (ρM) and infected red blood cells (ρY), the immunosensitivity of the gametocytes

(kG), the RBC infection rate (β), and the production rate of RBC from the bone marrow (πX). We once again see a significant negative effect of the immunogenicity of gametocytes (ρG) for the coupled model only, as we did when we used the number of merozoites as the response function, but now we also see a significant negative effect of the death rate of infected RBC (μY) for the coupled model. Both the immunosensitivity of merozoites (kM) and of infected red blood cells (kY) positively affect both models, while the death rate of immune cells (μB) shows a significantly positive PRCC for the coupled model only. Thus, this study identifies the most important parameters that drive the transmission mechanism of the disease within the human host. The identification of these key parameters is vital to the formulation of effective control strategies for combating the spread of the disease. In other words, the results of this sensitivity analysis suggest that a strategy that reduces the infection rate (reduce β), increase the production of RBC (increase πX), the immunogenicity of the infected red blood cells (increase ρY), the merozoites (increase ρM), the gametocytes (increase ρG), and the immunosensitivity of the infected RBCs (increase kY), the merozoites (increase kM) and the gametocytes (increase kG), decrease the death rate of immune cells (reduce μB), increase the death rate of infected RBCs (increase μY) will be effective in curtailing the spread of malaria within an infected human in the community. Effective control of the parasite in an infected individual can be achieved by administration of highly sensitive drugs or vaccine. However, the mechanism for sure control is beyond the scope of this current study. Next, we determine the PRCC values for each parameter of the population level model (2.4)-(2.8) and coupled model (2.10), using as response function the total number of infected humans (IH); the parameter ranges and baseline values are tabulated in Table 2 (see Fig. 6(a) and (b)). The PRCC cutoff points are set at ± 0.0818 and ± 0.0827 respectively for the two models, corresponding to a p-value of 0.01 (as in Marino et al. (2008); cut-off values are indicated with vertical lines in Fig. 6(a) and (b)). The significance of the non-zero PRCC of the population level model (2.4)-(2.8) with p-value of 0.01 was determined using Eq. (J-1) with 988 degrees of freedom and N = 1000, the number of runs (simulations). There is considerably less overlap in the parameters with 41

Ecological Complexity 38 (2019) 31–55

F.B. Agusto, et al.

Fig. 6. PRCC values of the population-level model (2.4)-(2.8) (panels (a) and (c)) and coupled model (2.10) (panels (b) and (d)). (a) Total number of infected humans using populationlevel model; (b) Total number of infected humans using the coupled model; (c) Total number of infected mosquitoes using population-level model; (d) Total number of infected mosquitoes using the coupled model. Parameter values and ranges used are as given in Table 2. Vertical lines give PRCC values significantly different from zero with p-value = 0.01; negative values smaller than these values and positive values greater than these values are significant.

significantly negative or positive PRCC values when comparing the population-level model and the coupled model. From Fig. 6(a) and (b), we note that the key parameters from both models with negative PRCC values that most influence the disease transmission dynamics are the death rate of mosquitoes (μV) and the infectious human recovery rate (γH). In addition, for the population-level model, the human recruitment rate (πH) also has a significantly negative PRCC value. In contrast, there are seven parameters that uniquely have negative effects when considering the coupled model. These are the malaria induced death rate (δh), the sporozoites clearance rate (μS), the immunogenicity of the gametocytes (ρG), the merozoites (ρM), and the infected red blood cells (ρY), the immunosensitivity of the gametocytes (kG), and the red blood cell infection rate (β). Three parameters have a significantly positive affect on the total number of infected humans for both models. These are the scaling factor of transmission probability per contact in mosquitoes (βm) and in humans (βh), and the mosquito biting rate (bM). In addition, for the coupled model, we observe an additional 5 parameters with significantly positive PRCC; these are the mosquito and human recruitment rates (πM and πH), the sporozoite production rate (πS), the immunosensitivity of the merozoites (kM), and the immune cells death rate (μB). We observed from Fig. 6(a) and (b) that for the populationlevel model, the magnitude of the PRCC is reduced when the models are coupled, unlike in the within-human host model (2.2). The sensitivity analysis was also implemented using the total number of infected mosquitoes (IM) as the response for the populationlevel model (2.4)-(2.8) and coupled model (2.10). As with the previous analyzes, we set the PRCC cutoff point for the population-level and coupled models at ± 0.0818 and ± 0.0827 respectively. Use of this response parameter caused the two models to differ greatly in which parameters showed significant PRCC (Fig. 6(c) and (d)). Only one parameter with a significantly negative PRCC was shared by both models, the death rate of mosquitoes (μV). There were two additional parameters with negative PRCC for the population-level model, the human recovery rate (γH) and the recruitment rate in humans (πH). In contrast, there were 7 additional parameters with negative PRCC values for the coupled model. These were the immunogenicity of the

gametocytes (ρG), the merozoites (ρM), and the infected red blood cells (ρY), the immunosensitivity of the gametocytes (kG), the death rate of infected RBCs (μY), the RBC infection rate (β), and the production rate of RBC from the bone marrow (πX). There were no parameters in common for the two models with positive PRCC. For the populationlevel model, parameters with significantly positive effects on the total number of infected mosquitoes included the scaling factor of transmission probability per contact in humans (βh) and in mosquitoes (βm), the mosquito recruitment rate (πM) and the mosquito biting rate (bM). For the coupled model, there were three parameters with significantly positive PRCC values; these were the rate of host loss of malaria immunity (ωH), the immunosensitivity of the merozoites (kM), and the death rate of immune cells (μB). Similarly, we observed from Fig. 6(c) and (d), a significant reduction in the magnitude of the PRCC values in the coupled model compared with the population-level model (2.4)(2.8). Hence, from sensitivity analysis described above, we identified the most important parameters that drive the disease transmission at the population level. Based on our models, the results of these analyzes suggest that a strategy that increases human recovery rate (increase γH), the death rate of the mosquitoes (increase μV), and human recruitment rate (increase πH), decreases the mosquito biting rate (reduce bM), recruitment rate (decrease πM), the scaling factor of transmission probability per contact in mosquitoes (decrease βm) and in humans decrease βh) and decrease the sporozoites production and clearance rates (decrease πS and μS) will be effective in curtailing the spread of malaria in the community. These results further suggest that a strategy that increases the production rate of the red blood cells from the bone marrow (increase πX), decreases red blood cells infection (decrease β), the merozoites (increase kM), the gametocytes (increase kG) and the immunogenecity of the infected red blood cells (increase ρY), the merozoites (increase ρM) and the gametocytes (increase ρG) will be effective in curtailing the spread of malaria within an infected human in the community. In summary, the sensitivity analysis of the within-human host model (2.2), the population level model (2.4)-(2.8) and the coupled model 42

Ecological Complexity 38 (2019) 31–55

F.B. Agusto, et al.

(2.10) suggests the following:

host dynamics. For proper comparison, we consider the cases without immune response (B (t ) = 0 ) and with immune response B(t) ≠ 0. The values used for most parameters are the same as the baseline values in Table 2 except X = 5, = 0.09, c = 0.017, u = 0.0013, µY = 1/200, µB = 1/5, kY = 0.00015, kG = 0.00015, kM = 0.0001, Y = 0.001, G = 0.001, M = 0.001, and the initial conditions are taken as X (0) = 300, Y (0) = 20, M (0) = 30, G (0) = 15, B (0) = 15, SP (0) = 500, SH (0) = 100, IH (0) = 20, RH (0) = 20, SM (0) = 405, IM (0) = 20 . It should be pointed out that the values in these simulations are only of theoretical sense to illustrate the results in this paper. Thus, in Fig. 7(a), we observed an increase in the merozoites population in the absence of the immune response as we increased θM(bM) (i.e., we increased the biting rate (bM)). The increase in the merozoite population is due to the feed back from the infected mosquito population which increases as the biting rate increases (see merozoite part of model (2.10)). On the other hand, in the presence of immune response (see Fig. 7(b)) we observed an oscillation in the merozoites levels due to the actions of immune response which proliferate with increasing merozoites level. However, as the scaling factor θM(bM) increases, the increase in the merozoites level dissipates and the oscillation dampens. This damping is due to the increase in the feedback scaling factor θM(bM) (which increases with increasing biting rate bM). In Fig. 7(c) we observed a negligible change in the gametocytes level in the absence of immune response with the changing biting rate bM. However, with increasing proliferation of the immune response, we see a decrease in the gametocyte level as the scaling factor θM(bM) increases with increasing biting rate bM (see Fig. 7(d)). Notice that more immune cells are activated as the merozoites population increase with increased biting rating (see immune response part of model (2.10)), leading to elimination of more gametocytes population. On the other hand fewer gametocytes are eliminated at low biting rate when there are fewer merozoites present. Also observed in Fig. 7(d) is the damping oscillation in the gametocyte level with increasing biting rate bM. Next, we consider the impact of varying the scaling factor θM(bM) (which is defined as a linear function of the biting rate (bM)) on the population level dynamics. We observed in Fig. 8(a) a considerable reduction in the number of infected humans owing to the number of death induced by the disease. Recall that the diseased-induced mortality is a function of the merozoites (i.e., H (M (t )) = rH M (t ) ); hence the merozoites population (M(t)) increases as the scaling factor (θM(bM)) increases, leading to an overall increase in mortality of the infected human population. In Fig. 8(b), we observed similar decrease in the presence of immune response and in addition, more infected humans at low biting rate. This is due to the fact that at low biting rate there are fewer merozoites present as a result of the actions of the human host immune response which then lead to fewer death in the infected human population. It should be noted that the oscillation observed in the merozoites from the within-human host dynamics feed back into the population level dynamics and also the damping effect which occurs with increasing θM(bM). For the infected mosquitoes in the absence of immune response (see Fig. 8(c)), we observed an increase in the number of infected mosquitoes as the biting rate increases; the biting rate serves as a modification parameter that amplify the number infected (see the mosquito part of model (2.10)). In the presence of immune response we also observed an increase in the number of infected mosquitoes as the biting rate increases. In addition, we observe that the initial amplitude in the oscillations is greatest with the largest biting rate, while further oscillations respond in a complex way to the biting rate (see Fig. 8(d)). After the initial high amplitude, further oscillations are greatly dampened for the largest biting rate considered (bM = 0.9 ); however, these initial high amplitude is lesser for intermediate biting rates (that is, bM = 0.3 and bM = 0.7 ). Lastly, we consider the impact on the within-vector dynamics, of the scaling factor (θH(bM)) obtained from the feedback function H(IH) from the population level dynamics into the within-mosquito host dynamics

(i) The key parameters that influences the total number of merozoites using the within-human host model (2.2) and the coupled model (2.10) (see Fig. 5(a) and (b)) are the production rate of the red blood cells from the bone marrow (πX), the infection rate of the red blood cells (β), the immunogenecity of the infected red blood cells (ρY), the merozoites (ρM) and the gametocytes (ρG), and the gametocytes (kG). (ii) The dominant parameters impacting the total number of gametocytes for the within-human host model (2.2) and coupled model (2.10) (see Fig. 5(c) and (d)) are the production rate of the red blood cells (πX), the infection rate of the red blood cells (β), the death rate of immune cells (μB), the death rate of infected RBCs (μY), the immunosensitivity of the infected red blood cells (kY), the merozoites (kM) and the gametocytes (kG), and the immunogenicity of the infected red blood cells (ρY), the merozoites (ρM) and the gametocytes (ρG). (iii) The key parameters influencing the total number of infected humans using the population level model (2.4)-(2.8) and coupled model (2.10) (see Fig. 6(a) and (b)) are the infectious human recovery rate (γH), death rate of the mosquitoes (μV), human recruitment rate (πH), the scaling factor of transmission probability per contact in humans (βh) and mosquitoes (βm), mosquito biting rate (bM), mosquitoes recruitment rate (πM), malaria induced death rate (δh) and the sporozoites production and clearance rates (πS and μS). The dominant parameters also include the immunogenicity of the infected red blood cells (ρY), the merozoites (ρM) and the gametocytes (ρG), the immunosensitivity of the merozoites (kM) and gametocytes (kG), and the red blood cells infection rate (β), and the immune cells death rate (μB). (iv) The dominant parameters impacting the total number of infected mosquitoes using the population-level model (2.4)-(2.8) and coupled model (2.10) (see Fig. 6(c) and (d)) are the human recovery rate (γH), the host loss of malaria immunity (ωH), death rate of the mosquitoes (μV), mosquito biting rate (bM), scaling factor of transmission probability per contact in humans (βh) and mosquitoes (βm), recruitment rates in human (πH) and mosquitoes (πM). The dominant parameters also include the immunosensitivity of the merozoites (kM) and gametocytes (kG), the immunogenecity of the infected red blood cells (ρY), merozoites (ρM) and the gametocytes (ρG), the red blood cells production and infection rates (πX and β), the death rate of immune cells (μB), and the death rate of infected RBCs (μY). (v) The magnitude of the PRCC of the within-human host model (2.2) is amplified when the model coupled with the population-level model (2.4)-(2.8). However, this is not the case for the populationlevel model, as the magnitude of the PRCC is reduced due to the coupling. 5. Numerical simulations In this section we carry out a numerical exploration of the coupled model (2.10). We start by defining the scaling factor θM of the feedback function L(IM) from the population level dynamics into the withinhuman host dynamics via the merozoite sub-population as a linear function of the mosquito biting rate bM. Hence,

L (IM ) =

M (bM ) IM

=

m bM IM .

Similarly, we define the scaling factor θH of the feedback function H(IH) from the population level dynamics into the within-vector host dynamics as a linear function of the mosquito biting rate bM. Thus,

H (IH ) =

H (bM ) IH

=

h bM IH .

We then vary the biting rate bM and observe its impact on the dynamics of the coupled system by first noting its impact on the within-human 43

Ecological Complexity 38 (2019) 31–55

F.B. Agusto, et al.

Fig. 7. Simulation of the coupled model (2.10) as a function of time while varying θM(bM). (a) Merozoites without immune response; (b) Merozoites with immune response. (c) Gametocytes without immune response; (d) Gametocytes with immune response.

Fig. 8. Simulation of the coupled model (2.10) as a function of time while varying θM(bM). (a) Infected humans in the absence of immune response; (b) Infected humans in the presence of immune response. (c) Infected mosquitoes in the absence of immune response; (d) Infected mosquitoes in the presence of immune response.

44

Ecological Complexity 38 (2019) 31–55

F.B. Agusto, et al.

Fig. 9. Simulation of the coupled model (2.10) as a function of time while varying θH(bM). (a) Sporozoites in the absence of immune response; (b) Sporozoites in the presence of immune response.

via the infected humans. We note that as θH(bM) increases there is an observed initial increase in the sporozoites level Sp(t) in the absence of the human immune response (see Fig. 9(a)). However, in Fig. 9(b) we observed oscillations which amplify with increasing values of the scaling factor θH(bM). Note, that the oscillatory effect of the human immune response from the within-human host dynamics feedbacks on the within-mosquito dynamics. We summarized below some of the findings of the numerical exploration of the coupled model (2.10). First, the scaling factors θM and θH of the feedback functions L(IM) and H(IH) were defined as a linear function of the mosquito biting rate bM;

6. Discussion and conclusion In this paper, we developed a novel deterministic model linking the dynamics of a within-host model to population-level dynamics for malaria transmission. The coupling and feedback for the model were created using the various life history stages of the malaria parasite (sporozoites, merozoites, and gametocytes) both in the human host and the mosquito vector; consideration of both levels of dynamics in a twohost system is a novel feature of this model. At the population level, we assume that individuals in the community are homogeneous both in their mixing and in their response to the infection. Relaxing these assumptions will be important in considering important aspects of the transient behavior of the within-host systems. Our future work expanding the models developed here will explore a time-since-infection model that allows the assumption of within host homogeneity to be dropped. The study shows that the parasite-free equilibrium of the withinhuman host model (2.2) is locally asymptotically stable (LAS) when the associated reproduction number is less than one and unstable otherwise. The study further shows that the parasite-present equilibrium can be LAS or unstable depending the parameter set used. The disease-free equilibrium of the population level model (2.4)-(2.8) is locally asymptotically stable when the associated reproduction number is less than one, and the model exhibits backward bifurcation, where the stable disease-free equilibrium coexists with a stable endemic equilibrium. The coupled model (2.10) has a locally asymptotically stable infection-free equilibrium (IFE). Note that the presence of backward bifurcation in the population level model complicates the prospect of malaria eradication, particularly in malaria endemic regions due to malaria induce mortality (Agusto et al., 2013; Agusto and Lenhart, 2013; Chitnis et al., 2006; Niger and Gumel, 2008). Since the requirement that the reproduction number be less than one no longer guarantees disease elimination, the reproductive number needs to be lower than a threshold quantity to ensure that malaria is eliminated in the community; this threshold quantity, in turn, is less than one (Agusto et al., 2013; Agusto and Lenhart, 2013; Chitnis et al., 2006; Niger and Gumel, 2008). Fig. 2 depict the simulation of the within-host model (2.2). The figure shows an oscillatory dynamics due to the response of the individual’s immune system to the infection and the development of merozoites and gametocytes in the infected human. Fig. 2 shows that at every peak cycle (corresponding to the peak level in both merozoites and gametocytes) there is a sharp increase in the immune cells level. These peaks and time to peak match the lowest point of the healthy red

(a) A comparison of the numerical simulations considered the cases without immune response (B (t ) = 0 ) and with immune response B (t) ≠ 0; (b) At the within-human host (see Fig. 7), we observed that increasing θM(bM) (i.e., increasing the biting rate bM) leads to the following: (i) An increase in the merozoites population level in the absence of immune response; (ii) Oscillations in the merozoite levels in the presence of immune response; (iii) The observed oscillations in the merozoites level dampens with increasing biting rate; (iv) A reduction in the gametocyte level as the proliferation of the immune response increases due to increasing biting rate bM; (c) At the population level dynamics (see Fig. 8), increasing θM(bM) leads to the following: (i) A reduction in the number of infected humans in the absence of immune response due to an increase in the cumulative mortality; (ii) The oscillation and damping effect observed in the merozoites from the within-human host dynamics feed back into the population level dynamics; (iii) In the presence of immune response a reduction in the number of infected mosquitoes is observed with negligible oscillation; (d) At the within-vector dynamics (see Fig. 9), we observed that increasing θH(bM) leads to the following: (i) The sporozoites level increases as θH(bM) increases; (ii) The oscillatory effect of the human immune response from the within-human host dynamics feed back into the within-mosquito dynamics amplifying the oscillations in the sporozoites level as scaling factor θH(bM) increases.

45

Ecological Complexity 38 (2019) 31–55

F.B. Agusto, et al.

blood cells (see Fig. 2(e)). Note that the time unit in this study is in days (the units of all our model parameters, see Table 2) and not hours. As part of our future work, we will consider ways to relate our results here with the much faster oscillations in malarial symptoms observed in malaria patients due to the synchronous development of the malarial parasite within the host (Baron, 1996; Hoffman et al., 2012; Maier et al., 2009; Wiser, 2011). In other words, we see oscillations due to our inclusion of an immune response in malarial patients; we know there are also oscillations due to timing of development of the parasite, which is not included in our model. Using a global sensitivity analysis, this study identifies (for the given models) the most important parameters driving the disease transmission dynamics both at the within-human host level and in the coupled model by considering the model outcomes such as the total number of merozoites and the total number of gametocytes. The identification of these parameters is vital to the formulation of effective control strategies for combating the spread of the disease. In other words, the results of this sensitivity analysis suggest (under the assumptions of the models formulated in this study) that strategies that increase the production of RBC, reduce the infection rate, increase the immunogenecity of the infected red blood cells, merozoites and the gametocytes, increase the immunosensitivity of the merozoites, reduce the immunosensitivity of the gametocytes, decrease the death rate of immune cells, and increase the death rate of infected RBCs will be effective in curtailing the spread of malaria within an infected human host. This study also identified the most important parameters at the population level and also in the coupled model using model outcomes such as the total number of infected humans and the total number of infected mosquitoes. The results of these analyzes suggest that strategies that increase the human recovery rate, the death rate of the mosquitoes, and the recruitment rate of uninfected susceptible humans into the population, will be effective in curtailing the spread of malaria between humans at the population level. Increasing the production rate of the red blood cells from the bone marrow, the immunosensitivity of the infected red blood cells (RBC), the merozoites, and the gametocytes will all likewise decrease the spread of malaria at the between human host level. The same effect is seen by decreasing red blood cells infection and increasing the immunogenicity of infected RBC, merozoites and the gametocytes. Finally, decreasing the mosquito biting rate, the mosquito recruitment rate, and the scaling factor of transmission probability per contact in mosquitoes and in humans will decrease between host transmission. The scaling factor of transmission probability per contact in mosquitoes is in effect the effectiveness of infected human transmitting the parasite (gametocytes) to a susceptible mosquito; since within the mosquito’s midgut, the parasite has to overcome the actions of both the human immune response still acting in the blood meal (Bousema et al., 2013; Da et al., 2015; Gouagna et al., 2004), and the vectorâs innate immune response (Blandin and Levashina, 2004; Da et al., 2015; Yassine and Osta, 2010) and the low efficiency of gamete fertilization (Alavi et al., 2003; Da et al., 2015; Gouagna et al., 1998; Vaughan, 2007). The scaling factor of transmission probability per contact in humans measures the effectiveness of an infected mosquito in transmitting the malaria parasite (sporozoites) in it’s salivary glands to a susceptible human (Baton and Ranford-Cartwright, 2005; Beier, 1998; Centers for Disease Control and Prevention, 2015; Kappe et al., 2003; Sreenivasamurthy et al., 2013; Wiser, 2011). Unlike in the case of human, the mosquito innate immune response does not act on the parasite at this stage; thus, the parasite has a higher probability of being transmitted to a susceptible human. Thus, a strategy that decreases the sporozoites production and clearance rates will be effective in curtailing the spread of malaria in the community. A look at the result of the global sensitivity analysis for the coupled system with the different response functions indicate some common parameters that are vital to the disease burden. We have indicated above the implication or possible impact of different control measures on these parameters, without reference to which of these parameters

can actually be used for control. Furthermore, it is vital to know how much control efforts will be required both at the in-human host level and at the population level to curtail the disease and transmission in the community effectively. That is, should more effort be at controlling the disease at the population level, as is currently is, with the use of insecticide treated bednet? Or, our effort should be geared towards the use of effective treatment and prevention and how much of this is necessary to guide against the development of drug resistance malaria? Or what global control strategy targeting the two levels is necessary for effective disease control? In future work, we will aim to address these vital questions in the context of our coupled system. The model (2.10) requires the coupling of the dynamics at the hosts (humans and mosquitoes) level into the population level, and we assumed the crucial fast-slow framework in the theoretical analysis. This assumption requires that the death rate of the infected red blood cells are greater than the recovery rate in humans (i.e., μY > γH) and that the sporozoites clearance rate in the mosquito salivary glands be greater than the mosquitoes death rate (i.e., μS > μV). These assumptions are biologically feasible since infection turnover in humans occurs at a faster rate than the turnover at the population level. Similarly, the sporozoites turnover in mosquitoes occurs at a faster rate than the turnover at the population level, thus indicating the fast time scale within the two hosts compared to the much slower time scale at the population level. However, model (2.10) also includes a feedback function from the population level into the respective within host, the fast-slow assumption here requires that the death rate of the infected red blood cells be greater than the natural death rate in mosquitoes (i.e., μY > μV) for the feedback from infected mosquitoes into withinhuman dynamics and the sporozoites clearance rate in the mosquito salivary glands be greater than the recovery rate in humans (i.e., μS > γH). The novel feature of the coupled model in this study is the presence of feedback functions linking the dynamics of the epidemic process at the population level to the within-host (human and mosquito) level. Other studies have considered nested models where the interaction between the selective levels depends, for example, on the specific function describing the dependence of the host immune response to the within-host pathogen replication rate (Gilchrist and Sasaki, 2002) or on how shedding of pathogens affects both within-host pathogen levels and between-host transmission (Barfield et al., 2015). In our study, the functions linking the two dynamic levels are assumed to be linear functions of the mosquito biting rate. In future studies of this coupled model, we will consider more complex and perhaps more realistic functions linking the dynamic levels, such as a saturating feedback function. The findings of the numerical analyzes show that as the biting rate increases (in other words, as the feedback function increases), the merozoite population increases in the absence of a human host immune response. The presence of an immune response leads to the development of oscillations in the merozoite population since the immune cells add an internal feedback to the within-host system. However, these observed oscillations in the merozoite population dampen as the biting rate increases. This may be due to a similar effect as seen in Orive et al. (2005), where oscillations in the numbers of cells and virus within the host were seen to dampen with increasing viral migration between compartments within a structured host. Here, the biting rate acts as a type of âǣmigrationâǥ rate between hosts, and the multiple hosts are taking the place of the multiple within-host compartments. At the population level, increasing the biting rate leads to the same kind of oscillation and damping effect observed in the merozoites, as the withinhuman host dynamics feedback into the population level dynamics. At the within-vector dynamics level, we see that increasing the feedback function leads to an increase in the sporozoite level. Increasing the feedback function also results in the same kind of oscillatory effect earlier observed, again due to a feedback between the within-human host dynamics and the within-vector dynamics. However, here the oscillations in the number of sporozoites within the vector are amplified 46

Ecological Complexity 38 (2019) 31–55

F.B. Agusto, et al.

as we increase the feedback function, whereas the oscillations in the number of merozoites within the human host are dampened with increases in the feedback function. This amplifying oscillatory effects may be due to the fact that the within-vector model does not include any of the mosquito vector innate immune response (Blair, 2011; Blandin and Levashina, 2004; Richman et al., 1997; Yassine and Osta, 2010), and so represents an overly simplified model of the actual within-vector dynamics. In future work, we will consider in more detail the withinvector dynamics, adding an innate immune response to the mosquito vector.

Acknowledgments FBA began this work as the 2015 Langston Hughes visiting Professor at the Department of Ecology and Evolutionary Biology, University of Kansas. FBA acknowledges the useful discussion with Eric Numfor during the model formulation phase. FBA, acknowledges the supported of the Strategic Environmental Research and Development Program under grant RC-2639. MEO acknowledges support from NSF grant DEB1354754. The authors acknowledges with thanks the several reviewers for their comments and suggestions that have help improve this final version.

Appendix A. Proof of Lemma 2 Lemma 2: Let the initial data X(0) > 0, Y(0) ≥ 0, G(0) ≥ 0, M(0) ≥ 0 and B(0) ≥ 0. Then the solutions (X, Y, G, M, B) are all non-negative for all t > 0. Proof. Let 1= sup {τ > 0: X(0) > 0, Y(0) > 0, G(0) > 0, M(0) > 0, B(0) > 0 ∈ [0, τ]}. Thus τ1 > 0. It follows from the first equation of model (2.2) that

X (t ) =

X (t ) M (t )

X

µX X (t )

which can be written as:

d d

1

X (t ) exp

0

( M (s ) + µX ) d

=

1

X exp

0

( M (s ) + µX ) d .

Hence, 1

X ( 1) exp

0

( M ( s ) + µX ) d 1

1

X (0) =

p

X exp

0

0

( M (s ) + µX ) dp dp

so that 1

X ( 1) = X (0) exp

0

1

×

X

0

exp

( M (s ) + µX ) d 1 + exp p 0

1

0

( M ( s ) + µX ) d 1

( M (s ) + µX ) dp dp > 0.

It can similarly be shown that Y(t) > 0, G(t) > 0, M(t) > 0 and B(t) > 0 for all τ > 0. □ Appendix B. Proof of Lemma 3 Lemma 3: The region

WH

5 +

=

is positively-invariant for the model (2.2) with non-negative initial conditions in

Proof. Let (X (t ), Y (t ), M (t ), G (t ), B (t )) (2.2), gives

NWH (t )

X

5 +,

5 +.

be the solution of model (2.2) with non-negative initial conditions. Adding the first 2 equations of

µmin NWH (t ).

Thus, X

lim sup NWH (t )

µmin

.

From the next two equations of model (2.2)

lim sup M (t )

rµY X µM µmin

and

lim sup G (t )

c X . µG µmin

Similarly, the equations for the immune cells yields,

lim sup B (t ) <

,

and it follows that,

NWH (t )

0 if

NWH (0)

X

µmin G (t )

, M (t ) 0 if

0 if M (t )

G (t )

r X , µM µmin

c X , B (t ) < 0 if B (t ) > B¯ . µG µmin

Thus the vector fields corresponding to the equations in model (2.2) points inward into Φ. Hence all feasible solutions of the RBCs, IRBCs, the parasites of model (2.2) and all feasible solutions of the immune cells enter the region Φ (Niger and Gumel, 2011). Hence, Φ is positively invariant. 47

Ecological Complexity 38 (2019) 31–55

F.B. Agusto, et al.

Therefore, all possible solutions of model (2.2) enter the region

WH

5 +,

=

as τ → ∞. Hence, ΦWH is positively invariant. □

Appendix C. Proof of Lemma 5 Lemma 5: The PPE

is locally asymptotically stable (LAS) if ˜ WH | B =0 > 1.

1 WH

1 Proof. The local stability of the parasite-present equilibrium (EWH ) in the absence of immune response and feedback L (IH ) = 0, is established using the eigenvalues of the Jacobian matrix of system (2.2) evaluated at the equilibrium 1WH

J

M ** µX 0 M ** k1 u M ** rµY 0 c 0 0

=

1 WH

µM

X ** X ** u X ** 0 0

0 0 0 µG 0

where k1 = µY + c . The eigenvalues of the matrix J given as

J 11

M ** µX 0 M ** k1 u M ** rµY

=

WH

µM

0 kY Y ** kM M ** kG G** Y Y ** + M M ** + G G ** 1 WH

, µB

are determined from the characteristics polynomials of matrices J 11

WH

and J 41 . Matrix J 11 WH

WH

is

X ** X ** u X **

The characteristics polynomial is expressed as 3

+ b1

2

+ b 2 + b3 ,

where 2 ˜ WH µX2 + ˜ WH k1 µX + ˜ WH µM µX + u µ ˜ WH

b1 =

X

X

2 ˜ WH µX (k1 + µM ) + u ˜ WH

b2 =

b3 = k1 µM µX ( ˜ WH

X

1).

The coefficients b1, b2, b3 > 0 provided ˜ WH > 1 and b1b2 > b3, thus it follows from Routh-Hurwitz criteria that the matrix J with negative real parts. The Matrix J 41 is given as

1 WH

have eigenvalues

WH

J 41

WH

µG 0

=

Y Y **

+

kG G** + G G **

M M **

The eigenvalues of matrix J 41

WH

˜

Y µX µM ( WH

µG , and

(rµY

Thus matrix J 41

WH

˜

Y µX µM ( WH

(rµY

are given as

+

M µX

( ˜ WH

1)

+

c

˜

G µX µM ( WH

(rµY

1) k1 u ) µG

µB .

have eigenvalues with negative real parts provided

1)

k1 u)

1)

k1 u )

µB

M µX

+

( ˜ WH

1)

+

c

˜

G µX µM ( WH

1) < µB . k1 u ) µG

(rµY

1 Thus, it follows that the parasite-present equilibrium (EWH ) in the absence of immune response is local stability whenever ˜ WH > 1 and ˜ Y µX µM ( WH (rµY

k1 u)

1)

+

˜ M µX ( WH

1)

+

c G µX µM ( ˜ WH 1) (rµY k1 u) µG

Appendix D. Existence of the PPE

< µB . □

2 WH

The existence of the parasite-present equilibrium of Eq. (2.2) to zero, we obtain the following

0 =

X

X (t ) M (t )

0 =

X (t ) M (t )

µY Y (t )

0 = B (t )[

µG G ( t ) YY

(t ) +

in the absence of feedback L (IM ) = 0, is established as follows. Setting the right-hand side

µX X (t ),

0 = L (IM ) + rµY Y (t ) 0 = cY (t )

2 WH ,

cY (t )

µM M (t )

kY B (t ) Y (t ), kM B (t ) M (t )

u X (t ) M (t ),

kG B (t ) G (t ), M M (t )

+

G G (t )]

(D-1)

µB B (t ).

From the fifth equation of (D-1), with B ≠ 0, we have Y ** = ⁎⁎

(µB

G G ** Y

48

M M **)

. Substituting Y

⁎⁎

into the remaining equations of (D-1) gives:

Ecological Complexity 38 (2019) 31–55

F.B. Agusto, et al.

X **M **

0 =

X

0 =

X **M ** +

µX X **, G G **

k1 (

+

M M **

µB )

+

kY B ** (

G G **

Y G G **

rµY (

0 =

M M **

+

c(

G G **

+

M M **

M M **

µB )

,

Y

µB )

µM M **

Y

0 =

+

µB )

µG G **

Y

u X **M **

kM B **M **,

kG B **G **,

where k1 = µY + c . This give the equilibrium

M ** =

b2 + 4ac , 2a

b+

X ** =

Y ** =

G G **

(µB

M M **)

,

G** =

Y

rµ Y µB (µG + kG B **) (µG + kG B **)*(µM + kM B ** + u X **) Y + rµY (µG + kG B **)

M

c

G

c (µB + µG

Y

M M **) + kG B **

, Y

+ c *(µM + kM B ** + u X **)

, G

where

a = µX u (kG B ** Y + c G + µG Y ), b = c (µX kM B ** + µX µM X u ) G + µX rµY (µG + k G B **) M + (µX kM B ** + µX µM X u )(µG + kG B **) Y + r µY µB (µG + kG B **), c = X (µM + kM B **)(µG + kG B **) Y + X rµY (µG + kG B **) M + X c (µM + kM B **)

G.

Substituting M into the expression for G gives the equilibrium expression for G⁎⁎ in terms of B⁎⁎; also, substituting these into Y⁎⁎ and X⁎⁎ 2 respectively gives their respective expressions at equilibrium in terms of B⁎⁎. Thus, the equilibrium WH exist provided B⁎⁎ ≠ 0 and ** ** M + G < µ , from the fifth equation of (D-1). M G B ⁎⁎

⁎⁎

Appendix E. Proof of Lemma 6 Lemma 6:The PPE

2 WH

is locally asymptotically stable (LAS) if the Jacobian matrix J

2 WH

have negative eigenvalues.

Proof. To examine the local stability of the parasite-present equilibrium ( in the presence of immune response, we follow the approach in Tumwiine et al. (2008) by using the following simplified form of system (2.2): 2 WH )

X **M ** X = M ** + µX , = k1 + kY B **, Y ** X ** ** rµY Y cY ** = µM + u X ** + kM B **, = µG + kG B **, M ** G ** 0 = Y Y ** + G G** + M M ** µB , with B ** 0, and k1 = µ Y + c. Hence, the Jacobian matrix of (2.2), evaluate at X

0

X **

XM ** Y **

M ** J

2 WH

=

u M ** 0

c

0

B **

0

X

0

kY Y **

0

kM M **

M **

0 B **

Y

M

The characteristics polynomial of matrix J 5

+ c1

4

+ c2

3

+ c3

2

+ c4

1

0

cY ** G **

B **

(E-1)

and using the simplified system (E-1) is given as

X **

rµY Y **

rµY

2 WH

G 2 WH

kG G ** 0 is given as

+ c5,

where the coefficients are expressed as

49

Ecological Complexity 38 (2019) 31–55

F.B. Agusto, et al.

rµ Y ** M **X ** cY ** + Y + , Y ** M ** G** Y Y ** + kM M M ** + G k G G **)

X + X ** = B ** (kY

c1 = c2

[(Y **) 2rµY + (M **)2 X **] c [(Y **) 2rµY X + (M **) 2 (X **) 2 + X Y **M **] , + M **X **Y ** G **M **X ** = B **Y **kY (c G + rµ Y M ) + M **X ** (rµY + Y B kM B **) B ** X (Y ** Y kY + M ** M kM + G** G kG ) X **M **B ** (M ** M kM + G ** G kG ) + + X ** Y ** +

c3

u 2M **X **

+ +

X

c [(Y **) 2B **

Y kY

+ M **Y **B ** G**

rµY Y ** [X **Y **G**B **

Y kY

M kM

+ M **

+ (G **) 2B **

X]

G kG X

+c

X Y **]

X **M **G** u 2X **M ** [c (Y **) 2 + X **M **G ** ] , Y **G **

+ M ** (X **) 2c G kM + X **M ** Y kM X ] X ** B ** X [c (Y **)3 Y kY + cM ** (Y **) 2 M kM + (M **) 2 M kM G **X ** + (G**)2M ** G kG X **] + X **Y **G ** Y **B **rµY [(G**)2 G kG X + Y **G** Y kY X + X **Y **G**c G kY + X ** (Y **)2c Y kY ] + X **M **G ** c (X **Y **M **r 2µ Y + X ** (M **) 2B ** M kM + X **Y **M **B ** Y kM + (Y **) 2B **r M µY kY ) + G** X **M ** 2 [(G**)2B **u G kG + X **M **cu + Y **G **B **u Y kY + M **G**B ** Y kM ] , G**

c4 =

B ** [Y **r

c5 = M **B ** +

G

M µY kY X

+ Y **c

(X **G**rkG µY

+

X kM c )

+

B (Y **)2crkY µ Y X ( G G ** + X **M **G**

Y Y **)

+ X **Y **M ** Y kM + (Y **)2r M µ Y kY ] X **G** X **M **B ** G 2 [(Y **)2cukY + M **ckM Y ** + X **M **G**ukG ] Y ** X **Y **M **B **c Y 2 (ukY Y ** + kM M **) . G ** B **c

X

[X ** (M **) 2

G kY X

M kM

If the coefficients c1, c2, c3, c4, c5 > 0 and c1 c2 c3 > c32 + c12 c4 and (c1 c4 c5)(c1 c2 c3 c32 c12 c4 ) > c5 (c1 c2 c3) 2c1 c52. Then, by Routh-Hurwitz criteria 2 2 the matrix J WH will have eigenvalues with negative real parts. Hence, the matrix J WH will be locally stable. □ Appendix F. Existence of the PPE

3 WH (IM )

The existence of the parasite-present equilibrium 3WH (IM ), in the presence of feedback L (IM ) = right-hand side of Eq. (2.2) to zero, we obtain the following

0 =

X

0 =

X (t ) M (t )

0 =

M IM

0 = cY (t ) 0 = B (t )[

X (t ) M (t )

µG G ( t ) YY

(t ) +

cY (t )

µM M (t )

kY B (t ) Y (t ), kM B (t ) M (t )

u X (t ) M (t ),

kG B (t ) G (t ), M M (t )

With B ≠ 0, we have Y ** = ⁎⁎

0 is established as follows. Setting the

µX X (t ),

µY Y (t )

+ rµ Y Y (t )

M IM

+ (µB

G G (t )] G G ** Y

(F-1)

µB B (t ). M M **)

, from the fifth equation of (F-1). Substituting Y

50

⁎⁎

into the remaining equations of (F-1) gives:

Ecological Complexity 38 (2019) 31–55

F.B. Agusto, et al.

X **M **

0 =

X

0 =

X **M **

µX X **, G G **

k1 (µB

M M **)

kY B ** (µB

G G **

Y

0 =

H IH

G G **

rµY (µB

+

M M **)

µM M **

Y

0 =

G G **

c (µB

M M **)

µG G**

Y

M M **)

,

Y

u X **M **

kM B **M **,

kG B **G**,

where k1 = µY + c . This leads to the equilibrium points

Y ** (IM ) = M ** (IM ) = G ** (IM ) =

b2 + 4ac , 2a

b+

X ** (IM ) =

(µG + kG B **)(µB kM B ** + µB µM + µB u X X ** M M IM ) (µG + kG B **)(µM + kM B ** + u X **) Y + rµY (µG + kG B **) M + c (µM + kM B ** + u X **) (µG + kG B **)

Y M IM

(µG + kG B **)(µM + kM B ** + u X **)

Y

+c

G M IM

, G

+ rµY µB (µG + kG B **)

+ rµY (µG + kG B **)

c (µB kM B ** + µB µM + µB u X X ** (µG + kG B **)(µM + kM B ** + u X **) Y + rµY (µG + kG B **)

M

+ c (µM + kM B ** + u X **)

, G

M M IM ) M

+ c (µM + kM B ** + u X **)

, G

where

a = µX u (kG B ** Y + c G + µG Y ), b = c (µX kM B ** + µX µM X u + X M IM ) G + µX rµY (µG + k G B **) M + (µX kM B ** + µX µM X u + X M IM )(µG + k G B **) Y + r µY µB (µG + k G B **), c = X (µM + kM B **)(µG + kG B **) Y + X rµY (µG + kG B **) M + X c (µM + kM B **)

G.

Substituting X into the expression for Y , M and G gives their respective expressions at the equilibrium 3WH (IM ) in terms of B⁎⁎ and the IM. Thus, from the fifth equation of (F-1), we have that the equilibrium 3WH (IM ) exist provided B⁎⁎ ≠ 0 and M M ** (IM ) + G G ** (IM ) < µB . ⁎⁎

⁎⁎

⁎⁎

⁎⁎

Appendix G. Proof of Lemma 3.1.4 Lemma 3.1.4: The PPE

3 WH (IM )

is locally asymptotically stable (LAS) provided B* ≠ 0 and

Proof. To establish the local stability of the parasite-present equilibrium ( Tumwiine et al. (2008) by using the following simplified form of system (2.2):

3 WH (IM ) )

Y Y ** (IM )

+

Hence, the Jacobian matrix of (2.2), evaluate at X

M ** (IM ) J

3 WH (IM )

=

X ** (IM )

0

X ** (IM )

u M ** (IM )

XM ** (IM ) Y ** (IM )

X ** (IM )

rµ Y

rµY Y ** (IM ) M ** (IM )

0

c

0

B **

B **

The characteristics polynomial of matrix J 5

+ d1

4

+ d2

3

+ d3

2

+ d4

1

0

0

kY Y ** (IM )

0

kM M ** (IM )

cY ** (IM ) G ** (IM )

0 Y

and using the simplified system (G-1) is given as

0

M

B **

3 WH (IM )

is given as

G

kG G** (IM ) 0

+ d5,

where the coefficients are expressed as

51

+

G G ** (IM )

< µB .

in the presence of immune response, we follow the approach in

X ** (IM ) M ** (IM ) X = M ** (IM ) + µX , = k1 + kY B **, X ** (IM ) Y ** (IM ) rµY Y ** (IM ) cY ** (IM ) = µM + u X ** (IM ) + kM B **, = µG + kG B **, M ** (IM ) G ** 0 = Y Y ** (IM ) + G G ** (IM ) + M M ** (IM ) µB , with B ** 0, and k1 = µY + c. 3 WH (IM )

M M ** (IM )

(G-1)

Ecological Complexity 38 (2019) 31–55

F.B. Agusto, et al.

rµ Y ** (IM ) M **X ** (IM ) cY ** (IM ) X , + + Y + X ** (IM ) Y ** (IM ) M ** (IM ) G ** (IM ) d2 = B ** (kY Y Y ** (IM ) + kM M M ** (IM ) + G kG G** (IM )) u 2M **X ** (IM ) d1 =

+ + d3 = + + + +

d4 =

[(Y ** (IM ))2rµY + (M ** (IM ))2 X ** (IM )] M **X ** (IM ) Y ** (IM ) c [(Y ** (IM )) 2rµ Y X + (M ** (IM )) 2 (X ** (IM ))2 + X Y ** (IM ) M ** (IM )] , G ** (IM ) M ** (IM ) X ** (IM ) B **Y ** (IM ) kY (c G + rµY M ) + M ** (IM ) X ** (IM ) (rµ Y + Y B kM B **) B ** X (Y ** (IM ) Y kY + M ** (IM ) M kM + G ** (IM ) G kG ) X ** (IM ) X ** (IM ) M ** (IM ) B ** (M ** (IM ) M kM + G** (IM ) G kG ) Y ** (IM ) c [(Y ** (IM )) 2B ** Y kY + M ** (IM ) Y ** (IM ) B ** M kM + M ** (IM ) X ] G** (IM ) rµY Y ** (IM )[X ** (IM ) Y ** (IM ) G** (IM ) B ** Y kY + (G** (IM )) 2B ** G kG X + c X ** (IM ) M ** (IM ) G** (IM ) u 2X ** (IM ) M ** (IM )[c (Y ** (IM )) 2 + X ** (IM ) M ** (IM ) G** (IM ) ] , Y ** (IM ) G** (IM ) [B ** [Y ** (IM ) r M µY kY X + Y ** (IM ) c G kY X + M ** (IM )(X ** (IM )) 2c G kM X

+ X ** (IM ) M ** (IM ) + (M ** (IM ))2

X ]]/[X ** (IM )]

kM

Y

+ [B **

X

[c (Y ** (IM ))3

2 M kM G **X ** (IM ) + (G ** (IM )) M ** (IM )

+ [Y ** (IM ) B **rµ Y [(G** (IM )) 2 + X ** (IM )(Y ** (IM

))2c

G kG X

+ Y ** (IM ) G ** (IM )

Y kY X

M kM

+ X ** (IM ) Y ** (IM ) G** (IM ) c

+ [c (X ** (IM ) Y ** (IM ) M ** (IM ) r

+ X ** (IM ) Y ** (IM ) M ** (IM ) B **

M kM

+ cM ** (IM )(Y ** (IM )) 2

G k G X ** (IM )]]/[X ** (IM ) Y ** (IM ) G ** (IM )]

Y kY ]]/[X ** (IM ) M ** (IM ) G ** (IM )]

+ X ** (IM )(M ** (IM ))2B **

Y kY

X Y ** (IM )]

Y

2µ

G kY

Y

kM

+ (Y ** (IM ))2B **r M µY kY )]/[G ** (IM )] [X ** (IM ) M ** (IM ) 2 [(G** (IM )) 2B **u G kG + X ** (IM ) M ** (IM ) cu + Y ** (IM ) G** (IM ) B **u Y kY + M ** (IM ) G** (IM ) B ** Y kM ]]/[G** (IM )], d5 = M ** (IM ) B ** (IM )

(X ** (IM ) G ** (IM ) rkG µY

G

B (Y ** (IM ))2crkY µ Y X ( G G ** (IM ) + + X ** (IM ) M ** (IM ) G ** (IM ) +

X kM c )

+ X ** (IM ) Y ** (IM ) M ** (IM ) Y kM + (Y ** (IM )) 2r M µY kY ] X ** (IM ) G** (IM ) X ** (IM ) M ** (IM ) B ** G 2 [(Y ** (IM ))2cukY + M ** (IM ) ckM Y ** (IM ) + X ** (IM ) M ** (IM ) G** (IM ) ukG ] Y ** (IM ) X ** (IM ) Y ** (IM ) M ** (IM ) B **c Y 2 (ukY Y ** (IM ) + kM M ** (IM )) . G** (IM ) B **c

X

[X ** (IM )(M ** (IM ))2

+

Y Y ** (IM ))

M kM

The coefficients d1, d2, d3, d4, d5 > 0 and d1 d2 d3 > d32 + d12 d 4 and (d1 d 4 d5)(d1 d2 d3 d32 d12 d 4) > d5 (d1 d2 d3)2d1 d52 . Thus, by Routh-Hurwitz 2 2 criteria the matrix J WH have eigenvalues with negative real parts. Hence, the matrix J WH is locally stable. □ Appendix H. Proof of Lemma 9 Lemma 9: The region

BH

=

H

3 +

M

×

2 +

is positively-invariant for the model (2.4)-(2.8) with non-negative initial conditions in

5 +.

Proof. The following steps are followed to establish the positive invariance of ΦBH (i.e., solutions in ΦBH remain in ΦBH for all t > 0). The rate of change of the population is obtain by adding the equations of the model (2.4)-(2.8) and this gives

dNBH (t ) = dt

µH NBH (t )

H

H (M (t )) IH (t ).

(H-1)

and it follows that

dNBH (t ) dt

H

dNM (t ) = dt

M

µH NBH (t ), µV NM (t ).

The solutions are given as

NBH (t )

NBH (0) e

µH t

+

H

µH

(1

e

µH t

)

and NM (t ) = NM (0) e

µV t

+

M

µV

(1

thus, it follows that

52

e

µM t

),

Ecological Complexity 38 (2019) 31–55

F.B. Agusto, et al. H

limNBH (t )

M

and lim NM (t )

µH

t

µV

t

.

Thus, the region ΦBH is positively-invariant. Hence, it is sufficient to consider the dynamics of the flow generated by (2.4)-(2.8) in ΦBH. In this region, the model can be considered as being epidemiologically and mathematically well-posed (Hethcote, 2000). Thus, every solution of the model (2.4)(2.8) with initial conditions in ΦBH remains in ΦBH for all t > 0. Therefore, the ω-limit sets of the system (2.4)-(2.8) are contained in ΦBH. □ Appendix I. Proof of Lemma 11 Lemma 11: The infection-free equilibrium (IFE) of the coupled model (2.10), given by C

=

( X X * + µM g1)2 + 4rµY X X * X X * + µM + g1

J

0 C

=

X X*

0 0 0 µG 0 0 0 0 0 0 0

X X*

µM

u

X X*

0 0 0 0 0 0 0 0

0 0 0 0 µB 0 0 0 0 0 0

0 0 0 0 0 kp 0 0 0 0 0

0 0 0 0 0 0 µH 0 0, 0 0

where g1 = µY + c, kp = µS + S , k1 = H + µH , k2 = infection-free equilibrium (IFE), are given as: u X X * + µM + g1 2

is

locally C

=

+

(u X X * + µM

g1)2 + 4rµY X X *

2

asymptotically

is locally asymptotically stable (LAS) provided

< 1.

Proof. To establish the local stability of the infection-free equilibrium ( equilibrium (IFE), this is stated below as:

µX 0 0 g1 0 rµY 0 c 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 C,

stable

(u X X * + µM g1)2 + 4rµY X X * u X X * + µM + g1

,

provided

H

0 0 0 0 0 H

0 k1 H (M **) 0 0

H

0 k2 0 0

0 0 0 0 0 0 0 0 0 µV 0

we linearize the coupled model evaluated at the infection-free

0 0 M

0 0 0 0 0 0 0 µV

+ µH . The eigenvalues of Jacobian matrix J

u X X * + µM + g1 2

the

0 0 0 0 0 0

0 C ),

first

(u X X * + µM

g1)2 + 4rµY X X * 2 u X X * + µM + g1 2

eigenvalue

,

µX , +

µG ,

(u X X * + µM

0 C

µB ,

of the coupled model, evaluated at the

k p,

k1,

g1)2 + 4rµY X X * 2

k2, < 0.

µH , This

µV ,

µV . The IFE

holds

provided

< 1. □

Appendix J. Latin hypercube sampling and PRCC J1. Latin hypercube sampling (LHS) The Latin hypercube sampling (LHS) (Blower and Dowlatabadi, 1994; Marino et al., 2008; McKay et al., 2000; Sanchez and Blower, 1997) is a stratified sampling without replacement technique. The LHS technique takes k random parameter distributions, divides them into N equal probability intervals, and then samples from each interval. N is the sample size determined a priori and it is at least k + 1, where k is the number of model parameters varied, but it should be much larger to ensure accuracy. The sampling is performed by randomly selecting parameter values exactly once (without replacement) from each parameter intervals, thereby exploring the entire range for each parameter. The LHS technique performs the sampling independently for each parameter. Following the sampling, an N × k matrix (LHS matrix) is generated, where N is the number of simulations (sample size) and k is the number of parameters varied. N model solutions are then simulated, using each combination of parameter values, these correspond to each row of the LHS matrix. The model response (or output) is collected for each model simulation. J2. Partial rank correlation coefficient (PRCC) The partial rank correlation coefficient (PRCC) analysis, determines the sensitivity of the model to each parameter. Correlation provides a measure of the strength of a linear association between an input (parameters) and an output (the model response functions). The correlation coefficient (CC) between xj and y is calculated as follows:

rxj y =

Cov (xj , y ) Var (xj ) Var (y )

N i=1

=

N i=1

(x ij

(x ij

x¯)(yij x¯)2

N i=1

y¯) (yij

y¯)2

,

j = 1, 2,

,k

with rxj y [ 1, + 1]. Cov(xj, y) represents the covariance between xj and y, while Var(xj) and Var(y) are the variance of xj and the variance of y respectively, and x¯ and y¯ are the sample means. The significance of a non-zero PRCC can be tested by computing the following statistic (see Eq. (7) in Marino et al. (2008) and p. 143 in Anderson (2003)) for each PRCC (γ):

T =

N

2 1

p 2

tN

2 p

(J-1)

53

Ecological Complexity 38 (2019) 31–55

F.B. Agusto, et al.

where T follows a student’s t distribution with (N 2 p ) degrees of freedom. N is the sample size and p is the number of inputs/parameters whose effects are discounted when the PRCC is calculated. The significance tests determines if a PRCC is significantly different from zero or if two PRCC values are significantly different from each other.

Gebru, T., Ajua, A., Theisen, M., Esen, M., Ngoa, U., Issifou, S., Adegnika, A., Kremsner, P., Mordmüller, B., Held, J., 2017. Recognition of plasmodium falciparum mature gametocyte-infected erythrocytes by antibodies of semi-immune adults and malariaexposed children from gabon. Malar. J. 16 (1), 176. Gilchrist, M., Coombs, D., 2006. Evolution of virulence: interdependence, constraints, and selection using nested models. Theor. Popul. Biol. 69 (2), 145–153. Gilchrist, M., Sasaki, A., 2002. Modeling host–parasite coevolution: a nested approach based on mechanistic models. J. Theor. Biol. 218 (3), 289–308. Gouagna, L., Bonnet, S., Gounoue, R., Verhave, J., Eling, W., Sauerwein, R., Boudin, C., 2004. Stage-specific effects of host plasma factors on the early sporogony of autologous plasmodium falciparum isolates within anopheles gambiae. Trop. Med. Int. Health 9 (9), 937–948. Gouagna, L.C., Mulder, B., Noubissi, E., Tchuinkam, T., Boudin, C., 1998. The early sporogonic cycle of plasmodium falciparum in laboratory-infected anopheles gambiae: an estimation of parasite efficacy. Trop. Med. Int. Health 3 (1), 21–28. Hellriegel, B., 1992. Modelling the immune response to malaria with ecological concepts: short-term behaviour against long-term equilibrium. Proc. R. Soc. Lond. B 250 (1329), 249–256. Hethcote, H., 2000. The mathematics of infectious diseases. SIAM Rev. 42 (4), 599–653. Hoffman, R., Benz Jr, E., Silberstein, L., Heslop, H., Weitz, J., Anastasi, J., 2012. Hematology: Diagnosis and Treatment. Elsevier Health Sciences. Iggidr, A., Kamgang, J.-C., Sallet, G., Tewa, J.-J., 2006. Global analysis of new malaria intrahost models with a competitive exclusion principle. SIAM J. Appl. Math. 67 (1), 260–278. Kappe, S., Kaiser, K., Matuschewski, K., 2003. The plasmodium sporozoite journey: a rite of passage. Trends Parasitol. 19 (3), 135–143. Kelly, J., Williamson, S., Orive, M., Smith, M., Holt, R., 2003. Linking dynamical and population genetic models of persistent viral infection. Am. Nat. 162 (1), 14–28. Laperriere, V., Brugger, K., Rubel, F., 2011. Simulation of the seasonal cycles of bird, equine and human west nile virus cases. Prev. Vet. Med. 98 (2), 99–110. Legros, M., Bonhoeffer, S., 2016. A combined within-host and between-hosts modelling framework for the evolution of resistance to antimalarial drugs. J. R. Soc. Interface 13 (117), 20160148. Li, Y., Ruan, S., Xiao, D., 2011. The within-host dynamics of malaria infection with immune response. Math. Biosci. Eng. 8 (4), 999–1018. Lou, Y., Zhao, X., 2010. A climate-based malaria transmission model with structured vector population. SIAM J. Appl. Math. 70 (6), 2023–2044. Maier, A., Cooke, B., Cowman, A., Tilley, L., 2009. Malaria parasite proteins that remodel the host erythrocyte. Nat. Rev. Microbiol. 7 (5), 341–354. Marino, S., Hogue, I., Ray, C., Kirschner, D.E., 2008. A methodology for performing global uncertainty and sensitivity analysis in systems biology. J. Theor. Biol. 254 (1), 178–196. McKay, M., Beckman, R., Conover, W., 2000. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 42 (1), 55–61. McQueen, P., Williamson, K., McKenzie, F., 2013. Host immune constraints on malaria transmission: insights from population biology of within-host parasites. Malar. J. 12 (1), 206. Mideo, N., Alizon, S., Day, T., 2008. Linking within-and between-host dynamics in the evolutionary epidemiology of infectious diseases. Trends Ecol. Evol. 23 (9), 511–517. Mideo, N., Nelson, W., Reece, S., Bell, A., Read, A., Day, T., 2011. Bridging scales in the evolution of infectious disease life histories: application. Evolution 65 (11), 3298–3310. Murall, C., McCann, K., Bauch, C., 2012. Food webs in the human body: linking ecological theory to viral dynamics. PLoS ONE 7 (11), e48812. Naotunne, T.d.S., Karunaweera, N., Mendis, K., Carter, R., 1993. Cytokine-mediated inactivation of malarial gametocytes is dependent on the presence of white blood cells and involves reactive nitrogen intermediates. Immunology 78 (4), 555. Niger, A., Gumel, A., 2008. Mathematical analysis of the role of repeated exposure on malaria transmission dynamics. Differ. Equ. Dyn. Syst. 16 (3), 251–287. Niger, A., Gumel, A., 2011. Immune response and imperfect vaccine in malaria dynamics. Math. Popul. Stud. 18 (2), 55–86. Orive, M., Stearns, M., Kelly, J., Barfield, M., Smith, M., Holt, R., 2005. Viral infection in internally structured hosts. I. Conditions for persistent infection. J. Theor. Biol. 232 (4), 453–466. Park, M., Loverdo, C., Schreiber, S., Lloyd-Smith, J., 2013. Multiple scales of selection influence the evolutionary emergence of novel pathogens. Philos. Trans. R. Soc. Lond. B 368 (1614), 20120333. Portugal, S., Carret, C., Recker, M., Armitage, A., Gonçalves, L., Epiphanio, S., Sullivan, D., Roy, C., Newbold, C., Drakesmith, H., et al., 2011. Host-mediated regulation of superinfection in malaria. Nat. Med. 17 (6), 732. Richman, A., Dimopoulos, G., Seeley, D., Kafatos, F., 1997. Plasmodium activates the innate immune response of anopheles gambiae mosquitoes. EMBO J. 16 (20), 6114–6119. Robert, V., Macintyre, K., Keating, J., Trape, J., Duchemin, J., Warren, M., Beier, J., 2003. Malaria transmission in urban sub-Saharan Africa. Am. J. Trop. Med. Hyg. 68 (2), 169–176. Rubel, F., Brugger, K., Hantel, M., Chvala-Mannsberger, S., Bakonyi, T., Weissenböck, H., Nowotny, N., 2008. Explaining usutu virus dynamics in Austria: model development

References Agusto, F., Del Valle, S., Blayneh, K., Ngonghala, C., Goncalves, M., Li, N., Zhao, R., Gong, H., 2013. The impact of bed-net use on malaria prevalence. J. Theor. Biol. 320, 58–65. Agusto, F., Lenhart, S., 2013. Optimal control of the spread of malaria superinfectivity. J. Biol. Syst. 21 (04), 1340002. Agusto, F., Marcus, N., Okosun, K., 2012. Application of optimal control to the epidemiology of malaria. Electron. J. Differ. Equ. 2012 (81), 1–22. Alavi, Y., Arai, M., Mendoza, J., Tufet-Bayona, M., Sinha, R., Fowler, K., Billker, O., Franke-Fayard, B., Janse, C., Waters, A., et al., 2003. The dynamics of interactions between plasmodium and the mosquito: a study of the infectivity of plasmodium berghei and plasmodium gallinaceum, and their transmission by anopheles stephensi, anopheles gambiae and aedes aegypti. Int. J. Parasitol. 33 (9), 933–943. Anderson, R., May, R., Gupta, S., 1989. Non-linear phenomena in hostparasite interactions. Parasitology 99 (S1), S59–S79. Anderson, T.W., 2003. An introduction to multivariate statistical analysis, third ed. Wiley Series in Probability and Statistics. Wiley-Interscience, Hoboken, NJ. Aneke, S., 2002. Mathematical modelling of drug resistant malaria parasites and vector populations. Math. Methods Appl. Sci. 25 (4), 335–346. Babiker, H., Schneider, P., Reece, S., 2008. Gametocytes: insights gained during a decade of molecular monitoring. Trends Parasitol. 24 (11), 525–530. Barfield, M., Orive, M., Holt, R.D., 2015. The role of pathogen shedding in linking withinand between-host pathogen dynamics. Math. Biosci. 270, 249–262. Baron, S., 1996. Epidemiology–Medical Microbiology. University of Texas Medical Branch at Galveston. Baton, L., Ranford-Cartwright, L., 2005. Spreading the seeds of million-murdering death: metamorphoses of malaria in the mosquito. Trends Parasitol. 21 (12), 573–580. Beauchemin, C., 2006. Probing the effects of the well-mixed assumption on viral infection dynamics. J. Theor. Biol. 242 (2), 464–477. Beier, J., 1998. Malaria parasite development in mosquitoes. Annu. Rev. Entomol. 43 (1), 519–543. Blair, C., 2011. Mosquito RNAi is the major innate immune pathway controlling arbovirus infection and transmission. Future Microbiol. 6 (3), 265–277. Blandin, S., Levashina, E., 2004. Mosquito immune responses against malaria parasites. Curr. Opin. Immunol. 16 (1), 16–20. Blower, S., Dowlatabadi, H., 1994. Sensitivity and uncertainty analysis of complex models of disease transmission: an HIV model, as an example. Int. Stat. Rev. 229–243. Bousema, T., Churcher, T., Morlais, I., Dinglasan, R., 2013. Can field-based mosquito feeding assays be used for evaluating transmission-blocking interventions? Trends Parasitol. 29 (2), 53–59. Cen, X., Feng, Z., Zhao, Y., 2014. Emerging disease dynamics in a model coupling withinhost and between-host systems. J. Theor. Biol. 361, 141–151. Centers for Disease Control and Prevention, 2015. Malaria. http://www.cdc.gov/malaria/ about/biology/. Chitnis, N., Cushing, J., Hyman, J., 2006. Bifurcation analysis of a mathematical model for malaria transmission. SIAM J. Appl. Math. 67 (1), 24–45. Chiyaka, C., Garira, W., Dube, S., 2008. Modelling immune response and drug therapy in human malaria infection. Comput. Math. Methods Med. 9 (2), 143–163. Cunnington, A., Riley, E., Walther, M., 2013. Stuck in a rut? Reconsidering the role of parasite sequestration in severe malaria syndromes. Trends Parasitol. 29 (12), 585–592. Da, D., Churcher, T., Yerbanga, R., Yaméogo, B., Sangaré, I., Ouedraogo, J., Sinden, R., Blagborough, A., Cohuet, A., 2015. Experimental study of the relationship between plasmodium gametocyte density and infection success in mosquitoes; implications for the evaluation of malaria transmission-reducing interventions. Exp. Parasitol. 149, 74–83. Dawes, E., Churcher, T., Zhuang, S., Sinden, R., Basáñez, M.-G., 2009. Anopheles mortality is both age-and plasmodium-density dependent: implications for malaria transmission. Malar. J. 8 (1), 228. Day, T., Alizon, S., Mideo, N., 2011. Bridging scales in the evolution of infectious disease life histories: theory. Evolution 65 (12), 3448–3461. Diekmann, O., Heesterbeek, J., Metz, J., 1990. On the definition and computation of the basic reproduction ratio r0 in models for infectious diseases in heterogeneous populations. J. Math. Biol. 28, 503–522. Feng, Z., Velasco-Hernandez, J., Tapia-Santos, B., 2013. A mathematical model for coupling within-host and between-host dynamics in an environmentally-driven infectious disease. Math. Biosci. 241 (1), 49–55. Feng, Z., Velasco-Hernandez, J., Tapia-Santos, B., Leite, M., 2012. A model for coupling within-host and between-host dynamics in an infectious disease. Nonlinear Dyn. 68 (3), 401–411. Ferguson, H., Read, A., 2004. Mosquito appetite for blood is stimulated by plasmodium chabaudi infections in themselves and their vertebrate hosts. Malar. J. 3 (1), 12. Frost, S., Dumaurier, M., Wain-Hobson, S., Brown, A.L., 2001. Genetic drift and withinhost metapopulation dynamics of HIV-1 infection. Proc. Natl. Acad. Sci. 98 (12), 6975–6980. Funk, G., Jansen, V., Bonhoeffer, S., Killingback, T., 2005. Spatial models of virus-immune dynamics. J. Theor. Biol. 233 (2), 221–236.

54

Ecological Complexity 38 (2019) 31–55

F.B. Agusto, et al. and calibration. Prev. Vet. Med. 85 (3), 166–186. Sanchez, M., Blower, S., 1997. Uncertainty and sensitivity analysis of the basic reproductive rate: tuberculosis as an example. Am. J. Epidemiol. 145 (12), 1127–1137. Santos, G., Torres, N., 2013. New targets for drug discovery against malaria. PLoS ONE 8 (3), e59968. Sreenivasamurthy, S., Dey, G., Ramu, M., Kumar, M., Gupta, M., Mohanty, A., Harsha, H., Sharma, P., Kumar, N., Pandey, A.o., 2013. A compendium of molecules involved in vector-pathogen interactions pertaining to malaria. Malar. J. 12 (1), 1. Sutherland, C., 2009. Surface antigens of plasmodium falciparum gametocytesa new class of transmission-blocking vaccine targets? Mol. Biochem. Parasitol. 166 (2), 93–98. Teboh-Ewungkem, M.I, Yuster, T., 2010. A within-vector mathematical model of plasmodium falciparum and implications of incomplete fertilization on optimal gametocyte sex ratio. J. Theor. Biol. 264 (2), 273–286. Thieme, H., 2003. Mathematics in Population Biology. Princeton University Press.

Tumwiine, J., Mugisha, J., Luboobi, L., 2008. On global stability of the intra-host dynamics of malaria and the immune system. J. Math. Anal. Appl. 341 (2), 855–869. Van den Driessche, P., Watmough, J., 2002. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math. Biosci. 180 (1), 29–48. Vaughan, J., 2007. Population dynamics of plasmodium sporogony. Trends Parasitol. 23 (2), 63–70. White, N., 2017. Malaria parasite clearance. Malar. J. 16 (1), 88. Wiser, M., 2011. Malaria. http://www.tulane.edu/~wiser/protozoology/notes/malaria. html. World Health Organization, 2017. Malaria. http://www.who.int/mediacentre/ factsheets/fs094/en/. Yassine, H., Osta, M., 2010. Anopheles gambiae innate immunity. Cell. Microbiol. 12 (1), 1–9.

55