# A Stochastic Model for Malaria Transmission Dynamics

Authors: Rachel Waema Mbogo; Livingstone S. Luboobi; John W. Odhiambo
Journal: Journal of Applied Mathematics
 (2018)

---

## Abstract

Malaria is one of the three most dangerous infectious diseases worldwide (along with HIV/AIDS and tuberculosis). In this paper we compare the disease dynamics of the deterministic and stochastic models in order to determine the effect of randomness in malaria transmission dynamics. Relationships between the basic reproduction number for malaria transmission dynamics between humans and mosquitoes and the extinction thresholds of corresponding continuous-time Markov chain models are derived under certain assumptions. The stochastic model is formulated using the continuous-time discrete state Galton-Watson branching process (CTDSGWbp). The reproduction number of deterministic models is an essential quantity to predict whether an epidemic will spread or die out. Thresholds for disease extinction from stochastic models contribute crucial knowledge on disease control and elimination and mitigation of infectious diseases. Analytical and numerical results show some significant differences in model predictions between the stochastic and deterministic models. In particular, we find that malaria outbreak is more likely if the disease is introduced by infected mosquitoes as opposed to infected humans. These insights demonstrate the importance of a policy or intervention focusing on controlling the infected mosquito population if the control of malaria is to be realized.


---

## Body

## 1. Introduction

Malaria is an infectious disease caused by the Plasmodium parasite and transmitted between humans through bites of femaleAnopheles mosquitoes. Approximately half of the world’s population is at risk of malaria. It remains one of the most prevalent and lethal human infections throughout the world. An estimated 40% of the world’s population lives in malaria endemic areas. Most cases and deaths occur in sub-Saharan Africa. It causes an estimated 300 to 500 million cases and 1.5 to 2.7 million deaths each year worldwide. Africa shares 80% of the cases and 90% of deaths [1]. According to the website of the World Health Organization [2] there were approximately 214 million new cases of malaria and 438,000 deaths worldwide in 2015. Most cases were reported in the African region.Recently, the incidence of malaria has been rising due to drug resistance. Various control strategies have been taken to reduce malaria transmissions. Since the first mathematical model of malaria transmission was introduced by Ross [3], quite a number of mathematical models have been formulated to investigate the transmission dynamics of malaria. Xiao and Zou [4] used mathematical models to explore a natural concern of possible epidemics caused by multiple species of malaria parasites in one region. They found that epidemics involving both species in a single region are possible.Li and others [5] considered fast and slow dynamics of malaria model with relapse and analyzed the global dynamics by using the geometric singular perturbation theory. They suggested that a treatment should be given to symptomatic patients completely and adequately rather than to asymptomatic patients. On the other hand, for the asymptomatic patients, their results strongly suggested that, to control and eradicate the malaria, it is very necessary for governments to control the relapse rate strictly. Relapse is when symptoms reappear after the parasites had been eliminated from blood but persist as dormant hypnozoites in liver cells [6]. This commonly occurs between 8 and 24 weeks and is commonly seen withP. vivax andP. ovale infections. Other papers also consider the influence of relapse in giving up smoking or quitting drinking; please see [7].Chitnis et al. [8] and Li et al. [5] assumed that the recovered humans have some immunity to the disease and do not get clinically ill but they still harbour low levels of parasite in their blood streams and can pass the infection to mosquitoes. After some period of time, they lose their immunity and return to the susceptible class. Unfortunately, Li and others did not consider that the recovered humans will return to their infectious state because of incomplete treatment.Stochasticity is fundamental to biological systems. In some situations the system can be treated as a large number of similar agents interacting in a homogeneously mixing environment, and so the dynamics are well-captured by deterministic ordinary differential equations. However, in many situations, the system can be driven by a small number of agents or strongly influenced by an environment fluctuating in space and time [9–12].Stochastic models incorporate discrete movements of individuals between epidemiological classes and not average rates at which individuals move between classes [13–15]. In stochastic epidemic models, numbers in each class are integers and not continuously varying quantities [13]. A significant possibility is that the last infected individual can recover before the disease is transmitted and the infection can only reoccur if it is reintroduced from outside the population [16]. In contrast, most deterministic models have the flaw that infections can fall to very low levels—well below the point at which there is only one infected individual only to rise up later [17]. In addition, the variability introduced in stochastic models may result in dynamics that differ from the predictions made by deterministic models [16].For a large population size and a large number of infectious individuals, the deterministic thresholdR0>1 provides a good prediction of a disease outbreak. However, this prediction breaks down when the outbreak is initiated by a small number of infectious individuals. In this setting, Markov chain (MC) models with a discrete number of individuals are more realistic than deterministic models where the number of individuals is assumed to be continuous-valued [18].Motivated by these works, in this paper, we propose a model which is an extension of the model formulated by Huo and Qui (2014), who assumed that the pseudorecovered humans can recover and return to the susceptible class or relapse and become infectious again. Using the extended model, we will formulate the basic reproductive numberR0 and use it to compare the disease dynamics of the deterministic and stochastic models in order to determine the effect of randomness in malaria transmission dynamics.This paper is organized as follows; in Section2, we present a malaria transmission deterministic model with relapse, which is an extension of the model in [6]. We compute the basic reproduction number, R0, of the malaria transmission deterministic model using the next-generation matrix approach. The stochastic version of the deterministic model and its underlying assumptions necessary for model formulation are presented and discussed in Section 3. In this section, we also compute the stochastic threshold for disease extinction or invasion by applying the multitype Galton-Watson branching process. In Section 4, we show the relationship between reproductive number of the deterministic model and the thresholds for disease extinction of the stochastic version; we also illustrate our results using numerical simulations. We conclude with a discussion of the results in Section 5.


## 2. The Malaria Deterministic Model

The interaction of the Mosquito-host is shown in Figure1.Figure 1
The mosquito-human transmission cycle. An infectious mosquito bites a susceptible, uninfected human and transmits the virus via saliva. Once the human becomes infectious (usually accompanied by symptoms), the human can transmit the pathogen to an uninfected mosquito via the blood the mosquito ingests.Source. Figure 1 is reproduced from Manore and Hyman (2016) [20] [under the Creative Commons Attribution License/public domain].


### 2.1. Model Formulation

In this section, we introduce a deterministic model of malaria with relapse, an extended form of the model in [6]. The model proposes a more realistic mathematical model of malaria, in which it is assumed that the pseudorecovered humans interact with infected mosquitoes and may acquire more parasites causing reinfection or, due to incomplete treatment, the infection my reoccur (relapse) and return to infectious class. Also the pseudorecovered human may lose their immunity and return to the susceptible class. Given that humans might get repeatedly infected due to not acquiring complete immunity, then the population dynamics are assumed to be described by the SIRS model; hence we consider a deterministic compartmental model which divides the total human population size at time t, denoted by N(t), into susceptible individuals Sh(t) (those who are not currently harbouring the parasite but are liable to be infected), infectious individuals Ih(t) (those already infected and are able to transmit the disease to mosquitoes), and pseudorecovered individuals Rh(t) (those who are treated from the disease but with partial recovery and hence can transmit the disease to mosquitoes). Mosquitoes are assumed not to recover from the parasites so the mosquito population can be described by the SI model.The structure of model is shown in Figure2.Figure 2
Schematic representation of the mosquito-human malaria transmission dynamics.


### 2.2. Variable and Parameter Description for the Model

The variables for the model are summarized in description of variables for the malaria transmission model in the Notations.The parameters for the model are described as in description of parameters for the malaria transmission model in the Notations.


### 2.3. Equations of the Model

Assuming that the disease transmits in a closed system which translates into the simplifying assumption of a constant population size with the birth rate equal to death rate, that is,λ=μ (see [21]) so that N(t)=N, hence, the above assumptions lead to the following system of differential equations which describe the interaction between mosquitoes and humans: (1)dShtdt=λN-βShImN+ρ1Rh-μSh,dIhtdt=βShImN+ρ2RhImN-γ+μ1Ih,dRhtdt=γIh-ρ1+ρ2ImN+μ1Rh,dSmtdt=δM-α1SmIhN-α2SmRhN-ηSm,dImtdt=α1SmIhN+α2SmRhN-η1Im,where Sh, Ih, Rh, Sm, and Im represent the number of susceptible humans, infectious humans, recovered humans, susceptible mosquitoes, and infectious mosquitoes, respectively.


### 2.4. Computation ofR0 Using the Next-Generation Matrix Approach

The basic reproduction number,R0, is defined as the secondary infections produced by one infective agent that is introduced into an entirely susceptible population at the disease-free equilibrium. The next-generation matrix approach is frequently used to compute R0. The original nonlinear system of ODEs including these compartments can be written as ∂Xi/∂t=F-V, where F=(Fi) and V=(Vi) represent new infections and transfer between compartments, respectively [22–24].The Jacobian matrices ofF(X) and V(X) at the disease-free equilibrium E0 are, respectively, (2)F=FE0=00β000α1α20V=VE0=γ+μ1-ρ20-γμ1+ρ1+ρ2000η1.The matrix J=F-V is the Jacobian matrix evaluated at the DFE. (3)J=F-V=-γ-μ1ρ2βγ-μ1-ρ1-ρ20α1α2-η1.The inverse matrix of V is(4)V-1=μ1+ρ1+ρ2γ+μ1μ1+ρ1+μ1ρ2ρ2γ+μ1μ1+ρ1+μ1ρ20γγ+μ1μ1+ρ1+μ1ρ2γ+μ1γ+μ1μ1+ρ1+μ1ρ20001η1.FV-1=00βη000γα2+α1μ+ρ1+ρ2γ+μμ+ρ1+μρ2γ+μα2+α1ρ2γ+μμ+ρ1+μρ20.The matrix FV-1 is called the next-generation matrix. The (i,k) entry of FV-1 indicates the expected number of new infections in compartment i produced by the infected individual originally introduced into compartment k. The model reproduction number, R0, which is defined as the spectral radius of FV-1 and denoted by ρ(FV-1), is given by(5)R0=ρFV-1=βγα2+α1ρ1+ρ2+μ1η1γ+μ1ρ1+μ1+μ1ρ2.Here R0 is associated with disease transmission by infected humans as well as the infection of susceptible humans by infected mosquitoes.Simplifying (5) we find an equivalent form, which is defined as the product of the transmission from mosquito to human and from human to mosquito as follows: (6)R0=βη1α1γ+μ1+α2ρ1+ρ2+μ1.In (6), R0 is associated with disease transmission by infected mosquitoes as well as infection of susceptible mosquitoes by infected humans. The term β/η1 represents the number of infected humans generated by infectious mosquito in its life span. The term α1/(γ+μ1) represents the number of infected mosquitoes generated by infectious human during the infectious period of the individual while the term α2/(ρ1+ρ2+μ1) represents the number of infected mosquitoes generated by pseudorecovered human during his/her infectious period. Susceptible humans acquire infection following effective contacts with infected mosquitoes. Susceptible mosquitoes acquire malaria infection from infected humans in two ways, namely, by infectious human or pseudorecovered humans.Also, from (6), we can rewrite R0 as infection between mosquitoes and infectious humans or infection between mosquitoes and pseudorecovered humans. (7)R0=R01+R02,where (8)R01=βη1α1γ+μ1,(9)R02=βη1α2ρ1+ρ2+μ1.R01 in (8) represents the product of individuals generated by infectious human and infected mosquito, respectively, during their infectious time.R
02 in (9) represents the product of individuals generated by pseudorecovered human and infected mosquito, respectively, during their infectious time.For model (1), the disease dies out if R0<1 and the disease persists if R0>1. Hence R0<1, iff R01<1 and R02<1.


## 2.1. Model Formulation

In this section, we introduce a deterministic model of malaria with relapse, an extended form of the model in [6]. The model proposes a more realistic mathematical model of malaria, in which it is assumed that the pseudorecovered humans interact with infected mosquitoes and may acquire more parasites causing reinfection or, due to incomplete treatment, the infection my reoccur (relapse) and return to infectious class. Also the pseudorecovered human may lose their immunity and return to the susceptible class. Given that humans might get repeatedly infected due to not acquiring complete immunity, then the population dynamics are assumed to be described by the SIRS model; hence we consider a deterministic compartmental model which divides the total human population size at time t, denoted by N(t), into susceptible individuals Sh(t) (those who are not currently harbouring the parasite but are liable to be infected), infectious individuals Ih(t) (those already infected and are able to transmit the disease to mosquitoes), and pseudorecovered individuals Rh(t) (those who are treated from the disease but with partial recovery and hence can transmit the disease to mosquitoes). Mosquitoes are assumed not to recover from the parasites so the mosquito population can be described by the SI model.The structure of model is shown in Figure2.Figure 2
Schematic representation of the mosquito-human malaria transmission dynamics.


## 2.2. Variable and Parameter Description for the Model

The variables for the model are summarized in description of variables for the malaria transmission model in the Notations.The parameters for the model are described as in description of parameters for the malaria transmission model in the Notations.


## 2.3. Equations of the Model

Assuming that the disease transmits in a closed system which translates into the simplifying assumption of a constant population size with the birth rate equal to death rate, that is,λ=μ (see [21]) so that N(t)=N, hence, the above assumptions lead to the following system of differential equations which describe the interaction between mosquitoes and humans: (1)dShtdt=λN-βShImN+ρ1Rh-μSh,dIhtdt=βShImN+ρ2RhImN-γ+μ1Ih,dRhtdt=γIh-ρ1+ρ2ImN+μ1Rh,dSmtdt=δM-α1SmIhN-α2SmRhN-ηSm,dImtdt=α1SmIhN+α2SmRhN-η1Im,where Sh, Ih, Rh, Sm, and Im represent the number of susceptible humans, infectious humans, recovered humans, susceptible mosquitoes, and infectious mosquitoes, respectively.


## 2.4. Computation ofR0 Using the Next-Generation Matrix Approach

The basic reproduction number,R0, is defined as the secondary infections produced by one infective agent that is introduced into an entirely susceptible population at the disease-free equilibrium. The next-generation matrix approach is frequently used to compute R0. The original nonlinear system of ODEs including these compartments can be written as ∂Xi/∂t=F-V, where F=(Fi) and V=(Vi) represent new infections and transfer between compartments, respectively [22–24].The Jacobian matrices ofF(X) and V(X) at the disease-free equilibrium E0 are, respectively, (2)F=FE0=00β000α1α20V=VE0=γ+μ1-ρ20-γμ1+ρ1+ρ2000η1.The matrix J=F-V is the Jacobian matrix evaluated at the DFE. (3)J=F-V=-γ-μ1ρ2βγ-μ1-ρ1-ρ20α1α2-η1.The inverse matrix of V is(4)V-1=μ1+ρ1+ρ2γ+μ1μ1+ρ1+μ1ρ2ρ2γ+μ1μ1+ρ1+μ1ρ20γγ+μ1μ1+ρ1+μ1ρ2γ+μ1γ+μ1μ1+ρ1+μ1ρ20001η1.FV-1=00βη000γα2+α1μ+ρ1+ρ2γ+μμ+ρ1+μρ2γ+μα2+α1ρ2γ+μμ+ρ1+μρ20.The matrix FV-1 is called the next-generation matrix. The (i,k) entry of FV-1 indicates the expected number of new infections in compartment i produced by the infected individual originally introduced into compartment k. The model reproduction number, R0, which is defined as the spectral radius of FV-1 and denoted by ρ(FV-1), is given by(5)R0=ρFV-1=βγα2+α1ρ1+ρ2+μ1η1γ+μ1ρ1+μ1+μ1ρ2.Here R0 is associated with disease transmission by infected humans as well as the infection of susceptible humans by infected mosquitoes.Simplifying (5) we find an equivalent form, which is defined as the product of the transmission from mosquito to human and from human to mosquito as follows: (6)R0=βη1α1γ+μ1+α2ρ1+ρ2+μ1.In (6), R0 is associated with disease transmission by infected mosquitoes as well as infection of susceptible mosquitoes by infected humans. The term β/η1 represents the number of infected humans generated by infectious mosquito in its life span. The term α1/(γ+μ1) represents the number of infected mosquitoes generated by infectious human during the infectious period of the individual while the term α2/(ρ1+ρ2+μ1) represents the number of infected mosquitoes generated by pseudorecovered human during his/her infectious period. Susceptible humans acquire infection following effective contacts with infected mosquitoes. Susceptible mosquitoes acquire malaria infection from infected humans in two ways, namely, by infectious human or pseudorecovered humans.Also, from (6), we can rewrite R0 as infection between mosquitoes and infectious humans or infection between mosquitoes and pseudorecovered humans. (7)R0=R01+R02,where (8)R01=βη1α1γ+μ1,(9)R02=βη1α2ρ1+ρ2+μ1.R01 in (8) represents the product of individuals generated by infectious human and infected mosquito, respectively, during their infectious time.R
02 in (9) represents the product of individuals generated by pseudorecovered human and infected mosquito, respectively, during their infectious time.For model (1), the disease dies out if R0<1 and the disease persists if R0>1. Hence R0<1, iff R01<1 and R02<1.


## 3. Malaria Stochastic Epidemic Model

For the mosquito-human malaria dynamics, since time is continuous and the disease states are discrete, we derive the stochastic version of the deterministic model (1) using a continuous-time discrete state Galton-Watson branching process (CTDSGWbp), which is a type of stochastic process. The malaria CTDSGWbp model is a time-homogeneous process with the Markov property. The model takes into account random effects of individual birth and death processes, that is, demographic variability. A stochastic process is defined by the probabilities with which different events happen in a small time interval Δt. In our model there are two possible events (production and death/removal) for each population. The corresponding rates in the deterministic model are replaced in the stochastic version by the probabilities that any of these events occur in a small time interval of length Δt [25, 26].


### 3.1. Model Formulation

Let time be continuous,t∈[0;∞), and let Sh(t), Ih(t), Rh(t), Sm(t), and Im(t) represent discrete random variables for the number of susceptible humans, infectious humans, recovered humans, susceptible mosquitoes and infectious mosquitoes, respectively, with finite space.(10)Sht,Iht,Rht,Smt,Imt∈0,1,2,3,…,H,where H is positive and represents the maximum size of the populations.If a disease emerges from one infectious group withR0>1 and if i infective agents are introduced into a wholly susceptible population, then the probability of a major disease outbreak is approximated by 1-(1/R0)i while the probability of disease extinction is approximately (1/R0)i [13]. However, this result does not hold if the infection emanates from multiple infectious groups [27]. For multiple infectious groups, the stochastic thresholds depend on two factors, namely, the number of individuals in each group and the probability of disease extinction for each group. Further, the persistence of an infection into a wholly susceptible population is not guaranteed by having R0 greater than one.For CTDSGWbp models, the transition from one state to a new state may occur at any timet. If the process begins in a state G(0), after a random time period τ, it transits to a new state G(τ). The process remains in state G(τ) for a random time t after which it moves through to the new state G(p), with p=τ+t [27]. This process continues throughout the model. The state transitions and rates for the stochastic model are presented in Table 1.Table 1
State transitions and rates for the CTDSGWbp malaria model.
Event
Population components att
Population components att+Δt
Transition probabilities
Birth of humans
S
h
,
I
h
,
R
h
,
S
m
,
I
m
S
h
+
1
,
I
h
,
R
h
,
S
m
,
I
m
λ
Δ
t
Death of susceptible humans
S
h
,
I
h
,
R
h
,
S
m
,
I
m
S
h
-
1
,
I
h
,
R
h
,
S
m
,
I
m
μ
S
h
Δ
t
Infection of humans
S
h
,
I
h
,
R
h
,
S
m
,
I
m
S
h
-
1
,
I
h
+
1
,
R
h
,
S
m
,
I
m
β(Sh/N)ImΔt
Recovery rate of humans
S
h
,
I
h
,
R
h
,
S
m
,
I
m
S
h
+
1
,
I
h
,
R
h
-
1
,
S
m
,
I
m
ρ
1
R
h
Δ
t
Relapse rate of humans
S
h
,
I
h
,
R
h
,
S
m
,
I
m
S
h
,
I
h
+
1
,
R
h
-
1
,
S
m
,
I
m
ρ
2
R
h
Δ
t
Treatment rate
S
h
,
I
h
,
R
h
,
S
m
,
I
m
S
h
,
I
h
-
1
,
R
h
+
1
,
S
m
,
I
m
γ
I
h
Δ
t
Death of infected humans
S
h
,
I
h
,
R
h
,
S
m
,
I
m
S
h
,
I
h
-
1
,
R
h
,
S
m
,
I
m
μ
1
I
h
Δ
t
Death of recovered humans
S
h
,
I
h
,
R
h
,
S
m
,
I
m
S
h
,
I
h
,
R
h
-
1
,
S
m
,
I
m
μ
1
R
h
Δ
t
Birth of mosquitoes
S
h
,
I
h
,
R
h
,
S
m
,
I
m
S
h
,
I
h
,
R
h
,
S
m
+
1
,
I
m
δ
Δ
t
Death of susceptible mosquitoes
S
h
,
I
h
,
R
h
,
S
m
,
I
m
S
h
,
I
h
,
R
h
,
S
m
-
1
,
I
m
η
S
m
Δ
t
Infection of mosquitoes from humans
S
h
,
I
h
,
R
h
,
S
m
,
I
m
S
h
,
I
h
,
R
h
,
S
m
-
1
,
I
m
+
1
α
1
I
h
+
α
2
R
h
(
S
m
/
N
)
Δ
t
Death of infected mosquitoes
S
h
,
I
h
,
R
h
,
S
m
,
I
m
S
h
,
I
h
,
R
h
,
S
m
,
I
m
-
1
η
1
I
m
Δ
t


### 3.2. The Branching Process Approximation

We use branching process to analyze the malaria dynamics near the disease-free equilibrium. Since infectious human, pseudorecovered humans, and infectious mosquitoes are the only sources of infection, we apply the multitype branching process in the three variablesIh(t), Rh(t), and IM(t). The susceptible humans and mosquitoes are assumed to be at the disease-free state [28]. We use the multitype Galton-Watson branching process (GWbp) to determine disease invasion and extinction probabilities. More review on the GWbp branching theory process can be accessed through [18, 27, 28]. We now define the offspring pgfs for the three variables, where each offspring pgf has a general form (11)Gis1,s2,…,sn=∑kn=0∞∑kn-1=0∞⋯∑k1=0∞Pik1,k2,…,kns1k1s2k2⋯snkn,where Pi(k1,k2,…,kn)=prob(Xi1=k1,Xi2=k2,…,Xin=kn) is the probability that one infected individual of type i gives birth to kj individuals of type j; see [29].For one malaria infectious human, there are three possible events: infection of a mosquito, recovery of the infectious human, or death of the infectious human. The offspring pgf for infectious humans define the probabilities associated with the “birth” of secondary infectious mosquito or the “death” of the initial infectious human, given that the process started with only one infectious human; that isIh(0)=1, Rh(0)=0, and IM(0)=0.The offspring pgf forIh is given by(12)G1s1,s2,s3=∑k1=0∞∑k2=0∞∑k3=0∞P1k1,k2,k3s1k1s2k2s3k3,where P1(k1,k2,k3)=prob(X11=k1,X12=k2,X13=k3) is the probability that one infectious human through infection produces another infectious human k1 or a pseudorecovered human k2 or an infectious mosquito k3.Similarly, the offspring pgf forRh is given by(13)G2s1,s2,s3=∑k1=0∞∑k2=0∞∑k3=0∞P2k1,k2,k3s1k1s2k2s3k3,where P2(k1,k2,k3)=prob(X21=k1,X22=k2,X23=k3) is the probability that one pseudorecovered human through infection produces an infectious human k1 or another pseudorecovered human k2 or an infectious mosquito k3.Lastly, the offspring pgf forIm is given by(14)G3s1,s2,s3=∑k1=0∞∑k2=0∞∑k3=0∞P3k1,k2,k3s1k1s2k2s3k3,where P3(k1,k2,k3)=prob(X31=k1,X32=k2,X33=k3) is the probability that one infectious mosquito through infection produces an infections human k1 or a pseudorecovered human k2 or another infectious mosquito k3.The power to whichsj is raised is the number of infectious individuals generated from one infectious individual. If an individual recovers or dies, then no new infections are generated, hence (sj0).The offspring pgfs forIh, Rh, and Im are used to calculate the expected number of offsprings produced by a single infectious human or by pseudorecovered human or infectious mosquito. They are also used to calculate the probability of disease extinction.The specific offspring pgfs forIh, Rh, and Im are defined using the rates in description of parameters for the malaria transmission model in the Notations, when the initial susceptible populations are near disease-free equilibrium, Sh(0)≈N and Sm(0)≈M.From (12), the offspring pgf for infectious human, given Ih(0)=1, Rh(0)=0, and IM(0)=0, is given by (15)G1s1,s2,s3=α1s1s3+γ+μ1α1+γ+μ1.For G1, one infectious human dies or is treated with probability (μ1+γ)/(α1+γ+μ1); this means an infectious human dies before infecting a susceptible mosquito. The term γ/(α1+γ+μ1) represents the probability that the infectious human is treated and moves to the pseudorecovered class, this results in “X11=0, X12=01, and X13=0, though there is movement of an infectious human to pseudorecovered state due to partial treatment (this is not new offspring).” The infectious human infects a mosquito with probability α1/(α1+γ+μ1). This means an infectious human infects a susceptible mosquito and remains infectious, which results in “X11=1, X12=0, and X13=1.” Note the term s1s3 in (15) means one infectious human generates one infectious mosquito (s3 raised to power one) and remains infectious (s1 raised to power one).For one pseudorecovered human, there are four events: infection of a mosquito, relapse to infected class, successful treatment of the pseudorecovered human, or death of the recovered host. Similarly, from (13), the offspring pgf for recovered humans given Ih(0)=0, Rh(0)=1, and IM(0)=0 is given by (16)G2s1,s2,s3=α2s2s3+ρ1+ρ2+μ1α2+ρ1+ρ2+μ1.For G2, one pseudorecovered host dies or relapses or is fully treated with probability (ρ1+ρ2+μ1)/(α2+ρ1+ρ2+μ1) or infects a mosquito with probability α2/(α2+ρ1+ρ2+μ1). This means a pseudorecovered human infects a susceptible mosquito and remains infectious, which results in “X21=0, X22=1, and X23=1.” Note the term s2s3 in (16) means one recovered human generates one infectious mosquito (s3 raised to power one) and remains infectious (s2 raised to power one).For one infectious mosquito, there are only two events: infection of a susceptible human or death of the mosquito. From (14), the offspring pgf for infected mosquito given Ih(0)=0, Rh(0)=0, and IM(0)=1 is given by (17)G3s1,s2,s3=βs1s3+η1β+η1.For G3, one infectious mosquito dies with probability η1/(β+η1) or infects a human with probability β/(β+η1). This means an infectious mosquito infects a susceptible human and remains infectious, which results in “X31=1, X32=0, and X33=1.” Note the term s1s3 in (17) means one infectious mosquito generates one infectious human (s1 raised to power one) and remains infectious (s3 raised to power one).


### 3.3. The Relationship betweenR0 and the Stochastic Threshold S0

The offspring pgfs, evaluated at(1,1,1), gives the expectation matrix with elements mji. Below are the offspring pgfs evaluated at (1,1,1). (18)∂G1∂s1=α1α+γ+μ1,∂G1∂s2=γ,∂G1∂s3=α1α+γ+μ1,∂G2∂s1=ρ2,∂G2∂s2=α2α2+ρ1+ρ2+μ1,∂G2∂s3=α2α2+ρ1+ρ2++μ1,∂G3∂s1=ββ+η1,∂G3∂s2=0,∂G3∂s3=ββ+η1.The expectation matrix of the offspring pgfs, evaluated at (1,1,1), is given by (19)M=∂G1∂s1∂G2∂s1∂G3∂s1∂G1∂s2∂G2∂s2∂G3∂s2∂G1∂s3∂G2∂s3∂G3∂s3=α1γ+μ1+α1ρ2γ+μ1+α1α1α+γ+μ1γμ1+α2+ρ1+ρ2α2μ1+α2+ρ1+ρ2α2μ1+α2+ρ1+ρ2ββ+η10ββ+η1.The entries m11, m21, and m31 represent the expected number of infectious humans, pseudorecovered humans, and infectious mosquitoes, respectively, produced by one infectious human. Similarly, the entries m12, m22, and m32 represent the expected number of infectious humans, pseudorecovered humans, and infectious mosquitoes, respectively, produced by one pseudorecovered human. Lastly, the entries m13, m23, and m33 represent the expected number of infectious humans, pseudorecovered humans, and infectious mosquitoes, respectively, produced by one infectious mosquito.The matrixC=WM-I is the Jacobian matrix evaluated for the stochastic version.(20)C=WM-I=-γ-μ1γα1ρ2-μ1-ρ1-ρ2α2β0-η1,where W=diagγ+μ1+α1,μ1+α2+ρ1+ρ2,β+η1 is a diagonal matrix and I is the identity matrix.From (3) and (20), we show that (21)F-V=WM-IT=-γ-μ1ρ2βγ-μ1-ρ1-ρ20α1α2-η1.The spectral radius of matrix M obtained by finding the eigenvalues of matrix M is given by (22)S0=ρM=max⁡βγ+μ1+2β+η1α1β+η1γ+μ1+α1,α2α2+ρ1+ρ2+μ1.Taking the first expression of S0=ρ(M) in (24), we have (23)S0=ρM=βγ+μ1+α1+α1β+η1β+η1γ+μ1+α1=ββ+η1+α1γ+μ1+α1.This gives the probability of malaria transmission by either infectious human or by infectious mosquito. From (23), (24)S0=ρM=βγ+μ1+α1+α1β+η1β+η1γ+μ1+α1.The probability of disease extinction is one if ρ(M)<1. Hence from (24) we have (25)βγ+μ1+α1+α1β+η1<β+η1γ+μ1+α1.Expanding and simplifying the inequality, we have(26)βα1<η1γ+μ1which reduces to (27)βη1α1γ+μ1<1⟹R01<1.When the infection is between infectious humans and infectious mosquitoes, (27) is true. The result in (27) agrees with the deterministic reproduction number for disease elimination. Hence we conclude that the probability of disease elimination in the CTDSGWbp model is one iff (28)ρM<1⟹R0<1.


### 3.4. Deriving Probability of Disease ExtinctionP0 Using Branching Process Approximation

To find the probability of extinction (no outbreak), we compute the fixed points of the system(q1,q2,q3)∈(0,1) of the offspring pgfs for the three infectious stages; that is, we solve GIh=q1, GRh=q2, and GIm=q3. The solutions of these systems are (1,1,1) and (q1,q2,q3); see [19]. Equating G1(s1,s2,s3) in (15) to q1, then letting s1=q1 and s3=q3, and solving for q1, we have (29)q1=β+η1γ+μ1βγ+μ1+α1=γ+μ1α1+γ+μ1+α1α1+γ+μ11R01.Equating G2(s1,s2,s3) in (16) to q2, letting s2=q2 and s3=q3, and solving for q2, we have (30)q2=β+η1α1μ1+ρ1+ρ2-η1γ+μ1α2+α1βα2+β+η1μ1+ρ1+ρ2.Equating G3(s1,s2,s3) in (17) to q3, then letting s1=q1 and s3=q3, and solving for q3, we have (31)q3=η1γ+μ1+α1β+η1α1=η1β+η1+ββ+η11R01.The expression for q1 in (29) has a biological interpretation. Beginning from one infectious human, there is no outbreak if the infectious human recovers or dies with probability (γ+μ1)/(α1+γ+μ1) or if there is no successful transmission to a susceptible mosquito with probability (α1/α1+γ+μ1)(1/R01). This implies that if there is successful contact, then the probability of successful transmission from infectious human to susceptible mosquito is 1-1/R01.The expression forq3 in (31) has a biological interpretation. Beginning from one infectious mosquito, there is no outbreak if the infectious mosquito dies with probability η1/(β+η1) or if there is no successful transmission to a susceptible human with probability (β/(β+η1))(1/R01).There are some other important relationships; if disease transmission is by infectious mosquitoes as well as infection of susceptible mosquitoes by infectious humans, then(32)q1∗q3=1R01.From (32), we see that, in both mosquito and human populations, the probability of no successful transmission from infectious human to susceptible mosquito and from infected mosquito to susceptible human is 1/R01.If disease transmission is by infectious mosquitoes as well as infection of susceptible mosquitoes by pseudorecovered humans, then(33)q2∗q3=1R02.From (33), we see that, in both mosquito and human populations, the probability of no successful transmission from pseudorecovered humans to susceptible mosquito and from infected mosquito to pseudorecovered humans is 1/R02.To compute the probability of disease extinction and of an outbreak for our malaria model, we recall that, for multiple infectious groups, the stochastic thresholds depend on two factors, namely, the number of initial individuals in each group and the probability of disease extinction for each group. Usingq1, q2, and q3 in (29)–(31) and assuming initial individuals for infectious humans, pseudorecovered humans and infected mosquitoes are Ih(0)=h0, Rh(0)=r0, and Im(0)=m0, respectively. Then the probability of malaria clearance is given by (34)P0=q1h0∗q2r0∗q3m0=β+η1γ+μ1βγ+μ1+α1h0∗α1β+η1μ1+ρ1+ρ2-η1α2γ+μ1+α1+α1β+η1α2+μ1+ρ1+ρ2r0∗η1γ+μ1+α1β+η1α1m0.


## 3.1. Model Formulation

Let time be continuous,t∈[0;∞), and let Sh(t), Ih(t), Rh(t), Sm(t), and Im(t) represent discrete random variables for the number of susceptible humans, infectious humans, recovered humans, susceptible mosquitoes and infectious mosquitoes, respectively, with finite space.(10)Sht,Iht,Rht,Smt,Imt∈0,1,2,3,…,H,where H is positive and represents the maximum size of the populations.If a disease emerges from one infectious group withR0>1 and if i infective agents are introduced into a wholly susceptible population, then the probability of a major disease outbreak is approximated by 1-(1/R0)i while the probability of disease extinction is approximately (1/R0)i [13]. However, this result does not hold if the infection emanates from multiple infectious groups [27]. For multiple infectious groups, the stochastic thresholds depend on two factors, namely, the number of individuals in each group and the probability of disease extinction for each group. Further, the persistence of an infection into a wholly susceptible population is not guaranteed by having R0 greater than one.For CTDSGWbp models, the transition from one state to a new state may occur at any timet. If the process begins in a state G(0), after a random time period τ, it transits to a new state G(τ). The process remains in state G(τ) for a random time t after which it moves through to the new state G(p), with p=τ+t [27]. This process continues throughout the model. The state transitions and rates for the stochastic model are presented in Table 1.Table 1
State transitions and rates for the CTDSGWbp malaria model.
Event
Population components att
Population components att+Δt
Transition probabilities
Birth of humans
S
h
,
I
h
,
R
h
,
S
m
,
I
m
S
h
+
1
,
I
h
,
R
h
,
S
m
,
I
m
λ
Δ
t
Death of susceptible humans
S
h
,
I
h
,
R
h
,
S
m
,
I
m
S
h
-
1
,
I
h
,
R
h
,
S
m
,
I
m
μ
S
h
Δ
t
Infection of humans
S
h
,
I
h
,
R
h
,
S
m
,
I
m
S
h
-
1
,
I
h
+
1
,
R
h
,
S
m
,
I
m
β(Sh/N)ImΔt
Recovery rate of humans
S
h
,
I
h
,
R
h
,
S
m
,
I
m
S
h
+
1
,
I
h
,
R
h
-
1
,
S
m
,
I
m
ρ
1
R
h
Δ
t
Relapse rate of humans
S
h
,
I
h
,
R
h
,
S
m
,
I
m
S
h
,
I
h
+
1
,
R
h
-
1
,
S
m
,
I
m
ρ
2
R
h
Δ
t
Treatment rate
S
h
,
I
h
,
R
h
,
S
m
,
I
m
S
h
,
I
h
-
1
,
R
h
+
1
,
S
m
,
I
m
γ
I
h
Δ
t
Death of infected humans
S
h
,
I
h
,
R
h
,
S
m
,
I
m
S
h
,
I
h
-
1
,
R
h
,
S
m
,
I
m
μ
1
I
h
Δ
t
Death of recovered humans
S
h
,
I
h
,
R
h
,
S
m
,
I
m
S
h
,
I
h
,
R
h
-
1
,
S
m
,
I
m
μ
1
R
h
Δ
t
Birth of mosquitoes
S
h
,
I
h
,
R
h
,
S
m
,
I
m
S
h
,
I
h
,
R
h
,
S
m
+
1
,
I
m
δ
Δ
t
Death of susceptible mosquitoes
S
h
,
I
h
,
R
h
,
S
m
,
I
m
S
h
,
I
h
,
R
h
,
S
m
-
1
,
I
m
η
S
m
Δ
t
Infection of mosquitoes from humans
S
h
,
I
h
,
R
h
,
S
m
,
I
m
S
h
,
I
h
,
R
h
,
S
m
-
1
,
I
m
+
1
α
1
I
h
+
α
2
R
h
(
S
m
/
N
)
Δ
t
Death of infected mosquitoes
S
h
,
I
h
,
R
h
,
S
m
,
I
m
S
h
,
I
h
,
R
h
,
S
m
,
I
m
-
1
η
1
I
m
Δ
t


## 3.2. The Branching Process Approximation

We use branching process to analyze the malaria dynamics near the disease-free equilibrium. Since infectious human, pseudorecovered humans, and infectious mosquitoes are the only sources of infection, we apply the multitype branching process in the three variablesIh(t), Rh(t), and IM(t). The susceptible humans and mosquitoes are assumed to be at the disease-free state [28]. We use the multitype Galton-Watson branching process (GWbp) to determine disease invasion and extinction probabilities. More review on the GWbp branching theory process can be accessed through [18, 27, 28]. We now define the offspring pgfs for the three variables, where each offspring pgf has a general form (11)Gis1,s2,…,sn=∑kn=0∞∑kn-1=0∞⋯∑k1=0∞Pik1,k2,…,kns1k1s2k2⋯snkn,where Pi(k1,k2,…,kn)=prob(Xi1=k1,Xi2=k2,…,Xin=kn) is the probability that one infected individual of type i gives birth to kj individuals of type j; see [29].For one malaria infectious human, there are three possible events: infection of a mosquito, recovery of the infectious human, or death of the infectious human. The offspring pgf for infectious humans define the probabilities associated with the “birth” of secondary infectious mosquito or the “death” of the initial infectious human, given that the process started with only one infectious human; that isIh(0)=1, Rh(0)=0, and IM(0)=0.The offspring pgf forIh is given by(12)G1s1,s2,s3=∑k1=0∞∑k2=0∞∑k3=0∞P1k1,k2,k3s1k1s2k2s3k3,where P1(k1,k2,k3)=prob(X11=k1,X12=k2,X13=k3) is the probability that one infectious human through infection produces another infectious human k1 or a pseudorecovered human k2 or an infectious mosquito k3.Similarly, the offspring pgf forRh is given by(13)G2s1,s2,s3=∑k1=0∞∑k2=0∞∑k3=0∞P2k1,k2,k3s1k1s2k2s3k3,where P2(k1,k2,k3)=prob(X21=k1,X22=k2,X23=k3) is the probability that one pseudorecovered human through infection produces an infectious human k1 or another pseudorecovered human k2 or an infectious mosquito k3.Lastly, the offspring pgf forIm is given by(14)G3s1,s2,s3=∑k1=0∞∑k2=0∞∑k3=0∞P3k1,k2,k3s1k1s2k2s3k3,where P3(k1,k2,k3)=prob(X31=k1,X32=k2,X33=k3) is the probability that one infectious mosquito through infection produces an infections human k1 or a pseudorecovered human k2 or another infectious mosquito k3.The power to whichsj is raised is the number of infectious individuals generated from one infectious individual. If an individual recovers or dies, then no new infections are generated, hence (sj0).The offspring pgfs forIh, Rh, and Im are used to calculate the expected number of offsprings produced by a single infectious human or by pseudorecovered human or infectious mosquito. They are also used to calculate the probability of disease extinction.The specific offspring pgfs forIh, Rh, and Im are defined using the rates in description of parameters for the malaria transmission model in the Notations, when the initial susceptible populations are near disease-free equilibrium, Sh(0)≈N and Sm(0)≈M.From (12), the offspring pgf for infectious human, given Ih(0)=1, Rh(0)=0, and IM(0)=0, is given by (15)G1s1,s2,s3=α1s1s3+γ+μ1α1+γ+μ1.For G1, one infectious human dies or is treated with probability (μ1+γ)/(α1+γ+μ1); this means an infectious human dies before infecting a susceptible mosquito. The term γ/(α1+γ+μ1) represents the probability that the infectious human is treated and moves to the pseudorecovered class, this results in “X11=0, X12=01, and X13=0, though there is movement of an infectious human to pseudorecovered state due to partial treatment (this is not new offspring).” The infectious human infects a mosquito with probability α1/(α1+γ+μ1). This means an infectious human infects a susceptible mosquito and remains infectious, which results in “X11=1, X12=0, and X13=1.” Note the term s1s3 in (15) means one infectious human generates one infectious mosquito (s3 raised to power one) and remains infectious (s1 raised to power one).For one pseudorecovered human, there are four events: infection of a mosquito, relapse to infected class, successful treatment of the pseudorecovered human, or death of the recovered host. Similarly, from (13), the offspring pgf for recovered humans given Ih(0)=0, Rh(0)=1, and IM(0)=0 is given by (16)G2s1,s2,s3=α2s2s3+ρ1+ρ2+μ1α2+ρ1+ρ2+μ1.For G2, one pseudorecovered host dies or relapses or is fully treated with probability (ρ1+ρ2+μ1)/(α2+ρ1+ρ2+μ1) or infects a mosquito with probability α2/(α2+ρ1+ρ2+μ1). This means a pseudorecovered human infects a susceptible mosquito and remains infectious, which results in “X21=0, X22=1, and X23=1.” Note the term s2s3 in (16) means one recovered human generates one infectious mosquito (s3 raised to power one) and remains infectious (s2 raised to power one).For one infectious mosquito, there are only two events: infection of a susceptible human or death of the mosquito. From (14), the offspring pgf for infected mosquito given Ih(0)=0, Rh(0)=0, and IM(0)=1 is given by (17)G3s1,s2,s3=βs1s3+η1β+η1.For G3, one infectious mosquito dies with probability η1/(β+η1) or infects a human with probability β/(β+η1). This means an infectious mosquito infects a susceptible human and remains infectious, which results in “X31=1, X32=0, and X33=1.” Note the term s1s3 in (17) means one infectious mosquito generates one infectious human (s1 raised to power one) and remains infectious (s3 raised to power one).


## 3.3. The Relationship betweenR0 and the Stochastic Threshold S0

The offspring pgfs, evaluated at(1,1,1), gives the expectation matrix with elements mji. Below are the offspring pgfs evaluated at (1,1,1). (18)∂G1∂s1=α1α+γ+μ1,∂G1∂s2=γ,∂G1∂s3=α1α+γ+μ1,∂G2∂s1=ρ2,∂G2∂s2=α2α2+ρ1+ρ2+μ1,∂G2∂s3=α2α2+ρ1+ρ2++μ1,∂G3∂s1=ββ+η1,∂G3∂s2=0,∂G3∂s3=ββ+η1.The expectation matrix of the offspring pgfs, evaluated at (1,1,1), is given by (19)M=∂G1∂s1∂G2∂s1∂G3∂s1∂G1∂s2∂G2∂s2∂G3∂s2∂G1∂s3∂G2∂s3∂G3∂s3=α1γ+μ1+α1ρ2γ+μ1+α1α1α+γ+μ1γμ1+α2+ρ1+ρ2α2μ1+α2+ρ1+ρ2α2μ1+α2+ρ1+ρ2ββ+η10ββ+η1.The entries m11, m21, and m31 represent the expected number of infectious humans, pseudorecovered humans, and infectious mosquitoes, respectively, produced by one infectious human. Similarly, the entries m12, m22, and m32 represent the expected number of infectious humans, pseudorecovered humans, and infectious mosquitoes, respectively, produced by one pseudorecovered human. Lastly, the entries m13, m23, and m33 represent the expected number of infectious humans, pseudorecovered humans, and infectious mosquitoes, respectively, produced by one infectious mosquito.The matrixC=WM-I is the Jacobian matrix evaluated for the stochastic version.(20)C=WM-I=-γ-μ1γα1ρ2-μ1-ρ1-ρ2α2β0-η1,where W=diagγ+μ1+α1,μ1+α2+ρ1+ρ2,β+η1 is a diagonal matrix and I is the identity matrix.From (3) and (20), we show that (21)F-V=WM-IT=-γ-μ1ρ2βγ-μ1-ρ1-ρ20α1α2-η1.The spectral radius of matrix M obtained by finding the eigenvalues of matrix M is given by (22)S0=ρM=max⁡βγ+μ1+2β+η1α1β+η1γ+μ1+α1,α2α2+ρ1+ρ2+μ1.Taking the first expression of S0=ρ(M) in (24), we have (23)S0=ρM=βγ+μ1+α1+α1β+η1β+η1γ+μ1+α1=ββ+η1+α1γ+μ1+α1.This gives the probability of malaria transmission by either infectious human or by infectious mosquito. From (23), (24)S0=ρM=βγ+μ1+α1+α1β+η1β+η1γ+μ1+α1.The probability of disease extinction is one if ρ(M)<1. Hence from (24) we have (25)βγ+μ1+α1+α1β+η1<β+η1γ+μ1+α1.Expanding and simplifying the inequality, we have(26)βα1<η1γ+μ1which reduces to (27)βη1α1γ+μ1<1⟹R01<1.When the infection is between infectious humans and infectious mosquitoes, (27) is true. The result in (27) agrees with the deterministic reproduction number for disease elimination. Hence we conclude that the probability of disease elimination in the CTDSGWbp model is one iff (28)ρM<1⟹R0<1.


## 3.4. Deriving Probability of Disease ExtinctionP0 Using Branching Process Approximation

To find the probability of extinction (no outbreak), we compute the fixed points of the system(q1,q2,q3)∈(0,1) of the offspring pgfs for the three infectious stages; that is, we solve GIh=q1, GRh=q2, and GIm=q3. The solutions of these systems are (1,1,1) and (q1,q2,q3); see [19]. Equating G1(s1,s2,s3) in (15) to q1, then letting s1=q1 and s3=q3, and solving for q1, we have (29)q1=β+η1γ+μ1βγ+μ1+α1=γ+μ1α1+γ+μ1+α1α1+γ+μ11R01.Equating G2(s1,s2,s3) in (16) to q2, letting s2=q2 and s3=q3, and solving for q2, we have (30)q2=β+η1α1μ1+ρ1+ρ2-η1γ+μ1α2+α1βα2+β+η1μ1+ρ1+ρ2.Equating G3(s1,s2,s3) in (17) to q3, then letting s1=q1 and s3=q3, and solving for q3, we have (31)q3=η1γ+μ1+α1β+η1α1=η1β+η1+ββ+η11R01.The expression for q1 in (29) has a biological interpretation. Beginning from one infectious human, there is no outbreak if the infectious human recovers or dies with probability (γ+μ1)/(α1+γ+μ1) or if there is no successful transmission to a susceptible mosquito with probability (α1/α1+γ+μ1)(1/R01). This implies that if there is successful contact, then the probability of successful transmission from infectious human to susceptible mosquito is 1-1/R01.The expression forq3 in (31) has a biological interpretation. Beginning from one infectious mosquito, there is no outbreak if the infectious mosquito dies with probability η1/(β+η1) or if there is no successful transmission to a susceptible human with probability (β/(β+η1))(1/R01).There are some other important relationships; if disease transmission is by infectious mosquitoes as well as infection of susceptible mosquitoes by infectious humans, then(32)q1∗q3=1R01.From (32), we see that, in both mosquito and human populations, the probability of no successful transmission from infectious human to susceptible mosquito and from infected mosquito to susceptible human is 1/R01.If disease transmission is by infectious mosquitoes as well as infection of susceptible mosquitoes by pseudorecovered humans, then(33)q2∗q3=1R02.From (33), we see that, in both mosquito and human populations, the probability of no successful transmission from pseudorecovered humans to susceptible mosquito and from infected mosquito to pseudorecovered humans is 1/R02.To compute the probability of disease extinction and of an outbreak for our malaria model, we recall that, for multiple infectious groups, the stochastic thresholds depend on two factors, namely, the number of initial individuals in each group and the probability of disease extinction for each group. Usingq1, q2, and q3 in (29)–(31) and assuming initial individuals for infectious humans, pseudorecovered humans and infected mosquitoes are Ih(0)=h0, Rh(0)=r0, and Im(0)=m0, respectively. Then the probability of malaria clearance is given by (34)P0=q1h0∗q2r0∗q3m0=β+η1γ+μ1βγ+μ1+α1h0∗α1β+η1μ1+ρ1+ρ2-η1α2γ+μ1+α1+α1β+η1α2+μ1+ρ1+ρ2r0∗η1γ+μ1+α1β+η1α1m0.


## 4. Numerical Simulations

In this section, we illustrate numerically the disease dynamics of model (1) using parameter values in Table 2. The numerical simulations are done using Maple codes.Table 2
Model parameter values.
Parameter
Description
Units
Parameter values
Source
λ
Birth rate of humans
Per day
0.000039
[5]
μ
Death rate of susceptible humans
Per day
0.000039
[5]
β
Infection rate of humans
Per day
0.02
[19]
μ
1
Death rate of infected humans
Per day
0.00039
[5]
ρ
1
Recovery rate of humans
Per day
0.01
[19]
ρ
2
Relapse rate of humans
Per day
0.002
[5]
γ
Treatment rate
Per day
0.037
[5]
μ
1
Death rate of recovered humans
Per day
0.00034
[5]
δ
Birth rate of mosquitoes
Per day
0.143
[5]
η
Death rate of susceptible mosquitoes
Per day
0.143
[5]
α
1
Infection rate from infectious human
Per day
0.072
[5]
α
2
Infection rate from recovered human
Per day
0.0072
[5]
η
1
Death rate of infected mosquitoes
Per day
0.143
[19]


### 4.1. Effects of Model Parameters onR0

Using parameter values in Table2, we identified how different input parameters affect the reproduction number R0 of our model as shown in Figure 3.Figure 3
Tornado plot for parameters that influenceR0.From Figure3, we see that increase in γ, ρ1, μ1, and η1 will decrease R0. Also decrease in β, α1, α2, and ρ2 will decrease R0.From the Tornado plot, the infection of susceptible humans by infected mosquitoes (denoted byβ) is a major factor in the malaria transmission dynamics. Reducing β would reduce R0 significantly hence reducing the possibility of disease outbreak. Vector control is the main way to prevent and reduce malaria transmission. If coverage of vector control interventions within a specific area is high enough, then a measure of protection will be conferred across the community. WHO recommends protection for all people at risk of malaria with effective malaria vector control which agrees with the information from Figure 3. Two forms of vector control insecticide-treated mosquito nets and indoor residual spraying are effective in a wide range of circumstances. Another factor which affects the malaria transmission significantly is the death rate of infected mosquitoes (denoted by η1). Increasing η1 would decrease R0 hence reducing the chances of disease outbreak.


#### 4.1.1. Graphical Representation of Parameters Effects onR0

Figure4 shows the effects of relapse and recover rates on the basic reproduction number. From the simulations, we find that R0 is increasing with increase in relapse rate, while it is decreasing with increase in recovery rate. To control and eradicate the malaria epidemic, it is important and necessary for governments of endemic countries to decrease the relapse rate and increase the treatment and recovery rate.Figure 4
Effects of relapse and recovery rate onR0.
(a)
Effects of relapse onR0
(b)
Effects of treatment rate onR0
(c)
Effects of recovery rate onR0


### 4.2. Probability of Disease Extinction

Using parameter values in Table2, we compute numerically the probability of disease extinction P0 and of an outbreak 1-P0 for our malaria model using different initial values for the infectious classes.The probability of disease extinction is high if the disease emerges from infected humans. It is very low if the infection emerges from infected mosquitoes. However, as the initial number of infected humans grows largely, there is a high probability of disease outbreak as illustrated in Table3. The probability of disease extinction is significantly low if the disease emerges from infected mosquitoes. Therefore, the disease dynamics for model system (1) at the beginning of the epidemics are being driven by initial number of infected mosquitoes.Table 3
Probability of disease extinctionP0 and of an outbreak 1-P0 for the malaria model.
h
0
r
0
m
0
P
0
1
-
P
0
R
01
R
02
R
0
1
0
0
0.6838
0.3162
R
01
>
1
R
02
<
1
R
0
>
1
0
1
0
0.9834
0.0166
R
01
>
1
R
02
<
1
R
0
>
1
0
0
1
0.1365
0.8635
R
01
>
1
R
02
<
1
R
0
>
1
1
1
0
0.6724
0.3276
R
01
>
1
R
02
<
1
R
0
>
1
1
1
1
0.0918
0.9082
R
01
>
1
R
02
<
1
R
0
>
1
10
0
0
0.0223
0.9777
R
01
>
1
R
02
<
1
R
0
>
1
0
10
0
0.8459
0.1541
R
01
>
1
R
02
<
1
R
0
>
1
0
0
10
0
1
R
01
>
1
R
02
<
1
R
0
>
1
10
10
0
0.0188
0.9812
R
01
>
1
R
02
<
1
R
0
>
1
10
10
10
0
1
R
01
>
1
R
02
<
1
R
0
>
1A female mosquito will continue to bite and draw blood until her abdomen is full. If she is interrupted before she is full, she will fly to the next person. After feeding, the mosquito rests for two or three days before laying her eggs, and then it is ready to bite again. Since mosquitoes do not recover from the infection, then a single infected mosquito is capable of biting and infecting so many susceptible humans in their lifespan which may in turn infect so many uninfected mosquitoes there by reducing the probability of disease clearance and increasing the probability of a major disease outbreak. Moreover, mosquitoes are the reservoir host for the parasites that cause malaria and it takes a long time for them to be malaria-free, hence the high probability of malaria outbreak if the parasite is introduced by infected mosquito.Table3 depicts that, at the beginning of malaria outbreak, any policy or intervention to control the spread of malaria should focus on controlling the infected mosquito population as well as the infected humans. If more effort to control the disease is only focused on the infected humans, then it is very unlikely that malaria will be eliminated. This is an interesting insight from the stochastic threshold that could not be provided by the deterministic threshold.


### 4.3. Numerical Simulation of Malaria Model

Using parameter values in Table2, we numerically simulate the behavior of model (1). Initial conditions are Sh(0)=99, Ih(0)=1, Rh(0)=0, Sm(0)=999, and Im(0)=1.


#### 4.3.1. Numerical Simulation WhenR0<1

From Figure5, the analysis shows that when R0<1 and S0<1, then the probability of disease extinction is P0=0.9476≃1, although this agrees with (30), which points out that the probability of disease elimination in the CTDSGWbp model is one iff (35)S0<1⟹R0<1.There is still a small probability 1-P0=0.0524 of disease outbreak. R0 has been widely used as a measure of disease dynamics to estimate the effectiveness of control measures and to inform on disease management policy. However, from the analysis in Figure 5, it is evident that R0 can be flawed and disease can persist with R0<1 depending on the kind of disease being modeled.Figure 5
Malaria dynamics whenR0=0.4564 and S0=0.5661 and P0=0.9476.
(a)
Human population dynamics whenR0<1 and S0<1
(b)
Mosquito population dynamics whenR0<1 and S0<1


#### 4.3.2. Numerical Simulation WhenR0>1

Using parameter values in Table2, we numerically simulate the behavior of model (1) when R0>1. Initial conditions are Sh(0)=99, Ih(0)=1, Rh(0)=0, Sm(0)=999, and Im(0)=1.The analysis from Figure6 suggests that when R0>1, there is still some probability of disease extinction (P0=0.3713). Although the probability is low, there is still a chance to clear the disease. Therefore, from the analysis in Figure 6, disease can be eliminated with R0>1 and hence the R0 threshold should not be the only parameter to consider in quantifying the spread of a disease.Figure 6
Malaria dynamics whenR0=2.4629 and P0=0.3713.
(a)
Human population dynamics whenR0>1
(b)
Mosquito population dynamics whenR0>1


#### 4.3.3. Effects of Relapse on Human Population

Figure7 shows changing effects of relapse on the human populations.Figure 7
Changing effects of relapse on the human populations.
(a)
Human population dynamics when relapse rate is increased;ρ2=0.4
(b)
Human population dynamics when relapse rate is reduced;ρ2=0.001To control and eradicate the malaria epidemic, it is important and necessary for governments of endemic countries to decrease the relapse rate as can be seen in Figure7.


## 4.1. Effects of Model Parameters onR0

Using parameter values in Table2, we identified how different input parameters affect the reproduction number R0 of our model as shown in Figure 3.Figure 3
Tornado plot for parameters that influenceR0.From Figure3, we see that increase in γ, ρ1, μ1, and η1 will decrease R0. Also decrease in β, α1, α2, and ρ2 will decrease R0.From the Tornado plot, the infection of susceptible humans by infected mosquitoes (denoted byβ) is a major factor in the malaria transmission dynamics. Reducing β would reduce R0 significantly hence reducing the possibility of disease outbreak. Vector control is the main way to prevent and reduce malaria transmission. If coverage of vector control interventions within a specific area is high enough, then a measure of protection will be conferred across the community. WHO recommends protection for all people at risk of malaria with effective malaria vector control which agrees with the information from Figure 3. Two forms of vector control insecticide-treated mosquito nets and indoor residual spraying are effective in a wide range of circumstances. Another factor which affects the malaria transmission significantly is the death rate of infected mosquitoes (denoted by η1). Increasing η1 would decrease R0 hence reducing the chances of disease outbreak.


### 4.1.1. Graphical Representation of Parameters Effects onR0

Figure4 shows the effects of relapse and recover rates on the basic reproduction number. From the simulations, we find that R0 is increasing with increase in relapse rate, while it is decreasing with increase in recovery rate. To control and eradicate the malaria epidemic, it is important and necessary for governments of endemic countries to decrease the relapse rate and increase the treatment and recovery rate.Figure 4
Effects of relapse and recovery rate onR0.
(a)
Effects of relapse onR0
(b)
Effects of treatment rate onR0
(c)
Effects of recovery rate onR0


## 4.1.1. Graphical Representation of Parameters Effects onR0

Figure4 shows the effects of relapse and recover rates on the basic reproduction number. From the simulations, we find that R0 is increasing with increase in relapse rate, while it is decreasing with increase in recovery rate. To control and eradicate the malaria epidemic, it is important and necessary for governments of endemic countries to decrease the relapse rate and increase the treatment and recovery rate.Figure 4
Effects of relapse and recovery rate onR0.
(a)
Effects of relapse onR0
(b)
Effects of treatment rate onR0
(c)
Effects of recovery rate onR0


## 4.2. Probability of Disease Extinction

Using parameter values in Table2, we compute numerically the probability of disease extinction P0 and of an outbreak 1-P0 for our malaria model using different initial values for the infectious classes.The probability of disease extinction is high if the disease emerges from infected humans. It is very low if the infection emerges from infected mosquitoes. However, as the initial number of infected humans grows largely, there is a high probability of disease outbreak as illustrated in Table3. The probability of disease extinction is significantly low if the disease emerges from infected mosquitoes. Therefore, the disease dynamics for model system (1) at the beginning of the epidemics are being driven by initial number of infected mosquitoes.Table 3
Probability of disease extinctionP0 and of an outbreak 1-P0 for the malaria model.
h
0
r
0
m
0
P
0
1
-
P
0
R
01
R
02
R
0
1
0
0
0.6838
0.3162
R
01
>
1
R
02
<
1
R
0
>
1
0
1
0
0.9834
0.0166
R
01
>
1
R
02
<
1
R
0
>
1
0
0
1
0.1365
0.8635
R
01
>
1
R
02
<
1
R
0
>
1
1
1
0
0.6724
0.3276
R
01
>
1
R
02
<
1
R
0
>
1
1
1
1
0.0918
0.9082
R
01
>
1
R
02
<
1
R
0
>
1
10
0
0
0.0223
0.9777
R
01
>
1
R
02
<
1
R
0
>
1
0
10
0
0.8459
0.1541
R
01
>
1
R
02
<
1
R
0
>
1
0
0
10
0
1
R
01
>
1
R
02
<
1
R
0
>
1
10
10
0
0.0188
0.9812
R
01
>
1
R
02
<
1
R
0
>
1
10
10
10
0
1
R
01
>
1
R
02
<
1
R
0
>
1A female mosquito will continue to bite and draw blood until her abdomen is full. If she is interrupted before she is full, she will fly to the next person. After feeding, the mosquito rests for two or three days before laying her eggs, and then it is ready to bite again. Since mosquitoes do not recover from the infection, then a single infected mosquito is capable of biting and infecting so many susceptible humans in their lifespan which may in turn infect so many uninfected mosquitoes there by reducing the probability of disease clearance and increasing the probability of a major disease outbreak. Moreover, mosquitoes are the reservoir host for the parasites that cause malaria and it takes a long time for them to be malaria-free, hence the high probability of malaria outbreak if the parasite is introduced by infected mosquito.Table3 depicts that, at the beginning of malaria outbreak, any policy or intervention to control the spread of malaria should focus on controlling the infected mosquito population as well as the infected humans. If more effort to control the disease is only focused on the infected humans, then it is very unlikely that malaria will be eliminated. This is an interesting insight from the stochastic threshold that could not be provided by the deterministic threshold.


## 4.3. Numerical Simulation of Malaria Model

Using parameter values in Table2, we numerically simulate the behavior of model (1). Initial conditions are Sh(0)=99, Ih(0)=1, Rh(0)=0, Sm(0)=999, and Im(0)=1.


### 4.3.1. Numerical Simulation WhenR0<1

From Figure5, the analysis shows that when R0<1 and S0<1, then the probability of disease extinction is P0=0.9476≃1, although this agrees with (30), which points out that the probability of disease elimination in the CTDSGWbp model is one iff (35)S0<1⟹R0<1.There is still a small probability 1-P0=0.0524 of disease outbreak. R0 has been widely used as a measure of disease dynamics to estimate the effectiveness of control measures and to inform on disease management policy. However, from the analysis in Figure 5, it is evident that R0 can be flawed and disease can persist with R0<1 depending on the kind of disease being modeled.Figure 5
Malaria dynamics whenR0=0.4564 and S0=0.5661 and P0=0.9476.
(a)
Human population dynamics whenR0<1 and S0<1
(b)
Mosquito population dynamics whenR0<1 and S0<1


### 4.3.2. Numerical Simulation WhenR0>1

Using parameter values in Table2, we numerically simulate the behavior of model (1) when R0>1. Initial conditions are Sh(0)=99, Ih(0)=1, Rh(0)=0, Sm(0)=999, and Im(0)=1.The analysis from Figure6 suggests that when R0>1, there is still some probability of disease extinction (P0=0.3713). Although the probability is low, there is still a chance to clear the disease. Therefore, from the analysis in Figure 6, disease can be eliminated with R0>1 and hence the R0 threshold should not be the only parameter to consider in quantifying the spread of a disease.Figure 6
Malaria dynamics whenR0=2.4629 and P0=0.3713.
(a)
Human population dynamics whenR0>1
(b)
Mosquito population dynamics whenR0>1


### 4.3.3. Effects of Relapse on Human Population

Figure7 shows changing effects of relapse on the human populations.Figure 7
Changing effects of relapse on the human populations.
(a)
Human population dynamics when relapse rate is increased;ρ2=0.4
(b)
Human population dynamics when relapse rate is reduced;ρ2=0.001To control and eradicate the malaria epidemic, it is important and necessary for governments of endemic countries to decrease the relapse rate as can be seen in Figure7.


## 4.3.1. Numerical Simulation WhenR0<1

From Figure5, the analysis shows that when R0<1 and S0<1, then the probability of disease extinction is P0=0.9476≃1, although this agrees with (30), which points out that the probability of disease elimination in the CTDSGWbp model is one iff (35)S0<1⟹R0<1.There is still a small probability 1-P0=0.0524 of disease outbreak. R0 has been widely used as a measure of disease dynamics to estimate the effectiveness of control measures and to inform on disease management policy. However, from the analysis in Figure 5, it is evident that R0 can be flawed and disease can persist with R0<1 depending on the kind of disease being modeled.Figure 5
Malaria dynamics whenR0=0.4564 and S0=0.5661 and P0=0.9476.
(a)
Human population dynamics whenR0<1 and S0<1
(b)
Mosquito population dynamics whenR0<1 and S0<1


## 4.3.2. Numerical Simulation WhenR0>1

Using parameter values in Table2, we numerically simulate the behavior of model (1) when R0>1. Initial conditions are Sh(0)=99, Ih(0)=1, Rh(0)=0, Sm(0)=999, and Im(0)=1.The analysis from Figure6 suggests that when R0>1, there is still some probability of disease extinction (P0=0.3713). Although the probability is low, there is still a chance to clear the disease. Therefore, from the analysis in Figure 6, disease can be eliminated with R0>1 and hence the R0 threshold should not be the only parameter to consider in quantifying the spread of a disease.Figure 6
Malaria dynamics whenR0=2.4629 and P0=0.3713.
(a)
Human population dynamics whenR0>1
(b)
Mosquito population dynamics whenR0>1


## 4.3.3. Effects of Relapse on Human Population

Figure7 shows changing effects of relapse on the human populations.Figure 7
Changing effects of relapse on the human populations.
(a)
Human population dynamics when relapse rate is increased;ρ2=0.4
(b)
Human population dynamics when relapse rate is reduced;ρ2=0.001To control and eradicate the malaria epidemic, it is important and necessary for governments of endemic countries to decrease the relapse rate as can be seen in Figure7.


## 5. Discussions and Recommendations

In this study, we investigated the transmission dynamics of malaria using CTDSGWbp model. The disease dynamic extinction thresholds from the stochastic model were compared with the corresponding deterministic threshold. We derived the stochastic threshold for disease extinctionS0 and showed the relationship that exists between R0 and S0 in terms of disease extinction and outbreak in both deterministic and stochastic models.Our analytical and numerical results showed that both deterministic and stochastic models predict disease extinction whenR0<1 and S0<1. However, the predictions by these models are different when R0>1. In this case, deterministic model predicts with certainty disease outbreak while the stochastic model has a probability of disease extinction at the beginning of an infection. Hence, with stochastic models, it is possible to attain a disease-free equilibrium even when R0>1. Also we noticed that initial conditions do not affect the deterministic threshold while the stochastic thresholds are affected. Thus, the dynamics of the stochastic model are highly dependent on the initial conditions and should not be ignored.The probabilities of disease extinction for different initial sizes of infected humans and infected mosquitoes were approximated numerically. The results indicate that the probability of eliminating malaria is high if the disease emerges from infected human as opposed to when it emerges from infected mosquito at the beginning of the disease. The analysis has shown that any policy or intervention to control the spread of malaria at the beginning of an outbreak must focus not only on infected humans but also on controlling the infected mosquito population. Also our results strongly suggest that, to eradicate malaria, governments of endemic countries should increase the recovery rate and reduce the relapse rate.In conclusion, to achieve the WHO vision 2030 strategic global targets, which include (i) reducing malaria case incidence by at least90% by 2030 and (ii) reducing malaria mortality rates by at least 90% by 2030, the governments of endemic countries should embark on increasing the coverage of vector control interventions and reduce the malaria relapse rate as well as controlling the infected mosquito population.For future work, the study can be extended by considering vertical transmission of the mosquitoes. Another extension of the study can be the inclusion of immigration and emigration of individuals in both populations in order to investigate the effect of movement on the disease transmission dynamics on the two thresholds. Also the study can be extended by considering the effect of climatic conditions, such as rainfall patterns, temperature, and humidity on the disease transmission dynamics.



---

*Source: 2439520-2018-02-11.xml*