Research ArticleCancer therapy

Modeling chemotherapy-induced stress to identify rational combination therapies in the DNA damage response pathway

See allHide authors and affiliations

Science Signaling  24 Jul 2018:
Vol. 11, Issue 540, eaat0229
DOI: 10.1126/scisignal.aat0229

Modeling drug-chemotherapy combinations

Pharmacological inhibitors of DNA repair pathways can enhance the efficacy of chemotherapy. Finding those combinations that are most effective in a certain cancer type or exploit vulnerabilities therein is key to improving clinical outcomes. Alkan et al. incorporated various molecular, biochemical, and cellular parameters into a model of cancer cell responses to various chemotherapies to predict and then test the efficacy of combined DNA repair inhibition. Among the expected synergy of chemotherapy with an inhibitor of the kinase ATR, both in culture and in vivo using delayed-release liposome technology, there were some unexpected differences in outcome when applying low-dose versus high-dose chemotherapy regimens. These findings show how modeling in cells can both help inform therapeutic development and reveal new biology to investigate further.


Cells respond to DNA damage by activating complex signaling networks that decide cell fate, promoting not only DNA damage repair and survival but also cell death. We have developed a multiscale computational model that quantitatively links chemotherapy-induced DNA damage response signaling to cell fate. The computational model was trained and calibrated on extensive data from U2OS osteosarcoma cells, including the cell cycle distribution of the initial cell population, signaling data measured by Western blotting, and cell fate data in response to chemotherapy treatment measured by time-lapse microscopy. The resulting mechanistic model predicted the cellular responses to chemotherapy alone and in combination with targeted inhibitors of the DNA damage response pathway, which we confirmed experimentally. Computational models such as the one presented here can be used to understand the molecular basis that defines the complex interplay between cell survival and cell death and to rationally identify chemotherapy-potentiating drug combinations.


Chemotherapy drugs commonly target rapidly dividing cells, preferentially cancer cells, by inducing replicative stress (1). Upon genotoxic stress, cells trigger complex signaling pathways that connect cell cycle arrest and DNA repair in G1, S, or G2/M phase with phenotypic cell fate decisions made between survival, cell cycle re-entry and proliferation, senescence, and cell death (2). The complexity arises from the observation that the DNA damage response simultaneously activates prodeath mechanisms, such as apoptosis and necroptosis, and prosurvival mechanisms including DNA repair and autophagy, often via the same signaling proteins and pathways. Within this context, thresholds appear to be important, but our mechanistic understanding of how cells determine cell fate upon genotoxic stress remains limited (3, 4).

The DNA damage signaling pathway has been extensively reviewed in the literature (58). Although complex, many of the signaling proteins that compose this pathway are well described, and their protein-protein interactions are largely known. A major regulatory function is assigned to the atypical protein kinases (PKs) of the phosphatidylinositol 3-kinase (PI3K)–related kinase family members: DNA-dependent PK (DNA-PK), ataxia telangiectasia–mutated (ATM) PK, and ataxia telangiectasia and Rad3-related (ATR) PK (9). They are associated with the recruitment of repair proteins to the site of damage and are implicated in various repair pathways, such as homologous recombination (HR) and nonhomologous end-joining (NHEJ) for double-strand break (DSB) repair, base excision repair (BER) or nucleotide excision repair (NER) for single-strand break (SSB) repair, and replication fork stabilization (10). In addition to their role in various repair processes, the signals from these master kinases are transmitted through additional kinases and phosphatases, such as cell cycle checkpoint (CHK1/2) and various cell division cycle proteins. Eventually, cyclin-dependent kinases and/or p53 transcriptionally induced proteins, such as p21, regulate the cell cycle. Arresting the cell cycle is necessary for allowing successful completion of DNA repair (11). In the case of sustained DNA damage, cells can undergo programmed cell death by inducing apoptotic signaling (12).

Synthetic lethal approaches to cancer therapy development have provided novel strategies to specifically target cancer cells while sparing noncancer cells and thereby reducing toxicity. The U.S. Food and Drug Administration approvals of three poly(adenosine 5′-diphosphate–ribose) polymerase inhibitors—olaparib (Lynparza, AstraZeneca), rucaparib (Rubraca, Clovis Oncology), and niraparib (Zejula, Tesaro Inc.)—in BRCA mutant high-grade ovarian cancer prove that the DNA damage repair pathway can successfully be exploited to treat certain cancer types. Inhibitors targeting the DNA damage response pathway exhibit activity in combination with chemotherapies but might also have benefits as monotherapy in cancers that exhibit replication stress (RS) and thus may decrease normal tissue toxicity (13, 14). Targeted inhibitors against almost any kinase within the DNA damage repair pathway (ATR, CHK1, ATM, CHK2, and DNA-PK) have been created and are currently in clinical development with mixed success (15). Currently, there are two ATR inhibitors, VX-970 (recently acquired by Merck KGaA from Vertex Pharmaceuticals) and AZD6738 (AstraZeneca), in multiple phase 1 clinical trials for the treatment of solid tumors in combination with DNA-damaging agents (16). In addition, there are three CHK1 inhibitors: SRA737 (previously known as CCT245737, Sierra Oncology Inc.), MK-8776 (previously known as SCH-900776, Merck and Co.), and prexasertib (LY2606368, Lilly Oncology) in early clinical development. The field is expanding rapidly, but the challenge moving forward will be how to identify patients who will benefit from combination regimens containing a specific DNA damage response pathway inhibitor and DNA-damaging agent (chemotherapy) regimen and, thus, to maximize the clinical benefit to patients.

Computational models provide a tool to quantitatively understand the complex interplay between survival and cell death. Tentner et al. (17) developed a quantitative signal response data set for doxorubicin-induced DNA damage and used relational modeling techniques to understand the interplay between signaling pathways and phenotypic response. They found that the cell fate choice between apoptosis and cell cycle arrest in the case of doxorubicin is dose-dependent and modulated by tumor necrosis factor–α. They also demonstrated that extracellular signal–regulated kinases 1 and 2 are required for the maintenance of a G1/S arrest and subsequent apoptotic cell death induced after exposure to low doses of doxorubicin. Other computational modeling efforts concentrated on the pharmacokinetic and pharmacodynamic (PK/PD) effects of targeted inhibitors of the DNA damage response to inform the translation of preclinical findings into the clinic (18, 19).

Here, we applied a mechanistic computational modeling approach to quantitatively characterize the cell-intrinsic regulation of cell cycle and DNA damage response signaling with respect to phenotypic cell fate decisions. Our model not only predicted therapeutic targets within the DNA damage response that can potentiate the cytotoxicity of several common chemotherapies but also identified combinations that unexpectedly led to attenuation of the chemotherapy-mediated effects. Both cases were experimentally validated.


Construction of a multiscale computational model to predict cell fate in response to chemotherapy

Using ordinary differential equations (ODEs), we developed a quantitative computational model that describes protein phosphorylation events, changes in gene expression, and cell fate in response to stress induced by chemotherapies on different time scales (Fig. 1). DNA damage response signaling occurs within minutes to hours; gene expression changes occur within hours to days; and cellular responses such as DNA damage repair, cell cycle arrest, and apoptosis may happen in the order of days. The model is composed of a cell cycle and a canonical DNA damage response signaling pathway module. As opposed to models that resolve the detailed molecular events of the cell cycle (20), we included a simplified model that describes the cell cycle phases and their relationship to DNA damage. Our cell cycle model distinguishes G1, S, and G2 phases, as well as different types of DNA damage: We included SSBs and DSBs. Depending on the type of DNA damage induced by chemotherapy, cells transition to the respective damaged state and trigger a signaling response characteristic for the type of damage. From the signaling response, processes such as repair, cell cycle arrest, and apoptosis feed back into the cell cycle model. Our model is trained against population-level data and describes signal transduction, repair, arrest, and apoptosis of cells in three different cell cycle stages. This is accomplished by explicitly modeling these three cell populations, including differential action of the chemotherapy treatment on the various cell cycle stages and subsequent divergent signaling events. Thus, the model describes the time evolution of a population of cells in different cell cycle states. Similarly, the signaling module of the computational model independently describes the average signaling response for cells in different cell cycle states.

Fig. 1 Illustration of the computational model describing the DNA damage response signaling, gene expression changes, cell cycle stages, and cellular responses.

The model includes a representation of the cell cycle (stages G1, S, and G2) and the DNA damage repair signaling pathways (gray boxes). Chemotherapy dose and the cell cycle composition before treatment are the model inputs. Doxorubicin leads to DSBs in G2 phase, whereas SN38 and gemcitabine induce SSBs in S phase. DSBs and SSBs trigger distinct individual branches of the signaling model. Stage-dependent cell cycle checkpoints (stop signs) and different DNA repair pathways (NHEJ, HR, and NER/BER) link the signaling model back to the cell cycle model. The model includes p53-dependent apoptosis and drug effects on proliferation and survival signaling. Node-labeled “m_” indicates mRNA. The reactions are labeled with ν1 to ν102; the complete mathematical description of the computational model is available in data file S1. DDR, DNA damage response. qRT-PCR, quantitative reverse transcription polymerase chain reaction.

Three DNA-damaging agents were the focus of this study: gemcitabine, irinotecan in its active form SN38, and doxorubicin. Gemcitabine is a nucleotide analog that is incorporated into DNA during S phase. It blocks DNA polymerases, stalls replication forks, and leads to RS and SSBs (21, 22). SN38 stabilizes the complex between topoisomerase-I and DNA and leads to SSBs mostly in S phase (2325). Upon collision with moving DNA replication forks, those SSBs can become DSBs by collapse of the replication forks. Doxorubicin has a diverse set of mechanisms of action (26); we were particularly interested in its ability to stabilize topoisomerase-II and induce DSBs (27). By inclusion of a separate reaction pathway that only operates at high concentrations of doxorubicin, we also accounted for its additional non-DNA damage response–related mechanisms of action (28). In a simplified way, this additional reaction can represent additional cytotoxicity by the generation of reactive oxygen species (29) or the inhibition of transcription (30). Those mechanisms lead to additional cytotoxicity at doses higher than about 0.1 μM in vitro. As a result, our model harbors two pathways that can lead to apoptosis: One pathway is triggered by DNA damage signaling–related mechanisms, and the other pathway is triggered by non-DNA damage response–related mechanisms. The two types of DNA damage, SSBs and DSBs, trigger different parts of the signaling network by phosphorylating DNA-PK, ATM, and ATR. We used the chemotherapy dose as input into our computational model to simulate the activation of the DNA damage response pathways, cell cycle distribution, and ultimately cell fate (Fig. 1).

Different components of the signaling network are linked to the cellular responses and are triggered proportionally. For instance, the higher the phosphorylation level of CHK1, the lower the transition rates from S to G2 phase. Feedback loops to the cell cycle model include both the induction of cell cycle checkpoints and the activation of DNA damage response (Fig. 1). We implemented an early S-phase checkpoint that is triggered by transcriptional induction of p21, an S-phase checkpoint that is triggered by phospho-CHK1–dependent signaling, and a G2 checkpoint that is triggered by phospho-CHK2 (pCHK2)–dependent signaling. The repair processes shift the distribution of cells from the damaged, arrested cell cycle states back to the corresponding cycling states. Activity of phospho–DNA-PK (pDNA-PK) increases NHEJ repair, the activity of phospho-ATM (pATM) and phospho-ATR increases HR, and the activity of phospho-ATR increases SSB repair or replication fork stabilization.

To calibrate the model, an experimental data set using various experimental methods has been generated, which captures the time course and dose dependency of the response to chemotherapy at the mRNA, protein, phosphoprotein, and cellular levels. The cellular level consists of tracking proliferation and apoptosis by live-cell imaging. The data set was tailored to the scope and complexity of the computational model, which allowed us to identify the 95% confidence intervals for half of the parameters and lower or upper bounds for the remaining parameters. As a model system, we chose the U2OS osteosarcoma cell line because a large body of data already exists in the literature. For model calibration, more than 3000 individual data points have been measured for different time points and different chemotherapy treatment conditions. A variety of the experimental methods and data was used to calibrate the computational model (Fig. 2). Time-lapse microscopy data captured the dose-dependent effects of chemotherapies on the dynamics of cell proliferation and death (Fig. 2A). A red fluorescent reporter was used to count the number of living cells over time, and a green live cell dye reporting Caspase 3/7 activity was used to enumerate apoptotic cells. Flow cytometry was used to determine the cell cycle distribution before drug treatment as input into the computational model (Fig. 2B). The signaling dynamics of eight proteins of the DNA damage response pathway (pATM, pDNA-PK, pCHK1, pCHK2, phospho- and total p53, total p21, and yH2AX) were measured over a 24-hour time course at different drug concentrations by quantitative Western blotting (Fig. 2C).

Fig. 2 Experimental data sets used to train the computational model.

(A) Proliferation and cell death: exemplary microscopy images of untreated U2OS cells at the 24-hour time point and of cells treated with 1 μM doxorubicin at the 48-hour time point (additional images shown in fig. S1). (B) Distribution of cell cycle stages at baseline quantified by flow cytometry measuring DNA content [DAPI (4′,6-diamidino-2-phenylindole)] and DNA synthesis rate [EdU (5-ethynyl-2′-deoxyuridine)]. (C) DNA damage response signaling in U2OS cells was exposed to 0, 0.01, 0.05, and 0.1 μM gemcitabine, SN38, and doxorubicin; the abundance of phosphoproteins was assessed over the subsequent 24 hours. pATM, pCHK1, p53, and p21 were measured by Western blot; representative examples are shown. GAPDH, glyceraldehyde-3-phosphate dehydrogenase.

The experimental data used to calibrate the computational model consisted of time course data of all phosphorylated signaling proteins as a function of three concentrations, the full dose response at the 6-hour time point as measured by quantitative Western blotting, and the time-lapse data enumerating viable versus apoptotic cells over time at five doses (Fig. 3). The computational model was calibrated against the complete experimental training data set simultaneously. The model calibration involves tuning of the model parameters (for example, reaction rate constants or transition rates) to find the best possible fit of the model outputs (for example, dynamics of population growth or phospho-signaling) to the experimentally measured counterparts. Common challenges, especially when considering larger-scale models, are that the numerical algorithms that perform the model fitting can get trapped in local optima and that the identified parameter values retain high uncertainty (that is, nonidentifiability problem). We have analyzed the optimality of the model parameters and how well they are constrained by the experimental data (identifiability) and have found that the optimal set of parameters that we identified satisfies quality measures (figs. S2 and S3). The simulation results for the optimal set of parameters were compared to the experimental data (Fig. 3). Additional model simulations and data, including the two transcriptional targets of the DNA damage response pathway (p21 and Wip1 mRNA) measured by qRT-PCR, are shown in data file S1. Overall, the model captures the data well considering the scope and the relative simplicity of the model. We also quantified the degree of misfit for individual experimental conditions to identify parts of the model that could be further improved (table S1). Small-molecule inhibitors against the master kinases DNA-PK, ATM, and ATR often have off-target effects on other family members such as mammalian target of rapamycin or on the structurally similar PI3K kinases (31, 32). To represent those possible polypharmacological effects, we implemented on-target effects of the small-molecule inhibitors on the DNA damage repair signaling pathway and off-target effects on proliferation and survival pathways. Off-target effect can be statistically significant at higher concentrations and may involve many possible targets. We chose to implement a broad effect of the inhibitors in our model by directly inhibiting the initiation of S phase and inducing apoptosis of all cell populations independent of cell cycle stage.

Fig. 3 Trained computational model describes signaling and proliferation data as a function of time and dose.

(A and B) Quantification of DNA damage response protein phosphorylation by Western blot data (represented in Fig. 2C), analyzed as response to 0.01, 0.05, and 0.1 μM chemotherapy over 24 hours (A), and dose-response curves at the 6-hour time point (B). (C) Imaging-based quantification of cell proliferation dynamics with and without chemotherapy exposure by counting nuclear reporter–positive cells. Data are normalized to time point zero. Dots represent data of separate experiments; each replicate was individually used to calibrate the computational model. Lines indicate the corresponding simulation trajectories. Further experimental details are described in Methods.

Mathematically, the model consists of a set of ODEs. Each equation represents the dose-time behavior of an entity in our model (nodes in Fig. 1), for instance, protein levels in phosphorylated and nonphosphorylated form, mRNA levels, or amount of viable, damaged or apoptotic cells in the different cell cycle stages. Transition reactions between equations (coupling) are used to model transitions between the different cell populations (for example, moving viable cells to apoptotic cells) and protein species (for example, moving nonphosphorylated protein to phosphorylated protein promoted by the level of a kinase).

Identification of chemotherapy-potentiating drug combinations

It is well established that agents targeting the DNA damage response pathway can potentiate chemotherapies. We applied our computational model to perform a sensitivity analysis and identify which targets lead to the maximal inhibition of proliferation or cell death in the presence of doxorubicin, gemcitabine, or SN38. To this end, we simulated the effect on proliferation by combining each chemotherapy with any possible target within the DNA damage response pathway in silico (Fig. 4). Our model can quantitatively distinguish between proliferation, stasis, and cell death. We selected two concentrations for the chemotherapies: a subcytostatic dose of 0.001 μM (Fig. 4A) and a cytostatic or, in the case of doxorubicin, even cytotoxic dose of 1 μM (Fig. 4B). The effect of modulating the DNA damage response pathway was considered by simulating the inhibition of eight molecular targets within the DNA damage response pathway (ATM, ATR, DNA-PK, CHK1, CHK2, p53, p21, and a combination of CHK1 and CHK2). We assumed hypothetical potent targeted inhibitors with a subnanomolar dissociation constant (KD of 0.5 nM at a dose of 1 μM). We evaluate the fold change in cell number by comparing the effect of the drug combination after 3 days of continuous drug exposure versus chemotherapy alone.

Fig. 4 Sensitivity analysis of cellular responses to perturbation of the DNA damage response signaling in the presence of chemotherapies.

(A and B) For the chemotherapies, a low dose of 0.001 μM (A) and a high dose of 1 μM (B) were simulated. The inhibition of eight molecular targets within the DNA damage response pathway (ATM, ATR, DNA-PK, CHK1, CHK2, p53, p21, and CHK1/2) was simulated. The effect of chemotherapies is simulated as monotherapy or in combination (bars) and shown on the x axis. The outcome is quantified as log10 fold change of cell number after 3 days of the drug exposure. The vertical dashed line is aligned to the effect of the chemotherapy alone (first bar in each and dashed lines), which is the comparator for the effect of the combinations.

Upon consideration of inhibitor effects at high- and low-dose monotherapy and combination chemotherapy, the strongest predicted potentiation effect is observed for the combination of an ATR or CHK1 inhibitor with either gemcitabine or SN38 (or with both chemotherapies combined). At low and high chemotherapy concentrations, these combinations potentiate stasis and cell death, respectively (Fig. 4B). The model can explain this effect by the strong induction of SSBs and activation of the ATR/CHK1 signaling pathway by both gemcitabine and SN38. In this damaged state, the model predicts that cells rely on cell cycle arrest via the CHK1 pathway. Inhibition of ATR or CHK1 removes this checkpoint while cells are still damaged and thus drives the cells into cell death. With doxorubicin, the combination effect of ATR/CHK1 inhibition is attenuated because of the induction of DSBs in addition to SSBs, in which case cells are arrested predominantly by the ATM/CHK2 checkpoint and the inhibition of the CHK1 pathway is less effective. The model also predicts that, for low concentrations of doxorubicin, a DNA-PK inhibitor would increase stasis (Fig. 4A). At high concentrations, doxorubicin alone dominates the response, and no potentiation effect is predicted (Fig. 4B). The addition of an ATM inhibitor to high chemotherapy concentrations does not improve the effect but appears to attenuate the effect at low-dose chemotherapy.

Dose-response profiles of potentiating drug combinations and experimental validation

To further investigate the identified combinations, we performed simulations over a range of concentrations for an ATR or CHK1 inhibitor with gemcitabine and DNA-PK and an ATM inhibitor with doxorubicin to obtain drug response profiles. We compared model simulations to the experimental validation for combinations of the small-molecule inhibitors—VE-822 (ATR inhibitor), MK-8776 (CHK1 inhibitor), KU-0060648 (DNA-PK inhibitor), and KU-60019 (ATM inhibitor)—and chemotherapies and varied the drugs over a wide concentration range (Fig. 5 and figs. S4 and S5). The validation data set is an independent data set that was not used for the model building. Overall, the model predictions are in good qualitative agreement with the experimental data. For the ATR or CHK1 combination with gemcitabine, potentiated cytotoxicity is predicted by the model over a wide concentration range (Fig. 5A) and validated by the experimental data (Fig. 5B). For the ATR inhibitor, there is a difference between model simulation and validation data at high doses (Fig. 5B, white dashed boxes). At concentrations above 1 μM, the ATR inhibitor VE-822 is cytotoxic as single agent. This is potentially due to off-target effects such as multikinase inhibition. For the model simulations, we assume perfect target specificity of the in silico inhibitors. This explains the observed differences. The ATR inhibitor VE-822 and the CHK1 inhibitor MK-8776 differ in the quantity of their off-target effects, with the former leading to stronger cytotoxicity. The predicted combination of doxorubicin with the DNA-PK inhibitor shows an increase of cytostasis in the lower concentration range of chemotherapy. Notably, only potentiation of a cytostatic effect is achieved by this combination. It does not lead to potentiation of cell death as the ATR or CHK1 combinations. For the highest concentrations (Fig. 5B), this combination also shows increased cytotoxic effects that the model simulation fails to capture. At concentrations above 0.1 μM, doxorubicin is known to induce a non–topoisomerase II–mediated form of cell death in addition to its DNA damage–related effect (29). Our computational model is focused on the DNA damage response pathway and does not capture this aspect of doxorubicin mechanistically. Only at the highest concentration can an increase in cell death be achieved with this combination, which matches the experimental data. As single drugs, both the DNA-PK inhibitor KU-0060648 and the ATM inhibitor KU-60019 used here did not lead to statistically significant cytotoxic effects even at doses above 1 μM.

Fig. 5 Drug response profiles and experimental validation.

(A and B) Drug response profiles for gemcitabine and doxorubicin predicted by the computational model (A) and experimental validation using small-molecule inhibitors of ATR, CHK1, DNA-PK, and ATM (B). For the analysis, the drug and inhibitor concentrations were varied over a wide range, and outcome is quantified after 3 days. The black line indicates the stasis line as predicted by the computational model. Blue stars represent the low-dose chemotherapy, and the red stars indicate the high-dose chemotherapy in combination with a targeted inhibitor as chosen for the single-dose combinations (Fig. 4). The dashed white boxes indicate areas where the model prediction and experimental data are discordant.

Mechanism of action inferred from proliferation dynamics of drug combinations

Next to the dose, the choice of the time point is essential to observe the maximal effect of a drug combination. To further validate the model, we compared simulated proliferation/cell death dynamics for selected combinations and concentrations to experimental validation (Fig. 6). The experimental data sets containing inhibitors were not used during model building. Overall, the model and experimental data are in good agreement. We treated U2OS cells with 0.04 μM gemcitabine and 1 μM VE-822 and observed that, whereas gemcitabine alone only leads to cytostatic effects over the 3-day time course, the combination with an ATR inhibitor shows a decrease in cell number after 24 hours of treatment, which is in concordance with the observed increase in apoptotic cells during the same period (Fig. 6A). Within the first 24 hours of drug exposure, there is a small but significant increase in the number of live cells detected for the combination. The model simulation replicates this behavior and can therefore help to mechanistically explain this observation. Gemcitabine leads to cytostatic effect with a strong induction of the S-phase cell cycle checkpoint by the ATR/CHK1 pathway. Combining gemcitabine with ATR (or CHK1) inhibition releases this S-phase checkpoint while cells are still damaged. This forces cells into progressing their cell cycle prematurely and to succumb to the consequences of the sustained DNA damage. The premature cell division leads to a temporary increase in the experimentally counted number of live (nuclear reporter positive) cells by live-cell imaging, followed by a rapid reduction of the overall cell population. This mechanism of action is captured by time-lapse microscopy movies available in movie S1. Similar observations have been described earlier (33). After 48 hours, the number of apoptotic cells declined. The model was not able to capture this rapid decline. This is due to the limited ability of the current model to produce a sharp, time-delayed induction of apoptosis. To facilitate a delay in apoptosis in the model, we implemented a linear chain of reactions that can approximate a delay and chose a chain length of five. To increase the “sharpness” of the delay, a longer chain would need to be used. This would increase the model size substantially and make subsequent computations more time-consuming. Another possibility involves the use of time-delay differential equations, but they are not often used in the field because of more difficult numerical handling. We feel that the misfit is not statistically significant enough to justify a different model that would result in more problematic handling. Please note also that the number of apoptotic cells is relatively small compared to the number of viable cells at any given time. We further used the combination of 0.0016 μM doxorubicin and 0.5 μM KU-0060648 and observed that, whereas the single agents are not effective at these concentrations, the combination with a DNA-PK inhibitor leads to a significant induction of cytostasis (Fig. 6B). This is achieved by sustaining the DNA damage as initially predicted using the sensitivity analysis (Fig. 4A).

Fig. 6 Dynamics of live and apoptotic cells treated with potentiating drug combinations.

(A and B) Proliferation and apoptosis assessed in U2OS cells cultured with 0.04 μM gemcitabine and 1 μM ATR inhibitor VE-822 (A), or 0.0016 μM doxorubicin and 0.5 μM DNA-PK inhibitor KU-0060648 (B). The experimental data (dots, representative of separate experiments) are compared to the response simulated by the computational model (solid lines) together with their associated uncertainties (shades).

In vitro to in vivo translation of combining ATR inhibitor VE-822 with irinotecan in the MS-751 cervical cancer model

To test whether our in vitro findings in the U2OS cell line have a potentiating effect with an ATR inhibitor and chemotherapy and translates to increased efficacy in vivo, we selected the MS-751 cervical cancer model (the U2OS model does not grow in vivo). For testing a drug combination in vivo, the exposure profile of both drugs is an important consideration. Therefore, we attempted to mimic the in vitro experiment with the in vivo dosing as closely as possible. Both gemcitabine and SN38 had a comparable potentiating effect in combinations with an ATR inhibitor (Figs. 5 and 7A). The reported drug half-life in mice is about 0.28 hours for gemcitabine (34) and 1.1 to 3 hours for irinotecan (prodrug of SN38 that is used in vivo) (35). The potentiating combination that we identified relies on both drugs being present simultaneously at the site of therapeutic effect. Liposome-encapsulated drugs have the advantage of slowly releasing their payloads at the site of the tumor due to a combination of the enhanced permeability and retention effect and engineered controlled release rates, and this may prolong the exposure at the tumor site (36, 37). For our in vivo study, we therefore selected a sustained release formulation of SN38 [nanoliposomal irinotecan (nal-IRI) (38, 39)] as the chemotherapy partner and an experimental liposomal formulation of the ATR inhibitor VE-822 (Ls-VE-822) that we have developed. However, the translation of the effective dose in vitro to an in vivo dosing schedule is challenging. The dosing schedule used here was selected on the basis of the previous experience with the MS-751 model. First, we validated that SN38 in combination with an ATR inhibitor shows a potentiating effect in MS-751 in vitro (Fig. 7B). Despite harboring the HPV E6 protein that has been associated with increase degradation of p53 (40), we observe that MS-751 cells treated with SN38 alone arrest their cell cycle. This requires that either p53 remains functional in these cells or the S-phase checkpoint is not under direct control of p53. When treated with SN38 or gemcitabine, MS-751 cells also show induction of p53-S15 phosphorylation (fig. S6). Consequently, the combination still has the desired potentiating effect. Similarly, the drug combination significantly enhances the antitumor response in animals bearing the MS-751 xenografts in vivo (Fig. 7C). As opposed to the in vitro experiments, where the combinations were applied simultaneously, Ls-VE-822 was administered 24 hours after nal-IRI. The chemotherapy that produces DNA damage, ATR/Chk1 activation, and consequently cell cycle arrest was injected first. The ATR inhibitor was injected second to inhibit ATR/Chk1 and release the cell cycle checkpoint while cells were still damaged. Because of the slow-release formulation, both drugs will still be present in the tumor microenvironment at the same time. Ls-VE-822 alone does not alter tumor growth. Similarly, nal-IRI alone only shows moderate activity at this dose level. The combination data, however, suggest that the cytotoxicity potentiating combination that we identified with our modeling analysis can qualitatively translate into potentiated efficacy in vivo. Because the pharmacokinetics and tissue distribution of drugs, as well as the behavior of the tumor model, differ considerably between in vivo and in vitro, a more accurate translation cannot be expected. We further acknowledge that the differences could be even further amplified in human.

Fig. 7 Response to SN38 and VE-822 (ATR inhibitor) in the U2OS and MS-751 cervical cancer lines in vitro and MS-751 cervical cancer model in vivo.

(A and B) Time course data of cell number measurements in U2OS cells (A) and MS-751 cells (B) after treatment with 0.01 μM SN38, 1 μM VE-822, or in combination (n = 4), overlaid with the corresponding model simulation. Drugs were added at the 24-hour time point. (C) Growth of MS-751 xenografts in mice treated as indicated on days 17, 24, 31, and 38 (dashed lines) after tumor implantation. nal-IRI was given at a dose of 5 mg/kg and liposomal ATR inhibitor Ls-VE-822 at 20 mg/kg intravenously. For the combination treatment, Ls-VE-822 was injected 24 hours after the nal-IRI injection (n = 10 mice per group). PBS, phosphate-buffered saline.


Although we have a detailed understanding of the DNA damage response signaling pathway wiring, it is yet insufficiently understood how the signaling components trigger cell fate decisions between survival, cell cycle re-entry and proliferation, senescence, and cell death. Here, we have presented a computational model that links the DNA damage response signaling to phenotypic cellular responses such as DNA damage repair, cell cycle arrest, and apoptosis and incorporates multiple time scales. Previous work combining computational models and experiments has focused on a detailed understanding of the wiring of the signal transduction pathway (17), DNA damage sensing (41), or PK/PD models (18). To construct the model, we linked a simple computational model of the cell cycle to a canonical signaling model representing the DNA damage response pathway. The model was trained and validated with an extensive signaling data set and time-lapse microscopy tracking the cell fate by enumerating the number of live cells and apoptotic cells. Although we aimed at integrating multiple time scales, we also paid attention to the fact that the overall size of the computational model stays within practical limits. As opposed to developing a model that includes a high degree of mechanistic detail, we decided on a model complexity that is matched to the amount and quality of the available experimental data. Computational models that exhibit greater complexity bear the risks of overfitting, nonidentifiability, and consequently less accurate predictions (42). Our model focuses on the canonical DNA damage signaling pathway that includes targets of interest for cell cycle control. Signaling pathways connected to other cellular responses, such as apoptotic signaling and repair pathways, have been included in a less detailed manner. Because not all cell lines may behave like the U2OS cell line, future work will be required to include the characterization and prediction of heterogeneous responses across multiple cell lines (fig. S7). Understanding under which circumstances a given combination potentiates a chemotherapy will be important for biomarker and indication identification. We recently demonstrated that, for the ErbB pathways, mechanistic models can help to characterize heterogeneous responses in cell lines and have the potential to improve patient stratification (43). The current DNA damage model can serve as a backbone to understand where cellular heterogeneity will most likely influence the response to drug combinations. Model extension to incorporate additional regulatory mechanisms of the DNA damage response pathway might be necessary in the future.

Our integrated approach allows us to study and predict the response of U2OS cells to combination treatments of chemotherapies with inhibitors of the DNA damage response pathways. We applied our model to predict and understand the potentiating combination regimens. Because of our computational approach, the impact of modulators of the DNA damage repair pathway can be studied on multiple time scales: DNA damage signaling—which occurs within minutes—or cell cycle arrest, repair, and cell death—which occur over hours or days. Therefore, conclusions can be made about the dependence of drug combinations on cell cycle stage or the quality of the response, that is, cytostatic versus cytotoxic. Nonintuitive system behavior can be studied most comprehensively with the help of a computational model because of the various nonlinear feedback loops of the DNA damage response system.

Using the computational model, we systematically investigated potentiating drug combinations in vitro and in vivo between DNA damage–inducing chemotherapy and DNA damage signaling modulators. We showed that the efficacy of signaling modulators is highly dependent on the respective chemotherapy combination. We found that the dependence of the effect of the chemotherapy on the cell cycle alone is an important factor in this consideration. We further predicted cytotoxic versus cytostatic effects of various drug combinations and their dose dependency. Experimental validation data in the cell line where we trained the model were in line with model predictions. These results confirm that the ATR/CHK1 pathway in combination with SSB damage–inducing chemotherapy such as SN38 or gemcitabine is a promising drug regimen. Quantitative predictions of the population dynamics of cell growth and death in conjunction with time-lapse imaging data helped to study the mechanism of action of this drug combination. We speculate that the mechanism of action is based on a premature release of damaged cells from the ATR-dependent S-phase cell cycle checkpoint that leads to rapid cell death after mitosis. Finally, the proposed drug combination was tested in vivo using a highly stabilized liposomal formulation of irinotecan and a novel liposomal ATR inhibitor and showed similarly enhanced activity as observed with the in vitro data. The computational model presented here represents a framework that can be used to drive further translational research in the DNA damage response field and may help to design combination regimens that specifically target cancer cells and protect normal cells from chemotherapy treatment but still sensitize cancer cell damage.


Cell culture

U2OS and all other cells were grown in RPMI 1640 supplemented with 10% fetal bovine serum and antibiotics. Imaging and time lapse–compatible NucLight red cell lines were generated as recommended by Essen BioScience Inc. using lentiviral particles, infection, and puromycin selection protocols.

Chemotherapy agents, small-molecule inhibitors, and cell cycle modulators

All drugs, such as chemotherapy agents and small-molecule inhibitors if not mentioned separately, were purchased from Sigma-Aldrich, SelleckChem, Tocris Bioscience, etc. nal-IRI (Onivyde) was a clinical-grade material manufactured and provided by Merrimack Pharmaceuticals Inc. The formulation was described previously (37). For details, see table S2.

Antibodies and Western blots

U2OS cells were exposed to doxorubicin (0, 0.1, 0.5, and 1 μM), SN38 (0, 0.01, 0.1, and 1 μM), and gemcitabine (0, 0.05, 0.1, and 1 μM) over a period of 24 hours. Controls (no drug) and corresponding chemotherapy-treated samples were lysed with 2% SDS containing lysis buffer and stored at −80°C until all samples were collected and run at the same time for the Western blot analysis. The corresponding labeled target specific bands were selected and quantified using the Image Li-Cor Studio software and analysis tool. All target and phospho-specific antibodies were purchased from Cell Signaling Technology and/or Epitomics and were used at 1:1000 dilutions. A list of all antibodies is available in table S3. Standard antibody-based Western blot and immunohistochemistry protocols were used.

Whole-cell lysis protocol

Cells were grown and treated in a 6-cm dish scale. To harvest cells, medium was removed and quickly replaced by ice-cold PBS. PBS was then replaced with 250 μl of 2% SDS lysis buffer. After 5 min of incubation, the lysed cells were scraped off and loaded onto a cell lysate homogenizer microcentrifuge spin-column (catalog no. 79656, QIAshredder, Qiagen). The filtrate was then loaded onto a 0.2-m centrifugal filter column [catalog no. ODM02C34, Nanosep MF (0.2 m), Pall]. The lysates were then stored at −80°C. SDS-containing lysis buffer (2%) was modified after (44).


Live-cell imaging was performed by using the IncuCyte Zoom (Essen BioScience). Fluorescence-activated cell sorter (FACS)–based experiment was performed by using FACSCalibur and/or Aria (BD Biosciences).

Cell cycle stage analysis (FACS)

FACS-based experiment was performed by using Click-iT Plus Assay Kits for Flow Cytometry (Life Technologies) and FACSCalibur and/or Aria (BD Biosciences). Methods and protocols were performed as described by the manufacturer and the Massachusetts Institute of Technology Flow Cytometry Core Facility (

Quantitative gene expression analysis by qRT-PCR

Cells were cultured and treated accordingly in a 96-well plate scale. RNA was extracted using the Qiagen RNeasy Kit (catalog no. 73404, RNeasy Plus Mini Kit, Qiagen) following the manufacturer’s protocol. Nondiluted RNA (~100 ng) was used for complementary DNA synthesis with a commercial kit (catalog no. 4368814, Applied Biosystems) following the manufacturer’s protocol. The gene expression was measured by using and following the manufacturer’s guidelines for the gene-specific TaqMan probes and ViiA Real-Time PCR system and software (Life Technologies).

Computational analysis and model building

Computational analysis and mechanistic model building were performed using the MATLAB software (MathWorks) and the Data2Dynamics modeling environment (45) that is tailored to high-end model calibration and parameter estimation. Mathematical details about the model structure and equations are available in data file S1. The model was calibrated simultaneously to the entirety of the available experimental data. We validated the performance of the parameter calibration procedure by the quality control suggested in (46) by running multiple independent estimations from randomized initial locations and also analyzed identifiability of the model parameters using the profile likelihood method (47). Additional details are available in data file S1.

Preparation of a liposomal ATR inhibitor (VE-822)

Cholesterol (Chol) was purchased from Avanti Polar Lipids, and hydrogenated soy phosphocholine (HSPC) and methoxy-poly(ethylene glycol)-1,2-distearyl-sn-glyceryl (PEG2000-DSG) were obtained from Lipoid GmbH. Chol, HSPC, and PEG2000-DSG were codissolved in 100% ethanol (200 proof, catalog no. 459828, Sigma-Aldrich) at a molar ratio of 3:2:0.03 at 65°C. Sucrose octasulfate (SOS) was acquired from Molekula, and the triethylamine salt was formed by titration of the acidic form of SOS with triethylamine as previously described (37). The solution of premade 1.1 M triethylammonium (TEA)–SOS was mixed with the lipid solution at 60° to 65°C. This suspension was extruded three times through five stacked polycarbonate track–etched filters (Corning Nuclepore) with an average pore size of 100 nm using an argon pressure extruder (Lipex Biomembranes) at 60° to 65°C, and resulting unilamellar liposomes were quickly chilled in ice and then stored at 4°C. The concentration of phospholipid was measured using a phosphate assay, and particle diameter was recorded using a Malvern Nanosizer. Before drug loading, the TEA-SOS gradient was created by removing excess of untrapped TEA-SOS using 500-kDa polysulfone hollow fiber filters (lot no. D06-E500-10-N, Spectrum Labs). Osmolarity of liposomes was balanced by mixing with 50% dextrose solution to a final dextrose concentration of 15% (w/v).

VE-822 (lot no. SSC40215, MedKoo Biosciences) was dissolved in 15% dextrose solutions in deionized water by titration with 1 M HCl and heating at 45°C and then filtrated with 0.2-μm Nalgene 13-mm syringe filters. Drug concentration in the solution was detected by high-performance liquid chromatography. A stock solution of VE-822 containing drugs (8 to 10 mg/ml) was added to the liposomes at a drug/lipid ratio of 800 g/mol phospholipid, and the pH was adjusted to pH 6.5 using 1M Hepes buffer (pH 6.5) and 0.1 N NaOH.

The liposome-drug mixture was incubated with occasional agitation for 30 min at 65°C. The incubation mixture was quickly cooled down and incubated for 10 min at 0°C and then allowed to attain ambient temperature. Unencapsulated drug was removed by 500-kDa polysulfone hollow fiber filters eluted with HBS-6.5 buffer [5 mM Hepes, 250 mM NaCl (pH 6.5)]. Liposome fractions eluted in the void volume were combined, sterilized by 0.2-μm filtration, and stored at 4° to 6°C before use.

Combination therapy on MS-751 cervical cancer xenografts in mice

Antitumor efficacy of nal-IRI in combination with liposomal ATR inhibitor (Ls-VE-822) was studied in the MS-751 cervical xenograft model. The cells were obtained from American Type Culture Collection and propagated in RPMI 1640 medium supplemented with 10% fetal bovine serum, penicillin G (50 U/ml), and streptomycin sulfate (50 μg/ml) at 37°C with 5% CO2 as recommended by the supplier. NCR nu/nu homozygous athymic female nude mice (4 to 5 weeks old; weight, at least 16 g) were obtained from Tacoma. The mice were inoculated subcutaneously in the right flank with 0.1 ml of the suspension containing 5 × 106 cells suspended in PBS supplemented with 30% Matrigel. When tumors achieved the size between 150 and 350 mm3, the animals were assigned to the treatment groups according to the following method. The animals were ranked according to the tumor size and divided into six categories of decreasing tumor size. Four treatment groups of 10 animals per group were formed by randomly selecting one animal from each size category, so that, in each treatment group, all tumor sizes were equally represented. The animals received four tail vein injections, at the intervals of 7 days, of the following preparations: control (Hepes-buffered saline, pH 6.5) or nal-IRI at a dose of 5 mg/kg per injection. The liposomal ATR inhibitor was administrated 24 hours later at a dose of 20 mg/kg per injection. Liposomes for injections were prepared as described above. The animal weight and tumor size were monitored twice weekly. The tumor progression was monitored by palpation and caliper measurements of the tumors along the largest (length) and smallest (width) axes twice a week. The tumor sizes were determined twice weekly from the caliper measurements using the formula: tumor volume = [(length) × (width)2]/2. To assess treatment-related toxicity, the animals were also weighted twice weekly. When the tumors in the group reached 10% of the mouse body weight, the animals in the group were euthanized.


Fig. S1. Additional images from live-cell imaging of U2OS cell line.

Fig. S2. Results of 500 independent model fits from randomized starting parameters.

Fig. S3. Identifiability analysis of model parameter by analysis of likelihood profiles.

Fig. S4. Drug response profiles for all chemotherapy and targeted inhibitors predicted by the computational model.

Fig. S5. Additional, experimentally measured drug response profiles for chemotherapy and targeted inhibitors.

Fig. S6. Western blot data for U2OS and MS-751 cell lines treated with chemotherapies, ATR inhibitor, or combinations.

Fig. S7. Additional drug response profiles comparing the potentiating effect of gemcitabine and ATR inhibition between the U2OS cell line and a panel of lung cancer cell lines.

Table S1. Quantification of model misfit for individual experimental conditions.

Table S2. Chemotherapy drugs and targeted inhibitors used.

Table S3. Antibodies used for Western blotting.

Movie S1. Time-lapse video of U2OS cells treated with a combination of gemcitabine and ATR inhibitor.

Data file S1. Model parameters and computational analyses.


Acknowledgments: We thank M. Lee, A. Hsu, S. Morandell, and A. Tentner for initial consulting on DNA damage signaling and for providing a list of previously evaluated antibodies and cell lines. We also thank the internal review team at Merrimack Pharmaceuticals for reviewing and helping improve the manuscript. Funding: No external funding was used to finance this study. Author contributions: M.S., T.H., and A.R. created the computational model and performed all in silico analyses. O.A. performed all in vitro experiments. A.K. performed the in vivo experiments. B.S., D.C.D., M.B.Y., O.A., and A.R. designed the study. All co-authors wrote the manuscript and created the figures. Competing interests: O.A., B.S., M.S., A.K., T.H., D.C.D., and A.R. were employed by Merrimack Pharmaceuticals Inc. when they contributed to this work. M.B.Y. was a scientific advisor for Merrimack Pharmaceuticals Inc. when he contributed to this work and owns stock in the company. Data and materials availability: All experimental data, the computational model, and MATLAB code to run model simulations are freely available for download ( The computational model was also made available on the BioModels Database ( All other data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials.
View Abstract

Navigate This Article