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.
1. 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 to 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–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 threshold 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 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, , 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 Figure 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 , 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.
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 , 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.
2.4. Computation of Using the Next-Generation Matrix Approach
The basic reproduction number, , 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 . The original nonlinear system of ODEs including these compartments can be written as , where and represent new infections and transfer between compartments, respectively [22–24].
The Jacobian matrices of and at the disease-free equilibrium are, respectively, The matrix is the Jacobian matrix evaluated at the DFE. The inverse matrix of isThe matrix is called the next-generation matrix. The entry of indicates the expected number of new infections in compartment produced by the infected individual originally introduced into compartment . The model reproduction number, , which is defined as the spectral radius of and denoted by , is given byHere 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), is associated with disease transmission by infected mosquitoes as well as infection of susceptible mosquitoes by infected humans. The term represents the number of infected humans generated by infectious mosquito in its life span. The term represents the number of infected mosquitoes generated by infectious human during the infectious period of the individual while the term 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 as infection between mosquitoes and infectious humans or infection between mosquitoes and pseudorecovered humans. where in (8) represents the product of individuals generated by infectious human and infected mosquito, respectively, during their infectious time.
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 and the disease persists if . Hence , iff and .
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 . 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].
3.1. Model Formulation
Let time be continuous, , and let , , , , and represent discrete random variables for the number of susceptible humans, infectious humans, recovered humans, susceptible mosquitoes and infectious mosquitoes, respectively, with finite space.where is positive and represents the maximum size of the populations.
If a disease emerges from one infectious group with and if infective agents are introduced into a wholly susceptible population, then the probability of a major disease outbreak is approximated by while the probability of disease extinction is approximately [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 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 a state , after a random time period , it transits to a new state . The process remains in state for a random time after which it moves through to the new state , with [27]. This process continues throughout the model. The state transitions and rates for the stochastic model are presented in Table 1.
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 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 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 , , and .
The offspring pgf for is given bywhere is the probability that one infectious human through infection produces another infectious human or a pseudorecovered human or an infectious mosquito .
Similarly, the offspring pgf for is given bywhere is the probability that one pseudorecovered human through infection produces an infectious human or another pseudorecovered human or an infectious mosquito .
Lastly, the offspring pgf for is given bywhere is the probability that one infectious mosquito through infection produces an infections human or a pseudorecovered human or another infectious mosquito .
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 ().
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, and .
From (12), the offspring pgf for infectious human, given , , and , is given by For , one infectious human dies or is treated with probability ; this means an infectious human dies before infecting a susceptible mosquito. The term represents the probability that the infectious human is treated and moves to the pseudorecovered class, this results in “, , and , 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 . This means an infectious human infects a susceptible mosquito and remains infectious, which results in “, , and .” Note the term in (15) means one infectious human generates one infectious mosquito ( raised to power one) and remains infectious ( 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 , , and is given by For , one pseudorecovered host dies or relapses or is fully treated with probability or infects a mosquito with probability . This means a pseudorecovered human infects a susceptible mosquito and remains infectious, which results in “, , and .” Note the term in (16) means one recovered human generates one infectious mosquito ( raised to power one) and remains infectious ( 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 , , and is given by For , one infectious mosquito dies with probability or infects a human with probability . This means an infectious mosquito infects a susceptible human and remains infectious, which results in “, , and .” Note the term in (17) means one infectious mosquito generates one infectious human ( raised to power one) and remains infectious ( raised to power one).
3.3. The Relationship between and the Stochastic Threshold
The offspring pgfs, evaluated at , gives the expectation matrix with elements . Below are the offspring pgfs evaluated at . The expectation matrix of the offspring pgfs, evaluated at , is given by The entries , , and represent the expected number of infectious humans, pseudorecovered humans, and infectious mosquitoes, respectively, produced by one infectious human. Similarly, the entries , , and represent the expected number of infectious humans, pseudorecovered humans, and infectious mosquitoes, respectively, produced by one pseudorecovered human. Lastly, the entries , , and represent the expected number of infectious humans, pseudorecovered humans, and infectious mosquitoes, respectively, produced by one infectious mosquito.
The matrix is the Jacobian matrix evaluated for the stochastic version.where is a diagonal matrix and is the identity matrix.
From (3) and (20), we show that The spectral radius of matrix obtained by finding the eigenvalues of matrix is given by Taking the first expression of 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 . Hence from (24) we have Expanding and simplifying the inequality, we havewhich 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
3.4. Deriving Probability of Disease Extinction Using Branching Process Approximation
To find the probability of extinction (no outbreak), we compute the fixed points of the system of the offspring pgfs for the three infectious stages; that is, we solve , , and . The solutions of these systems are and ; see [19]. Equating in (15) to , then letting and , and solving for , we have Equating in (16) to , letting and , and solving for , we have Equating in (17) to , then letting and , and solving for , we have The expression for in (29) has a biological interpretation. Beginning from one infectious human, there is no outbreak if the infectious human recovers or dies with probability or if there is no successful transmission to a susceptible mosquito with probability . This implies that if there is successful contact, then the probability of successful transmission from infectious human to susceptible mosquito is .
The expression for in (31) has a biological interpretation. Beginning from one infectious mosquito, there is no outbreak if the infectious mosquito dies with probability or if there is no successful transmission to a susceptible human with probability .
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 .
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 .
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 , , and in (29)–(31) and assuming initial individuals for infectious humans, pseudorecovered humans and infected mosquitoes are , , and , respectively. Then the probability of malaria clearance is given by
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.
4.1. Effects of Model Parameters on
Using parameter values in Table 2, we identified how different input parameters affect the reproduction number of our model as shown in Figure 3.
From Figure 3, we see that increase in , , , and will decrease . Also decrease in , , , and will decrease .
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 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 ). Increasing would decrease hence reducing the chances of disease outbreak.
4.1.1. Graphical Representation of Parameters Effects on
Figure 4 shows the effects of relapse and recover rates on the basic reproduction number. From the simulations, we find that 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.
(a) Effects of relapse on
(b) Effects of treatment rate on
(c) Effects of recovery rate on
4.2. Probability of Disease Extinction
Using parameter values in Table 2, we compute numerically the probability of disease extinction and of an outbreak 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.
4.3. Numerical Simulation of Malaria Model
Using parameter values in Table 2, we numerically simulate the behavior of model (1). Initial conditions are , , , , and .
4.3.1. Numerical Simulation When
From Figure 5, the analysis shows that when and , then the probability of disease extinction is , although this agrees with (30), which points out that the probability of disease elimination in the CTDSGWbp model is one iff There is still a small probability of disease outbreak. 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 can be flawed and disease can persist with depending on the kind of disease being modeled.
(a) Human population dynamics when and
(b) Mosquito population dynamics when and
4.3.2. Numerical Simulation When
Using parameter values in Table 2, we numerically simulate the behavior of model (1) when . Initial conditions are , , , , and .
The analysis from Figure 6 suggests that when , there is still some probability of disease extinction (). 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 and hence the threshold should not be the only parameter to consider in quantifying the spread of a disease.
(a) Human population dynamics when
(b) Mosquito population dynamics when
4.3.3. Effects of Relapse on Human Population
Figure 7 shows changing effects of relapse on the human populations.
(a) Human population dynamics when relapse rate is increased;
(b) Human population dynamics when relapse rate is reduced;
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.
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 extinction and showed the relationship that exists between and 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 and . However, the predictions by these models are different when . 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 . 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 least by 2030 and (ii) reducing malaria mortality rates by at least 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.
Notations
Description of Variables for the Malaria Transmission Model| Susceptible humans | |
| Infectious humans | |
| Pseudorecovered humans | |
| Susceptible mosquito | |
| Infectious mosquito. |
| Natural birth rate of humans | |
| Natural death rate of humans | |
| Disease induced death rate of humans | |
| Transmission rate from an infectious mosquito to a susceptible human | |
| Recovery rate | |
| Relapse rate | |
| Treatment rate | |
| Natural birth rate of mosquitoes | |
| Natural death rate of mosquitoes | |
| Parasite induced death rate of mosquitoes | |
| Transmission rate from an infectious human to a susceptible mosquito | |
| Transmission rate from a pseudorecovered human to a susceptible mosquito | |
| The total size of human population | |
| The total size of mosquito population. |
Conflicts of Interest
All authors declared no conflicts of interest regarding the publication of this paper.
Authors’ Contributions
Rachel Waema Mbogo produced the first draft of the manuscript. All authors critically reviewed the paper and approved the final version.