Analysis of ex vivo drug response data of Plasmodium clinical isolates: the pros and cons of different computer programs and online platforms

Background In vitro drug susceptibility testing of malaria parasites remains an important component of surveillance for anti-malarial drug resistance. The half-maximal inhibition of growth (IC50) is the most commonly reported parameter expressing drug susceptibility, derived by a variety of statistical approaches, each with its own advantages and disadvantages. Methods In this study, licensed computer programs WinNonlin and GraphPad Prism 6.0, and the open access programs HN-NonLin, Antimalarial ICEstimator (ICE), and In Vitro Analysis and Reporting Tool (IVART) were tested for their ease of use and ability to estimate reliable IC50 values from raw drug response data from 31 Plasmodium falciparum and 29 P. vivax clinical isolates tested with five anti-malarial agents: chloroquine, amodiaquine, piperaquine, mefloquine, and artesunate. Results The IC50 and slope estimates were similar across all statistical packages for all drugs tested in both species. There was good correlation of results derived from alternative statistical programs and non-linear mixed-effects modelling (NONMEM) which models all isolate data simultaneously. The user-friendliness varied between packages. While HN-NonLin and IVART allow users to enter the data in 96-well format, IVART and GraphPad Prism 6.0 are capable to analyse multiple isolates and drugs in parallel. WinNonlin, GraphPad Prism 6.0, IVART, and ICE provide alerts for non-fitting data and incorrect data entry, facilitating data interpretation. Data analysis using WinNonlin or ICE took the longest computationally, whilst the offline ability of GraphPad Prism 6.0 to analyse multiple isolates and drugs simultaneously made it the fastest among the programs tested. Conclusion IC50 estimates obtained from the programs tested were comparable. In view of processing time and ease of analysis, GraphPad Prism 6.0 or IVART are best suited for routine and large-scale drug susceptibility testing. Electronic supplementary material The online version of this article (doi:10.1186/s12936-016-1173-1) contains supplementary material, which is available to authorized users.

beneficial, it is important to consider that resistance to one component of the therapy can be masked by a partner drug which retains high anti-malarial efficacy. In vitro assays also provide an opportunity to assess drug susceptibility of parasites to individual drugs, thereby allowing preventive measures to be taken before clinical treatment failure occurs. In addition, in vitro assays enable the measurement of drug sensitivity without the confounding effects of clinical efficacy such as host immunity and the pharmacokinetics of the drug [2][3][4].
A variety of assays are available to measure drug susceptibility in Plasmodium falciparum. In most, the parasites' drug susceptibility is defined by measuring growth (i.e., schizont maturation) or replication (i.e., re-invasion assays) in the presence of varying concentrations of antimalarial compounds. Although these assays are highly informative, the comparison of data between laboratories or field sites is often problematic as numerous variations exist between protocols, such as the initial parasitaemia, incubation time, culture haematocrit, and the use of alternative media and supplements [5,6].
Standard analysis of in vitro (culture-adapted Plasmodium strains) and ex vivo (fresh Plasmodium clinical isolates) assay data is commonly conducted by using nonlinear regression to fit a sigmoid E max model to each sample's concentration-effect data. The sigmoid E max model is comprised of four parameters: minimum and maximum effects, steepness (=slope) of the dose-response curve, and concentration of the drug required to inhibit growth to 50 % of that observed in the absence of the drug (IC 50 ). In this approach, the definition of IC 50 is based on the assumptions that there is: (a) a monotonic relationship between the dose of the compound tested and the response in the assay; and, (b) a consistent 50 % response can be defined clearly [7]. IC 50 values can also be estimated using alternative statistical models, such as polynomial regression [8] and sigmoid inhibition models based on non-linear regression [9,10]. Determination of reproducible IC 50 values remains a challenge for investigators as it is time-consuming and often subjective as the process itself involves visual inspection of individual dose-response curves [10]. Moreover, laboratories have different criteria for accepting or rejecting assay data, leading to a variety of selection biases, such as excluding assays from highly resistant parasites for which IC 50 estimates tend to be less precise [11].
Different computer programs and online platforms are available, each with its own algorithms and features, to facilitate the processing of raw drug assay data. In this study, five of the most commonly used analytical platforms were assessed by applying them to the same dataset of ex vivo drug response data collected from patients with clinical malaria and comparing the derived output parameters, their utility, robustness, and consistency of results.

Clinical isolates
From April 2012 to January 2013, Plasmodium species isolates were collected from patients with malaria attending an outpatient clinic in Timika, Papua Province, Indonesia. In this region, clinical trials have confirmed high levels of multidrug-resistant P. vivax and P. falciparum [12][13][14]. Patients with symptomatic malaria were recruited into the study if they had a single species infection with P. falciparum or P. vivax, with a parasitaemia of between 2000 and 80,000 μL −1 as determined by microscopy, and with synchronous infection with at least 70 % parasites at ring stage. Patients were excluded from the study if they had taken any anti-malarials in the preceding month. After written informed consent was obtained, 5 mL venous blood was collected by venipuncture. Host white blood cells (WBC) were removed using cellulose column filtration and packed infected red blood cells (IRBC) used for the ex vivo drug susceptibility assay.

Ex vivo drug susceptibility assay
Plasmodium drug susceptibility was measured using a protocol modified from the World Health Organization (WHO) microtest as described previously [12,15]. Twohundred microlitres of a 2 % haematocrit blood medium mixture (BMM), consisting of RPMI 1640 medium plus 10 % AB+ human serum (P. falciparum) or McCoy's 5A medium plus 20 % AB+ human serum (P. vivax) were added to each well of pre-dosed drug plates containing 11 serial concentrations (twofold dilutions) of the antimalarials (maximum concentration shown in parentheses) chloroquine (2992 nM), piperaquine (769 nM), mefloquine (338 nM), artesunate (49 nM), and amodiaquine (158 nM). Pre-dosed drug plates were made by diluting the compounds in 50 % methanol, followed by lyophilization, and stored at 4 °C until use. All antimalarials were obtained from the World Wide Antimalarial Resistance Network (WWARN) Malaria Drug Reference Material Scheme [16].

Dose-response analysis
The raw dose-response data for each drug were visually checked to ensure that they were amenable to regression modelling. Five programs and/or online tools were evaluated.

WinNonlin 4.1 (Pharsight Corporation)
IC 50 values were estimated using the program WinNonlin 4.1 (Pharsight Corp) by fitting a sigmoid E max model (Eq. 1) to each individual isolate's concentration-effect data by non-linear regression E represents the percentage of schizonts observed after normalization to the control wells using the following equation: E raw is the percentage of schizont observed in the drug well and divided by the percentage of schizont observed in the control well (E control ). Maximum inhibition occurs when there is no observable schizont (i.e., inhibition of growth), that is, E equals 0. In Eq. 1, E 0 represents minimum per cent growth, E max maximum per cent growth, IC 50 the concentration of the drug required to inhibit 50 % of the control parasites' schizont growth, C the drug concentration, and γ the steepness of the curve (slope).
Initial raw data were transposed into an Excel converter worksheet to fulfil the program's template requirements. The converted data were then transposed into the program's data table and analysed using the pharmacodynamic inhibitory sigmoid E max model presented above (see Eq. 1). Key outputs and additional parameters are presented in Table 3. Warning messages were shown when potential inaccuracies were identified, such as steep response without observed values on the slope, and dose-response curves were provided for visual inspection.
Isolates with estimated E max and E 0 values above 15 % or below 100 and 0, respectively, were re-analysed. The majority of these results were caused by either resistant isolates in which the response was deemed inadequate at the highest drug concentration tested, or lower schizont counts at the lowest concentration tested than those of untreated controls. In such cases, additional extreme data points were added assuming E = 1 with no drug and 0 when very high drug concentrations were present.

GraphPad Prism (version 6.01; GraphPad Software, Inc.)
Dose-response analysis of the log-transformed drug concentrations (represented by X in Eq. 2) and response (represented by Y in Eq. 2) was performed using the 'log (inhibitor) versus response-variable slope' equation (Hill equation): Bottom represents the value of Y for the minimal curve asymptote (response in the absence of drug), Top is the 1 + 10 LogIC 50 10 x γ value of Y for the maximal curve asymptote (response produced by an infinitely high concentration of drug). LogIC 50 is the logarithm of the drug concentration required to inhibit 50 % of growth, and γ represent the slope of the curve. The upper and lower limits of γ were constrained to 10 and 0, respectively. Raw data were first transferred into an Excel converter worksheet to fulfil the program's template requirements. The converted raw data were then transposed into the program's Y-axis data table and the drug concentrations entered into the X-axis data table; the data were logtransformed prior to analysis. In GraphPad Prism 6.0, the results were available immediately after data entry. Key and additional outputs are outlined in Table 3. Doseresponse curves were also produced for visual inspection. The program showed 'ambiguous' text warnings when it was unable to find the combination of parameters that best fit the data. A warning text was shown for the following reasons: (i) data were not collected over a wide enough range of drug concentrations; (ii) the model had redundant parameters (i.e., the model could be fitted with multiple sets of parameters); or, (iii) the model resulted in a biphasic curve, leading to different IC 50 values [17]. In the case of 'ambiguous' results due to reason (i), additional extreme data points were added, assuming Y equals to 0 when very high drug concentrations were present, and data were re-analysed.

HN-NonLin
HN-NonLin is a freeware program [18] that uses logtransformed data to fit a polynomial regression model with the following equation (Eq. 3): The goal of polynomial regression is to find the values of the parameters A, B, C, etc., that are going to generate the best-fit curve for the observed data points. When using a polynomial regression model, investigators need to specify the order of the polynomial (i.e., the numbers of parameters to be fitted into the equation) [19]. In HN-NonLin, users can choose the order of the polynomial (i.e., 1-7) that will determine the shape of the doseresponse curve: Polynomial '1' will give a straight line, while polynomial '7' will produce a perfect interpolation (i.e., the regression line will pass through all data points). In this study, polynomial '3' was chosen since polynomials at higher order will produce perfect correlation, but undesirable oscillations caused by outliers [20].
To conduct data modelling using HN-NonLin, raw data were transferred into an Excel converter worksheet in order to fit the 96-well plate format of the program. After the polynomial was set, data modelling was performed and the key parameters, as well as additional outputs generated as presented in Table 3. As with the other programs, additional extreme data points were added in cases where drug-response data were deemed inadequate, assuming that no schizont was observed at very high drug concentrations.

ICEstimator 1.2
The ICEstimator is a freeware online program that carries out a non-linear regression on the relative concentration-effect points by using an inhibitory sigmoid E max model [9].
Raw data were transferred into an Excel converter worksheet to comply with the program's requirements. Converted raw data were then entered via the ICEstimator 1.2 website [21]. For the analysis of each isolate, drug concentration (x-axis) and response data (y-axis) needed to be entered individually. In the current study, the unit for concentration was set to 'nM' and for response 'other' was chosen. In addition, 'schizont counting' was chosen as the method for analysis. Data modelling was performed online generating the outputs outlined in Table 3. Dose-response curves were available for individual inspection. Warning signs were displayed if: (i) the ratio of the upper and lower limits of the 95 % CI of the IC 50 was greater than two; (ii) convergence was unsuccessful, even when the slope was set at 10; (iii) the parasite's growth ratio (mean effect at drug-free control/ mean effect at maximum drug concentration) was less than two; (iv) the lower limit of the IC 50 95 % CI was less than 0; and/or, (v) IC 99 > maximum concentration [9]. The results, as well as dose-response curves were displayed directly on the website and were available for download in portable data file (pdf ) and comma-separated value (csv) format.

In Vitro Analysis and Reporting Tool (IVART)
IVART is a freeware online program that was developed by WWARN. IVART is a program that is based on the code of the ICEstimator program and described in detail on the WWARN website [10].
Raw data were transferred into an Excel converter worksheet and converted into a 96-well format to comply with the program's requirements. Converted raw data were then entered into the WWARN pre-build Excel worksheet as per the developer's instructions [22]. This worksheet was then uploaded through the WWARN website [23] and data analysis performed online. The results took approximately 1 min to be generated following data upload and could be downloaded as pdf and/ or csv files. Key and additional outputs are outlined in Table 3. In the pdf output, warning messages were shown when (a) drug concentration ranges were incorrect; or, (b) core criteria, as described elsewhere [10], were not met. If warning signs were shown, additional drug concentrations and data points were added as described above and the data re-analysed.

Non-linear mixed-effects modelling using the package NONMEM
The estimates obtained by all the programs tested were compared with the estimates produced by NONMEM. NONMEM is a commercial software package that analyses data from all isolates simultaneously using nonlinear mixed-effects models to take into account the two sources of variability, within-and between-isolate variability. NONMEM allows the specification of any non-linear regression equation to describe the effectconcentration curves and can incorporate other confounding factors of the drug assay, such as parasite stage composition at the start of the assay, assay duration, and the use of different drug plate batches. For this analysis, the sigmoid E max model (i.e., Eq. 1) was used and no additional confounding factors were included. The derivation of IC 50 values using NONMEM is described in detail elsewhere [24].

Data analysis
Results from each program were pooled and statistical analysis conducted using GraphPad Prism 6.0. Geometric mean and 95 % reference range of IC 50 and slope parameters were calculated for each drug using all available programs. NONMEM was used as the reference method and parameter estimates compared with each of the other methods. Bland-Altman plots using the derived IC 50 estimates were created to compare agreement between programs.
In addition, a qualitative analysis was conducted to assess the ease of use and application of each methodology. This included how the raw data were processed, data conversion and/or transformation, data upload, data analysis, processing time, capability to perform batch analysis, capability of offline analysis, and costs.

Ethical approval
Ethical approval for this project was obtained from the Human Research Ethics Committee of the Northern Territory Department of Health and Families and Menzies School of Health Research (HREC 2010-1396), Darwin, Australia, and the Eijkman Institute Research Ethics Commission (EIREC-47), Jakarta, Indonesia.

Overview and comparison of outputs
Ex vivo drug susceptibility assays to chloroquine, amodiaquine, piperaquine, mefloquine, and artesunate were assessed in field isolates from 60 patients presenting with single-species infections of P. falciparum (n = 31) and P. vivax (n = 29). Drug susceptibility data were successfully derived from all samples. Baseline characteristics of the isolates processed are presented in Table 1 and the derived drug susceptibility for both species depicted in Fig. 1. The population geometric mean IC 50 s were similar for the five different methodologies with no significant differences when compared to the NONMEM reference estimates (see Tables 2 and 4). With the exception of HN-NonLin, all programs generated estimates of the slope (γ) and corresponding 95 % confidence intervals (Table 3).
For resistant isolates, both the IC 50 and the shape of the dose response curve shift, and hence both parameters, are required to calculate the inhibitory responses at other concentrations, such as the IC 90 or IC 99 . All programs produced dose-response curves for visual inspection, but only WinNonlin, GraphPad Prism 6.0, and HN-NonLin generated statistics for the goodness of fit of the final model. Although this statistic is not shown in ICE and IVART, warning messages are displayed when the regression fails to produce an acceptable fit.
Visual inspection of both raw data and dose-response curves is inevitable for final interpretation of IC 50 estimates. When data cannot be modelled for reasons such as similar schizont growth at each drug concentration compared to the drug-free control, or non-sigmoidal curves, none of the software packages generated IC 50 values or drug response curves. In these cases, ICE, IVART and WinNonlin also displayed error messages. For assays with a large discrepancy of response between duplicate concentrations on the slope, GraphPad Prism 6.0 and WinNonlin tended to produce higher IC 50 estimates and lower slope values compared to the other packages. This shift in IC 50 s was most prominent when  one of duplicates was zero on the slope (i.e., GraphPad Prism 6.0 and WinNonlin automatically excluded the zero value as outlier). The IVART, as described by Woodrow and colleagues [10], was designed primarily for ≥48 h in vitro re-invasion assays, rather than the schizont maturation assay used in the current analysis. IVART generates a growth ratio, defined as the uninhibited growth divided by maximally inhibited growth. In the schizont maturation assay, maximum inhibition leads to zero schizont count, thus producing infinite values at these concentrations [10]. Nevertheless, IVART could still be successfully applied to the schizont maturation data, producing similar estimates as the other programs.

Agreement with NONMEM
Bland-Altman plots were produced to measure the agreement between the different programs and NON-MEM for every anti-malarial according to the species of infection. Overall, there was good concordance between all programs and NONMEM. The biases (mean differences in estimates) were close to zero for all drugs assessed (range 3.80-1.76 nM; Additional files 1, 2, 3, 4 and 5). In total, 80 out of 1485 (5.4 %) assays had differences outside the 95 % limits of agreement of the mean difference and were termed outliers. The majority of these outliers (30/80, 37.5 %) were derived from assays with a large discrepancy in duplicate microscopy readings, possibly caused by human error, leading to responses with biphasic curves [9,25]. In these outliers, the difference in IC 50 values with the NONMEM value was higher in GraphPad Prism 6.0 compared to other methods. Data with a very steep slope (i.e., slope estimate = 10) accounted for 30 % (24/80) of outliers. For all these outliers, NONMEM consistently produced lower slope values compared to other methods, ranging from 2.9 to 4.9 (Table 4). For the assays with steep slopes, GraphPad Prism 6.0, IVART and ICE tended to produce higher IC 50 estimates than NONMEM. In total, 23.8 % (19/80) of the outliers had parasite growth at the highest concentration of the drug tested (i.e., 100 % inhibition was never reached). Whereas NONMEM interpreted these data within the context of the whole population, for the other programs, E min was anchored to zero, assuming zero growth was eventually reached with very high drug concentrations. However, this approach needs to be interpreted with caution. Data with 'plateau responses' (i.e., the same schizont counts were recorded for two or more different drug concentrations before it reached the maximum effect) were present in 5 % (4/80) assay outliers. The rest of the outliers (3.7 %; 3/80) were generated by IVART where a very wide confidence interval for the IC 50 values was generated.

Qualitative analysis
The key features of each package are summarized in Table 3.

Data conversion and transformation
All of the programs had individual template requirements. To fulfil these, raw data were converted into an Excel file first before transposing into the appropriate template. This process was critical and error-prone for WinNonlin, ICE and GraphPad Prism 6.0, where raw data needed to be entered manually in either row or column format. Data conversion for HN-NonLin and IVART was easier, since the data could be transferred in a 96-well format, reducing the risk for transcription errors. The in vitro assay relies on a process of normalizing the response at each concentration against a control well with no drug exposure. With the exception of WinNonlin, the normalization algorithm is included in the analysis for all the statistical packages assessed in this study.

Plate format and templates
Drug susceptibility assays generally use drug plates in a 96-well format. Although the eight-well rows of a 96-well plate are sufficient to generate valid dose-response curves and IC 50 estimates [4], most investigators prefer to use the 12-well columns of a 96-well plate in order to have a broader drug concentration range. In WinNonlin, GraphPad Prism 6.0, HN-NonLin, and IVART, users can enter pre-made drug concentration templates prior to processing, without the need to reset the template for every single analysis. This feature also allows the analysis of multiple drugs for the same isolate. This is not possible with ICE, which requires the drug concentration and response to be entered manually for each drug per isolate.
The number of drugs that can be assessed also varies. GraphPad Prism 6.0 allows the user to enter a userdefined number of drug concentrations as long the drug concentrations entered into the program are logtransformed. The drug template in HN-NonLin allows a maximum of 12 different drugs with seven different concentrations and one control well in a 96-well format. In IVART, the layout is more flexible and allows testing of eight to 12 different drugs with seven different concentrations plus one control well per drug, or 11 different concentrations and one control well per drug. Different 96-well plate templates employed by these programs can be suited to accommodate the investigator's needs. With horizontal drug dilutions, investigators are able to obtain more data points and a wider dynamic range with a specific drug; this is particularly useful when testing novel anti-malarials with unknown IC 50 s. With vertical dilutions, investigators are able to test more anti-malarials in parallel, but with a lower dynamic range, an approach more suited for anti-malarials for which the IC 50 range is well defined.

High throughput analysis
The ability to perform multiple plate analysis provides significant advantages for high throughput processing of large numbers of samples and a range of different antimalarials. This feature is only available for IVART and GraphPad Prism 6.0.

Data processing time
Outputs from all programs can be obtained as an Excel worksheet (HN-NonLin), a csv file (ICE and IVART), pdf file (IVART), or can be exported to Word or Excel (Win-Nonlin and GraphPad Prism 6.0). The process from raw data conversion until final outputs was timed for each package.
Since some programs are capable of analysing multiple drug plates and/or isolates, timing was performed by running a single isolate for five drugs on two 96-well plates (=single run). A single run could be completed on average in 40 s with GraphPad Prism 6.0 since multiple drug plates can be analysed in parallel and the results are updated instantly when new data are entered. HN-NonLin took 45 s, with the delay attributable to data from different plates having to be entered separately. Results from IVART can be obtained in approximately 90 s because raw data need to be uploaded to, and results downloaded from a server. A single run for ICE took 4-5 min, the main reason being that response data need to be entered per drug and per isolate.

Offline capabilities
For investigators working in resource-limited field settings with unreliable internet connection, the capability of offline data analysis has significant advantages. IVART and ICE require internet connection for uploading raw data to a server and calculations are generated online. WinNonlin, GraphPad Prism 6.0 and HN-NonLin allow offline data analysis.

Conclusions
The lack of standardized approaches to the analysis of in vitro drug susceptibility continues to confound antimalarial resistance surveillance and comparison of data generated from different laboratories. Reassuringly, the current study demonstrates that the various statistical programs produce good overall concordance for IC 50 estimates across different drugs and for both Plasmodium species tested. However, there were some important differences. In 5.4 % of drug assays, there were significant differences between approaches. In general, these were attributable to how the packages dealt with noisy data.
Since serial dilutions of drug concentrations are almost universally adopted, the confidence of estimates for highly resistant isolates can be wide. No matter which analytical platform is used, these important isolates need to be individually scrutinized. For routine and medium to high throughput drug susceptibility testing, GraphPad Prism 6.0 and IVART offered distinct advantages, including their capability to process multiple plates in parallel and the availability of all necessary output parameters to assess the accuracy of data modelling and hence, facilitate interpretation of the final data. For smaller scale drug testing, all statistical packages tested generated comparable IC 50 estimates. Investigators working in remote areas with poorer infrastructure, including limited internet connection, may prefer the freely available HN-NonLin for data analysis.