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.

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).

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.

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.

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.

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₁ 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:

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:

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:

The Appendix provides derivations of Equations (10) and (11).

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:

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:

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.

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).

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.

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.

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, σ₁ 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 %).

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).

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.

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.

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.

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.

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).

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.

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.

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.