The antiproliferative response to anticancer treatment is the result of concurrent responses in all cell cycle phases, extending over several cell generations, whose complexity is not captured by current methods. In the proposed experimental/computational approach, the contemporary use of time-lapse live cell microscopy and flow cytometric data supported the computer rendering of the proliferative process through the cell cycle and subsequent generations during/after treatment. The effects of treatments were modelled with modules describing the functional activity of the main pathways causing arrest, repair and cell death in each phase. A framework modelling environment was created, enabling us to apply different types of modules in each phase and test models at the complexity level justified by the available data. We challenged the method with time-course measures taken in parallel with flow cytometry and time-lapse live cell microscopy in X-ray-treated human ovarian cancer cells, spanning a wide range of doses. The most suitable model of the treatment, including the dose-response of each effect, was progressively built, combining modules with a rational strategy and fitting simultaneously all data of different doses and platforms. The final model gave for the first time the complete rendering in silico of the cycling process following X-ray exposure, providing separate and quantitative measures of the dose-dependence of G1, S and G2M checkpoint activities in subsequent generations, reconciling known effects of ionizing radiations and new insights in a unique scenario.
Anticancer research spans a wide range of scales, from the microscopic/molecular up to the macroscopic level of clinical assessment of treatment efficacy. On an intermediate scale of preclinical testing in vitro and in vivo, cancer research maintains the need to evaluate the antiproliferative activity of old and new treatment options in cancer cell populations. In fact, even for new drugs, assessment of activity against a molecular target does not guarantee antiproliferative activity and control of expansion – if not eradication – of cancer cells, which is clearly the primary goal of any treatment. However, the impressive advances in the molecular understanding of cancer (s) have not been accompanied by a deeper comprehension of proliferation in treated samples, often simplistically estimated with scores based on treated/control ratios. These do not catch the complexity of the phenomena in play, leaving unanswered several questions, from the contribution of inter-cell heterogeneity to the origin of dose-dependence, from the dynamics of the response of checkpoint controls to the overlapping of cytostatic and cytotoxic phenomena [1]–[3]. In particular, biomedical research is now acknowledging the importance of cell-to-cell diversity and dynamic variability for protein levels [4]–[7] and response to environmental interactions [8] or treatments [9], even within populations of genetically identical cells. This calls for a probabilistic view of the single cell response representing a change of perspective compared to exclusively molecularly-driven research, often dominated by a deterministic-mechanistic view according to which a drug makes a given effect or not, without specifying the probabilities of the events and their dose-dependence. Mathematical modelling is increasingly adopted to tackle biological complexity, together with experimental procedures producing appropriate quantitative data, for in silico rendering of biological structures and processes in different fields and scales, from X-ray crystallography to medical imaging [10]–. The question is normally tackled by adopting a computational model of the biological phenomenon, whose inputs are meaningful biological parameters and outputs are measurable quantities. For instance, in the crystallography field a model of the diffraction keeps the 3D structure of a molecule as input and gives as output the data that a molecule' s crystal would produce when challenged in X-ray diffraction experiments. The model can be used in two ways: to infer the 3D structure from experimental data (optimization problem) or to simulate the expected data from hypothetical 3D structures (simulation) [15]. Adopting a conceptually similar approach, we present here a mixed experimental/computational method (Figure 1) to render the process of proliferation at the cell population level, using a computational model whose input parameters are simple descriptors of the functional activities of the main intracellular molecular controls of the cell cycle and whose outputs can be directly fitted to data obtained by time-lapse live cell microscopy (TL) and DNA flow cytometry (FC). The two platforms convey complementary information, FC focusing on cell distributions in G1, S, G2M cell cycle phases, TL on lineage trees following cells in subsequent generations. Valuable examples of modelling proliferation over subsequent generations can be found in the immunological field, were data from cells labelled with suitable trackers are exploited to follow lymphocyte expansion after introduction of stimuli triggering their entry in cycle [16], [17]. These models usually do not include cell cycle phases. On the other hand, models based on cell cycle phases were used [18], [19] and adapted to model treatment [20] not distinguishing generations, although effects of treatment were observed not only in cells directly exposed but also to their descendants [21], [22]. In this paper we attempt to build a model deciphering the effects of treatment including both cell cycle phases and subsequent generations, exploiting together FC and TL data. By modelling proliferation at the cell population level it was possible to incorporate different kinds of heterogeneity in the structure of the model, both for untreated and treated cell populations. For untreated cells, intercell heterogeneity was included by distinguishing G1, S and G2M phases and introducing frequency distributions of the duration of each phase. After treatment, observations demonstrate that some cells die other survive, some are blocked in a given phase, and among them some die and others cycle and reach subsequent generations and so on. In this case the heterogeneity is rendered by the model parameters themselves, which are probabilities or rates of occurrence of each phenomenon in play and were applied to the cell cohorts arriving at each checkpoint at a given time. A multi-phase and multi-generation model for proliferation in the absence of treatment was already presented [23] showing the feasibility to reproduce proliferation at a deeper level than models explaining simply the time course of the overall cell number or equivalent measures, or cell phase percentages alone or cell generations alone. This model was here modified to include the effects of a treatment. Our previous attempts in treatment modelling were based on “single-cycle” models, where cells after division re-enter in the same cycle and generations were mixed, and on FC data, alone or coupled with absolute cell counts, with simplified descriptions of the time-dependence of checkpoint activities, and trial-and-error simulations [24]. These treatment models were revised and modules rendering the main effects of treatments were included in each phase and generation of the multi-phase and multi-generation model. This eventually enabled us to explain most of time-dependences required in the previous treatment models as different effects in subsequent generations. Moreover, a non-linear fitting procedure was included and, correspondingly, the experimental design was strongly enriched, with the inclusion of TL data, for different treatment doses. Setting the method obliged us to study and solve a number of related subproblems, concerning the data structure, the cell cycle model and the optimization (Figure 1). A robust rendering of the proliferation process up to the fourth generation was achieved, designing a formal workflow for model building and fitting, considering contemporaneously data from treatments at different doses and exploiting cross-validation between FC and TL datasets. This approach was used here to disclose the details of the proliferation of a cancer cell population after X-ray exposure, leading to a measure of the dose- and time-dependence of the main checkpoint activities and of the balance between recycling and death over three cell generations after treatment. Reconstruction of the process was then used for in silico experiments, evaluating the consequences of default or potentiation of specific checkpoints on the overall effect of treatment.
Human ovarian cancer cells (IGROV-1) in exponential growth were exposed to different doses of X-rays (0 h). The response to the treatment was monitored over time using TL and FC (Figure 2). Two FC techniques were applied (Figure 2A): i) monoparametric cell cycle analysis, which gave FC DNA histograms, ii) BrdU pulse-chase, which gave biparametric DNA-BrdU histograms at different times after treatment. Through a first-level analysis (supplementary Text S1), we obtained time-courses of FC cell percentages in cell cycle phases and BrdU subsets, which were collected in the FC database. TL movies (see representative examples as supplementary Videos 1 and 2) were analysed tracking all the cells in the fields of view at 0 h, defined as “generation 0” (gen0), and their descendants, constituting gen1, gen2 etc. Observations were collected in a “lineage database” reporting for each cell: i) a cell identification code, including information about treatment, field, lineage and generation number; ii) time of birth; iii) time and kind of outcome (mitosis (M), death (D), survival at 72 h (S), re-fusion of two newborn cells to make a polyploid one (R) or loss from the field of view (FL) ) (Figure 2B). Other events, like anomalous mitoses (resulting in three or four offspring), occurred with a frequency lower than 1% and were neglected in this analysis. Lineage data were then analysed, deriving cell population statistics. For each dose we calculated: i) the frequency distribution of intermitotic times; ii) the time course of the number of cells in each generation and of polyploid cells, relative to the number of cells present at 0 h and corrected for FL cells; iii) the frequency of the cell events in each generation. These data made up the three sections of the TL database. Qualitative visual inspection of FC and TL data after this first level of analysis already enabled us to depict a preliminary view of the effects of X-ray exposure (supplementary Text S2). TL data (Figure S1) indicated the presence of cell cycle delays even at 0. 5 Gy, and cell death was detectable from 2. 5 Gy, and up to gen2; some sibling cells eventually re-fused, producing a subset of polyploid cells within which mitosis was a rare event. FC data provided complementary information, indicating the activation of blocks/delays in all cell cycle phases, but without distinguishing different generations (Figures S2, S3). The overall process of proliferation of cancer cell populations involves the passage of cell cohorts through cycle phases in subsequent generations. Treatment changes the cycling (i. e. the “unperturbed” cycling of untreated controls), imposing delays and blocks and killing cells with different mechanisms in each phase. This is reflected in the structure of the modelling framework implemented in silico (Figure 3A). Response to treatment was modelled with modules for G1, S and G2M checkpoint response in gen0, gen1, gen2 etc. superimposing and modifying the flow of the cells through the cycle. An additional module governs polyploidization at the end of each generation. The user interface of the software enables the researcher to build the desired model, including the modules of choice for G1, S and G2M in gen0, gen1 etc. (see supplementary Dataset S2). The rendering of the proliferation works on cohorts of cells whose positions within G1, S or G2M are tracked by their age (time from entering) in that phase. Cells are grouped in each phase within age compartments of length Δ and straightforward balance equations connect the age distributions of cells at time t to those at time t+Δ (Figure S4). When a group of cells ends the cycle and divides, their offspring enter the next generation at the beginning of G1. This modelling of the cell flow through the cycle was extensively tested in previous studies [24]–[27] and is consistent with the general mathematical theory of age-structured cell populations [24]–[30], with the simple assumption of the balance of the number of cells entering, lost and exiting in each compartment. Recently, we introduced a separate simulation for cell populations of different generations, using this model to render the proliferation of a cell population by fitting at once FC and TL data [23]. The section “Modelling proliferation of untreated cells” recapitulates the model, including a further generalization for future uses. Using this model, cycling of the untreated IGROV-1 cell population was simulated (Figure S5) on the basis of average and coefficient of variation of the time spent in each phase and refined including quiescence characteristics for the cell line under study, calculated in preliminary TL and FC experiments, during asynchronous exponential growth. In addition to measuring cell cycle percentages and the number of cells in each generation, these experiments included the measure of intermitotic times, a detailed time course after BrdU pulse labelling and evaluation of quiescent cells by continuous BrdU labelling. In order to render proliferation in the presence of an anticancer treatment, the basic proliferation model was nested with a model of the response to treatment, considering that the response can be very different not only in each cell cycle phase but also in subsequent generations. For this purpose, cycling of treated samples was simulated by perturbing the basic flow of untreated cells with modules that can be applied to any phase. Module types for G1 phase are shown in Figure 3, panels B to D (see Computational Methods for details). Modules render complex biological phenomena, e. g. the operation of a checkpoint in a specific cell cycle phase, with the minimum choice of parameters that enables quantifying the main antiproliferative effects (block, block recovery or death) occurring in that phase. The cycling cells death module (Figure 3B) mimics an immediate cell loss, possibly occurring at high doses with no detectable previous cell cycle arrest. The mitosis/polyploidy module (Figure 3C) deviates polyploid cells to a specific pool. Four different types of modules render checkpoint activity with increasing levels of complexity (Figure 3D). At the lowest level (type I), the effect is simply a delay of the progression in the phase where it is located. Type II acts permanently arresting a fraction of the cells transiting to the next phase in a specific compartment of blocked cells, which then may die or not, according to a “death probability” parameter. Type II can be applied together with a type I module, modelling a situation where some cells are permanently arrested and eventually die, while the others are simply delayed. Type III aims at rendering the competition between repair and cell death in blocked cells, adding a new parameter describing recycling to the type II module. Type IV adds further complexity to type III, introducing time-dependence for the onset of one or more of the effects (block, recycling or death). Parameters of each module are probabilities (e. g. a “block probability”) or rates (e. g. a “death rate”) of a specific event in the phase and generation where the checkpoint is located. For any set of values of the input parameters of all modules, the software (supplementary Text S4) gives the simulation of the entire time course of the progression of the cell population within the phases and in subsequent generations, of both BrdU− and BrdU+ cycles, and plots the results. Outputs are simulated data equivalent to those obtainable with different experiments, like cell cycle flow cytometry, time lapse imaging, Coulter counter and potentially any proliferation-based test (see Computational Methods). A rational workflow (Figure S6 and supplementary Text S4) was devised to build progressively the most suitable model of X-ray treatment, starting from the simplest and adding complexity when required, reaching a balance between the desired detail of rendering the process and the available data, avoiding over-simplification on one hand and over-parameterization on the other. Testing of tentative models was based on non-linear fitting of individual doses first, then performing a multi-dose fitting considering all data together. For multi-dose fitting, single dose parameters were constrained to obey simple dose-response relationships, using Hill or gamma functions. The workflow was successfully applied to fit simultaneously all FC and TL data of all doses in the experiment of X-ray exposure of IGROV-1 cells (Figure 4). The resulting best fit with the final model gives a dynamic (0–72 h) rendering of the flow of cells through the phases of the cell cycle and through subsequent generations, for controls and each X-ray dose (Videos S1 to S5). Figure 5 shows the cell cycle distribution at representative times and doses taken from the Movies 1 (control), 2 (0. 5 Gy) and 4 (5 Gy). In particular the 0. 5 Gy panels show that cells arrested in gen1 at 24 h were no longer there at 72 h, and most cells (85%) were cycling in gen3 at 72 h. Instead, at 72 h, 5 Gy-treated cells are distributed in the G1 and G2M blocks in gen1 and gen2 (40%) and in the polyploidy compartments (17%), with few (10%) cycling cells in gen3. The details of the dynamics of the process can be appreciated in the respective Supplementary Videos. Supplementary Dataset S3 includes FC and TL data of all X-rays doses and the user interface to simulate the final model, together with the instructions to download and run the simulation program as a Matlab standalone executable file (PaoSim_MultiDose). Rendering the proliferation in silico enabled us to retrieve information not available by direct observation and to quantify the processes in play. In particular, the best fit provided a quantitative evaluation of the activity of each checkpoint, separately in irradiated cells (gen0) and their descendants (Figure 6). We found that the G1 checkpoint was activated in cells directly exposed to radiation at all doses (Figure 6A, gen0), inducing a partial G1 block immediately after exposure. Blocking activity was relatively weak at 0. 5 Gy and high at 2. 5 Gy, with no further increase at 5 and 10 Gy. G2M block probability was similar to that in G1 at 2. 5 Gy, but was stronger and continued to increase at higher doses. Among cells reaching G2M in gen0 (i. e. not only cells originally in G2M but also those irradiated in G1 and S that progressed in the cycle) >95% were intercepted there at 10 Gy. The fate of cells blocked in G1 and G2M was different. All G1-blocked cells were able to recycle only with 0. 5 Gy, then gradually the recycling decreased and death increased, and at 10 Gy the majority of these cells eventually died (Figure 6B, gen0). In contrast, the balance between repair and death among G2M-blocked cells in gen0 was completely on the side of repair, as at least 95% of those cells eventually succeeded in dividing and entering gen1 also at the highest dose. After the first division, cells irradiated with 0. 5 Gy were minimally perturbed in gen1 and gen2,25% of them being temporarily arrested or delayed in G1 (Figure 6A, gen1 and gen2). At higher doses, cell cycling was strongly hampered in all phases. 50% of the cells traversing G1 were intercepted and blocked there with 2. 5 Gy or higher doses (Figure 6A, gen1 and gen2). Most of these G1-blocked cells died or were still arrested at 72 h (Figure 6B, gen1 and gen2). A dose-dependent G2M block was detected, reaching 0. 8 probability at 10 Gy (Figure 6A, gen1 and gen2), with different outcome in gen1 and gen2. The majority of cells recovered from G2M block up to 5 Gy in gen1, while in gen2 most G2M-blocked cells died (Figure 6C, gen1 and gen2). At 10 Gy mortality in the G2M block was high in both generations. The activity of S checkpoint is shown in Figure 6D, reporting the reduction of the rate of DNA synthesis. Gen0 cells exposed to X-rays while in S phase (BrdU+) were only moderately delayed in that phase and eventually reached G2M. The delay was stronger at 10 Gy for cells entering S after having overcome the G1 checkpoint (BrdU−), requiring more than 20 h on average to complete DNA synthesis. Cells in gen1 were then strongly delayed while traversing S phase at 2. 5 Gy and higher doses, and no S checkpoint activity was observed in gen2. Cell re-fusion was measured by the polyploidization module, acting on cells entering gen1 and gen2. Polyploidization was dose-dependent (Figure 6E), with about one third of re-fusions among cells that divided after 10 Gy. Most polyploid cells survived without dividing up to the end of observation, at 72 h, so that the fate of these cells remained unknown. The final model can be used for in silico experiments, addressing specific questions. In a first experiment in silico we explored the impact on the overall proliferation of checkpoint activity in gen1 and gen2. We ran the simulation assuming the effects of radiation would be limited to the cells directly exposed to X-rays, with no effects on their descendants. In this scenario, we found that cancer cell expansion was not stopped even at the highest dose (Figure 7A), demonstrating that the events in gen1 and gen2 are crucial for the response to treatment. At the lowest dose, in the absence of effects on descendants, the growth curve almost overlaps the growth of untreated cells, so most of the modest growth inhibition has to be attributed to perturbations in gen1 and gen2. As the dose increases, the perturbations in gen0 become more obvious, but cell cycle checkpoint activities on the descendants of irradiated cells still make a major contribution to the overall effect. Other in silico experiments were run strengthening the G1 block or preventing the exit from G2M block, miming hypothetical co-treatments with agents eliciting these specific effects in combination with X-rays. Figure 7B shows the effects of co-treatments with two hypothetical drugs: drug A acting in G1, resulting in potentiation of the G1 block, drug B acting in G2M, inhibiting recycling from that phase, possibly acting on one of the stages of the repair process. At 0. 5 Gy the effect of drug B is irrelevant, in view of the limited importance of the G2M checkpoint at this dose, while drug A, acting on G1, would potentiate the X-ray effects. At 2. 5 Gy drug B would potentiate X-rays more than drug A in the short term, but in a long term the differences tend to disappear and eventually re-growth would be similar in both instances. At very high doses, the effects of drug A would become irrelevant, because the G1 checkpoint is already strong with X-rays alone. Drug B would be more effective in the short term, but less on a longer term, as it would prevent cells reaching the G1 checkpoint in gen1 and gen2 where cytolethal activity was greater. The results suggest that strengthening the G1 block would potentiate the effects of radiation, but only at intermediate doses. A drug preventing the exit from G2M block would be more effective in the short term, but less on a longer term.
Mathematical modelling of biological or biomedical systems is now facing the challenge of quantitative rendering of the biological phenomena, taking into consideration the vast amounts of data produced by experimental biomedicine, with appreciation of the biological complexity, with its dynamics, and the processes of measure, with their strengths and limitations. Also in cancer research the demand for proper quantitative interpretation of the processes in play is rising [31], [32]. The method we propose enables a detailed study of the mechanisms concurring to build the response to anticancer treatment at the “cell population” level, where measures were made on single cells and on a statistically representative number of them, accounting for inter-cell variability. This provides a different view respect to that obtained with bulk measures on cell homogenates, representing “average” quantities, unpredictably weighted by the relative numerosity of cell subsets with different amounts of the measured quantity. At the cell population level, we may exploit a very rich amount of information, coming from established or new techniques. It remains a challenge to find the way to organize all that we can measure, and possibly put together complementary information from different techniques, in order to give sense to the observed heterogeneity. This challenge must be addressed, if we really aim to reach an unbiased picture of the behaviour of a multicellular system, from in vitro cultures to in vivo organs or tumours. We demonstrated that it is possible to integrate FC and TL platforms via in silico rendering of the dynamics of the underlying cell cycle progression. This was achieved with a flexible model-building environment, which enabled us to challenge a variety of tentative models in a manageable time. Once the model of unperturbed proliferation was established, we progressed to study how the system responds to a hypothetical treatment. Different types of modules were implemented in the model, miming the activity of cell cycle controls at increasing levels of complexity. The user can select a specific effect and look at its consequences in simulated experiments. Wishing to build a model for a real treatment, we strengthened the parameter estimation by fitting data from different treatment levels (four X-ray doses in our example) simultaneously and optimizing the dose-response for each parameter, with simple regularity assumptions. This combination of data from different doses and platforms provided internal cross-validation to exclude the existence of multiple solutions in models including several variable parameters. Our optimization procedure initially establishes a relatively simple model with a unique best fit solution, which was then cautiously refined, up to the definition of a final model, by which the data of all doses were fitted with average errors 1. 6% (FC cell cycle percentages) and 3. 6% (TL cell generation numbers) (Figure 4). The final model provided a unifying view of the effects of radiation that puts in order and defines the relative importance at the different doses of the main phenomena in play, which are described separately in the literature [21], [22], [33], [34] (supplementary Text S5). Study of the final model disclosed the dynamics of checkpoint activities and their dose-dependence, showing up basic differences between G1 and G2M checkpoint responses that were not perceivable by simple data inspection. In particular, G2M-blocked cells were eventually able to recover from the block and divide even after treatment with 10 Gy, while at this dose the majority of G1-blocked cells eventually died there. The response was not restricted to gen0 cells, but cell cycle progression was also hampered in gen1 and gen2. Our results strongly support the concept of a limited efficiency of a single checkpoint to process the damage [35], suggesting that cycling of individual cells may be subsequently hampered in different phases. The model was also used as a basis to explore, with in silico experiments, the consequences of default or potentiation of specific checkpoints, which in principle could become feasible with new, molecularly targeted co-treatments, or may occur in tumours with particular genetic defects. Obviously, the real response to treatment is even more complex than the proposed model. It does not include, for instance, possible correlations between sibling cells for probabilities of specific outcomes, or signalling between cells or spatial issues driving neighbour cells towards a common fate. Still our model works at a quite high complexity level, including details of the effects in each phase and generation at different treatment doses. This was made possible avoiding over-parametrization only through fitting of a very rich and complex data set. Further details would require more parameters and additional data, which are not precluded in future studies.
Asynchronous proliferation of untreated cell populations in the biological system is achieved by a balance of the cell cycling process, quiescence and cell death [37]. The multi-generation cell cycle model, with variable phase durations but without quiescence and death, is shown in Figure S4 and Supplementary Text S3. State variables are G1 (k, t, geni), S (k, t, geni), G2M (k, t, geni), giving the number of cells in generation geni, with age k at time t in the respective cell cycle phases. Model parameters are (, CVG1, , CVS, , CVG2M), i. e. the average and coefficient of variation of the phase durations, on which basis the exit probability for each compartment () were calculated (see supplementary Text S3). Because quiescence and spontaneous cell death are often not negligible even in untreated cell populations growing in vitro in optimal environmental situations, quiescence modules were included in G1, S and G2M to render these phenomena. Simulation is intended to start at a laboratory time (0 h) corresponding to the start of a treatment in the experimental setting. Thus the model requires the input of the age distributions at 0 h (G1 (k, 0,0), S (k, 0,0) and G2M (k, 0,0) ) from which the time evolution is simulated according to the previous equations. The user can provide any starting distribution, but we chose the asynchronous distribution in all our studies aiming at rendering experiments with non-synchronised cell populations. Any cell population, under constant environmental conditions, reaches asynchronous exponential growth, and good laboratory practice requires that cells are in this condition (or as near as possible) before an in vitro treatment, for data reproducibility. In this condition the probability distribution of cells in the cell cycle phases and ages is time-independent. A desynchronization routine of our program automatically calculates the asynchronous distribution, by running the cell cycling procedure from an approximated initial distribution until the cell cycle percentages varied by less than the desired precision (0. 1% in the present study). Then the G1, S, and G2M age distributions were normalized in order to start the cell cycle simulation with a given number of cells at 0 h (e. g. 1000 cells), providing G1 (k, 0,0), S (k, 0,0) and G2 (k, 0,0). A program performing the desynchronization routine is freely available as specified in supplementary Dataset S1 (see also the Software section in Supplementary Text S4). To simulate pulse-chase BrdU experiments, we run two cell cycle simulations independently for BrdU− cells and BrdU+ cells. The effect of pulse labeling was reproduced by splitting the starting distribution, assigning G1 (k, 0,0) and G2 (k, 0,0) (plus quiescent cells) to the BrdU− cycle and S (k, 0,0) to the BrdU+ cycle. Thereafter the two cell cycle simulations proceeded independently. The program can also reproduce the effect of labeling of any duration (if required for specific experiments, e. g. with continuous BrdU labeling), transferring cells entering S phase from the BrdU− to the BrdU+ cycle up to the time when BrdU is removed. These features guarantee maximum flexibility to the program, to simulate different experimental plans. At the end of a simulation run, the program derives the time courses of several kinds of simulated data with simple calculations from the time courses of the state variables. First, the proportions of cells in G1, S and G2M are trivially calculated by summing the number of cells in all compartments of that phase and dividing by the total number of cells. These predictions are compared with experimental FC percentages at the measurement times. Then the program calculates the time courses of the number of cells in each generation (Ngen (i) (t) ), summing up the number of cells in all compartments of all phases of the generation. Predicted Ngen (i) (t) s are directly compared with the corresponding TL data. The time course of overall cell number (N (t) ) is also calculated, summing the number of cells in all generations, and can be compared with experimental measures obtained by TL or other techniques (e. g. Coulter counter). Similarly, simulated percentages of cells in the BrdU− and BrdU+ subsets are calculated summing the number of cells in the appropriate subset compartments and plotted with experimental counterparts. The time course of the percentages of undivided BrdU+ cells (%Und+ (t) ) is obtained by summing the numbers of cells in all gen0 compartments of the BrdU+ cell cycle and dividing by the overall cell number. Another valuable item is the the distribution of Tc, measurable by TL. It comes from frequency measures of intermitotic times in a sampling of the whole population (we had a total of 2480 cells), which were normalized and assumed to be representative of the distribution of Tc in the whole population (experimental F (Tc) ). While the simulated average is easily calculated directly in our model as, and compared with the corresponding datum, a specific routine was designed to calculate the whole F (Tc), on the basis of FG1 (k), FS (k) and FG2M (k). Considering two phases, e. g. S and G2M, the joint probability that a cell traverses S and G2M with times kS and kG2M is FS (kS) ×FG2M (kG2M). Thus the probability that the duration of S+G2M is kSG2M is given by the formula: F (kSG2M) = ∑FS (kS) ×FG2M (kG2M) where the summation includes all combinations of kS and kG2M, such that kS+kG2M = kSG2M. Then F (Tc) is calculated as F (kG1SG2M) = ∑ FG1 (kG1) ×F (kSG2M) where the summation includes all combinations of kG1 and kSG2M, such that kG1+kSG2M = kG1SG2M. Simulated data are derived in the same way when perturbation modules are included except in the case of F (Tc), which was calculated only for untreated cells. Summarizing the simulation procedure, upon input of the parameter values, the program: i) calculates Fph (k) and βph (k), ii) runs the subroutine to produce the asynchronous starting cell distribution, iii) simulates the movement of the cells through the cell cycle in subsequent generations, calculating the time course of the state variables up to the desired end-time, and iv) calculates the time course of all derived quantities (such as %G1 (t), %S (t) and %G2M (t), %BrdU+ (t), %Und+ (t), , F (Tc), Ngen (i) (t) etc.) and updates their plots with the corresponding experimental data. Perturbations of cell cycling induced by treatment are modelled with modules for G1, S and G2M response in the gen0, gen1, gen2 etc. , superimposing and modifying the “unperturbed” flow of the cell through the cycle. An additional module governs polyploidization at the end of each generation. The user can choose among several modules types, emulating different biological effects, and locate them in any phase and generation, so that it is possible to deal with the different antiproliferative effects of anticancer treatments without modifications of the software. Supplementary Dataset S2 includes the user interface to build a single-dose model, together with the instructions to download and run the simulation program as a Matlab standalone executable file (PaoSim_SingleDose). The interested reader can also test the model for untreated cells using the PaoSim_SingleDose program with no perturbation modules. During the optimization procedure, the choice of the module and module type for each phase was guided by the law of parsimony: in general we applied the module type that satisfactorily reproduced the data with the smallest number of parameters. The optimization procedure therefore involved several steps, as described in supplementary Text S4.
The antiproliferative response to anticancer treatment is the result of concurrent effects in all cell cycle phases, where molecular control pathways (checkpoints) are activated and cells may be arrested to repair DNA damage or killed if not able to succeed in the repair process. The complexity and inter-cell variability of these phenomena are not captured by the available methods, and the origin of the dose-dependence of the response remains elusive. In this work, we present an experimental-computational method that discloses and measures the individual responses of cell cycle controls in each phase and generation. We demonstrate that the method, exploiting jointly data sets obtained by flow cytometry and time-lapse in vivo imaging with a suitable experimental design, is able to achieve a full reconstruction in silico of the actual movement of cell cohorts following X-ray exposure, providing separate and quantitative measures of the dose-dependence of G1, S and G2M checkpoint activities in subsequent generations. Best fit parameters values are actual measures of the probability of activation of the specific pathways of arrest, repair or death within the cell population, linking the molecular scale to the “macroscopic” response, with full appreciation of its dynamics and inter-cell heterogeneity.
Dynamic Rendering of the Heterogeneous Cell Response to Anticancer Treatments