The dynamics of multiple species and strains of malaria

This paper presents a deterministic SIS model for the transmission dynamics of malaria, a life-threatening disease transmitted by mosquitos. Four species of the parasite genus Plasmodium are known to cause human malaria. Some species of the parasite have evolved into strains that are resistant to treatment. Although proportions of Plasmodium species vary considerably between geographic regions, multiple species and strains do coexist within some communities. The mathematical model derived here includes all available species and strains for a given community. The model has a disease-free equilibrium, which is a global attractor when the reproduction number of each species or strain is less than one. The model possesses quasi-endemic equilibria; local asymptotic stability is established for two species, and numerical simulations suggest that the species or strain with the highest reproduction number exhibits competitive exclusion.


Introduction
Malaria is an infectious disease of humans caused by a parasite, Plasmodium, which infects red blood cells. It is widespread in tropical and subtropical regions, including much of sub-Saharan Africa, Asia and the Americas. According to the World Health Organization, about 250 million clinical cases of malaria occur globally each year resulting in more than one million deaths, with sub-Saharan Africa being the worst affected region (WHO, 2013). The four common species of the parasite are: Plasmodium falciparum, Plasmodium vivax, Plasmodium ovale and Plasmodium malariae. P. falciparum is the most virulent and potentially lethal species to humans. Infections resulting from the species P. vivax, P. malariae and P. ovale are generally less serious and are usually not life-threatening. A relatively new species, Plasmodium knowlesi, causes malaria in macaques monkeys but can also infect humans. Given that more than one species can be found in some locations, it is not uncommon to be infected with more than one species of Plasmodium at the same time. Resistance to antimalarial drugs has been one of the main obstacles in the fight against malaria. Drug resistance has so far been documented in P. falciparum, P. vivax and P. malariae (WHO, 2000(WHO, -2010. Mathematical modelling is a powerful tool that has been instrumental in understanding and combating malaria infection (McKenzie, 2000;McKenzie & Samba, 2004). The classical works of Ross (1910) and MacDonald (1955) have led to various forms of compartmental models for malaria infection. Since the 1970s, malaria models have been expanded and brought closer to data (Dietz, Molineaux, & Thomas, 1974;Gosoniu & Vounatsou, 2011). There are several models of malaria, each of which describes certain aspects of the disease (Teboh-Ewungkem, Ngwa, & Ngonghala, 2013). Some models consider interacting populations of humans and mosquitos (Dietz, 1988;Mandal, Sinha, & Sarkar, 2013). Others are solely focused on the dynamics within the host (Mitchell & Carr, 2010;Molineaux & Dietz, 1999) or vector (Ngwa, Niger, & Gumel, 2010;Ngwa, Wankah, Fomboh-Nforba, Ngonghala, & Teboh-Ewungkem, 2014) and some do consider drug resistance (Agusto, 2014;Aneke, 2002). The few models devoted to multiple species have been limited to two species (Xiao & Zou, 2013) P. falciparum and P. vivax (Pongsumpun & Tang, 2009) or two strains of a single species (Esteva, Gumel, & de Leon, 2009). This paper presents a deterministic model for monitoring the transmission dynamics of all possible circulating malaria parasites in a community. This study is motivated by the fact that multiple species and strains of Plasmodium parasite do coexist in some communities (Gnémé et al., 2013;WHO, 2000WHO, -2010Xiao & Zou, 2013). Our primary objective is to construct a model with n infectious compartments that includes different species and strains. For the purpose of simplicity, the model will treat each new strain as a separate species. The model derived then can be used to study (i) multiple species of plasmodium, (ii) multiple strains of the same species and (iii) a combination of these two.

Mathematical model
We assume that there is no immunity to the disease and every new strain is treated as another species. The schematics of the SIS model that takes into account all n species (that includes strains) of malaria in a local community is given in Figure 1. All the parameters are assumed to be positive: S h and S m are the susceptible human and mosquito populations, respectively; θ and ν are the natural birth/death rates for humans and mosquitos, respectively. We assume that all new offspring are susceptible in all populations. I hk and I mk are the infectious human and mosquito populations respectively for each species where k = 1, . . . , n; γ k and β k are the human recovery and death rates for species k. The average biting rate of the mosquitos is a, ρ j is the probability that a mosquito bite will lead to an infection in a human and σ j is the probability that mosquitos will be infected as a result of taking blood meals from infectious humans. The total human and mosquito populations are N h and N m , respectively.
The following equations, which summarize the model, are based on the assumptions provided in the schematic diagram, Figure 1: for k = 1, . . . , n. We remark here that in those communities (for example the Americas) where malaria-induced death is very small, the β k parameters can be set to zero.
In the remainder of this paper, while focusing our analysis on the infective equations of the system (1), we will assume that the total population sizes for both humans (N h ) and mosquitos (N m ) are constant and normalized to one. By defining the new variables u h = S h /N h , u m = S m /N m , x k = I hk /N h , and y k = I mk /N m , for k = 1, . . . , n, the infective equations of the system becomeẋ where c k = γ k + β k + θ, α k = (aρ k N m )/(c k N h ), δ k = aσ k /ν and for k = 1, . . . , n.
The basic reproduction number R 0 , is an established epidemiologic statistic that provides an indication of the average number of secondary infections produced when one infected individual is introduced into a naive host population. Let R 0k denote the reproduction number of secondary infections that we expect to be produced by introducing a single human infected with species (or strain) k. This will lead to the infection of α k cases in mosquitos, each of which will lead to δ k cases in humans. The reproduction number for each species (or strain) is therefore R 0k = α k δ k , which is the product of α k and δ k . The basic reproduction number of the model is then given by the sum R 0 = n j=1 R 0j .

Stability analysis
In this section, we study the long term behaviour of the disease-free equilibrium and the endemic equilibria of the proposed model. Suppose that no deaths result from the infections (i.e. β k = 0 for all k). The region D = {(x 1 , . . . , x n , y 1 , . . . , y n ) ∈ R 2n + : n j=1 x j ≤ 1, n j=1 y j ≤ 1} is positively invariant and contains all solutions of the system (2).
Using that u h = 1 − n j=1 x j and u m = 1 − n j=1 y j , system (2) becomeṡ for k = 1, . . . , n. The equilibria of system (3) satisfy the equations Multiply the x k -y k pairs for every k: For any k value, x k = 0 if and only if y k = 0 (this is immediate from (4) and (5)); if this is the case for all k values, we obtain the disease-free equilibrium. We remark that the reproduction number for species k is R 0k = α k δ k . If we assume that R 0k = R 0l for k = l, then by the previous equation, the only other equilibria are given by the quasi-endemic states because if we assume that x * k = 0 and y * k = 0 for some k, then from (6) we get which implies the uniqueness of a non-zero pair (x * k , y * k ) in case of R 0k = R 0l for k = l. We called the equilibrium quasi-endemic because only one of the species is present. We observe that the quasi-endemic equilibria are biologically relevant only when R 0k > 1. Theorem 1: The disease-free equilibrium of system (3) is locally asymptotically stable if R 0k < 1 for all k = 1, . . . , n and unstable if there exists an R 0k > 1 for some k.
Proof: The Jacobian matrix, J 0 , of the system (3) at the origin is given by the sparse matrix which upon manipulation yields The matrix J 0 therefore has n pairs of eigenvalues given by Clearly, if R 0k < 1 for all k = 1, . . . , n, all n pairs of eigenvalues are negative and therefore the disease-free equilibrium is locally asymptotically stable. If R 0k > 1 for some k, then there is at least one positive eigenvalue and the disease-free equilibrium is unstable.
The following result is stronger than Theorem 1 and provides a sufficient condition for the eradication of the k th species or strain. Theorem 2: If R 0k < 1 for some k, then x k (t) → 0 and y k (t) → 0 as t → ∞. Proof: Since S h ≤ N h and S m ≤ N m , then u h ≤ 1 and u m ≤ 1. For some k, system (2) leads to the differential inequalityẇ . The eigenvalues of the matrix M k are given by the pair If R 0k < 1 for a given k, then λ k,1 < 0 and λ k,2 < 0. Let v 1 and v 2 be eigenvectors associated with the eigenvalues λ k,1 and λ k,2 respectively, then the systeṁ , has a solution of the formw k (t) = v 1 e λ k,1 t + v 2 e λ k,2 t . Clearly, it can be seen thatw k (t) → 0 as t → ∞. Observing that x k (t) ≤x k (t) and y k (t) ≤ỹ k (t) concludes the proof. The immediate consequence of Theorem 2 is the following result. Theorem 3: If α k < 1 and δ k < 1 (note that in this case R 0k < 1) for all k = 1, . . . , n, then the disease-free equilibrium of the system (3) is globally asymptotically stable.
We next turn our attention to the quasi-endemic equilibrium. Due to the complexities involved in the analysis, we will only present the case of two species. Theorem 4: If n = 2, the quasi-endemic equilibrium (x * 1 , 0; y * 1 , 0) of the system (3) is locally asymptotically stable if 1 < R 01 and R 02 < R 01 . Proof: Jacobian of the system (3) at the equilibrium (x * 1 , 0; y * 1 , 0) is Substituting the values of x * 1 and y * 1 , after simplification this yields The eigenvalues of the preceding matrix are given by the conjugate pairs The first pair is negative when This gives which is the same as α 2 δ 2 < α 1 δ 1 , i.e. R 02 < R 01 . The second pair is negative when This is the same as

Results and discussion
In this section, we present numerical simulations to illustrate some of the theoretical results established in the previous section for the model (3). We apply the model to two different cases in a given community. The first case simulates the transmission dynamics of two species, P. falciparum and P. vivax. The second case simulates the transmission dynamics of two strains of P. falciparum. Note that n = 2 in both cases. In each case, we will study the effect of varying the density of vectors in relation to human hosts and varying the vector biting rate. We assumed that there are no disease-induced deaths in all our simulations. The values of the parameters that we used in our simulations were taken from related works (Aron & May, 1982;Esteva, Gumel, & de Leon, 2009;Pongsumpun & Tang, 2009) and others were chosen to explore the behaviour of the model (see Table 1).
The results reported here are based on the initial conditions: x 1 (0) = .1 and x 2 (0) = .1. We ran the simulation for different sets of initial conditions chosen by giving a small perturbation to the disease-free equilibrium and the qualitative form of the final finiteamplitude steady-state solutions was the same for all of them.
The initial and long-term behaviour of the model for infectious populations of the species, P. falciparum and P. vivax, are shown in Figure 2. The results in Figure 2 show that the disease dies out in both species. For both population ratios, the calculated reproductive number for each species was less than one and thus established that the disease-free equilibrium is a global attractor.
By varying the average biting rate of mosquitos, the results in Figure 3 show that the disease persists in P. falciparum, but dies out in P. vivax. Let R 0f and R 0v denote the reproduction numbers of P. falciparum and P. vivax, respectively. We found that R 0f > 1 and R 0v < 1 for the results in Figure 3(a), and R 0f > R 0v > 1 for the simulations in Figure 3(b). The model clearly predicts competitive exclusion when R 0f = R 0v , where the species with the higher reproduction number eventually displaces the other. This is why P. vivax does not survive irrespective of the fact that R 0v > 1. It is clear that the model has a quasi-endemic equilibrium that is a global attractor.
The simulations in Figure 4 exhibit the general behaviour of two strains of the species, P. falciparum (one strain resistant to treatment and the other not). For this example, let R 0r and R 0n denote the reproduction numbers of the resistant and non-resistant strains, respectively. Figure 4(a) shows both strains dying out. This is expected since the calculated values of R 0r and R 0n are both less than one. However, when we varied the total mosquito to human population ratio from 5 to 10 while keeping the average biting rate constant, we observed that R 0r > 1 and R 0n < 1 as depicted in Figure 4(b). The outcomes of the simulations in the cases of the lower and higher mosquito to human population ratio clearly underscore the fact that controlling the vector population is a key aspect in the fight against malaria, and also validate the seasonal prevalence of malaria, although the present model does not incorporate seasonality.
An increase in the average biting rates of mosquitos from .25 to .5 yields R 0r > R 0n > 1 for both population ratios as illustrated in Figure 5. Similar to the previous case, the resistant  strain dominates the non-resistant strain because it has a higher reproduction number. Figure 5(b) shows that under favourable conditions (high vector population and biting rates), the non-resistant strain may persist for some time but is eventually displaced by the resistant strain. Here again, we see that the model possesses a quasi-endemic equilibrium which is a global attractor.

Conclusion
In this paper, we have developed a susceptible-infectious model for the dynamics and transmission of malaria which can be used to study multiple species (and/or strains) of the disease. The model is based on the assumption that there is no immunity to the disease. It is formulated mathematically as a system of differential equations for the proportions of susceptible and infectious individuals in all populations (vector and host). The existence of a disease-free equilibrium that is globally asymptotically stable under certain conditions is confirmed. Numerical simulations for two species/strains established that if the reproduction number of at least one species or strain is greater than one, then there is a quasi-endemic equilibrium, which is a global attractor. The model predicts that a species or strain will be eradicated if its reproduction number is less than one. The model further predicts that all species or strains may persist for some time if each has a reproduction number greater than one; however, the species or strain with the highest reproduction number will eventually displace the others.
There are various factors that we did not consider in this model, such as seasonality, age structure of both humans and mosquitos, incubation period and spatial distribution. Some of these are currently under investigation by the authors.

Disclosure statement
No potential conflict of interest was reported by the authors.