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

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 EC50 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 EC50 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 predosed 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 50 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 3 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 50 ) 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: 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 50 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 50 values in the parasite population is unimodal (and this is a critical point) then, given that the EC 50 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.
But the fit, and thus estimated EC 50 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 50 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 50 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 50 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 50 , EC 90 , 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.
Standard dose-response curve in antimalarial drug suscepti-bility testing (arbitrary units) 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. Drug concentration Growth/uptake/production (%)

Estimating EC 50 based on the population estimates
If, from the in vitro drug susceptibility experiment, it is known that the EC 50 lies in the interval (a to b), where b = 2a if doubling solutions were used, the Bayesian estimate of EC 50 in this interval can be constructed in the following way: 1. estimate distribution of EC 50 values for this batch of isolates (using 'good' estimates coming from 'good' fits). This involves estimating EC 50 for each isolate.
2. using distribution of EC 50 , calculate median value of EC 50 in the interval (a,b) 3. this summary measure value is the Bayesian estimate of EC 50 in the interval (a,b).
Characterization of the distribution of EC 50 values is required, and this is not always possible, but if the EC 50 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 50 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 cyto-Illustrative example of a distribution of antimalarial susceptibility concentration-growth inhibition curves fitted using standard models 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 50 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.

Normal distribution
Probability density f(x) of a normally distributed random variable x is given by the expression 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: where Φ(t) = 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 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

Response (%)
Concentration (ng/mL) (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)

Box -Cox family of transformations
The Box-Cox transformation is defined as: 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 50 , 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 50 is established, a large number of points (for example 10 6 ) from this distribution can be readily generated and then the estimate of the EC 50 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 50 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 For each of the isolates which were most resistant to mefloquine, halofantrine, chloroquine and artesunate the EC 50 values were estimated using three methods: (a) by fitting a standard 3-parameter sigmoid curve using Win-Nonlin ® (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 50 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 50 values was estimated as described above, while noninformative 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 50 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 50 values for each drug. All of them could be assumed to be lognormally distributed. Overall, after the selection described in points 1, 3 and 4, about half of the isolates contributed to the estimation.  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 50 but the EC 50 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 50 estimates (<10%) were observed with lower bound changes ( Table 4).

Calculation of the Bayesian EC50 Halofantrine
In total 197 isolates were tested, of which 85 satisfied the criteria (see above, points 1-5
In Table 5 Bayesian EC 50 estimates (c) are compared with the standard estimates obtained using WinNonlin ® (b). The standard method oversestimated the EC 50 of the most resistant isolate by between 15% (artesunate) and 43% (mefloquine).

Discussion
Antimalarial drug resistance is widely monitored using invitro 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 50 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 50 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 50 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 50 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 doseresponse measurements should be repeated with selection of higher concentrations covering the likely EC 50 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.