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 G₂/M phase. Over time, cells in the control group accumulated in G₀/G₁ 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 G₁, G₀/G₁, S, and G₂/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; 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^®, 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₂/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₁, G₀/G₁, S, and G₂/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.

Materials and Methods

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 Culture

The pancreatic cancer cell lines MiaPaCa-2 and BxPC-3 were purchased from American Type Culture Collection (ATCC). MiaPaCa-2 cells were grown in DMEM (Cellgro, Manassa, VA) supplemented with 10% (v/v) fetal bovine serum (Cellgro). BxPC-3 cells were cultured in RPMI (Cellgro), 10% (v/v) fetal bovine serum and 1% (v/v) sodium pyruvate (Gibco, Grand Island, NY). Cells were maintained in 5% CO₂ at 37°C with 95% humidity and grown as monolayers in T75 tissue culture flasks (BD Biosciences, Bedford, MA).

Cell Proliferation Assay

To enable exponential cell growth for the duration of the experiment, 1.5 − 2 × 10⁶ 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.

TABLE 1

Table 1. Concentrations of gemcitabine (G) and trabectedin (T) tested as single agents and in combination in cell culture studies.

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₀/G₁, S, and G₂/M phases. Apoptosis was measured by quantifying the sub G₁ 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 step-wise 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.

FIGURE 1

Figure 1. Schematic of the pharmacodynamic models for single and combination drug effects on the cell cycle. (A) Cell cycle base model. (B) Alternative cell cycle base model. (C) Cell cycle model of gemcitabine for MiaPaCa-2 cells. (D) Cell cycle model of gemcitabine for BxPC-3 cells. (E) Cell cycle model of trabectedin for MiaPaCa-2 and BxPC-3 cells. (F) Cell cycle model of drug combination for MiaPaCa-2 cells. (G) Cell cycle model of drug combination for BxPC-3 cells.

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:

$\begin{array}{ll}\frac{d{G}_1}{dt}=2k_3{G}_2-k_1{G}_1{G}_1(0)=G_1^0& (1)\end{array}$

$\begin{array}{ll}\frac{dS}{dt}=k_1{G}_1-k_{2}SS(0)=S^0& (2)\end{array}$

$\begin{array}{ll}\frac{d{G}_2}{dt}=k_{2}S-k_3{G}_2{G}_2(0)=G_2^0& (3)\end{array}$

$\begin{array}{ll}{R}_{tot}=G_1+S+G_2& (4)\end{array}$

where G₁ is the number of cells in G₀/G₁ phase, S is the cell number in S phase, and G₂ are cells in G₂/M phase. In our study, neither G₀ and G₁ nor G₂ and M were distinguishable experimentally. The total cell number is R_tot and the first-order transition rates among the consecutive phases are k₁, k₂ and k₃. 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/k₂ + 1/k₃.

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:

$\begin{array}{ll}\frac{d{G}_1}{dt}=2k_3{G}_2\left(1-\frac{R_{tot}}{{2R}_{ss}}\right)-k_1{G}_1{G}_1(0)=G_1^0& (5)\end{array}$

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₀/G₁ was modified by slowing the transition rate k₁. The final cell cycle model (Figure 1B) without drug effects is then:

$\begin{array}{l}\frac{d{G}_1}{dt}=2k_3{G}_2\left(1-\frac{R_{tot}}{2R_{ss}}\right)-k_1\left(1-\frac{I_{max}{t}^{\gamma }}{I{T}_{50}^{\gamma }+t^{\gamma }}\right)G_1\\ G_1(0)=G_1^0& (6)\end{array}$

$\begin{array}{ll}\frac{dS}{dt}=k_1\left(1-\frac{I_{max}{t}^{\gamma }}{I{T}_{50}^{\gamma }+t^{\gamma }}\right)G_1-k_{2}SS(0)=S^0& (7)\end{array}$

$\begin{array}{ll}\frac{d{G}_2}{dt}=k_{2}S-k_3{G}_2{G}_2(0)=G_2^0& (8)\end{array}$

$\begin{array}{ll}{R}_{tot}=G_1+S+G_2& (9)\end{array}$

where I_max is the maximal growth inhibition effect of cell contact, IT₅₀ 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:

$\begin{array}{ll}{G}_1^{*}=\frac{R_{ss}}{1+\frac{k_1}{k_2}\left(1-I_{max}\right)+\frac{k_1}{k_3}\left(1-I_{max}\right)},& (10)\end{array}$

$\begin{array}{ll}{S}^{*}=\frac{k_1}{k_2}{\left(1-I_{max}\right)G}_1^{*}and{G}_2^{*}=\frac{k_1}{k_3}{\left(1-I_{max}\right)G}_1^{*}.& (11)\end{array}$

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 time-dependent 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:

$\begin{array}{cc}\begin{array}{l}\frac{d{G}_1}{dt}=2k_3{e}_{Z,3}(C_{Gem},C_{Tbd})G_2(1-\frac{R_{tot}}{2R_{ss}})G_1(0)=G_1^0\\ -k_1{e}_{Z,\text{1}}(C_{Gem},C_{Tbd})(\text{1}-\frac{I_{max}{t}^{\gamma }}{I{T}_{50}^{\gamma }\text{+}{t}^{\gamma }})G_1\\ -a_{Gem}(i_{G1}^{Gem})G_1-a_{Tbd}(i_{G1}^{Tbd})G_1\end{array}& (12)\end{array}$

$\begin{array}{cc}\begin{array}{l}\frac{dS}{dt}=k_1{e}_{Z,\text{1}}(C_{Gem},C_{Tbd})(\text{1}-\frac{I_{max}{t}^{\gamma }}{I{T}_{50}^{\gamma }\text{+}{t}^{\gamma }})G_{1}S(0)=S^0\\ -k_2{e}_{Z,\text{2}}(C_{Gem},C_{Tbd})S\\ -a_{Gem}(i_S^{Gem})S-a_{Tbd}(i_S^{Tbd})S\end{array}& (13)\end{array}$

$\begin{array}{cc}\begin{array}{l}\frac{d{G}_2}{dt}=k_2{e}_{Z,\text{2}}(C_{Gem},C_{Tbd})S-k_3{e}_{Z,\text{3}}(C_{Gem},C_{Tbd})G_2\\ G_2(0)=G_2^0\\ -a_{Gem}(i_{G2}^{Gem})G_2-a_{Tbd}(i_{G2}^{Tbd})G_2\end{array}& (14)\end{array}$

$\begin{array}{cc}\begin{array}{l}\frac{dA}{dt}=a_{Gem}(i_{G\text{1}}^{Gem})G_1\text{+}{a}_{Tbd}(i_{G\text{1}}^{Tbd})G_1\text{+}{a}_{Gem}(i_S^{Gem})S\\ \text{+}{a}_{Tbd}(i_S^{Tbd})S+a_{Gem}(i_{G\text{2}}^{Gem})G_2\text{+}{a}_{Tbd}(i_{G\text{2}}^{Tbd})G_2\\ -k_{apop}AA(\text{0})\text{= }\text{0}\end{array}& (15)\end{array}$

$\begin{array}{cc}{R}_{tot}=G_1\text{+}S\text{+}{G}_2\text{+}A& (16)\end{array}$

where C_Gem and C_Tbd are the concentrations of gemcitabine and trabectedin, A is the cell number in the sub G₁ 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:

$a_X(i)=\{\begin{array}{c}{K}_{4,X}ifi=1\\ 0ifi=0\end{array}(17)$

where i = 1 indicates that a delayed apoptotic effect occurs as described by:

$\begin{array}{cc}{K}_X=\frac{K_{max,X}{C}_X^{{\gamma }_X}}{K{C}_{50,X}^{{\gamma }_X}+C_X^{{\gamma }_X}}& (18)\end{array}$

$\frac{d{K}_{1,X}}{dt}=\frac{1}{{\tau }_X}(K_X-K_{\text{1},X})K_{\text{1},X}(0)=0(19)$

$\frac{d{K}_{j,X}}{dt}=\frac{1}{{\tau }_X}(K_{j-1,X}-K_{j,X})K_{j,X}(0)=0j=2,\dots ,4(20)$

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₀/G₁, 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₁ increases with gemcitabine concentrations as described by:

$\begin{array}{ll}{e}_{M,1}(C_{Gem},0)=exp({\alpha }_{k1}(C_{Gem}-C_{ref}))& (21)\end{array}$

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₁ will increase exponentially with increasing gemcitabine concentrations. Gemcitabine-induced cell cycle arrest in S phase is modeled by inhibition of transition rate k₂ with:

$\begin{array}{l}{e}_{M,2}(C_{Gem},0)\\ =\left(1-\frac{I_{max,Gem}{(1-exp(-k(C_{Gem}-C_{ref})))C}_{Gem}^{{\gamma }_{Gem1}}}{{IC}_{50,Gem}^{{\gamma }_{Gem1}}+C_{Gem}^{{\gamma }_{Gem1}}}\right)& (22)\end{array}$

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, $\frac{1}{2}I{C}_{50}$, IC₅₀ and 2IC₅₀ 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₃ is not affected by time and concentration, and we set:

$\begin{array}{ll}{e}_{M,3}\left(C_{Gem},0\right)=1.& (23)\end{array}$

Apoptotic effects are assumed to occur in the S phase only, and therefore we have $i_{G1}^{Gem}=i_{G2}^{Gem}=0$, $i_S^{Gem}=1$ and $i_{G1}^{Tbd}=i_S^{Tbd}=i_{G2}^{Tbd}=0$. The overall model utilizes Equations (12)–(20) with Equations (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₁ and we set

$\begin{array}{ll}{e}_{B,1}\left(C_{Gem},0\right)=1.& (24)\end{array}$

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₂ is given by:

$\begin{array}{l}{e}_{B,2}(C_{Gem},0)=(1-exp(-{\alpha }_{k2}t))\\ \left(1-\frac{I_{max,Gem}(1-exp(-k(C_{Gem}-C_{ref})))C_{Gem}^{{\gamma }_{Gem1}}}{{IC}_{50,Gem}^{{\gamma }_{Gem1}}+C_{Gem}^{{\gamma }_{Gem1}}}\right).& (25)\end{array}$

The exponential term was applied to k₂ to describe k₂ increases with rate constant α_k2. This is based on the assumption that gemcitabine perturbs k₂ by changing cell cycle regulation such as cyclin-CDK protein expression over time (Yip-Schneider et al., 2001).

In addition, we assume that k₃ increases as gemcitabine concentrations increase, as described by:

$\begin{array}{ll}{e}_{B,3}\left(C_{Gem},0\right)=\text{ exp}\left({\alpha }_{k3}\left(C_{Gem}-C_{ref}\right)\right).& (26)\end{array}$

Apoptosis occurs in the G₁ and S phases, and we set $i_{G1}^{Gem}=i_S^{Gem}=1$, $i_{G2}^{Gem}=0$ and $i_{G1}^{Tbd}=i_S^{Tbd}=i_{G2}^{Tbd}=0$. The overall model includes Equations (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₂/M phases and from G₂/M to G₀/G₁ 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₁ increases with increasing trabectedin concentration by:

$\begin{array}{ll}{e}_{M,1}\left({0,C}_{Tbd}\right)=exp\left({\beta }_{k1}{C}_{Tbd}\right).& (27)\end{array}$

The transition rate k₂ increases over time by:

$\begin{array}{ll}{e}_{M,2}\left({0,C}_{Tbd}\right)=\left(1-\frac{I_{max,Tbd,S}exp\left(-k_{S}t\right)C_{Tbd}^{{\gamma }_{Tb{d}_1}}}{I{C}_{50,Tbd,S}^{{\gamma }_{Tbd1}}+C_{Tbd}^{{\gamma }_{Tb{d}_1}}}\right)& (28)\end{array}$

whereas k₃ is inhibited over time:

$\begin{array}{ll}{e}_{M,3}\left(0,C_{Tbd}\right)=\left(1-\frac{I_{max,Tbd,G2}\left(1-exp\left(-k_{G2}t\right)\right)C_{Tbd}^{{\gamma }_{Tb{d}_2}}}{I{C}_{50,Tbd,G2}^{{\gamma }_{Tbd2}}+C_{Tbd}^{{\gamma }_{Tb{d}_2}}}\right).& (29)\end{array}$

Additionally, apoptosis can initiate in all phases by setting $i_{G1}^{Tbd}=i_S^{Tbd}=i_{G2}^{Tbd}=1$ and $i_{G1}^{Gem}=i_S^{Gem}=i_{G2}^{Gem}=0$. The overall model entails Equations (12)–(20) with Equations (27)–(29).

Model Equations for Trabectedin on BxPC-3 Cells

Transition rates k₁ and k₂ were modeled with Equations (27) and (28). We assumed that the inhibition of transition rate k₃ increases over time and k₃ exhibits concentration-dependency:

$\begin{array}{l}{e}_{B,3}({0,C}_{Tbd})=\left(1-\frac{I_{max,Tbd,G2}\left(1-exp\left(-k_{G2}t\right)\right)C_{Tbd}^{{\gamma }_{Tb{d}_2}}}{I{C}_{50,Tbd,G2}^{{\gamma }_{Tbd2}}+C_{Tbd}^{{\gamma }_{Tb{d}_2}}}\right)\\ exp\left({\beta }_{k3}{C}_{Tbd}\right).& (30)\end{array}$

Apoptosis can initiate in all phases, described by $i_{G1}^{Tbd}=i_S^{Tbd}=i_{G2}^{Tbd}=1$ and $i_{G1}^{Gem}=i_S^{Gem}=i_{G2}^{Gem}=0$. The overall model includes Equations (12)–(20) with Equations (27), (28), and (30).

Cell Cycle Model of Drug Combinations

The time- and concentration-dependent drug effects of gemcitabine and trabectedin as single agents on k₁ and k₃ 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, ψ₁ and ψ₂ (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 ψ₁ or ψ₂ equals 1, the effect of the combination is additive, based on the applied combination effect term. A ψ₁ or ψ₂ 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₁ was modeled with multiplication of single drug effects. For k₁ we obtain:

$\begin{array}{ll}{e}_{M,1}(C_{Gem},C_{Tbd})=exp\left({\alpha }_{k1}(C_{Gem}-C_{ref})\right)exp\left({\beta }_{k1}{C}_{Tbd}\right).& (31)\end{array}$

Because both gemcitabine and trabectedin act on the transition from S to G₂/M, we assumed a competitive behavior and multiplied the interaction parameter ψ₁ by the IC₅₀ of gemcitabine:

$\begin{array}{ll}E\left(C_{Gem},C_{Tbd}\right)=\left(1-\frac{I_{max,Gem}\left(1-exp\left(-k\left(C_{Gem}-C_{ref}\right)\right)\right){\left(\frac{C_{Gem}}{{{\psi }_{1}IC}_{50,Gem}}\right)}^{{{\gamma }_{Gem}}_1}+I_{max,Tbd,S}exp\left(-k_{S}t\right){\left(\frac{C_{Tbd}}{{IC}_{50,Tbd,S}}\right)}^{{\gamma }_{Tbd1}}}{1+{\left(\frac{C_{Gem}}{{{\psi }_{1}IC}_{50,Gem}}\right)}^{{{\gamma }_{Gem}}_1}+{\left(\frac{C_{Tbd}}{{IC}_{50,Tbd,S}}\right)}^{{\gamma }_{Tbd1}}}\right)& (32)\end{array}$

and set

$\begin{array}{ll}{e}_{M,2}\left(C_{Gem},C_{Tbd}\right)=E\left(C_{Gem},C_{Tbd}\right).& (33)\end{array}$

In Equation (18), ψ₂ is multiplied by the KC_{50, Gem} of gemcitabine. Because gemcitabine has no effect on k₃, (compare Equation 23), the final equations for the combination model are Equations (12)–(20) with Equations (29), (31)–(33) where we set $i_{G1}^{Gem}=i_{G2}^{Gem}=0$, $i_S^{Gem}=1$ and $i_{G1}^{Tbd}=i_S^{Tbd}=i_{G2}^{Tbd}=1$.

Model Equations for Combinations on BxPC-3 Cells

For k₁ we apply Equation (27). The transit rate k₂ is influenced by the competitive effect Equations (32) and a time-dependent effect, (compare Equation 25):

$\begin{array}{ll}{e}_{B,2}\left(C_{Gem},C_{Tbd}\right)=\left(1-exp\left(-{\alpha }_{k2}t\right)\right)E\left(C_{Gem},C_{Tbd}\right).& (34)\end{array}$

For k₃:

$\begin{array}{l}{e}_{B,3}\left(C_{Gem},C_{Tbd}\right)\\ =\left(1-\frac{I_{max,Tbd,G2}{\left(1-exp\left(-k_{G2}t\right)\right)C}_{Tbd}^{{\gamma }_{Tbd2}}}{{IC}_{50,Tbd,G2}^{{\gamma }_{Tbd2}}+C_{Tbd}^{{\gamma }_{Tbd2}}}\right)\\ exp\left({\alpha }_{k3}(C_{Gem}-C_{ref})\right)exp\left({\beta }_{k3}{C}_{Tbd}\right).& (35)\end{array}$

The model equations are Equations (12)–(20) with Equations (27), (32), (34), (35) and $i_{G1}^{Gem}=i_S^{Gem}=1$, $i_{G2}^{Gem}=0$, $i_{G1}^{Tbd}=i_S^{Tbd}=i_{G2}^{Tbd}=1$.

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:

$\begin{array}{ll}{V}_i={\left({\sigma }_1+{\sigma }_2{Y}_i\right)}^2& (36)\end{array}$

where V_i is the variance of the ith time point, σ₁ and σ₂ 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₀/G₁ phase with time (Figures 2A,C). Numerous factors can cause accumulation of the cell population in the G₀/G₁ 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₀/G₁ phase (Hayes et al., 2005; Choresca et al., 2009; Dalman et al., 2010).

FIGURE 2

Figure 2. Cell cycle phase distributions in the vehicle control groups for MiaPaCa-2 and BxPC-3 cells. (A,C) Percent of cells in G₀/G₁, S, and G₂/M phases at 0, 24, 48, 72, and 96 h (MiaPaCa-2) and 0, 16, 26, 48, 72, and 96 h (BxPC-3) for vehicle control groups. Columns, mean of triplicate; bars, SD. N = 3. (B,D) Data and model fitting of vehicle control groups. Numbers of cells in G₀/G₁, S, and G₂/M phases, as well as total cells over time, are represented by symbols. Model fittings based on the model of Figure 1B are indicated by lines.

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⁶ cells (fixed to a value observed in the study, data not shown) for MiaPaCa-2 and 7.98 × 10⁵ 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 2

Table 2. Pharmacodynamic model parameter estimates for two cell lines. Parameters were obtained after fitting to the cell cycle base model shown in Figure 1B.

Gemcitabine Effects upon Cell Cycle and Model Prediction

Figure 3 shows the effects of gemcitabine on the cell cycle distribution of MiaPaCa-2 (Figures 3A–C) and BxPC-3 cells (Figures 3D–F). As time increased, 45 nM gemcitabine induced S phase accumulation of MiaPaCa-2 cells (Figure 3A and Supplemental 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 concentration-dependent 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.

FIGURE 3

Figure 3. Gemcitabine effects on the cell cycle for MiaPaCa-2 and BxPC-3 cells. (A,D) Histograms showing temporal changes in cellular DNA content after incubation with 45 nM gemcitabine in MiaPaCa-2 cells and 34 nM in BxPC-3 cells. All data are in triplicate. (B,E) fraction of cells in subG₁, G₀/G₁, S, and G₂/M phases after incubating cells with concentrations of gemcitabine shown for 48 (MiaPaCa-2) or 76 (BxPC-3) h. Columns, mean of triplicate determinations; bars, SD. (C,F) Data and model fitting results of cells incubating cells with concentrations of gemcitabine shown for up to 120 h. Symbols represent data for cells numbers in G₀/G₁, S, G₂/M, subG₁ phase as well as total cell numbers; lines show model fittings based on PD models in Figure 1C (MiaPaCa-2) and Figure 1D (BxPC-3).

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⁻¹ for MiaPaCa-2 and 1 nM⁻¹ for BxPC-3 cells. The model assumed that gemcitabine-induced apoptosis occurred in S phase in MiaPaCa-2 cells, and in both G₀/G₁ 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₁ increased with rate of α_k1 (0.0612 nM⁻¹) with increased gemcitabine concentrations. For BxPC-3, it was assumed that k₂ and k₃ increased with rates of α_k2 (0.0635 h⁻¹) and α_k3(0.0605 nM⁻¹) as gemcitabine concentrations increased.

TABLE 3

Table 3. Pharmacodynamic model parameter estimates. Parameters were obtained after fitting to cell cycle models for either gemcitabine (Figures 1C,D) or trabectedin (Figure 1E) as single agents.

Trabectedin Effects upon Cell Cycle and Model Prediction

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₂/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 concentration-dependent (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₂/M phase.

FIGURE 4

Figure 4. Trabectedin effects on the cell cycle for MiaPaCa-2 and BxPC-3 cells. (A,D) histograms showing temporal changes in cellular DNA content after continuous exposure of cells to 0.8 nM (MiaPaCa-2) or 1.1 nM (BxPC-3) trabectedin. All data are in triplicate. (B,E) Fraction of cells in sub G₁, G₀/G₁, S, and G₂/M phase after incubating cells in the concentrations of trabectedin shown for 24, 48, and 72 h (MiaPaCa-2) or 26, 48, and 76 h (BxPC-3). Columns, mean of triplicate determinations; bars, SD. (C,F) Temporal change in number of cells in G₀/G₁, S, G₂/M, sub G₁ phase, as well as total cell numbers, after incubating cells with concentrations of trabectedin shown for up to 120 (MiaPaCa-2) and 96 h (BxPC-3). Symbols represent observed data and lines represent fitted curves generated from PD model shown in Figure 1E.

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₂/M phase transition k₂ and the G₂/M to G₀/G₁ phase transition k₃. In order to characterize early S phase arrest followed by subsequent G₂/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₂/M phase increased with time with rate k_G2. The k_S was 1 h⁻¹ for MiaPaCa-2 and 0.00204 h⁻¹ for BxPC-3, whereas k_G2 was 0.0270 h⁻¹ for MiaPaCa-2 and 0.0516 h⁻¹ 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₂/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₂/M phase arrest than BxPC-3 cells. In addition, we also assumed concentration-dependent trabectedin effects on k₁, with rates of β_k1 of 2.64 nM⁻¹ for MiaPaCa-2 and 0.637 nM⁻¹ for BxPC-3. For BxPC-3, k₃ also increased with trabectedin concentration, with a rate of β_k3 of 0.296 nM⁻¹. 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₂/M phase cells (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.

FIGURE 5

Figure 5. Drug combination effects and model-fittings based on individual drug mechanisms for MiaPaCa-2 and BxPC-3 cells. (A,B) Observed data and model fitting results for cells exposed to different combinations of drug concentrations. Symbols represent cell numbers in G₀/G₁, S, G₂/M, and sub G₁ phases; lines depict model fittings to the PD model shown in Figures 1F,G. (C) Shown are temporal cell proliferation profiles and model fitting results for single drug effects and four different drug combinations for total BxPC-3 cells. Symbols represent experimental data and lines indicate the corresponding model fittings.

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 ψ₁ was multiplied at IC_{50, Gem}. The ψ₁ interaction parameter was fixed to 1 for MiaPaCa-2 cells (because modeling of ψ₁ 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 ψ₂ was multiplied by KC_{50, Gem} to characterize drug interactions for induction of apoptosis. Estimates of ψ₂ were 0.499 for MiaPaCa-2 cells and 0.363 for BxPC-3 cells, indicating synergistic interactions in inducing apoptosis.

TABLE 4

Table 4. Pharmacodynamic interaction model parameter estimates.

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₀/G₁, S, G₂/M) and in apoptosis (sub G₁). 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₀/G₁ 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₁ phase. If cells are able to progress to S phase, normal cycling occurs. Otherwise, cells enter a resting G₀ 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₀/G₁ 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₁ 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. (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₀/G₁ and an inhibition of transition rates from S phase to G₂/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₂/M phase. In order to capture the time-dependency, equations were constructed such that an inhibition term on the transition from S to G₂/M decreased with time, and the inhibition term on the transition from G₂/M to G₀/G₁ 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₀/G₁ phases. The models were constructed with differences between cell lines. For example, based on the KC₅₀ of the two cell lines (Miao et al., 2016a), MiaPaCa-2 is more trabectedin-sensitive while BxPC-3 is more gemcitabine-sensitive. 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, ψ₁ was applied to represent drug interactions at the point of S phase arrest, and ψ₂ was applied to represent drug interactions at the point of induction of apoptosis. However, ψ₁ was approximately equal to 1, suggesting that drug interactions in inducing S phase accumulation were additive, whereas ψ₂ 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₂/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 time-dependent. 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 multi-stage 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.

Author Contributions

Participated in research design: XM and WJJ. Conducted experiments: XM. Performed data analysis: XM and GK. Wrote or contributed to the writing of the manuscript: XM, GK, SAO, RMS, WJJ.

Funding

This work was supported by National Institutes of Health Grants 24211 to WJJ and CA168454 and CA198096 to RMS, the National Research Fund, Luxembourg to GK, and co-funded under the Marie Curie Actions of the European Commission (FP7-COFUND).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors would like to acknowledge the assistance of Nancy Pyszczynski for facilitating data analysis using the FACSCalibur flow cytometer. The authors also thank Dr. Paul Wallace and Joseph D. Tario of Roswell Park Cancer Institute for their support and help with the project. Trabectedin was a gift from PharmaMar.

Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fphar.2016.00421/full#supplementary-material

Abbreviations

CCM, cell cycle model; DDI, drug-drug interaction; PK/PD, pharmacokinetics/pharmacodynamics.

References

Ariens, E. J., Van Rossum, J. M., and Simonis, A. M. (1957). Affinity, intrinsic activity and drug interactions. Pharmacol. Rev. 9, 218–236.

PubMed Abstract | Google Scholar

Basse, B., Baguley, B. C., Marshall, E. S., Joseph, W. R., van Brunt, B., Wake, G., et al. (2004). Modeling cell death in human tumor cell lines exposed to the anticancer drug paclitaxel. J. Math. Biol. 49, 329–357. doi: 10.1007/s00285-003-0254-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Chakraborty, A., and Jusko, W. J. (2002). Pharmacodynamic interaction of recombinant human interleukin-10 and prednisolone using in vitro whole blood lymphocyte proliferation. J. Pharm. Sci. 91, 1334–1342. doi: 10.1002/jps.3000

PubMed Abstract | CrossRef Full Text | Google Scholar

Choresca, C. H. Jr., Koo, O. J., Oh, H. J., Hong, S. G., Gomez, D. K., Kim, J. H., et al. (2009). Different culture conditions used for arresting the G0/G1 phase of the cell cycle in goldfish (Carassius auratus) caudal fin-derived fibroblasts. Cell Biol. Int. 33, 65–70. doi: 10.1016/j.cellbi.2008.09.015

PubMed Abstract | CrossRef Full Text | Google Scholar

D'Argenio, D. Z., Schumitzky, A., and Wang, X. (2009). ADAPT 5 User's Guide: Pharmacokinetic/Pharmacodynamics Systems Analysis Software. Los Angeles, CA: Biomedical Simulations Resource.

D'Incalci, M., Erba, E., Damia, G., Galliera, E., Carrassa, L., Marchini, S., et al. (2002). Unique features of the mode of action of ET-743. Oncologist 7, 210–216. doi: 10.1634/theoncologist.7-3-210

PubMed Abstract | CrossRef Full Text

D'Incalci, M., and Galmarini, C. M. (2010). A review of trabectedin (ET-743): a unique mechanism of action. Mol. Cancer Ther. 9, 2157–2163. doi: 10.1158/1535-7163.MCT-10-0263

PubMed Abstract | CrossRef Full Text | Google Scholar

D'Incalci, M., and Zambelli, A. (2015). Trabectedin for the treatment of breast cancer. Expert Opin. Investig. Drugs. 25, 105–115. doi: 10.1517/13543784.2016.1124086

PubMed Abstract | CrossRef Full Text | Google Scholar

Dalman, A., Eftekhari-Yazdi, P., Valojerdi, M. R., Shahverdi, A., Gourabi, H., Janzamin, E., et al. (2010). Synchronizing cell cycle of goat fibroblasts by serum starvation causes apoptosis. Reprod. Domest. Anim. 45, e46–e53. doi: 10.1111/j.1439-0531.2009.01520.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Deer, E. L., González-Hernández, J., Coursen, J. D., Shea, J. E., Ngatia, J., Scaife, C. L., et al. (2010). Phenotype and genotype of pancreatic cancer cell lines. Pancreas 39, 425–435. doi: 10.1097/MPA.0b013e3181c15963

PubMed Abstract | CrossRef Full Text | Google Scholar

Evan, G. I., and Vousden, K. H. (2001). Proliferation, cell cycle and apoptosis in cancer. Nature 411, 342–348. doi: 10.1038/35077213

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferrell, J. E. Jr., Tsai, T. Y., and Yang, Q. (2011). Modeling the cell cycle: why do certain circuits oscillate? Cell 144, 874–885. doi: 10.1016/j.cell.2011.03.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Florian, J. A. Jr., Eiseman, J. L., and Parker, R. S. (2005). Accounting for quiescent cells in tumour growth and cancer treatment. Syst. Biol. (Stevenage). 152, 185–192. doi: 10.1049/ip-syb:20050041

PubMed Abstract | CrossRef Full Text | Google Scholar

Gajate, C., An, F., and Mollinedo, F. (2002). Differential cytostatic and apoptotic effects of ecteinascidin-743 in cancer cells. Transcription-dependent cell cycle arrest and transcription-independent JNK and mitochondrial mediated apoptosis. J. Biol. Chem. 277, 41850–41859. doi: 10.1074/jbc.M204644200

PubMed Abstract | CrossRef Full Text | Google Scholar

Gardner, S. N. (2002). Modeling multi-drug chemotherapy: tailoring treatment to individuals. J. Theor. Biol. 214, 181–207. doi: 10.1006/jtbi.2001.2459

PubMed Abstract | CrossRef Full Text | Google Scholar

Hamed, S. S., and Roth, C. M. (2011). Mathematical modeling to distinguish cell cycle arrest and cell killing in chemotherapeutic concentration response curves. J. Pharmacokinet. Pharmacodyn. 38, 385–403. doi: 10.1007/s10928-011-9199-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Hamed, S. S., Straubinger, R. M., and Jusko, W. J. (2013). Pharmacodynamic modeling of cell cycle and apoptotic effects of gemcitabine on pancreatic adenocarcinoma cells. Cancer Chemother. Pharmacol. 72, 553–563. doi: 10.1007/s00280-013-2226-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Hayes, O., Ramos, B., Rodriguez, L. L., Aguilar, A., Badia, T., and Castro, F. O. (2005). Cell confluency is as efficient as serum starvation for inducing arrest in the G0/G1 phase of the cell cycle in granulosa and fibroblast cells of cattle. Anim. Reprod. Sci. 87, 181–192. doi: 10.1016/j.anireprosci.2004.11.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Jusko, W. J. (1973). A pharmacodynamic model for cell-cycle-specific chemotherapeutic agents. J Pharmacokin Biopharm. 1, 175–200. doi: 10.1007/BF01062346

CrossRef Full Text | Google Scholar

Koch, G., Schropp, J., and Jusko, W. J. (2016). Assessment of non-linear combination effect terms for drug-drug interaction. J. Pharmacokinet. Pharmacodyn. 43, 461–479. doi: 10.1007/s10928-016-9490-0

CrossRef Full Text | Google Scholar

Koch, G., Walz, A., Lahu, G., and Schropp, J. (2009). Modeling of tumor growth and anticancer effects of combination therapy. J. Pharmacokinet. Pharmacodyn. 36, 179–197. doi: 10.1007/s10928-009-9117-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Kozusko, F., Chen, P., Grant, S. G., Day, B. W., and Panetta, J. C. (2001). A mathematical model of in vitro cancer cell growth and treatment with the antimitotic agent curacin A. Math. Biosci. 170, 1–16. doi: 10.1016/S0025-5564(00)00065-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Miao, X., Koch, G., Straubinger, R. M., and Jusko, W. J. (2016a). Pharmacodynamic modeling of combined chemotherapeutic effects predicts synergistic activity of gemcitabine and trabectedin in pancreatic cancer cells. Cancer Chemother. Pharmacol. 77, 181–193. doi: 10.1007/s00280-015-2907-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Miao, X., Shen, S., Li, J., Shen, X., Wang, X., Koch, G., et al. (2016b). “Mathematical models of temporal proteomic profiles of pancreatic cancer cells in response to gemcitabine and trabectedin combinations,” in Poster Presentation at the 2016 AAPS Annual Meeting and Exposition (Denver, CO), Poster 08R1030.

Morgan, M. A., Parsels, L. A., Parsels, J. D., Mesiwala, A. K., Maybaum, J., and Lawrence, T. S. (2005). Role of checkpoint kinase 1 in preventing premature mitosis in response to gemcitabine. Cancer Res. 65, 6835–6842. doi: 10.1158/0008-5472.CAN-04-2246

PubMed Abstract | CrossRef Full Text | Google Scholar

Oberstein, P. E., and Olive, K. P. (2013). Pancreatic cancer: why it is so hard to treat? Therap. Adv. Gastroenterol. 6, 321–337. doi: 10.1177/1756283X13478680

PubMed Abstract | CrossRef Full Text | Google Scholar

Panetta, J. C., Evans, W. E., and Cheok, M. H. (2006). Mechanistic mathematical modeling of mercaptopurine effects on cell cycle of human acute lymphoblastic leukaemia cells. Br. J. Cancer. 94, 93–100. doi: 10.1038/sj.bjc.6602893

CrossRef Full Text | Google Scholar

Ribba, B., You, B., Tod, M., Girard, P., Tranchand, B., Trillet-Lenoir, V., et al. (2009). Chemotherapy may be delivered based on an integrated view of tumor dynamics. IET Syst. Biol. 3, 180–190. doi: 10.1049/iet-syb.2008.0104

CrossRef Full Text | Google Scholar

Robinson, H. M., Jones, R., Walker, M., Zachos, G., Brown, R., Cassidy, J., et al. (2006). Chk1-dependent slowing of S-phase progression protects DT40 B-lymphoma cells against killing by the nucleoside analogue 5-fluorouracil. Oncogene 25, 5359–5369. doi: 10.1038/sj.onc.1209532

PubMed Abstract | CrossRef Full Text | Google Scholar

Senderowicz, A. M. (2004). Targeting cell cycle and apoptosis for the treatment of human malignancies. Curr. Opin. Cell Biol. 16, 670–678. doi: 10.1016/j.ceb.2004.09.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Sherer, E., Hannemann, R. E., Rundell, A., and Ramkrishna, D. (2006). Analysis of resonance chemotherapy in leukemia treatment via multi-staged population balance models. J. Theor. Biol. 240, 6468–6661. doi: 10.1016/j.jtbi.2005.11.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Simoens, C., Korst, A. E., De Pooter, C. M., Lambrechts, H. A., Pattyn, G. G., Faircloth, G. T., et al. (2003). In vitro interaction between ecteinascidin 743 (ET-743) and radiation, in relation to its cell cycle effects. Br. J. Cancer. 89, 2305–2311. doi: 10.1038/sj.bjc.6601431

PubMed Abstract | CrossRef Full Text | Google Scholar

Yip-Schneider, M. T., Sweeney, C. J., Jung, S. H., Crowell, P. L., and Marshall, M. S. (2001). Cell cycle effects of nonsteroidal anti-inflammatory drugs and enhanced growth inhibition in combination with gemcitabine in pancreatic carcinoma cells. J. Pharmacol. Exp. Ther. 298, 976–985.

PubMed Abstract | Google Scholar

Zhu, X., Straubinger, R. M., and Jusko, W. J. (2015). Mechanism-based mathematical modeling of combined gemcitabine and birinapant in pancreatic cancer cells. J. Pharmacokinet. Pharmacodyn. 42, 477–496. doi: 10.1007/s10928-015-9429-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Appendix

Behavior of Cell Cycle Model with Saturation for Increasing Time

For t → ∞ we obtain from Equations (6)–(8):

$\begin{array}{ll}0=2k_3{G}_2^{*}\left(1-\frac{G_1^{*}+S^{*}+G_2^{*}}{2R_{ss}}\right)-k_1\left(1-I_{max}\right)G_1^{*}& (\mathrm{A1})\end{array}$

$\begin{array}{ll}0=k_1\left(1-I_{max}\right)G_1^{*}-k_2{S}^{*}& (\mathrm{A2})\end{array}$

$\begin{array}{ll}0=k_2{S}^{*}-k_3{G}_2^{*}.& (\mathrm{A3})\end{array}$

From Equation (A3) one finds:

$\begin{array}{ll}{S}^{*}=\frac{k_3}{k_2}{G}_2^{*}.& (\mathrm{A4})\end{array}$

From Equations (A2) and (A4) we have:

$\begin{array}{ll}{G}_2^{*}=\frac{k_1\left(1-I_{max}\right)}{k_3}{G}_1^{*}.& (\mathrm{A5})\end{array}$

Hence, with Equation (A5) we obtain:

$\begin{array}{ll}{S}^{*}=\frac{k_1}{k_2}\left(1-I_{max}\right)G_1^{*}.& (\mathrm{A6})\end{array}$

Substituting Equations (A5) and (A6) in (A1) results with b = 1−I_max in:

$\begin{array}{l}0=2k_{1}b{G}_1^{*}\left(1-\frac{G_1^{*}\left(1+\frac{k_1}{k_2}b+\frac{k_1}{k_3}b\right)}{2R_{ss}}\right)-k_{1}b{G}_1^{*}.\end{array}$

We set $a=1+\frac{k_1}{k_2}b+\frac{k_1}{k_3}b$, and cancel out the zero-solution to obtain:

$\begin{array}{ll}0=1-\frac{G_1^{*}a}{R_{ss}}\iff G_1^{*}=\frac{R_{ss}}{a}& (\mathrm{A7})\end{array}$

yielding Equations (10) and (11). Summation of $G_1^{*}$, S* and $G_2^{*}$ from Equations (10) and (11) fulfills $R_{tot}^{*}=R_{ss}^{*}$.