Integrated experimental and simulation study of the response to sequential treatment with erlotinib and gemcitabine in pancreatic cancer

The combination of erlotinib with gemcitabine is one of the most promising therapies for advanced pancreatic cancer. Aiming at optimizing this combination, we analyzed in detail the response to sequential treatments with erlotinib → gemcitabine and gemcitabine → erlotinib with an 18 h interval, adopting a previously established experimental/computational approach to quantify the cytostatic and cytotoxic effects at G1, S and G2M checkpoints. This assessment was achieved by contemporary fits of flow cytometric and time-lapse experiments in two human pancreatic cancer cell lines (BxPC-3 and Capan-1) with a mathematical model reproducing the fluxes of cells through the cycle during and after treatment. The S-phase checkpoint contributes in the response to erlotinib, suggesting that the G1 arrest may hamper S-phase cytotoxicity. The response to gemcitabine was driven by the dynamics of the progressive resumption from the S-phase arrest after drug washout. The effects induced by single drugs were used to simulate combined treatments, introducing changes when required. Gemcitabine → erlotinib was more than additive in both cell lines, strengthening the cytostatic effects on cells recovering from the arrest induced by gemcitabine. The interval in the erlotinib → gemcitabine sequence enabled to overcome the antagonist effect of G1 block on gemcitabine efficacy and improved the outcome in Capan-1 cells.


Computer simulation of cell cycle and drug effects
A thorough description and discussion of the approach has been published before [1][2][3][4][5], here we will sketch the outlines of the modeling.
First we modeled the basic cell proliferation process in untreated cells (controls) and then we modeled proliferation during and after treatment applying perturbations to basic cell cycling to reproduce the data of time course experiments with different drug concentrations and schedules.

Basic cell cycle model (untreated cells)
Wishing to consider not only the overall cell number, as in traditional pharmacodynamic models, but also the amount of information conveyed by the time courses of flow cytometry and time-lapse experiments, we introduced a proper cell cycle and age structure. The basic dynamics of progression inside each phase is described by continuity equations, in discrete time intervals Δ (from 0, i.e. the start of the experiment, to t end at the end of the experiment) on the age distributions of cells in each phase. The "state vectors" N ph (a i ,t) represent the number of cells at age a i in a phase ph (G 1 , S or G 2 M) at time t. Each phase is divided in k ph compartments (k ph Δ = T ph , where T ph is the duration of phase ph) grouping cells with ages respectively from 0 to Δ (a 1 ), from Δ to 2Δ (a 2 )… from (k ph -1)Δ to k ph Δ (a kph ). Model parameters are the durations of cell cycle phases (T G1 , T S , T G2M ). Thus the time evolution of N ph (a,t) within a phase is given by: while cells entering a phase are those at the end of the previous phase at the previous time step: where the factor 2 in the first equation accounts for the fact that two daughter cells are produced for each mitotic cell ending the previous cycle.
The overall number of cells is obviously N(t) = Σ ph Σ i N ph (a i ,t) while typical flow cytometric data like percentages of cells in cell cycle phases (%G 1 , %S and %G 2 M) can be easily derived, e.g. %G 1 (t)= 100 x Σ aG1 N G1 (a G1 ,t)/N(t).
When dealing with time-lapse data, a further generation structure can be trivially included, repeating the previous equations in each cell generation (gen h , starting from gen 0 of the cells at the start of the experiment): N ph,genh (a i + 1 ,t + Δ) = N ph,genh (a i , t) N G1,genh (a 1 , t + Δ) = 2 N G2M,gen(h-1) (a kG2M , t) N S,genh (a 1 , t + Δ) = N G1,genh (a kG1 , t) N G2M,genh (a 1 , t + Δ) = N S,genh (a kS , t) This structure enables to derive the overall cell number in each generation, and compare modeling prediction with the corresponding experimental results obtained by cell tracking in time-lapse experiments.
However this first basic model with fixed phase durations is not realistic and would be inconsistent with the direct time-lapse measure of intermitotic times, demonstrating a wide distribution of values among individual cells also in a homogenous environment, as it occurs in culture flasks in in vitro studies. In the cell lines used in the present work, intermitotic times varied from 12 h to 48 h (Supplementary Figure 1) with typical rightskewed distributions, well described by a log-normal or, better, by a reciprocal-normal function [5,6]. Thus, at the second level of complexity of the model, intercell variability was included, assuming that the duration of each phase is distributed according to a reciprocal-normal (F ph (a)), specified by its mean ( T ph ) and standard deviation or coefficient of variation (CV ph ). The distribution of intermitotic times F(Tc) is calculated as where the summation includes all combinations of a G1 ,a S , aG 2M such that a G1 + a S + a G2M = Tc F ph enables to calculate the probability that cells having reached the age a i will complete and exit the phase ph in the following time step as: Integrated experimental and simulation study of the response to sequential treatment with erlotinib and gemcitabine in pancreatic cancer

Supplementary Materials
with the following balance equations, connecting the cell cycle distribution at time t with that at t+Δ: Intercell variability (if constant) drives the cell population towards a unique asymptotic asynchronous steady state of growth, where the overall number of cells increases exponentially, but the percentage of cells in each age and phase is constant [7]. In the practice, the steady state of growth represents a reasonable approximation of the common exponential phase of growth in cell cultures (but also in vivo, as far as environmental conditions remain constant), where cell cycle distributions are timeindependent. By starting an experiment when the cell population is in this steady state, as it is usually done, the initial cell distribution was determined by a preliminary desynchronization run for each set of parameters' values.
In this way the models of steady state growth of BxPC-3 and Capan-1 were derived with the following parameters values: TG1 8.6 h; T S 11.0 h; TG2M 4.7 h; CV G1 85%; CV S and CV G2M 30% for BxPC-3 and T G1 10.3 h; T S 11.0 h; T G2M 5.0 h; CV G1 , CV S and CV G2M 30% for Capan-1. The computer program in use would enable to build models simulating cell cycling with higher complexity, e.g. introducing quiescent cells and loss and simulating a progressive approach to confluence of the cell population [5], but these effects were negligible and the study of these details was out of the aim of the present work.

Modeling response to treatment
The effects of treatment were described by parameters measuring delays, blocks and killing of cells in each phase. This is reflected in the structure of the modeling framework in silico, including modules for G 1 , S and G 2 M checkpoint responses in generations 0, 1, 2, etc., superimposing and modifying the flow of the cells through the cycle.
Modules render complex biological phenomena, such as the operation of a checkpoint in a specific cell cycle phase, with the minimal choice of parameters so that the main antiproliferative effects (block, block's recovery or death) occurring in that phase can be quantified, as previously described with full mathematical details [5].
Briefly, to test drug or radiation effects in a particular phase, the researcher can choose among different types of modules, rendering checkpoint activity with increasing levels of complexity. At the lowest level (type I), the effect is simply a delay of the progression in the phase where it is located, lengthening the phase. The balance equations inside a phase were modified by a "delay" parameter (Del) as follows: N ph,genh (a i + 1 ,t + Δ) = (1-Del ph ) (1-β ph (a i )) N ph,genh (a i ,t) + Del ph (1-β ph (a i + 1 )) N ph,genh (a i + 1 ,t) meaning that a fraction (1-Del ph ) of cells progresses in the cell cycle while a fraction Del ph remains in the same age compartment at each time step. Analogous straightforward modifications apply to other equations.
The associated "delay" parameter in S phase is equivalent to the fractional reduction of the average DNA replication rate, providing a quantitative measure of the inhibition of DNA synthesis, net of other confounding factors (like the overall number of cells or the number of dying cells) concurring in measures obtained with other tests. Cell death was regulated by a second, independent, "death rate" parameter (DR), simply representing the fraction of cells in the phase that die at each time step. Thus equations become: being (1-DR ph ) the fraction of cells surviving, while the fraction DR ph is lost.
A type II module acts by arresting permanently a fraction of the cells transiting to the next phase in a specific compartment of blocked cells ("block probability" parameter, pBL ph ). Cell death is measured by a second, independent parameter -"death rate in block" (DRBL ph ) -representing the fraction of cells in the blocked-cell compartment that die at each time step. We used type I or II modules to emulate G 1 and G 2 M checkpoint activities, while only type I was considered appropriate to describe the activity of S-phase checkpoint, resulting in a delay or complete inhibition of DNA replication and not in an arrest at a specific point within the S phase as assumed by type II. New equations introducing the compartments of blocked cells (B G1 and B G2M ) were added, while corresponding changes were applied to the first compartment in the subsequent phases, reached only by the fraction of non-blocked cells: Type III aims at rendering a possible resumption of blocked cells in cycle, adding a new parameter describing recycling -"recycling rate" (Rec ph ), the fraction of cells in the blocked-cell compartment that resume cycling, entering the next phase -to the type II module. Further complexity can be introduced as desired, including, in addition to phase and generation dependence, a timedependence for the onset or disappearance of one or more of the effects (block, recycling or death). The previous equations are accordingly modified: Notice that parameters of each module are probabilities (e.g. "block probability") or rates (e.g. "death rate") of a specific event in the phase and generation where the checkpoint is located, used in a first order approximation to determine the fraction of cells blocked or dead within the specific cohorts of cells arriving at the checkpoints. This enabled us to introduce the heterogeneity of the response to treatment in the model, avoiding overfitting by considering the rich data set of the time-course of flow cytometric and time-lapse data.

Building erlotinib and gemcitabine models
In a preliminary screening of biologically sound models we initially tested models including only G 1 delay in the case of erlotinib and S delay in the case of gemcitabine, optimizing parameter values, but these models were unable to predict the experimental time courses. The final models were reached by progressively increasing the number of parameters until simulation of the time course of cell cycling satisfactorily described the data for all drug concentrations and experimental platforms, avoiding over-fitting by the use of the likelihood ratio test. The final gemcitabine model included the following modules: 1) S-phase module type I in generations 0 and 1, with time dependence of the Del S parameter. Maximum reduction of the DNA synthesis rate (dependent on drug concentration) was reached immediately (i.e. within 0.5 h) when the drug was added and persisted several hours after drug removal. A sharply decreasing Hill function (sigmoidicity -10), with half-time (time to reach half the maximum value, or half the difference between maximum and minimum when minimum is not zero) of 18 h (generation 0) or 21 h (generation 1) was suitable in all treatment groups of both cell lines to describe the timedependence of the S-phase delay parameter. At higher gemcitabine concentrations data were incompatible with complete recovery of DNA synthesis and a minimum long-term delay was maintained through generations 0 and 1. S-phase death rate was observable in BxPC-3 only after 48 h (generation 0) or 72 h (generation 1) and was described by an increasing Hill function (sigmoidicity 10) with half-time 48 h or 72 h.
2) G 1 phase module type I in generation 0, with time dependence of Del G1 parameter in generation 0. Timedependence was described by a Hill function (sigmoidicity -10) with half-time depending on drug treatment (BxPC-3: 6 h, 15 h and 27 h for 20, 40 and 120 nM gemcitabine; Capan-1: 15 h for 100 nM; at 30 nM no G 1 delay was detected).
3) G 2 M phase module type II in generation 0, with time dependence of pBl G2M parameter. Time-dependence was described by an increasing Hill function (sigmoidicity 10), initially zero, with half-maximum at 12 h (BxPC-3) or 6 h (Capan-1), meaning that this effect did not involve cells initially in G 2 M (which all exit before that time) but only cells reaching G 2 M after recovery of the DNA synthesis rate in the previous S phase.
Thus seven variable parameters were optimized in each gemcitabine treatment group: maximum S delay (generations 0 and 1), long-term S delay in generation 0, long-term S delay in generation 1, S-phase death rate in generation 0, S-phase death rate in generation 1, maximum G 1 delay (generation 0), maximum G 2 M delay (generation 0).
In the erlotinb model we separated the response in the periods during (0-48 h) and after (> 48 h) treatment, applying distinct modules in the two intervals as follows: 1) G 1 phase module type II in generation 0, with time dependence of pBl G1 during treatment. The onset of the effect was not immediate after drug addition: timedependence was described by an increasing Hill function (sigmoidicity 1), initially zero, with half-maximum at 6 h.
2) G 1 phase module type II in generation 1 during treatment. The model cannot detect short-term time dependence for the onset of this effect, as a negligible amount of cells has time to complete G 1 in generation 1 in the first 6 h. Thus pBl G1 was kept time-independent.
3) S-phase module type I in generation 0, with time dependence of Del S (sigmoidicity 1, half-maximum at 6 h) during treatment. The S-phase death rate was assumed constant and equal in generations 0 and 1, as the data were insensitive to more detailed, time-dependent models. 4) S-phase module type I in generations 1 and 2 during treatment. Because the data were not sensitive to different values in these generations, a unique Del S parameter was used for both. 5) G 2 M phase module type II in generation 0, with time dependence of pBl G2M (sigmoidicity 1, half-maximum at 6 h) during treatment. 6) G 1 phase module type III in generations 0 and 1 after treatment (mainly). Although no more cells were blocked in G 1 after treatment, previously blocked cells exited the block with a recycling rate following an increasing Hill function with half-maximum at 48 h and sigmoidicity 10, meaning that the phenomenon actually started before the end of treatment. A unique pBl G1 parameter was used for both generations.
7) S-phase module type I in generations 1 and 2 after treatment. Because the data were not sensitive to different values in these generations, a unique Del S parameter was used for both. 8) G 2 M phase module type I in generations 1 and 2 after treatment. Because the data were not sensitive to different values in these generations, a unique pBl G2M parameter was used for both.
The following nine parameters were optimized in each erlotinib concentration, during treatment: maximum G 1 block probability in generation 0, G 1 block probability in generation 1, maximum S delay in generation 0, S-phase death rate in generations 0 and 1, S delay in generations 1 and 2, maximum G 2 M block probability in generation 0; and after treatment: maximum G 1 recycling rate in generations 0 and 1, S delay in generations 1 and 2, G 2 M delay in generations 1 and 2.
In Capan-1, with the concentrations used in this study, the data were explained with a subset of six parameters, setting the following parameters to zero: S-phase death rate in generations 0 and 1, S delay in generations 1 and 2 and G 2 M delay in generation 1 and 2.

Optimization and uncertainty analysis
Models were optimized by non-linear fitting of data for individual doses, simulating a complete time course with 0.5 h time step at each tentative set of parameters' values. The objective function to be minimized during the optimization was the overall negative log-likelihood (-log(L)), obtained by summing the contributions of all experimental procedures and assuming normally distributed errors with variance typical of each procedure [5].
We estimated the uncertainty of the best-fit parameters, calculating likelihood-based confidence intervals, giving the range for each parameter within which the fit remained not significantly worse than that obtained with the best-fit parameters, according to a likelihood ratio test. We took into account every parameter in the best scenario reproducing the experimental data and changed them, one by one. Likelihood-based 95% confidence intervals for each parameter were obtained by raising or lowering its value until log(L) was reduced from its bestfit value log(Lbest) to log(L) = log(Lbest) -χ 2 0.05,1 /2 [8].