Age-structure as key to delayed logistic proliferation of scratch assays

Scratch assays are in-vitro methods for studying cell migration. In these experiments, a scratch is made on a cell monolayer and recolonisation of the scratched region is imaged to quantify cell migration rates. Typically, scratch assays are modelled by reaction diffusion equations depicting cell migration by Fickian diffusion and modelling proliferation by a logistic term. In a recent paper (Jin, W. et al. Bull Math Biol (2017)), the authors observed experimentally that during the early stage of the recolonisation process, there is a disturbance phase where proliferation is not logistic, and this is followed by a growth phase where proliferation appears to be logistic. The authors did not identify the precise mechanism that causes the disturbance phase but showed that ignoring it can lead to incorrect parameter estimates. The aim of this work is to show that a non-linear age-structured population model can account for the two phases of proliferation in scratch assays. The model consists of an age-structured cell cycle model of a cell population, coupled with an ordinary differential equation describing the resource concentration dynamics in the substrate. The model assumes a resource-dependent cell cycle threshold age, above which cells are able to proliferate. By studying the dynamics of the full system in terms of the subpopulations of cells that can proliferate and the ones that can not, we are able to find conditions under which the model captures the two-phase behaviour. Through numerical simulations we are able to show that the resource concentration in the substrate regulates the biphasic dynamics.


Introduction
Cell migration and proliferation are central processes in the development and homoeostasis of many organisms. ey also play a key role in pathological processes such as cancer invasion, chronic in ammatory diseases and vascular diseases [1,2,3]. To understand the biochemical and physical cues that regulate cell migration and proliferation, in-vitro assays are used to quantify the ability of cell populations to migrate and proliferate under controlled situations [4,5,6]. Scratch assays are typically used to study cell migration [6,7]. A scratch assay involves: growing a cell monolayer to con uence; creating a "scratch" in the monolayer; and monitoring the cell dynamics as the scratch closes [7]. e resulting time-course data are then analysed to estimate the cell migration rate [8,9]. Other common cell-based assays are proliferation assays that focus on measuring the cell number or the proportion of cells that are dividing [10]. e simplest proliferation assay consists of growing a cell monolayer to low density on a two-dimensional substrate and measuring the cell number change in time [11,12,10]. See Fig. 1 (a) and (b) for a schematic representation of the typical time progression of these two in-vitro assays.
Many mathematical models have been developed to describe these in-vitro assays and to test hypotheses on the mechanisms that govern cell migration and proliferation [13,14,15,16,6,17,18,19,20]. e macroscale dynamics of scratch assays are o en described by reaction-di usion equations [15,17,18] in which cell migration is modelled as Fickian di usion and cell proliferation is modelled by a logistic source term. For the proliferation assay, where no scratch is performed, the logistic growth equation is widely used [14,16,6]. Model selection is commonly performed by matching the experimental time-course data with the model evolution and accepting the model if the residual error is small [21,22,23]. However, several authors have pointed out that this is not enough for model validation [24,25].
In [21] the ability of the logistic growth model to describe cell proliferation in scratch assays was studied. A series of scratch and proliferation assays using PC-3 prostate cancer cells were performed and the changes in cell density in two subregions located far from the scratch were quanti ed. ese two subregions were chosen so the changes in cell density were not due to cell migration. To assess the suitability of the logistic growth model, they analysed the model t to experimental data and the per capita growth rate of the experimental data, σ(N ) = dN dt 1 N , as a function of the cell density, N . Calibrating solutions of the logistic model to the experimental data showed a good t for both assays, however, analysing the per capita growth rate revealed di erent behaviours between the proliferation and the scratch assays. e authors observed that for proliferation assays, the per capita growth rate could be well described by a linearly decreasing function of the cell density (see Fig. 1 (e)), a result consistent with the logistic model. However, for the scratch assay data, during the rst 18 hours of the experiment, the per capita growth rate was found to increase with cell density, whereas for 8 ≤ t ≤ 48h, the per capita growth rate was found to decrease approximately linearly with the cell density (see Fig. 1 (f)).
Guided by their experimental observations, the authors in [21] proposed that cell proliferation in scratch assays involves two phases: an initial disturbance phase during which proliferation is not logistic and a growth phase during which proliferation is approximately logistic. Since this behaviour was not observed in the proliferation assays, the authors concluded that it was caused by the scratching procedure. ey hypothesised that scratching might create chemical or mechanical disturbances but did not test their hypothesis.
In this work, we show that a nonlinear age-structured model can account for the disturbance and growth phases observed in scratch assays. We refer to a cell's age as the cell's temporal position within the cell cycle. e cell cycle consists of four phases: rst gap G 1 , synthesis stage S, second gap G 2 and mitosis stage M . G 1 and G 2 consist of growth phases, DNA replication occurs in S phase and cell division in M phase [26]. e progression in the cell cycle is regulated by environmental cues and intracellular checkpoints involving various proteins, in particular cyclins and cyclin dependent kinases [26,27].
By considering an age-structured model we can investigate how heterogeneity in the cell population age distribution a ects the total cell population dynamics. Previous work has shown that heterogeneity in cell age distribution may still generate logistic growth at the population scale [28]. Other age-structured models have described how the cell age distribution can a ect the overall cell population dynamics in scratch assays: the speed with which the cells invade the vacant region [29,30], the e cacy of anti-cancer drugs, particularly phase-speci c drugs [31] and the in uence of growth factors [13]. Agent-based models have also been considered to study the e ect of heterogeneity in cell age distribution on the overall dynamics in in-vitro assays [19,20].
We consider an age-structured model that was rst introduced in [32]. e model considers the interplay between a single cell population and the resource concentration available in the substrate. e model assumes a resource-dependent G1/S transition age, above which, cells are able to proliferate. is critical age naturally divides the cell population into "immature" and "mature" cells. By studying the dynamics of the full system in terms of these subpopulations, we are able to nd conditions under which the cell population evolution is of logistic-type and the per capita growth rate follows a biphasic behaviour. We In the plots (c)-(f), we consider the experimental data from a proliferation and a scratch assay performed in [21]. In (c) and (d) the mean cell population is plo ed over the duration of the experiment. In (e) and (f) the per capita growth rate σ(N ) = dN dt 1 N is plo ed with respect to the cell population N . We calculate the per capita growth rate in the same way as in [21]. A biphasic trend can be observed in (f) (for the scratch assay) but not in (e) (for the proliferation assay). validate our predictions via numerical simulations and show that the resource concentration regulates the disturbance phase. e remainder of this paper is organised as follows. In Sect. 2, we describe the nonlinear age-structured model. In Sect. 3, we describe the dynamics of the full model in terms of the mature and immature subpopulation dynamics. We then derive conditions under which the per capita growth rate will exhibit biphasic behaviour and logistic-type proliferation. In Sect. 4, we investigate numerically these conditions and determine parameter regimes in which the resource concentration regulates the disturbance phase. We conclude in Sect. 5 with a discussion of the results and conclusions.
2 An age-structured model with resource-regulated proliferation To take into account the cell cycle heterogeneity of a cell population, we consider a McKendrick-von Foerster model [33,34]. e model considered in this work was rst used to describe the mean eld dynamics of a stochastic multiscale model of a cell population with oxygen-regulated proliferation [32]. Here the model describes the dynamics of an in-vitro cell population assumed homogeneously distributed in space, like cell cultures in proliferation assays or cells in far away regions from the wound in scratch assays, as in [21]. We are interested in the population dynamics with respect to the resources available in the medium. We denote by n(a, t), the number of cells of age a at time t. We consider a ∈ [0, ∞]. We denote by T > 0 the duration of the experimental observation. e model focuses on the cell population proliferation dependence on the growth factors, oxygen and nutrients available in the medium. We refer to these components as a single, generic resource and denote it by c(t).
We assume cells mature with constant speed 1, die with rate µ and proliferate at a rate b(a, c(t)) which we consider to be age and resource-dependent rate [35,36]. Combining the above assumptions, it is straightforward to show that the evolution of the cell density function n(a, t) We suppose further that when a cell divides it produces two daughter cells of age a = 0. By considering all possible division events, we deduce that We consider a coarse-grained description of the cell cycle by lumping S, G2 and M into one phase, so we consider a two phase model G1 and S-G2-M. We assume cells proliferate at a constant rate, τ −1 p , provided they successfully enter the S-G2-M phase. Entry to the S-G2-M phase is regulated by the G1/S checkpoint which has been shown to depend on the presence of resources needed for cell growth [37]. erefore, we assume that the proliferation rate is given by where H is the Heaviside function and a G1/S (c(t)) denotes the G1/S transition age, the age at which a cell passes from the G1 to the S phase. e G1/S transition age, a G1/S (c) is given as a function of c to capture the dependence of the checkpoint on the available resource concentration. e authors in [32] describe this dependence by a simple scaling which was derived from analysing an oxygen-dependent cell cycle progression model, where a − and β are positive constants. e positive quantity c cr is the critical resource concentration that allows cell proliferation. For high resource values, the transition age becomes smaller, so the cell requires less time to proliferate. On the other hand, for resource values slightly higher than c cr , the transition age becomes bigger and cells take longer to transition to the G1/S phase and be able to proliferate. For c ≤ c cr , cells are assumed to enter the quiescent state. If we assume that the resource, c(t), is supplied at a constant rate S and consumed by all cells at a constant rate k; then where N (t) = ∞ 0 n(a, t)da denotes the total number of cells at time t. e full model consists of Eqs. (1) -(4) and the initial conditions e long-time dynamics of this model was described in [32]. Before beginning our analysis, we summarize below those results from [32] that are relevant for our work.
• For large times, the solution of the age-structured model approaches a separable solution, Q(a, t) = A(a) exp(λt) i.e. n(a, t) ≈ Q(a, t) for t >> 1, where λ satis es the Euler-Lotka equation: in which c ∞ > 0 is the steady state value of the resource concentration. A(a) is referred as the stable age-distribution of the model. e separable solution is assumed to be in equilibrium when λ = 0 and therefore N (t) → N ∞ and dN dt → 0 as t → ∞.
• For the separable solution to be in equilibrium, the average waiting time to division a er the G1/S transition, τ p , must be smaller than the average cell life span, µ −1 .
• If Eq. (7) is satis ed, a er solving Eq. (6) with λ = 0, they obtained that: e transition age, a G1/S , reaches the steady-state value e steady-state resource concentration is given by e steady-state value of the total cell population density is given by: e above results will be useful for parametrising our model.

Model analysis
In this section, we derive necessary and su cient conditions for Eqs. (1) -(5) to exhibit logistic-type behaviour and biphasic behaviour in the per capita growth rate. ese conditions can be obtained by distinguishing two subpopulations: cells whose age is greater than the G1/S transition age and which can proliferate, and the cells whose age is below the threshold. We refer to these subpopulations as, respectively, mature and immature cells. We now analyse the full system with respect to these subpopulation dynamics.

Mature and immature subpopulation dynamics
Let us denote by X(t) and Y (t), respectively, the subpopulations of mature and immature cells so that where a G1/S (t) := a G1/S (c(t)).
With these de nitions, it is possible to obtain the following result: eorem 1. Let X and Y be de ned as in (11) and let n(a, t) and c(t) satisfy the model Eqs.
(1) -(5). en the dynamics of X and Y are given by the following system: and c(t) is formally given by Furthermore, if c(t) ≥ c cr ∀t ∈ [0, T ], then the number of cells with G1/S transition age a G1/S (t) at time t, n(a G1/S (t), t), is formally given by for t ≤ a G1/S (t) and for t > a G1/S (t). Proof.
Integrating Eq. (1) over the domains [0, a G1/S (t)] and [a G1/S (t), ∞], we obtain the dynamics of X and Y , respectively. e evolution of c(t) is given by solving the di erential equation (4) while assuming N (t) is a known function. Finally, the formula for the cell concentration with G1/S transition age a G1/S (t) at time t, n(a G1/S (t), t), is obtained by solving the age-structured model by the method of characteristics. e characteristic curves of Eq. (1) are the lines a = t + µ where µ is a constant. By solving Eq. (1) along the characteristic curves and taking into account the boundary condition (2) and initial condition (5), we obtain that n(a, t) is given by: By considering a = a G1/S (t) and n(0, t − a G1/S (t)) in terms of the mature subpopulation Y in Eq. (15), we derive the expressions (13) and (14) for n(a G1/S (t), t).
eorem 1 reduces the analysis of the full model to the analysis of the system (12). Given Eqs. (13) and (14), we notice that for t ≤ a G1/S (t), (12) consists of a non-homogeneous coupled linear system of di erential equations that depends on the initial age-distribution and for t > a G1/S (t), the dynamics of Y is described by a state-dependent delay di erential equation, in which a(t) Having the dynamics of the full model expressed in terms of these subpopulations, we can now describe the overall proliferation and the per capita growth rate in terms of the dynamics of these two subpopulations dynamics. Let us denote by σ(N (t)) := dN (t) dt 1 N (t) , the per capita growth rate. e following theorem gives the evolution of the per capita growth rate and the total population evolution in terms of the mature and immature cell subpopulations. eorem 2. e evolution of N (t) and σ(N (t)) for t ∈ [0, T ] is given by: from which we obtain Eq. (18).

Conditions for delayed logistic proliferation
In this section we identify conditions under which the per capita growth rate of the age-structured model exhibits biphasic dynamics and logistic-type growth. First, we interpret the experimental behaviour in mathematical terms: B1 e cell population undergoes a logistic-type behaviour, i.e, • e population increases in size monotonically, B2 e per capita growth rate exhibits biphasic behaviour, i.e. there exists t 1 > 0 such that the per capita growth rate, σ(N (t)), has the following behaviour: eorem 3. Necessary and su cient conditions for the age-structured model with resource-regulated proliferation, Eqs. (1) -(5), to exhibit the behaviour B1 and B2 are as follows and there exists t 1 > 0 such that Proof. Given the relationship in Eq. (17), it is clear that (20) is true if and only if the derivative of N is positive. To derive (21), we notice from the chain rule that dσ(N (t)) dt = dσ(N ) dN dN dt and given that dN dt > 0, the sign of d dt (σ(N (t))) is the same as dσ(N ) dN . Finally, the condition (19) needs to be ful lled so that the age-structured model has a stable steady-state solution as described in [32]. eorem 3 gives necessary and su cient conditions for the model to exhibit the behaviour described in B1 and B2. We note that the conditions for the biphasic behaviour are given by Eq. (21) which can be interpreted biologically as an initial increase in the fraction of mature cells followed by a stabilization phase in which the mature cell fraction decreases to its steady-state value. Since the fraction of mature cells is regulated by the transition age, a G1/S , and this critical age is regulated by the resource concentration, we deduce that changes in the resource dynamics may in uence the biphasic behaviour.

Numerical results
In this section, we perform numerical simulations to verify the predictions of eorem 3. First, we present the discretization scheme that we use to solve Eqs. (1) -(5).

Discretisation scheme
e model is composed of a hyperbolic partial di erential equation with a nonlocal boundary condition coupled to an ordinary di erential equation. We use a spli ing method to discretise the coupled system. e ordinary di erential equation is discretised by an implicit Euler method. e partial di erential equation is discretised using the Rothe method: rst discretised in time via the implicit Euler method and then in space by linear nite elements [38]. e method is stabilised by considering the streamline upwind/Petrov Galerkin formulation [39]. e nite element discretisation is implemented in C++ using the so ware DEAL.II [40].
To discretise the age domain [0, ∞], we notice that the cell age is naturally bounded. Hence there exists a max > 0 such that n(a, t) = 0 for a ∈ (a max , ∞), and so we restrict our a ention to a ∈ [0, a max ]. For the nite element discretisation to be computationally e cient and not create spurious oscillations around the discontinuity point of the proliferation rate (given by the Heaviside function), we discretise the age domain [0, a max ] di erently in two subdomains: the subdomain [0, a * ], where a * is given by Eq. (8), is discretised with a step size of k = a * /2 10 ≈ 0.00469081, and the domain [a * , a max ] is discretised with a step size of k = (100 − a * )/2 10 ≈ 0.0442. We consider a smaller step size for the rst domain since for the simulations considered in this work, the evolution of a G1/S (t) occurs in this subdomain. For the time discretisation, we consider a time step of h = 0.1.

Model parametrisation
We consider as a case study, PC-3 prostate cancer cells. In order to parametrise our model, we use as a guideline the parameters that were estimated in [28] from the experimental time-course data of scratch assays. e authors estimated a proliferation rate ofλ = 0.053 ± 0.005 [hr] −1 . We assume their estimate corresponds in our model to τ −1 p , i.e.
Estimates of doubling times of PC-3 prostate cancer cells are in the range of 25-33 hours [41,42]. In our model, the doubling time, which we denote as DT, is given by When the cell population reaches the steady-state age distribution, DT is given by: By considering this formula, the estimated value for the proliferation rate τ −1 p (Eq. (22)) and the range of values for the doubling time (Eq. (23)), we obtain a range of values for the death term, µ ∈ [0.0233−0.0333]. We assume that µ = 0.0283 [hr] −1 (the midpoint value). e carrying capacity was also estimated in [28]  [hr] [cell] . Using this value for the consumption rate, k, and the estimate of the carrying capacity (Eq. (26)) in the formula of the steady-state value of the total population in our model (Eq. (10)), we determine the value of the resource ux,S = 26.63 × 10 −4 [nM ] [hr] . Regarding the parameters that describe the resource dependence of the G1/S transition age, we consider the parameter values that were derived from an intracellular model for an oxygen-regulated proliferation rate in [32] the dimensionless parameter values are β = 0.2 and c cr = 0.23. For determining the parameter value of a − , we notice that the admissible range of values for a G1/S (t) is 6 ≤ a G1/S (t) ≤ 14 given the range values of DT (Eq. (23)) and the estimate of τ p (Eq. (22)). We consider a − = 8.25 [hr], so a G1/S (t) is in this range for all our simulations.
We consider the following initial age-distribution that is a multiple of the steady-state age distribution: where a * is the steady-state value of the transition age (Eq. (8)) and κ is estimated so the initial number of cells corresponds to the average number of cells consider as initial condition in [28] for the most con uent initial condition, N (0) = 223. We set c 0 = c ∞ as the steady-state value of the resource concentration given by Eq. (9). e default parameters are listed in Table 1.

Reference dynamics
Under the default parameter values, the model satis es the conditions of eorem 3 and therefore we expect a logistic-type behaviour and a biphasic dynamics. In Fig. 2 (a) and (b) we plot the time evolution of the total cell population N (t) and of the resource concentration c(t), respectively. As expected, the total cell population follows a logistic-type behaviour: it grows exponentially and saturates. e resource concentration has an initial increase until it reaches a maximum value and then decreases until it reaches its steady-state value. In Fig. 2 (c) we plot the evolution of the transition age a G1/S and observe its dependence on the resource concentration: it decreases to a minimum value and then increases towards the steady-state value given by Eq. (8). In Fig. 2 (d) we plot the per capita growth rate σ(N ) as a function of the total cell population N ; as expected, we observe biphasic behaviour, with an initial increase in σ with respect to cell density followed by a linear decrease.

Sensitivity analysis
In this subsection we perform a sensitivity analysis focusing on those parameters that modulate the resource dynamics in order to understand their role in determining the biphasic behaviour of the per capita growth rate. We focus our sensitivity analysis on two parameters that can be experimentally manipulated: the resource ux, S, and the initial concentration of resource, c 0 . We rst vary the resource ux, S, in the domain [25.63 × 10 −4 − 28.63 × 10 −4 ]. We consider this domain so the cell population steady-state value, N ∞ , is within the con dence interval given by Eq. (25). As expected, from Eq. (10), the steady-state value at which the total cell population saturates, N ∞ , increases as the value of S increases (Fig. 3). We observe also that the maximum value of c(t) increases and the minimim values of a G1/S (t) decreases as S increases, see  1) and (4). In (a) we plot the total cell population evolution and observe it follows a logistic-type growth. In (b) we plot the resource evolution that follows a rapid increase and then a monotonic decrease to the steady-state value. In (c) we plot the evolution of the transition age a G1/S (t) and the inverse dependence of the resource concentration on the transition age can be observed. In (d) we plot the per capita growth rate against the total cell population evolution for which two proliferation phases can be observed (a rapid increase, then slower decrease). Parameter values as per Table 1.
We now vary the initial resource concentration c 0 . We focus only on increasing the value of the resource concentration above c ∞ , since considering values below c ∞ will make the population decrease and this will be in disagreement with experimental observations. In Fig. 4 we show that when we increase c 0 , from 0.034 to 0.114, the duration of the disturbance phase is shortened. e model predicts that by increasing c 0 , the total cell population evolution does not change signi cantly, however the resource evolution changes from biphasic to monotonic, i.e. the duration of the disturbance phase decreases and eventually disappears as c 0 increases (see Fig. 4 (b)). Increasing the value of the initial resource concentration shi s the observation time of the dynamics and now it is only possible to observe the resource concentration decline to the steadystate value. is e ect can be similarly observed in Fig. 4 (c) where the initial decrease of the transition age a G1/S evolution as c 0 increases. We observe that increasing the initial resource concentration decreases the duration of the disturbance phase where the dependence of the per capita growth rate with respect to the total cell population is not linear (see 4 (d)). We notice that during the disturbance phase, the per   1) and (4), as we vary the initial resource concentration value, c 0 . Increasing c 0 : (a) does not a ect the total cell population, increases the maximum resource concentration but reduces the time it takes to reach it, and (c) makes the initial decrease of the transition age a G1/S disappear. (d) Increasing the initial resource concentration, c 0 , a ects the plots of the per capita growth rate against the total cell population by shortening the duration of the disturbance phase where there is no monotonic decreasing dependence. Parameter values as per Table 1. capita growth rate does not increase monotonically with respect to the total cell population as we saw previously for the reference dynamics and when varying the resource ux, S. For the simulations with c 0 bigger than c ∞ = 0.034, the per capita growth rate initially decreases, then increases and nally follows a linear decreasing trend with respect to the total cell population.
Finally, we perform a sensitivity analysis on the rest of the model parameters and consider di erent initial age-distributions to show that our model predictions are robust. We observe that uncertainty of the parameters, on the range of values that we considered, has li le e ect on the overall dynamics. e two proliferation phases are present for all the set of parameters considered (see Supporting Information section 1). We then consider di erent inital age-distributions that ful ll eorem 3 and show that the biphasic dynamics is present in all cases considered (see Supporting Information section 2).

Discussion and conclusions
In this work we have presented an age-structured model with resource-dependent proliferation rate that captures for the rst time the biphasic behavior in the per capita growth observed experimentally in [21]. We analysed the full model in terms of two subpopulations: cells which are able to proliferate or not (mature and immature cells). We then derived necessary and su cient conditions under which the model presents a logistic type behavior and two phases of proliferation: an initial phase, which in [21] was named disturbance phase, in which proliferation does not follow a logistic growth and a growth phase, where proliferation is approximately logistic. e biphasic behaviour was demonstrated to be a result of an initial increase of the fraction of mature cells followed by a decline to their steady-state concentration. We then parametrised the model using PC-3 prostate cancer cells as a case study. Finally, through numerical simulations, we showed that varying the resource initial value the plot changes the dependence of the per capita growth rate against the total cell population. e model predicts that the duration of the disturbance phase is decreased as the initial resource concentration is increased.
In view of this model, the experimental observations in [21] can be explained as follows: the scratch procedure decreases the cell number in the plate and by replenishing the medium in the same quantity, as customary, the resource concentration is increased, therefore triggering the biphasic dynamics in the per capita growth rate. e predictions of the model can be experimentally tested by modifying the resource concentration in the substrate and examining the plots of per capita growth rate against the total cell population as was performed in [21]. e model predicts that by increasing the initial resource concentration, the disturbance phase would shorten. Varying the resource concentration has been shown to a ect the overall dynamics of the scratch assay in other situations [13].
ere are di erent ways that we can extend this work. e model was parametrised using estimates calculated in [21]. eir estimates were calculated with respect to the logistic equation and since we are considering a di erent model, the values are expected to di er. e model could be parametrised using more appropiate methods as has been done for similar models [44,13]. Other alternative hypothesis could be analysed: such as mechanical and chemical disturbances [45]. An inference-based modelling approach could be performed to test which model best explains the experimental data or design new experiments to discriminate between feasible alternatives.
In summary, this work identi es the interplay between the mature subpopulation and the resource concentration as the responsibles for the biphasic behaviour observed experimentaly in [21]. Experimental validation and further mathematical modelling will help elucidate the impact of heterogeneity in cell age distribution in the overall dynamics.