Research Article  Open Access
MultipleStrain Malaria Infection and Its Impacts on Plasmodium falciparum Resistance to Antimalarial Therapy: A Mathematical Modelling Perspective
Abstract
The emergence of parasite resistance to antimalarial drugs has contributed significantly to global human mortality and morbidity due to malaria infection. The impacts of multiplestrain malarial parasite infection have further generated a lot of scientific interest. In this paper, we demonstrate, using the epidemiological model, the effects of parasite resistance and competition between the strains on the dynamics and control of Plasmodium falciparum malaria. The analysed model has a trivial equilibrium point which is locally asymptotically stable when the parasite’s effective reproduction number is less than unity. Using contour plots, we observed that the efficacy of antimalarial drugs used, the rate of development of resistance, and the rate of infection by merozoites are the most important parameters in the multiplestrain P. falciparum infection and control model. Although the drugresistant strain is shown to be less fit, the presence of both strains in the human host has a huge impact on the cost and success of antimalarial treatment. To reduce the emergence of resistant strains, it is vital that only effective antimalarial drugs are administered to patients in hospitals, especially in malariaendemic regions. Our results emphasize the call for regular and strict surveillance on the use and distribution of antimalarial drugs in health facilities in malariaendemic countries.
1. Introduction
The emergence of parasite resistance [1–4] to antimalarial drugs has contributed significantly to human mortality and morbidity due to malaria infection, worldwide [5–7]. A global malaria control strategy of 1992 [8] that advocated for early diagnosis and prompt treatment has been heavily compromised by the emergence of parasite resistance to antimalarial drugs. The evolution of parasite resistance has been described in [9] as an example of a Darwinian evolution. Parasites undergo mutations in their genome in response to the drugtreated human host. These mutations reduce the rate of parasite elimination from the host and increase their survival chances [9]. The most extensively used antimalarial drugs against the deadly Plasmodium falciparum malaria are chloroquine (CQ) and sulfadoxinepyrimethamine (SP) [10, 11]. These drugs are cheap, easily available, and slowly eliminated from the human body [11]. However, the extensive use of CQ and SP has resulted in P. falciparum resistance. This has led to global increase in malaria cases and mortality [12]. In response, the World Health Organization (WHO) in 2006 recommended the use of artemisininbased combination therapies (ACTs) as a firstline treatment for uncomplicated P. falciparum malaria [13]. Resistance to ACTs which are currently the standard treatment for P. falciparum is likely to cause global health crisis especially in African regions where P. falciparum malaria is endemic [11].
The emergence of parasite resistance to malaria therapy dates back to the 19th century. Quinine (1963) was the firstline antimalarial drug against P. falciparum [14]. High mortality cases coupled with high parasite resistance led to the introduction of a second drug, chloroquine (CQ), in 1934 [15]. A decade later, CQ was considered the firstline antimalarial drug by several countries until 1957, when the first focus of P. falciparum resistance was detected along the ThaiCambodia border [16]. In Africa, P. falciparum resistance to CQ was first discovered among travelers from Kenya to Tanzania [17]. By 1983, CQ resistance had spread to Sudan, Uganda [18], Zambia [19], and Malawi [20]. Unlike Africa, CQ was replaced for the first time with sulfadoxinepyrimethamine (SP) as a firstline antimalarial drug in Thailand in 1967. Several other countries in Asia and South America followed thereafter [10]. Resistance to SP was, however, reported the same year [21] in the region. In 1988, CQ was replaced for the first time in Africa. KwaZuluNatal Province of South Africa replaced CQ with SP [22]. In 1993, the Malawian government changed the treatment policy from CQ to SP. Other African countries followed thereafter: Kenya, South Africa, and Botswana (in 1998); Cameroon and Tanzania (in 2001); and Zimbabwe (in 2000) [23]. The effectiveness of SP was equally undermined by resistance. Unlike CQ, P. falciparum resistance to SP was mainly attributed to the long halflife of the drug [24]. Confirmed resistance to the artemisinin derivatives was first reported in Cambodia and Mekong regions in 2008 [25].
To leverage on parasite resistance, cost of treatment, and burden of malaria infection to communities and governments, the WHO recommends the use of artemisininbased combination therapies (ACTs) as the first and secondline treatment drugs for uncomplicated P. falciparum malaria [25]. ACT is a combination of artemisinin derivatives and a partner monotherapy drug. Artemisinin derivatives include artemether, artesunate, and dihydroartemisinin. These derivatives reduce the parasite biomass within the first three days of therapy, while the partner drug, with longer halflife, eliminates the remaining parasites [26]. The WHO currently recommends five different ACTs: (1) artesunateamodiaquine (AS + AQ), (2) artesunatemefloquine (AS + MQ), (3) artesunate + sulfadoxinepyrimethamine (AS + SP), (4) artemetherlumefantrine (AMLM), and (5) dihydroartemisininpiperaquine (DHA + PPQ). Additionally, artesunatepyronaridine may be used in regions where ACT treatment response is low [26]. Access to ACT has been tremendous in the last 8 years, with a recorded increase of 122 million procured treatment courses for the period 2010–2016. However, resistance to currently used ACTs has important public health consequences, especially in the African region, where resistant P. falciparum is predominant.
Numerous crosssectional studies [27, 28] have revealed the possible impacts of multiple strains of P. falciparum on the development of resistance to ACTs. In [29] and citations therein, drugsensitive parasites are shown to strongly suppress the growth and transmission of drugresistant P. falciparum parasites. Although hightransmission settings such as subSaharan Africa account for about 90% of all global malaria deaths, resistance to antimalarial drugs has been shown to emerge from lowtransmission settings, such as Southeast Asia and South America [29]. Causes of parasite resistance to ACTs are diverse. Historical studies [30, 31] indicate that antimalarialresistant parasites could emerge from a handful of lineages. It is argued elsewhere [32, 33] that recombination during sexual reproduction in the mosquito vector could be responsible for the delayed appearance of multilocus resistance in hightransmission regions. Moreover, owing to repeated exposure for many years, individuals in hightransmission settings are likely to develop clinical immunity to malaria, leading to stronger selection for resistance [34]. Studies in [29] also support the hypothesis that inhost competition between drugsensitive and drugresistant parasites could inhibit the spread of resistance in hightransmission settings. Owing to their integral role in the recent success of global malaria control, the protection of efficacy of ACTs should be a global health priority [35].
Mathematical models of inhost malaria epidemiology and control constitute important tools in guiding strategies for malaria control [36, 37] and the associated financial planning [38]. While some researchers have focussed on probabilistic models [39, 40], others have investigated the effects of drug treatment and resistance development using dynamic models [41, 42]. A deterministic model by Esteva et al. [43] monitored the impact of drug resistance on the transmission dynamics of malaria in a human population. In [29], the impacts of withinhost parasite competition are shown to inhibit the spread of resistance [44, 45]. On the contrary, some models [39, 46] have suggested that withinhost competition is likely to speed up the spread of resistance in hightransmission settings due to a phenomenon called “competitive release.” In this paper, we provide theoretical insights using mathematical modelling of the impacts of multiplestrain infections on resistance, dynamics, and antimalarial control of P. falciparum malaria.
The rest of the paper is organized as follows: In Section 2, we formulate the withinhuman malaria model that has both the drugsensitive and drugresistant P. falciparum parasite strains subject to antimalarial therapy. In Section 3, we analyze the model based on epidemiological theorems. Withinhost competition between parasite strains and the effects of antimalarial drug efficacy on parasite clearance are discussed in Section 4. Sensitivity analysis and multiplestrain infection and its effects on resistance and malaria dynamics are demonstrated in Section 5. We conclude the paper in Section 6 by emphasizing the need for antimalarial therapy with the potential to eradicate multiplestrain infection due to P. falciparum.
2. Model Formulation
We present in this paper a deterministic model that describes the withinhumanhost competition and transmission dynamics of two strains of P. falciparum parasites during malaria infection. The compartmental model considers the coinfection and competition between the drugsensitive (dss) and the drugresistant (drs) P. falciparum strains in the presence of antimalarial therapy. The drs arise presumably from the dss. The rare mechanism here could possibly be due to single point mutation [47]. Both drs and dss initiate immune responses that follow densitydependent kinetics.
Our model is composed of eight compartments: susceptible/healthy/unparasitized erythrocytes (red blood cells) , parasitized/infected erythrocytes ( and ), merozoites ( and ), gametocytes ( and ), and immune cells . The healthy erythrocytes (RBCs) make up the resource for competition between the drugresistant and drugsensitive parasite strains. The infected red blood cells (IRBCs) and different erythrocytic parasite life cycles are categorized based on the strain of the infecting parasite. The merozoites are therefore categorized into drugsensitive and drugresistant strains, denoted by and , respectively. The merozoites invade the healthy erythrocytes during the erythrocytic stage, leading to formation of infected erythrocytes. The variable denotes the red blood cells (RBCs) infected with drugsensitive merozoites, whereas refers to the RBCs infected with drugresistant merozoites. Similarly, the variables and represent drugsensitive and drugresistant gametocytes, respectively. Owing to saturation in cell and parasite growth, we consider the nonlinear Michaelis–Mented–Monod function described in [48, 49] and used in [50–53] to model the reductive effects of the immune cells on the parasite and infectedcell populations.
The density of the healthy RBCs is increased at the rate per healthy RBC per unit time from the host’s bone marrow, and healthy RBCs die naturally at a rate . Following parasite invasion by free floating merozoites, the healthy erythrocytes get infected by both drugsensitive and drugresistant merozoite strains at the rates β and , respectively. The parameter (with ) accounts for the reduced fitness (infectiousness) of the resistant parasite strains in relation to the drugsensitive strains. The destruction of the healthy red blood cells is however limited by the adaptive immune cells W. This is represented by the term , where γ is a measure of the efficacy of the immune cells. The equation that governs the evolution of the healthy RBCs is hence given by
The parasitized erythrocytes are generated through mass action contact (invasion) between the susceptible healthy erythrocytes X and the blood floating merozoites ( and ). The merozoites subdivide mitotically, within the infected erythrocytes, into thousands of other merozoites, leading to cell burst and emergence of characteristic symptoms of malaria. Additionally, a single infected erythrocyte undergoes hemolysis at the rate to produce P secondary merozoites, sustaining the erythrocytic cycle. The drugsensitive IRBCs burst open to generate more drugsensitive merozoites or drugsensitive gametocytes at the rate . Similar dynamics are observed with the drugresistant IRBCs, where the drugresistant gametocytes are generated at the rate from IRBCs. Treatment with ACT is assumed to disfranchise the development of the merozoite within the infected erythrocyte. The druginfested erythrocytes are hence likely to die faster. This is represented by the term , where represents the antimalarialspecific treatment efficacy. In this paper and for purposes of illustration and simulations, corresponds to the efficacy of artemetherlumefantrine (AL), which is the recommended firstline antimalarial ACT drug for P. falciparum infection in Kenya. We assume that no treatment is available for erythrocytes infected with the resistant parasite strains. The time rate of change for and takes the following form:
The drugresistant merozoites and the drugresistant gametocytes die naturally at the rates and , respectively. It is further assumed that drugsensitive merozoites and gametocytes may develop into drugresistant merozoites and gametocytes at the rates and , respectively. The cost of resistance associated with AL is represented by the parameter . Parasite resistance to antimalarial drugs exacerbates the erythrocytic cycle and increases the cost of treatment [54, 55]. The higher the resistance to antimalarial therapy, the higher the density of malarial parasites in blood. We therefore model this decline in drug effectiveness by rescaling the density of merozoites produced per bursting parasitized erythrocyte P by the factor , where implies no resistance; that is, the ACT is highly effective in eradicating the parasites. If corresponds to maximum resistance, the used ACT drug is least effective in treating the infection. The converse of these descriptions applies to the drugresistant P. falciparum parasite strains. The equations that govern the rate of change of the infected red blood cells and the merozoites take the following form:
Antimalarial therapy increases the rate of elimination of drugsensitive merozoites and gametocytes. This is represented by the nonnegative enhancement parameters ζ and η, respectively.
Although the innate immunity is faster, it is often limited by the on and off rates in its response to invading pathogens [56, 57]. The adaptive immunity, on the contrary, is very slower at the beginning but lasts long enough to ensure no parasite growth in subsequent infections [27]. We assume an immune system that is independent of the invading parasite strain. For purposes of simplicity, we only consider the adaptive immune system, which is mainly composed of the cells [58]. We adopt the assumption that the background recruitment of immune cells is constant (at the rate ). Additionally, the production of the immune cells is assumed to be boosted by the infective and infected cells , , and at constant rates , , and , respectively. Circulating gametocytes, infective merozoites, and infected erythrocytes are removed phagocytotically by the immune cells at the rates , , and , respectively. The immune cells also get depleted through natural death at the rate . The equation for the immune cells takes the following form:
Following invasion by the merozoites, the IRBCs either produce merozoites or differentiate into gametocytes upon bursting. The total erythrocyte population at any time t, denoted by , is therefore given by
Similarly, the sum total of P. falciparum parasites, denoted by , within the host at any time t is described by the following equation:
The above dynamics can be represented by the schematic diagram in Figure 1. The list of model variables and model parameters is provided in Tables 1 and 2, respectively.


2.1. Model Equations
Based on the above model descriptions and schematic diagram shown in Figure 1, the model in this paper consists of the following nonlinear system of ordinary differential equations:subject to the following initial conditions:
3. Model Analysis
3.1. Positivity and Uniqueness of Solutions
The consonance between a formulated epidemiological model and its biological reality is key to its usefulness. Given that all the model parameters and variables are nonnegative, it is only sound that the model solutions be nonnegative at any future time within a given biological space.
Theorem 1. The region with solutions of system (7)–(14) is positively invariant under the flow induced by system (7)–(14).
Proof. We need to show that every trajectory from the region will always remain within it. By contradiction, assume (where refers to time) in the interval , such that , but for , , and , , , and for . Notice that, at , is declining from the original zero value. If such an X exists, then it should satisfy the differential equation (7). That is,We arrive at a contradiction, i.e., . This shows the nonexistence of such . This argument can be extended to all the remaining seven variables . The process of verification is however simpler. We can follow the steps as presented in [59, 60]. Let the total erythrocyte population evolve according to the following formulation:where . Similarly, the total density of malarial parasites is described bywhere .
The solutions of equations (14), (17), and (18) are, respectively, given aswhereHere, and represent the initial total populations of erythrocytes and malarial parasites, respectively. We observe that all the solutions of equations (14), (17), and (18) remain nonnegative for all future time, . Moreover, the total populations are bounded: , and . Thus, all the state variables of model system (7)–(14) and all their corresponding solutions are nonnegative and bounded in the phase space φ, whereIt is obvious that φ is twice continuously differentiable function. That is, . This is because its components , are rational functions of state variables that are also continuously differentiable functions. We conclude that the domain φ is positively invariant. It is therefore feasible and biological meaningful to study model system (7)–(14).
Theorem 2. The model system (7)–(14) has a unique solution.
Proof. Let so that and as presented in system (7)–(14). Similarly, let be a vector defined in . The model system (7)–(14) can hence be written aswhere denotes a column vector of state variables and represents the righthand side (RHS) of system (7)–(14). The result is as follows.
Lemma 1. The function g is continuously differentiable in x.
Proof. All the terms in g are either linear polynomials or rational functions of nonvanishing polynomials. Since the state variables are all continuously differentiable functions of t, all the elements of vector g are continuously differentiable. Moreover, let . By the mean value theorem,where denotes the mean value point and the directional derivative of the function g at m. However,where is the coordinate unit in . We can clearly see that all the partial derivatives of g are bounded and that there exists a nonnegative U such thatTherefore, there exists such thatThis shows that the function g is Lipschitz continuous. Since g is Lipschitz continuous, model system (7)–(14) has a unique solution by the uniqueness theorem of Picard [61].
3.2. Stability Analysis of the ParasiteFree Equilibrium Point (PFE)
The inhost malaria dynamics are investigated by studying the behaviour of the model at different model equilibrium points. Knowledge on model equilibrium points is useful in deriving parameters that drive the infection to different stability points. The model system (7)–(14) has a parasitefree equilibrium point given by
Using the nextgeneration operator method by van den Driessche and Watmough [62] and matrix notations therein, we obtain a nonsingular matrix Q showing the terms of transitions from one compartment to the other and a nonnegative matrix F of new infection terms as follows:where , , , , and , .
The effective reproduction number of model system (7)–(14) associated with the parasitefree equilibrium is the spectral radius of the nextgeneration matrix , where
It follows thatwhere
From equation (31), it is evident that, in a multiplestrain P. falciparum malaria infection, the progression of the disease depends on the reproduction number of different parasite strains. If the threshold quantity , the drugsensitive parasite strains will dominate the drugresistant strains and hence the driver of the infection. To manage the infection in this case, the patient should be given antimalarials that can eradicate the drugsensitive parasites. Conversely, if , the infection is mainly driven by the drugresistant parasite strains. In this scenario, the used antimalarial drugs should be highly efficacious and effective enough to kill both the drugresistant and drugsensitive parasite strains in the blood of the human host. This result is quite instrumental in improving antimalarial therapy for P. falciparum infections. The best antimalarials should be sufficient enough to eradicate both parasite strains within the human host.
Based on Theorem 2 in [63], we have the following lemma.
Lemma 2. The parasitefree equilibrium point is locally asymptotically stable if and unstable otherwise.
The Jacobian matrix associated with the inhost model system (7)–(14) at is given bywhere the terms are as defined in (30). It is clear from matrix (33) that the first four eigenvalues are (from column 1), (from column 8), (from column 7), and (from column 6). They are all negative. The remaining four eigenvalues are obtained from the roots of the following quartic equation:where
Due to complexity in the coefficients of the polynomial (34), we shall rely on the Routh–Hurwitz stability criterion [64], which provides sufficient condition for the existence of the roots of the given polynomial on the left half of the plane.
Definition 1. The solutions of the quartic equation (34) are negative or have negative real parts provided that the determinants of all Hurwitz matrices are positive [64].
Based on the Routh–Hurwitz criterion, the system of inequalities that describe the stability region is presented as follows:(i)(ii)(iii)(iv)
From (35), it is clear that . Upon simplifying in (36), we obtainwhere
Thus,
Similarly, the expression for can be rewritten as follows:
Lastly, upon simplifying equation (37), we obtain
Since all the coefficients of the quartic equation (34) are nonnegative, all its roots are therefore negative or have negative real parts. Hence, the Jacobian matrix (33) has negative eigenvalues or eigenvalues with negative real parts if and only if the effective reproduction number is less than unity. Equilibrium point is therefore locally asymptotically stable when (when both and ). This implies that an effective antimalarial drug would cure the costrain infected human host, provided that the drug reduces the effective reproduction number to less than 1.
Lemma 2 shows that P. falciparum malaria can be eradicated/controlled within the human host if the initial parasite and cell populations are within the basin of attraction of the trivial equilibrium point . To be certain to eradicate/control the infection irrespective of the initial parasite and cell populations, we need to prove the global stability of the parasitefree equilibrium point. This is presented in the following section.
3.3. Global Asymptotic Stability Analysis of the ParasiteFree Equilibrium Point
Following the work by Kamgong and Sallet [65], we begin by rewriting system (7)–(14) in a pseudotriangular form:where is a vector representing the densities of noninfective population groups (unparasitized erythrocytes and immune cells) and represents the densities of infected/infective groups (infective P. falciparum parasites and/or infected host cells) that are responsible for disease transmissions. For purposes of clarity and simplicity to the reader, we shall represent with and with in . We assume the existence of a parasitefree equilibrium in φ: . Thus,
We analyze system (43) based on the assumption that it is positively invariant and dissipative in φ. Moreover, the subsystem is globally asymptotically stable at on the projection of φ on . This implies that whenever there are no infective malarial parasites, all cell populations will settle at the parasitefree equilibrium point . Finally, in (43) is a Metzler matrix that is irreducible for any . We assume adequate interactions between and among different parasites and cell compartments in the model.
The matrices and are easily computed from subsystem in (43) so that we have
We can easily see that the eigenvalues of matrix are both real and negative (, ). This shows that the subsystem is globally asymptotically stable at the trivial equilibrium . Additionally, from subsystem , we obtain the following matrix:
Notice that all the offdiagonal entries of are nonnegative (equal to or greater than zero), showing that is a Metzler matrix. To show the global stability of the parasitefree equilibrium , we need to show that the square matrix in (46) is Metzler stable. We therefore need to prove the following lemma.
Lemma 3. Let K be a square Metzler matrix that is block decomposed:where and are square matrices. The matrix K is Metzler stable if and only if and are Metzler stable.
Proof. The matrix K in Lemma 3 refers to in our case. We therefore letResults from analytical computations based on Maple software givewhere , , and .
From equation (50), it is evident that all the diagonal elements of matrix are negative and the rest of the elements in the matrix are nonnegative. This shows that matrix is Metzler stable, and the parasitefree equilibrium point is globally asymptotically stable in the biologically feasible region φ of model system (7)–(14). Epidemiologically, the above result implies that when there is no malaria infection, different cell populations under consideration will stabilize at the parasitefree equilibrium. However, if there exists a P. falciparum infection, then an appropriate control in forms of effective antimalarial drugs would be necessary to clear the parasites from the human blood and restore the system to the stable parasitefree equilibrium state.
3.4. Coexistence of ParasitePersistent Equilibrium Point
The existence of a nontrivial equilibrium point represents a scenario in which the P. falciparum parasites are present within the host and the following conditions hold: . Upon equating the righthand side of system (7)–(14) to zero and solving for the state variables, we obtain the parasitepersistent equilibrium point , wherewhere .
Using Descartes’ “Rule of Signs” [66], it is evident that irrespective of the sign of in (53), in (56), in (59), in (60), in (63), and in (66), the state variables can only have one real positive solution. This shows that the model system (7)–(14) has a unique parasitepersistent equilibrium point .
3.5. Stability of the Coexistence of ParasitePersistent Equilibrium Point
Here, we shall prove that the coexistence of parasitepersistent equilibrium is locally asymptotically stable when . We shall follow the methodology by Esteva and Vargus presented in [67], which is based on the Krasnoselskii technique [68]. This methodology requires that we prove that the linearization of system (7)–(14) about the coexistence of parasitepersistent equilibrium does not have a solution of the formwhere , , and the real part of ξ is nonnegative (). Note that is a set of complex numbers.
Next, we substitute a solution of the form (69) into the linearized system (7)–(14) about the coexistence of parasitepersistent equilibrium. We obtain