AGE-STRUCTURED CELL POPULATION MODEL TO STUDY THE INFLUENCE OF GROWTH FACTORS ON CELL CYCLE DYNAMICS. FRÉDÉRIQUE BILLY

Cell proliferation is controlled by many complex regulatory networks. Our purpose is to analyse, through mathematical modeling, the effects of growth factors on the dynamics of the division cycle in cell populations. Our work is based on an age-structured PDE model of the cell division cycle within a population of cells in a common tissue. Cell proliferation is at its first stages exponential and is thus characterised by its growth exponent, the first eigenvalue of the linear system we consider here, a growth exponent that we will explicitly evaluate from biological data. Moreover, this study relies on recent and innovative imaging data (fluorescence microscopy) that make us able to experimentally determine the parameters of the model and to validate numerical results. This model has allowed us to study the degree of simultaneity of phase transitions within a proliferating cell population and to analyse the role of an increased growth factor concentration in this process. This study thus aims at helping biologists to elicit the impact of growth factor concentration on cell cycle regulation, at making more precise the dynamics of key mechanisms controlling the division cycle in proliferating cell populations, and eventually at establishing theoretical bases for optimised combined anticancer treatments.

the DNA. For instance, in G 1 the cell can check its environmental conditions and decide whether to undergo division, quiescence or apoptosis. The time of this 'decision' is called the restriction point [19,36,47]. This decision is irreversible: once the cell has passed the restriction point, the transition from G 1 to S normally occurs and the cell continues its progression in the cell cycle whatever the environmental conditions. Biological experiments have shown that growth factors were generally responsible for the synthesis of cyclin D, a protein whose binding to cyclin-dependent kinases (CDKs) 4 and 6 is known to regulate the transition from G 1 to S [38,39]. Thus, by improving the cell environmental conditions, an increase in the growth factor concentration is in most cases supposed to reduce the duration of the G 1 phase dynamics, i.e., to increase the probability that the transition from G 1 to S occurs. Although the progression of the cell in phases S, G 2 and M is known to occur whatever the growth factor concentration, it is not yet very clear whether this concentration influences the dynamics of these phases or not. The dynamics of the cell cycle is certainly affected by the level of growth factors in its environment, but much remains to be done to precisely understand the mechanisms that underlie this dynamics.
Proliferation occurs in healthy and tumour tissue. However in cancer cells, various regulation mechanisms are inefficient, which results in uncontrolled tissue growth. In particular, checkpoints, that should induce cancer cell death, are no more efficient and let cancer cells divide. The circadian clock, that exists and is physiologically active in all nucleated cells, is also assumed to play a role on cell cycle regulation in a 24-hour periodic manner. The molecular functioning of this clock is quite complex since more than 15 genes are involved in positive and negative feedback loops. In humans, cell circadian clocks of the whole body are synchronised by a central pacemaker located in the hypothalamus [27,28]. Some biological and clinical experiments have shown that disrupted circadian clocks may enhance tumour growth [15,17,18,21]. Although it has clearly been observed that circadian rhythms play a role in cell cycle regulation, this role has still to be investigated to understand their mechanisms of action, and how circadian control could be used to prevent or treat cancer. In this way, several experiments led by one of us (C. Feillet) resulted in amounts of available biological data related to the cell division cycle with and without circadian control.
In [8], some of us based a theoretical study on biological data to investigate, through mathematical modelling, the population dynamics of healthy and cancer cells and set an optimisation problem for cancer chronotherapeutics, solving it numerically. With the same goal to base ourselves on biological data and eventually help improve cancer treatments, the aim of the present study was, temporarily leaving aside the influence of circadian clocks, to analyse the effects of growth factors on cell cycle regulation and thus on cell proliferation by means of an age-structured mathematical model of cell population dynamics. This paper is organised as follows. Section 2 presents the mathematical model we used and the method we adopted to determine the model parameters. Numerical results are presented in Section 3. A conclusion, together with a discussion and consideration of future works, constitutes Section 4. 2. Mathematical model.  [1,2,3,4,5,6,7,8,10,11,20,23,24,26,35,40,44,46]. The availability of amounts of biological data enables to build more relevant population models.
We consider here age-structured cell cycle models, age referring to age of cells in phases of the cell cycle. The main interest of considering an age-structured model is in distinguishing, in a representation of the cell division cycle, between physiological time (age, taking into account by a structure variable relevant biological variability in a proliferating cell population) and external time. Such models may be relevant in the perspective of controlling the cell cycle by drugs that act on it, as do most anticancer drugs, and in particular cytotoxic drugs.
The cell division cycle being divided into I phases (classically I=4: G 1 , S, G 2 and M ), the evolution of the densities n i (t, x) of cells of age x at time t in phase i is well described by the following McKendrick model [33]: The model for I phases was first introduced in [12]. The particular case I = 1 has received attention from the authors [10,11]. In each phase i of this model, the cells are ageing with speed v i (transport term), they may die (with rate d i ) or proceed to next phase (with rate K i→i+1 ) in which they start with age 0.
Solutions to (1), if the coefficients are time-periodic, or stationary, satisfy n i (t, x) ∼ C 0 N i (t, x)e λt , asymptotically in a L 1 sense with respect to time [37], where N i are defined by (for T −periodic coefficients): (2) The growth exponent λ, first eigenvalue of the system, thus governs the long-time behaviour of the population, since the N i are bounded, according to the normalisation condition of the last equation. The study of λ is therefore of crucial importance. More details about the asymptotic behaviour of (1) can be found in [37].
In the sequel, we will focus on the case of stationary phase transition coefficients (K i→i+1 (t, x) = K i→i+1 (x)), with no death rates (d i (t, x) = 0) and constant velocities (v i (x) = v i ).
As shown in [12], the first eigenvalue λ is then given as the only positive solution to the following equation, which in population dynamics is referred to, in the 1-phase case (I = 1) with no death term, as Lotka's (or Euler-Lotka) equation: In the stationary case with no death rate, integrating Equation (1) along its characteristics, we can obtain the formula: This formula can be interpreted in the following form: the probability that a cell which entered phase i at time t stays for at least a duration x in phase i is given by: where τ i represents the time spent by the cell in the phase i. This leads to the natural consideration of the time τ i spent in phase i as a random variable on [0, +∞[, with probability density function f i : Notice that it is necessary for this interpretation to be coherent, and f i to be a probability density function, to impose +∞ 0 K i→i+1 (x)dx = +∞, which physiologically means that all cells leave any phase in finite time, without any hypothetical maximum age in phase to be introduced in the cell cycle model.

Identification of model parameters.
In this section we present the method we used to identify the model parameters from recent imaging data on individual cells that enabled us to assess the variability of cell cycle phase durations in populations of cells.
2.2.1. FUCCI reporters to identify model parameters. From a biological point of view, as said in Section 1, the cell cycle is classically considered as composed of four phases named G 1 (gap 1), S (DNA synthesis), G 2 (gap 2) and M (mitosis). One challenge of our modelling study was to determine the expression of the parameters K i→i+1 mentioned in the model equations (1) for each phase of the cell cycle (i = 1 . . . 4) (recalling that we assumed d i = 0 and v i constant for all i = 1 . . . 4). To get such an expression, we needed to have access to the distribution of the duration of the phases of the cell cycle within a cell population.
FUCCI is the acronym of fluorescent ubiquitination-based cell cycle indicator. This is a recently developed technique that allows tracking progression within the cell cycle of an individual cell with a high degree of contrast [41,42]. The FUCCI method consists in developing two fluorescent probes indicating whether a tracked cell is in the G 1 phase or in one of the phases S, G 2 or M of the cell cycle. Sakaue-Sawano et al. [41,42] fused red-and green-emitting fluorescent proteins to proteins called Cdt1 and Geminin. Cdt1 and Geminin oscillate reciprocally: Cdt1 level is highest in the G 1 phase and falls down when the cell enters the S phase, whereas Geminin level is highest in the S, G 2 and M phases and falls when the cell enters the G 1 phase. Let us mention that Cdt1 and Geminin are degraded due to the process of ubiquitination, which is what is referred to ("U") in the name of the reporter method. Consequently, the nucleus of a FUCCI cell fluoresces in red when this cell is in the G 1 phase of the cell cycle, and in green when it is in phases S, G 2 or M . This method allows to measure the time a tracked cell spends in the phase G 1 and in the combined phase S/G 2 /M of the cell cycle. Thus, by tracking each cell in a population, we can get the distributions of the duration of these phases within the population, and so we can deduce the probability density functions of the random variables representing the duration of these phases (see below Section 2.2.3 for details).

Identification procedure.
For the parameter identification procedure we used FUCCI data transmitted to us within the C5Sys EU project by F. Delaunay's team, Institute of Biology Valrose, CNRS UMR 7277, INSERM U1091, Université Nice-Sophia-Antipolis, France. The cell lines were obtained by one of us (C. Feillet) by recloning cell cycle phase and circadian clock (through Rev-Erb α) markers and generating a stable NIH 3T3 cell line (mouse embryonic fibroblasts). To get fluorescences that were compatible with the one of the circadian clock marker, that was a yellow-green fluorescence, red-and far-red-(instead of red-and green-, cf Section 2.2.1) emitting fluorescent proteins (mKO2 and E2-Crimson) were respectively fused to Cdt1 and Geminin. So the nucleus of each of our FUCCI cells fluoresced in red when this cell was in the G 1 phase of the cell cycle, and in far-red when it was in phases S, G 2 or M . In the sequel, for historical reasons and to avoid confusion between red and far-red fluorescences, we will consider the initial color code according to which the nucleus of a FUCCI cell fluoresces in red when this cell is in the G 1 phase of the cell cycle and in green when it is in phases S, G 2 or M .
These NIH 3T3 cells were proliferating in a liquid medium (Dulbecco's Modified Eagle Medium, DMEM) with foetal bovine serum (FBS). FBS is known to contain several growth factors such as cytokines and proteins. Two experiments were performed on non confluent (i.e., proliferating) cells according to the concentration of FBS in the medium: 10% (usually considered as the "reference" concentration) and 15%. To avoid synchronisation of cells due to changing culture medium, cells were kept proliferating three days in the chosen medium before the measurements began. Thus, as measurements began, cells were at different stages of the cell cycle.
To compare numerically equivalent populations of cells in the two qualitative groups, 10% and 15% FBS medium, and as cells had also a fluorescent marker for their circadian clock (through the expression of the protein RevErb-α), we decided to keep only within the (obviously more numerous) 15% FBS group only those cells that showed the most robust circadian rhythms. This choice possibly enhanced the effect of FBS concentration increase, not only on cell cycle phase duration times, but also on possible clock synchronising effects, which was the searched-for effect. In the present study, this is the only way we will take care of the circadian clock data. A more detailed analysis will be provided in a forthcoming study (cf Section 4).
The data processed in the identification procedure thus consisted of time series of intensities recording the red and green fluorescences emitted by individual NIH 3T3 cells proliferating within an in vitro homogeneous population, in a medium composed of a given FBS concentration. Only cells that were alive during the whole measurement period were tracked, so that our assumption of a zero death rate in the model is in accordance with these experimental conditions. The intensities were recorded every fifteen minutes, over 63 hours. A graph representing such a time series is presented on Figure 1.
We considered only data with at least the duration of a complete cell cycle, and measured the duration of the phases G 1 and S/G 2 /M within this cell cycle. The end of a cell cycle is characterised by the birth of two daughter cells that were also labelled and tracked. Consequently we measured the cell cycle duration as the time between cell birth and the FIGURE 1. Example of a time series of the intensity of red (deep grey, related to G 1 ) and green (light grey, related to S/G 2 /M ) fluorescences obtained by using the FUCCI method on a NIH 3T3 cell within a population in liquid medium.
division of this cell into two daughter cells. This method is equivalent to the one used in [8] since cell division is also characterised by a fast disappearance of the green fluorescence. During the transition from G 1 to S, red and green fluorescences overlap, so that it was not trivial to determine the duration of phase G 1 . As in [8], we decided to define the end of phase G 1 as the time at which red fluorescence was maximum before decreasing. The duration of phase S/G 2 /M was obtained by subtracting the duration of phase G 1 from the duration of the cell cycle. This method is summarised on Figure 2.

2.2.3.
Expression of the transition rates. Using the identification method presented above (Section 2.2.2), we identified 117 data for the duration of the phases G 1 and S/G 2 /M on individual cells cultured in a liquid medium containing 10% of foetal bovine serum (FBS), and 150 such data on individual cells cultured in a liquid medium containing 15% of FBS. The mean value of these durations and the corresponding standard deviation are given in Table 1 As expected, we could notice that the standard deviation of G 1 was higher than the one of S/G 2 /M , which reflects the known fact that the duration of G 1 is more variable than the one of S/G 2 /M in cell populations, a phenomenon that is common knowledge among biologists. More surprisingly, we could notice that increasing the concentration of FBS from 10% to 15% decreased not only the duration of G 1 but also the one of S/G 2 /M , leading to an overall higher growth rate of the population by a shortening of both phases.
We rounded each duration to the nearest hour. As in [8], the distributions of the durations of G 1 and of S/G 2 /M within the population were fitted to experimental data by using Gamma laws.
For all x ≥ 0, we thus used the following probability density functions: where Γ is the Gamma function, 1 [γi;+∞[ the indicator function of interval [γ i ; +∞[ and where the parameters α i , β i and γ i are given in Table 2. The corresponding curves are presented on Figure 3.   As in [8], we chose Gamma laws because they allowed a good (phenomenological) fit to our experimental data while keeping a reasonable number of parameters to be estimated. Moreover, there is a clear physiological basis to this choice: if the parameter α is an integer, the Gamma distribution is the law often used to represent probabilities of waiting times of the sum of α i.i.d. random variables representing waiting times, each one of them following an exponential law with the same parameter β. In our context, within G 1 or S/G 2 /M , these waiting times could be times between triggerings of crucial switches in a cascade of protein expressions leading to a phase transition, e.g., G 1 /S. Such an explanation, or parts of it, has been proposed, in the literature, in this context or others dealing with gene or protein expression, for instance in [13,32,43].
The mean, m, and the standard deviation, sd, of a random variable that follows a shifted Gamma law with parameters α, β,γ are given by the following formulas: The parameters we found for Gamma distributions led to a mean duration and a standard deviation on R + respectively of 9.0h and 3.1h for the G 1 phase and of 12.0h and 1.9h for the S/G 2 /M phase in the case of 10% of FBS, and of 7.7h and 1.9h for the G 1 phase and of 10.5h and 1.6h for the S/G 2 /M phase in the case of 15% of FBS. These figures, summarised in Table 3, are very close to the ones mentioned in Table 1 3. Mean (m) and standard deviation (sd) (in hours) of the Gamma distributed duration of the phases G 1 and S/G 2 /M for two experimental conditions (culture medium composed of 10% of FBS or of 15% of FBS), according to the parameters mentioned in Table 2.
As the experimental data were performed in vitro in a liquid medium, and as cells were kept in the medium three days before the recordings in order not to be synchronised by a "serum shock" (meaning a steep increase in growth factors contained in FBS induced by sudden serum adjunction) from the beginning of recordings, we could consider that there was no synchronisation between cells, hence no time dependency of the control of the growth process at the cell population level. These facts were consistent with our assumption of stationary transition coefficients, i.e., the assumption under which transition rates from G 1 to S/G 2 /M (K 1→2 ) and from S/G 2 /M to G 1 (K 2→1 ) did not depend on time, but only on the age of cells in the two phases. From Equation (6), we deduced the expression of the cumulative distribution function [8]: and thus we obtained: where f i represents the experimentally determined probability density function of the random variable representing the duration the cell spent in phase i. One may note here that the right hand side is usually known among probabilistic scientists as the hazard rate.
3.1. Discretisation scheme. We discretized our equations by means of finite differences as we did in [8]. As we introduced in the present paper a velocity v i in the model that is not necessarily equal to 1, we had to choose time and age steps, respectively denoted by ∆t and ∆x, so as to satisfy the CFL condition. Thus we assumed: 3.2. Internal validation. To make sure that our numerical results were in agreement with the biological data that we used to build our model ("internal validation"), we performed simulations in the case of no time control, that is, where the K i→i+1 (x) were given by the expression (10). As we only wanted, in this internal validation stage, to compare with biological data numerical results that were performed on the basis of modelling assumptions that corresponded to the two experimental conditions (10% and 15% of FBS), without assessing velocities, we assumed here that velocities v 1 and v 2 were constant, both equal to 1. The graphs of the transition rates we obtained from formula (10) and experimental data are presented on Figure 4.   Figure 5 presents the time evolution of the percentage of cells in phases G 1 and S/G 2 /M over the duration of one cell cycle resulting from numerical and biological experiments (biological data were preliminarily synchronised "by hand", i.e., by deciding that all cells were at age nought in phase G 1 at the beginning of simulations). We notice that modelled numerical data are very close to the raw biological data in the two experimental conditions. Therefore, we can conclude that the model and the method used to represent the proliferation phenomenon and fit our experimental data may have led us close to biological likelihood.
As mentioned in Section 2.1, the growth exponent λ, first eigenvalue of System (2), can be computed from Lotka's equation (cf Equation (3) and [12] for details). In the particular case of our experimentally based study described above, we had: With the coefficients of the two Gamma distributions identified from FUCCI data (see Table 2) where v 1 = v 2 = 1, this yielded λ ≈ 0.033h −1 for the experiment with a medium composed of 10% of FBS and λ ≈ 0.038h −1 for the experiment with a medium composed of 15% of FBS, which corresponds to doubling times (T d = ln (2) Table 4.  Figure 6.
As expected, the oscillations of the percentages of cells in G 1 and in S/G 2 /M were damped. These percentages rapidly reached a steady state. This phenomenon reflects the progressive desynchronisation of cells in the population (asynchronous cell growth [1,2,3]). We could notice that, even if the cell population submitted to 15% of FBS grew faster than the one submitted to 10% of FBS, desynchronisation occurred more rapidly in the case of 10% of FBS than in the case of 15% of FBS, i.e., it took longer for cells to desynchronise in the culture medium composed of 15% of FBS than in the medium composed of 10% of FBS. This corresponds to a lower phase duration variability in the case of 15% of FBS than in the case of 10% of FBS, lower variability in phase duration meaning sharper phase transition (cf. Section 2.2.3).
As remarked in Section 2.2.3 and in Table 4, an increased concentration in FBS (hence in growth factors) in the medium resulted in shorter durations of phases G 1 and S/G 2 /M and thus in faster cell proliferation. To model the effects of higher growth factor concentration in the medium on cell population dynamics, we considered the experimental data for a concentration of 10% of FBS in the medium as reference data. Our aim was to recover the population dynamics for a concentration of 15% of FBS in the medium using this reference population as a basis. We hal-00730457, version 1 -10 Sep 2012 CELL POPULATION MODEL WITH GROWTH FACTORS 13 assumed that FBS influenced only cell ageing velocity and we wanted to check how this assumption was consistent with our biological data.
We thus calculated the probability density function deduced from the distribution of the durations of phases G 1 and S/G 2 /M fitted in the population of cells proliferating in a medium composed of 10% of FBS, using a shifted Gamma distribution model, whose parameters were found to be α 1 = 1.80, β 1 = 0.43h −1 , γ 1 = 4.83h, α 2 = 16.96, β 2 = 2.22h −1 , γ 2 = 4.37h. From these figures we deduced the transition rates K i→i+1,10% under the assumption v 1 = v 2 = 1 (cf Section 2.2.3) and subsequently used these parameters of the age-structured model, except that v 1 and v 2 were no longer fixed, to identify these velocities in the case of 15% of FBS in the medium.
For convenience in this preliminary study, we assumed that v 1 = v 2 = v, i.e., that 15% FBS cells were ageing with the same, constant, velocity in G 1 and in S/G 2 /M . This velocity v was then viewed as a parameter to fit the percentages of cells in G 1 and S/G 2 /M resulting from the biological experiment with a medium composed of 15% of FBS (and synchronised "by hand" as in Section 3.2). Thus the model we used was: where K i→i+1,10% corresponded to the transition rates determined in the case of 10% FBS (cf Section 2.2.3).
Our aim was thus to minimise the mean squared error that quantified the discrepancy between the percentages of cells in G 1 and in S/G 2 /M computed by Equations (13) and the experimental ones. We thus used the mean squared method by means of the fmincon Matlab function with the velocity v as parameter (the constraint being the positivity of v). This method led to a computed value of v equal to 1.095. The corresponding computed and experimental percentages of cells in G 1 and S/G 2 /M are presented on Figure 7. We can notice that modelled numerical data were very close to raw biological data.
This value v = 1.095 led to an exponential growth exponent λ (computed by means of Equation (3)) equal to 0.045h −1 and so to a doubling time T d equal to 15.4h. These figures tend to demonstrate that cells would proliferate approximately 10% faster in a medium composed of 15% of FBS than in a medium composed of 10% of FBS. 4. Conclusion and discussion. Growth factors, such as those contained in foetal bovine serum (FBS) are known to influence cell cycle dynamics. A lot of studies were and are still performed by scientists aiming at fully understanding the underlying mechanisms, see for instance [14,16,25,29,30,31,45]. We proposed an age-structured mathematical model whose parameters were computed on the basis of biological data, and that enabled us to analyse the cell synchronisation in the cell cycle according to the FBS concentration. This model also enabled us to recover the effects of an increase in the FBS concentration in the culture medium.
In our biological data, FBS decreased both the durations of G 1 and S/G 2 /M . These results are in agreement with the results presented in the literature, for instance in [16,31]. In [16], the authors studied the effects of several concentrations of epidermal growth factor (EGF) on the proliferation and cell cycle regulation of cultured human amnion epithelial cells, using flow cytometry and gene expression measurements. They concluded that EGF increase resulted in reduced expression of the cell cycle control genes, which explained the increased concentration of cells in S and G 2 /M phases they observed when using EGF supplementation.
The value of the predicted exponential growth exponent λ for the case of 15% FBS (λ = 0.045h −1 ) was higher than the one we computed on the basis of the raw experimental data (λ = 0.038h −1 ). Nevertheless it seemed to us that the experimental data were "visually" better fitted by the predictive model (Equations (13)) than by the experimentally-based model (Equations (1)). This difference could be due to the fact that we directly fitted results of the model (percentages of cell in G 1 and S/G 2 /M ) instead of fitting it indirectly, i.e., through the Gamma distributions. Another way to improve the predictive character of our model would be to recur to another minimisation method (e.g., CMAES [22]).
Moreover, in this preliminary study, we assumed that FBS only impacted the cell ageing velocities, in the same manner in the phase G 1 as in the phase S/G 2 /M . In fact, the mechanisms leading to changes in the cell cycle dynamics due to an increase (or decrease) of the FBS concentration in the culture medium might be much more complex.
Indeed, one may naturally question the way in which we investigated the effect of an increase in growth factors on the cell cycle model parameters: we focused on the velocity parameters v 1 and v 2 with which phases G 1 and S/G 2 /M are cruised by proliferating cell populations in the division cycle, which implicitly assumes that protein synthesis mechanisms in phases G 1 and G 2 , the major biological phenomena occurring in these phases, are the main targets of growth factors, and this may be an oversimplified modelling assumption. Indeed, as mentioned in the introduction, it is known that growth factors affect in particular Cyclin D and its inhibitors at the restriction point of phase G 1 , which could be differently represented in the model, namely by a direct action on phase transition rates K i→i+1 . This has been done only graphically (Figure 4) in this paper, where the choice made here was to focus on velocity parameters v 1 and v 2 .
We based our study on two FBS concentrations. Biological FUCCI data performed on other cell types (healthy and/or cancer cells) with other varying FBS concentrations in the culture medium would help us to further analyse the effects of FBS on the cell cycle dynamics and thus to improve the model predictability.
In a forthcoming work, we plan to introduce time control in the model presented in this paper in order to analyse the effect of increasing the FBS concentration of the culture medium on cell circadian clocks. Such an analysis could help us to better understand the role of the circadian clock in cell cycle regulation of healthy and cancer cells in the presence of growth factor receptor (GFR) inhibitors. This should prove helpful to design theoretically optimised cancer treatments, that combine cytotoxic drugs, that exert their main action by arresting the cell cycle at phase transitions via DNA damage, ATM triggering and p53 control, and GFR blockers. Such therapeutic combinations, currently in use in the clinic ( [34], see also [9] and references therein), should indeed benefit from better knowledge of interactions between growth factors and the cell cycle from a modelling point of view.