Growth rate inhibition metrics correct for confounders in measuring sensitivity to cancer drugs.

Drug sensitivity and resistance are conventionally quantified by IC50 or Emax values, but these metrics are highly sensitive to the number of divisions taking place over the course of a response assay. The dependency of IC50 and Emax on division rate creates artefactual correlations between genotype and drug sensitivity, while obscuring valuable biological insights and interfering with biomarker discovery. We derive alternative small molecule drug-response metrics that are insensitive to division number. These are based on estimation of the magnitude of drug-induced growth rate inhibition (GR) using endpoint or time-course assays. We show that GR50 and GRmax are superior to conventional metrics for assessing the effects of small molecule drugs in dividing cells. Moreover, adopting GR metrics requires only modest changes in experimental protocols. We expect GR metrics to improve the study of cell signaling and growth using small molecules and biologics and to facilitate the discovery of drug-response biomarkers and the identification of drugs effective against specific patient-derived tumor cells.

INTRODUCTION

The quantification of drug response is fundamental to discovering therapeutic molecules, investigating their mechanisms of action[1-3] and studying signal transduction, cell division and other biological processes using chemical biology approaches[4,5]. In the case of anti-cancer drugs, cells are typically exposed to drugs over a range of concentrations, and the number of viable cells (or surrogates such as ATP using CellTiter-Glo®) are measured several days later. Data comprising cell counts in the presence of drug divided by counts for untreated controls are fitted to a sigmoidal curve to compute (i) the concentration of drug at which the cell count is half the control (IC), (ii) the fraction of viable cells at the highest drug concentration (E), and (iii) the area under the dose-response curve (AUC)[6,7]. Dose-response and genomic datasets are often combined to discover drug response biomarkers[1-3,8,9], but it has recently been found that large-scale drug response data vary from one study to the next[10] for reasons that remain poorly understood[11-13].

We show here that for dividing cells, traditional drug response metrics such as IC suffer from a fundamental flaw when they are estimated from cell counts made at the end of the experiment (the standard approach): if cells undergo different numbers of divisions during the course of an assay due to natural differences in proliferation rate, variation in growth conditions, or changes in the duration of an experiment, IC, and AUC values will vary dramatically, independent of any changes in the underlying biology. Thus biomarkers that predict sensitivity under one (potentially arbitrary) set of assay conditions may not predict sensitivity under slightly different conditions. We therefore propose a new method for parameterizing drug response, the normalized growth rate inhibition (GR), which is based on comparing growth rates in the presence and absence of drug. Parameterization of GR data yields GR, GR, GR, and h (Hill Slope) values that are largely independent of cell division rate and assay duration (we use “area over the curve”, GR, rather than AUC for reasons discussed in online methods). GR metrics can be determined with modest changes in experimental procedures and we propose that they replace IC and E values in assessing cellular response to drugs, RNAi, and other perturbations in which control cells divide over the course of the assay.

RESULTS

Definition of normalized growth rate inhibition (GR)

We used computer simulation to model drug response by three idealized cell lines having identical sensitivity to a cytostatic drug (i.e. a drug that arrests but does not kill cells) and different division times (T = 1.8, 2.4 or 3.9 d). These division times correspond to the lower quartile, median, and upper quartile for breast cancer cell lines[3] and are similar to those of NCI-60 cells[14]. In the slowly dividing cell line (T = 3.9 d), the total number of cells did not double in a typical three-day assay, thus E ≥ 0.5 and IC was undefined. In the case of the two faster growing cell lines, IC and E values fell as division rate increased (Fig. 1a) simply because cell number (or CTG value) was normalized to a drug-naïve control in which cell number increases as division time fell (compare curves in panels of Fig. 1a).

Figure 1

Modeling drug response and the dependence of drug response metrics on division time

(a) Simulation of a simple drug-response model yields relative cell counts across a concentration range for a cytostatic drug for a slow- (left), medium- (middle), and fast-growing cell line (right; Td: division time). Black lines correspond to untreated control samples and red lines denote 50% growth inhibition. Black marks show where IC and E are evaluated. (b–d) Methods for evaluating GR value: (b) Conceptual approach based on growth rates (k and k(c)), (c) fixed interval approach based either on cell number at the start (x) and end of the experiment (x and x(c)), or (d) time-dependent value based on cell count before and after a time interval 2 × Δt (x(c, t ± Δt)). (e) Simulated data showing relative cell count (green lines) and GR value (purple lines) for a cytostatic drug assayed over three days. The darker the line the longer the division time. IC and GR are projected onto the x-axis, E and GR onto the y-axis. (f) IC or E (green) and GR or GR (purple) computed from a theoretical three-day assay with cells dividing at different rates (AUC and GR values in Supplementary Fig. 1c).

We can compensate for the confounding effects of division rate on drug response measurements by computing “normalized growth rate inhibition” at time t in the presence of drug at concentration c: where k(c,t) is the growth rate of drug-treated cells and k(0) is the growth rate of untreated control cells (Fig. 1b). The GR value is simply the ratio between growth rates under treated and untreated conditions normalized to a single cell division. The sign of the GR value relates directly to response phenotype: it lies between 0 and 1 in case of partial growth inhibition, equals 0 in the case of complete cytostasis and lies between 0 and −1 in case of cell death. Given GR values for a range of drug concentrations, GR is the concentration at which GR(c) = 0.5, GR is the maximal measured GR value, and h is the slope of the sigmoidal fit; GR is calculated by integrating the GR curve over a range of concentrations (see online methods).

In practice, GR values can be estimated from endpoint measurement of cell number in treated and treated samples, given the initial cell number (Fig. 1c; this is related to the procedure for GI50 determination, see Supplementary Note 1). Alternatively, the doubling time for untreated cells can be measured under the same conditions in parallel experiments and used in place of the initial cell number (see online methods). A time-dependent GR value can be evaluated given cell count measurements at two or more time points. Time-dependent GR values capture adaptive responses, varying kinetics of drug-target interaction, drug efflux, etc (Fig. 1d). Introducing time as a variable makes it possible to relate drug-induced changes in cell states to dynamic measures of drug response at a molecular level (equations for all calculations are provided in online methods with links to scripts).

To compare GR dose-response curves to conventional curves we created synthetic data for cells with T ranging from one to four days exposed to a drug that is partially cytostatic, fully cytostasic, or cytotoxic (models are described in online Methods). For all three drugs IC and E values were strongly correlated with division time and assay duration, but this was not true of corresponding GR metrics (Fig. 1e–f, Supplementary Fig. 1a–b and methods). In the case of drugs that kill cells rapidly, independently of cell cycle state, GR still varies with growth rate and time-dependent GR values are preferable (Supplementary Fig. 2 and methods). AUC combines IC or E data and is often less sensitive to experimental noise[1,6,11] but it suffers from the same dependence on cell division time; GR corrects for this (Supplementary Fig. 1c).

GR metrics are robust to experimental variability

To study how changing cell division affects IC values we expressed BRAFV600E in hTERT-immortalized retinal pigment epithelial (RPE) cells and then exposed them to etoposide, a topoisomerase II inhibitor that has a cytostatic effect in RPE cells (and whose mechanism of action is independent of BRAF). Oncogene over-expression is known to slow the growth of non-transformed cells[15,16] and expression of BRAFV600E in RPE cells under the control of a doxycycline-inducible promoter was observed to increase division time 3-fold as the doxycycline concentration increased from 0 to 60 nM. Under these conditions the estimated IC value for etoposide increased 100-fold and E increased from 0.25 to 0.6 (Fig. 2a and Supplementary Fig. 3a). In contrast, GR and GR values varied only slightly under these conditions because the effects of etoposide per cell division were unchanged.

Figure 2

GR values are independent of both the length of the assay and the division time

(a) Effect of altering division time of hTERT RPE-1 cells on metrics of etoposide sensitivity. Relationship between division time (T) and doxycycline (DOX) concentration (left) to cells expressing BRAFV600E under the control of a DOX-regulated promoter. Values for IC and GR (middle), and E and GR (right) were evaluated at 48 h at different concentration of DOX. Large and undefined values for IC were set at 30 µM for purposes of illustration. (b) Effect of altering division time of MCF 10A cells on metrics of etoposide sensitivity. Relationship between division time (T) and concentration of EGF in a serum-free media (left). Values for IC and GR (middle), and E and GR (right) were evaluated at 59 h at different concentrations of EGF in a serum-free media. Large and undefined values for IC were set at 30 µM for purposes of illustration. (c) Evaluation of relative cell count (top) and GR values (bottom)for a drug concentration close to the GR value (left) and computed response metrics at different time points (middle and right) as estimated from live-cell imaging of MCF 10A cells exposed to one of five drugs with different mechanisms of action. (d) Time-dependent GR values estimated over 18 h intervals for MCF 10A cells grown as spheroids treated with omipalisib (left panel) and corresponding time-dependent GR values (right). Dashed line shows the GR value evaluated at 90 h. In all panels, the data shown derive from one of three biological replicates.

In a second experiment, we measured etoposide sensitivity in MCF 10A cells grown in serum-free medium supplemented with different concentrations of epidermal growth factor (EGF). The division time of MCF 10A cells varied from >10 d at 0.02 ng ml−1 EGF to <2 d at 1 ng ml−1 EGF (Fig. 2c). Across this range, IC and E for etoposide were substantially more variable than GR or GR values (Fig. 2b and Supplementary Fig. 3b). Thus, potentially arbitrary and unintended differences in cell culture conditions can lead to large variation in IC or E values that are unlikely to have a biological basis.

To test the GR approach in a typical small-scale drug sensitivity screen we exposed MCF 10A and BT-20 cells, expressing H2B-mCherry to facilitate automated cell counting, to five different drugs chosen for diverse mechanisms of action. We measured cell number every ~8 h over 3 d using an automated microscope. Estimated IC and E values converged only after three divisions (~60 h) whereas GR metrics stabilized by the first division (<20 h for MCF 10A and ~36 h for BT-20; Fig. 2c and Supplementary Fig. 3c). This confirms that GR metrics are substantially less dependent on assay duration than IC and E. In the case of very slow and uneven growth, by primary human tumor cells for example, the stabilization of GR values within one cell division is likely to be a real advantage in obtaining reliable estimates of drug sensitivity.

When grown in 3-D culture, MCF 10A cells are known to adapt to PI3K/mTOR inhibitors (such as omipalisib), becoming less sensitive over time[17]. End-point measures of drug sensitivity do not report on such adaptive responses. However, when we monitored spheroids by live-cell imaging, we observed that the time-dependent GR and GR values were lowest ~20 h after omipalisib addition and increased ~10-fold by day 4 (Fig. 2d); endpoint GR values lay midway in this range (Fig. 2d). Thus, time-dependent GR data directly capture the decreasing effectiveness of omipalisib in MCF 10A spheroids, enabling detection and further analysis of adaptive mechanisms.

Analysis of high-throughput data using GR metrics

The majority of large-scale drug response datasets published to date neither report cell division rates nor allow us to estimate them. An exception is a study by Heiser et al.[3], which recorded cell numbers for breast cancer cell lines before and after three-day exposure to a panel of anti-cancer drugs. For many drugs in this dataset, IC correlates with division rate[6] (e.g. for cell cycle inhibitors, regression coefficient of –0.54, Spearman’s p-value < 10−66, N=2956, Supplementary Fig. 4a). However, we could reproduce this correlation using an idealized model that does not assume any biological connection between drug sensitivity and division rate and by repeatedly simulating drug responses using random parameters (Supplementary Fig. 4b and methods). This finding suggests that the correlation between drug sensitivity and division rate found in experimental data is spurious; the correlation was also absent when drug response was measured using GR (Spearman’s p-value = 0.31, Supplementary Data 1) suggesting that it is an artifact of the way IC is calculated.

Paclitaxel, a microtubule inhibitor widely used in chemotherapy, is one drug exhibiting a strong negative correlation with division rate and also substantial variation (100-fold) in IC across cell lines (Fig. 3a). This relationship has been previously described and is thought to arise because paclitaxel acts primarily on mitotic cells and the faster a cell line divides, the more likely it is to be in mitosis[18,19]. However, reanalysis of data in Heiser et al.[3], show that GR values for paclitaxel actually span a narrow range centered around 10 nM[14] (Fig. 3b and Supplementary Fig. 4c), close to the estimated affinity of paclitaxel for assembled microtubules in vitro[20]. In contrast, GR values for paclitaxel vary considerably across cell lines and subtypes. In HER2-amplified (HER2amp) and triple-negative breast cancer (TNBC) lines GR values are negative, which is indicative of a cytotoxic response, whereas hormone receptor-positive (HR+) and non-malignant lines generally exhibit cytostatic responses (Supplementary Fig. 4c). This is consistent with data showing that HER2amp and TNBC (basal-like) human tumors are more taxane-sensitive that other types of breast cancer[21]. We propose that future attempts to find biomarkers predictive of paclitaxel response focus on variation in GR rather than IC[22,23].

Figure 3

Evaluation of GR metrics in a high-throughput dataset

(a–b) Fitted dose-response curves for (a) relative cell count and (b) GR values for paclitaxel in breast cancer cell lines (3-day assay data from Heiser et al.[3]). Red denotes cytotoxic response and blue denotes cytostatic response; darker curves denote cell lines with fewer divisions. Marginal distributions (below and to the right of the dose-response plots) show the relation between estimated IC or GR and E or GR values and number of divisions for that line over three days; ρ values are Spearman’s correlation coefficients. (c) Number of divisions over three days for cell lines in the dataset grouped by subtype: HER2-amplified (HER2amp), triple-negative breast cancer (TNBC), hormone receptor-positive (HR+), and non-malignant (NM). (d–e) Distribution of IC50 (green) and GR50 (purple) for all drugs in Heiser et al.[3] (d) or for ErbB2 inhibitors only (e) grouped by clinical subtype. P-values were derived from a ranksum test for IC or GR values for NM cells (e) or for HER2amp cells (f) lines versus all other subtypes.

The average division rate of breast cancer cell lines in culture differs with clinical subtype: among cells studied by Heiser et al.[3], HR+ and HER2amp subtypes divided most slowly (median T = 3.2 d), TNBC cells divided faster (median T = 2.3 d) and non-malignant cells were faster again (median T = 1.8 d, Fig. 3c). As measured by IC values, non-malignant cells were on average more sensitive to anti-cancer drugs than tumor cells, whereas GR values showed that the mean and range of drug sensitivity was similar (Fig. 3d). Focusing on HER2amp lines, IC values for EGFR/ErbB2 inhibitors were similar across breast cancer subtypes, even though HER2amp human tumors are preferentially sensitive to such drugs and ErbB2 inhibitors are front-line therapy for this disease[24,25]. GR data for HER2amp cell lines in Heiser et al. show that this subtype of breast cancer is in fact ~10-fold more sensitive to EGFR/ErbB2 inhibitors in vitro than other subtypes (Fig. 3e; p = 1.3×10−4). The failure of IC values to show the preferential sensitivity of HER2amp cell lines to EGFR/ErbB2 inhibitors arises because their relatively slow growth rate is a hidden confounder in IC calculation. From these data we conclude that artefactual dependency of IC on cell division rate obscures meaningful associations between genotype and drug sensitivity and also create associations where none exist.

Effects of cell density on drug sensitivity

GR metrics allow us to quantify how drug sensitivity changes in the face of variables that affect division rate. One such variable is seeding density: increasing density has widely been reported to promote drug resistance[26-29]. To investigate this, we cultured six breast cell lines representative of different subtypes at six seeding densities over a 32-fold range. Twenty-four hours after plating, cells were exposed to eleven drugs with diverse mechanism of action. Growth rates, IC, E, and GR metrics were estimated by imaging and counting fixed cells 72 h after treatment (Supplementary Data 2). Overall, IC and E correlated positively with seeding number for MCF 10A, MDA-MB-231, and Hs 578T cells and negatively for BT-20 and SK-BR-3 cells (p < 0.01, Fig. 4a). In the first three cell lines, division rate decreased with density, presumably due to contact inhibition, depletion of essential media components, and other effects (Fig. 4b)[30-32]. In BT-20 and SK-BR-3 cells, division rate increased with density, presumably due to conditioning of the growth medium[31,33]. Correlations between density and IC and E were the weakest for MCF7 cells, which grew equally well across all densities. Thus, the effect of density on number of divisions, and consequently on IC and E, varied dramatically with cell line (Supplementary Fig. 5). However, only a small number of cell line-drug pairs exhibited a statistically significant association between plating density and GR or GR values. We wondered whether these were situations in which the biology of drug response was altered by density.

Figure 4

Plating density affects division rate and drug sensitivity

(a) Spearman’s correlation between estimated values for IC, E, GR, and GR and seeding number for the six breast cancer cell lines shown. Data derive from an experiment in which cells were plated at a range of six densities and treated with eleven drugs having diverse mechanisms of action (see Online Methods and data 2). Significance: * for p<0.05, ** for p<0.01, and *** for p<0.001. (b) Relationship between seeding number and the number of divisions after 3 days. Error bars are the SEM of three biological replicates.

In MCF 10A cells, six drugs were associated with a significant variation in GR or GR (or both) across seeding densities; RNA-seq revealed that gene expression in these cells also varied with cell density with significant enrichment for genes involved in catabolic processes, cellular respiration, and wounding responses[30-32,34] (Supplementary Fig. 6a). When the data was analyzed by principle component analysis, PC1 (which captured 32% of variance) was strongly correlated with the number of cells at the time of collection (Spearman’s ρ = 0.98) and less so with time in culture (ρ= 0.57) suggesting that cell density at the time of assay and not the history of the culture was the primary determinant of transcriptional state (Supplementary Fig. 6b). Follow up studies with methotrexate showed that the ratio between cell number and media volume was the key variable for sensitivity to this drug, not cell density per se (this was also true of oligomycin, an inhibitor of ATP Synthase; Fig. 5a and Supplementary Fig. 7); since cells grew in a constant volume of media, GR50 was therefore time-dependent. Density (and time) dependent variation in GR50 was also observed for linsitinib, an IGF1R inhibitor in Phase II clinical trials (Fig. 5b). Variation in GR50 with time and density was reduced by co-treatment with the metalloprotease inhibitor batimastat, suggesting a role for autocrine conditioning of the microenvironment in drug response.

Figure 5

Time-dependent GR metrics reveal diverse mechanisms of drug sensitivity and resistance

Evaluation of three selected drugs for which seeding density had a significant effect on GR values in MCF 10A cells based on data summarized in Fig. 4. (a) Influence of ratio of cell number to media volume for drugs related to metabolism. Time-dependent GR values for MCF 10A cells at different seeding numbers following exposure to 0.1 µM methotrexate (MTX) (left). Time-dependent GR values for MCF 10A cells for MTX or oligomycin (middle); values are capped at 30 µM for illustration purposes. GR dose-response curves at 3 days for MCF 10A treated with MTX with a constant number of seeded cells but in different volumes of media (right). Data shown derive from one of two biological replicates. (b) Media conditioning and adaptive resistance to linsitinib, an IGF1R inhibitor. GR dose-response curves at 3 days for MCF 10A treated with linsitinib at different cell densities (left). Error bars are the SEM of three biological replicates. Time-dependent GR values for MCF 10A at different seeding numbers following exposure to 20 µM linsitinib in the absence (middle) or presence (right) of 10 µM batimastat, a matrix metalloproteinase inhibitor. Data shown are from one of three biological replicates. (c) Paclitaxel-induced cytostasis at low cell density and apoptosis at high density. GR values for MCF 10A cells seeded at different densities and treated with paclitaxel for 3 days (left). Error bars are the SEM of three biological replicates. Fraction of MCF 10A cells positive for cleaved caspase 3 at different time points and cell densities following exposure to 0.03 µM paclitaxel (middle). Time-dependent GRmax values for paclitaxel in MCF 10A cells seeded at different densities (right). Data shown are from one of two biological replicates.

In the case of paclitaxel, GR values remained at ~5–10 nM across plating densities but GR was strongly density-dependent, varying from ~0 (cytostatic) at low cell densities to negative values (cytotoxic) at higher densities (Fig. 5c, left). Time-dependent GR reached its greatest negative value at 24 h concomitant with the an increase in the fraction of taxol-treated cells that contain cleaved caspase-3, confirming that negative GR values corresponded to elevated apoptosis (Fig. 5c, middle and right). The molecular basis of this effect is unknown, but it may be one reason for the discrepancy between cell-killing dynamics in low-density culture and in xenograft tumors in which density is high[35].

DISCUSSION

Accurate measurement of drug sensitivity and resistance is the cornerstone of cancer biology, pharmacology, and many fundamental studies on cell signaling and cell division. In this paper we demonstrate theoretically and experimentally that variation in division rate seriously confounds existing drug response metrics. We also show that division rate varies with cell type, media composition and seeding density, often in unpredictable ways. Cell division rate slows down in some cell lines as density increases while speeding up in others. Such variation can change apparent IC 100-fold or more and therefore introduces artificial correlations in data, obscures the true effects of drug action, and introduces unknown complications into biomarker discovery. As an alternative we propose GR metrics that are computed by comparing growth rates in the presence and absence of drug. GR and GR are robust to variation in cell division rate and should replace IC and E in studies in which control cells divide (including the study of drugs, gene depletion/overexpression and variation in the extracellular environment). GR quantifies the potency of a drug on a per division basis, ensuring that fast- and slow-growing cells with similar biochemical responses to drug are scored equivalently. GR captures the maximal drug effect on growth rate and differs from E in that it falls between 1 and –1, where negative values denote cell death, 0 denotes cytostasis, and positive values denote partial inhibition. GR, and h (Hill Slope) values can also be calculated from GR curves; the former is often the most robust metric in the face of experimental noise and the later quantifies an important relationship between dose and response that often neglected[6].

In the simplest version of our method, GR metrics are computed after measuring cell number or a surrogate (e.g. CTG) before and after exposure of a culture to varying concentrations of drug or other perturbation for a fixed time (Fig. 1c). When cell division rates under similar culture conditions are known from previous data, only the final cell number is needed (Fig. 1d; although we believe that before and after data on cell number is a valuable control to ensure data quality). Time-course data makes it possible to compute time-dependent GR values and quantify phenomena such as delayed response, drug adaptation, or variation in the kinetics of drug-target interaction (Fig. 1e). To facilitate the use of GR metrics by others we provide MATLAB and python routines and an online GR calculator (www.grcalculator.org). Moreover, the data in this paper, including the GR/GR values for the Heiser et al. dataset, are available online (http://lincs.hms.harvard.edu/).

Large-scale drug-response studies based on existing response metrics[1,2,8,9] are discrepant for poorly understood reasons[5,11,12], raising concerns about the value of drug response biomarkers[10]. We speculate that comparison of datasets across centers (or even within a center) might be confounded by differences in plating density, growth media, and other factors that affect cell division rate. It might be possible to correct for this in existing data by computing GR values post facto, but this will require recreating the original assay conditions followed by measurement of division rates.

Based on the results in this paper, we believe that use of GR metrics in lieu of traditional IC, E, or AUC values will improve our ability to identify genes and biological processes responsible for drug sensitivity and resistance. GR metrics decouple any effect that genotype or microenvironment have on division rate from their effects on drug sensitivity. Cell biology studies involving the modification of genes or the microenvironment often result in changes in cell division rate, leading to potentially spurious correlations with drug sensitivity (as illustrated here by oncogene overexpression and changes in EGF levels) and should also benefit from the use of GR metrics. By analogy with antimicrobial susceptibility testing for bacterial infections, it has recently been suggested that cancer therapy might be personalized by screening primary human tumor cells against panels of drugs[36,37]. Such cells grow slowly and unevenly in culture, making division number a poorly controlled variable. Accounting for such differences using GR metrics should create drug response data that are more reproducible and useful for optimizing patient therapy.

ONLINE METHODS

Metrics of drug response

Determining relative cell counts

Cell counts in the presence of drug are normalized to DMSO-treated controls grown on the same plate under the same conditions. In the current study, relevant conditions include seeding density, the concentration of exogenous growth factors (e.g. EGF in the case of MCF 10A cells), and the concentration of a second drug such as doxycycline or batimastat. For each cell line, drug, and drug concentration, we define the relative cell count as x(c)/x, where x(c) is the count in the presence of drug and x is the 50%-trimmed mean of the count for control cells. Technical replicates are averaged to yield an average relative cell count. We typically collect data from three technical replicates (on three separate plates).

Calculating GR values using end-point drug response data

Normalized growth rate inhibition is calculated according to the formula: where x(c) and x are as described above and x is the 50%-trimmed mean of the cell count from a sample grown in parallel and measured just prior to drug exposure.

Alternatively, the untreated division time T=ln(2)/k(0) can be measured in independent experiments and used in place of the initial cell number (x0 = x · 2−): where T is the duration of the assay.

Calculating time-dependent GR values

GR values can be evaluated over a time interval (2 × Δt) around any time point t based on the equation: The time-dependent GR values in the current paper were computed with 2 × Δt = 12 h to 18 h, which corresponds to about half a cell division time.

Curve fitting and estimating drug response metrics

GR data are fitted to a sigmoidal curve as follows (Supplementary Fig. 8a–c): where the fitted parameters are: If the fit of the curve is not significantly better than that of a flat curve (i.e. GR(c) ≡ GR) based on an F-test with cutoff of p = 0.05, the response is considered flat and the parameter GEC is set to 0 (Supplementary Fig. 8d).

: the effect of the drug at infinite concentration (GR ≡ GR(c → ∞)). GR lies between −1 and 1; negative values correspond to cytotoxic responses (i.e. induction of cell death, Supplementary Fig. 8a), a value of 0 corresponds to a fully cytostatic response (Supplementary Fig. 8b), and a positive value to partial growth inhibition (Supplementary Fig. 8c).

: the Hill coefficient of the fitted curve which reflects how steep the dose response curve is. In practice, we typically constrain h to a value between 0.1 and 5.

: the concentration at half-maximal effect. To avoid artefacts in curve fitting we constrain GEC to be two orders of magnitude higher and lower than the experimentally tested concentration range (in practice, this is usually about 10−7 to 103 µM).

Inferred drug response metrics

The value is the concentration of drug at which GR(c = GR) = 0.5. If the value for GR is above 0.5, the GR value is not defined and is therefore set to +∞ (Supplementary Fig. 8c). By extension other thresholds can be defined in a similar manner. For example GR corresponds to the concentration at which a drug is fully cytostatic: GR(c = GR) = 0.

The is the maximum effect of the drug at the highest tested concentration; it lies between −1 and 1; a value of 0 corresponds to a fully cytostatic response and negative values correspond to cytotoxic response. GR can be estimated from the fitted curve or obtained directly from experimental data; we often do the later.

For time course data, all metrics are evaluated at each time point individually.

Area under the curve and over the curve (GRAOC)

Another common metric for quantifying dose response is the area under the response curve (AUC), which is based on integrating the dose response curve over the range of tested concentrations. In the case of GR curves, which can have negative values, is more intuitive to use the area over the curve: where GR(c) are measured GR values are discrete concentrations c. GR has the benefit that, in the case of no response, it has a value of 0. It is important to note that GR values (like conventional AUC) can only be used to compare responses evaluated across the same drug concentration range.

The GR value captures variation in potency and efficacy at the same time. The calculation of GR at discrete (experimentally determined) concentrations has the advantage that it does not require curve fitting and is therefore free of fitting artifacts. This is especially useful for assays where fewer than five concentrations are measured and curve fitting is unreliable. GR values are also more robust to experimental noise than metrics derived from curve fitting; GR values are particularly sensitive to outlier values when directly obtained from data.

Drug concentration range

The drug concentrations used for fitting drug response curves need to span a sufficiently wide range and have sufficiently intermediate values, to obtain reliable estimates for GEC, h, GR. Denser sampling provides more precise estimates, especially in the case of steep dose response curves. Optimal design of dose-response curves has been discussed elsewhere[7]. In practice we suggest using nine doses spanning four orders of magnitude from 1 nM to 10 µM. This range can be shifted to lower concentrations for more potent drugs, or higher values for less potent drugs with the caveat that GR values for different drugs can only be compared over the same concentration range. We suggest discarding any GR value which is more than an order of magnitude above the highest tested concentration because values extrapolated from the fitted curves are more subject to fitting artefacts than interpolated values. Similarly, GR value is not properly constrained if the GR(c) dose-response curve does not reach a plateau at the highest measured concentration. In such cases, GR is the most reliable metric.

Computing GR metrics

Source code for computing GR metrics are provided (Supplementary Software). To facilitate the computation of GR metrics we provide updated source code available under an open source software license and MATLAB and python scripts at https://github.com/sorgerlab/gr50_tools. We also provide an on-line calculator at www.grcalculator.org. This web site contains a user guide, various tutorials and explanatory materials and example datasets, including all of the data in the current manuscript.

Theoretical model of drug response

To simulate the effect of division time on GR and conventional drug response metrics under different assumptions about the degree of cytostasis or cell killing, we developed a theoretical model of drug response. To the first approximation, cell growth can be considered exponential, with drugs either decreasing the division rate or killing cells in a cell cycle-dependent manner: where x is the cell count, k is the untreated growth rate (per day), c is the drug concentration, S is the maximal inhibitory effect, SC is the concentration at half-maximal effect of drug, and h the Hill coefficient. The growth rate k corresponds to the division rate k as k = ln(2)·k = ln(2)/T, where T is the division time. S can be larger than 1 to account for drugs inducing cell death at a specific phase of the cell cycle. The model can also be generalized to account for drugs that induce cell death independent of the cell cycle: where k is the maximal killing rate (per day), LC is the concentration of drug that produces half-maximal cell killing. Integrating these equations for an assay of t days yields the cell count at concentration c: where x ≡ x(t = 0) is the cell number at the time of treatment. Thus, the relative cell count is: and the normalized growth rate inhibition (GR value) is: This equation for GR(c) is independent of the length of the assay t and thus, the metrics GR, GR, GR, and h are also independent of t. For cases where drug action is mainly related to the cell cycle (k =0), GR values are also independent of the untreated growth rate k. As shown in Supplementary Fig. 2 and 3a–b, the impact of k>0 is minimal on GR and relatively small on GR. This is also illustrated by the analytical formula for GRGR = 2^(1 − S − k/k) − 1. Note that the metrics derived from growth inhibition (GI) values used in Heiser et al.[3] such as GI are more robust than traditional metrics, but still depend on both division time and assay length (see Supplementary Note).

Parameters for simulations in figure panels

Model parameters used in the numerical simulations shown in Fig. 1b, Supplementary Fig. 1, 2, and 7 were as follows: For Supplementary Fig. 3b, the parameters were randomly sampled within the following distribution:

Division rate k: normal distribution around 0.9 divisions per day with standard deviation of 0.46 and a lower bound of 0.14.

Hill coefficient h: uniform distribution between 1.5 and 2.5

Maximum inhibition SM: uniform distribution between 0 and 2

Half inhibition concentration S50: log-uniform distribution ranging from 0.31 and 10

Maximum toxic effect TM: 0 for 50% of values, uniform distribution between 0 and 0.5 for the other 50% of values.

Half inhibition concentration T50: log-uniform distribution ranging from 0.56 and 5.6

Experimental methods and data processing

Cell lines and tissue culture

MCF 10A, Hs 578T, MDA-MB-231, MCF7, SK-BR-3, and BT-20 were obtained from the ATCC and grown according to ATCC recommendations. For time-lapse experiments MCF 10A and BT-20 cells were modified by inserting an H2B-mCherry expression cassette (gift of R Benezra, Addgene plasmid # 20972)[38] that comprised AAVS1 homology arms, the hPGK promotor, and SV40 polyA terminator (gift from R Jaenisch, Addgene plasmid # 22072)[39] into the AAVS1 safe harbor genomic locus using CRISPR/Cas9. MCF 10A-H2B-mCherry cells were grown in the same manner as the parental strain, with the exception that traditional DMEM was replaced with FluoroBrite DMEM (Thermo Fisher Scientific) for imaging. The modified hTERT RPE-1 cells (gift from J Chen) were created by inserting the full length BRAF(V600E) expression cassette (Addgene plasmid # 15269)[40] drive by a tet-inducible promotor (Addgene plasmid # 41394). Cell identity was confirmed by short tandem repeat (STR) profiling at the Dana-Farber Cancer Institute and all cells were tested with the MycoAlert PLUS mycoplasma detection kit (Lonza) and found to be free of Mycoplasma prior to analysis.

Drugs and dyes

Drugs were obtained from commercial vendors and tested for purity in-house as described in detail in the HMS LINCS drug collection database (http://lincs.hms.harvard.edu/db/sm/). Drugs and reporter dyes were dispensed directly onto multiwall plates using a D300 Digital Dispenser (Hewlett-Packard). To stain cells, YOYO-1 and TOTO-3 (Thermo Fisher Scientific) were used at 250 nM and 100 nM respectively; NucView 488 caspase 3 substrate (Biotium) was used at 200 nM.

Manipulating cell growth rate to determine effects on drug sensitivity

RPE-1 or MCF 10A-H2B-mCherry cells were plated in 384-well plates using the Multidrop Combi Reagent Dispenser (Thermo Scientific) at 250 and 500 cells per well respectively. To modulate the growth rate in RPE-1 cells, we induced expression the BRAF(V600E) oncogene by treating with indicated doses of doxycycline using a D300 Digital Dispenser (Hewlett-Packard). To modulate the growth rate in MCF 10A-H2B-mCherry we serum-starved cells twice with DMEM/F12 media supplemented with 0.1% bovine serum albumin and 1% penicillin-streptomycin. Media changes and cell washing was performed using a EL406 Microplate Washer Dispenser (BioTek). Cells were treated with indicated doses of human epidermal growth factor (EGF, eprotech) using a D300 Digital Dispenser. After 24 hours the cells were treated with a dilution series of etoposide using a D300 Digital Dispenser. RPE-1 cells were stained and fixed for analysis at the time of drug treatment and after 72 hours. MCF 10A-H2B-mCherry were imaged in an IncuCyte ZOOM live cell imager (Essen Bioscience) starting at the time of EGF treatment and drug sensitivity was evaluated 72 hours after drug addition.

Evaluating drug response metrics in MCF 10A and BT-20 over time

MCF10 A-H2B-mCherry and BT-20-H2B-mCherry at 1250 and 2500 cells per well respectively in 384-well plates using the Multidrop Combi Reagent Dispenser (Thermo Scientific) and grown for 24 hours. Cells were treated with a dilution series of the indicated drugs using a D300 Digital Dispenser (Hewlett-Packard) and imaged after drug addition in an Operetta (Perkin Elmer) for high content imaging system equipped with a live-cell chamber over a period of 96 hours. For these experiments, we used the following drugs:

Etoposide, Topiosomerase inhibitor

Linsitinib, IGF1R inhibitor

Omipalisib/GSK2126458, panPI3K/mTOR inhibitor

PLX4720, B-RAF inhibitor

Tanespimycin/17-AAG, HSP90 inhibitor

Evaluating drug sensitivity in MCF 10A spheroids over time

MCF10 A-H2B-mCherry cells were plated at 200 cells per well into ultra-low attachment coated, flat bottom, 384-well plates (Corning) with the addition of 2.5% growth factor reduced basement membrane matrix Matrigel (Corning) to the growth media. After 48 hours cells were treated with a dilution series of omipalisib using a D300 Digital Dispenser (Hewlett-Packard) and imaged in an IncuCyte ZOOM live cell imager (Essen Bioscience) for an additional 96 hours. Drug sensitivity was evaluated by computing GR values every 18 hours over a 10–90 hour time frame.

Evaluating drug sensitivity in breast cancer cells plated at different seeding densities

MCF 10A, Hs 578T, MDA-MB-231, MCF7, SK-BR-3, and BT-20 were plated at densities ranging from 156 to 5000 cells per well in 384-well plates using the Multidrop Combi Reagent Dispenser (Thermo Scientific) and grown for 24 hours. Cells were treated with a dilution series of the indicated drugs using a D300 Digital Dispenser (Hewlett-Packard). Cells were stained and fixed for analysis at the time of drug treatment and after 72 hours of incubation with drug. For these experiments, we used the following drugs:

Erlotinib, EGFR inhibitor

Etoposide, Topiosomerase inhibitor

Lapatinib, EGFR/ErbB2 inhibitor

Linsitinib, IGF1R inhibitor

Methotrexate, dihydrofolate reductase inhibitor

Omipalisib/GSK2126458, panPI3K/mTOR inhibitor

Paclitaxel, target microtubules

Palbociclib, CDK4/6 inhibitor

PLX4720, B-RAF inhibitor

TAE684, ALK inhibitor

Tanespimycin/17-AAG, HSP90 inhibitor

Investigating density-dependent drug effects

MCF 10A-H2B-mCherry cells were plated at densities that ranged from 156 to 5000 cells per well in 384-well plates using the Multidrop Combi Reagent Dispenser (Thermo Scientific) and grown for 24 hours. Cells were treated with a dilution series of drugs using a D300 Digital Dispenser (Hewlett-Packard) and imaged after drug addition in an Operetta (Perkin Elmer) for high content imaging system equipped with a live-cell chamber over a period of 72 hours.

In the case of methotrexate and oligomycin, 1250 cells were plated in 20–120 µl of media per well, treated with a dilution series of drug, and imaged for 72 hours.

In the case of linsitinib, cells were treated with a dilution series of linsitinib either with or without 10µM batimastat using a D300 Digital Dispenser and imaged in an IncuCyte ZOOM live cell imager (Essen Bioscience) for an additional 72 hours.

In the case of paclitaxel, cells were treated with a dilution series of paclitaxel and 200 nM of NucView 488 caspase 3 substrate (Biotium) using a D300 Digital Dispenser (Hewlett-Packard) and imaged after drug in an IncuCyte ZOOM live cell imager (Essen Bioscience) for an additional 72 hours. For immunofluorescence experiments, cells were grown for 24 hours and then treated with a dilution series of paclitaxel using a D300 Digital Dispenser (Hewlett-Packard) and incubated for 3, 6, 12, and 24 hours. Cells were fixed for 30 min in 3% formaldehyde, permeabilized for 30 min in phosphate buffered saline (PBS) with 0.3% Triton X-100 (Sigma-Aldrich), washed twice in PBS with 0.1% Tween 20 (Sigma-Aldrich; PBS-T), and blocked for 60 min with Odyssey blocking buffer. Anti-active Caspase-3 antibody (BD Biosciences) was diluted 1:1000 in Odyssey blocking buffer and incubated for 16 h at 4°C. Cells were washed three times in PBS-T for 5 min and incubated with Alexa Fluor 488 conjugated goat anti-rabbit secondary antibody for 60 min at room. Cells were washed two times in PBS-T, once with PBS, and stained for 30 min with whole cell stain (Thermo Fisher Scientific) and Hoechst (Thermo Fisher Scientific), and washed three times in PBS.

Fixed-cell endpoint assays

After drug treatment cells were stained at the indicated timepoints with 2 µM Hoechst 33342 (Thermo Fisher Scientific) and LIVE/DEAD Far Red Dead Cell Stain (Thermo Fisher Scientific) for 30 minutes and fixed with 3% formaldehyde (Sigma Aldrich) for 30 minutes. Fixed cells were imaged using an Operetta microscope and analyzed using the Columbus image data storage and analysis system (Perkin Elmer).

Live-cell time course assays

Cells expressing H2B-mCherry were imaged at the indicated timepoints after drug treatment in an Operetta high content imaging system equipped with a live-cell chamber or an IncuCyte ZOOM live cell imager. Dead cells were identified by counter-staining cells with YOYO-1 or TOTO-3 (Thermo Fisher Scientific) and apoptotic cells were identified with the NucView 488 caspase 3 substrate (Biotium). Live, dead, or apoptotic cell lines were identified using the Columbus image data storage and analysis system or the IncuCyte analysis software. Spheroid growth was estimated by measuring the sum of the spheroid area per well using IncuCyte analysis software.

mRNA analysis

Cells plated at different densities in 384-well plates were harvested at the indicated times and RNA was extracted using the RNeasy mini kit (Qiagen). To ensure sufficient amounts of RNA, wells having a low number of cells were pooled. Libraries were prepared by the Broad Technology Labs (BTL) following the protocol for SCRB-Seq described in[41]. Transcripts were quantified by the BTL computational pipeline using Cuffquant version 2.2.1[42]. Analysis of the transcriptional data was performed in MATLAB using standard libraries and in-house scripts. Gene set enrichment analysis was performed using the GSEA software v2.1.0 from the Broad Institute with the MSigDB v4.0 GO biological process set[43].

• Cytostatic:	SM=1, S50=1.5,	TM=0,	h=1.6
• Partial response:	SM=0.65, S50=1.2,	TM=0,	h=1.6
• Cytotoxic:	SM=2.6, S50=2,	TM=0,	h=1.6
• Partial toxic response:	SM=0.45, S50=1.2,	TM=0.05, T50=5,	h=1.6
• Mixed response:	SM=0.65, S50=1.2,	TM=0.1, T50=3,	h=1.6
• Complete cytotoxic:	SM=2.6, S50=1.2,	TM=1, T50=3,	h=1.6