UNRAVELING THE COMPLEXITY OF CELL CYCLE EFFECTS OF ANTICANCER DRUGS IN CELL POPULATIONS

. Cell cycle perturbations occur after treatment with all anticancer drugs. The perturbations are usually classiﬁed as cytostatic (cell cycle arrest) or cytotoxic (cell killing). Our approach for analysis of cell cycle perturbations in vitro was to consider all the data provided by diﬀerent experimental tests and interpret them through a mathematical formulation of the problem. The model adopted for data analysis and interpretation is the result of merging two models, one for the cell cycle and the other for the drug eﬀects. The ﬁrst exploits the results of the theory of age-structured cell population dynamics while the second is based on distinct parameters (”eﬀect descriptors”) directly linked to cell cycle arrest, damage repair or cell death in G 1 and G 2 M and to inhibition of DNA synthesis and death in S . The set of values of the eﬀect descriptors which are coherent with all experimental data are used to estimate the cytostatic and cytotoxic eﬀects separately. Applying the procedure to data from in vitro experiments, we found complex but biologically consistent patterns of time and dose dependence for each cell cycle eﬀect descriptor, opening the way for a link to the parallel changes in the molecular pathways associated with each eﬀect.


1.
Introduction. Mathematical theories of the behaviour of cell populations rely on variables associated with cell cycle progression. Cell volume, for instance, obviously increases in the normal cell cycle from birth to division. However, it does not represent cell cycle progression well. Cell age is probably the variable of choice, if we want to make kinetics considerations, and eventually outline the time evolution of certain global quantities of a cell population.
If we want to base the theory on experimental grounds, we will probably be disappointed that the quantities that can be measured differ from the mathematical variables. In addition, there is no measurable marker of cell age and biological measures are not as precise, reproducible, and independent of the "biological model" (i.e. a given cell line) as in engineering, physics or chemistry.
In our approach, we aimed to bridge the gap between theory and experiment. The starting point is the experimental measures, in particular the DNA content, which so far is the measurable quantity nearest to the mathematical "cell age" variable. Thus, in the first sections we describe the experimental techniques and the relationships between the measures and the basic kinetics quantities. Then, 324 PAOLO UBEZIO we outline a mathematical theory based on the measurable quantities, for both "unperturbed" and "perturbed" (by drugs or radiation) cell populations.
2. Experimental Techniques for Studying Cell Populations. The technical evolution of microscopy, as a way to look at cells, followed two directions in the 1960s. The first, aimed at overcoming the subjective nature of microscopic observation, led to computerised image analysis; the second, aimed at real measures of particular cell features and representativeness of the whole population, led to the design of Coulter counters and flow cytometers. Those are both known as "zeroresolution" instruments, because they do not reconstruct an image of the cells, but measure one parameter (or a few) on a large number of individual cells, suspended in a fluid and driven to a sensing zone at high speed.
Coulter counters measure the number of cells growing, for example, in a tissue culture flask, so cell population growth can be studied in normal conditions and under the influence of drugs. Flow cytometry measures the cell content of any substance for which there is a specific fluorescent marker, after cytochemical treatment of the cells with that probe. We can record the fluorescence intensity of each cell and build up a frequency distribution for it on thousands of cells analysed, statistically representative of the whole population. The distribution of the fluorescence intensity is then interpreted as the frequency distribution of the cell content of the substance recognised by the marker.
The number of substances that can be measured by flow cytometry is enormous, exploiting antibody technology. Still important is the analysis of DNA content, using specific fluorochromes. Studies of cell cycle and cell population dynamics still rely on measurement of DNA content, as a marker of the cells position in the cycle. Figure 1 shows schematically the procedure for measuring DNA by flow cytometry.
The histogram shown in Figure 1 is readily interpreted on the basis of the cell cycle. Since 1951 [5] we have known that DNA synthesis does not start immediately in newborn cells, but is confined to a subsequent period ("phase S") which ends some time before the onset of mitosis. Thus, the mitotic cell cycle is divided into four phases: i) "phase G 1 " from birth to the start of DNA synthesis, when cells have DNA content "x"; ii) "phase S", when cells progressively increase their DNA content from "x" to "2x"; iii) "phase G 2 " from the end of DNA synthesis to the start of mitosis, when cells have DNA content "2x"; iv) "phase M ", the mitotic phase, still with DNA content "2x". This is easily confirmed by flow cytometric DNA histograms of proliferating cells, where a first peak (around the value x = 100 in Figure 1) is detected, grouping all G1 cells (usually this is the biggest peak), followed by a region of distributed cells, with no peaks around any particular value ("phase S" cells), and ending with a second peak (around the value "2x") grouping G 2 and M cells.
DNA histograms can then be further analysed to retrieve the percentages of cells in the various phases (%G 1 , %S, %G 2 M ). Unfortunately, this may call for sophisticated fitting procedures, because S-phase cells are not detected separately from G 1 cells, as the two subpopulations are partially superimposed due to instrumental broadening, and the same occurs in the region between S and G 2 M cells.  2) The flow cytometer sucks the cells from the tube, aligns and carries them, one at a time, to the sensing zone, where the excitation light beam is focused. As each cell passes through the light beam, the excited PI molecules generate a fluorescence signal. Several detectors collect this light pulse; each connected to an optical path selecting a particular wavelength range. One optical path (red) is optimised for detection of PI fluorescence.
3) The instrument converts the optical pulse into an electrical signal, which is amplified, measured, digitalized and recorded. In a few minutes thousands of cells are analysed and the frequency distribution of DNA content is visualised.
In the example shown there is a first peak of G 1 cells, a flat (S-phase) region and a second peak of G 2 M cells.
In mathematical terms, the relationships between the measured fluorescence distribution f (y) and the unknown DNA distribution g(x) is given by the equation where p(y|x) is the probability density of fluorescence y for a cell with a DNA content x, and e(y) is an error term.
3. Characterisation of Asynchronous Exponential Growth by Flow Cytometry. Once we have measured cell phase percentages from DNA content analysis by flow cytometry, the next question is how these data can be used to study the kinetics of a cell population. We need a link between flow cytometric quantities and the basic kinetic features of a cell population, such as the durations of cell cycle phases. This link can be established in asynchronous cell populations.
In the experimental practice, a proliferating cell population is said to be asynchronous when the percentage of cells in every phase of the cell cycle remains constant in time (within the limits of experimental error). At the same time, the cell population grows exponentially (again within the range of experimental error). This is a frequent biological observation, when environmental and nutritional changes are slow within the observation period. Mathematically, in the framework of age-structured cell populations, an asynchronous population presents a constant age density distribution, with perfectly exponential cell number growth.
We can bridge the gap between experimental data and mathematical theory by including the phase structure in the theory. A very simple mathematical model of the cell cycle can be built (see [8]), assuming no inter-cell variability (i.e. the duration of G 1 is T G1 in all cells, and T S and T G2M for S and G 2 M length), no quiescence, and no cell loss. This gives the following one-to-one relationship between the percentages of cells in G 1 , S and G 2 M (measurable quantities) and the fraction of T D (doubling time) spent in the various phases: Steels formulae estimates T G1 , T S and T G2M , since T D = T C in this model and can be independently measured from the growth curve (e.g. using a Coulter counter).
Starting from this basic model several elements of complexity can be added, bringing the mathematical description closer to reality.
First, inter-cell variability of phase durations can be introduced, because without this variability the asynchronous status will never be reached. The mathematical equations become much harder to manage, but fortunately computer simulations (see [10]) show that Steel's formulae still hold good, re-interpreting T G1 , T S , T G2M as mean values for the distributions of phase transit times. Thus, inter-cell variability acts as a potent motor that makes a cell population asynchronous in a short time, allowing us to apply the theory of asynchronous populations to real situations. When asynchronous growth is reached, we can measure DNA content and apply the simplest cell cycle theory to estimate the mean cell cycle times. The procedure has been applied to data obtained in vitro but in order to approach the in vivo situation we have to introduce quiescence and cell loss.
To introduce a "quiescent" status into our mathematical model, we might assume that a fraction of the newborn cells will enter it; this will not necessarily happen immediately after birth, but for the sake of simplicity we can consider as quiescent any cell that will become quiescent. Bertuzzi et al. (see [3] and [4]) extended this analysis, admitting (regular) re-cycling of quiescent cells. In the asynchronous steady state, the fraction of proliferating cells (growth fraction, GF ) is constant, meaning that the quiescent and proliferating cell numbers increase the same way as the whole population. Steel's formulae still hold but the (mean) cell cycle time T C is no longer equal to the doubling time, and is also related to the growth fraction, by the formula: This means that T C (or GF ) must be calculated independently, in addition to T D , in order to obtain the value of T G1 , adding a new experimental step. Cell loss normally occurs in vivo and helps keep the cell number dynamically constant in normal cell populations; but it is also present in growing tumor cell populations, where the control of cell number is lost. The rate of cell loss may change in the phases of the cell cycle. A detailed mathematical model including quiescent cells and different rates of cell loss in each phase has been published ( [3] and [4]). A first approach considers uniform cell loss (i.e. not dependent on phase or quiescence). In this case, a new measurable quantity replaces T D in Steel's formulae, namely the potential doubling time (T P OT ), which is the doubling time the population would have if there were no cell loss (see [8], [3] and [4]).
In this way, an originally simple mathematical theory was made more and more complex by introducing new but measurable quantities, improving our modelling in parallel with the experimental requirement. The theory suggested that GF (or T C ) and T P OT should be measured, and indicated how to interpret the data to get a complete picture of the kinetic features of a cell population in a steady state of growth.
The limits of Steel's formulae were eventually wider than the restricted hypotheses used to obtain them analytically. The error introduced by intercellular variability of phase durations is negligible, and the procedure can be adapted to take account of quiescent cells or cell loss.

Multiparametric Flow Cytometry.
A convenient way to measure quantities related to the cell cycle is multiparametric flow cytometry, where additional parameters are measured with DNA in the same cell. The resulting kinetic characterisation is not only of asynchronous but also of perturbed (by treatment) cell populations. The power of these measures lies in the fact that we do not simply obtain frequency distributions of the single parameters, but have full information about how they are correlated. For instance, we can measure DNA content together with protein content. We will get a biparametric DNA/protein content histogram where the protein content of G 1 , S and G 2 M cells can be measured separately, within the same data set.
A multiparameter method very useful in kinetics studies is based on bromodeoxyuridine (BrdU) labelling. BrdU is a thymidine analogue, which is incorporated by growing cells during DNA synthesis in place of thymidine itself. The cells are exposed to BrdU for a while in their environment (in vitro or in vivo), then collected, fixed, double stained (e.g. DNA is marked by PI -red fluorescence-and BrdU by a specific antibody conjugated with fluorescein -green fluorescence-) and analysed by flow cytometry. This gives a biparametric histogram correlating DNA content and BrdU content (Figure 2).
In a pulse-chase experiment, cells are pulse-labelled with BrdU and collected at intervals. In this way, BrdU+ cells, which were in S-phase at the time of labelling, have time to flow out and are detected in G 2 M , G 1 , etc., in the next cycle (incorporated BrdU is transmitted to the descendants) (Figure 3).
A pulse-chase experiment with many time points provides a lot of information about the kinetic properties of the cell population (see [9]), including not only mean phase durations but also their variability. However, this complexity makes analysis of the data anything but straightforward. Nevertheless, a simplified pulse-chase design is used to estimate T S and T P OT in human tumors ( [1] and [2]), allowing the first direct measures of kinetic quantities in vivo in humans.

Flow Cytometric Evaluation of the Effects of Drugs on the Cell Cycle.
Anticancer drugs perturb the growth of normal and tumor cells in several ways, with either "cytostatic" or "cytotoxic" effects. A cytostatic drug should stop cell proliferation, and a cytotoxic drug should kill cells. However, classical anticancer drugs exert both effects, with a strong dependence on the dose and on the schedule of treatment. Nowadays, basic research is focused on molecular differences between normal and tumor cells, on how these differences influence the effect of a drug and on the design of new drugs against specific molecular targets. The molecular origins of the cytostatic and cytotoxic effects have been clarified and the importance of checkpoints in G 1 , G 2 and M is accepted. At these checkpoints, the cell cycle can be interrupted; the cells are blocked until the drug-induced damage is repaired or they are committed to programmed death (apoptosis). Other drugs interfere with DNA synthesis. In this scenario, flow cytometry plays an important role, allowing prompt evaluation of cell cycle perturbations, at least in vitro.
The easiest way to detect cell cycle perturbations is by one-parameter DNA flow cytometry. In a typical in vitro experiment, an asynchronous cell population of tumor cells is treated with a drug, following a schedule of interest, and then cells are collected and analysed at various times. Simple cell cycle perturbations like a complete block of exit from G 1 or G 2 M phases are readily detected by simple inspection of DNA histograms. Most drugs, however, perturb all cell cycle phases in some way.  treated with various concentrations of the drug cisplatin for 1h (see [6]). Cells were collected and analysed at 24h, i.e. the duration of a cell cycle in control cells. It is immediately clear that the effects strongly depend on the dose. At lower doses, the G 2 M peak becomes higher and higher, but the G 1 peak never disappears. This suggests that G 2 M is not the only phase affected by the drug. Moreover, this kind of experiment tells us that blocks or delays in cell cycle phases are not absolute and equal for all cells, even in genetically homogeneous populations such as lines growing in vitro. Only a fraction of cells may fall into a block, often depending on the drug dose. Similarly, blocked cells may exit from the block (after repairing the drug-induced damage), or not, and the exit kinetics is expected to be heterogeneous as well. At higher doses, no accumulation in G 2 M is detected. In these cases, the limits of one-parameter analysis become evident. Whether %G 1 , %S, %G 2 M change or not depends on many factors. First of all, the accumulation of cells in a given phase (for instance, G 2 M ) not only requires a drug effect in that phase (i.e. a block in G 2 M ), but also an influx of cells (i.e. cells have to exit G 1 and S to reach G 2 M ). Therefore, measured cell percentages depend on the history of the cells in the other phases and are the resulting balance between the influx and efflux history in that phase.
Before interpreting our data, we would like a more dynamic picture, giving us a first idea of the fluxes of cells through cell cycle phases. This can be provided by DNA-BrdU biparametric flow cytometry. Figure 5 shows the effect of 50 µM and 200 µM cisplatin in a pulse-chase experiment, 6h after treatment and BrdU incorporation. The progression of BrdU+ cells treated with 50 µM is delayed: they do not shift towards G 2 M as in controls, and none have divided; thus, cells reaching G 2 M have not completed that phase yet. However, BrdU-cells initially in G 2 M are no longer detected there at 6h, meaning they have left this phase. We have therefore confirmation that phases other than G 2 M (namely S phase) are strongly affected by the treatment and the influx to G 2 M is diminished. We also reach a new degree of complexity of cell cycle perturbations: cells may behave differently depending on which phase they were in at the time of treatment. BrdU+ cells in S-phase at the time of treatment are blocked in G 2 M , while BrdU-cells in G 2 M at the time of treatment are not blocked.
The scenario is quite different for the cells treated with 200 µM cisplatin. Here the DNA-BrdU dot plot is different from the control at 6h, but similar to the control at 0h, meaning the cell cycle is frozen in all phases. That is why the %G 2 M does not increase in those cells: not because there is no G 2 M block, but because there is a block in G 2 M , S and G 1 . What happens to these cells? We can bet that at least part will die. In the 24h DNA histogram of 200 µM treated cells, many cell fragments are detected -to the left of the G 1 peak. We can also do other kinds of analysis (flow or static cytometry) to verify whether cell death occurs by apoptosis; the answer is yes in our cisplatin experiment. The important thing is that cell cycle percentages measured by flow cytometry are obviously affected by cell loss. For instance, selective death of G1 and S cells would increase %G 2 M , as would a G 2 M block.
We can in part be helped in the data interpretation by coupling the cell percentages obtained by flow cytometry with a measure of the absolute cell number, for + - Control, 6h Control, 0h instance with a Coulter counter. This was done in the cisplatin experiment. More experimental data were collected in this experiment than are shown here. With such a large amount of data (varying both the cisplatin concentration and the observation time), a coherent interpretation applying intuition to visual inspection of the plots becomes impossible. Without a quantitative description of the cell cycle effects, we cannot catch the complexity of these phenomena. In other words, we really need a mathematical model to read our data. Actually, we need two mathematical models, nested: one for cell cycle fluxes and one for drug effects. We have attempted this approach.

Computer Simulation of Cell Population Growth and Drug Effects.
What we present here is essentially the state of the art of studies started in our laboratory several years ago and still continuing. A more complete description of the method can be found in our publications listed below.

A cell cycle model
To build our cell cycle model, we wanted to keep the essential features we have seen before, allowing the link with experimental flow cytometric data. Thus, we considered three cell cycle phases, G 1 , S and G 2 M , which have a given mean duration with its standard deviation. As we are interested in the behaviour of populations of thousands of cells, and not the fate of single cells, we chose a "deterministic" model, by fixing the probability distribution of phase durations and applying it to each cohort of cells entering a phase at a given time. Starting from the probability Cells completing G 1

Block probability
Killing rate Recycling rate Dead cells Figure 6. Scheme of the three main drug effects in phase G 1 . The parameters "Block probability," "recycling rate" and "killing rate" are time-and dose-dependent.
distribution of duration of a phase, we can calculate the fixed fraction (b(k)) of cells of age "k" completing that phase at any time, while 1 − b(k) progresses in that phase to age k + 1 (see [11], [12]). The alternative would have been a cellby-cell simulation using a Montecarlo approach, which would have enabled us to catch statistical fluctuations of phase transit times, but without an experimental counterpart and with enormously longer computing time. Two parallel cycles can be used simultaneously, one for BrdU+ cells and the other for BrdU-cells, but the same reasoning holds for any two distinguishable subpopulations. At the highest level of complexity, the model comprises separate compartments of definitively quiescent cells in each phase and a "G0-like" phase within G 1 with a first-order exit to G 1 . Distinct cell death rates exist for cells in each phase, quiescent or proliferating, BrdU+ or BrdU-. This allows enormous flexibility in designing any reasonable kinetic scenario for either in vitro or in vivo studies. Once the scenario and a starting age distribution are chosen, the program gives the corresponding time course of cell numbers and percentages in each compartment as output. Any quantity measurable by DNA and DNA-BrdU flow cytometry can be calculated: that is why this program was initially used to study how flow cytometric quantities are related to underlying kinetic characteristics (see [10], [11]).

A model of drug effects
The model of drug effects was designed on the basis of our current biological knowledge, introducing appropriate parameters (see [6]), as descriptors of drug effects in each cell cycle phase. From the cell population point of view, drug response is heterogeneous so we are obliged to describe it in terms of probabilities. The model of the main G 1 responses is represented in Figure 6.
The molecular machinery of the cells may respond to the drug challenge activating a block at a given point (checkpoint) of the cycle inside G 1 . This is described by the G 1 "block probability," which is the fraction of cells (among all cells completing G1 at a given time) that enter the G 1 block, instead of proceeding to S. Thus, this parameter represents the probability of being intercepted at the G 1 checkpoint and blocked there. Blocked cells are collected in a single compartment (without a specific age) and may subsequently either re-enter the cycle or die in the block, depending on the "recycling" and "killing" rates.
The G 1 "Recycling rate" is the proportion of blocked cells re-entering the cycle, in the unit of time. This is indicative of recovery of cells blocked at a checkpoint. The reciprocal of the recycling rate is the mean duration of the block, which gives the time needed to repair, in a broad sense, the damage that caused the block.
The G 1 "Killing rate" is the proportion of cells removed from the group of G 1 blocked cells, in the unit of time. For general application, the present version of the program also includes independent killing rates that act on cycling cells in each phase.
The response in G 2 is modelled the same way as the G 1 response, using a G 2 "block probability," representing the probability of being intercepted at the G 2 checkpoint and blocked there, and G 2 recycling and G 2 killing rates. The main drug effect in S phase is not a checkpoint block, but a delay of progression in the phase, caused by inhibition of DNA synthesis or a reduction of the DNA synthesis rate. This led us to consider a different modelling strategy, loosening the concept of "age" by reducing the fraction of cells (1 − b(k)) expected to mature from age k to k + 1 in the step-time by a factor (Figure 7). This S "Delay rate" is the fraction of cells in each age compartment of S unable to progress to the next age in the unit of time. The overall effect of this parameter is equivalent to a reduction of the average rate of DNA synthesis (reduced by the exact value of the parameter itself), resulting in a lengthening of phase S. In the extreme situation (delay rate = 1) age progression does not occur, indicating a complete "freezing" of the cell cycle.
Interactive optimisation procedure From data for untreated cells (control), the best scenario of cell cycle progression is outlined. This gives the starting age distribution (usually asynchronous) and the reference growth of controls, the baseline on which the drug acts. Then the data for treated cells are examined. For a given input of drug-effect descriptors, the simulation program gives the time course of the output data coherent with that hypothetical scenario. These results are compared directly with all experimental data (absolute overall cell numbers, %G 1 , %S, %G 2 M , % BrdU+, % undivided BrdU+, etc.). As the experimental errors with flow cytometric percentages and cell counts are about 3% and 10% respectively, the fitting is considered satisfactory when all experimental data are reproduced with the same precision. It takes only a few minutes to run a single simulation in any scenario of cell cycle perturbations. Thus, at this stage, many scenarios can be tested with the aim of finding, by trial and error, all -or at any rate as many as possible -the scenarios fitting the data within the experimental range of error. Then, instead of choosing the scenario that fits the data best, we can accept as equally probable all the scenarios fitting the data within the given error and challenge the alternatives by specific additional measures, as suggested by the simulation itself. Eventually, a single set of values for the "block" and "delay" probability, "recycling" and "killing" rates in the G 1 , . How the "delay rate" mimics an effect on the synthesis of DNA. In the absence of treatment, a fraction b(k) of the cells at age "k" in S completes the phase and goes to G 2 , at each steptime. The remaining 1 − b(k) should continue progression in S, reaching age k + 1. When a drug delays DNA synthesis, the latter fraction is reduced by the "delay rate" by an amount equal to the reduction of the DNA synthesis rate. Killing rate is not shown.

Do se ( M)
t=0h-6h t=6h-24h t=24h-48h t=48h-72h G 1 block probability G 1 recycling rate G 1 killing rate  Figure 8. Outlines of the parameters describing the activity of the G 1 checkpoint (block percentage, recycling rate and killing rate), after treatment for 1h with cisplatin, at the concentrations indicated.
S and G 2 M phases provides a coherent picture of the main effects of various drug concentrations over time.
As an example, Figure 8 shows the estimated activity at the G 1 checkpoint (not immediately evident from rough inspection of data) of the cisplatin experiment, in terms of our descriptors of drug effects (see [6]). G 1 block activity was strongest immediately after treatment and then decreased. Recycling after repair mainly occurred during the 6-48h interval, at intermediate doses. The death rate was high at the highest doses immediately after treatment, and late G 1 killing was lower even at intermediate doses.
In conclusion, the method was proved feasible and is now being employed to study both old and new anticancer drugs [7].