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 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

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 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

Separating the negative terms, we obtain the following system:wherewith

Note that and all the elements in the square matrix H are nonnegative. The coordinates of are all positive, and the coordinate of the vector is described by the notation for . Additionally, the equilibrium satisfies . If we assume, for example, that system (73) has a solution of the form , then there exists a small positive real number ϵ, such that , where . Note also that is a norm in the field of complex numbers.

Next, we show that . To do so, we apply proof by contradiction. We let and . For the case when , the determinant () of (70) is given bywhere the positive terms are as defined in matrix (29).

It is clear that the above determinant is nonnegative (). Consequently, the system (70) can only have the trivial solution . On the contrary, for , we assume and define . This implies that and . The minimality of ϵ means that . While considering the nonnegativity property of H, if we assume the norms on the two sides of (73), we shall have

This implies that , which is a contradiction. Therefore, and is locally asymptotically stable when .

4. Numerical Simulations

4.1. Boundary Equilibrium Points

In this section, we show by means of numerical simulation the existence and stability of a positive parasite-persistent equilibrium point that involves only one of the parasite strains under study.

4.1.1. Drug-Sensitive-Only Persistent Equilibrium Point

This is an equilibrium point where only the drug-sensitive parasite strains are present in the infected human host. That is, the populations . This steady state is only feasible if no resistant parasites emerge from infected red blood cells and the use of antimalarial treatment does not lead to resistance development; that is, . The original model (7)–(14) is thus reduced to

Numerically, this equilibrium point is illustrated, as shown in Figure 2.

4.1.2. Drug-Resistant-Only Persistent Equilibrium Point

In this case, the population of the drug-sensitive parasite strains declines to zero as the density of the resistant strains grows and stabilizes at an optimal population size. This is also illustrated numerically, as shown in Figure 3.

4.2. Within-Host Competition between Parasite Strains

We investigate the competitive exclusion principle by simulating the model system (7)–(14) under different values of the threshold quantities and in (31). Model (7)–(14) is simulated so that and , and we achieve a convergence to the drug-sensitive-only endemic equilibrium point , as shown in Figure 4(a). Again, using the parameter values in Table 3 with and (, ), the solutions of and converge to the drug-resistant-only endemic equilibrium point (Figure 4(b)).

Provided that both and are greater than 1 (as shown in Figure 4(c)), the parasite-infected red blood cells remain persistent in the host. This implies that the merozoites (both drug-sensitive and drug-resistant) continue to multiply in the absence of antimalarial therapy, , or in the presence of ineffective antimalarial drugs. Similar results are observed in the dynamics of merozoites ( and ), as shown in Figure 5. It should be noted that the dominant merozoite strains are likely to drive the infection under these conditions. As the density of one strain increases, the population of the other strain is likely to decrease due to a phenomenon known as competitive exclusion principle. The most fit parasite strain survives as the weaker competitor dies out, as shown in Figure 5(a). Both drug-sensitive and drug-resistant merozoites would remain persistent if poor-quality antimalarial drugs are administered to P. falciparum malaria patients. Thus, in the absence of efficacious antimalarial drugs like ACTs with the potential to eradicate resistant merozoites, we are likely to experience an exponential growth in the density of drug-resistant merozoites, as displayed in Figure 5(b). This may lead to severe malaria and eventual death of the patient.

The bifurcation analysis of both scenarios is presented in Figure 6 (with and without competition between the parasite strains). When there is competition between the parasite strains, as shown in Figure 6(a), we observe that the strain with a higher threshold quantity would exclude the other strain. A decrease in the population of the drug-sensitive strain would pave way for a surge in the population of the drug-resistant strains, and vice versa. This is despite the fact that some drug-resistant strains emerge from the drug-sensitive strains as a result of mutation [77]. In Figure 6(b), we observe coexistence of the strains that do not compete with each other. Like the resistance strains, the sensitive strains are only present when their threshold quantity, , is greater than unity. Both strains are however present when and . Additionally, when and , we arrive at the parasite-free equilibrium (PFE) point, as shown in Figures 6(a) and 6(b).

4.3. Antimalarial Drug Effects and Parasite Clearance

The effects of antimalarial drug treatment are monitored by establishing first and foremost that

Thus, is a decreasing function of (the efficacy of the antimalarial drug used). Therefore, using a highly efficient antimalarial drug could lead to a scenario where and (disease-free state shown in Figure 7(c)). In Figure 7(a), model system (7)–(14) is simulated by varying the efficacy of the antimalarial drug and other model parameters chosen such that and . The higher the efficacy of the used antimalarial, the lower the density of infected erythrocytes. Thus, governments and ministry of health officers should only roll out or permit the administration of antimalarials or ACTs that can eradicate (totally) both the drug-resistant and the drug-sensitive strains of P. falciparum parasites.

The rate of development of resistance by the drug-sensitive merozoites, , is shown to have very minimal impact on the dynamics of infected red blood cells as long as and (Figure 7(b)). Nevertheless, analytical results indicate that the higher the rate of development of resistance, the lower the severity of future malaria infections. This is presented as

Other parameters that have direct negative impacts on the progression of malaria infection are the efficacy of the immune effectors, γ, and the rate of therapeutic elimination of drug-sensitive merozoites, ζ:

Further simulations based on contour plots (see [78] for theory on contour plots) are used to ascertain the relational effects of selected pairs of model parameters on the disease threshold quantities and . In Figure 8(a), both β and increase the reproduction number due to drug-sensitive P. falciparum parasite strains. A direct relationship exists between the two parameters: the higher the decay rate of the immune cells, the higher the rate of infection of healthy erythrocytes.

In Figure 8(b), we observe the least increase in with respect to an increase in relative to . Antimalarial therapy is shown to be very effective in reducing the severity of P. falciparum infection. Conversely, the number of merozoites produced per dying blood schizont, P, is shown in Figure 8(c) to have a very high positive impact on and hence on the severity of malaria infection due to drug-sensitive parasite strains. Clinical control should target and eradicate infected red blood cells to diminish the erythrocytic cycles of infections.

We observe in Figure 9(b) that the rate at which merozoites develop resistance due to treatment failure has no resultant effects on the rate of formation of gametocytes that undergo sexual reproduction within the mosquito vector. The higher the value of , the higher the cost of resistance, as shown in Figure 9(a). The higher the density of drug-resistant parasite strains, the higher the level of resistance and hence the cost of disease control. Unfortunately, highly effective antimalarial drugs (such as ACTs) that can eradicate both parasite strains are slightly expensive in several P. falciparum malaria-endemic regions [79]. Like the parameter P, the parasite infection rate β is shown to have a direct positive effect on the threshold quantity (Figure 9(c)) due to drug-resistant parasite strains. Effective antimalarials should hence target new cell infections and eliminate recrudescence (by killing already infected erythrocytes).

5. Effects of Multiple-Strain Infection and Fitness Cost on Parasite Clearance

Numerous studies [27, 80] have suggested the negative impacts of drug resistance on the fitness and ability of the parasite to dominate the P. falciparum infection. Resistance to antimalarial drugs imposes fitness cost on the drug-resistant parasite. The drug-resistant parasite strains are thought to experience impaired growth within the human host [29]. The cost of resistance is further exacerbated due to the competition between parasite strains within an infected human host. In Figure 10(a), the area under the curve for the drug-resistant strain or the number of infected erythrocytes is lower than that of the drug-sensitive strains. However, in a multiple-strain infection (Figure 10(b)), the area difference is much bigger. This implies that competition between the parasite strains within the human host could result in elimination of one of the parasite strains provided that both and are less than unity.

The presence of multiple strains of P. falciparum parasites is likely to complicate and worsen the severity of malaria disease infection in humans. Figures 11 and 12 show the simulated model (7)–(14) for single- and multiple-strain infections, in the absence of preexisting immunity and antimalarial drugs. The persistence of gametocytes in Figures 11(b) and 12(b) is consistent with the actual observations of human malaria infection in the absence of antimalarial therapy [81]. Acquired immunity is shown in Figure 11(c) to increase and eventually level-off at higher levels to contain future infections.

Although the aspect of timing is key in these multiple-strain infections, we assumed here that the two strains are introduced at the same time. In the long run, it is evident in Figures 10 and 12 that the sensitive strain overtakes the resistant strain. We argue that this could be as a result of strain-specific adaptive responses that symmetrically affect the sensitive parasites.

Unlike single-strain P. falciparum parasite infections, data on multiple-strain infections are not readily available. Nevertheless, a multiple-strain infection (drug-sensitive and drug-resistant) as presented in this paper is biologically reasonable and consistent with that of P. Chabaudi described in [82].

5.1. Sensitivity Analysis

In this paper, the primary model output of interest for the sensitivity analysis is the infected erythrocytes (). However, the effective reproduction number is a threshold quantity which represents on overage the number of secondary infected erythrocytes due to merozoite invasions. We can therefore measure the sensitivity indices of the effective reproduction number of model system (7)–(14) relative to model parameters. For example, the sensitivity of relative to the parameter is given by the following formulation:

Using the parameter values in Table 3, the expressions for sensitivity for all the parameters in are evaluated and presented in Table 4. The higher the numerical value of the sensitivity index (S.I), the greater the variational impact of the parameter on the disease progression. A parameter with a negative index decreases the model when they are increased. On the other hand, a parameter with a positive index would generate a proportional increase in when they are magnified. Results shown in Table 4 indicate that the rate of infection of healthy erythrocytes by the merozoites β, the density of merozoites generated from each of the bursting schizonts P, the efficacy of antimalarial drug used , and the rate at which drug-sensitive merozoites develop resistance are the four most influential parameters, in determining the disease dynamics as presented in model system (7)–(14).

Results from sensitivity analysis emphasize the use of highly efficacious antimalarial drugs such as ACTs in malaria-endemic regions. This would mitigate the many cases of malaria in the region and further help to reduce emerging cases of parasite resistance to existing therapies. Drugs with a higher parasite clearance rate would greatly reduce resistance, which is associated with longer parasite exposure to antimalarial drugs. It is imperative, therefore, that governments and ministry of health personnel in malaria-endemic countries enforce the use of efficient antimalarial drugs that not only cure infected malaria patients but also eliminate the chance of P. falciparum parasites to develop resistance to existing therapy.

6. Conclusion

In this paper, a deterministic model of multiple-strain P. falciparum malaria infection has been formulated and analysed. The parasite strains are categorized as either drug-sensitive or drug-resistant. The infected erythrocytes and the malaria gametocytes are similarly grouped according to the strain of the parasite responsible for their existence. The immune cells are incorporated to reduce the invasive characteristic of the malaria merozoites. Antimalarial therapy is applied to the model but only targets red blood cells infected with drug-sensitive merozoites. Based on the next-generation matrix method, we computed the effective reproduction number of the formulated model. Based on , it is evident that the success of P. falciparum infection in the presence of multiple-parasite strains is directly dependent on the ability of the individual parasite strains to drive the infection. The parasite strain with a higher threshold value, , is likely to dominate the infection. Prescribed antimalarial drugs should therefore be effective enough to eradicate both drug-sensitive and drug-resistant parasite strains in vivo. Linearization of the model at the parasite-free equilibrium reveals the local asymptotic stability of the trivial equilibrium point.

By rewriting the model in the pseudotriangular form, the parasite-free equilibrium is also shown to be globally asymptotically stable. Although the parasite-persistent equilibrium exists, its expression based on a single-model variable proved to be mathematically intractable. The use of antimalarial treatment may eradicate one parasite strain so that we arrive at either a drug-sensitive-only persistent equilibrium point or a drug-resistant-only persistent equilibrium point.

To assess the impacts of the different parasite strains to disease dynamics, the model is simulated for different values of the threshold quantities and . We observed that when and , then both parasite strains are persistent and the infection becomes severe. If and , then the drug-sensitive parasites would decline to zero as the drug-resistant strains continue to multiply and remain persistent, increasing the severity of infections. On the other hand, if and , then the drug-resistant parasite strains would be eradicated. Moreover, provided that the threshold quantities and are less than unity, the use of an efficacious antimalarial drug would help eradicate P. falciparum infection.

The efficacy of antimalarial drug is shown to have direct negative impact on the density of infected red blood cells. The higher the efficacy of administered antimalarial drug, the lower the population of infective merozoites and the smaller the density of infected erythrocytes. This ensures prompt recovery from malaria infections. This result is consistent with that in [72, 83]. The efficacy of antimalarial drug is however shown to have least effect on the population of drug-resistant infected erythrocytes. The rate of development of resistance by drug-sensitive parasites is also shown to drive the infection due to resistant parasite strains. Using contour plots and results from sensitivity analysis, we observe that the efficacy of antimalarial drug used , the density of blood floating merozoites produced per infected erythrocyte P, the rate of development of resistance , and the rate of infection by merozoites β are the most important parameters in the disease dynamics and control.

Finally, 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. To improve malaria therapy and reduce cases of parasite resistance to existing therapy, our results call for regular and strict surveillance on antimalarial drugs in clinics and hospitals in malaria-endemic countries.

Data Availability

All data used in this study are included in this published article.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this article.

Authors’ Contributions

All authors contributed to all sections of this manuscript.

Acknowledgments

The authors are thankful to the anonymous referees for their constructive comments. The authors would also like to thank the Strathmore Institute of Mathematical Sciences for its support in the production of this manuscript. The authors acknowledge with gratitude the financial support from the German Academic Exchange Service (DAAD) (ST32-PKZ:91711149) and the Kenya National Research Fund (NRF-Kenya) (NRF-Phd Grant Titus O.O), in the production of this manuscript.