Correction: Measurement and models accounting for cell death capture hidden variation in compound response

An amendment to this paper has been published and can be accessed via a link at the top of the paper.


Introduction
Quantifying cellular response to therapeutic compounds is essential to understanding their mechanisms of action and assessing therapeutic efficacy [1][2][3] . In the case of cancer treatments, and often with other diseases, drug activities are evaluated by quantifying the number of live cells after a short period using direct or surrogate measurements 4,5 . However, quantities beyond the number and viability of cells provide valuable information about the cellular response. Along with altering cell proliferation, promoting cell death is another important index of drug efficacy 6,7 . Incomplete eradication of drug-susceptible malignant cells allows the survival of drug-tolerant persister cell populations that can develop resistance by multiple routes [8][9][10] . Moreover, cell death can occur via a variety of mechanisms, including apoptosis and necroptosis, and selection among these outcomes can potently modulate cancer immunogenicity 11 . Limited understanding of these underlying cellular responses further complicates the assessment of drug combinations. Combination treatments are typically evaluated for their ability to enact greater effects than either compound alone 12 , but typically only by quantifying viability.
Here, we show that directly measuring both cell growth and death can provide valuable information for interpreting the response of cells to single and combination treatments. We propose a framework for quantifying drug response that accounts for the compound-induced changes in rates of cell growth and death. This approach reveals extensive differences in cell response, otherwise hidden by simply quantifying cell number. Of course, trade-offs exist for the breadth versus depth of analysis that can be performed to characterize cellcompound response. We show that end-point analysis preserves much of the distinct outcomes we observe for kinetic measurements while allowing similarly highthroughput analysis to those of live-cell number surrogates. These results demonstrate the need for and an approach to more precisely quantify the nature of cellcompound response and interactions.

Results
Viability alone is insufficient to distinguish cell growth and death effects To test whether growth and death are confounded in livecell measurements, we first explored the uncertainty in a model using only these measurements (Fig. 1a). We fit typical dose-response measurements of H1299 non-smallcell lung cancer cells to the chemotherapy doxorubicin (Fig. 1b) to a model incorporating both cell growth and death. We assumed no cell death in the absence of drug to show the best-case scenario of assessing drug response. This model was identifiable for the live-cell number (Fig. 1c), and the IC 50 and E max of compound effect on cell viability were narrowly defined as 21.4 ± 1.0 nM and 26.7 ± 3.4% (90% confidence interval, CI), respectively. In contrast, the model showed large uncertainty in inferred cellular growth or death rates (Fig. 1d, e). At the maximum dose, the predicted growth rate ranged 0.33-0.71 day −1 and death rate 0.02-0.40 day −1 (90% CI). The large uncertainty in outcome was due to a strong correlation in the fit values of drug effect on the growth and death rates (Fig. 1f). This shows that given only live-cell number one is incapable of distinguishing between reductions in cellular growth rate and increases in cell death. The number of divisions and cell deaths can vary largely while similarly fitting live-cell measurements. Moreover, the number of cell divisions and cumulative dead cells can differ drastically while resulting in the same cell viability (Fig. 1g).
High-throughput measurements of cell death quantify compound response To quantify pharmacologic response, we extended our experimental measurements to those of cell death. While quantifying the number of cells over time by phase, we used an Annexin V probe to measure phosphatidylserine exposure during apoptosis and a membrane-impermeable DNA dye, YOYO-3, to measure permeabilized apoptotic and dead cells (Fig. 2a) 13,14 . The areas occupied by cells, Annexin V, and YOYO-3 signal in each image were then analyzed to determine the total, apoptotic, and dead cells relative to the whole image area. We evaluated this realtime imaging method by measuring the response to doxorubicin (DOX) in H1299 cells. DOX strongly reduced the number of cells (Fig. 2b, top), as seen before (Fig. 1b). At the same time, we observed an increase in Annexin V while YOYO-3 increased minimally throughout the assay. Fitting these data to a model of cell growth and death (Fig. 2c), we observed a strong decrease in the inferred growth rate (div), and only a modest increase in cell death rate (deathRate; Fig. 2d, top). We next compared these measurements to those with another chemotherapy, vinorelbine (NVB), again observing a dose-dependent decrease in the number of live cells (Fig. 2b, bottom). At the same time, we observed a large increase in both Annexin V and YOYO-3 signal. This was reflected in our subsequent analysis, inferring an increase in the death rate (Fig. 2d, bottom). The fraction of cells dying through apoptosis (apopfrac) was also inferred to be lower in NVB as compared with DOX treatment (Fig.  2d). While kinetic measurements provide a wealth of information, end-point measurement is more amenable to high-throughput experiments. We confirmed that our analysis provided qualitatively similar results using only the first and last measurements in each experiment, demonstrating that cell growth and death can be quantified using either kinetic or endpoint measurements ( Supplementary  Fig. S1).
To independently verify these opposing outcomes upon DOX or NVB treatment, we used a carboxyfluorescein succinimidyl ester (CFSE)-based proliferation assay to verify the distinct growth rate effects inferred by our analysis. We measured CFSE intensity of untreated cells every 24 h starting from a day after cell labeling ( Fig. 2e, left). The detected intensity over time was used to estimate the number of times each cell had divided. Consistent with the inferred growth rates of our model, DOX-treated cells distributed more toward higher CFSE intensity in a dose-dependent manner, implying fewer cell divisions than with non-treated cells (Fig. 2e, right). In contrast, the distribution of NVB-treated cells remained more similar to non-treated cells. Our inferred cell growth rates and CFSE measurements were overall consistent (Fig. 2f).
To validate the inferred cell death rates, we measured the induction of cleaved caspase-3, an apoptotic marker, after treatment with DOX or NVB (Fig. 2g). After 24-h treatment, cleaved caspase-3 was detected in NVB-treated cells, but not DOX-treated cells. Both drugs induced the caspase-3 cleavage by 48-h treatment. These observations were consistent with the inferred cell death rates, indicating that NVB induces faster cell death than DOX. Interestingly, the level of cleaved caspase-3 after 72-h treatment was lower in NVB-treated cells compared with DOX-treated cells. This may be explained by the lower inferred apoptosis fraction in NVB treatment while the cell death rates increased ( Fig. 2d) implying an increase in non-apoptotic cell death, such as caspase-independent cell death. Alternatively, a large fraction of NVB-treated cells are dead after 72 h. Phosphatidylserine (to which Annexin V binds) is irreversibly externalized to the cell surface by caspase-dependent scramblases during apoptosis, in contrast, it is not exposed by caspase-independent cell death 15,16 , and NVB may induce caspase-independent cell death as well as apoptosis. DOX and NVB are known to operate through distinct mechanisms-inducing double-stranded breaks or preventing microtubule polymerization, respectively. As a result, each compound has a differing dependency on p53 status and leads to arrest in distinct cell cycle phases, supporting that each might engage distinct cell death programs [17][18][19] . Collectively, measuring and analyzing phase, Annexin V and YOYO-3 signals can quantify both the growth and death rate effects of drugs on cells.

Targeted compounds also display distinct phenotypic consequences
We next evaluated whether the growth-death model can dissect the response of cancer cells to targeted compounds as well. We treated a non-small cell lung cancer cell line PC9 with seven different targeted drugs whose effects on cell growth and death we expected to vary according to their mechanisms of action. We also used paclitaxel, a chemotherapeutic drug that is widely known to interfere with the cell cycle resulting in cell death and reduced proliferation, for comparison. The total, apoptotic and dead cell measurements from each compound treatment were diverse ( Supplementary Fig. S2), and these differences were reflected in the inferred cell division and death rates ( Fig. 3; Supplementary Fig. S3).
We were able to classify the tested compounds into three types: a compound that (1) both inhibits cell division and induces cell death, (2) only inhibits cell division, and (3) only induces cell death. As expected, paclitaxel fell into the first type by simultaneously exhibiting strong cell growth suppression and death. Similar to paclitaxel, the PI3Kɑ inhibitor BYL719, the pan-PIM kinase inhibitor PIM447, the EGFR tyrosine kinase inhibitor erlotinib, and CDK7/12 inhibitor THZ1 were grouped into the first type. In contrast, both OSI-906, a dual IGF-1R/INSR tyrosine kinase inhibitor, and binimetinib, a MEK1/2 inhibitor, mainly decreased cell division rates. OSI-906 inhibited the division rate as effective as paclitaxel within the range of tested doses. LCL161 is a small molecule SMAC mimetic that antagonizes multiple inhibitor of apoptosis (IAP) family proteins and augments apoptosis induction. Consistent with its mechanism of action, LCL161 showed a minimal effect on cell division while strongly enhancing cell death. Taken together, the growth-death model allows us to interpret the response of cancer cells to different targeted drugs in terms of cell division and death rates, which otherwise would not be revealed from overall phenotypic changes without a panel of experiments.
Compounds with disparate phenotypic outcomes can appear synergistic when only analyzed by viability Based on the changes in rate parameters by targeted drugs in Fig. 3, we wondered how drugs with divergent phenotypic effects might influence cell behavior when combined. We selected one compound from each of the three groups identified from Fig. 3; PIM447 affected both cell division and death rates, OSI-906 affected only cell division, and LCL161 affected only cell death. Combination treatment between OSI-906 and LCL161 or PIM447 was quantified. Each compound's effect closely matched those we expected from the single-agent treatments (Figs. 3 and 4a, f; Supplementary Fig. S4a, b).
Intriguingly, the nature of each drug interaction was dependent upon whether we took cell death into account. A Bliss additivity model using just cell confluence indicated that both combinations led to a synergistic interaction (Fig. 4b, g). However, we also used our combination treatment data to fit a model in which both growth and death rates respond based on a Hill curve dose-response relationship ( Supplementary Fig. S4a, b). A model assuming Bliss additivity for the growth rate and an additive interaction for the death rate fit our measurements much more closely, despite also having to account for both phase and cell death measurements (Fig. 4c, h). Investigating the source of this discrepancy, we noted that the perceived synergy arose with an increase in cell death ( Fig. 4b/g vs. d/i). We were surprised by the difference in outcomes between each model as the IC 50 and Hill coefficient of the latter model (Fig. 4c, h) are assumed to be equivalent for both phenotypes. However, we expect that the difference arises from the very distinct E max values for growth and death phenotypes (Fig. 4e, j). These terms do not simplify into one maximal effect term when entered into an exponential growth model accounting for both growth and death (methods), indicating a complex interaction between each phenotype's effects. Indeed, our model closely matched the results of just analyzing viability for co-treatment with binimetinib and OSI-906, which both preferentially inhibit growth (Fig. 4k-o). Therefore, we conclude that not only are measures of cell death an important component of pharmacologic response but that ignoring cell death can lead to spurious conclusions of drug synergy.

Discussion
More than half of preclinical agents entering clinical trials fail to be approved due to lack of efficacy 20,21 , and the success rate for clinical approval of oncology drugs remains low as 13% 22 . These statistics reveal the limitation in current preclinical experimental models and analysis in predicting clinical outcomes. Recently, new strategies of evaluating drug response, such as patientderived xenografts and organoids, have been introduced to more closely match patient's tumors to effective compounds [23][24][25] . Two major phenotypic changes we can observe as cellular responses to drugs are cell growth and death. Each response differs according to the administered drug and contributes to the overall effect. In this paper, we demonstrate that cell viability measurements alone cannot distinguish these phenotypic outcomes (Fig. 1). These responses can strongly influence clinical outcome 11 . Thus, we propose a high-throughput assay that can be paired with analysis to quantify both the growth and death effects of drug response for better estimation of drug efficacy. Applying this analysis, we identify that compounds can have distinct outcomes by specifically driving growth, death, or mixed effects (Fig. 3). Cell death can also take distinct forms, which can easily be quantified as relatively more or less apoptotic (Fig. 2). Furthermore, additive drug interactions in terms of cell growth or death can appear synergistic when only assessed through cell viability (Fig. 4). Overall, these results show that cytotoxic drug response should be assessed by the distinct phenotypes of cell growth and death, and then demonstrate an approach to do so.
Separating these phenotypic outcomes provides future opportunities for cancer treatment. A more detailed view of phenotypic drug response should enable treatment optimization both in a population-and patient-specific manner 26 . For maximally effective cancer treatment, single or combination treatments should likely modulate both phenotypes. Purely cytostatic therapy leaves drugpersister cells able to undergo genetic or epigenetic changes giving rise to resistance 9 . On the other hand, cell death may not overcome the replenishment of tumor cells at toxicity-limiting doses. Among the targeted therapies we evaluated, an EGFR tyrosine kinase inhibitor erlotinib exhibited strong inhibition on cell division while modest induction on cell death. Consistent with our result, recent studies on mechanism of erlotinib resistance reported that drug-tolerant colonies maintain throughout drug exposure in quiescent state and expand over long period of time acquiring resistance 8,10 . A combination of erlotinib with a cell death-inducing agent may therefore prevent the survival of drug-tolerant subpopulation cells.
Our analysis results in drug combination also provided a plausible cause for frequent failure of preclinical regimens in clinical trials. The insulin-like growth factor-1 receptor (IGF-1R) signaling pathway is a wellcharacterized pathway involved in cancer cell survival and promoting drug resistance 27 . IGF-1R inhibitors have been validated for their efficacy in various cancer types as mono-and combined therapies in preclinical studies [28][29][30] . Despite promising preclinical evidence, only limited responses were observed in clinical trials [31][32][33][34] . From our drug combination experiments with the IGF-1R inhibitor OSI-906, we discovered that the analysis by cell viability alone found a strong synergistic interaction, while our analysis accounting cell growth and death identified an additive response. Thus, synergy observed in preclinical assays may have been due to the nature of the assay moreso than molecular synergy.
This more detailed analysis of phenotypic response can also help to optimize the microenvironment and host immune response to cancer. For example, host response is potently influenced by the route of cell death 35,36 . These methods can be applied not only in vitro but also in more complex, translationally relevant models, such as organoids and in vivo 37 . By doing so, treatments could be tailored to not just maximally reduce cell viability but optimize cytotoxicity and cell death programs to mount a host response 38 .
The approach here will benefit from improvements in the single-cell resolution of drug response and more exactly distinguishing cell death programs. Cell-to-cell variability impacts both drug response and the development of resistance 39 , and single-cell technologies have enabled in-depth molecular analyses of this heterogeneity 40,41 . Within a tumor population, drug treatment can dynamically shift the balance of variability present, and so dynamic information will likely be critical 8,42 . Automated cell tracking would enable drug-response quantification, with single-cell resolution, while preserving the lineage relationship of cells 43,44 . Finally, while we distinguish relatively more or less apoptotic cell death here, various cell death programs exist 45 . Improved methods to visualize distinct forms of cell death in populations of cells will allow distinct forms of cell death to be separately quantified 46 . The diverse outcomes observed in our study suggest that deeper phenotypic characterization will uncover more differences between compounds in their elicited drug responses.

End-point cell viability assay and time-lapse microscopy
For the end-point cell viability assay in Fig. 1, cells were seeded at 1.5 × 10 3 cells per well in 96-well plates and cultured overnight. Then, cells were treated with doxorubicin. After 72 h, CellTiter Glo reagent (Promega, Madison, WI) was added to each well, and luminescence was detected according to the manufacturer's instructions.
Cells were seeded as indicated above for the time-lapse measurements. The next day, each indicated treatment was added, along with IncuCyte Annexin V Green Reagent (Essen BioScience, Ann Arbor, MI) and 300 nM YOYO-3 in media containing 1.25 mM CaCl 2 . Cells were then cultured and imaged within the IncuCyte Zoom or S3 (Essen BioScience) every 3 h. Four fields of view were taken per well. Fluorescence images were manually thresholded, and the fraction of image area with Annexin V and/or YOYO-3 signal was quantified using IncuCyte Zoom or S3 software (Essen BioScience). Finally, the fraction of area occupied by cells was analyzed by brightfield analysis.
Hill curve identifiability model related to Fig. 1 A model of exponential growth along with death was fit to viability measurements assuming a Hill dose-response relationship. For comparing the model to the data, the fit residuals were assumed to be normally distributed. In the absence of drug, the growth rate was measured and experimentally set to be 0.0315 1/h, and cells were assumed to not undergo cell death. The minimum growth rate (at infinite concentration of drug) was fit using a uniform prior between 0.0 and the growth rate in the absence of drug. The maximal death rate (at infinite concentration of drug) was fit using a log-normal prior of −2.0 ± 2.0 1/h (log 10 scale). The Hill slope was fit using a log-normal prior of 0.0 ± 1.0 (log 10 scale). Both the IC 50 and Hill slope were assumed to be the same for growth and death rates.

Growth model structure
Cell behavior was modeled using a series of kinetic equations incorporating cell growth and death. We represent the overall state of a population of cells as v ¼ ½L; E; D a ; D n , respectively indicating the number of live cells, cells within early apoptosis, dead cells via apoptosis, and dead cells via a non-apoptotic process. Using such notation, the time derivative was defined as: where R g (or div) is the rate of cell division, R d (or deathRate) is the rate of cell death, f (or apopFrac) is the fraction of dying cells which go through apoptosis, and τ (or d) determines the rate of conversion from early to late apoptosis.
integrating these equations provides the solution:

Growth model inference
Predicted cell numbers were fit to experimental measurements using Markov chain Monte Carlo 47 . The percent area positive for cell confluence, Annexin V stain, or YOYO-3 stain was quantified and assumed to be proportional to the number of cells positive for each marker. Cell confluence was assumed to be the total of cells in all states. Apoptotic cells were assumed to be positive for Annexin V signal, then positive for both signals after late apoptosis. Non-apoptotic cells were assumed to just be positive for YOYO-3 signal after dying. Each rate parameter was fit to the corresponding measurements within a single drug condition over time. An entire experiment, corresponding to a set of different compounds and concentrations, was fit simultaneously, allowing for a background offset and conversion factor of each quantity to be fit across the experiment. div was set to have a uniform prior of 0.0-0.35 1/h. deathRate, and d were set to have log-normal prior distributions of mean 0.01 1/h with standard deviation 0.5 (log 10 scale). By inspecting a calibration experiment and manually counting the cells within a field, we measured the conversion between number of cells and area of signal for the confluence, Annexin V, and YOYO-3 images. In addition, we quantified the ratio of positive area for each pair of signals when a single cell was positive for both. Each of these were set as log-normal prior distributions on the conversion values between number of cells and positive area. Finally, we observed appreciable background in the Annexin V and YOYO-3 signal, leading to signal in the absence of cells. Therefore, we set log-normal priors for the background levels with mean 0.1% of area and standard deviation of 0.1 (log 10 scale). Each data point was assumed to have independent, normally distributed error around the model prediction.
Sampling convergence was verified by checking that two independent runs generated insignificant differences, checking for ergodicity through the Geweke criterion comparing the first and second half of each run, and verifying an effective sample size of greater than 200. Sampling failures were solved by increasing the number of tuning samples.

CFSE-based cell proliferation analysis
Cell division was measured using carboxyfluorescein diacetate succinimidyl ester (CFSE) dilution analysis. Cells were labeled with 5 μM CFSE (Invitrogen, Carlsbad, CA) according to the manufacturer's protocol. The stained cells were seeded overnight in 60-mm dishes at a density of 2 × 10 5 cells per dish, and then treated with indicated drugs next day. For 72 h at 24-h intervals, cells were collected and fixed in 4% paraformaldehyde prior to acquisition on a BD LSRFortessa flow cytometer (BD Biosciences, San Jose, CA). CFSE signal intensity of 1 × 10 4 cells was recorded and analyzed to measure cell divisions. The same cell line was labeled the day of the analysis to determine initial labeling.

Western blot analysis
Cells were seeded at a density of 2 × 10 5 cells per 60-mm dish 24 h prior to drug treatment then treated with the indicated conditions for 24, 48, and 72 h. After incubation, cells were lysed in 10 mM Tris-HCl pH 8.0, 1 mM EDTA, 1% Triton-X 100, 0.1% Na deoxycholate, 0.1% SDS, and 140 mM NaCl, freshly supplemented with protease and phosphatase inhibitor (Boston Bio Products, Ashland, MA). Protein concentration was measured by a bicinchoninic acid assay. In total, 10 μg of protein from each cell lysate was subjected to SDS-PAGE, and then transferred to a polyvinylidene difluoride membrane. Each membrane was incubated overnight with antibody against cleaved caspase-3 (Cell Signaling Technology, Danvers, MA, #9664) or 1.5 h with HRP-conjugated β-actin antibody (Cell Signaling Technology, #12262). β-actin was used as a loading control for western blot analysis.

Drug interaction fitting
Drug interaction was assumed to follow the Bliss independence model 48 . Where indicated, this was taken to be defined as a proportional decrease in the viability of cells. That is, cell viability was normalized to 1.0 for the control condition, and then the proportional decrease in cell viability was calculated by 1.0 minus cell viability. Synergy or antagonism was identified by a greater or lesser decrease in viability than predicted, respectively.
Alternatively, Bliss additivity was defined in conjunction with a model incorporating cell death. d and apopfrac were assumed to be constant across drug concentration or combination and fit using the same priors as before. The growth rate in the absence of drug was fit using the lognormal prior of −1.5 ± 0.1/h (log 10 scale) based on experimental growth measurement. Cells were assumed to undergo no cell death in the absence of drug. An E max of growth inhibition was fit using a Beta prior (ɑ = 1.0, β = 1.0), where 1.0 indicates complete growth inhibition and 0.0 no growth inhibition. The E max of death effect was fit using a log-normal prior of −2.0 ± 0.5/h (log 10 scale), where the value indicates the maximal death rate. The half-maximal concentration (EC 50 or IC 50 ) and Hill coefficient of each compound were fit using the same priors as before for these quantities and assumed to be the same for both growth and death effects.