Quantifying cell cycle regulation by tissue crowding

The spatiotemporal coordination and regulation of cell proliferation is fundamental in many aspects of development and tissue maintenance. Cells have the ability to adapt their division rates in response to mechanical constraints, yet we do not fully understand how cell proliferation regulation impacts cell migration phenomena. Here, we present a minimal continuum model of cell migration with cell cycle dynamics, which includes density-dependent effects and hence can account for cell proliferation regulation. By combining minimal mathematical modelling, Bayesian inference, and recent experimental data, we quantify the impact of tissue crowding across different cell cycle stages in epithelial tissue expansion experiments. Our model suggests that cells sense local density and adapt cell cycle progression in response, during G1 and the combined S/G2/M phases, providing an explicit relationship between each cell cycle stage duration and local tissue density, which is consistent with several experimental observations. Finally, we compare our mathematical model predictions to different experiments studying cell cycle regulation and present a quantitative analysis on the impact of density-dependent regulation on cell migration patterns. Our work presents a systematic approach for investigating and analysing cell cycle data, providing mechanistic insights into how individual cells regulate proliferation, based on population-based experimental measurements.


Introduction
The coordination of cell proliferation across space and time is crucial for the emergence of collective cell migration, which plays a fundamental role in development, including tissue formation and morphogenesis, and also at later stages for tissue regeneration and homeostasis.Cells adapt their division rates in response to mechanical constraints within tissues [1,2], allowing cell populations to self-organise and eventually form and maintain tissues and complex structures.Moreover, disruptions in the control of cell proliferation often result in tumour formation [3][4][5].Although significant experimental efforts have been devoted to understand the mechanical regulation of cell proliferation [6] and its interplay with collective cell migration, existing mathematical models have failed to describe these constraints and how they affect cell cycle progression [7][8][9].
In order to understand cell proliferation regulation, numerous experimental studies have explored how spatial and mechanical constraints within tissues affect different stages of the cell cycle.The cell cycle consists of four main stages, namely: the G1 phase, where cells grow and prepare for DNA replication; the S phase, during which DNA synthesis occurs; the G2 phase, characterised by further cell growth and preparation for mitosis; and finally, the M phase, where cell division takes place.Cells can also exit the cell cycle and enter G0, where they become quiescent.The experimental visualisation of cell cycle stages can be achieved via the widely used FUCCI cell-cycle marker [10], which consists of red and green fluorescent proteins that are fused to proteins Cdt1 and Geminin, respectively.Cdt1 exhibits elevated levels during the G0/G1 phase and decreased levels throughout the remaining cell cycle stages, whereas Geminin shows high expression during the S, G2, and M phases; allowing thus to distinguish between these different stages -see Fig. 1.Several extensions of the FUCCI system exist now [11]; for instance FUCCI4 allows for the simultaneous visualisation of the G1, S, G2, and M phases [12].
Experimental studies of cell migration are often performed in epithelia due to their strong cell-cell adhesion which gives rise to collective and cohesive motion.Moreover, they play a fundamental role in multicellular organisms as they serve as protective layers for various body surfaces and organs.Epithelial cell proliferation is regulated by mechanical forces, which can accelerate, delay, arrest, or re-activate the cell cycle.In particular, extensive research has focused on the G1-S boundary, revealing that intercellular tension can favour this transition [13], while tissue pressure can halt progression based on crowding [2].
More generally, the extracellular regulation of switches from G0 and G1, and within substages of G1, has been well-known for many years [14].However, and contrary to initial assumptions, cells also have the ability to regulate progression through stages of the cell cycle following the G1-S transition in response to external cues.These external signals might involve not only mechanical forces [15], but also nutrients and growth factors [16,17].In epithelia, this question was explored recently by Donker et al. [18], revealing a mechanical checkpoint in G2 which controls cell division.In particular, this checkpoint allows cells to regulate progression through G2, via sensing of local density, explaining why dense regions in epithelia contain groups of cells that are temporarily halted in G2.
Experimental studies employing FUCCI and variations of it have thus successfully linked mechanical constraints to cell cycle progression.These studies have employed qualitative analysis, direct measurements of cell cycle stage durations     < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 1 o z F 3 / i h f 7 C E a 0 P a 3 q g 4 0 o i s S g = " > A A A B 7 X i c b V B N S w M x E J 3 U r 1 q / q h 6 9 B I v g q e w W U Y 9 F L x 4 r 2 A 9 o l 5 J N s 2 1 s N l m S r F C W / g c v H h T x 6 v / x 5 r 8   [19] at different time points (adapted).Initial tissue diameter ∼ 3.4 mm.Scale bars correspond to 1 mm.(C) Segmented data showing G1 (red), S/G2/M (green), and post-mitotic (gray) cells.Note that in the model we combine post-mitotic cells and cells in G1. (D) Zoomed-in segmented data at the tissue edge and centre, corresponding to the black squares in (C).(E) Fraction of cell-cycle state cells in the tissue centre and in the tissue edge -defined as regions extending ∼ 200 µm from the tissue center and tissue edge, respectively.(F) Density profiles in polar coordinates at t = 40 h, showing cells in G1, S/G2/M, and post-mitotic cells.(E) and (F) show the average of eleven independent tissue expansions with the same experimental initial condition, with shaded regions indicating one standard deviation with respect to the mean.[18], or metrics associated with cell cycle progression, such as cell area [2], and Geminin/Cdt1 or EdU signals [20,21].However, these approaches omit a quantitative comparison between model and data, hence limiting the depth of mechanistic insights that can be derived.
Here, we present a quantitative investigation into the mechanical regulation of cell cycle progression by sensing of local tissue density First, we construct a mathematical model of cell cycle dynamics that accurately captures the impact of tissue crowding on cell cycle progression.By combining minimal mathematical modelling, Bayesian inference, and recent experimental data [19], we provide further evidence, consistent with previous experimental studies [2,18], that density-dependent effects operate throughout the cell cycle and together serve as a regulating mechanism for the growth of epithelial tissues.Our work thus constitutes a systematic approach towards the quantification of densitydependent effects regulating cell cycle progression.Moreover, the obtained parameter estimates reveal an explicit relation between the duration of different cell cycle stages and tissue density, which is consistent with the experimental measurements of Donker et al. [18].

Methods
Mathematical models of cell cycle dynamics.We build on the model proposed by Vittadello et al. [7] to describe two cell populations, ρ 1 (x, t) and ρ 2 (x, t), in different stages of the cell cycle.We represent by ρ 1 the density of cells that are in G0/G1, while ρ 2 gives the density of cells in the S/G2/M phases of the cell cycle -see Fig. 1.In the model, cell motility is described via linear diffusion, with a diffusion constant D > 0 for both cell populations [22].In order to effectively capture density-dependent effects controlling cell cycle progression, we assume that the transitions between different cell cycle stages are regulated by two crowding functions, f (ρ) and g(ρ), which depend on the total cell density ρ = ρ 1 + ρ 2 .In particular, the transition rate from G1 to S is given by k 1 f (ρ), while the division rate (from S/G2/M to G1) is given by k 2 g(ρ), where k 1 , k 2 > 0 are intrinsic rates of cell cycle progression.With this, the model reads where the factor of two in the equation for ρ 1 represents cell division into two daughter cells, and ∆ = d i=1 ∂ 2 xi is the Laplacian operator in dimension d.These equations are solved first in polar coordinates (assuming radial symmetry in two spatial dimensions, d = 2) to describe epithelial tissue expansion experiments, and then in one spatial dimension (d = 1) to study travelling wave behaviour, and the impact of tissue crowding on cell migration phenomena.
In order to accurately capture density-dependent effects regulating cell cycle progression, we assume that f and g are non-increasing functions of the total density ρ.Again, this is motivated by the experimental observations of Streichan et al. [2] and Donker et al. [18].Furthermore, we assume f (0) = g(0) = 1, so that k 1 and k 2 represent density-independent transition rates.Note that setting f = g ≡ 1 gives rise to an exponential growth model (i.e.no dependence on density).On the other hand, choosing f ≡ 1 and g(ρ) = (1−ρ/K) + we recover the Vittadello et al. model [7].Here, we assume that f (ρ) and g(ρ) decrease linearly with the total cell density so that where K 1 , K 2 > 0 are constants controlling the duration of G1 and the S/G2/M phases, respectively, and (z) + = max(z, 0).The specific form of these crowding functions is chosen here for simplicity, although other functions sharing the same properties show similar qualitative behaviour.We follow a Bayesian approach [8,[23][24][25] to calibrate the model given in Eqs.(1).In particular, given experimental measurements of the cell densities {ρ D k (x i , t j )} i,j for k = 1, 2, and a vector of model parameters θ = (D, k 1 , k 2 , K 1 , K 2 ), we estimate the posterior probability distribution p(θ|ρ D ), which gives the probability density for the model parameters taking specific values.The posterior distribution, thus, can be used to quantify the uncertainty associated with specific parameter values, given the experimental observation.We refer the reader to the Supplementary Information for more details on Bayesian inference.

Results
Tissue expansion experiments.We compare our model predictions to the experiments performed by Heinrich et al. [19] studying the expansion and growth dynamics of a single circular epithelial tissue -see Fig. 1B.In these experiments, MDCK cells expressing the FUCCI markers are cultured in a silicone stencil for 18 hours and, after stencil removal, the cell population is allowed to freely expand for 46 hours.Given that the average cell cycle duration for MDCK cells is around 16 hours, this enables each cell to potentially undergo 2-3 cell divisions during the experiment.Local densities are then quantified by segmenting the fluorescence images in ImageJ and counting the number of nucleus centroids -Fig.1C.Note that post-mitotic cells do not fluoresce and appear dark, which makes the FUCCI system unreliable for cell counting.To quantify the density of post-mitotic cells, Heinrich et al. used a convolutional neural network to identify nuclei from phase contrast images [26] -see [19] for more details.Moreover, and in line with previous work [24], the model takes as initial condition the quantified density profile ten hours after stencil removal, so that the impact of the stencil on the dynamics is reduced.Note that after this time, cell densities near the tissue centre are relatively high (∼ 3500 cells/mm 2 , which corresponds to around 50-70% of Squares represent the estimated cell density obtained by averaging eleven experimental realisations, which we use to calibrate the model.Shaded regions denote one standard deviation with respect to the mean -see Supplementary Fig. 2 for confidence intervals in the model predictions.Numerical simulations in polar coordinates were obtained by using the posterior modes as parameter values, and no-flux boundary conditions -for details on the numerical scheme we refer to the Supplementary Information.In order to minimise the effects of the stencil removal on cell behaviour, the initial condition corresponds to the experimental density profile ten hours after stencil removal. the maximum saturation density for MDCK [19,24]) and a fraction of cells in this region are likely to be found in a quiescent state due to contact inhibition of locomotion and proliferation [27].
The experiments by Heinrich et al. [19] reveal a higher density of cells in G0/G1 at the centre of the tissue, where the total cell density is also higher -see Fig. 1D-E.The tissue edge, in contrast, is characterised by a larger number of cells which are preparing to divide (green) or are directly post-mitotic (gray).This agrees with previous observations of epithelial cells, which are known to control progression from G1 to S in response to spatial constraints [2].Note, however, that the density of cells in S/G2/M in the tissue centre is low but non-zero, even at later times in the experimentsee Fig. 1D-F -as observed also by Donker et al. [18].
For the sake of simplicity, here we consider post-mitotic cells (gray in Fig. 1) and cells in G0/G1, as one single cell population.Quantifying post-mitotic cell density is crucial in order to estimate both K 1 and K 2 in Eqs.(2), given that these parameters are measures of contact inhibition of proliferation, typically associated with regions of higher cell density [28].
To calibrate the model, we fit to the estimated cell density obtained by averaging eleven experimental realisations.We show the univariate marginal posterior distributions corresponding to the model parameters in Fig. 2A, confirming that all model parameters are practically identifiable.In particular, all marginal posteriors show well-defined and unimodal distributions, with a relatively narrow variance.For more details on model calibration we refer to the Supplementary Information (Supplementary Fig. 1).
The posterior distributions in Fig. 2A are not only useful to inform further model predictions, but also give insights into the fundamental mechanisms underlying cell proliferation.In particular, given the intrinsic transition rate from G1 to S, k 1 , and the constant K 1 in Eqs.(2), we can estimate the average duration of the combined G1/post-M phase, for a given fixed density ρ, as Analogously, the estimated average duration of the S/G2/M phases is given by 1/k 2 g(ρ) = 1/(k 2 (1 − ρ/K 2 ) + ).Note, however, that these are only estimates of the timescales associated with different cell cycle stages.Put together, these estimates predict for a range of densities between 4000-4500 cells/mm 2 , a population doubling time of 14-20 hours.In Fig. 2B we plot these timescales as a function of the density ρ, observing how the duration of the different cell cycle stages increases with density.These results confirm again, in line with previous experimental measurements [2,18], that cell cycle dynamics are tightly regulated by density-dependent effects.In particular, these estimates are consistent with the experimental measurements of Donker et al. [18] -taking into account that the initial cell densities in our datasets are around ρ ∼ 3500 cells/mm 2 .At very low densities, however, our estimates predict a relatively short cell cycle duration.This suggests that the shape of the crowding functions f and g might be closer to a constant function in this regime.
In Fig. 2C we show numerical solutions of the model (Eqs.( 1) and ( 2)), taking the posterior modes as parameter values.These confirm that the model can describe cell cycle dynamics inside expanding epithelial tissues.Notably, the model captures the tissue expansion speed, as well as the S/G2/M density peak near the edge of the tissue, which results from density-dependent effects regulating the cell cycle.We also note that this type of density profile is possible in the model when crowding-dependent effects are stronger in the early stages of the cell cycle (G1/post M) and weaker in the latter ones (S/G2/M).In terms of Eqs.(2) this requires having K 1 < K 2 , which is correctly identified from the data.
We observe that the model overestimates the experimental density for early times of the experiment and, as a result of the model fit, underestimates it at later times.This is likely due to the transient behaviour that cells exhibit immediately after stencil removal [29,30], which could have an impact on cell behaviour even after the first ten hours of expansion, as suggested also in previous studies [24].However, we emphasise that tissue edge motion can be well described by the model.
A similar behaviour is reported when the model is compared to a second set of experiments performed by Heinrich et al. [19].In this case, we use the obtained parameter estimates to describe the expansion of initially smaller epithelial monolayers (initial diameter ∼ 1.7 mm).We highlight that the mathematical model can capture the expansion dynamics near the tissue edge as well as the expansion speed (see Supplementary Fig. 3), even though model parameters were inferred from the large tissue expansions.
Tissue colonisation experiments.Our model, together with the experiments of Heinrich et al. [19], reveals the intrinsic connection between tissue crowding and cell cycle progression, showcasing how this interplay can give rise to spatiotemporal patterns of cell proliferation in growing tissues.Next, we show how the model can be used to study and describe similar patterns observed in several other experimental studies using FUCCI and variants of it.
Streichan et al. [2] show, using a tissue barrier assay, how the cell cycle can be reactivated by allowing cells to migrate and colonise free space -see top row in Fig. 3A.These experiments are initialised by growing MDCK-2 FUCCI cells in the G0/G1 phase within a removable barrier.After barrier removal, the tissue quickly colonises the available space, and cells behind the barrier, which were initially in G0/G1, reactivate their cycle by entering S phase.On the other hand, cells located further behind the barrier remain at high density and do not progress through the cell cycle.
By solving numerically Eqs. ( 1) on a one-dimensional domain -see bottom row in Fig. 3A -we immediately observe how a model accounting for density-dependent regulation predicts similar behaviour to that observed experimentally 1 .In particular, and as inferred from the experimental data of Heinrich et al. [19], the calibrated model predicts that crowding-dependent effects have a greater impact at the G1-S transition, compared to the S/G2/M phases.In terms of the model and the choice of crowding functions (Eqs.( 2)), this once again requires Cell cycle regulation and cell migration.Given that assuming K 1 < K 2 seems necessary in order to obtain biologically realistic model predictions, what role does tissue crowding play in shaping cell migration patterns?We explore this question by varying the values of K 1 and K 2 in Eqs.(2) -Fig.3B.First, we observe that the density of S/G2/M peaks near the tissue edge when K 2 > K 1 > 0 and remains low in the tissue bulk as long as K 2 ≫ K 1 .However, for K 2 ∼ K 1 , the height of this peak decreases and the fraction of S/G2/M cells in the tissue bulk increases.On the other hand, when we assume a higher influence of density during S/G2/M relative to G1/post-M (K 2 < K 1 ), we observe that the tissue centre shows a higher fraction of cells in S/G2/M, in contrast with previously reported observation of contact inhibition of proliferation [27].
The numerical solutions in Fig. 3B suggest that lowdensity initial conditions lead to travelling wave solutions in one spatial dimension: ρ 1 (x − ct), ρ 2 (x − ct), with c > 0 being the wave speed, and x denoting the spatial coordinate.Standard arguments -see Supplementary Information -predict the existence of a minimum travelling wave speed in terms of only three model parameters Interestingly, this suggests that the invasion speed is inde-  1) for different values of K1 and K2 and at time points t = 50, 100, 150, 200 h.Units of K1 and K2 are cells/mm 2 .Initial conditions: ρ1(x, 0) = ρ2(x, 0) = 500 cells/mm 2 for x < 850 µm, and ρ1(x, 0) = ρ2(x, 0) = 0 cells/mm 2 otherwise.In all cases, all parameters except for K1 and K2 are fixed (taken from posterior modes).
pendent of cell cycle regulation, and only depends on cell motility (D), and the intrinsic, density-independent growth rates (k 1 and k 2 ).However, we highlight that, as shown in the figure, crowding constraints play an important role in shaping collective migration patterns.
The expression for the minimum travelling wave speed facilitates a comparison between the two-stage model proposed here (Eqs.(1)), and conventional single-population models of cell migration of the form where F is a non-increasing function satisfying The intrinsic growth rate of the population, r, is related to the intrinsic rates of cell cycle progression, k 1 and k 2 , via r 3) can be approximated by which agrees with the prediction of the well-known Fisher-Kolmogorov-Petrovsky-Piskunov (FKPP) equation (F (ρ) = 1 − ρ/K for a maximum cellular density K > 0) in one spatial dimension.Using the estimated parameter values we obtain 4r/(k 1 + k 2 ) ∼ 0.98, and in this case Eq. ( 3) predicts a minimum travelling wave speed of c min ∼ 33 µm/h, while the FKPP approximation yields c min ∼ 26 µm/h; both of them within the measured values by Heinrich et al. [19].
A better comparison with the two population model can be obtained by setting rF (ρ) = λ(ρ), where λ(ρ) is the dominant eigenvalue of the growth matrix , as given by Eqs.(1).In this case, , agreeing with the prediction from Eq. (3).
More generally, we noted that c min does not depend on the choice of crowding functions f and g; however, crowding constraints have an impact on the observed migration patterns (Fig. 3B).To understand how growth and cell cycle regulation lead to the patterns observed experimentally, we investigate travelling wave solutions in a simplified version of our model, utilising the same parameters as in Eqs. ( 1) and ( 2) (see Supplementary Information).In particular, we set f (ρ) = H(K 1 − ρ), and g(ρ) = H(K 2 − ρ), where H(•) denotes the Heaviside function.This reduced model does not accurately approximate the model presented in Eqs. ( 1) and ( 2), but nonetheless it captures the same qualitative < l a t e x i t s h a 1 _ b a s e 6 4 = " B 4 e q U b 9 < l a t e x i t s h a 1 _ b a s e 6 4 = " z o t f t 7 z 6 O T v p A R A 4 N J q V d u R 7 N 6 w = " > A A A C A X i c d V D L S g N B E J y N r x h f U S + C l 8 E g e A q 7 q 9 H k I A S 9 e I x g H p D E M D u Z T Y b M P p j p F c M S L / 6 K F w + K e P U v v P k 3 z i Y r q G h B Q 1 H V T X e X E w q u w D Q / j M z c / M L i U n Y 5 t 7 K 6 t r 6 R 3 9 x q q C C S l N behaviour, and hence we expect that the relevant phenomena show similar dependence with respect to model parameters (see Supplementary Fig. 4).The analysis of travelling wave solutions for this simpler model suggests that the density of S/G2/M cells in the tissue bulk, ρ bulk

2
, only depends on the ratio of cell cycle progression rates, κ = k 1 /k 2 , and on the ratio of densities associated with crowding constraints, K 1 /K 2 , (Fig. 4).In particular, we obtain where α(κ) = 2/( √ κ 2 + 6κ + 1−κ−1).A similar dependence with respect to the model parameters is observed numerically for the model given by Eqs. ( 1) and (2) (Supplementary Fig. 4).For our estimated parameters, the expression above predicts ρ bulk 2 /K 2 ∼ 0.3, which is consistent with experimental observations.We also highlight that, as long as K 1 < K 2 , and k 1 and k 2 are of a similar order of magnitude, this expression predicts that the number of cells in S/G2/M in the tissue bulk will be small in comparison to the number of cells in G1/post-M (Fig. 4B).In particular, note that ρ bulk Interestingly, the travelling wave analysis also reveals that, when ρ bulk 2 > 0, the difference in S/G2/M cell density between the tissue bulk, ρ bulk 2 , and the tissue edge, ρ edge 2 , depends only on the difference of densities associated to crowding constraints at the G1-S and G2-M boundaries (Fig. 4C), where For our estimated parameters, we obtain ρ edge 2 − ρ bulk 2 ∼ 500 cells/mm 2 , again consistent with the experimental observations.These analytical expressions confirm the impact of density-dependent effects on cell migration and suggest that differences in the regulation of cell cycle stages contribute to the emergence of cell proliferation patterns.

Density-dependent effects and experimental design.
We have demonstrated that our model can be calibrated to experimental data collected by Heinrich et al. [19] to provide confident estimates of all parameters and, from there, used to extract and quantify crowding constraints regulating the cell cycle.An obvious question to ask is whether the model parameters could also be confidently estimated from other datasets, in particular where the cell density remains much lower and the impact of tissue crowding is reduced.To explore this question, we attempt to estimate the model parameters (including K 1 and K 2 ) using data from a lowdensity scratch assay with 1205Lu melanoma cells -see Fig. 5.We highlight that melanoma cells are highly metastatic and often display uncontrolled and invasive migration, in contrast to the highly collective and regulated movement exhibited by epithelial cells.In this experiment, tissues are seeded at an initial density of ∼ 400 cells/mm 2 (5% of the theoretical maximum packing density [7]), and data is collected every 16 hours, over two full days, allowing cells to potentially undergo 1-2 cell cycles.The posterior distributions obtained for the different model parameters reveal estimates for D, k 1 and k 2 that are consistent with previous studies [8].However, the low experimental densities do not allow for the quantification of density-dependent effects -the parameters K 1 and K 2 cannot be estimated with any degree of confidence (see Supplementary Fig. 5).This non-identifiability of K 1 and K 2 suggests the use of a simpler model, which assumes that cell cycle progression is independent of density-dependent effects (f (ρ) = g(ρ) = 1) and hence is only valid in the low-density regime.Indeed, when calibrated to data from the low-density scratch assay it x (µm) Figure 5: Absence of densitydependent effects in a low-density scratch-assay experiment (1205Lu melanoma cells).In this case, an exponential growth model can reproduce the experimental data.Density is normalised by using the theoretical maximum density corresponding to hexagonal close packing of cells [7].Top row is adapted from [8], with scale bars corresponding to 200 µm.Numerical solutions of Eqs.
provides accurate parameters estimates (Supplementary Fig. 6), and an excellent agreement with the experimental data -see Fig. 5.This result clearly illustrates both the key role that mathematical modelling can play in the experimental design process, and the importance of considering parameter identifiability in the process of model construction.

Discussion
In this work, we have presented a new mathematical model of cell migration with cell cycle dynamics which captures and quantifies cell cycle regulation by sensing levels of tissue crowding.In line with previous experimental studies, by combining minimal modelling and Bayesian inference, we confirm that cell cycle progression is monitored via crowding constraints [2,18], and present a systematic approach towards the quantification of interactions regulating cell proliferation.Our model is capable of quantifying cell cycle data from experiments using the FUCCI system, and enables the extraction of mechanistic insights into how individual cells regulate proliferation based on population-level measures.
The model presented here offers several applications to further our understanding of cell-cell interactions in cell proliferation.In particular, our model presents a systematic way to quantify the impact of drugs and gene knockouts/knockdowns interfering with cell proliferation.By using parameter estimation techniques, applied to different experimental datasets, we can gain insights into the regulatory roles of specific genes in the cell cycle.Another possible application concerns the study of cell migration in biomaterials incorporating cadherin proteins, which have recently been shown to slow down cell cycle dynamics [20].Furthermore, generalisations of the FUCCI system could allow for a finer representation of the different cell cycle stages -for instance, FUCCI4 [12] allows for the simultaneous visualisation of the four stages of the cell cycle.In line with these methodologies, extensions of our model (Eqs.( 1)) to multi-stage cell populations are straightforward, and could enable a more exhaustive explanation of the role of spatial constraints across all four cell cycle stages [18].
In the case of Heinrich et al.'s experiments [19], the excellent imaging quality allowed us to perform an accurate quantification of the cellular density profiles.This, in turn, facilitated model development and the subsequent inference of model parameters from the data, with the estimated parameters showing a low uncertainty.While the parameter identifiability of such mathematical models can be evaluated a priori under the assumption of infinite ideal data [31,32], biologically realistic datasets are finite, and often contain a significant level of noise, which can, in certain instances, constrain the ability to confidently estimate model parameters.More generally, and as we have illustrated, practical constraints in the experimental data often relate to the level of model complexity which can be inferred from experiments and the confidence in model parameter estimates.
Continuum models are a widely adopted approach for describing cell migration.However, these models come with limitations: they tend to neglect local structure, especially in situations involving multiple cell populations.Such local structure can be observed in Fig. 1; (C) and (D) show some degree of local correlation in the cell phases, however this phenomenon is lost when averaging radially to obtain the density profiles in (E).Agent-based models [21,33,34] can help mitigate some of these issues, by providing more understanding of the generation and maintenance of spatial structure, but at the cost of increased computational times for simulation and inference, additional model parameters, and limited analytical tractability.We emphasise, however, that cell cycle dynamics appear to be globally desynchronised, as observed in previous studies [35], and so our differential equation-based model remains appropriate for this study, where the data is generated by averaging over a number of experimental replicates.
The model presented here is minimal in the sense that it assumes that cell movement is random, and it ignores basic cell-cell interactions which are typical of epithelial cell migration such as cell-cell adhesion.While local cell density is likely to have an impact on cell motility [19], previous work shows that for individual expanding epithelial tissues, the linear diffusion model provides a good approximation [24].Note, however, that it is important to account for population pressure and its impact on cell movement when considering tissue-tissue interactions [36].Additional research is needed to determine whether more complicated models [37,38], incorporating cell-cell adhesion and other basic interactions offer deeper mechanistic understanding.Moreover, the model given in Eqs. ( 1) assumes that at low densities, the duration of each of the cell cycle stages follows an exponential distribution.While this assumption contradicts experimental observations [39,40] and can be mitigated by representing the cell cycle as a multi-stage process [9,41], such models break the cell cycle into a very large number of stages, limiting the potential for calibration to experimental data.Additional investigation is required to understand the extent to which more complicated models can provide further insights into how cells coordinate proliferation and migration to give rise to complex collective behaviours.For example, in the context of the cell cycle, an option is to explicitly incorporate cell cycle stage via the use of an age-structured model [42] that includes density-dependent regulation.Our results indicate that adopting a quantitative approach [43], that carefully examines quantitative data through the lens of mathematical modelling and Bayesian inference, can help provide answers to this question.

S1 Bayesian parameter estimation
All experimental datasets [2,8] consist of direct measurements of the density of cells in the G1/post-M, and S/G2/M phases of the cell cycle.We denote these measurements by {ρ D 1 (x i , t j ), ρ D 2 (x i , t j )} i,j .The next step in order to estimate the different model parameters is to assume a so-called error model, which relates experimental measurements with the model predictions given by the solutions of the model: ρ 1 (x, t), ρ 2 (x, t).For simplicity, here we assume that the residuals are independent and normally distributed where σ 1 and σ 2 are parameters to be estimated from the data.The white noise assumption has the advantage that simple likelihood-based methods can be used for inference.In particular, the log-likelihood of observing the data, given specific model parameters θ, can be written as While we assume a simplistic noise model to perform the parameter inference, model misspecification and the temporal resolution of the measured data are likely to introduce correlations between residuals [3].In particular, some degree of correlation might be expected given the model parameters are very well-determined with a relatively small variance -see Fig. 2 in the main text.More recently, a binomial measurement error model has been suggested in order to mitigate some of the inconsistencies of the white noise assumption [9].Other commonly used error models assume different forms of multiplicative noise, which preserve the positivity of the data [5,6].A more comprehensive study of the error model is left as a subject for future investigation.
Maximising the log-likelihood function would give a set of parameters θ * that we could use to generate further model predictions.However, this approach does not give any information on the associated uncertainty, which here is of particular interest given that the parameters are estimated from noisy experimental data.To explore parameter identifiability, we follow a Bayesian approach, in which uncertainty associated with model parameters (θ) is quantified in a posterior distribution . This posterior distribution can be calculated from Bayes' theorem where P (ρ D | θ) = exp D (θ) is the likelihood of observing the measured data, and π(θ) is the prior distribution of the parameter vector θ.For the tissue expansion experiments [2], we assume a log-uniform prior on D, k 1 , k 2 , with bounds: 10 1 µm 2 /h < D < 10 4 µm 2 /h, This assumption allows us to consider a broad range of orders of magnitude for these parameters, although simpler uniform priors could also be used.For the parameters K 1 , K 2 , σ 1 , σ 2 , uniform priors with the following bounds were used: For the scratch assay data, we follow [8] and assume σ 1 = σ 2 = σ, and a uniform prior in all model parameters with the following conservative bounds: We use a Metropolis-Hastings MCMC (Markov chain Monte Carlo) sampler with adaptive proposal covariance to infer the posterior distributions.This is implemented in the parameter estimation toolbox pyPESTO [7].In the MCMC algorithm, a Markov Chain starts at position θ and accepts a potential move to θ * with probability q = min{1, P (θ | ρ D )/P (θ * | ρ D )}.In this way, the Markov chain tends to move towards high values of the posterior distribution, while still allowing for transitions to regions of lower probability in order to move away from local maxima.Figures S1 and S4 show typical MCMC iterations for both sets of experimental data and the corresponding data.We show the obtained stationary posterior distributions in the main text, and in Fig. S4.

S2 Outline of the numerical scheme
We briefly explain the numerical scheme used to solve our model in polar coordinates.For the tissue expansion experiments we assume solutions with radial symmetry: ρ 1 (x, t) = ρ 1 (r, t), ρ 2 (x, t) = ρ 2 (r, t), where r denotes the distance from the tissue centre.Hence, we can write We use a finite-volume scheme [2] and discretise the domain into a small circle C 0 of radius r 1/2 = δr/2, and concentric annuli C i with inner radii r i−1/2 = (i − 1/2)/δr, for i = 1, 2, . . ., N .If x ∈ C i , we approximate where |C i | denotes the volume of C i .In particular, by integrating the equation for ρ 1 over C 0 we obtain where |C 0 | = πr 2 1/2 denotes the area of C 0 .The first integral can be calculated exactly to obtain The last two integrals can be approximated to obtain where ρ 0 = ρ 0 1 + ρ 0 2 .Similarly, we integrate the equation for ρ 1 over C i , i ≥ 1, to obtain where An analogous set of equations can be obtained for ρ 2 following the same arguments.Finally, we approximate the derivatives ∂ r ρ as We solve the resulting differential equations using a fourth-order Runge-Kutta method implemented in the scipy.integrate.odeclass in Python.

S3 Minimum travelling wave speed
We look for travelling solutions in the model given by Eqs.(2) in the main text, in one spatial dimension.We assume that the crowding functions f (ρ) and g(ρ) are non-increasing with ρ, and non-negative.In the comoving reference frame, we can write: ρ 1 (x, t) = U 1 (z), ρ 2 (x, t) = U 2 (z), where z = x − ct, and c ≥ 0 denotes the wave speed.By denoting where the primes indicate differentiation with respect to z.
The set of steady states of system (1) consists of the origin (U 1 , V 1 , U 2 , V 2 ) = (0, 0, 0, 0) and any state of the form (α, 0, U * − α, 0), with f (U * ) = g(U * ) = 0 and 0 ≤ α ≤ U * .Note that whenever f (ρ), g(ρ) > 0 for all ρ ≥ 0, the latter does not exist.As usual with linear diffusion models, the stability of the origin gives a lower bound on the wave speed c.In particular the Jacobian of system (1) at the origin reads The eigenvalues λ i of the linearized system about this point satisfy the polynomial equation By defining the roots of this quartic polynomial can be expressed as We seek biologically realistic solutions with U 1 , U 2 ≥ 0, and hence the eigenvalues must be real.In particular, this demands γ ± ≥ 0, which establishes the minimum travelling wave speed found in [10] By writing we observe that when 4k 1 k 2 /(k 1 + k 2 ) 2  1, the minimum travelling wave speed can be approximated by which agrees with the minimum speed predicted by the Fisher-Kolmogorov-Petrovsky-Piskunov (FKPP) equation [4].

S4 Study of travelling wave solutions
A commonly used approach to obtain approximate solutions for travelling waves is the socalled Canosa's method [1].This procedure is a standard singular perturbation technique, and consists of a transformation y = −z/c, where D/c 2 := ε is treated as a small parameter.The first-order perturbation in ε approximates, within a small error, travelling solutions of the well-known FKPP equation, even though in this case ε is not necessarily small [4].
In our case, by using the esimated parameters and a wave speed of 30 µm/h, we obtain ε ∼ O(1).We highlight, however, that the lowest order approximation in ε provides an excellent approximation of the travelling wave -see Fig. 1.By using the transformation y = −z/c, we can write system (1) as Observe that, given the sign of the transformation y = −z/c, we need to impose the following  1) approximation (dashed lines), obtained from solving the ordinary differential equations ( 5) and (6).Model parameters are taken from posterior distribution modes.Although the analysis as ε → 0 looks like a singular perturbation problem, setting ε = 0 gives a valid first-order approximation.This is due to the fact that the nonlinear terms in Eqs. ( 3) and (4) vanish at both boundaries [4].Hence, we can look for a regular perturbation expansion in both U 1 and U 2 .By denoting the order O(1) solutions as u 1 and u 2 we obtain where u = u 1 + u 2 .In the figure above (Fig. 1), we compare the approximate solutions obtained by solving this system with the full travelling wave solutions.We highlight that the lowest order approximation provides an excellent approximation of the travelling wave shape.

S4.1 A simplified model
In order to make analytical progress we set f and g to be Heaviside functions: f (u) = H(K 1 −u) and g(u) = H(K 2 −u).This model is not an approximation of the model presented in the main text, but a simplification which preserves the same qualitative behaviour.Hence, we expect that the observed phenomena show similar dependence on the model parameters; this will be numerically confirmed later.In particular, note that this simplified model also describes two density checkpoints, at the G1-S boundary, and during the G2/M phases.The parameters K 1 and K 2 , respectively, quantify the cell density associated with these checkpoints.We rewrite Eqs. ( 5) and.(6) in terms of the variables (u, u 2 ), du dy = k 2 u 2 g(u) , Depending on the relative values of the total cell density, u, and the density checkpoints parameters K 1 , K 2 , we distinguish three possible cases.As inferred from the experimental data, we assume K 1 < K 2 .
t e x i t s h a 1 _ b a s e 6 4 = " D Z H t D 5 S i r o Y X 8 G F M O j F E 9 F J u + 3 k = " > A A A B 8 X i c b V B N S w M x E J 3 4 W e t X 1 a O X Y B H q p e y K q M e i F 4 8 V 7 A e 2 S 8 m m 2 T Y 0 m y x J V i h L / 4 U X D 4 p 4 9 d 9 4 8 9+ Y t n v Q 1 g c D j / d m m J k X J o I b 6 3 n f a G V 1 b X 1 j s 7 B V 3 N 7 Z 3 d s v H R w 2 j U o 1 Z Q 2 q h N L t k B g m u G Q N y 6 1 g 7 U Q z E o e C t c L R 7 d R v P T F t u J I P d p y w I C Y D y S N O i X X S 4 6 j n R 5 W u H q q z X q n s V b 0 Z 8 D L x c 1 K G H P V e 6 a v b V z S N m b R U E G M 6 v p f Y I C P a c i r Y p N h N D U s I H Z E B 6 z g q S c x M k M 0 u n u B T p / R x p L Q r a f F M / T 2 R k d i Y c R y 6 z p j Y o V n 0 p u J / X i e 1 0 X W Q c Z m k l k k 6 X x S l A l u F p + / j P t e M W j F 2 h F D N 3 a 2 Y D o k m 1 L q Q i i 4 E f / H l Z d I 8 r / q X V f / +o l y 7 y e M o w D G c Q A V 8 u I I a 3 E E d G k B B w j O 8 w h s y 6 A W 9 o 4 9 5 6 w r K Z 4 7 g D 9 D n D 5 4 W k D s = < / l a t e x i t > k1f (⇢) < l a t e x i t s h a 1 _ b a s e 6 4 = " + D x P G e G M b a b 2 p r / 4 b b / 4 b s + 0 e t P X B w O O 9 G W b m B T F n 2 r j u t 1 N Y W 9 / Y 3 C p u l 3 Z 2 9 / Y P y o d H b S 0 T R W i L S C 5 V N 8 C a c i Z o y z D D a T d W F E c B p 5 1 g c p v 5 n S e q N J P i w U x j 6 k d 4 J F j I C D Z W e p w M 6 q N o l y 7 y e M o w D G c Q A V 8 u I I a 3 E E d G k B B w j O 8 w h s y 6 A W 9 o 4 9 5 6 w r K Z 4 7 g D 9 D n D 5 4 W k D s = < / l a t e x i t > k1f (⇢) < l a t e x i t s h a 1 _ b a s e 6 4 = " + D x P G e G M b a b 2 p r / 4 b b / 4 b s + 0 e t P X B w O O 9 G W b m B T F n 2 r j u t 1 N Y W 9 / Y 3 C p u l 3 Z 2 9 / Y P y o d H b S 0 T R W i L S C 5 V N 8 C a c i Z o y z D D a T d W F E c B p 5 1 g c p v 5 n S e q N J P i w U x j 6 k d 4 J F j I C D Z W e p w M 6 q N B D e r A 4 B G e 4 R X e P O 2 9 e O / e x 7 x 1 x c t n j u A P v M 8 f S X K O 8 Q = = < / l a t e x i t > ⇢ 1 l e 4 Q 0 p 9 I L e 0 c e i t Y D y m W P 4 A / T 5 A 0 r 2 j v I = < / l a t e x i t >

Figure 1 :
Figure 1: (A) Schematics of the FUCCI cell cycle marker system and model conceptualisation.Transitions in the model given by Eqs.(1) are regulated by the crowding functions f (ρ) and g(ρ), dependent on the total cell density ρ = ρ1 + ρ2.(B) FUCCI fluorescence images from the experiments of Heinrich et al. [19] at different time points (adapted).Initial tissue diameter ∼ 3.4 mm.Scale bars correspond to 1 mm.(C) Segmented data showing G1 (red), S/G2/M (green), and post-mitotic (gray) cells.Note that in the model we combine post-mitotic cells and cells in G1. (D) Zoomed-in segmented data at the tissue edge and centre, corresponding to the black squares in (C).(E) Fraction of cell-cycle state cells in the tissue centre and in the tissue edge -defined as regions extending ∼ 200 µm from the tissue center and tissue edge, respectively.(F) Density profiles in polar coordinates at t = 40 h, showing cells in G1, S/G2/M, and post-mitotic cells.(E) and (F) show the average of eleven independent tissue expansions with the same experimental initial condition, with shaded regions indicating one standard deviation with respect to the mean.

Figure 2 :
Figure2: Density-dependent effects regulate cell cycle dynamics in epithelial tissue expansion experiments[19].Parameter estimation and model-data comparison for the model given by Eqs.(1) and 2. (A) Univariate marginal posterior distributions for the model parameters.Posterior modes are given by (D, k1, k2, K1, K2) = (1300 ± 66 µm 2 /h, 0.612 ± 0.015 h −1 , 0.457 ± 0.011 h −1 , 4965 ± 38 cells/mm 2 , 5435 ± 45 cells/mm 2 ), where errors correspond to one standard deviation.(B) Estimated duration of the G0/G1/post M (red) and S/G2/M (green) phases, as well as the whole cell cycle (blue), as a function of cell densities.Solid lines correspond to posterior modes and shaded regions are obtained sampling from the posterior distribution.(C) Comparing data and model predictions.Squares represent the estimated cell density obtained by averaging eleven experimental realisations, which we use to calibrate the model.Shaded regions denote one standard deviation with respect to the mean -see Supplementary Fig.2for confidence intervals in the model predictions.Numerical simulations in polar coordinates were obtained by using the posterior modes as parameter values, and no-flux boundary conditions -for details on the numerical scheme we refer to the Supplementary Information.In order to minimise the effects of the stencil removal on cell behaviour, the initial condition corresponds to the experimental density profile ten hours after stencil removal.

2 )
4 y 3 T 5 h 0 l v 8 = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " D 1 3 n b 9 1 4 9 w w o n Q m C G l u o 4 d a y 9 F U l P M y L T Q S x S J E R 6 h A e k a F I g T 5 a W z o 6 f w y D h 9 G E b S P K H h z P 0 + k S K u 1 I Q H p p M j P V S / a 5 n 5 V 6 2 b 6 P D c S 6 m I E 0 0 E n i 8 K E w Z 1 B L M E Y J 9 K g j W b G E B Y U n M r x E M k E d Y m p 4 I J 4 e u n 8 H 9 o u 2 X n r O z c n J b q l 4 s 4 8 u A A H I J j 4 I A a q I N r 0 A Q t g M E 9 e A B P 4 N k a W 4 / W i / U 6 b 8 1 Z i 5 l 9 8 E P W 2 y f / x p D 0 < / l a t e x i t > ⇠ K 2 K 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " 4 1 1 z J 1 n i z X N O x E u 6 1 z 5 P z S 8 K I Z A = " > A A A B / X i c d Z D J S g N B E I Z 7 X G P c x u X m p T E I n s J M Y o z e g l 4 8 R j A L Z M a h p 9 O T N O l Z 6 K 4 R 4 x B 8 F S 8 e F P H q e 3 j z b e w s g o r + 0

Figure S1 : 2 Figure S5 :
Figure S1: MCMC iterations for the large tissue expansions experimental data.Parameters D, k 1 , k 2 , K 1 , K 2 correspond to the model presented in the main text, and σ 1 , σ 2 are error model parameters.