An R Package for Assessing Drug Synergism/Antagonism

Synergistic and antagonistic drug interactions are important to consider when devel-oping mixtures of anticancer or other types of drugs. Boik, Newman, and Boik (2008) proposed the MixLow method as an alternative to the Median-Eﬀect method of Chou and Talalay (1984) for estimating drug interaction indices. One advantage of the MixLow method is that the nonlinear mixed-eﬀects model used to estimate parameters of concentration-response curves can provide more accurate parameter estimates than the log linearization and least-squares analysis used in the Median-Eﬀect method. This paper introduces the mixlow package in R , an implementation of the MixLow method. Results are reported for a small simulation study.


Introduction
In many branches of medicine, drugs are administered in the form of mixtures. In order to develop new and more effective mixtures it is useful to analytically assess, with respect to the mixture's intended effect, the interactions that occur within a mixture. Such drug interactions can be of an additive nature, or can be antagonistic or synergistic (sub-additive or supraadditive, respectively). The approaches used to assess drug interactions vary depending on the design of the experiment, the expected shape of the concentration-response curve, and other factors. For reviews see Berenbaum (1989); Greco, Bravo, and Parsons (1995); Merlin (1994); Tallarida (2001). The approach described here is applicable to the common in-vitro situation where within-unit and between-unit measurements are available, responses follow a sigmoidal pattern, and ratios between drugs in a mixture are fixed (that is, various dilutions of the mixture and its component drugs are tested). While such data could be generated in many types of experiments, the data discussed in this paper are obtained from in-vitro cytotoxicity experiments, where cancer cells are exposed to a drug for a specified length of time (typically 72 hours) and then cell viability is indirectly measured, usually via fluorescence readings after addition of a suitable dye. Such cytotoxicity assays use multi-well incubation trays, where each tray receives one drug or mixture, each column of the tray might receive a different drug concentration, and replicate trays are tested for each drug and mixture. Responses are typically modeled as a sigmoidal function of the drug concentration.
Several indices for assessing drug interactions have been proposed, but perhaps the most popular one is the Loewe index. The Loewe index for a mixture is a function of the mixture's concentration-response curve as well as the curves of its component drugs. Therefore, estimating the parameters of these curves is a first step in estimating the index. A common method to estimate both the curve parameters and the index is the Median-Effect method developed by Chou and Talalay (1984). This method linearizes the sigmoidal response curve and then estimates concentration-response parameters using ordinary least squares regression. The interaction index is estimated based on regression results. The Median-Effect method has been criticized on several points, including the linearization step and the need for extensive preprocessing of data Greco et al. (1995). Recently,  proposed the MixLow method, which is an improvement over the Median-Effect method. The MixLow method utilizes a nonlinear mixed-effects model to more accurately estimate concentration-response parameters. In addition, data preprocessing is minimal.
The three components of the MixLow method are a nonlinear mixed-effects model, the Loewe index, and a method to calculate confidence intervals for the index. See  for complete details. In some applications it is not necessary to estimate a Loewe index and only curve parameter estimates are desired (IC50 values, for example). In such cases, estimation of the Loewe index can be omitted. Note that the R package drm by Ritz and Streibig (2005) also analyzes concentration-response data, but does not use a mixed-effects model to do so.
Concentration-effect data are modeled using a sigmoidal function. The utility of the sigmoidal function stems from its simplicity (as few as two parameters) and its empirical usefulness in modelling cytotoxicity data. For these reasons, sigmoidal functions are often used for modelling cytotoxicity data. The MixLow method is novel in its use of a nonlinear mixedeffects model for estimating sigmoidal curve parameters from concentration-response data, for the purpose of estimating drug interaction indices and associated confidence intervals. This paper describes the mixlow package, an implementation of the MixLow method in R (R Development Core Team 2009). An overview of the mathematical model is given, components of the package are summarized, and results from a small simulation study are reported. The package is available from the Comprehensive R Archive Network at http://CRAN.R-project. org/package=mixlow.

MixLow method
Let the random variable fa signify the fraction of cells affected by a drug concentration, and define φ = E[fa]. In some contexts, φ is estimated based on the data and in other contexts a concentration is estimated that produces a fixed value of φ. Denote the φ-effective log concentration of drug d by ψ d,φ . This is the concentration that produces a fraction affected equal to a fixed φ. For example, exp(ψ d,0.2 ) is the concentration that inhibits proliferation of a population by 20 percent. By convention, this is referred to as the IC20.

Basic MixLow model
The MixLow method utilizes a nonlinear mixed-effects model to estimate parameters of the sigmoidal concentration-response curve. Responses {Y d,t,w } are modeled as: where and the subscripts d, t, w refer to the d-th drug, t-th tray, and w-th well. In addition, γ d indexes the steepness of the curve for drug d at the IC50, and c d,t,w is the log concentration of drug d in well w and tray t. Here drug d could refer to a specific mixture or could refer to a component drug. Model (1) can be used to estimate parameters for single drugs or a mixture, or can be used to estimate parameters simultaneously for all component drugs and the mixture. For treatment control wells (those that receive media, cells, and a drug concentration of zero), φ d,t,w = 0. For other wells, 0 < φ d,t,w < 1. The expected value of exp(µ + b t ) refers to the expected mean response of control wells over all trays assessed, where b t is a random deviate specific to tray t (i.e., a "tray" effect). Values {b t } are independently distributed as b t ∼ N(0, σ 2 b ). The error terms are independently distributed as d,t,w ∼ N(0, f (·)), where f (·) is a function discussed later.
The Loewe index for a mixture of n drugs is given by: where m d,φ is an unknown constant signifying the log concentration of drug d in the mixture when the mixture is at its φ-effective log concentration, τ d is the fraction of the mixture composed of drug d, and ψ m,φ is the φ-effective log concentration of the mixture. The mixture is synergistic if L φ < 1, it is additive if L φ = 1, and it is antagonistic if L φ > 1. EstimateL φ is obtained by using estimatesψ m,φ andψ d,φ in Equation (3). To obtain expressions forψ m,φ andψ d,φ , note that if c d,t,w in (2) is equated toψ d,φ , thenφ d,t,w becomes φ. That is, To write the Loewe index as an explicit function ofψ d,φ , first solve (4) forψ d,φ to obtain Using (3) and (5), the estimator of the index is given bŷ Confidence intervals based on SE log(L φ ) are obtained using the Delta method as follows: where the superscript refers to transpose. The termv ar ψ γ is obtained from the observed information matrix, which is produced when Model (1) is fit by maximizing the likelihood function. The confidence interval forL φ is where t df,1− α 2 is a multiplier obtained from the t distribution with df degrees of freedom and α significance level.
Equation (8) is a piecewise confidence interval. In a typical linear regression setting, a simultaneous interval can be constructed using the Working-Hotelling procedure. If directly applied to the Loewe index, the Working-Hotelling interval would be where p is the number of parameters used to construct the interval and df is the degrees of freedom. The interval (9) has unknown properties, however. Model (1) is nonlinear and contains random effects. The standard error is obtained at φ using a first-order approximation to a linear model. Nevertheless, we can use the Working-Hotelling procedure to obtain a rough estimate of the difference in widths between the pointwise and continuous intervals. In a typical two-drug experiment, the interval will be based on six parameters (ψ and γ for each drug and the mixture). With three replicate 96-well trays per drug (see discussion below), the degrees of freedom will be roughly 540. Therefore, the multiplier in the continuous 95 percent interval will be larger than that in the piecewise interval by a factor of roughly or about 1.81.

Modified MixLow model
In some data sets, responses follow a sigmoidal curve but the curve appears to have a non-zero asymptote. An example is illustrated in Figure 1. The solid curve in the figure is the sum of Figure 1: Example concentration-response curve with non-zero asymptote.
the dotted sigmoidal curve and the offset from zero. This type of pattern can be modeled by altering Model (1) to the following: where Here, 0 ≤ λ d ≤ 0.5 is a drug-dependent weight, and φ , ψ , and γ all refer to the sigmoidal curve with zero asymptote. When λ d = 0, Model (11) collapses to Model (1). To write the Loewe index as a function of φ,ψ d,0.5 , andγ d , an expression forψ d,φ is obtained in a manner similar to (5). Specifically,ψ Using Model (11), and without additional assumptions, the Loewe index can only be estimated reads data from file mixlowData <-prepareData(trayData) prepare raw data cellLines <-getCellLines(mixlowData) get a vector of cell line names drugs <-getDrugs(mixlowData) get a vector of drug names trays <-getTrays(mixlowData) get a vector of tray names parameterDefaults <-getNlsParameterDefaults(trays) make data frame of parameter default values for NLS nlsData <-doNls(mixlowData, parameterDefaults) run nonlinear least-squares (NLS) model NlmePrintVarFunctions() print variance function options for nonlinear mixed-effects (NLME) model nlmeData <-doNlme(mixlowData, nlsData) run NLME model loeweData <-doLoewe(mixlowData, nlmeData) run Loewe analysis Table 1: Functions in the mixlow package in order of use.
In Figure 1, the maximal response (E[exp(µ |b t )]) is 900, the half-maximal response (at concentration ψ 0.5 ) is 450, and the asymptote occurs at a response of 100. The total response (solid curve) is equal to the asymptotic response plus the sigmoidal response (dotted curve). The half-maximal response of the dotted curve is 400, which occurs at a concentration ψ 0.5 .

Summary of the mixlow package
The mixlow package contains ten functions, listed in Table 1 and ordered by chronological use. In addition, there are plot, print, and summary methods for objects trayData, mixlowData, nlsData, nlmeData, and/or loeweData. Detailed information on the functions is available via the help commands in R. A summary of their use and actions is provided below. To use the mixlow package, the following six steps are performed in sequence: 1. Design the experiment, collect the data, and prepare the data file.
2. Read the data file using the function readDataFile.
3. Prepare and subset the data using the function prepareData.
4. Obtain starting values for the fixed-effects parameters of the nonlinear mixed-effects (NLME) model by using nonlinear least squares (NLS) analysis. This is done using the function doNls.
5. Use the NLS parameters as starting values for the NLME model. The mixed-effects model is estimated using the function doNlme.
6. Use the NLME parameters to estimate the Loewe index. This is done using the function doLoewe.

Experimental design and preparation of the data file
The appropriate experimental design for use with the mixlow package is the "ray" design, which has the following two characteristics: all drugs are tested individually at a variety of concentrations a mixture is made of two or more drugs and this mixture is tested at a variety of concentrations.
As a result, ratios between component drugs of a mixture are constant over all concentrations tested. The MixLow method will not work for checkerboard designs where ratios between drugs in a mixture vary. Each tray should test only one drug or one mixture, and each drug or mixture should be tested in replicate trays.
In a typical synergism experiment, two drugs and their mixture might be tested, although the MixLow model can be used on larger mixtures. Three replicate trays might be used for each drug and the mixture, or nine trays in total. Within each tray, 12 of the wells might be treatment-control wells that receive cells and media but no drug, and 60 might be treatment wells that receive drug, media, and cells. Ten different concentrations of the drug might be tested in the treatment wells (6 wells per concentration). Lastly, 24 of wells might be optical control wells (called blanks here for convenience). These can be of two types: Type bbt: optical control wells contain only media Type bbc: optical control wells contain media plus drug, and each drug concentration tested in treatment wells is also tested in optical control wells.
The first type is referred to here as bbt, or blanks-by-tray, and the second type is referred to as bbc, or blanks-by-concentration. The use of blanks-by-concentration is recommended when drugs can induce concentration-dependent responses (e.g., autofluorescence). For blanks-bytray, treatment-and treatment-control-well responses are adjusted by subtracting the mean response of all bbt wells in a tray. For blanks-by-concentration, a 1st to 4th degree polynomial is fit to the bbc responses and predicted values for each concentration are subtracted from treatment-and treatment-control-well responses.
The only preprocessing of data used in the Mixlow method is subtraction of an optical correction factor for each well based on bbt or bbc responses. Relative to the variance of treatmentwell responses at each drug concentration, the variance of optical control well responses is typically very small. Therefore, only a few bbt wells are needed per tray, or a few bbc wells are needed per concentration. Because of the low variance, use of optical correction factors is not likely to introduce significant uncertainty into the model.
In many cytotoxicity experiments, a reasonable arrangement for a 96-well tray is to use the last two rows of every column as bbc wells, use the first six rows of the first two columns as treatment-control wells, and use the first six rows of the third to twelfth columns as treatment wells, with each column receiving a different drug concentration. This layout is shown in  (e.g., responses in the center wells are different from responses in the peripheral wells), then a randomized concentration assignment should be employed.
To use the readDataFile function, data must be stored in a specially formatted ASCII file. An industry standard XML format for cytotoxicity data has not yet been developed, and the interim format described here was designed to facilitate cut-and-paste operations from plate reader data displayed in spreadsheet format. If a XML standard format is developed in the future, then the package can be modified to take advantage of this. Two sample data files are included with the package and can be used as templates for constructing new data files, as well as for testing out the capabilities of the program. One of these files contains simulated data, and one contains data from a synergism experiment described in .
The data file can be created in a spreadsheet (in which case it must be saved in tab-delimited plain text format without automatic enclosing quotations). Row labels (the first column in a data file if a spreadsheet is used) are required, and the order of row labels must not change from that given in the sample files. The required layout for the first three columns is shown in Table 3. The first eight row labels (up to and including col) occur only once. Other row labels are repeated for each tray. To conserve space, only three rows and two columns are shown for the conc, label, and resp entries. Entries for the following labels are optional, but all row labels, including the space labels, must still be present: global_notes, name, Institution, email, assay, date, drug_name_full, seed_density, assay_time, tray_notes, media, drug_solvent, solvent_drug_ratio.
The first eight rows contain general information about the experiments and the remainder of the file occurs in blocks, one block for each tray. Entries for tray_label and drug_name_short should contain only digits, letters, and the underscore character. The drug_name_short entry should consist of a unique, abbreviated drug name containing no more than 20 characters. This will be the drug name that appears in plots and in text output. The composition entry specifies the fraction of each drug in a mixture. For each drug, a string such as "drug_A = 0.5" must be given (for a hypothetical drug_A). The values provided should sum to 1.0 over all drugs, so if a single drug is tested in a tray the string would look like "drug_A = 1". The conc entries specify the drug concentration in each well and the label entries specify the type of well (rx for treatment, bbt or bbc for blanks). Treatment-control wells are also given the label rx. The resp entries specify the responses measured in each well.
The user does not need to manually alter the raw data in any way. The only exception is when there is good reason to believe that data from certain wells in a tray are anomalies. In any other notes on this tray media type of media used drug_solvent solvent used (DMSO, for example) solvent_drug_ratio ratio of solvent to drug space conc drug conc. in row 1, col 1 drug conc. in row 1, col 2 conc drug conc. in row 2, col 1 drug conc. in row 2, col 2 conc drug conc. in row 3, col 1 drug conc. in row 3, col 2 label label for row 1, col 1 label for row 1, col 2 label label for row 2, col 1 label for row 2, col 2 label label for row 3, col 1 label for row 3, col 2 resp response in row 1, col 1 response in row 1, col 2 resp response in row 2, col 1 response in row 2, col 2 resp response in row 3, col 1 response in row 3, col 2 this case the responses for these wells can be replaced by periods (".") in the resp entries of the data file, which will cause these wells to be ignored. If desired, the excludeWells argument of the readDataFile function can be used to exclude a list of wells from all trays. The row and column of each cell to be excluded need to be specified in the excludeWell list.

Generating starting values for fixed-effects parameters
Once a file is read using the readDataFile function, and data are subset and prepared using the prepareData function, the next step is to generate rough estimates of concentrationresponse curve parameters (µ, ψ d,0.5 , log γ d , and λ d ). These rough estimates will be used as starting values for fixed effects parameters in the nonlinear mixed-effects model. The estimates are generated using the doNls function, which is a wrapper for R's nls function.
The doNls function also produces output that can be plotted, and these plots can be useful for understanding the data and estimating parameter values graphically. An nls analysis is performed separately on each tray.
By default, only Model (1) is estimated. If the user provides a starting value for λ d , then Model (11) will also be estimated. If Model (11) is estimated, the results for Model (1) will be compared to those for Model (11) and the one that produced the lowest BIC score will automatically be selected. If Model (11) is estimated, the user can specify a threshold value (lower bound) for λ d . If the nls function estimates a value for λ d that is below this threshold, Model (1) will be selected. Final parameter estimates are obtained from the selected model.
The nls function requires starting values. As a convenience, the getNlsParameterDefaults function can be used to generate a data frame of parameter defaults for each drug. Parameter values of NaN indicate the following default values will be used: ψ d,0.5 will be set at the closest tested log concentration corresponding to 1 /2 the mean of control wells std(ψ d,0.5 ) will be set at the absolute value of ψ d,0.5/2 log(γ d ) will be set to zero std(log(γ d ) will be set to one µ will be set to the log of the mean of all control wells λ d will be set to NaN, which indicates that a lambda term will not be used. If a lambda term is used, its value should be set between zero and 0.5.
In order to increase the chance of model convergence, the doNls function will sample ψ d,0.5 and γ d starting values from a normal distribution with mean and standard deviation equal to the user-supplied or default-generated values. The user can specify the number of random samples (numRandomTries) to draw. After fitting, the models are ranked by BIC scores. The starting values for the top few models are averaged and used as starting values for a final model.
If verbose is set equal to TRUE, information on the various models constructed will be printed, including a summary table that contains estimates for γ d , ψ d,0.5 , µ, and λ d . These parameters are noted as g1, p1, u1, and lambda in the printout, respectively. (Values for λ d will be included only if a λ-model was estimated). Parameters γ d , ψ d,0.5 , and µ are in log scale. The print and summary methods for the data objects returned from the prepareData, doNls, doNlme, and doLoewe functions can also be used to examine the results.

Running the mixed-effects model
The function doNlme reads output from the NLS model and runs the NLME model. The doNlme function is a convenience wrapper for the nlme function from the nlme package Pinheiro, Bates, DebRoy, Sarkar, and R Development Core Team (2009). The user supplies an analysis argument, which can be either single or multiple. If single, each drug will be Variance Function Analysis Type d,t,w ∼ N (o, σ 2 d,t,w ), where σ d,t,w equals: analyzed separately. If multiple, all drugs will be analyzed together. The user must also specify the variance function to be used from the list given in Table 4. If analysis is multiple, the user can specify more than one variance function in the form of a vector, in which case one model will be estimated for each variance function and the best model will be identified using the BIC criterion. If analysis is single, the user must provide a list, named by drug, which specifies the vector of variance functions to be used for each drug. For a given drug, if more than one variance function is specified, one model will be estimated for each variance function and the best model will be identified using the BIC criterion. For heteroscedastic errors, variance function 3 (multiple drugs) or 2 (single drugs) may be appropriate in many cases. Simpler error functions could be tried if these do not allow convergence.
As with the doNls function, if the verbose argument is set equal to TRUE, information on the models constructed will be printed. Parameter values and standard errors are printed, along with other information. The rc values given at the bottom of the summary for each set of drugs is a goodness-of-fit statistic similar to the r 2 -statistic (it is the average model concordance correlation Vonesh, Chinchilli, and Pu 1996).
Lastly, for testing purposes, if only one tray of data is available for a drug, the function will duplicate concentrations and responses in that tray to make two (equivalent) replicates so that analysis can be completed. This can be used, for example, to obtain rough IC50 estimates when only one tray is available.

Loewe analysis
After parameter estimates for the concentration-response curves have been made by the doNlme function, these can be passed to the doLoewe function so that Loewe indices can be estimated. The primary result from the Loewe analysis is a data frame of fraction affected values and their corresponding Loewe index estimates. Standard errors and upper and lower confidence limits for the index at each fraction affected value are also produced. By default, fraction affected values range from 0.02 to 0.98 (in steps of 0.01).
The output of the doLoewe function also includes a measure of the degree of statistically significant synergism and antagonism at each fraction affected value. For synergism, this value is zero at any fraction affected value where the upper confidence interval limit is greater than one. Otherwise, the value is one minus the upper confidence interval limit. For antagonism, this value is zero at any fraction affected value where the lower confidence interval limit is less than one. Otherwise, the value is the lower confidence interval limit minus one. These values can be useful for summarizing the index over a range of fraction affected values. For an example application, see .

Analysis of example data
The following code illustrates use of mixlow functions on a data set included with the package. These data are for vincristine (vin), topotecan (topo), and their 1:1 mixture (mixture) tested against A549 human lung cancer cells. The response curve for vincristine suggests a non-zero asymptote.

Simulation study
Mixtures that contain more than a few drugs can be problematic in that the mixed-effects model may not converge or may take a long time to converge if all drugs and the mixture are analyzed together. If data are poorly behaved, convergence problems could occur for any size mixture. If convergence becomes a problem, one alternative is to estimate the concentrationresponse curve parameters separately for each drug and mixture. The user can instruct the doNlme function to perform this type of single analysis. The disadvantage of the singledrug approach is that control-well information is obtained from the trays of one drug only. In contrast, control-well data from all trays are used when all drugs and the mixture are analyzed together (multiple drug analysis). To determine if single-drug analysis makes a large difference in coverage of the confidence intervals of the Loewe index, a simulation was performed that compared the two approaches.
The simulation is based on Simulation C discussed in . Model (1) was used to generate the simulation data using the parameters listed in Table 5. These parameters are typical for cytotoxicity assays. One thousand experiments were simulated. Each experiment involved three drugs tested in three trays, or nine trays in total. Each tray used 14 control  wells, two bbt wells, and 80 treatment wells (8 treatment wells for each of 10 drug concentrations). The experiment simulated a sham mixture of a drug with itself. Therefore, the parameters used to generate data for drug 1, drug 2, and drug 3 were identical. Data for drug 3 were taken to represent the sham mixture of a 50/50 mixture of drugs 1 and 2. As a sham mixture, the true Loewe index is 1.0 for all fraction affected values.
Average coverage of the confidence intervals for the Loewe index over all 1,000 experiments is shown in Figure 3  for all fraction affected values, regardless of the type of analysis used. At fraction affected values near 0.58, coverage for multiple-drug analysis was closer to nominal values than for single-drug analysis (0.950 vs. 0.938, respectively). In conclusion, both types of analysis produced excellent coverage. Multiple-drug analysis is recommended, however, if convergence is not an issue.

Summary
This paper has introduced the mixlow package in R. The package should be useful for students, researchers, and technicians who need to obtain accurate parameter estimates for sigmoidal concentration-response curves based on in-vitro experiments. In addition, it should be useful to those who must estimate a drug interaction index to determine if a mixture is antagonistic, additive, or synergistic.