A Mathematical Analysis of Competitive Dynamics and Aggressive Treatment in the Evolution of Drug Resistance in Malaria Parasites

: Experimental evidence supports the counterintuitive notion that rapid eradication of pathogens within a host, infected with both drug-sensitive and -resistant malaria parasites, can actually accelerate the evolution of drug-resistant pathogens. This study aims to analyze the competitive dynamics between these two strains through a mathematical model and evaluate the impact of aggressive treatment on the spread of drug resistance. We conducted equilibrium, uncertainty, and sensitivity analyses to assess the model, identifying and measuring the influence of key factors on the outcome variable (the population of drug-resistant parasites). Both equilibrium and local sensitivity analyses concurred that the density of drug-resistant parasites is notably affected by genetic instability, the production rate of red blood cells, the number of merozoites, and competition factors. Conversely, there is a negative relationship between genetic instability and one of the competition coefficients. Global sensitivity analysis offers a comprehensive examination of the impact of each input parameter on the temporal propagation of drug resistance, effectively accounting for the interplay among parameters. Both local and global sensitivity analyses underscore the continuous impact of drug treatment on the progression of drug resistance over time. This paper anticipates exploring the underlying mechanisms of drug resistance and providing theoretical support for developing more effective drug treatment strategies.


Introduction
Drug resistance in malaria refers to the ability of strains of Plasmodium parasites to persist and reproduce even when exposed to recommended or higher doses of antimalarial drugs [1].This resistance arises due to the widespread use of common antimalarial medications, leading to resistance in severe forms of the disease.Despite intensified efforts for control, progress towards malaria elimination has stalled in recent years [2,3].The World Health Organization has recommended several artemisinin-based combination therapies (ACTs) for the treatment of uncomplicated falciparum malaria [4].While ACTs have significantly reduced malaria morbidity and mortality rates globally, the challenge of drug resistance persists, particularly in regions such as Southeast Asia [5], South America [6], Papua New Guinea [7], and Eastern Africa [8].
In theory, drug treatment is a crucial intervention for patients infected with malaria, as it can alleviate symptoms and reduce mortality rates.However, during drug treatment, sensitive parasites are typically eliminated.If the concentration achieved is insufficient to kill all drug-resistant mutants, populations of surviving resistant mutants will rapidly expand, enabling their survival and reproduction [9].This scenario increases the likelihood of drug resistance spreading to the next generation, potentially leading to its widespread prevalence.Importantly, experimental evidence suggests that, when an aggressive treatment regimen is employed, sensitive parasites are swiftly eradicated, potentially giving a competitive advantage to drug-resistant parasites [10][11][12].Therefore, the competition between sensitive and drug-resistant parasites is crucial in this context [13][14][15].In contrast, a more moderate treatment approach allows some sensitive parasites to persist and compete with drug-resistant strains, thereby suppressing the latter and delaying the onset of resistance [16].Based on this analysis, it is evident that drug treatment plays a pivotal role in influencing the development of drug resistance.Consequently, there is an urgent need to explore the relationship between drug-resistant parasites and drug concentration, which is a topic of paramount importance.
Mathematical modeling is considered a crucial tool for understanding the evolution of drug resistance.This tool not only elucidates interactions between different parasite strains, but also provides a theoretical foundation for managing resistance.In recent years, various epidemiological models have emerged for studying drug-resistant parasites.One study indicated that high transmission intensity may hinder the evolution of resistance by affecting host population immunity [17].Similarly, another study found that early transmission of drug-sensitive parasites and asymptomatic infections can delay the onset of drug resistance [18].Additionally, models have focused on the relative size of the host population infected with resistant parasites and its impact on the spread of drug resistance [19].Mathematical models for malaria have integrated treatment and varying levels of resistance to analyze the role of treatment in the development of drug resistance [20].Concurrently, research has investigated the potential contribution of co-infection competitive release to the emergence of resistance [21].While some studies aim to devise optimal strategies for controlling resistance spread through model analysis, research simultaneously addressing drug treatment and competition is relatively scarce [22,23].Therefore, it is necessary to construct a comprehensive model that integrates these key elements and draws upon experimental findings from multiple sources [15,16,24].
One of the primary objectives of this paper is to identify influential factors and explore the relationship between the density of drug-resistant parasites and input parameters.Typically, this goal can be achieved by directly computing the mathematical expression representing the state variable of drug-resistant parasites in our model.However, due to the complexity of our model, deriving a specific formulation for the state variable of drug-resistant parasites proves to be a challenging task.Consequently, we are compelled to seek alternative approaches.Uncertainty and sensitivity analyses serve as important tools for gaining insights into the mechanisms governing system behaviors.In particular, uncertainty analysis helps quantitatively assess biases in the results of epidemiological studies, enhancing our understanding of the inherent uncertainties in the data and methodologies used [25].By conducting sensitivity analysis, we can not only identify the parameters significantly influencing the outcome variables of interest, but also quantify how input uncertainties impact the overall results of the model [26].
This paper focuses on a thorough investigation of the model proposed in [27], which considers the dynamics of drug treatment and within-host competition between drugsensitive and drug-resistant parasites.Built upon empirical observations [15,16,24], these observations highlight a counterintuitive phenomenon: rapid pathogen clearance in a host co-infected with both drug-sensitive and drug-resistant malaria parasites may paradoxically accelerate the emergence of drug resistance.This paradox arises from the competitive interactions between the two parasite strains.To validate and refine the model, we aligned it with experimental data to facilitate parameter estimation for subsequent statistical analysis.Subsequently, we conducted a series of analytical procedures, including uncertainty analysis, local sensitivity analysis, and global sensitivity analysis, to investigate the behavior of the model under various conditions.We particularly focused on the state variable representing the density of drug-resistant parasites, depicting it as a function of various input parameters through uncertainty analysis.Sensitivity analysis helps clarify how changes in the input parameters influence the evolution of drug resistance within the parameter space defined by the parameter ranges.In conclusion, the outcomes of this study hold the potential to provide insights for formulating optimal policies to manage drug resistance and control disease transmission.

Model
In this paper, we considered the model proposed in [27]. ( The model considers the dynamic changes in three distinct cell populations: uninfected red blood cells denoted as S(t), drug-sensitive malaria parasites represented by I s (t), and drug-resistant malaria parasites indicated as I r (t).Uninfected red blood cells are continually supplied at a constant net rate Λ from the bone marrow and have a natural life expectancy of 1/d 1 days.Malaria parasites infect red blood cells at a rate proportional to the contact rate associated with their population size.The transmission rates for the wild-type and drug-resistant parasites are denoted as β 1 and β 2 , respectively.Following a series of multiplication cycles, each infected red blood cell produces α infective merozoites (in the case of P. falciparum, each infected cell yields a substantial number of merozoites, typically ranging from 8 to 32 in humans).The parameter p represents the proportion of drug-sensitive parasites released from a red blood cell infected by a drug-resistant strain.The competition between the two parasite strains is characterized by the terms γ 1 I s (t)I s (t) and γ 1 I s (t)I s (t), which align with the formulation used in numerous intra-host competition models.Both parasite strains experience mortality at per capita rates denoted as d 2 and d 3 .Additionally, the clearance rate of the wild-type strain is represented by µ.

Parameter Estimation
Parameter estimation was carried out within the dynamic framework of malaria infection, encompassing both competition and drug treatment.The model parameter values were estimated using the least squares method to fit a dataset on malaria infection from [24].The least squares problem involves minimizing the discrepancy between the sum of squared differences and the observed data over a parameter vector space constrained by a predefined feasible region.It is crucial to recognize that the standard least squares formulation comprises two fundamental components: mathematical and statistical [28].The mathematical model utilized in this study is the malaria model (1) we developed.The statistical model is built on the assumption that the output of the model and associated random deviations (measurement errors) are governed by the random variables: where θ = (α, β 1 , β 2 , γ 1 , γ 2 , d 2 , d 3 , µ, p) and I r (t j ; θ) denotes the output of the mathematical model (1).The random variables ϵ j model the random deviations away from I r (t, θ) and are assumed to satisfy that all ϵ j are independent and identically distributed random variables, E(ϵ j ) = 0 for every j and var(ϵ j ) = σ 2 < ∞ for every j.The objective function is given by where θ OLS is the estimator, the set of parameter vectors θ constrained by a prespecified feasible region Θ.Θ is primarily determined based on the biological significance of the parameter.j represents the time point at which the time series data are observed, and n is the number of data points available for inference.Hence, the model solution I r (t; θ) yields the best fit to the time series data Y j .The lsqcurve f it function in Matlab 2015b (Mathworks, Inc., Natick, MA, USA) is used to find the least squares best fit to the data.Moreover, in Matlab, two numerical optimization methods are available to solve the nonlinear least squares problem: The trust-region-reflective algorithm and the Levenberg-Marquardt algorithm.

Equilibrium Analysis
In order to explore the relationship between the drug-resistant parasites and each parameter, computing the interior equilibrium is of primary importance.To determine the interior equilibrium of system (1), we denote P * = (S * , I * s , I * r ) and let the right-hand side of (1) be zero.It then follows from the last two equations of (1) that and Substituting I * r and I * s into the first equation of (1), we obtain the S * that satisfies where It is noteworthy that not all I * r values can keep positive due to the complex existence condition of I * r .According to Theorem 3 of [27], we can obtain the existence condition of the interior equilibrium P * .Model (1) also exhibits a fixed-point bifurcation.For details, please refer to Reference [27].
Moreover, the interior equilibrium exists by further assuming that S * (one of S * ± ) satisfies (H1)

Uncertainty Analysis
Uncertainty analysis was formally used to quantify the variability in the outcome variable, which was due to the uncertainty in estimating the values of the input parameters [29,30].This type of analysis provides insight into the range of potential results that the model can generate.To initiate the uncertainty analysis, it was essential to define the prior distribution of the input parameters under consideration.Subsequently, a Latin Hypercube Sampling (LHS) method was applied.The outcomes were visualized through a histogram, shedding light on how the distribution changes when different assumptions and prior distributions for the bias parameters are altered [25].To gain a comprehensive understanding of the outcome variable, various statistical measures such as the mean, variance, median, quantiles, and others were employed as valuable approaches.

Local Sensitivity Analysis (LSA)
Local sensitivity refers to the sensitivity of parameters with respect to a given parameter set, where only one parameter varies in a certain range and the other parameters are fixed.If the output variable y can be defined as follows: where x = (x 1 , x 2 , . . ., x n ) represent the input parameters, the usual method to compute the sensitivity of the i-th input parameter on the outcome variable is by the partial derivative at the nominal value of parameter x i .Derivative-based sensitivity analysis finds its rationale in the Taylor series expansion [31].Note that it was hard to quantify the sensitivity if x i and y varied in different units of measurement; therefore, normalizing the sensitivity indices was necessary, which took the form: where σ x i and σ y represent the standard deviation of input parameter x i and outcome variable y in Monte Carlo sampling, respectively.The normalized sensitivities of the model output are then given as L 1 and L 2 measure the impact of the input parameter on the model outcome variable to give a preference for tuning the parameter with a high sensitivity.Derivative-based sensitivity methods are widely performed due to the efficient calculations, but this also has the drawback that they provide information about local sensitivity only.

Global Sensitivity Analysis (GSA)
Compared to the local sensitivity analysis, GSA is performed over the entire parameter space, where all parameter variations are simultaneously considered.In addition, the interaction between the input parameters we are concerned with would be calculated.In this paper, we mainly apply the variance-based GSA method, i.e., the Sobol sensitivity method.The essence of Sobol sensitivity analysis is to measure sensitivity through variance decomposition [32,33].According to Equation ( 6), x = (x 1 , x 2 , . . ., x n ) represent the input parameters, which are defined on an n-dimensional unit cube H n = {x|0 ≤ x i ≤ 1, i = 1, 2, . . ., n}.Based on the Sobol method, y = f (x) can be decomposed into single model parameters and the subitem function of the parameter interaction: The number of all subitems is 2 n , and the subitem function is obtained by calculating the following multiple integrals: Similarly, other subitem functions of high order can be achieved.f 0 is a constant, and the integral of every summand over any of its own variables is zero: The subitem functions in Equation (10) satisfy Equation (11), and it can be inferred that every subitem function in Equation ( 9) is orthogonal, which is: ..,j l (x j 1 , . . ., x j l )dx = 0, for (x i 1 , . . ., x i s ) ̸ = (x j 1 , . . ., x j l ).(12) Based on the above properties, the variance of output variable V(y) can be decomposed as follows: where V i 1 ,...,i s = 1 0 f 2 i 1 ,...,i s (x i 1 , . . ., x i s )dx i 1 , . . ., x i s . is the partial variance corresponding to the subitem function of Equation ( 9).The Sobol global sensitivity indices are defined by We mainly considered first-order Sobol indices and total Sobol indices, which are and Equation ( 15) is the first-order Sobol index.Equation ( 16) is the total-order sensitivity index, as an extension of the Sobol sensitivity indices, which is defined as the ratio of the sum of the related sensitivity indices.Equations ( 15) and ( 16) explain the reason why the higher value of S i , the more importance of parameter x i on the outcome variable, but this is not applicable in reverse.In other words, a low value of the first-order index cannot necessarily stand for having no influence on the outcome variable due to the parameters' interactions, which lead to a larger total-order index.Hence, total-order indices have a great significance.Parameter x i has no impact on the outcome variable in the case ST i = 0 and vice versa.In summary, first-order Sobol indices refer to the influence of each parameter on the outcome variable alone, which measure the importance of each parameter itself, whereas the total Sobol indices not only characterize the contribution of the concerned parameters, but also their interactions.The calculation of Sobol sensitivity indices involves the computation of multiple integrals, which is very complicated and difficult, especially for complex nonlinear models.Therefore, the Monte Carlo method is employed to approximate the multiple integral solution.

Parameter Estimates
Optimal fits were obtained for the dataset in Figure 1.The model has a goodness-of-fit ratio R 2 = 0.9436.In this paper, we fixed the parameter value of the production rate of RBCs (Λ) and the decay rate of RBCs (d 1 ) based on the research of Anderson [34].The estimated parameter values are shown in Table 1, where d 3 < d 2 accorded with the fact that the resistant mutants have a higher death rate than wild-type strains, reflecting the biological cost of resistance [19].In addition, the solution curve was obtained on the basis of the parameter set we estimated (Figure 2), and the three cell populations involved in this model finally reached a steady state.

Equilibrium Analysis
Equilibrium analysis is employed to observe the relationship between the parameters and the drug-resistant parasites in this paper.Calculating the interior equilibrium of the dynamical system is a necessary step, then we primarily focus on the part of I r in the equilibrium.Equilibrium analysis results are obtained in Figure 3.The results show a strong positive linearity between α and I * r as α varies in the interval [8,32] (Figure 3a), indicating the more merozoites were produced, the higher the value of I * r was.As the parameters β 1 and β 2 increase, I * r exhibits the tendency of decreasing and increasing, respectively (Figure 3b,c), but in general, their impact on I * r is not obvious, such that their effective intervals are much smaller.The trends of γ 1 (rising) and γ 2 (falling) on I * r are opposite, which exactly shows the importance of competition, that drug-resistant parasites could be controlled by competition to some extent, and the resistant parasites could not be suppressed once they obtain a strong competitive ability.Moreover, the results of the upward tendency of µ and d 2 just coincide with the research of Hansen [16] that aggressive treatment (a high level drug concentration) can fail and threaten patient health and lifespan.I * r reduces with the increasing d 1 and d 3 ; the reason is obvious: that d 3 , as the decay rate of drug-resistant parasites, will inevitably lead to the reduction of I * r with its increasing.The rising death rate of uninfected RBCs indicates the reduction of the resources of drugresistant parasites, thereby leading to their decreasing.Therefore, on the contrary, Λ, as the production rate of uninfected RBCs, has a positive effect on I * r .In addition, the negative linear relationship between p and I * r demonstrated that the gene stability can affect I * r to a great extent.

Uncertainty Analysis Results
Uncertainty analysis is performed in order to describe the range of possible outcomes given a set of inputs (where each input parameter has some uncertainty).Therefore, in this paper, uncertainty analysis would be performed to describe the range in the density of drug-resistant parasites in steady states, given that a range of input parameters is simulated.There were 50,000 input parameters samples extracted to obtain 11,888 effective samples that yielded a positive equilibrium.Descriptive statistics (mean, var, etc.) were calculated to describe the distribution of the outcome variables I * r that we were concerned with.The LHS uncertainty technique is a formal approach to quantify the effect the uncertainty of the input variables on the prediction precision of the output variables for which we aimed.First, we sampled all parameters by using the LHS method (Figure 4).The descriptive statistics of I * r is shown in Table 2.The acquired frequency distribution in Figure 4

Sensitivity Analysis Results
Sensitivity analysis is performed in order to describe how sensitive the outcome variables are to the variation of individual input parameters.Since there may be multiple input parameters, sensitivity analysis could determine which ones drive the majority of the variation in the outcome.

Local Sensitivity Analysis
Local sensitivity refers to the sensitivity of parameters with respect to a given parameter set.This type of analysis is usually carried out by partial derivatives of the output functions at a fixed and predefined point of the parameter space.We are mainly concerned with the local sensitivity of the parameters we estimated in Table 1.The LSA results offered a sensitivity order, which is p, α, β 2 , Λ, d 3 , γ 1 , γ 2 , β 1 , d 1 , d 2 , µ from large to small by the measurement of L 1 and L 2 .The mean of the sensitivity index has the same direction as the equilibrium analysis results of each parameter.Table 3 and Figure 5 illustrate that d 2 and µ have a little effect on I r in the local parameter space, and p affect I r most near the estimated parameter, indicating that the proportion of drug-sensitive parasites released from an infected RBC by drug-resistance parasites is fairly important in a local parameter space.On the contrary, the direct effect of the decay rate of drug-sensitive parasites and the drug concentration was not significant may due to ignoring the interaction between these parameters.
i )/n are the L 1 -and L 2 -norm, respectively.b The Mean column is the mean of the sensitivity functions.Local sensitivity of state variable I r (drug-resistant parasites) on all parameters.The light grey region represents from the minimum value to the maximum value; the dark grey region is from the mean value minus the standard deviation to the mean value plus the standard deviation.The light blue region ranges from the 5th percentile to the 95th percentile; the dark blue region ranges from the 25th percentile to the 75th percentile.

Global Sensitivity Analysis
When the sensitivity index for an output variable is examined over the time course of the simulation, the relative importance of the parameters is shown to vary widely depending on the time point of inspection.The global sensitivity analysis results were obtained by employing the Sobol method, which provided a factor-based decomposition of the output variance.The first-order indices only account for the direct contribution of each parameter.The total Sobol indices quantified the fraction of the variance to explain the variability of each parameter full consideration of their interaction effects.Figure 6 depicts the first-order Sobol indices and the total Sobol indices; both of them range from 0 to 1 for all parameters.The first-order Sobol index of γ 2 is high in the first few days, but then decreases sharply, and other parameters remain zero in the concerned time course.The total Sobol index of γ 2 is high during the time interval, and p, Λ, β 2 , α, µ, and d 2 increase to a high value with time varying, indicating their relative importance to the drugresistant strains.The other remaining parameters including β 1 , d 1 , and γ 1 are relatively less sensitive to resistant strains.Combining the two figures, we could deduce that, between those parameters, an interaction may exist.

Discussion and Conclusions
We have studied model (1) with the aim of exploring the impact of within-host competition and drug treatment on drug-resistant parasites.The results of our analysis have not only provided predictions, but have also painted a clearer picture of the influential parameters, laying the groundwork for the formulation of effective strategies in resistance management.
One of the key parameters, denoted as α, exhibits a significant correlation with drugresistant parasites (I * r ), underscoring the importance of the infected red blood cell's capacity to produce merozoites.The infection process is characterized by repetitive cycles of invasions and merozoite production, as detailed by Anderson [34].Consequently, inhibiting merozoite invasion of erythrocytes emerges as a viable strategy to reduce the number of drug-resistant parasites.The rate of malaria parasite infection plays a pivotal role in determining the invasive and infectious abilities of parasites [35].As a result, variations in β 1 and β 2 can directly impact the population of drug-resistant parasites, as illustrated in Figure 3b,c.It is worth noting a phenomenon wherein increasing β 1 and β 2 leads to diminishing changes in I * r , suggesting that, within a certain range, adjusting the infection rate can effectively manage the evolution of drug resistance.
There is a growing interest in the potential of drug-based strategies beyond the initial treatment of malaria cases to further reduce morbidity and mortality and advance towards malaria elimination.Nevertheless, it is crucial to acknowledge the harsh reality that the emergence of resistance to anti-malarial drugs or combination therapies exposes the vulnerability of current strategies.The rapid elimination of drug-sensitive parasites can create a favorable environment for the survival of resistant parasites [24,36], a phenomenon well explained by the trends in I * r seen in Figure 3g,i.Furthermore, aggressive elimination of sensitive strains can lead to competitive release of drug-resistant parasites, given the absence of competitors (drug-sensitive parasites) to suppress them, highlighting the significance of competition between the two parasite types [16].This relationship is depicted in Figure 3d,e, illustrating how the competitive ability of resistant strains directly affects their population.Another essential factor in the control of drug resistance is the genetic stability of resistant parasites, as evidenced by the decline in their numbers in Figure 3j.Hence, efforts to enhance the genetic stability of resistant parasites hold promise in resistance management.Finally, the tendencies of I * r in response to varying parameters such as Λ, d 1 , and d 3 are both reasonable and evident.
To further elucidate the sources of uncertainty, sensitivity analysis emerges as a fundamental technique, encompassing both local sensitivity analysis (LSA) and global sensitivity analysis (GSA).Typically, conducting LSA involves calculating the derivative of the outcome variable concerning each input parameter, making slight perturbations near the best point estimates.Numerous studies have delved into LSA, often by analyzing the basic reproduction number (R 0 ) to examine disease equilibria in various epidemiological models [37,38].In some cases, additional conditions are necessary because R 0 alone may not determine the existence of interior equilibria.Consequently, we have opted to directly focus on the interior equilibrium for our sensitivity analysis, utilizing numerical methods to gain deeper insights.Our LSA findings emphasize the significance of genetic instability and the number of merozoites in influencing the outcomes of the model.Intriguingly, drug concentration appears to be the least sensitive parameter to resistant strains.This observation may be attributed to the inherent limitations of LSA, which does not account for interactions between inputs, potentially making it challenging to establish a definitive ranking of the key parameters.
Global sensitivity analysis, on the other hand, excels in addressing this limitation by considering interactions among input parameters and varying all of them simultaneously within a certain range for epidemiological models.In contrast to LSA, GSA highlights the importance of µ, possibly due to its strong interactions with other parameters.Additionally, γ 2 maintains the highest sensitivity throughout the analysis, suggesting that drug therapy indirectly induces changes in the population of resistant parasites, potentially through competition dynamics.The implementation of aggressive treatment could trigger competitive release of drug-resistant parasites, as sensitive strains diminish.Hence, the competitive aspect should not be underestimated, and enhancing competition between the two strains may offer a viable approach to controlling drug-resistant parasites.The temporal variation in the sensitivity of state variables to parameters is a valuable aspect of infectious disease modeling studies [39].Notably, the sensitivity of µ declines around the 40th day, possibly resulting from the elimination of drug-sensitive parasites by treatment, leaving fewer sensitive parasites to make a diminishing contribution.Parameters p, α, Λ, and β 2 continue to play a similarly crucial role compared to the LSA and maintain a high level of sensitivity throughout the simulation period.
Generally, Sobol global sensitivity analysis allows for accurate assessment of the parameters' main effects and interactions and is applicable to various parameter space dimensions and model types, albeit with high computational costs.In contrast, MeFAST [40] incorporates the use of multiple statistical tests for the significance of the parameters; however, its performance may be suboptimal in highly nonlinear systems, and its evaluation of high-order interactions among parameters is limited.The Morris method [41] incurs lower computational costs, making it suitable for high-dimensional parameter spaces and complex models; nevertheless, its evaluation of interactions among parameters is limited.Derivative-based global sensitivity [40] assessment methods can provide rapid and accurate results, but are only applicable to differentiable functions and may have limited sensitivity evaluation capabilities for highly nonlinear systems.Therefore, when selecting methods, it is necessary to consider a combination of factors such as computational costs, applicability, and evaluation accuracy.
The model effectively aligns with the available data, yielding a set of estimated parameter values that form the foundation for subsequent analysis.In our investigation, we conducted equilibrium analysis, uncertainty analysis, and sensitivity analysis to explore the parameters that significantly influence and are chiefly responsible for the development of drug resistance.It is important to note that our study does not aim to directly propose specific drug resistance management strategies, but rather, intends to offer recommendations derived from mathematical findings and logical reasoning for more effective drug resistance control.The scope of this research solely considers the impact of drug treatment on sensitive strains.In future work, we plan to explore an extended dynamical system that incorporates the immune response into the model [42][43][44].The immune response can play a pivotal role in preventing the re-invasion of merozoites or eliminating infected red blood cells [45].

Figure 1 .
Figure 1.Drug-resistant parasites' data (red dot) with model best fit (blue line).Each data point is assumed to have the same measurement error, estimated from competitive release of drug resistance data in the research ofRead [24].The y-axis is plotted using a log scale, reflecting the drug-resistant parasite density.

Figure 2 .
Figure 2. The solution curve of the model with the estimated parameters according to the data.The solid line, the dashed line, and the dotted line represent red blood cells, drug-sensitive parasites, and drug-resistant parasites, respectively.

Figure 3 .
Figure 3. Relationship between eleven sampled input parameters and equilibrium of I r .These results were obtained using a sample size of 1000.(a) I r vs. α.(b) I r vs. β 1 .(c) I r vs. β 2 .(d) I r vs. γ 1 .(e) I r vs.γ 2 .(f) I r vs. d 1 .(g) I r vs. d 2 .(h) I r vs. d 3 .(i) I r vs µ.
illustrates a wide range of I * r for a great uncertainty in estimating the values of eleven input parameters, with a range from 4.61 × 10 −5 to 200.3142.Accordingly, 5% to 95% of I * r fell in the interval [1.5219, 51.7296] due to the left-skewed character of the frequency distribution.The variance was fairly large, which indicated a high volatility of I r , due to the uncertainty of all parameters.

FrequencyFigure 4 .
Figure 4. Uncertainty analysis of equilibrium of I r on all parameters.

Figure 6 .
Figure 6.Sobol sensitivity indices for state variable I r (drug-resistant parasites).The left graph depicts the tendency of the first-order Sobol index, and the right graph describes the total sobol index of all model parameters.

Table 1 .
Parameters in the model.

Table 2 .
Descriptive statistics from the uncertain analysis of Ir in Figure4.

Table 3 .
Value of local sensitivity function.
a Here, L 1