Pierre Gabriel
email: gabriel@ann.jussieu.fr

Shawn P Garbett

Darren R Tyson

Glenn F Webb
email: glenn.f.webb@vanderbilt.edu

The contribution of age structure to cell population responses to targeted therapeutics

Keywords: Cell cycle, intermitotic time, renewal equation, exponentially modified gaussian

Cells grown in culture act as a model system for analyzing the effects of anticancer compounds, which may affect cell behavior in a cell cycle position-dependent manner. Cell synchronization techniques have been generally employed to minimize the variation in cell cycle position. However, synchronization techniques are cumbersome and imprecise and the agents used to synchronize the cells potentially have other unknown effects on the cells. An alternative approach is to determine the age structure in the population and account for the cell cycle positional effects post hoc.

Here we provide a formalism to use quantifiable lifespans from live cell microscopy experiments to parameterize an age-structured model of cell population response.

Introduction

When examined individually in time lapse microscopy experiments, cells grown in culture display variability in the length of their cell cycles (between mitotic events), and this variability is represented by intermitotic time (IMT) distributions [START_REF] Dowling | Modelling cell lifespan and proliferation: is likelihood to die or to divide independent of age?[END_REF][START_REF] Lee | Modeling t cell proliferation and death in vitro based on labeling data: Generalizations of the smithmartin cell cycle model[END_REF][START_REF] Liu | A randomperiods model for expression of cell-cycle genes[END_REF]. These distributions are usually obtained from asynchronously dividing populations of cells, which achieve a steady-state age structure when the population is growing exponentially. In experimental studies that examine the effects of perturbations on cellular proliferation, it is often desirable to determine whether the perturbation is affecting cells in particular stages within their cell cycle, i.e. whether the perturbation has cell cycle-specific effects. Since we do not know where a cell is in the cycle, an alternative is to define the age of a cell as the time elapsed from its last division and deal with this measurable quantity instead of the cell cycle phase. In this paper we adopt this definition and provide a formalism to convert data obtained as IMT distributions to parameterize an age-structured population model and, thus, identifying the contribution of age structure to the response of the cell population to perturbation. Models of age-structured populations using partial differential equations, such as those originally discribed by A. Lotka, A.G. McKendrick, W.O. Kermack, and F. von Foerster, are well adapted to model the dynamical features of experimental cultures of cells transiting the cell cycle with variable IMT. These models have been widely studied from a mathematical perspective [START_REF] Ainseba | An innovative multistage, physiologically structured, population model to understand the european grapevine moth dynamics[END_REF][START_REF] Ainseba | Cml dynamics: Optimal control of age-structured stem cell population[END_REF][START_REF] Calsina | A model of physiologically structured population dynamics with a nonlinear individual growth rate[END_REF][START_REF] Clément | One-parameter semigroups[END_REF][START_REF] Diekmann | On the stability of the cell size distribution[END_REF][START_REF] Engel | One-parameter semigroups for linear evolution equations[END_REF][START_REF] Metz | The dynamics of physiologically structured populations[END_REF][START_REF] Perthame | Transport equations in biology[END_REF][START_REF] Picart | Optimal control problem on insect pest populations[END_REF][START_REF] Webb | Theory of nonlinear age-dependent population dynamics[END_REF], but the application of the model to experimental data has been hampered by an inability to determine the age-dependent model parameters. The usual approach for parameter estimation is to solve numerically an inverse problem (see [START_REF] Ackleh | Parameter identification in size-structured population models with nonlinear individual rates[END_REF][START_REF] Banks | Estimation of dynamic rate parameters in insect populations undergoing sublethal exposure to pesticides[END_REF][START_REF] Banks | A new model for the estimation of cell proliferation dynamics using cfse data[END_REF][START_REF] Doumic | Numerical solution of an inverse problem in sizestructured population dynamics[END_REF][START_REF] Doumic | On the calibration of a size-structured population model from experimental data[END_REF][START_REF] Engl | A regularization scheme for an inverse problem in agestructured populations[END_REF][START_REF] Gyllenberg | The inverse problem of linear age-structured population dynamics[END_REF][START_REF] Picart | Parameter identification in multistage population dynamics model[END_REF][START_REF] Pilant | Determining a coefficient in a first-order hyperbolic equation[END_REF][START_REF] Pilant | Determining the initial age distribution for an age structured population[END_REF][START_REF] Rundell | Determining the birth function for an age structured population[END_REF][START_REF] Rundell | Determining the death rate for an age-structured population from census data[END_REF][START_REF] Perthame | On the inverse problem for a size-structured population model[END_REF] on this question for structured population models), but this requires extensive input data and is specific to a given situation. A much more convenient approach is to assume that the distributed parameters lie in a class of functions with only a few constants (power laws or other forms [START_REF] Gabriel | The shape of the polymerization rate in the prion equation[END_REF]) and obtain the parameters from the fit of the functions to experimental data.

One function that often provides a better fit to IMT data than others, such as log-normal, inverse normal or gamma functions, is an exponentially modified Gaussian (EMG) [START_REF] Golubev | Exponentially modified Gaussian (EMG) relevance to distributions related to cell proliferation and differentiation[END_REF][START_REF] Tyson | A method to quantify population-level cell proliferation dynamics from single-cell data[END_REF]. Under conditions in which the IMT distribution can be explained by an EMG model, we submit that the age-dependent division rate can be identified as an error function. Starting from this observation, we present a simple method to recover the parameters of this error function from parameters fitted to the experimental IMT data. Once reliably parameterized, the age-structured model can be used to make predictions about cell age-dependent effects of perturbations, for example, whether cells arrest during their cell cycle in response to treatment with antiproliferative compounds.

Individual based models (IBM) have also been used to track individual cell behavior in proliferating cell populations [START_REF] Alarcón | A cellular automaton model for tumour growth in inhomogeneous environment[END_REF][START_REF] Mallet | A cellular automata model of tumorimmune system interactions[END_REF][START_REF] Gerlee | A hybrid cellular automaton model of clonal evolution in cancer: The emergence of the glycolytic phenotype[END_REF][START_REF] Patel | A cellular automaton model of early tumor growth and invasion: The effects of native tissue vascularity and increased anaerobic tumor metabolism[END_REF][START_REF] Powathil | Modelling the effects of cell-cycle heterogeneity on the response of a solid tumour to chemotherapy: Biological insights from a hybrid multiscale cellular automaton model[END_REF][START_REF] Quaranta | Invasion emerges from cancer cell adaptation to competitive microenvironments: Quantitative predictions from multiscale mathematical models[END_REF][START_REF] Quaranta | Mathematical modeling of cancer: The future of prognosis and treatment[END_REF][START_REF] Ramis-Conde | Modeling the influence of the e-cadherin--catenin pathway in cancer cell invasion: A multiscale approach[END_REF][START_REF] Ribba | The use of hybrid cellular automaton models for improving cancer therapy[END_REF]. IBM have advantages in their direct connection to discrete events and to cell-cell interactions. IBM are usually readily implementable to computer simulations, but may require many repeated runs to access an average outcome. IBM are usually difficult to analyze theoretically with respect to parametric input, particularly with highly sensitive parameters. The advantage of continuum differential equations models, such as the cell age structured models we develop here, are their tractability for theoretical analysis and their reproducibility for specific parametric input. The practical difficulties of determining this parametric input is the problem we focus upon here. IBM and differential equations models should be viewed as companion approaches, which complement and compare their insights and predictions for proliferating cell population behavior, when individual cell variation is a primary consideration.

The age-structured model

As with human populations, we can associate an age to each individual cell in a cell population. We define the age of a cell as the time elapsed from its last mitosis. Starting from this definition, we derive an age-structured model giving the dynamics of populations of proliferating and quiescent cells. This model is adapted to investigate the effect of treatment by the drug erlotinib on in vitro PC-9 cancer cell lines. It has been recently shown in [START_REF] Tyson | A method to quantify population-level cell proliferation dynamics from single-cell data[END_REF] that the main effect of erlotinib on cancer cells is to induce entry into quiescence. In [START_REF] Tyson | A method to quantify population-level cell proliferation dynamics from single-cell data[END_REF] a system of ordinary differential equations (without age structure) is used to model these experiments. We hypothesize here that erlotinib induced quiescence is linked to the age of the cells involved. Many age-structured models with quiescence can be found in the literature (see [START_REF] Arino | Necessary and sufficient conditions for asynchronous exponential growth in age structured cell populations with quiescence[END_REF][START_REF] Brikci | Analysis of a molecular structured population model with possible polynomial growth for the cell division cycle[END_REF][START_REF] Brikci | An age-and-cyclin-structured cell population model for healthy and tumoral tissues[END_REF][START_REF] Gyllenberg | Age-size structure in populations with quiescence[END_REF][START_REF] Gyllenberg | A nonlinear structured population model of tumor growth with quiescence[END_REF][START_REF] Gyllenberg | Quiescence in structured population dynamics: applications to tumor growth[END_REF] for examples). Here we present a simple age-structured model with quiescence which allows us to explain observed delays in response to erlotinib. We start from the observation that there is no effect of the treatment on total population growth during the first twenty hours (see [START_REF] Tyson | A method to quantify population-level cell proliferation dynamics from single-cell data[END_REF] and Figure 7). This time corresponds almost exactly to the minimal age of division observed for PC-9 cells. It suggests that erlotinib acts only during a specific phase of the cell cycle, which based on its biological activity would be expected to be in G1. The model we present is based on this idea and considers a fractional rate f of cells that become quiescent in an age-dependent manner, where the fraction is assumed to reflect the dose of erlotinib used for the treatment.

We start from the the McKendrick-Von Foerster's model which is widely used to model the cell cycle (see [START_REF] Metz | The dynamics of physiologically structured populations[END_REF][START_REF] Perthame | Transport equations in biology[END_REF][START_REF] Webb | Theory of nonlinear age-dependent population dynamics[END_REF] and references). The partial differential equation in this model provides the evolution of the density p(t, a) of cells with age, or "cell cycle phase", a at time t. To take into account the quiescent cells, we introduce the quantity Q(t) which representsthe density of quiescent cells at time t. It evolves according to an ordinary differential equation coupled to the equation on p(t, a). We obtain the system

                               ∂ ∂t p(t, a) + ∂ ∂a p(t, a) + β(a)p(t, a) + µ p(t, a) = 0, t ≥ 0, a > 0, p(t, 0) = 2(1 -f ) ∞ 0 β(a)p(t, a) da, p(0, a) = p 0 (a), d dt Q(t) = 2f ∞ 0 β(a)p(t, a) da -µ Q(t), Q(0) = Q 0 . (1)
In this model, the proliferaing cells age one-to-one with time at speed da dt = 1, and divide with rate β(a) ≥ 0. When a cell divides at mitosis, it produces two daughters which are either proliferating with age a = 0, or quiescent. This is taken into account by the boundary condition at a = 0 and the first term in the equation for Q(t), where the parameter f ∈ [0, 1] represents the fraction of proliferating cells which become quiescent. The coefficient µ is a death rates. The number of proliferating cells at time t with age between a 1 and a 2 is a 2 a 1 p(t, a) da, and the total number of cells at time t is

N (t) = ∞ 0 p(t, a) da + Q(t).
The division rate β has a probabilistic interpretation: the probability that a cell did not divide by age a is given by

P(a) = e -a 0 β(a ′ ) da ′ .
Since all the proliferating cells must divide at some time by definition, the division rate has to satisfy

lim a→+∞ a 0 β(a ′ ) da ′ = +∞. ( 2
)
The model is completed with intial data p 0 (a) and Q 0 . We choose the time t = 0 to be the beginning of the erlotinib treatment, so at this time there are no quiescent cells (Q 0 = 0). The age distribution of the proliferating cells is assumed to be at equilibrium (i.e. p 0 (a) = const p(a), see Section 3 for the mathematical definition of p). The experimental values of the total population N (t) along time are ploted after normalization by the initial value N (0) on a log-scale (see Figure 7). Because of this normalization, we consider an initial distribution such that ∞ 0 p 0 (a) da = 1 which leads to p 0 (a) = p(a) because of the definition of p.

We want to compare the solutions of model (1) to the experimental observations presented in [START_REF] Tyson | A method to quantify population-level cell proliferation dynamics from single-cell data[END_REF]. The first step is to estimate the different parameters of the model. To estimate the value of the fraction f for different doses of treatment, we use data available in [START_REF] Tyson | A method to quantify population-level cell proliferation dynamics from single-cell data[END_REF]. The fraction F of quiescent cells is estimated by examination of whether cells treated with erlotinib divide or not before the end of the experiment. We want to use this experimental fraction F to estimate the coefficient f of the model. For the sake of simplicity, assume that the death rate µ is 0 in [START_REF] Ackleh | Parameter identification in size-structured population models with nonlinear individual rates[END_REF]. In this case, at the end of the labeling period t 0 , the quantity of labeled quiescent cells corresponds to Q(t 0 ) and the quantity of proliferating cells corresponds to t 0 0 p(t, 0) dt. Thus, the fraction F of quiescent cells at t 0 is

F = Q(t 0 ) Q(t 0 ) + t 0 0 p(t, 0) dt . (3)
Now we compute this quantity from model [START_REF] Ackleh | Parameter identification in size-structured population models with nonlinear individual rates[END_REF], keeping in mind that we have assumed no mortality. We have

Q(t 0 ) = t 0 0 dQ dt (t) dt = 2f t 0 0 ∞ 0 β(a)p(t, a) da dt and t 0 0 p(t, 0) dt = 2(1 -f ) t 0 0 ∞ 0 β(a)p(t, a) da dt.
Finally we obtain

Q(t 0 ) Q(t 0 ) + t 0 0 p(t, 0) dt = 2f 2f + 2(1 -f ) = f. (4)
So when there is no death rate, the experimental fraction F corresponds exactly to the coefficient f.

In our simulations we consider positive a death rate µ as suggested by the results in Section 4.2. But because the numerical value we recover for µ is very small (µ ≪ 1 in Figure 6), we can consider that f is still well-approximated by F and we use the fractions available in [START_REF] Tyson | A method to quantify population-level cell proliferation dynamics from single-cell data[END_REF] to parameterize f. The identification of the coefficients β and µ is much more delicate because β(a) is a distributed function. The two following sections are devoted to presenting a method to recover these coefficients from IMT distributions. Then these coefficients are used to compare the solutions of Equation ( 1) to experimental data in Section 5.

Modelling the intermitotic time

Since the age of a cell is defined as the time elapsed from its last mitosis, the IMT of a cell is its age at division. This definition allows us to interpret the IMT distributions in terms of a dynamic age-structured population model.

Experimental IMT distributions can be seen as histograms which represent, for a given population, the density of cells with a certain age of division (see Figure 1 for an example taken from [START_REF] Tyson | A method to quantify population-level cell proliferation dynamics from single-cell data[END_REF]). The age of division is distributed into N a bins of width ∆a. For all i between 1 and N a , the height H i of the i th bar represents the density of cells with an age of division in the window [i∆a, (i + 1)∆a]. The histogram (H i ) 1≤i≤Na is normalized to represent a density ∆a Na i=1

H i = 1.
(5) [START_REF] Tyson | A method to quantify population-level cell proliferation dynamics from single-cell data[END_REF] and the bin width is ∆a = 10 8 .

We briefly explain here the method used in [START_REF] Tyson | A method to quantify population-level cell proliferation dynamics from single-cell data[END_REF] to build IMT histograms (we refer to this paper for more details). The data are obtained using extended temporally resolved automated microscopy (ETRAM) in which cell nuclei are fluorescently labeled, imaged by automated time lapse fluorescence microscopy, and tracked as individual cells from the resultant image stacks. Cells are subjected to various microenvironmental conditions (such as the addition of a drug) at a specific time during image acquisition (set to time t 0 ) and the effect of the perturbation from that point in time is followed across the entire population and individual cells within it. The duration of observation (T ) is chosen large enough to observe that almost all cycling cells divide before this final time. The data are organized into the bins to obtain a histogram, which is then normalized to represent a density of cells. The next step is to describe the IMT distributions in terms of an age-structured model and use this interpretation to parameterize Equation [START_REF] Ackleh | Parameter identification in size-structured population models with nonlinear individual rates[END_REF]. We want to estimate the coefficients β and µ, and for this we fix f = 0 in [START_REF] Ackleh | Parameter identification in size-structured population models with nonlinear individual rates[END_REF]. This corresponds to the situation when there is no erlotinib treatment and thus all the cells are proliferating. We recover in this case the original McKendrick-Von Foerster's model

               ∂ ∂t p(t, a) + ∂ ∂a p(t, a) + β(a)p(t, a) + µ p(t, a) = 0, t ≥ 0, a > 0, p(t, 0) = 2 ∞ 0 β(a)p(t, a) da, p(0, a) = p 0 (a). (6)
       λp(a) + ∂ a p(a) + β(a)p(a) + µ p(a) = 0, p(0) = 2 ∞ 0 β(a)p(a) da, p(•) > 0, p(a) da = 1, (7)
This problem has a unique solution given by

p(a) = p(0)e -a 0 (β(a ′ )+µ+λ) da ′
where λ > 0 is the unique value such that

1 = 2 ∞ 0 β(a)e -a 0 (β(a ′ )+µ+λ) da ′ da (8)
and

p(0) = ∞ 0 e -a 0 (β(a ′ )+µ+λ) da ′ da -1
.

Then we can prove that, for large times, p(t, a) ∼ const p(a)e λt (see Appendix A for more details and references). If the population of cells proliferates over a sufficiently long time, we can assume that this asymptotic behavior is reached and use it to investigate the IMT distributions. An experimental observation that the total population (independent of age structure) is growing exponentially is an indicator that the population has effectively reached the equilibrium age distribution, which can be checked using ETRAM for other time series data collection (see Figure 2). Now we give a continuous expression of the IMT distribution in terms of the age-structured model. The age distribution of the cells relative to time t 0 (time of perturbation) is given by a truncation of the equilibrium age distribution p t (in hour)

p0 (a) = ρ e -a 0 (β(a ′ )+µ+λ) da ′ if 0 ≤ a ≤ t 0 , 0 if a > t 0 .
Figure 2: The evolution of the total population N (t) is plotted on log-scale (blue) (the data is from [START_REF] Tyson | A method to quantify population-level cell proliferation dynamics from single-cell data[END_REF]). The data are well-fitted with a line of slope λ = 0.022 (correlation coefficient: R 2 = 0.9967). This is evidence that the total population is growing exponentially fast with exponential constant λ. Thus, N (t) = exp λt N (0), and if t ⋆ is the population doubling time, then N (t ⋆ ) = 2N (0), and

t ⋆ = ln 2/λ.
where ρ is a scaling constant. We then follow this age distribution along time and obtain for t > 0 p(t, a) = ρ e -a 0 (β(a ′ )+µ+λ) da ′ e λt if t ≤ a ≤ t + t 0 , 0 otherwise.

According to the age-structured model, and because of the normalization (5), the IMT distribution satisfies

I T (a) := C -1 T T 0 β(a)p(t, a) dt, C T := ∞ 0 T 0 β(a)p(t, a) dt da, (9)
by definition of the division rate β. The fact that no cell can divide in a time less than t 0 means

∀ a ≤ t 0 , β(a) = 0. ( 10
)
Under this condition, for T large, the function I T is close to the function

I ∞ (a) := C -1 ∞ β(a)e -a 0 (β(a ′ )+µ) da ′ (11)
where

C ∞ := ∞ 0 β(a)e -a 0 (β(a ′ )+µ) da ′ da. ( 12
)
This convergence is made mathematically precise by the following claim (the proof is given in Appendix B for a more general case in which it is not assumed that p(t, a) is close to the equilibrium age distribution):

Claim 1. Under Assumptions (2) and (10), we have the convergence

∞ 0 |I T (a) -I ∞ (a)| da ----→ T →∞ 0.
Remark 2. We can easily prove that, under the additional assumption

lim a→+∞ β(a)e -a 0 β(a ′ ) da ′ = 0, ( 13
)
we also have

sup a≥0 |I T (a) -I ∞ (a)| -----→ T →+∞ 0.
Condition (13) is satisfied from Assumption (2), if for example, β bounded or monotonic.

Because in the experimental protocol T is chosen large enough to observe no dividing cells at the end, we can approximate I T by I ∞ which has a simple expression in terms of the division and death rates [START_REF] Brikci | Analysis of a molecular structured population model with possible polynomial growth for the cell division cycle[END_REF].

From the intermitotic time to the division rate

In this section we explain how we can recover the parameters β and µ from the IMT distribution I ∞ . We first present the method in the case when there is no death (µ = 0). This simplification is useful to give an inversion formula which gives the rate β in terms of I ∞ . Then we extend the method to include possible death rates.

The case µ = 0

In the case µ = 0, the constant C ∞ is equal to 1 and we have

I ∞ (a) = β(a)e -a 0 β(a ′ ) da ′ .
This expression can be inverted to recover the division rate from the IMT distribution (see [START_REF] Hinow | A mathematical model separates quantitatively the cytostatic and cytotoxic effects of a HER2 tyrosine kinase inhibitor[END_REF][START_REF] Billy | Synchronisation and control of proliferation in cycling cell population models with age structure[END_REF])

β(a) = I ∞ (a) ∞ a I ∞ (a ′ ) da ′ . (14)
We start from the fitting of the experimental IMT distributions, which are observed to be positively skewed. From a fitting procedure for I ∞ , we use Equation ( 14) to recover the division rate of the McKendrick-Von Foerster's equation.

First consider as in [START_REF] Hinow | A mathematical model separates quantitatively the cytostatic and cytotoxic effects of a HER2 tyrosine kinase inhibitor[END_REF] that the IMT distribution is a shifted gamma function (see also [START_REF] Bernard | Analysis of cell kinetics using a cell division marker: Mathematical modeling of experimental data[END_REF] and references therein). Setting

I ∞ (a|m, σ) = 0 if 0 ≤ a ≤ m, a -m σ 2 e -a-m σ if a ≥ m, (15)
we can solve explicitly Equation ( 14) and we find

β(a) =    0 if 0 ≤ a ≤ m, a -m σ(σ + a -m) if a ≥ m. (16)
So by fitting an experimental IMT distribution with a shifted gamma function, we obtain two parameters m and σ which allows reconstruction of the division rate of the renewal equation.

To have a smoother transition at the minimum age of division m, one can consider a second shifted gamma function

I ∞ (a|m, σ) =    0 if 0 ≤ a ≤ m, (a -m) 2 2σ 3 e -a-m σ if a ≥ m. (17)
Then the corresponding β is

β(a) =    0 if 0 ≤ a ≤ m, 1 σ (a -m) 2 2σ 2 + 2σ(a -m) + (a -m) 2 if a ≥ m. ( 18
)
The different functions [START_REF] Calsina | A model of physiologically structured population dynamics with a nonlinear individual growth rate[END_REF] to [START_REF] Douglas | Numerical methods for a model of population dynamics[END_REF] are plotted in Figure 3 It has been observed that an exponentially modified Gaussian (EMG) is often a better model for IMT distributions than the gamma function [START_REF] Golubev | Exponentially modified Gaussian (EMG) relevance to distributions related to cell proliferation and differentiation[END_REF][START_REF] Tyson | A method to quantify population-level cell proliferation dynamics from single-cell data[END_REF] . An EMG is defined as the convolution of a Gaussian with a decreasing exponential, but after solving it can be written with three parameters as

I ∞ (a|β 0 , m, σ) = β 0 Erf c m -a σ e -2β 0 β 0 σ 2 2 -m+a . ( 19
)
where the (complementary) error function is defined by

Erf c(z) = 1 - 2 √ π z 0 e -t 2 dt.
Replacing I ∞ by an EMG in Equation ( 14), we cannot compute explicitly the expression for β. But by numerical comparison, we obtain a division rate β that is essentially indistinguishable from an error function (see Figure 4). 14) by using an EMG with coefficients m = 22, σ = 2 and β 0 = 0.2 to fit I ∞ [START_REF] Doumic | On the calibration of a size-structured population model from experimental data[END_REF]. Right: The formula for β, as an error function with these same three parameters, is numerically indistinguishable from the numerically obtained β (R 2 = 1).

Instead of fitting the IMT distribution by an EMG and then fitting the corresponding β by an error function, we may directly assume that β is an error function

β(a) = β 0 Erf c m -a σ . ( 20
)
We can then explicitly derive a new fitting formula for I ∞ due to Equation ( 11)

I ∞ (a|β 0 , m, σ) = β 0 Erf c m -a σ e -a 0 β 0 Erf c( m-a ′ σ ) da ′ ( 21
)
where the integral in the exponential can be computed as

a 0 Erf c m -a ′ σ da ′ = m Erf c m σ - σ √ π e -( m σ ) 2 -(m -a)Erf c m -a σ + σ √ π e -( m-a σ ) 2 . ( 22
)
Using this formula for the IMT instead of an EMG formula, the fitting parameters provide immediately the division rate β.

When we know the Malthusian growth parameter (population growth rate λ) from experimental data (see Figure 2), it is possible to recover a division rate β such that relation [START_REF] Arino | Necessary and sufficient conditions for asynchronous exponential growth in age structured cell populations with quiescence[END_REF] is satisfied. This is an important step in order to parameterize and apply an age-structured model to experimental data. We notice the relationship of β in terms of

I ∞ (a) is 2 ∞ 0 I ∞ (a)e -λa da = 2 ∞ 0 β(a)e -a 0 β(a ′ ) da ′ e -λa da = 1.
To take advantage of this relationship between the IMT data, β(a), and the malthusian parameter λ at the equilibrium distribution, we define a new histogram ( Hi ) 1≤i≤N by

∀ i, Hi := 2H i e -λa i ∆a N i=1 2H i e -λa i (23)
where a i := (i + 1 2 )∆a is the mean age of the i th bar. Thus defined, the new histogram incorporates information about λ and satisfies the relation ∆a

∞ i=1 Hi = 1. ( 24
)
We then fit this new histogram ( Hi ) with the model

Ĩ∞ (a) := 2I ∞ (a)e -λa (25)
instead of fitting (H i ) with I ∞ (a). Because of the normalization ( 24), we expect that the fitting provides parameters such that ∞ 0 Ĩ∞ (a) da ≈ 1, and this relation can be checked numerically a posteriori.

Method and example. We divide the process in three steps and illustrate it by an example. Each required fitting step can be performed using the freely available Ezyfit Matlab toolbox [< http://www.fast.u-psud.fr/ezyfit/ >].

Step 1: Determine equilibrium IMT distribution a) Obtain a histogram (H i ) for the experimental IMT distribution of control (untreated) cells.

Here we use an IMT distribution obtained in [START_REF] Tyson | A method to quantify population-level cell proliferation dynamics from single-cell data[END_REF] for PC-9 lung cancer cells using ETRAM (see Figure 1).

b) Plot the time evolution of the total population on a log e -scale, verify a linear fit and obtain the slope as the experimental value for the Malthusian parameter λ of the cell population (see Figure 2 for an example).

c) Construct the new histogram ( Hi ) from (H i ) and λ by using the definition [START_REF] Engl | A regularization scheme for an inverse problem in agestructured populations[END_REF].

Step 2: Obtain parameters from model fit to IMT distribution a) Choose a form for I ∞ as a gamma function [START_REF] Calsina | A model of physiologically structured population dynamics with a nonlinear individual growth rate[END_REF] or [START_REF] Diekmann | On the stability of the cell size distribution[END_REF], or as the new EMG form [START_REF] Dowling | Modelling cell lifespan and proliferation: is likelihood to die or to divide independent of age?[END_REF].

b) Fit the histogram ( Hi ) with the corresponding form Ĩ∞ from definition [START_REF] Gabriel | High-order WENO scheme for polymerization-type equations[END_REF]. For PC-9 cancer cells we choose the form [START_REF] Dowling | Modelling cell lifespan and proliferation: is likelihood to die or to divide independent of age?[END_REF], because in [START_REF] Tyson | A method to quantify population-level cell proliferation dynamics from single-cell data[END_REF] it was observed that the IMT distribution appeared to be an EMG (see Figure 5 for the example).

c) Verify that the correlation coefficient R 2 of the IMT data and the chosen form Ĩ∞ (a) is close to 1.

Step 3: Parameterize age-structured model

If numerical integral

∞ 0 Ĩ∞ (a) da is close to 1, the fitting parameters provide a good approximation of the division rate β, which can then be propagated through the population, based on the choice of I ∞ . In the example, β is given by Equation (20) with the numerical paramaters of Figure 5.

Introducing a death rate

If we include a death rate µ > 0 in the model, then there does not exist an inversion formula as Equation ( 14) to recover β analytically. So we consider that β is a rational fraction given by Equation ( 16) or [START_REF] Douglas | Numerical methods for a model of population dynamics[END_REF], or an error function given by Equation ( 20), and we derive the corresponding model I ∞ from Equation [START_REF] Brikci | Analysis of a molecular structured population model with possible polynomial growth for the cell division cycle[END_REF]. The death rate µ will appear in the expression of I ∞ as an additional fitting parameter and will be determined together with the parameters of β by the fitting of experimental IMT distributions. Notice that a fixed rate µ based on measured data can also be used if available to avoid the addition of a free parameter.

Unlike the case µ = 0 for which C ∞ = 1, we have for a positive death rate that C ∞ < 1. This constant is a function of the fitting parameters (see Equation ( 12)), but we do not have an analytic expression of this function in general. A solution to this problem is to use the Malthusian parameter λ, which incorporates information about µ and β(a), to define a form Ĩ∞ (a) as in Section 4.1. Indeed we obtain, because of relation [START_REF] Arino | Necessary and sufficient conditions for asynchronous exponential growth in age structured cell populations with quiescence[END_REF], Ĩ∞ (a) := 2β(a)e -a 0 (β(a ′ )+µ) da ′ e -λa which does not involve C ∞ . Then, the division rate β(a) is obtained from Ĩ∞ (a) as before using the modified histogram ( Hi ) defined from (H i ) by [START_REF] Engl | A regularization scheme for an inverse problem in agestructured populations[END_REF].

Example. We fit the same distribution as in Section 4.1 still considering that β is an error function.

With a constant death rate µ, we obtain the four parameters model

Ĩ∞ (a|β 0 , m, σ, µ) = 2β 0 Erf c m -a σ e -a 0 (β 0 Erf c ( m-a ′ σ )+µ+λ) da ′ (26)
where the integral a 0 β 0 Erf c m-a ′ σ da ′ is given by Equation [START_REF] Engel | One-parameter semigroups for linear evolution equations[END_REF]. The fitting provides new parameters for β and a positive death rate µ (see Figure 6). The correlation R 2 is slightly better than in Figure 5 and the integral ∞ 0 Ĩ∞ (a) da is significantly closer to 1. So we can assume that mortality has to be considered for this cell line.

Numerical simulation and discussion

Once the parameters accurately estimated, we use them to compare the behavior of the solution of Equation ( 1) to the experimental data in [START_REF] Tyson | A method to quantify population-level cell proliferation dynamics from single-cell data[END_REF]. Since the model does not have analytic solutions, we need to perform numerical simulations. Many numerical methods are available for structured population equations (see for instance [START_REF] Ackleh | An implicit finite difference scheme for the nonlinear size-structured population model[END_REF][START_REF] Angulo | Numerical integration of nonlinear size-structured population equations[END_REF][START_REF] Angulo | Numerical integration of fully nonlinear size-structured population models[END_REF][START_REF] Douglas | Numerical methods for a model of population dynamics[END_REF][START_REF] Gabriel | High-order WENO scheme for polymerization-type equations[END_REF][START_REF] Kostova | An explicit third-order numerical method for size-structured population equations[END_REF]). Here we use a scheme based on the method of characteristics, as in [START_REF] Angulo | Numerical integration of nonlinear size-structured population equations[END_REF][START_REF] Angulo | Numerical integration of fully nonlinear size-structured population models[END_REF][START_REF] Kostova | An explicit third-order numerical method for size-structured population equations[END_REF], for its anti-dissipative properties.

We can see in Figure 7 that the numerical simulations are very similar to the experimental curves. In particular, the delay of twenty hours before the effect of the treatment on the growth of the population is apparent, and this twenty hour period pulses two more times as the population approaches a new equilibrium distribution during the total time of the experiment. We have thus developed an agestructured model that can explain the dynamic effects of erlotinib on PC-9 cells, which are intrinsically dependent on the age of proliferating cells. The model is also accurate in a quantitative point of view. The quiescent fractions f which are used in the numerical simulations are not chosen arbitrarily to obtain the adequate behavior. They are linked to the treatment dose of erlotinib and estimated from experimental data in a rigorous way so that they are realistic.

The model has been chosen as simple as possible to be able to recover all the parameters from experimental IMT distributions. This simplicity leads to qualitative differences between observed data and simulations in Figure 7. Experimental observation suggests the population reaches equilibrium when treated with drug whereas simulated population size is still increasing at the end of the experiment. An age-structure for quiescent cells together with an age-dependent rate of death should be consider to explain this plateau effect. Experimental data also suggest that drug treatment increases population growth rate above that of untreated cells briefly, but model does not capture this. This could be obtained by considering a lower death rate for quiescent cells than for proliferating cells. But for such more complex models, additional experimental data and a new parameter estimation method would be necessary to estimate the death rates. t (h) t (h) Figure 7: Left: experimental data without erlotinib, with a 50nM dose and with a 5000nM dose. The total quantity N (t) is plotted on a log-scale. Right: numerical simulation of model ( 1) for f = 0, f = 0.6 and f = 0.84. The curves represent the evolution of ln(P (t) + Q(t)) where P (t) := ∞ 0 p(t, a) da is the total quantity of proliferating cells at time t. They are obtained by solving Equation (1) numerically with the parameters of Figure 6 and with ν = 0.004.

Conclusion

Linking experimental observation of cell behavior between the single-cell and population scales has recently been described using newly developed mathematical models [START_REF] Tyson | A method to quantify population-level cell proliferation dynamics from single-cell data[END_REF]. However, this approach does not take into account the possible cell age-dependent effects of a perturbation on cell behavior, such as would be expected if the effects occur in cells at a specific position in the cell cycle. Since these studies were performed with asynchronously dividing cell populations, it is evident that the mathematical models of these experiments should have age structure as a primary feature. In fact, to fit the data, the authors had to use an artificial time offset to account for the age-structured effects. Here we provide a formalized approach to accurately account for cell age-dependent effects on cellular behavior. A major difficulty in the parameterization of age-structured models is the determination of the age dependent division rate. Our study provides a method for the quantitative recovery of this rate by fitting experimental IMT distributions to special forms, such as gamma functions or exponential modified Gaussians. This model, once successfully parameterized, is very useful for simulating and analyzing age-dependent phenomena in cell population dynamics. We have presented one such application for the in vitro treatment of cancer cells by erlotinib. This example shows the utility of age-structured population models in explaining the connection of drug therapy to phenomena such as cell cycle phase entry into quiescence. The method we have presented can be implemented readily to many issues in cell population behavior when there is experimental data based on cell age.

Aknowledgments

The research stays of P. Gabriel at Vanderbilt University have been financially supported by a grant of the Fondation Pierre Ledoux Jeunesse Internationale.

A Convergence to the equilibrium

The long-time behavior p(t, a) ∼ const p(a)e λt can be proven by using semi-groups methods [START_REF] Metz | The dynamics of physiologically structured populations[END_REF][START_REF] Webb | Theory of nonlinear age-dependent population dynamics[END_REF] or General Relative Entropy techniques [START_REF] Michel | General relative entropy inequality: an illustration on growth models[END_REF][START_REF] Perthame | Transport equations in biology[END_REF]. We detail here a result provided by General Relative Entropy.

Consider φ the unique solution to the adjoint eigenvalue problem For the sake of simplicity, the proof of Theorem 3 is given in the case µ = 0. But it can be easily adapted to the case with a death rate.

Proof of Theorem 3 in the case µ = 0. First we extend the initial distribution and the division rate β(a) to negative ages by setting p0 (a) = β(a) = 0 for a < 0. Since the daughters of the labeled cells are not tracked, the age distribution p(t, a) of labeled cells satisfies a transport equation without boundary condition. Thus it writes, using the characteristic method, ∀ t, a ≥ 0, p(t, a) = p0 (a -t)e -t 0 β(a-s) ds .

Introducing this expression in the definition of I T (a) we obtain, with changes of variables, if 0 ≤ a ≤ t 0 , P (t 0 ) if t 0 ≤ a ≤ T, P (t 0 ) -P (a -T ) if T ≤ a ≤ T + t 0 , 0 if a ≥ T + t 0 .

I
Because β(a) = 0 for a ≤ t 0 , we obtain by integration of I T (a) on R +

Figure 1 :

 1 Figure1: Plot of an histogram (H i ) 1≤i≤Na representing the IMT distribution of a population of PC-9 lung cancer cells. The data is from[START_REF] Tyson | A method to quantify population-level cell proliferation dynamics from single-cell data[END_REF] and the bin width is ∆a =10 8 .

Figure 3 :

 3 Figure 3: The two different gamma functions and their corresponding division rate β(a) are plotted for coefficients m = 17 and σ = 2.

Figure 4 :

 4 Figure 4: Left: The division rate β(a) is obtained numerically from (14) by using an EMG with coefficients m = 22, σ = 2 and β 0 = 0.2 to fit I ∞[START_REF] Doumic | On the calibration of a size-structured population model from experimental data[END_REF]. Right: The formula for β, as an error function with these same three parameters, is numerically indistinguishable from the numerically obtained β (R 2 = 1).

Figure 5 :

 5 Figure 5: Fitting of experimental data in [53] with the model (21). The fitting parameters are β 0 = 0.14204, m = 24.456 and σ = 3.3451. The correlation coefficient is R 2 = 0.95363 and the integral of Ĩ∞ (a|β 0 , m, σ) is ∞ 0 Ĩ∞ (a) da ≈ 1.0983. The formula for the age dependent division rate is β(a) = β 0 Erf c( m-a σ ).

Figure 6 :

 6 Figure 6: Fitting of experimental data with the model (26). The fitting parameters are β 0 = 0.17879, m = 25.007, σ = 3.6141 and µ = 0.00333. The correlation coefficient is R 2 = 0.95611 and the integral of Ĩ∞ (a|β 0 , m, σ) is ∞ 0 Ĩ∞ (a) da ≈ 1.0132. The formula for the age dependent division rate is β(a) = β 0 Erf c( m-a σ ).

  t) + Q(t))

Theorem 3 .

 3 λφ(a) -∂ a φ(a) + β(a)φ(a) + µ φ(a) = 2φ(0)β(a), φ(•) ≥ 0, p(a)φ(a) da = 1. Then the General Relative Entropy method allows us to prove that ∞ 0 |p(t, a)e -λt -ρ 0 p(a)|φ(a) da ---→ )p 0 (a) da. B Convergence of the IMT distribution Suppose that Assumptions (2) and (10) are satisfied and consider an initial distribution p0 such that p0 (a) = 0 for all a > t 0 . Then we have the convergence ∞ 0 |I T (a) -I ∞ (a)| da ----→ T →∞ 0 where I T and I ∞ are defined in (9) and (11).

= C - 1 T T 0 β

 10 )p 0 (a -t)e -t 0 β(a-s) ds dt (a)p 0 (a -t)e -a a-t β(a ′ ) da ′ dt= C -1 T a a-T β(a)p 0 (u)e -a u β(a ′ ) da ′ du = C -1 T β(a)e -a 0 β(a ′ ) da ′ a a-T p0 (u)e u 0 β(a ′ ) da ′ du.Since the support of p0 is included in [0, t 0 ] and β(a) = 0 for a ∈ [0, t 0 ] due to Assumption (2), we have for all u ∈ R p0 (u)e u 0 β(a ′ ) da ′ = p0 (u) so I T (a) writesI T (a) = C -1 T β(a)e -a 0 β(a ′ ) da ′This function is nondecreasing and satifies P (a) = 0 for a ≤ 0 and P (a) = P (t 0 ) = ∞ 0 p0 for a ≥ t 0 , so we have a a-T p0 (u)du = P (a) -P (a -T ) =

1 = C - 1 T P (t 0 ) T 0 β 1 T T +t 0 T 1 T 0 β+ T +t 0 T C - 1 T

 10101001 (a)e -a 0 β(a ′ ) da ′ da + C -β(a)e -a 0 β(a ′ ) da ′ (P (t 0 ) -P (a -T )) da which givesC T -P (t 0 ) T 0 β(a)e -a 0 β(a ′ ) da ′ da = T +t 0 T β(a)e -a 0 β(a ′ ) da ′ (P (t 0 ) -P (a -T )) da ≤ P (t 0 ) T +t 0 T β(a)e -a 0 β(a ′ ) da ′ da -----→ )e -a 0 β(a ′ ) da ′ da = 1, a) -I ∞ (a)| da ≤ C -1 T P (t 0 ) -(a)e -a 0 β(a ′ ) da ′ da [P (t 0 ) -P (a -T )] -1 β(a)e -a 0 β(a ′ ) da ′ da + ∞ T +t 0 β(a)e -a 0 β(a ′ ) da ′ da ≤ C -1 T P (t 0 ) -1 T 0 β(a)e -a 0 β(a ′ ) da ′ da + K ∞ T β(a)e -a 0 β(a ′ ) da ′ da ----→T →∞ 0 and it ends the proof of Theorem 3.

  It is a transport equation which satisfies a maximum principle, namely if the initial distribution p 0 (a) is nonnegative (positive) then the distribution p(t, a) remains nonnegative (positive) for all time t > 0. The solutions of (6) have a remarkable behavior as time evolves, in that the age structure equilibrates no matter what the initial age structure of the population may be. This effect is known as asynchronous exponential growth, and its interpretation is that the population disperses over time to a limiting asymptotic equilibrium age structure where the fraction of the population in any age range [a 1 , a 2 ] satisfies independent of the initial age structure[START_REF] Webb | An operator-theoretic formulation of asynchronous exponential growth[END_REF]. Moreover, the solutions to this equation, as t → ∞, are known to behave like a separated variables solution, that is, p(t, a) = a function of time only × a function of age only. More precisely, consider the eigenvalue problem

	lim t→∞	a 2 a 1 p(t, a) da ∞ 0 p(t, a) da
	= a constant