The emergence of parasite resistance to antimalarial drugs has contributed significantly to global human mortality and morbidity due to malaria infection. The impacts of multiple-strain 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 multiple-strain P. falciparum infection and control model. Although the drug-resistant 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 malaria-endemic regions. Our results emphasize the call for regular and strict surveillance on the use and distribution of antimalarial drugs in health facilities in malaria-endemic countries.

1. Introduction

The emergence of parasite resistance [14] to antimalarial drugs has contributed significantly to human mortality and morbidity due to malaria infection, worldwide [57]. 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 drug-treated 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 sulfadoxine-pyrimethamine (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 artemisinin-based combination therapies (ACTs) as a first-line 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 first-line 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 first-line antimalarial drug by several countries until 1957, when the first focus of P. falciparum resistance was detected along the Thai-Cambodia 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 sulfadoxine-pyrimethamine (SP) as a first-line 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. KwaZulu-Natal 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 half-life 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 artemisinin-based combination therapies (ACTs) as the first- and second-line 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 half-life, eliminates the remaining parasites [26]. The WHO currently recommends five different ACTs: (1) artesunate-amodiaquine (AS + AQ), (2) artesunate-mefloquine (AS + MQ), (3) artesunate + sulfadoxine-pyrimethamine (AS + SP), (4) artemether-lumefantrine (AM-LM), and (5) dihydroartemisinin-piperaquine (DHA + PPQ). Additionally, artesunate-pyronaridine 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 cross-sectional 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, drug-sensitive parasites are shown to strongly suppress the growth and transmission of drug-resistant P. falciparum parasites. Although high-transmission settings such as sub-Saharan Africa account for about 90% of all global malaria deaths, resistance to antimalarial drugs has been shown to emerge from low-transmission settings, such as Southeast Asia and South America [29]. Causes of parasite resistance to ACTs are diverse. Historical studies [30, 31] indicate that antimalarial-resistant 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 high-transmission regions. Moreover, owing to repeated exposure for many years, individuals in high-transmission settings are likely to develop clinical immunity to malaria, leading to stronger selection for resistance [34]. Studies in [29] also support the hypothesis that in-host competition between drug-sensitive and drug-resistant parasites could inhibit the spread of resistance in high-transmission 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 in-host 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 within-host parasite competition are shown to inhibit the spread of resistance [44, 45]. On the contrary, some models [39, 46] have suggested that within-host competition is likely to speed up the spread of resistance in high-transmission settings due to a phenomenon called “competitive release.” In this paper, we provide theoretical insights using mathematical modelling of the impacts of multiple-strain 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 within-human malaria model that has both the drug-sensitive and drug-resistant P. falciparum parasite strains subject to antimalarial therapy. In Section 3, we analyze the model based on epidemiological theorems. Within-host competition between parasite strains and the effects of antimalarial drug efficacy on parasite clearance are discussed in Section 4. Sensitivity analysis and multiple-strain 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 multiple-strain infection due to P. falciparum.

2. Model Formulation

We present in this paper a deterministic model that describes the within-human-host competition and transmission dynamics of two strains of P. falciparum parasites during malaria infection. The compartmental model considers the coinfection and competition between the drug-sensitive (dss) and the drug-resistant (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 density-dependent 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 drug-resistant and drug-sensitive 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 drug-sensitive and drug-resistant 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 drug-sensitive merozoites, whereas refers to the RBCs infected with drug-resistant merozoites. Similarly, the variables and represent drug-sensitive and drug-resistant 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 [5053] to model the reductive effects of the immune cells on the parasite and infected-cell 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 drug-sensitive and drug-resistant 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 drug-sensitive 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 drug-sensitive IRBCs burst open to generate more drug-sensitive merozoites or drug-sensitive gametocytes at the rate . Similar dynamics are observed with the drug-resistant IRBCs, where the drug-resistant 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 drug-infested erythrocytes are hence likely to die faster. This is represented by the term , where represents the antimalarial-specific treatment efficacy. In this paper and for purposes of illustration and simulations, corresponds to the efficacy of artemether-lumefantrine (AL), which is the recommended first-line 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 drug-resistant merozoites and the drug-resistant gametocytes die naturally at the rates and , respectively. It is further assumed that drug-sensitive merozoites and gametocytes may develop into drug-resistant 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 drug-resistant 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 drug-sensitive 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 right-hand 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 Parasite-Free Equilibrium Point (PFE)

The in-host 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 parasite-free equilibrium point given by

Using the next-generation 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 parasite-free equilibrium is the spectral radius of the next-generation matrix , where

It follows thatwhere

From equation (31), it is evident that, in a multiple-strain P. falciparum malaria infection, the progression of the disease depends on the reproduction number of different parasite strains. If the threshold quantity , the drug-sensitive parasite strains will dominate the drug-resistant 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 drug-sensitive parasites. Conversely, if , the infection is mainly driven by the drug-resistant parasite strains. In this scenario, the used antimalarial drugs should be highly efficacious and effective enough to kill both the drug-resistant and drug-sensitive 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 parasite-free equilibrium point is locally asymptotically stable if and unstable otherwise.

The Jacobian matrix associated with the in-host 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


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 parasite-free equilibrium point. This is presented in the following section.

3.3. Global Asymptotic Stability Analysis of the Parasite-Free 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 parasite-free 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 parasite-free 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 off-diagonal entries of are nonnegative (equal to or greater than zero), showing that is a Metzler matrix. To show the global stability of the parasite-free 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 parasite-free 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 parasite-free 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 parasite-free equilibrium state.

3.4. Coexistence of Parasite-Persistent 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 right-hand side of system (7)–(14) to zero and solving for the state variables, we obtain the parasite-persistent 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 parasite-persistent equilibrium point .

3.5. Stability of the Coexistence of Parasite-Persistent Equilibrium Point

Here, we shall prove that the coexistence of parasite-persistent 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 parasite-persistent 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 parasite-persistent equilibrium. We obtainwhere , , and .

Upon simplifying the equations in (70), we obtainwhere