Overestimating resistance in field testing of malaria parasites: simple methods for estimating high EC₅₀ values using a Bayesian approach

Abstract

Conventional methods of assessing in-vitro antimalarial drug-concentration effect relationships in field testing of fresh isolates assess each parasite isolate individually. This leads to systematic overestimation of EC₅₀ values for the most resistant isolates, and thus overestimation of the degree of resistance. In antimalarial drug-susceptibility studies conducted on the north-western border of Thailand the overestimation of EC₅₀ for the most resistant isolate ranged from 15% for artesunate to 43% for mefloquine. If isolates cannot be stored for re-testing, more accurate estimations of the degree of resistance can be obtained using a Bayesian approach to data analysis which is described here.

Background

The development of resistance to antimalarial drugs poses one of the greatest threats to malaria control and is the main cause of recent increases in malaria morbidity and mortality. The precise quantitation of resistance is therefore of prime importance. Initially the only way of assessing resistance to antimalarials was by inference from clinical treatment failures. Occasionally these were supported by measurement of antimalarial drug concentrations in the patients' blood. Once methods of culturing Plasmodium falciparum became established, in vitro methods for measuring the effect of the antimalarial drug directly on malaria parasite were developed [1]. This allowed resistance (i.e. reduced susceptibility) to be differentiated from poor adherence or unusual pharmacokinetics as the cause of treatment failure.

In in vitro susceptibility tests, blood samples from malaria patients are obtained and the infecting malaria parasites are cultured ex-vivo in the presence of stepwise increases in the concentrations of antimalarial drugs. Some methods call for adaptation of parasites to culture first, while others put blood directly from patients into the test system. Field testing, where blood is taken and malaria parasites are cultured directly in 96-well plastic plates pre-dosed with antimalarial at different concentrations, is now widely used. These freshly obtained parasites are usually not cryo-preserved and so there is only the one opportunity to assess the parasite drug susceptibility.

Problems with estimating EC₅₀ for resistant isolates

The three most commonly used methods of in vitro susceptibility testing of malaria parasites are (i) the micro test which assesses inhibition of parasite growth to the schizont stage microscopically, (ii) the radioisotope test which measures uptake of ³H-hypoxanthine, and (iii) ELISA based methods which measure the production of lactate dehydrogenase or histidine rich protein by the parasite [2, 3]. Results of in vitro tests are expressed as the percentage parasite growth/viability plotted against the antimalarial drug concentration (or log concentration) to give a dose-response (concentration-effect) curve. The range of drug concentrations used are chosen because they are thought a priori to reflect the likely range of parasite susceptibilities and they are usually prepared as serial doubling dilutions The resulting dose-response curve is usually sigmoid (Figure 1) and can be described by a number of parameters: minimum growth/uptake/production, maximum growth/uptake/production, the slope, and importantly the midway point (EC₅₀) describing the concentration of drug which is required for 50% of the maximum inhibitory effect in the test system. These variables can be read directly from the plot, assessed by probit analysis, or they require fitting of a statistical model to data and then derivation from the equation of the model. A commonly used model is the sigmoid Emax model:

Figure 1

Standard dose-response curve in antimalarial drug susceptibility testing (arbitrary units). In this case the response or effect is inhibition of growth, uptake, or synthesis.

E(C) = E _max - (E _max - E _min )C^γ/(C^γ+EC ₅₀ ^γ)

where E_max is the maximum effect (e.g. minimum growth), E_min is the minimum effect (maximum growth), C is the concentration of drug, γ is the slope of the linear part of the curve, and EC₅₀ is defined as before. Parameter E_min sometimes is set to 0 to avoid overparameterization.

The most resistant isolate(s), by definition, represents the extreme value in one tail of a distribution of values. These values lie at the lower end of the dilution range (Figure 2). As concentrations are taken usually as serial dilutions, there are less measurements of drug effect at this end of range (because most of the concentrations in the test system have no effect). The relatively large intervals between the highest drug concentrations in the test system mean that often the middle range of inhibitory effect is not observed at all; for one dilution there may be no (or very little) drug effect and for the next dilution the maximum (or close to maximum) effect is measured (Figures 1 to 3). If the distribution of all EC₅₀ values in the parasite population is unimodal (and this is a critical point) then, given that the EC₅₀ lies at an unknown point between two values (two dilutions), the prior probability suggests that the true value lies closer to the median value of the population than to the more extreme value.

Figure 2

Illustrative example of a distribution of antimalarial susceptibility concentration-growth inhibition curves fitted using standard models. The most resistant isolate's curve fit is poor as inhibition was obtained only with the highest concentration tested (64 ng/mL), and so the EC₅₀ derived using standard curve fitting lies at the mid-point (48 ng/mL) between the two highest concentrations tested (32 and 64 ng/mL). The Bayesian approach (dotted line) provides a lower estimate closer to the rest of the population.

Figure 3

Dose-Response relationship; real data from ³hypoxanthine uptake inhibition in-vitro susceptibility tests for the most resistant isolates in the series (6) for (A) Artesunate, (B) Halofantrine, (C) Mefloquine, (D) Mefloquine, (E) Chloroquine, (F) Chloroquine.

But the fit, and thus estimated EC₅₀ value, is based only on the observations for that isolate. It is independent of the distribution of values from other isolates in the series and, in case of the sigmoid model, a symmetrical curve is fitted. The fit is poor but the midpoint value will lie close to half way between the two concentrations tested (Figure 2).

Standard fitting procedures, therefore, systematically overestimate EC₅₀ values for the most resistant parasites. Moreover, when curve fitting is attempted for resistant isolates and only responses close to 100% or 0% are observed, often either computational problems occur and the curve cannot be fitted at all, or the curve is fitted but the standard errors of the estimates are large. In this case, a priori information about the population could be very useful in planning the future experiments and deciding where the EC₅₀ is expected to be and where the concentrations measurements should be taken in order to pinpoint the 50% response. In this report a Bayesian approach to the estimation of EC₅₀ in the outlying most resistant isolates together with worked examples from field data is presented.

Methods

These suggestions apply to field methods of assessing in vitro susceptibility to malaria parasites where the inhibitory values (EC₅₀, EC₉₀, etc.) are extrapolated from observations over a preset concentration range. This range is set based on the earlier experience with field testing. The criteria for "goodness of fit" have not been standardized and, therefore, vary between investigators. The variation will affect the results, but this important issue is not discussed here.

Estimating EC₅₀ based on the population estimates

If, from the in vitro drug susceptibility experiment, it is known that the EC₅₀ lies in the interval (a to b), where b = 2a if doubling solutions were used, the Bayesian estimate of EC₅₀ in this interval can be constructed in the following way:

• 1. estimate distribution of EC₅₀ values for this batch of isolates (using 'good' estimates coming from 'good' fits). This involves estimating EC₅₀ for each isolate.

• 2. using distribution of EC₅₀, calculate median value of EC₅₀ in the interval (a,b)

• 3. this summary measure value is the Bayesian estimate of EC₅₀ in the interval (a,b).

Characterization of the distribution of EC₅₀ values is required, and this is not always possible, but if the EC₅₀ values can be shown to have a normal or lognormal distribution or can be transformed to normal using Box-Cox transformation [4], then the estimation of the corrected EC₅₀ value can be done using a pocket calculator, as described below. The corrections suggested apply to continuous distributions and are not appropriate for clearly discontinuous distributions (for example, those seen where one or more parasite isolates contain the cytochrome b mutations which confer high level atovaquone resistance in P. falciparum).

Normal distribution

Probability density f(x) of a normally distributed random variable x is given by the expression

f(x) = 1/(σ$\sqrt{2}$π)·exp(-(x-μ)²/2σ²)

where exp(z) = e^z, μ is the expectation or mean value of x and σ is the standard deviation of x.

The median value of x on the interval (a, b) is calculated (see Appendix II) as:

Median (x) = Φ^-1 ((Φ((a-μ)/σ) + Φ((b-μ)/σ))/2)·σ + μ     (1)

where Φ(t) = ${\displaystyle {\int }_{-\infty }^{\text{t}}\text{f}(\text{z})}$ dz and z is a standard normal distribution N(1,0). Values of Φ(z) are tabulated and can be found in any book with statistical tables.

Lognormal distribution

A random variable y is log-normally distributed if x = log(y) is normally distributed with log denoting the natural logarithm. The general formula for the probability density function of the lognormal distribution is

f(y) = 1/((y-θ)σ$\sqrt{2}$π)·exp(-(log(y-θ)/m)²/2σ²)

where σ is the shape parameter, θ is the location parameter and m is the scale parameter (equal to median). The calculation of the median value in the interval (a,b) may be performed using the log-transformed data. Firstly, the median of transformed x needs to be found in interval (log(a), log(b)). As transformed x is normally distributed, methods described in the paragraph above apply and equation (1) can be used. Then the median needs to be back-transformed (using the exponent function) to the original scale. Because of monotonic transformation, the median of transformed x in interval (log(a), log(b)) corresponds to the median of x in the interval (a,b)

Median (x) = exp (Φ^-1 ((Φ((log(a)-μ)/σ) + Φ((log(b)-μ)/σ))/2)·σ + μ)     (2)

Box -Cox family of transformations

The Box-Cox transformation is defined as:

T(y) = (y^λ-1)/λ

where y is the random variable and λ is the transformation parameter. For λ = 0, the natural log of the data is taken instead of using the above formula. Similarly as before, if T(y) has an median value m in interval (T(a), T(b)) then median of y in interval (a,b) is m^1/λ.

Simple method using computer simulation

An easier way to obtain the Bayesian estimate of the EC₅₀, for data where the effect goes from 100% to 0% over a single dilution, which does not require any equations, is to use computer simulations. They can be done in Excel^®, or with any statistical software. Once the distribution of EC₅₀ is established, a large number of points (for example 10⁶) from this distribution can be readily generated and then the estimate of the EC₅₀ for a resistant isolate can simply be found as a median value for all data points which are in the interval (a,b) between the two dilutions.

Data

The data used to illustrate these issues come from the Shoklo Malaria Research Unit (SMRU) in Thailand. They are in-vitro susceptibility data of malaria parasites isolated from patients recruited in two camps for displaced persons of the Karen ethnic minority situated in an area of forested hills on the north-western border of Thailand. Antimalarial drug susceptibility data were obtained using the hypoxanthine uptake inhibition assay and have been described and analysed elsewhere [5]. In total, 268 fresh isolates of P. falciparum from primary infections were assayed for in vitro drug susceptibilities to a wide range of antimalarials including chloroquine diphosphate, quinine citrate, mefloquine hydrochloride, halofantrine hydrochloride, artesunate, dihydroartemisinin, artemether, lumefantrine, and atovaquone. To illustrate this particular issue we selected data for four drugs; mefloquine, chloroquine, halofantrine and artesunate.

Statistical Analysis

For each of the examined drugs, the distribution of EC₅₀ values was estimated in the following way:

• 1. alll isolates which had at least one inhibitory response value between 30 and 70% of the maximum response were identified and selected as this was considered a priori to be a minimum requirement for the data points to give stable estimates

• 2. the E_max model using Stata^® software (StataCorp. 2005, ver 9) was fitted to each of the selected isolates

• 3. isolates for which the model did not reach convergence were excluded

• 4. isolates which had negative estimates of EC₅₀ were excluded

• 5. the distribution of EC₅₀ values in the remaining isolates (N) was examined

• 6. If possible, EC₅₀ values were transformed (logarithmic or Box-Cox transformation) to normality. Normality of the distributions was tested using the Kolmogorov-Smirnov test.

For each of the isolates which were most resistant to mefloquine, halofantrine, chloroquine and artesunate the EC₅₀ values were estimated using three methods: (a) by fitting a standard 3-parameter sigmoid curve using WinNonlin^® (Pharsight, ver 4.1); (b) by fitting the 3-parameter sigmoid curve using Gibbs sampling using WinBUGS^®; (c) by the "Bayesian" method based on the population distribution of EC₅₀ values proposed in this paper (equation 1).

The standard sigmoid model (a) was fitted using Nelder-Mead minimisation, with no specified boundaries for parameters and WinNonlin^® generated starting values. In the Gibbs sampling estimation (b) the distribution of EC₅₀ values was estimated as described above, while non-informative uniform priors with appropriate bounds were used for other parameters:

Emax ~ uniform(0.5,1.5)

γ ~ uniform(1,100)

The effect of the bounds of parameter γ on estimates of EC₅₀ were also examined by restricting the upper bound to 30, and then to 50. Parameter estimates were summarized as median and 95% range of the posterior estimates over 90,000 iterations.

Results

The measured responses of isolates most resistant to mefloquine (two series), halofantrine, chloroquine (two series) and artesunate are presented in Figure 3.

Table 1 describes estimated distributions of the EC₅₀ values for each drug. All of them could be assumed to be log-normally distributed. Overall, after the selection described in points 1, 3 and 4, about half of the isolates contributed to the estimation.

Table 1 Estimated distribution of EC₅₀ values

Table 2 lists parameter estimates obtained from the WinNonlin analysis (a) for the resistant isolates. Only for artesunate and chloroquine do the estimates have reasonable precision, for the other drugs the standard errors (CV%) are very large. This means that the model does not estimate the parameters accurately. In fact, when refitting the model with different starting values, different estimates of EC₅₀ (within ± 10% of the quoted value) were obtained each time (data not shown). In Table 3, corresponding parameter estimates obtained from fitting the three-parameter sigmoid curve using the Gibbs sampling method (b) are listed. Obviously the slope of the concentration-effect relationship is a critical determinant of EC₅₀ but the EC₅₀ estimates were relatively insensitive to these prior slope estimates (γ); change of the upper bounds for parameter γ did not affect substantially the estimated values of parameters whereas slight increases in the EC₅₀ estimates (<10%) were observed with lower bound changes (Table 4).

Table 2 Parameter estimates obtained from fitting the 3-parameter sigmoid curve to the most resistant isolates using standard individual data analysis (fitted using WinNonlin^®).

Table 3 Parameter estimates obtained from fitting the 3-parameter sigmoid curve to the data presented in the example using a Bayesian approach (Gibbs sampling) (fitted using WinBUGS^®)

Table 4 Estimates of EC₅₀ obtained from fitting the 3-parameter sigmoid curve to the data presented in the Example using WinBUGS, with upper boundary of γ set to 30 or 50.

Calculation of the Bayesian EC₅₀

Halofantrine

In total 197 isolates were tested, of which 85 satisfied the criteria (see above, points 1–5). EC₅₀ values estimated in these samples were found to have a lognormal distribution (p = 0.100, Kolmogorov-Smirnov test) with mean of 1.84 and standard deviation of 0.979; i.e. log(EC₅₀) ~ N(1.84, 0.979) (Table 1).

From the equation (2) the EC₅₀ was calculated as 41.2 ng/mL in the highest concentration interval (32.6 – 65.32 ng/mL).

Artesunate

Of 191 isolates, 89 satisfied the criteria and for them log(EC₅₀) ~ N(0.75, 0.888), p = 0.074 (Table 1) In the highest concentration interval (16.73–33.46 ng/mL) an EC₅₀ of 20.3 ng/mL is obtained.

Mefloquine

Out of 216 isolates, 124 satisfied the prespecified criteria and for these the EC₅₀ values have a log-normal distribution (p = 0.366), log(EC₅₀) ~ N(3.68, 0.633) (Table 1). In the highest concentration interval (115–230 ng/mL) an EC₅₀ of 137.8 ng/mL was calculated.

Chloroquine

Out of 183 samples, 117 satisfied the criteria and for them log(EC₅₀) ~ N(4.9, 0.680), p = 0.076 (Table 1). In the highest concentration interval (341.5–683 ng/mL) EC₅₀ of 420.4 ng/mL was obtained.

In Table 5 Bayesian EC₅₀ estimates (c) are compared with the standard estimates obtained using WinNonlin^® (b). The standard method oversestimated the EC₅₀ of the most resistant isolate by between 15% (artesunate) and 43% (mefloquine).

Table 5 Methodological differences in the estimation of antimalarial resistance; the highest EC₅₀ (ng/mL) (Brockman et al, 2000)

Discussion

Antimalarial drug resistance is widely monitored using in-vitro susceptibility testing. There are sentinel sites throughout the malaria affected world monitoring for drug resistance. A variety of methods have been developed and the results have provided valuable information in the assessment and mapping of antimalarial drug resistance [7]. Evaluation of stored isolates in reference centres allows proper standardisation of methodologies and repeated tests on single isolates. But most of this testing in the field is a "one-off" microtest on freshly obtained blood samples. When antimalarial drug susceptibility tests are reported the highest observed values observed are naturally of greatest interest as they may represent emerging drug resistance. Drug regimens should aim to cure all infections, and thus provide concentrations exceeding the inhibitory concentrations for the most resistant prevalent parasites. If only a single concentration range is evaluated in an in-vitro susceptibility assay using serial dilutions then, by definition, the true EC₅₀ of the most resistant isolates must lie above or between the largest concentration differences tested. Thus, unless the parasites are retested with a higher concentration range (which they usually cannot be), the precision of the estimated EC50 or EC90 value of the most resistant isolate will usually be the poorest of all the isolates assayed. Furthermore as curve fitting or probit analysis takes no account of other isolates in the series tested, then if there are two adjacent points with extremely different values (often zero and 100% inhibition), a curve will be fitted as lying symmetrically between the two. The EC₅₀ will be assessed as lying close to the mid-point between the two concentrations (Figure 2). But if the parasites can be shown to derive from a single distribution of susceptibilities, and this proviso is critical, the prior probability is that the true EC₅₀ value lies closer to the population mean value. Thus resistance is systematically overestimated. A better estimate is provided by the simple Bayesian analysis described above. Such a continuum conforming to a single distribution of susceptibilities is observed commonly for artesunate, dihydroartemisinin, artemether, artemisinin, chloroquine, desethylamodiaquine, quinine, quinidine, lumefantrine, piperaquine, pyronaridine and mefloquine. For example, the artesunate-mefloquine combination has been systematically deployed for over twelve years on the north-western border of Thailand. It has been reported in studies of antimalarial susceptibility that the most resistant isolate EC₅₀ values were 23.4 ng/mL for artesunate and 197 ng/mL for mefloquine. Reanalysis of the data using the Bayesian approach reduces this to 20.3 (15% less) ng/mL for artesunate and 138 ng/mL (43% less) for mefloquine. When antimalarial drug resistance is reported precise details of the concentration range tested and the analytical procedure used should always be provided.

For those resistance mechanisms in which single mutations confer large reductions in susceptibility, such as the Pfdhfr 164 mutation for pyrimethamine resistance or cyt b 268 mutations for atovaquone resistance there will clearly be discontinuous distributions of susceptibility and this method will not be appropriate. This emphasizes the importance of distributional assessments before using this Bayesian approach to analysis. Of course if parasites can be cryopreserved then the resistant parasites dose-response measurements should be repeated with selection of higher concentrations covering the likely EC₅₀ region. For freshly assayed isolates this will not be possible as the level of susceptibility is not known before the test. In this case a Bayesian adjustment should be made for extreme values if justified by the distributional assessment.

Conclusion

Conventional analytical methods for characterizing in-vitro antimalarial drug susceptibility assess each isolate independently and, consequently, overestimate the EC₅₀ and EC₉₀ for the most resistant isolates. Bayesian methods based on the distribution of EC₅₀ values in the whole series offer considerable improvement of the estimate. The method which is proposed here, does not require sophisticated software, nor does it make any assumptions about the dose-response relationship, and it provides more realistic estimates of the most resistant isolates' EC₅₀ values

Appendix I

Useful facts :

• 1. If f(x) is a density function for random variable x then

F(z) = P(x<z) = ${\displaystyle {\int }_{-\infty }^{\text{z}}\text{f}(\text{x})\text{ dx}}$; ${\displaystyle {\int }_{-\infty }^{\infty }\text{f}(\text{x})\text{ dx}}$;

P(a<x<b) = F(b)-F(a)

where P(x<z) denotes probability that random variable x is less then z.

• 2. f(x) = 1/(σ$\sqrt{2}$π)·exp(-(x-μ)²/2σ²) density function for N(μ, σ)

• 3. f(x) = 1/($\sqrt{2}$π)·exp(-x²/2) density function for N(0,1)

• 4. Φ(x) = 1/($\sqrt{2}$π) ${\displaystyle {\int }_{-\infty }^{\text{z}}{\text{exp(-x}}^{\text{2}}\text{/2) dx}}$

• 5. If x is random variable with distribution N(μ, σ) then

F(x) = Φ((x-μ)/σ))

Appendix II

• 1. Let assume x is N(μ, σ) with density function f(x) (as in point 2, Appendix I)

• 2. We want to find median value of x in the interval (a,b)

• 3. We need to define density function g(x) such that:

$\text{g}(\text{x})=\{\begin{array}{cc}0& \text{if x}<=\text{a}\\ \text{Af}(\text{x})& \text{if x in interval}(\text{a},\text{b})& \text{and}& {\displaystyle {\int }_{-\infty }^{\infty }\text{g}(\text{x})\text{ dx}=1.}\\ 0& \text{if x}>=\text{b}\end{array}$

$${\displaystyle {\int }_{-\infty }^{\infty }\text{g}(\text{x})}$$

dx = A (Φ((b-μ)/σ)) - Φ((a-μ)/σ))) = 1

Therefore A = 1/(Φ((b-μ)/σ)) - Φ((a-μ)/σ))).

• 4. Median value of x in interval (a,b) is such a z that

G(z) = P(x<z) = ${\displaystyle {\int }_{-\infty }^{\text{z}}\text{g}(\text{x})}$ dx = 0.5

$${\displaystyle {\int }_{-\infty }^{\text{z}}\text{g}(\text{x})}$$

dx = A ${\displaystyle {\int }_{\text{a}}^{\text{z}}\text{f}(\text{x})}$ dx = A (Φ((z-μ)/σ)) - Φ((a-μ)/σ))) =

(Φ((z-μ)/σ)) - Φ((a-μ)/σ)))/(Φ((b-μ)/σ)) - Φ((a-μ)/σ))) = 0.5

After transformation we get:

Φ((z-μ)/σ) = (Φ((a-μ)/σ) + Φ((a-μ)/σ))/2

and

z = Φ^-1 ((Φ((a-μ)/σ) + Φ((a-μ)/σ))/2)·σ + μ