A Stochastic Model for Malaria Transmission Dynamics

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 formalaria transmission dynamics between humans and mosquitoes and the extinction thresholds of corresponding continuous-time Markov chain models are derived under certain assumptions.The stochasticmodel 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 infectedmosquitoes 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.


Introduction
Malaria is an infectious disease caused by the Plasmodium parasite and transmitted between humans through bites of female Anopheles 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 with P. vivax and P. 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][10][11][12].
Stochastic models incorporate discrete movements of individuals between epidemiological classes and not average rates at which individuals move between classes [13][14][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 threshold  0 > 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 number  0 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 Section 2, we present a malaria transmission deterministic model with relapse, which is an extension of the model in [6].We compute the basic reproduction number,  0 , 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.

The Malaria Deterministic Model
The interaction of the Mosquito-host is shown in Figure 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 , denoted by (), into susceptible individuals  ℎ () (those who are not currently harbouring the parasite but are liable to be infected), infectious individuals  ℎ () (those already infected and are able to transmit the disease to mosquitoes), and pseudorecovered individuals  ℎ () (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 Figure 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.

Human population
Mosquito population : Schematic representation of the mosquito-human malaria transmission dynamics.

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 () = , hence, the above assumptions lead to the following system of differential equations which describe the interaction between mosquitoes and humans: where  ℎ ,  ℎ ,  ℎ ,   , and   represent the number of susceptible humans, infectious humans, recovered humans, susceptible mosquitoes, and infectious mosquitoes, respectively.

Computation of 𝑅 0 Using the Next-Generation Matrix
Approach.The basic reproduction number,  0 , 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  0 .The original nonlinear system of ODEs including these compartments can be written as   / = F − V, where F = (F  ) and V = (V  ) represent new infections and transfer between compartments, respectively [22][23][24].The Jacobian matrices of F() and V() at the diseasefree equilibrium  0 are, respectively, ) . ( The matrix J =  −  is the Jacobian matrix evaluated at the DFE. ) .
The inverse matrix of  is The matrix FV −1 is called the next-generation matrix.The (, ) entry of  −1 indicates the expected number of new infections in compartment  produced by the infected individual originally introduced into compartment .The model reproduction number,  0 , which is defined as the spectral radius of FV −1 and denoted by (FV −1 ), is given by Here  0 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: ).
In ( 6),  0 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  0 as infection between mosquitoes and infectious humans or infection between mosquitoes and pseudorecovered humans. where ) .
01 in (8) represents the product of individuals generated by infectious human and infected mosquito, respectively, during their infectious time.
02 in ( 9) represents the product of individuals generated by pseudorecovered human and infected mosquito, respectively, during their infectious time.

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 Δ.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 Δ [25,26].
where  is positive and represents the maximum size of the populations.
If a disease emerges from one infectious group with  0 > 1 and if  infective agents are introduced into a wholly susceptible population, then the probability of a major disease outbreak is approximated by 1 − (1/ 0 )  while the probability of disease extinction is approximately (1/ 0 )  [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  0 greater than one.
For CTDSGWbp models, the transition from one state to a new state may occur at any time .If the process begins in  1.

The Branching Process Approximation.
We use branching process to analyze the malaria dynamics near the diseasefree equilibrium.Since infectious human, pseudorecovered humans, and infectious mosquitoes are the only sources of infection, we apply the multitype branching process in the three variables  ℎ (),  ℎ (), and   ().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 where   ( 1 ,  2 , . . .,   ) = prob( 1 =  1 ,  2 =  2 , . . .,   =   ) is the probability that one infected individual of type  gives birth to   individuals of type ; 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 is  ℎ (0) = 1,  ℎ (0) = 0, and   (0) = 0.
Similarly, the offspring pgf for  ℎ is given by where is the probability that one pseudorecovered human through infection produces an infectious human  1 or another pseudorecovered human  2 or an infectious mosquito  3 .Lastly, the offspring pgf for   is given by where is the probability that one infectious mosquito through infection produces an infections human  1 or a pseudorecovered human  2 or another infectious mosquito  3 .
The power to which   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 ( 0  ).The offspring pgfs for  ℎ ,  ℎ , and   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 for  ℎ ,  ℎ , and   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,  ℎ (0) ≈  and   (0) ≈ .
From ( 12), the offspring pgf for infectious human, given  ℎ (0) = 1,  ℎ (0) = 0, and   (0) = 0, is given by For  1 , 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 " 11 = 0,  12 = 01, and  13 = 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 " 11 = 1,  12 = 0, and  13 = 1."Note the term  1  3 in ( 15) means one infectious human generates one infectious mosquito ( 3 raised to power one) and remains infectious ( 1 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  ℎ (0) = 0,  ℎ (0) = 1, and   (0) = 0 is given by For  2 , 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 " 21 = 0,  22 = 1, and  23 = 1."Note the term  2  3 in ( 16) means one recovered human generates one infectious mosquito ( 3 raised to power one) and remains infectious ( 2 raised to power one).For one infectious mosquito, there are only two events: infection of a susceptible human or death of the mosquito.
The expectation matrix of the offspring pgfs, evaluated at (1, 1, 1), is given by The entries  11 ,  21 , and  31 represent the expected number of infectious humans, pseudorecovered humans, and infectious mosquitoes, respectively, produced by one infectious human.Similarly, the entries  12 ,  22 , and  32 represent the expected number of infectious humans, pseudorecovered humans, and infectious mosquitoes, respectively, produced by one pseudorecovered human.Lastly, the entries  13 ,  23 , and  33 represent the expected number of infectious humans, pseudorecovered humans, and infectious mosquitoes, respectively, produced by one infectious mosquito.
The matrix C = W(M − ) is the Jacobian matrix evaluated for the stochastic version. where ) is a diagonal matrix and I is the identity matrix.From ( 3) and ( 20), we show that ) . ( The spectral radius of matrix M obtained by finding the eigenvalues of matrix M is given by Taking the first expression of  0 = (M) in ( 24), we have ) . ( This gives the probability of malaria transmission by either infectious human or by infectious mosquito.From (23), The probability of disease extinction is one if (M) < 1.
Hence from (24) we have Expanding and simplifying the inequality, we have which reduces to 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
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/ 01 ).This implies that if there is successful contact, then the probability of successful transmission from infectious human to susceptible mosquito is 1 − 1/ 01 .The expression for  3 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/ 01 ).
There are some other important relationships; if disease transmission is by infectious mosquitoes as well as infection of susceptible mosquitoes by infectious humans, then 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/ 01 .If disease transmission is by infectious mosquitoes as well as infection of susceptible mosquitoes by pseudorecovered humans, then 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/ 02 .
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.Using  1 ,  2 , and  3 in ( 29)-(31) and assuming initial individuals for infectious humans, pseudorecovered humans and infected mosquitoes are  ℎ (0) = ℎ 0 ,  ℎ (0) =  0 , and   (0) =  0 , respectively.Then the probability of malaria clearance is given by ( +  1 +  1 ) ) )

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.

Effects of Model
Parameters on  0 .Using parameter values in Table 2, we identified how different input parameters affect the reproduction number  0 of our model as shown in Figure 3.
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  0 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  0 hence reducing the chances of disease outbreak.

Graphical Representation of Parameters
Effects on  0 .Figure 4 shows the effects of relapse and recover rates on the basic reproduction number.From the simulations, we find that  0 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.

Probability of Disease Extinction.
Using parameter values in Table 2, we compute numerically the probability of disease extinction  0 and of an outbreak 1 −  0 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 Table 3.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.
A 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.
Table 3 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.2, we numerically simulate the behavior of model (1).Initial conditions are  ℎ (0) = 99,  ℎ (0) = 1,  ℎ (0) = 0,   (0) = 999, and   (0) = 1.5, the analysis shows that when  0 < 1 and  0 < 1, then the probability of disease extinction is  0 = 0.9476 ≃ 1, although this agrees with (30), which points out that the probability of disease elimination in the CTDSGWbp model is one iff

Numerical Simulation When 𝑅
There is still a small probability 1 −  0 = 0.0524 of disease outbreak. 0 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  0 can be flawed and disease can persist with  0 < 1 depending on the kind of disease being modeled.

Numerical Simulation
When  0 > 1.Using parameter values in Table 2, we numerically simulate the behavior of model ( 1) when  0 > 1.Initial conditions are  ℎ (0) = 99,  ℎ (0) = 1,  ℎ (0) = 0,   (0) = 999, and   (0) = 1.The analysis from Figure 6 suggests that when  0 > 1, there is still some probability of disease extinction ( 0 = 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  0 > 1 and hence the  0 threshold should not be the only parameter to consider in quantifying the spread of a disease.

Effects of Relapse on Human
Population. Figure 7 shows changing effects of relapse on the human populations.
To 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 Figure 7.

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 extinction  0 and showed the relationship that exists between  0 and  0 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 when  0 < 1 and  0 < 1.However, the predictions by these models are different when  0 > 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  0 > 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

Figure 1 :
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].

Figure 4 :
Figure 4: Effects of relapse and recovery rate on  0 .

Table 3 :
Probability of disease extinction  0 and of an outbreak 1 −  0 for the malaria model.