Measuring biological activity

The measure of effect in our screening is usually the inhibition of cellular response relative to the untreated level (vehicle alone). For untreated vehicle and treated levels V and T, we calculate a fractional inhibition I = 1–T/V. The inhibition ranges from 0% at the untreated level to 100% when T = 0. Inhibition levels are negative for agents that actually increase levels. Other effect measures, such as an activity ratio r = T/V may be more appropriate for some assays. When activity ratios (eg, fold increase over stimulated control) are being used, the effect can be measured using an induction I = ln(T/V). With this definition, all effect expressions are the same as for inhibition.

For oncology screens with time zero reference growth levels, we use a "Growth Inhibition" (GI) based on the measures proposed for the NCI60 screens [Shoemaker'06]. We define GI = 1 – [ T<V0 ? (T V0)/V0 : (T V0)/(V V0) ], using the conditional syntax [ condition ? result.if.true : result.if.false ], and where V0 is the vehicle-treated level at time zero. The descriptive levels for this measure are "silent" (GI=0), "active" (GI=0.5), "static" (GI=1), "toxic" (GI=1.5), and "lethal" (GI=2). As defined for the NCI60 screens [Shoemaker'06], the GI50, TGI, and LC50 are the concentrations that achieve active, static, and toxic levels respectively. Some treatments which actually increase growth will have T > V, with a corresponding "agonizing" reference level (GI= 0.5).

Error estimates for these measures result from comparing each treated level T to the median untreated level V ± σV, determined for each plate by finding the median level among the untreated control wells arranged across the plate. For example, applying error propagation σ2(f(xi)) = σ2(xi) (df/dx¬i)2 to the expression for fractional inhibition I, the estimated standard error is σI ~ (σV/V) sqrt(1–I). The error estimates are further increased to account for variations between replicate combination blocks as well as a minimum assumed fractional uncertainty of σmin ~ 3%. Thus for inhibition, the standard error estimate becomes σI ~ sqrt[ (σV/V)2 (1–I) + σrep2 + σmin2 ].

We use medians rather than averages to reduce the effect of occasional outliers on the consensus. Note that while medians are robust to outliers, they are more sensitive to statistical noise, yielding ~30% larger deviations for normally-distributed data. Standard deviations are estimated from the median absolute deviation (MAD), where for a normal distribution, the sample deviation σdat ~ 1.5 MAD [Filliben'05]. The error for the median is then σmed ~ σdat/sqrt(N 1), given N data values.

Identifying replicate dose series and dose matrices

Our cHTS screening produces phenotypic measurements for single agents or fixed ratio combinations as a series of concentrations (doses) or for pairwise combinations as dose matrices which contain all mixtures of of two single agents each sampled at a series of concentrations, including zero [Borisy '03; Lehár '09]. The information from replicate doses can be summarized by generating a consensus block of data (either dose series or dose matrix) out of the individual blocks from the replicates.

In the Chalice Analyzer, single agent dose series are identified by specifying the desired chemical’s name or C number, along with other search criteria specifying the assay details or other properties of the data set. The Analyzer will then collect all examples of the requested chemical that were tested as single agents and will either present them as separate blocks or merge them according to instructions from the user.

Combination data are identified by specifying pairs of chemical names or C numbers, along with other search criteria. Each combination dose matrix contains internal copies of the single agent curves which are used as the reference for combination effects. The properties of these single agent dose series are used to guide the merging process between replicate data blocks.

Merging replicate dose series

The single agent curve data are recorded at different doses (chemical concentrations) and arranged into a series with increasing dose. To merge dose series from replicate data sets, these series are resampled to generate a consensus dose series which will be used for calculating the median dose response values. Two algorithms are currently employed, basic and periodic (Fourier) resampling. For both algorithms, the data from all individual dose series are collected into a single data set and sorted by dose.

In basic resampling, a sliding window of a specified width in log-concentration (logarithm of the tested concentration) space is moved through the tested doses. Starting with the lowest nonzero tested concentration in the series, all data points within the specified window are merged, and the lowest tested data point outside the window serves as the starting dose for the next window. The window is advanced iteratively through the dose series until all data points have been merged, and for each such window the average log-concentration within the window is defined as the consensus dose. These consensus doses are then reconverted to concentrations through an inverse logarithm operation.

In periodic resampling, a Fourier algorithm is employed to define a regular series of doses with a fixed dilution factor covering the tested doses, phased to best align the consensus doses with the original sampling. For a set of N doses x covering M distinct log-concentrations (logarithm of the tested concentration), with δx and ∆x as the smallest and largest intervals between the nonzero sampled doses, a frequency interval δf is set to be 1/(2∆x). Then for all frequencies f covering the range from 2Mδf to 1/(δx) in steps of δf, the Fourier coefficients A(f) = ∑i cos(2πfxi) and B(f) = ∑i sin(2πfxi) are summed across all N doses xi, and the peak frequency fpeak is identified with the largest Fourier power A2+B2. The period Ppeak = 1/fpeak and optimal phase φpeak = Ppeak atan2(Bpeak,Apeak)/2π are calculated, and optimal doses ξk = φpeak + kPpeak are chosen for all integers k yielding ξk within the tested dose range. All data within a dosing window of ±Ppeak/2 from each ξk are then merged and assigned a consensus dose ξk, which are then reconverted to concentration through an inverse logarithm operation. Dosing windows in the periodic series that contain no data are assigned null response values.

After resampling, all data points within each consensus window are merged by taking the median of the corresponding response values, and assigning the consensus dose for that window. Error estimates are calculated for each consensus response value following the error propagation methods described in "Measuring biological activity" above.

Merging replicate dose matrices

Pairwise combinations are organized into dose matrices which contain all mixtures of two single agents each sampled at a series of concentrations, including zero [Borisy '03; Lehár '09]. When many such dose matrix blocks are collected for merging, the consensus block needs to account for the possibility of different single agent sampling between data blocks.

The process for merging dose matrices is as follows: (1) first, the single agents are aligned (put on the same vertical or horizontal axis) across all of the data blocks; (2) next, each single agent’s dose series is merged according to the single agent merging algorithms described above, yielding a separate consensus dose series for each single agent; (3) finally, all the combination points within each two-dimensional window defined by each of the single agents’ concentration windows are merged and assigned consensus concentration values corresponding to the single agent consensus doses. All data points within a dosing window are merged by taking the median data value with its associated uncertainty as described above in "Measuring biological activity". For those dosing points defined by the single agent sampling windows that contain no combination data, null response values are assigned.

Sigmoidal response models for Dose Response data

The single agent activity is characterized by fitting a sigmoidal function of the form I = ImaxCα/[Cα+EC50α]. Here, C is the concentration, EC50 is the EC50, and α is the sigmoidicity. We also use these parameters to calculate an inhibitory concentration ICcut, where the fitted curve crosses a chosen effect level Icut. The fitting is performed using a parameter search algorithm based on cost minimization. The cost function is a chi-squared between the fitted sigmoid and the measured data points with standard errors, with the last data point having double weight to ensure that compounds with activity only at the highest concentration are represented. Additional penalties are applied to restrict parameter ranges. Using fractional inhibitions, for example, hard limits on Imax are enforced by multiplying the cost function by ∆10, where ∆ is (Imax–Ilim)/0.001, where Ilim is the upper or lower limit allowable for that measure. For EC50, a softer multiplicative penalty of ∆2 is applied, where ∆ is the logarithm of the ratio between the fitted EC50 and the largest or smallest nonzero concentration in the data series. The sigmoidicity values are restricted between 1 and 10 using a penalty of ∆2, where ∆ is the distance from α outside the nearest limit. The curve fits were optimized using simplex cost minimization [Nelder&Mead'65; Press '97], with the starting simplex vectors connecting the center of each parameter’s range to its upper limit.

Combination effect reference models

Combination effects can be most readily characterized by comparing each data point’s inhibition to that of a combination reference model that was derived from the single agent curves [Greco '95]. Three models are generally used: (1) The highest single agent model IHSA(CX,CY) = max(IX,IY) is a simple reference model, where CX,Y are the concentrations of the X and Y compound, and IX,Y are the inhibitions of the single agents at CX,Y; (2) Bliss independence [Bliss'39] IBliss(CX,CY) = IX + IY – IXIY represents the statistical expectation for independent competing inhibitors; and (3) Loewe additivity [Loewe'28], where ILoewe(CX,CY) is the inhibition that satisfies (CX/ICx) + (CY/ICy) = 1, and ICx,y are the inhibitory concentrations at ILoewe for the single agent curves. Loewe additivity is the generally accepted reference for synergy [Greco '95], as it represents the combination response generated if X and Y are the same compound. Both IHSA and IBliss are easily calculated from IX,Y, but determining ILoewe requires interpolation and numerical root finding [Berenbaum'85].

To select promising candidates from large combination screens, it is most useful to design a numerical score that captures the qualities of desirable combinations, which can be used to prioritize the manual selection process. These scores should be designed for each screen to capture the experimental objectives, but there are usually three important considerations: (1) significant synergy over the Loewe additive model; (2) substantial activity where the synergy occurs; and (3) sufficient potency shifting.

Activity over Loewe additivity is most easily calculated using a simple volume score, where VLoewe = log fX log fY ∑ (Idata–ILoewe), summed over all non-single agent concentration pairs and where log fX,Y are the natural logarithm of the dilution factors used for each single agent. This effectively calculates a volume between the measured and Loewe additive response surfaces, corrected for varying dilution factors. This volume score emphasizes the overall synergistic or antagonistic effect of the combination, thus minimizing the effects of outlying data spikes and identifying combinations with a robust synergy across a wide range of concentrations and at high effect levels. VLoewe is positive for mostly synergistic combinations and negative for antagonism. The uncertainty σV can be calculated based on the measured errors σI and standard error propagation.

We usually use a "Synergy Score" S = fcov ln fX ln fY ∑ max(0,Idata) max(0,Idata–ILoewe), which is a positive-gated, inhibition-weighted volume over Loewe additivity. This provides an additional prioritization favoring combinations whose synergy occurs at high effect levels, ignoring antagonistic portions of the response surface. Here fX,Y are the the dilution factors used for each single agent and the coverage factor fcov accounts for missing data, scaling the score up by the ratio of total/tested combination dose matrix points. S is always positive, and its uncertainty σS can be calculated based on the measured errors σI and standard error propagation. An alternative to the synergy score is the "Hit Score" H = fCOV log fX log fY ∑ max(0,Idata) max(0,Idata–IHSA), which refers to the HSA model. The key distinctions between S and H lie in the different underlying models and also in how the single agents are used in the model calculations. In the Chalice Analyzer, the HSA model is calculated directly from the single agent responses at corresponding concentrations, while the Loewe additive model is derived from the sigmoidal fits to the single agent response curves.

To prioritize hits, we generally use distributions of our preferred score (eg, S) and its error to define an appropriate selection cutoff. For example, combinations with S > 3σS are "individually significant" at ~99% confidence, assuming normal errors. To estimate systematic experimental errors that are not tested by replicate plates, we sometimes use the distribution of synergy scores for any drug-with-itself combinations acquired during the experiment to determine a plausible range for non-detections. Alternatively, the score distribution for the whole experiment could be used to identify outliers at a chosen confidence level.

Displaying potency shifting using isobolograms

Potency shifting is usually shown using an isobologram [Greco '95] which shows how much less drug is required in combination to achieve a desired effect level, when compared to the single agent doses needed to reach that effect.

The choice of effect level for the isobologram display and combination index calculations is somewhat arbitrary, and can either be manually or automatically selected in the Analyzer. The automatic iso-level selection algorithm finds the observed Idata with the the largest Idata–ILoewe, excluding those points with Idata exceeding the lesser single agent’s Imax. This exclusion is applied to ensure that the isobologram reflects the best synergy at levels covered by both single agents. The chosen iso-level is then floored at a 20% inhibition level, beneath which inhibitions often are unreliable, and rounded to the nearest 5%.

Having selected an isobologram level Icut, the isobologram is drawn by identifying the locus of concentrations that correspond to crossing the chosen iso-level. In the Analyzer, this is done by finding the crossing point for each single agent concentration in a dose matrix across the concentrations of the other single agent. Practically, first each vertical concentration CY is held fixed while a bisection algorithm [Press '97] is used to identify the horizontal concentration CX in combination with that vertical dose that gives the chosen effect level in the response surface Z(CX,CY). These concentrations are then connected by linear interpolation to generate the isobologram display. The uncertainty for each crossing point is estimated from the response errors using bisection to find the concentrations where Z–σZ(CX,CY) and Z+σZ(CX,CY) cross Icut.

Scoring potency shifting using combination index

Potency shifting is scored a combination index [Chou&Talalay'84] CI. For a chosen iso-effect level Icut, CII = (CX/ECX)I + (CY/ECY)I, where (CX/ECX)I for a particular data point is the ratio of the X compound’s measured concentration to its effective concentration at the chosen inhibition level. The CI is a rough estimate of how much drug was needed in combination relative to the single agent doses required to achieve the chosen effect level. CI values in the range of 0.5-0.7 are typical for in vitro measurements of current clinical combinations [Greco '95]. The combination index error σCI is calculated using standard error propagation through the CI calculation based on the isobologram errors derived as described above in "Displaying potency shifting using isobolograms".

In the Chalice Analyzer, the best CI is reported from the many combination index values calculated for each Icut crossing concentration. Among all the measured CI values, the one with the largest signal-to-noise level (1-CI)/σCI is reported as the best combination index.

References

1. Berenbaum (1985). The expected effect of a combination of agents: the general solution. J Theor Biol 114(3): 413-31. #4021503
2. Bliss (1939). The toxicity of poisons applied jointly. Ann. Appl. Biol. 26: 585-615.
3. Borisy, Elliott, Hurst, Lee, Lehár, Price, Serbedzija, Zimmermann, Foley, Stockwell, Keith (2003). Systematic discovery of multicomponent therapeutics. PNAS 100(13): 7977-82. #12799470
4. Chou, Talalay (1984). Quantitative analysis of dose-effect relationships: the combined effects of multiple drugs or enzyme inhibitors. Adv Enzyme Regul 22: 27-55. #6382953
5. Filliben. (2005, 2005). Exploratory Data Analysis. Engineering Statistics Handbook Retrieved 10/14, 2007, from http://www.itl.nist.gov/div898/handbook/eda/section3/eda356.htm.
6. Greco, Bravo, Parsons (1995). The search for synergy: a critical review from a response surface perspective. Pharmacol Rev 47(2): 331-85. #7568331
7. Lehár, Krueger, Avery, Heilbut, Johansen, Price, Rickles, Short, Staunton, Jin, Lee, Zimmermann, Borisy (2009). Synergistic drug combinations tend to improve therapeutically relevant selectivity. Nat Biotechnol 27(7): 659-66. #19581876
8. Loewe (1928). Die quantitativen Probleme der Pharmakologie. Ergebn. Physiol. 27: 47-187.
9. Nelder, Mead (1965). A Simplex Method for Function Minimization. The Computer Journal 7(4): 308-313.
10. Press, Teukolsky, Vetterling, Flannery (1997). Numerical Recipes in C: the Art of Scientific Computing, 2nd Ed. Cambridge, CUP.
11. (2006). The NCI60 human tumor cell line anticancer drug screen. Nat Rev Cancer 6(10): 813-23. #16990858