Pharmacodynamic Modeling of Cell Cycle Effects for Gemcitabine and Trabectedin Combinations in Pancreatic Cancer Cells

Combinations of gemcitabine and trabectedin exert modest synergistic cytotoxic effects on two pancreatic cancer cell lines. Here, systems pharmacodynamic (PD) models that integrate cellular response data and extend a prototype model framework were developed to characterize dynamic changes in cell cycle phases of cancer cell subpopulations in response to gemcitabine and trabectedin as single agents and in combination. Extensive experimental data were obtained for two pancreatic cancer cell lines (MiaPaCa-2 and BxPC-3), including cell proliferation rates over 0–120 h of drug exposure, and the fraction of cells in different cell cycle phases or apoptosis. Cell cycle analysis demonstrated that gemcitabine induced cell cycle arrest in S phase, and trabectedin induced transient cell cycle arrest in S phase that progressed to G2/M phase. Over time, cells in the control group accumulated in G0/G1 phase. Systems cell cycle models were developed based on observed mechanisms and were used to characterize both cell proliferation and cell numbers in the sub G1, G0/G1, S, and G2/M phases in the control and drug-treated groups. The proposed mathematical models captured well both single and joint effects of gemcitabine and trabectedin. Interaction parameters were applied to quantify unexplainable drug-drug interaction effects on cell cycle arrest in S phase and in inducing apoptosis. The developed models were able to identify and quantify the different underlying interactions between gemcitabine and trabectedin, and captured well our large datasets in the dimensions of time, drug concentrations, and cellular subpopulations.


INTRODUCTION
Pancreatic cancer is a highly aggressive malignancy and shows resistance to almost all existing treatments (Oberstein and Olive, 2013). Gemcitabine (GEMZAR, Eli Lilly, Indianapolis, IN), a standard therapy for the treatment of advanced pancreatic cancer, can disrupt DNA replication and activate the S phase checkpoint (Yip-Schneider et al., 2001;Morgan et al., 2005; Abbreviations: CCM, cell cycle model; DDI, drug-drug interaction; PK/PD, pharmacokinetics/pharmacodynamics. Robinson et al., 2006). However, the benefits of gemcitabine monotherapy are limited, and combinations of other agents with gemcitabine may improve survival of pancreatic cancer patients. Trabectedin (YONDELIS R , Et-743; Johnson and Johnson Pharmaceutical Research and Development, Raritan, NJ, USA; PharmaMar S.A.U., Madrid, Spain) is a promising anticancer agent that has demonstrated clinical activity in many drug-resistant cancer cell lines, and has been approved by the US Food and Drug Administration for advanced soft tissue sarcoma. It has three tetrahydroisoquinoline rings. The A and B subunits bind covalently to the DNA minor groove and bend DNA toward the major groove, and the C ring protrudes to interact with adjacent macromolecules such as transcription factors (D'Incalci and Galmarini, 2010). Trabectedin was found previously to cause cell cycle arrest at S and G 2 /M phases in many human tumor cell lines (Gajate et al., 2002;Simoens et al., 2003). Because of its unique mechanisms of action (D'Incalci and Galmarini, 2010), trabectedin has been reported to exert anti-tumor activities in many malignancies, including soft-tissue sarcomas, ovarian carcinomas, and breast cancer (D'Incalci et al., 2002;D'Incalci and Zambelli, 2015).
Our previous report provided indications from the literature that gemcitabine and trabectedin have mechanisms that may interrelate to produce synergism in their chemotherapeutic effects and we demonstrated that the combination of gemcitabine and trabectedin exerts synergistic cytotoxic effects on pancreatic cancer cells (Miao et al., 2016a). Here we have extended the work, assessing cell cycle subpopulations in two pancreatic cancer cell lines to examine drug interactions, because asynchronous cancer cell cultures are composed of different subpopulations, and each may have different sensitivities to drugs. Previously we also developed a pharmacodynamic (PD) model that was able to characterize simultaneously 32 sets of data for single-agent and combined drug effects on pancreatic cancer cell lines (Miao et al., 2016a). Here we have expanded the model to integrate additional data regarding the temporal changes of cell numbers in sub G 1 , G 0 /G 1 , S, and G 2 /M phases, so as to determine how each subpopulation contributes to the observed effects of the drugs, as single agents or combined.
Cell cycle models have been developed previously to characterize cell cycle arrest and induction of apoptosis for drugs such as gemcitabine (Jusko, 1973;Hamed et al., 2013;Zhu et al., 2015). In this study, we extended a cell cycle model (Hamed et al., 2013) to integrate components of our previous model (Miao et al., 2016a) in order to characterize cell cycle effects of drug combinations. We measured cell proliferation as temporal changes in total cell numbers, as well as the fraction of cells in each phase of the cell cycle, and used the absolute cell number in each cell cycle phase as a PD endpoint for model fitting and qualification. The cell cycle models feature the dimensions of time, drug concentration, and drug effects on cell subpopulations. The application of mathematical modeling of cell subpopulation responses to combination therapy, and gaining an understanding of drug effects upon the transition rates between cell cycle phases, provides a greater insight into the molecular mechanisms underlying the synergistic effects of gemcitabine and trabectedin.

Reagents
Gemcitabine hydrochloride was purchased from Eli Lilly (Indianapolis, IN), dissolved in sterile double-distilled water, and stored at −20 • C at a stock concentration of 50 mM. Trabectedin, obtained as a gift from PharmaMar (Madrid, Spain), was prepared as a 1 mM stock solution in dimethylsulfoxide (DMSO) and stored at −20 • C.

Cell Proliferation Assay
To enable exponential cell growth for the duration of the experiment, 1.5 − 2 × 10 6 cells in 5 mL fresh medium were seeded in 6-well plates and allowed to adhere overnight. Cells then were treated in triplicate with 4 different concentrations of gemcitabine and trabectedin, alone or combined, for 6 different time intervals of up to 120 h ( Table 1). Equivalent volumes of water or DMSO were added as vehicle controls for the two drugs. At the appropriate times, triplicate wells of cells were washed twice with Dulbecco's phosphate buffered saline (Gibco) and suspended by incubating in 500 mL 1×Trypsin EDTA (Cellgro) for 5 min. Cells were counted using a Coulter Counter model Z2 (Beckman Coulter, Hialeah, FL), and analyzed by flow cytometry.

Cell Cycle Assay
Aliquots of the cell suspensions were fixed in 70% cold ethanol (Decan Laboratories, King of Prussia, PA) and stored at −20 • C until flow cytometry analysis, which was performed within 2 weeks. In brief, the ethanol was removed and cells were washed in cold Stain Buffer (BD Pharmingen, San Diego, CA) and treated for 30 min at room temperature in the dark with propidium iodide (PI) containing RNase (BD Pharmingen). Based upon analysis of the DNA content of each cell, histograms of cell distribution in the different cell cycle phases were obtained using a FACSort flow cytometer (Becton Dickinson, San Jose, CA). The CellQuest, WinList, and ModFit LT 4.0 software (Verity Software, Topsham, ME) were used for analysis of cell cycle distribution and determination of the fraction of cells in the G 0 /G 1 , S, and G 2 /M phases. Apoptosis was measured by quantifying the sub G 1 peak. Measurement accuracy was evaluated from the coefficient of variation (CV%). Each sample was assayed in triplicate. The number of cells in each cell cycle phase was obtained by multiplying the total cell numbers of each sample by the fraction of cells in each phase.

Mathematical Model
The models were premised on known cell cycle processes with components informed by observations of perturbations caused by the drugs. Mathematical models (Figure 1) were developed based on a series of ordinary differential equations using a stepwise modeling approach to characterize the experimental data.
In the first step, the cell cycle base model was constructed with saturation of cell numbers and contact inhibition of cell proliferation included to describe the vehicle control data. In the second step, the cell cycle base model was extended to include estimates of the cytostatic and cytotoxic effects of the single drugs. Model components for drug-and cell line-specific behaviors over time and concentration were included. Finally, the cell cycle models of gemcitabine and trabectedin were combined, and combination drug effects were incorporated in the model. Parameters estimated from previous steps were fixed when performing model fitting in the subsequent steps.

Cell Cycle Base Model
A cell cycle base model ( Figure 1A; Hamed and Roth, 2011;Hamed et al., 2013) was used previously to characterize the time course of three cell cycle phases in the absence of drug: where G 1 is the number of cells in G 0 /G 1 phase, S is the cell number in S phase, and G 2 are cells in G 2 /M phase. In our study, neither G 0 and G 1 nor G 2 and M were distinguishable experimentally. The total cell number is R tot and the first-order transition rates among the consecutive phases are k 1 , k 2 and k 3 . Equations (1)-(3) are a linear system and all phases have exponential behavior without saturation. The doubling time can be approximated by T d = 1/k 1 + 1/k 2 + 1/k 3 . Saturation of total cell growth toward a maximal cell count at steady-state R ss was introduced by slowing the doubling process in Equation (1) via: with this modification T d holds only for early times. Please note that multiplication of R ss by 2 is necessary to adjust for the doubling factor at G2, see Appendix for details. To account for cell contact inhibition, the outflow from G 0 /G 1 was modified by slowing the transition rate k 1 . The final cell cycle model ( Figure 1B) without drug effects is then: where I max is the maximal growth inhibition effect of cell contact, IT 50 is the time when the half-maximal inhibition effect is achieved, and γ is the Hill coefficient. For increasing times, the states of the three cell phases in Equations (6)-(9) converge toward: The Appendix provides derivations of Equations (10) and (11).

Cell Cycle Model of Gemcitabine and Trabectedin As Single Agents
Gemcitabine and trabectedin have concentration-and timedependent effects on the transition rates. In addition, the MiaPaCa-2 and BxPC-3 cell lines showed different behaviors, and appropriate model adjustments were necessary for these factors. To avoid repetitive representation of the specific model equations, we present a general model structure and list only the necessary adjustments for the different drugs and cell lines. The general model is: Frontiers in Pharmacology | www.frontiersin.org where C Gem and C Tbd are the concentrations of gemcitabine and trabectedin, A is the cell number in the sub G 1 phase, and k apop is the first-order rate constant for apoptosis. Drug and cell-line effects on the transition rates are modeled by the effect functions e Z,1 , e Z,2 , and e Z,3 for gemcitabine and trabectedin, where Z denotes either the MiaPaCa-2 (M) or BxPC-3 (B) cell line. Apoptotic drug effects are described by on/off functions of the form: where i = 1 indicates that a delayed apoptotic effect occurs as described by: where K X represents the non-linear cytotoxicity function, K j,X are transit steps, K max,X is the maximum killing rate, KC 50,X is the sensitivity constant, and γ X is the Hill coefficient. X indicates either Gem or Tbd. If no apoptotic effect exists, i = 0 in Equation (17). Four transit steps were previously found (Miao et al., 2016a) optimal to describe the time delay.

Gemcitabine
The general cell cycle models structure described above was refined to characterize gemcitabine effects on MiaPaCa-2 ( Figure 1C) and BxPC-3 ( Figure 1D) cells. Gemcitabine-induced cell cycle arrest was modeled as inhibition of the transition rate between S and G2/M phase. It was assumed that gemcitabine induced apoptosis from S phase for MiaPaCa-2 cells and both G 0 /G 1 , and S phases for BxPC-3 cells. This accounts for the greater sensitivity of BxPC-3 cells to gemcitabine compared to MiaPaCa-2 cells (Miao et al., 2016a).

Model Equations for Gemcitabine on MiaPaCa-2 Cells
The transition rate k 1 increases with gemcitabine concentrations as described by: where α k1 is a rate constant, and C ref is a reference concentration. In this study, C ref is always set to the lowest gemcitabine concentration in all subsequent terms. As a result, k 1 will increase exponentially with increasing gemcitabine concentrations. Gemcitabine-induced cell cycle arrest in S phase is modeled by inhibition of transition rate k 2 with: where I max, Gem represents gemcitabine maximum inhibition in S phase, IC 50, Gem is the gemcitabine concentration inducing 50% of cell cycle arrest, and γ Gem1 is the Hill coefficient. Of note, the cell cycle arrest in S phase is concentration-dependent, and thus an exponential term with rate of k is used to describe increased inhibition as gemcitabine concentrations increase. C ref is always set to the lowest gemcitabine concentration. The C ref was used for gemcitabine not for trabectedin. This is because the concentration-response curve for gemcitabine is more gradual compared to trabectedin, and the gemcitabine concentrations were chosen at about 0, 1 2 IC 50 , IC 50 and 2IC 50 values for both cell lines, thus the concentration-dependency is more prominent at high concentrations for gemcitabine. For trabectedin, the concentration-response curve is steeper, so the concentration-dependency is more prominent and there is no need to incorporate C ref . The exponential term has values between 0 and 1. The transition rate k 3 is not affected by time and concentration, and we set: Apoptotic effects are assumed to occur in the S phase only, and therefore we have i Gem  (21)-(23), and the schematic is depicted in Figure 1C.

Model Equations for Gemcitabine on BxPC-3 Cells
No effect of gemcitabine was observed on k 1 and we set e B,1 (C Gem , 0) = 1.
For early time points all drug-exposed cells are arrested in S phase. But this effect wanes with longer exposure times. Therefore, Equation (22) is extended by an exponential term, and the time-and concentration-dependent effect of gemcitabine on k 2 is given by: The exponential term was applied to k 2 to describe k 2 increases with rate constant α k2 . This is based on the assumption that gemcitabine perturbs k 2 by changing cell cycle regulation such as cyclin-CDK protein expression over time (Yip-Schneider et al., 2001).
In addition, we assume that k 3 increases as gemcitabine concentrations increase, as described by: Apoptosis occurs in the G 1 and S phases, and we set i Gem (12)-(20) with Equations (24)-(26). The schematic is depicted in Figure 1D.

Trabectedin
Trabectedin-induced cell cycle arrest was modeled as inhibition of the transition from S to G 2 /M phases and from G 2 /M to G 0 /G 1 phases, and cells in all cycle phases may commit to apoptosis. A diagram of the cell cycle model for trabectedin is shown in Figure 1E.

Model Equations for Trabectedin on MiaPaCa-2 Cells
The transition rate k 1 increases with increasing trabectedin concentration by: The transition rate k 2 increases over time by: whereas k 3 is inhibited over time:

Model Equations for Trabectedin on BxPC-3 Cells
Transition rates k 1 and k 2 were modeled with Equations (27) and (28). We assumed that the inhibition of transition rate k 3 increases over time and k 3 exhibits concentration-dependency:

Cell Cycle Model of Drug Combinations
The time-and concentration-dependent drug effects of gemcitabine and trabectedin as single agents on k 1 and k 3 are multiplied for the drug combination. Both drugs interact at the transition from S to G2/M phase and at the induction of apoptosis. Instead of multiplying the single effects, we assume competitive drug interactions for cell cycle arrest in S phase and apply the competitive combination effect term from Ariëns et al. (Ariens et al., 1957;Koch et al., 2016). For the induction of apoptosis, the effects of both cytotoxic drug effects are summed (Miao et al., 2016a). To account for remaining interaction effects not predicted by the model, two interaction parameters, ψ 1 and ψ 2 (Chakraborty and Jusko, 2002;Koch et al., 2009) were included at the points where the drugs interact on cell cycle inhibition and induction of apoptosis. If ψ 1 or ψ 2 equals 1, the effect of the combination is additive, based on the applied combination effect term. A ψ 1 or ψ 2 value smaller than 1 indicates synergistic interaction, and any value greater than 1 indicates antagonistic combination behavior.

Model Equations for Combinations on MiaPaCa-2 Cells
The drug interaction on k 1 was modeled with multiplication of single drug effects. For k 1 we obtain: (31) Because both gemcitabine and trabectedin act on the transition from S to G 2 /M, we assumed a competitive behavior and multiplied the interaction parameter ψ 1 by the IC 50 of gemcitabine: and set e M,2 (C Gem , C Tbd ) = E (C Gem , C Tbd ) .
In Equation (18), ψ 2 is multiplied by the KC 50,Gem of gemcitabine. Because gemcitabine has no effect on k 3 , (compare Equation 23

Model Equations for Combinations on BxPC-3 Cells
For k 1 we apply Equation (27). The transit rate k 2 is influenced by the competitive effect Equations (32) and a time-dependent effect, (compare Equation 25): For k 3 : The model equations are Equations (12)

Data Analysis
The modeling was performed using ADAPT 5 software (D'Argenio et al., 2009) with a naive pooled approach. Models were fitted to the data using the maximum likelihood (ML) estimation method. The variance model was: where V i is the variance of the ith time point, σ 1 and σ 2 are variance model parameters, and Y i is the predicted response at ith time point. The goodness of fit was assessed by visual inspection of model fittings, goodness-of-fit plots, the Akaike Information Criteria (AIC), and the coefficients of variation (CV %).

RESULTS
Cell Cycle without Perturbation (Baseline Model) Figure 2 shows the cell cycle distribution of MiaPaCa-2 (Figures 2A,B) and BxPC-3 cells (Figures 2C,D) in the absence of drug. Both cell lines showed progressive accumulation in the G 0 /G 1 phase with time (Figures 2A,C). Numerous factors can cause accumulation of the cell population in the G 0 /G 1 phase. Serum starvation, a reduction in nutrients, or contact inhibition as a result of increasing cell confluence activate cell growth checkpoints and cell cycle arrest occurs at the G 0 /G 1 phase (Hayes et al., 2005;Choresca et al., 2009;Dalman et al., 2010). Total cell numbers and the number of cells in different cycle phases, were modeled simultaneously with the cell cycle base model (Figure 1B). The model was able to characterize well the observed data (Figures 2B,D). Parameter estimates are shown in Table 2. These model estimates were then fixed for subsequent analyses. The approximate doubling time T d for MiaPaCa-2 was 31.9 and 16.8 h for BxPC-3 cells. Cells finally reached a steady-state R ss , which was 6.00 × 10 6 cells (fixed to a value observed in the study, data not shown) for MiaPaCa-2 and 7.98 × 10 5 cells for BxPC-3 cells. Under the experimental conditions used here, the time to reach half of the maximal cell contact inhibition was 45.4 h for MiaPaCa-2 and 61.0 h for BxPC-3 cells.  Table 1). A gradual increase of cells in S phase also was observed after exposure of BxPC-3 cells to 34 nM gemcitabine ( Figure 3D and Supplemental Table 1). Gemcitabine-induced S phase cell cycle arrest was concentrationdependent for both MiaPaCa-2 ( Figure 3B) and BxPC-3 cells ( Figure 3E). Apoptosis, measured as the percentage of  sub-diploid cells, appeared at the highest-tested gemcitabine concentrations for both cell lines. Based upon mechanisms of gemcitabine effects on the cell cycle, mathematical models were developed for MiaPaCa-2 ( Figure 1C) and BxPC-3 cells ( Figure 1D). Figures 3C,F show the model fittings. In general, the model well captured the trend of the observed data. Considering the large number and complexities of the datasets, model performance is acceptable. Table 3 shows parameter estimates for the drugs as single agents. Gemcitabine-induced S phase arrest was modeled as inhibition of the transition rate from S to G2/M phase. The maximum inhibition I max,Gem was fixed to 1 for both cell lines. The concentrations that induced half-maximal S phase arrest IC 50,Gem were 11.5nM for MiaPaCa-2 and 48.1 nM for BxPC-3. Because the maximal inhibition I max,Gem was concentration-dependent, inhibition increased proportional to the rate constant k, which was 0.0281 nM −1 for MiaPaCa-2 and 1 nM −1 for BxPC-3 cells. The model assumed that gemcitabine-induced apoptosis occurred in S phase in MiaPaCa-2 cells, and in both G 0 /G 1 and S phases in BxPC-3 cells, and the parameters related to gemcitabine killing effects were fixed to parameters derived in our previous study (Miao et al., 2016a). Minor differences between MiaPaCa-2 and BxPC-3 cells also included certain model assumptions; for MiaPaCa-2, it was assumed that k 1 increased with rate of α k 1 (0.0612 nM −1 ) with increased gemcitabine concentrations. For BxPC-3, it was assumed that k 2 and k 3 increased with rates of α k2 (0.0635 h −1 ) and α k3 (0.0605 nM −1 ) as gemcitabine concentrations increased. Figure 4 shows the effects of trabectedin on the cell cycle distribution for MiaPaCa-2 (Figures 4A-C) and BxPC-3 cells (Figures 4D-F). For both cell lines, trabectedin induced cell cycle arrest at S phase at early time points, followed by progression to G 2 /M phase at later time points. Figure 4A shows the cell cycle distribution of MiaPaCa-2 cells after exposure to 0.8 nM trabectedin for 0, 24, 48, 72, and 96 h. At 24 h, cell cycle arrest occurred in the S phase (Supplemental Table 1). Similar effects were observed in BxPC-3 cells (Figure 4D and Supplemental Table 1). Trabectedin effects on the cell cycle were concentrationdependent (Figures 4B,E), for both cell lines at early time points, and S phase accumulation increased with concentration. As exposure times increased (48-76 h) the accumulation shifted from S to G 2 /M phase.

Trabectedin Effects upon Cell Cycle and Model Prediction
Based upon the mechanistic effects of trabectedin on the cell cycle, mathematical models were developed ( Figure 1E). Trabectedin effects on cell cycle checkpoints were modeled with inhibition effects on both S to G 2 /M phase transition k 2 and the G 2 /M to G 0 /G 1 phase transition k 3 . In order to characterize early S phase arrest followed by subsequent G 2 /M phase arrest, models were constructed in such a way that S phase inhibition decreased with time with rate k S and inhibition in G 2 /M phase increased with time with rate k G2 . The k S was 1 h −1 for MiaPaCa-2 and 0.00204 h −1 for BxPC-3, whereas k G2 was 0.0270 h −1 for MiaPaCa-2 and 0.0516 h −1 for BxPC-3 cells. Trabectedin concentrations inducing half-maximal S phase inhibition (IC 50,Tbd,S ) were 0.222 nM for MiaPaCa-2 and 1.03 nM for BxPC-3 cells. Concentrations inducing half-maximal G 2 /M phase inhibition (IC 50,Tbd,G2 ) were 0.525 nM for MiaPaCa-2 and 0.681 nM for BxPC-3 cells, suggesting that MiaPaCa-2 cells are also more sensitive to trabectedin-induced G 2 /M phase arrest than BxPC-3 cells. In addition, we also assumed concentrationdependent trabectedin effects on k 1 , with rates of β k1 of 2.64 nM −1 for MiaPaCa-2 and 0.637 nM −1 for BxPC-3. For BxPC-3, k 3 also increased with trabectedin concentration, with a rate of β k3 of 0.296 nM −1 . Finally, trabectedin-induced apoptosis was assumed to occur in all cell cycle phases. For both cell lines, the parameters related to trabectedin killing effects were fixed to values obtained previously (Miao et al., 2016a).

Modeling Combined Drug Effects upon Cell Cycle
Analysis of experimental data suggested that trabectedin and gemcitabine exert different effects upon cell cycle progression. For both cell lines, gemcitabine induced cell cycle arrest in S phase (Figure 3), but with trabectedin a population of cells passed through S phase and accumulated in G 2 /M phase cells  Figure 1C (MiaPaCa-2) and Figure 1D (BxPC-3).
( Figure 4). In order to model combined drug effects upon the cell cycle, the PD models for the single agents were integrated (Figures 1F,G) and used to fit simultaneously the data for four different drug concentration combinations on MiaPaCa-2 ( Figure 5A) and BxPC-3 cells ( Figure 5B). Figure 5C shows the BxPC-3 cell growth model fitting of experimental data for BxPC-3 cells. The model captured the trend of the data well. Table 4 provides parameter estimates and statistics for the combined drugs. For both cell lines, a competitive interaction equation was used to model S phase inhibition by the two drugs. The interaction term ψ 1 was multiplied at IC 50,Gem . The ψ 1 interaction parameter was fixed to 1 for MiaPaCa-2 cells (because modeling of ψ 1 for S phase inhibition resulted in a value of about 1) and estimated as 0.975 for BxPC-3 cells, suggesting additive drug interaction on S phase inhibition. The ψ 2 was multiplied by KC 50,Gem to characterize drug interactions for induction of apoptosis. Estimates of ψ 2 were 0.499 for MiaPaCa-2 cells and 0.363 for BxPC-3 cells, indicating synergistic interactions in inducing apoptosis.

DISCUSSION
Mathematical models serve to help integrate information from complex studies, assess mechanisms of drug interactions, and guide selection and optimization of drug combinations. We previously introduced a mathematical framework based on a semi mechanism-based approach to model gemcitabine and trabectedin combination effects on pancreatic cancer cells (Miao et al., 2016a). That simple and reasonable approach could be applied to other drugs, and could also be useful to characterize time-and concentration-dependent drug-drug interactions. The present work goes beyond such semi-mechanistic models. It incorporates cell subpopulation information and enlarges the previous model in order to characterize drug combination effects of gemcitabine and trabectedin on cell cycle distribution.
Cancer cells are heterogeneous and asynchronous populations are composed of cell subpopulations in different cycle phases (Evan and Vousden, 2001). Drug sensitivity or mechanism may vary for cells in different cycle phases. Cell cycle arrest at different phases could influence cell sensitivity or resistance to drugs significantly and affect overall responses to single or multiple chemotherapeutic agents. To account for such drug effects, we developed four cell cycle models: a cell cycle base model lacking drug effects, cell cycle models for gemcitabine and trabectedin as single agents, and a cell cycle model for drug combinations. The models were applied to large experimental datasets that captured relevant pharmacodynamic endpoints such as cell proliferation, cell cycle phase, and induction of apoptosis. A step-wise modeling approach was used, first fitting the control (vehicle) data to the cell cycle base model (Figure 2), followed by fitting single drug effects to the cell cycle models of gemcitabine (Figure 3) and trabectedin (Figure 4). Ultimately, the combined data were fitted with a cell cycle model of the drug combination ( Figure 5). The total cell number was measured during drug exposure, as was the fraction of cells in each cycle phase (G 0 /G 1 , S, G 2 /M) and in apoptosis (sub G 1 ). The developed model incorporates not only changes in cell cycle distribution, but also the dynamics of the whole cell population. The dataset was large, consisting of ∼500 samples that each provided data for 5 cell populations, and was analyzed using ordinary differential equations.
The proposed cell cycle model framework was developed based on drug mechanisms observed in the study. For the vehicle-treated controls, the data show that the G 0 /G 1 fraction increased with time, reflecting the fact that when cells reach confluence, contact inhibition occurs, nutrients may become limited, and a cell growth checkpoint is activated toward the end of G 1 phase. If cells are able to progress to S phase, normal cycling occurs. Otherwise, cells enter a resting G 0 phase until they are able to resume cycling. For the vehicle control group, the logistic growth function was used to characterize saturation of total cells, while contact inhibition was used to characterize G 0 /G 1 phase cells increasing over time. The contact inhibition was modeled with a time-dependent Hill function. We have also tested contact inhibition with a feedback loop where the inhibition of k 1 is dependent on R tot . However, the time-dependent Hill type of equation was chosen as it allows us to estimate the time where the half-maximal effects of cell contact inhibition occurs and to simplify the convergence analysis as presented in the Appendix. Our model for gemcitabine was initially premised on that of Hamed et al. Parameters were obtained after fitting to cell cycle models for combined drug interaction effects ( Figures 1F,G).
(2013), but more extensive measurements provided clearer insights into drug effects on apoptosis. Further, for gemcitabine, time-and concentration-dependent cell cycle arrest occurred at the S phase, consistent with previous results (Yip-Schneider et al., 2001;Morgan et al., 2005). The gemcitabine-induced S phase arrest was modeled as both an increase in the transition rate from G 0 /G 1 and an inhibition of transition rates from S phase to G 2 /M phase. The S phase arrest increased with gemcitabine concentrations. Therefore, gemcitabine-induced cell death is modeled with a non-linear killing function, with the delay of killing effects modeled by four transit compartments. Trabectedin induced cell cycle arrest in S phase at earlier time points, but at later times, cells progressed through S to G 2 /M phase. In order to capture the time-dependency, equations were constructed such that an inhibition term on the transition from S to G 2 /M decreased with time, and the inhibition term on the transition from G 2 /M to G 0 /G 1 increased with time. In addition, we assumed that cells in all three cell cycle phases can commit to apoptosis after drug treatment. The models employed a similar basic structure for both MiaPaCa-2 and BxPC-3 cells, with only minor changes required to reflect the unique characteristics of the specific cell line. For example gemcitabine induced apoptosis in S phase for MiaPaCa-2 cells, but apoptosis occurred in BxPC-3 cells in both S and G 0 /G 1 phases. The models were constructed with differences between cell lines. For example, based on the KC 50 of the two cell lines (Miao et al., 2016a), MiaPaCa-2 is more trabectedin-sensitive while BxPC-3 is more gemcitabinesensitive. In addition, the two cells lines also have phenotypic and genotypic differences; for example, MiaPaCa-2 has mutant K-ras and p53 genes, whereas BxPC-3 has wild-type K-ras and p53 genes (Deer et al., 2010). The final model permitted a higher-resolution investigation of the basis of the synergistic interaction of these drugs than previously (Miao et al., 2016a). An interaction term ψ was incorporated into specific mechanistic model components to explore model behavior when interactions in those mechanisms were hypothesized to explain observed drug effects of the combination that exceeded model predictions. In the final developed model, ψ 1 was applied to represent drug interactions at the point of S phase arrest, and ψ 2 was applied to represent drug interactions at the point of induction of apoptosis. However, ψ 1 was approximately equal to 1, suggesting that drug interactions in inducing S phase accumulation were additive, whereas ψ 2 was <1, suggesting that drug interactions in inducing apoptosis were synergistic.
Cell-cycle-specific compartmental models have been applied previously for investigation of cancer chemotherapy agents. Simple models can be constructed with as few as 4-5 compartments that represent each cycle phase (Kozusko et al., 2001;Florian et al., 2005;Ribba et al., 2009). Compartments can be added to represent apoptosis or other death mechanisms (Basse et al., 2004;Panetta et al., 2006;Sherer et al., 2006;Zhu et al., 2015), and additional components may be included to account for more complex pharmacological mechanisms such as cell cycle regulation (Senderowicz, 2004;Ferrell et al., 2011). Most existing cell cycle models characterize single drug effects; few reports model drug combinations (Gardner, 2002;Zhu et al., 2015). Because of the layers of complexities in larger and richer data sets, there is a need for models to account for those complexities. The cell cycle models developed here allow investigation of the specific effects of the individual drugs in combination therapy upon cell cycle progression. For example, we identified that gemcitabine exerts influence in the S phase, whereas trabectedin exerts influence in both S and G 2 /M phases. In addition, the models have the flexibility to account for conditions under which drug effects upon specific cycle phases or phase transitions are concentration-or timedependent. For example, the degree of gemcitabine-induced S phase arrest increased with concentration. Finally, the developed model was able to account for cell line differences; the two cell lines investigated are different in phenotype and genotype (Deer et al., 2010), and in their sensitivity to gemcitabine and trabectedin (Miao et al., 2016a). Model development showed that to capture differences in cell line-specific behaviors, drug effects may be hypothesized to affect different, or even multiple cycle phase transition rates, and apoptosis may be induced from different cell cycle phases. Traditional approaches to model tumor growth often include a logistic or Gompertzian growth model. Such models are empirical and do not take into account complex biological mechanisms such as the cell cycle dynamics of subpopulations. Here we used models based upon cell cycle structure to describe cancer cell growth and drug effects. The final model characterized very well not only cell cycle dynamics, but also cancer cell proliferation ( Figure 5C). The model can characterize cancer cell growth without drug perturbation, and growth in response to single-or combined cell-cycle-specific agents.
In conclusion, we have modeled chemotherapeutic drug effects of single and combined agents on pancreatic cancer cell cycle dynamics and apoptosis. The study suggests a basis for the overall gemcitabine/trabectedin synergy observed previously (Miao et al., 2016a) when examines cell cycle subpopulations. The proposed model structure is not limited to gemcitabine and trabectedin as single or combined agents, and can be adapted to investigate efficacy of other cell-cycle specific agents. Even though the cell cycle models developed in this study are used to characterize in vitro data, the model structures may be kept and integrated with more model components accounting for cell and tissue heterogeneity when translating to in vivo work. However, the mechanisms underlying synergy within specific cell subpopulations remains undefined at the level of protein-and pathway interactions.
Further research, employing approaches such as high-resolution proteomic investigation of drug action, may provide the means to integrate information at that scale into the model structure proposed here.
The model components and complications that were introduced in this report for the dual drug effects on MiaPaCa-2 cells have been largely confirmed in our further studies that utilized proteomic and western blot methods to assess drug effects on various signaling pathways (Miao et al., 2016b). The present report reflects the second part of our multistage efforts to perform studies that employ modeling to assess relevance of diverse physiologic and pharmacologic processes and then employ more sophisticated methodology to delve into possible and identified mechanisms and complexities.