The role of the cell cycle in collective cell dynamics

Cells coexist together in colonies or as tissues. Their behaviour is controlled by an interplay between intercellular forces and biochemical regulation. We develop a simple model of the cell cycle, the fundamental regulatory network controlling growth and division, and couple this to the physical forces arising within the cell collective. We analyse this model using both particle-based computer simulations and a continuum theory. We focus on 2D colonies confined in a channel. These develop moving growth fronts of dividing cells with quiescent cells in the interior. The profile and speed of these fronts are non-trivially related to the substrate friction and the cell cycle parameters, providing a possible approach to measure such parameters in experiments.

Cell proliferation is controlled by the cell cycle, a genetic regulatory mechanism. To divide, cells must pass so-called "checkpoints" to ensure that the mother cell is properly prepared to divide. Fig. 1 (a) shows a schematic of the cell cycle, which is classified into four phases, each characterized by a set of separate events: growth and preparation for DNA replication (G1), DNA replication (S), preparation for mitosis (G2), and mitosis (M), immediately followed by cell division. There are three main checkpoints: The G1 checkpoint control mechanism ensures that the cell is ready for DNA synthesis; the G2 checkpoint ensures that the cell is ready to enter the M (mitosis) phase and divide; a checkpoint in the middle of mitosis (Metaphase Checkpoint) ensures that the cell is ready to complete cell division. Cells that have temporarily stopped dividing are said to have entered a quiescence state, also called the G0 phase. The switch into this state usually happens in * m.s.turner@warwick.ac.uk The phase of the cell cycle is parameterised by an angle θ with a parameter r a proxy for cell cycle activity, represented by a radial distance. The dynamics of r are controlled by a separate equation that is sensitive to physical stresses (see text and eq. (2)). Unstressed, the cell cycle orbits around the circle r = 1 in a space that can be thought of as a proxy for the biochemical/genetic concentrations, dividing every τ div when θ = θ div = 2nπ with n ∈ Z, just like in (a). The red line represents a (reversible) transition between a proliferating state and full quiescence at r = 0 with the rate of progress around the cycle f (r) a decreasing function of activity r, reflecting slowing of the cell cycle.
the G1 phase and is reversible, should conditions again become more favourable for growth and division. Cells that are densely packed, small and physically or chemically stressed tend to divide more slowly, corresponding to a slowing of the underlying cell cycle, sometimes to full quiescence [35,36]. The cell cycle involves an interplay between a range of genetic promoters and inhibitors, with analogues found across most eucaryotes. For example, upon receiving a pro-mitotic extracellular signal, G1 cyclins and cyclin-2 dependent kinases (CDKs) become active in preparing the cell for S phase, promoting the expression of transcription factors and enzymes required for DNA replication. In contrast, two families of genes, the cip/kip (CDK interacting protein/Kinase inhibitory protein) family and the INK4a/ARF (Inhibitor of Kinase 4/Alternative Reading Frame) family, prevent the progression of the cell cycle by binding to and inactivating the corresponding cell-cycle promoters [37]. Quiescent cells maintain a transcriptional state that is different from proliferating cells, achieved by restraining proliferation and cell-cycle progression genes [38]. A growing body of evidence suggests that quiescence is a non-terminal and tissue-specific state that can be initiated and sustained by mechanical factors, such as cell-to-cell friction and extracellular matrix (ECM) friction. Cells can sense external physical cues. Physical forces are transmitted via biochemical signaling pathways that regulate the cell cycle [22,39,40]. On the contrary, some recent experiments examining cellular protein dynamics have shown that processes such as protein expression and transcription are related to the cell cycle state. In other words, the cell cycle is not only a consequence of cell dynamics, but the cell cycle state itself influences cell behavior [41][42][43].
Underlying this complexity is a simple picture: the cell cycle can be viewed as a biochemical oscillator, producing almost regular division events if the cell is not otherwise perturbed, e.g. by lack of nutrient or physical stresses. Similar oscillators are known to produce circadian rhythms in plants [44,45] and there are also analogues in abiotic chemical systems in vitro, such as the well-known BZ reaction [46]. Such chemical oscillators require a minimum of two components, each influencing the production of the other, although in vivo there are many more than this. The course of these reaction can be visualised as a trajectory in the space of the concentrations of the various chemical species (on each axis). These trajectories form closed loops, signifying the presence of an oscillator 1 . A typical cartoon of the cell cycle, as shown in Fig 1, can also be seen as a simplified projection of these closed reaction loops. The network that controls the same process in prokaryotes is similar in function but differs in its molecular components and checkpoints, e.g. no nuclear envelope breakdown (or reformation) is necessary.
Cell biologists have been studying the cell cycle for many decades [47][48][49][50][51][52]. Historically, these studies sought to minimise the effect of mechanical heterogeneity and study the cell cycle "in isolation" as far as that is possible. This philosophy is now rapidly changing, with 1 In closed chemical systems a better description would be that of a spiral, as the reactive components are gradually depleted. a direct role for mechanics experimentally confirmed [36,[53][54][55]. Several key signalling pathways have been identified to correlate with mechanical-feedback capacity [56][57][58][59]. This is driving interest in the role of mechanical signalling, e.g. in cellular morphology, colony development and cell competition [60][61][62][63]. Alongside this there has been increasing interest in active tissues in the physics community, with the development of models of active, out of equilibrium materials that can be reminiscent of foams or soft glassy materials [25,26,38,64,65] and commonly employ continuum hydrodynamic analysis [24,[65][66][67][68][69][70][71]. Understanding the physical behaviour of the tissue is the main motivating factor and cell division is typically included in a fairly stylised form [72], if at all. The state of the art in this respect [25] incorporates cell elasticity and adhesive cell-cell interactions, as well as cell birth and death but can't really be said to incorporate any description of the cell cycle.

A. Motivation and outline
Our goal in the present work is a model that couples the physical and biochemical/genetic descriptions: the physics affects the biochemistry in the way mechanical stress slows the cell cycle, while the cell cycle affects the mechanical stresses because it controls cell proliferation and volume, in turn driving the flows that determine these stresses. The motivation for the present work is to combine a functional, rather than biochemically specific, model of the cell cycle with a realistic physical model. We aim to keep the structure of this model as simple as possible. In particular we assume (i) constant, uniform nutrient levels, (ii) waste or signalling compounds are rapidly removed/diluted and (iii) the role of both motility/migration and apoptosis (cell death) are initially neglected. Our motivation in neglecting these factors is to isolate the role of the cell cycle as far as possible in this study so that it can be studied in relative isolation. In section IV we outline how these other factors could be incorporated within the same framework to generalise our model.
We propose a minimal model that involves a cell cycle oscillator coupled to a physical model through the local stress (pressure) and cell volume. Our primary goal is to show that the cell cycle, and the parameters that underlie it, can non-trivially influence the dynamics of growing colonies, as well as the distribution of cell cycle activity, cell volume and stress within it. We propose the basic governing equations in section II which can be inferred to be reasonable starting points for a stylised model that reflect previous experimental studies: cell cycle progression control [73,74], mechanicalfeedback regulation of cell proliferation state [39,[75][76][77], and cell size regulation during cell division or 3 quiescence [12,78,79]. For simplicity we focus on colony expansion in a simple quasi 2D channel geometry and provide a mechanism for relating the speed and structure of the growth front to the parameters controlling the cell cycle. We first encode our model in a particlebased simulation in which each cell carries its own cellcycle oscillator controlling size and division events. We then compare this with a continuum analysis that offers substantial analytic insight and a number of scaling results, including an inverse square-root scaling of the front speed with the substrate friction. We obtain quantitative to semi-quantitative agreement between this continuum analysis and our simulations, validating the continuum model. We then compare our model with recent experiments on MDCK cells and find encouraging agreement with the cell area distributions and division rates, in spite of the fact that motility likely plays a role in the motion of these cells, except possibly at the lowest surface frictions.
We simulate colony expansion using a range of parameters to show how the speed and structure of the front depends on the cell cycle variables and show how this can be inverted to extract biophysical parameters from experimental measurements.

II. MODEL
We propose a mechanical-feedback cell cycle model, which in which cells can reversibly switch between proliferation and quiescence driven by local pressure. Here the progress around the cell cycle is measured by an angular variable θ, with division at θ = 2nπ with n ∈ Z, obeying with r a proxy for the cell cycle activity and take f (r) = ω 0 r for simplicity. An unstressed cell at r = 1 divides with a division time τ div = 2π/ω 0 . The cell cycle activity r is assumed to depend on the local pressure p according to and is therefore suppressed when p > 0, eventually approaching r = 0 if p ≥ p r . For simplicity we take g(r, p) = 1 τr (1 − r − p/p r ) with τ r a characteristic time on which the cell cycle responds to mechanical perturbation and where sensitivity is controlled by a reference pressure p r . There is a fixed point at r = 1 for vanishing pressure p = 0, corresponding to an unstressed cell cycle that undergoes circular orbits at r = 1 in a space that can be thought of as a proxy for the biochemical/genetic concentrations.
Except at division events the volume of any cell is assumed to be given by This encodes cellular growth with a rate assumed proportional to its rate of progress around the cell cycle ∂θ ∂t via a dimensionless growth rate c, but also that the cell volume reverts to a minimal quiescent volume V q if its cell cycle stalls, with a characteristic recovery time given by τ v . We nondimensionalise according tot = ω 0 t/(2π), τ i = ω 0 τ i /(2π),c = 2πc/ω 0 ,p = p/p r , andṼ = V /V q , hence all lengths are measured in units of V q in 2D. Using dot (˙) to indicate ∂ ∂t and then dropping all˜for convenience we haveθ = 2πr (4) Eq 6 only applies between division events and so we use the symbol v to denote the average volume, in the presence of division, as measured in our simulations and for consistency with later continuum analysis. In these dimensionless units the cell cycle is parameterised by only three variables: a growth rate c and recovery times for the cell cycle and volume, τ r and τ v respectively; the reference pressure p r will also prove important in providing a reference scale for hydrodynamic/frictional stresses, arising from flows. It's worth noting that τ r and τ v can be experimentally estimated by measuring the time length spend for a cell subjected to growthlimiting pressure to (i) stop dividing or (ii) finish shrinking to the quiescent volume [12,14]. In addition to the phenomenological description, we found the timescales can also be linked to molecular processes. For example, the way τ r leads to a lag in division response resembles the p53 tumor suppressor's activation can lead to later cell cycle arrest [80], while τ r reflects the period of inhibition of YAP/TAZ preventing uncontrolled cell volume growth [81].

A. Simulation methods
We employ a particle-based simulation in which each cell, with unique index i, carries its own local variables θ i (t), r i (t) and v i (t) with it. On division the mother cell volume is divided equally between two daughters that . (e) Growth fronts later develop at the left/right edges of the colony while the interior bulk contains smaller cells that rarely divide. (f) The instantaneous distribution of dimensionless pressure p and cell cycle activity r along the channel as shown in (e); in the bulk the pressure is falling towards the cell cycle reference (stalling) pressure driving the cell cycle activity down according to eq. (5), as the cells approach quiescence. The grey dashed lines connecting (e) and (f) mark the colony edges. In these, and all later figures unless noted otherwise, we employ a dimensionless surface friction ζ = 0.01, τr = τv = τ div , and a growth rate corresponding to a nominal average cell volume v1 = 5Vq five times the quiescent volume (see appendix for further details).
inherit their mother's cell cycle phase θ and activity r, see Fig 2a. As we will be focusing on 2D colony growth in the present work, cells are treated as elastic disks with a time-dependent natural radius to accommodate growth (and division), but the code can also be implemented in 3D with spherical cells. For cell-cell interactions, we choose a Hertzian contact mechanism in which a short-range attraction force simulates cell adhesion, and an elastic repulsive force mimics the cell-cell repulsion (see appendix for details). We performed control simulations in which we confirmed that the results did not depend on the strength of the attraction or repulsion or the functional form of the repulsion: The cells need merely be sufficiently repulsive that the system remains approximately incompressible and sufficiently adhesive that the tissue does not tear. We implement our simulations in an over-damped setting with a cellsubstrate friction that is proportional to the cell area: larger cells experience higher friction. Our simulations are mainly performed in channels resembling a strip of cells with periodic boundaries in the y-direction, expanding in the +/-x-direction. Periodic boundaries strictly equate to a cylindrical shell of cells but likely provide an excellent approximation to a planar channel, provided that the walls of the channel themselves do not provide significant friction. The symmetry direction, along x, has an open boundary, while the y-direction is given a periodic boundary condition in lieu of explicit walls to the channel. Simulations begin with one cell with an unstressed 2D volume (area) v 1 = 5V q , consistent with the nominal average volume at a corresponding growth rate of c = (v 1 − 1)τ v + log 2 (see section II B for derivation. This cell has a random initial phase, sampled uniformly from θ ∈ [0, 2π] 2 and is introduced into the centre of the channel. It then divides after time t ∼ τ div . As the cells repeatedly divide the colony expands roughly as a circle until it spans across the channel width, see Fig. 2 (b)-(d); (e) shows how the colony then develops two fronts, moving with speed s in the ±x-directions. Unless specified otherwise simulations are performed at τ v = τ r = τ div , and with a dimensionless surface friction, relating force per cell area to sliding velocity, ζ = 0.01 solved using a timestep ∆t = 10 −6 .
The pressure measurement of each cell is calculated by the virial stress tensor method [82], equivalent to summing all the inward forces and dividing by the cell area (perimeter). For more detail on methods see appendix A.

B. Continuum theory
To gain analytic insight we develop a continuum theory of this cellular material, noting that continuum models underlie many of the recent physics-based mod-5 els of cell cultures and tissues reviewed in the introduction.
Cells grow in volume during the interdivision process, as given by eq. (6), while the cell division process always acts to reduce the volume of each cell: at a division event the cell volume halves. We seek to approximate this by a process that is continuous in time. The effect of division on the mean cell volume can be incorporated into a revised version of eq. (6), where we now write the division-adjusted mean volume as v, to distinguish it from the inter-divisional volume V , The new term −r log 2 provides a mean field approximation for the reduction in mean volume due to cell division. To see this note thatv v = −r log 2 has solution v ∝ (1/2) rt with 1/r the time between volume-halving division events in units of τ div .
The mean cell volume at constant r is then given by the solutionv = 0, hence v = 1 + τ v (c − log 2)r. For an unstressed cell cycle r = 1 this defines a relationship between the unstressed division-adjusted mean volume v 1 and c so that eq. (7) becomes Eq. 8 is also used in the simulations to relate a nominal unstressed division-adjusted mean volume v 1 to growth rate c.

Continuity equation
We seek a continuity equation to establish the velocity u of the growing cell culture, treated as a continuum. To achieve this we exploit Gauss' theorem u · dS = ∇ · u dV evaluated around an infinitesimal volume dV containing dn = dV /v cells of volume v. The total rate of change of volume is the integrated outward flux of cells The change in the number of cells follows exponential growth according tȯ dn = log 2r dn in the mean field, while the average cell volume follows eq. (9). Combining these results and dividing by v dn we obtain the approximate relationship with c = log 2 + (v 1 − 1)/τ v from eq. (8).

Stokes-Darcy stress balance
Cells are assumed to experience friction with a substrate with a friction coefficient 3 ζ. In general there are also internal shear stresses, controlled by viscosities η and ξ, leading to the stress balance equation eqs. (10) and (11) are analogous to Stokes equations for a low Reynolds number fluid but here incorporate the role of cell division and a Darcy-like frictional force, characterised by ζ. The term in eq. (11) involving ∇ · u, with prefactor ξ depending on the bulk viscosity, usually vanishes for incompressible fluids [83]. In what follows, we set η = ξ = 0, assuming that the growth is substratefriction dominated.

Advected description
In the presence of the flow field u the lab-frame cell activity r and volume v are given by the advected forms of eq. (5) and eq. (9) In order to establish the dynamics we simultaneously solve for u, r, v and p using eqs. (10) to (13) 4. Channel flow (1D) For consistency with the particle based simulations we assume η = ξ = 0, noting that symmetry dictates that average motion is along the x-direction, i.e. u = ux and so shear stresses are absent in a coarse-grained theory. We seek solutions in which there is a steadystate front of dividing cells moving with constant speed s. It is convenient to transform to the co-moving frame using z = st − x, focussing on the right-moving front where (the arbitrary zero of time is chosen so that) the front is at z = 0 and the colony populates the space z > 0.
3ζ = ζ/pr, dropping the tilde; similarly for η and ξ With boundary conditions with s = u(0) shorthand for the front speed, to be determined. We solve eqs. (10) to (13) and (18) numerically, as outlined in appendix B 1. We can compute the corresponding boundary condition on the derivative of r from eq. (16) according to By L'Hopital's rule and using the boundary conditions eq. (18) this is Hence the (redundant) boundary condition Similarly for the derivative of v from eq. (17) Again using L'Hopital's rule and the boundary conditions eq. (18) we find leading to

Scaling relations
In order to better understand the structure of the growth front we can compute characteristic lengthscales, at least at the scaling level, by analysing the rate of change of the mechanochemical variables near the front.
1. The pressure at the front obeys p (0) = ζs from eqs. (15) and (18). A scaling estimate of the distance behind the front at which the pressure reaches the cell cycle reference (stalling) pressure is R p = 1/p (0), hence 2. The cell cycle activity obeys r (0) from eq. (21). A scaling estimate of the distance behind the front that the cell cycle will stall to r = 0 is R r = 1/|r (0)|, hence 3. The cell volume obeys v (0) from eq. (24). A scaling estimate of the distance behind the front at which the volume reaches the quiescent volume is 4. The local lab-frame cell velocity obeys u (0) = − log 2 from eqs. (14) and (18). A scaling estimate of the distance behind the front at which the cells become immobile is R u = s/|u (0)|, hence Where s can be calculated self-consistently or estimated by setting R u ∼ R p (say), revealing the scaling relations (that then holds for all variables) and III. RESULTS

A. Substrate friction
The friction forces experienced by cell colonies expanding along a 2D channel can be traced primarily to transient adhesion with the substrate [10,28,30,84]. While it remains difficult to measure the absolute substrate friction it is more feasible to make controlled variations, e.g. by preparing surfaces with different area fractions of surface coatings. Such an approach may allow a test of eq. (30) and ultimately may give a way to measure the substrate friction indirectly, e.g. by observing the spreading speed of pre-calibrated cell lines.  Fig 3a) shows the relationship between front velocity and surface friction, revealing fully quantitative agreement between our particle-based simulations and continuum theory and confirming the scaling result eq. (30) . Fig 3 also shows (b) the crossover from exponential growth to constant front speed and (c) an example of how the speed also depends on the cell cycle parameters.

B. Colony structure
During colony expansion, the organisational structure that emerges behind the growth front depends on the mechanochemical properties of individual cells, including the parameters controlling their cell cycle. We examine the steady state distributions of cell area, outward velocity, pressure, and cell cycle activity in a frame of reference co-moving with the front in which cells are active, growing and moving near the front and quiescent and stationary deep in the bulk, see Fig 4 and Supplemental Material at [URL will be inserted by publisher] for movies of the growth process at various surface frictions and cell cycle times. Fig 4 (a) shows that the stochastic variation in the velocity is larger than in the other variables. This may be associated with the relatively large velocities that are realised, immediately post-division. Two daughter cells have large overlap immediately after division, and the repulsive force drives rapid separation. Although the post-division velocity is large this has almost no effect on the average velocity: the daughters separate in opposite directions. We also predict the emergence of cells for which the sign of the velocity is reversed (located around z ≈ 100 for these parameter values). This indicates that cells can sometimes move inward, away from the front, and is a result of slow volume loss of the cells approaching quiescence "sucking" cells in this region away from the front. This phenomenon can only be realised within relatively sophisticated models in which delayed volume loss is encoded into the cellular response. Similar negative velocities have been observed in experiments on multicellular spheres [85]. Fig 4 (b) shows the pressure profile. Two features are worth highlighting. Firstly, the pressure in the bulk is systematically elevated slightly above p r (p = 1 in dimensionless units). This is due to the controlled loss of division activity, which means some cells continue to divide even when the material is near the stalling pressure, resulting in an overshoot. A consequence of the overshoot of the pressure is that the gradient of the pressure reverses its sign in the colony interior. This would provide a signature of contractility and is due to the shrinkage of the cells in this region. Secondly, two types of pressure distribution can arise. When the cell cycle volume adaptation time τ v is short enough the pressure  monotonously increases from the edge toward the bulk, stabilising near p r . Conversely, when τ v is large the pressure is non-monotonic, exhibiting a pressure peak near the front, with the pressure decreasing further into the bulk. This is a consequence of relatively rapid volume loss of cells behind the growing front. Both of these features are quite generic. There are hints that a pressure peak might also be observed in experiments on multicellular spheroids [86], although experimental measurements of pressure remain challenging. Fig 4 (c) and (d) shows that both the cell cycle activity and cell volume decrease monotonically, approaching single exponential decay deep in the bulk (in the case of the volume decaying to the quiescent value v = 1). Similar exponential decay emerges for u and p but this is not shown because of minor complications: the velocity can go negative and the bulk pressure is not fixed, see appendix B 2 for details.
For all mechanochemical variables there is good semiquantitative agreement between the continuum theory and the particle based simulations.

C. Cell cycle control of the growth front
The cell cycle parameters τ v and τ r affect the physical distribution of the mechanochemical variables across the growth front. In order to analyze this we define a half-decay length for each variable written as z u , z p , z r , z v respectively. These are defined as the smallest root of the following relations: For r and v we use r(z r ) = 0.5 and v(z v ) = (v 1 − 1)/2; for p and u we use p(z p ) = (p max − p min )/2, involving the empirically determined maximum and minimum values of pressure, similarly for u. See appendix A for details. Fig 5 shows the relationship between half-decay lengths with the two cell cycle times τ r and τ v of velocity u, pressure p, cell cycle activity r and cell area v respectively. This relationship can be inverted to relate the experimentally observable half-decay lengths to the underlying cell cycle parameters, e.g. τ r and τ v , see Fig 6. We choose to focus on the volume and velocity decay lengths here because they can be more easily visualised.

D. The cell cycle controls features inaccessible to purely "physics based" models
In this section we review how the incorporation of a stylised cell cycle introduces descriptive power not readily accessible to a phenomenological "physics only" model.
While modelling mechanical-feedback in cell dynam-ics is not new it is known that cellular behavior, very generally, is controlled by the cell cycle, a profoundly out-of-equilibrium chemical oscillator/sensor [51,52], see also the discussion in the introduction. This regulates the cell proliferation state and is involved in other decisions the cell makes about its development [36] except possibly under the most extreme stress when there is catastrophic (genuine physical) failure. A common understanding in cell biology is that the cell cycle is the master controller with physical forces providing input to this controller [50,[87][88][89][90]. This controller can also be manipulated in numerous other ways, e.g. temperature, chemical composition and genetic manipulation. It becomes increasingly contrived and clumsy to incorporate such factors into models that start from a purely physical description but relatively straightforward with (refinements to) models that incorporate a cell cycle explicitly. For these reasons we believe that it is important that models are able to reflect the underlying regulatory mechanisms of cell biology. However, the way in which a phenomenological "physics only" model differs from one with an explicit cell cycle regulation deserves to be addressed directly. The most obvious signature of this difference here is the fact our results depend on the timescales τ r and τ v that enter through the cell-cycle eqs. (5) and (6). These affect the emergent physical behaviour, including the front speed, velocity distribution, volume distribution, activity distribution and pressure distribution, as shown in Figs 3, 4, 5 and 6. Some of the qualitative features that are controlled by these parameters include the pressure overshoot and velocity sign change that are shown in Fig 4. In order to seek to underline this we identify a limiting case of our model, corresponding to a "physics only" limit. We take a fast response limit of our cell cycle oscillator in which τ r = τ v = 0; neither parameter then appears explicitly. Equations (5) and (9) then reduce to r = 1 − p and v = 1 + r(v 1 − 1), with r and v slave to the instantaneous p and r, in the spirit of direct physical control. eqs. (14) and (15) then reduce to u = −cr and p = ζu. Substituting the former into a derivative of the later we obtain p = cζ(p − 1). Hence p = 1 − e −z/l with l = 1/ √ cζ. All other variables (r, u, v) are similarly mono-exponential, with the same lengthscale. This universal mono-exponential behaviour is clearly quite distinct from that arising under more general values of τ r and τ v , as can be seen from Fig 4, e.g. the decay lengths vary and the pressure and velocity can be non-monotonic.
Finally, it is also well known that driven oscillators have distinct properties from non-oscillatory systems [91]. Given that a genuine oscillator is present in our model one might also anticipate improved descriptive power, e.g. under conditions where the system is driven in a time-varying fashion with a refined f (r, p, θ), introduced to capture the phase (checkpoint) sensitivity of physical stimulus.

E. Dynamics of area distribution
We compare our model with experimental data for the growth of 2D sheets of MDCK cells [12], see Fig 7. Our simulations show broadly similar trends in cell area distributions and division rates to those observed in MDCK cells, noting that these cells are known to be motile and so this level of agreement is unexpected. Rather surprisingly, we were unable to locate suitable experimental data on non-motile cells in 2D that are in the non-exponential growth phase and not limited by other factors, e.g. nutrients. That we provide motivation for such experiments in the future should be considered an objective of this work. Fig. 7 (b) shows the fastest division rate for the largest cells to be approximately τ div = 0.5 day. Hence we take τ r = τ v = τ div as nominal values for the cell cycle times and fix V q = 30µm 2 and v 1 = 15V q .
In our simulations the cells adjust their size naturally according to internal cell cycle activity r and hence local pressure p. The simulations are initialised with a single cell in the centre of 2D substrate, generating a roughly circular growth front for all times, see Supplemental Material at [URL will be inserted by publisher] for a movie of colony growth in this geometry. This is to align with the experiments and is different to the 2D channel considered in the rest of this study. As a steadily expanding colony emerges the cells in the bulk (interior) start switching to the quiescent state, with an area that approaches V q . As a result the median area of the colony decreases towards V q in the late stage of simulations, as Fig. 7 (a) and (b) both show. Time t = 0 in panel (b) is measured from a "morphological transition", close to the crossover to subexponential growth. The experiments reported cells that that were thinner in the early stages of growth, with the cell thickness roughly constant once the median area fell below 100 − 150µm 2 . This means the system only has a quasi 2D nature, with area a proxy for volume, after this point. Hence our theory can only reasonably be expected to apply in this regime.

Time/days
τ v = τ div , τ r = τ div Figure 7. Cell area distributions in MDCK cells are reproduced in our simulations. The three panels on the left (a,c,e) are simulation results; the panels on the right (b,d,f) are experimental results on spreading MDCK cells, reproduced with permission from [12]. Panels (a,b) show the evolution of colony-average median cell areaÃ as a function of time. In panel (a) we set t = 0 as the crossover from exponential to subexponential growth of the colony area, roughly 20τ div after initialisation. In panel (b) time t = 0 is measured from a "morphological transition" that is close to the crossover to subexponential growth, noting that substantial changes in cell thickness were reported before this time, see [12] for details. These thickness changes mean that comparison with quasi 2D simulations (a) only becomes appropriate for t 0. Panels (c,d) show the colony-wide probability distribution of individual cell area at different times (colors). Panels (e,f) show the distribution of the time between cell division events as a function of the premitotic area. In (e) every cell division event contributes a data point while in (f) data points with division times greater than 3 days are inferred indirectly. In (f) the five different datasets (colours) are from different experiments. See [12] for details. We choose simulation parameters to roughly match those of the MDCK cells: Vq = 30µm 2 , v1 = 15Vq and τ div = τr = τv = 0.5 day, see appendix A 5 for details.
We also record the distribution of cell areas at different times. Fig. 7 (c) and (d) both show the cell areas converging to a narrow, stationary size distribution centered around 30 − 35µm 2 . Because of the way we have incorporated cell cycle activity into our model the division time for individual cells can vary enormously, depending on local pressure, as both Fig. 7 (e) and (f) show. This leads to a range of division times from τ div up to ∼ 50τ div with large (unstressed) cells dividing faster than smaller cells. This degree of agreement with the experimental data would not be produced by models that do not incorporate a similar mechanism to control division rate and size.

IV. DISCUSSION
We have developed a model for confluent cells that incorporates a stylised cell cycle, regulating division and cell size. This both incorporates and, indirectly generates, mechanical feedback: the former via pressure sensitivity of the cell cycle activity and the latter via the growth and division of cells, driving flows and generating dynamical stresses. Our motivation is to develop a minimal model of this feedback, noting that the cell cycle model can be made more sophisticated, e.g. by alternative choices of f or g in eqs.
(1) and (2); the incorporation of apoptosis, cell motility or the supply and removal of nutrients, oxygen and waste byproducts; extensions to 3D or heterogeneity, e.g. mimicking different tissue or cell types. We developed a particle-based simulation to study this model, complemented with a simple continuum theory that is found to be in broad agreement with the simulations. In this work we have focussed mainly on colony development in a quasi 2D channel with adjustable substrate friction, neglecting the role of active motility.
This simple model may be useful as a reference tool to characterise substrate friction using pre-calibrated cell lines. We also studied the spatial structure of the growing front, where cells are proliferating and moving outwards. We focus on four fundamental variables: outward velocity, local pressure, cell cycle activity and cell area. Cells switch from proliferation to a quiescent state as they transition from the front into the interior (bulk). We propose a method to relate the parameters underlying the cell cycle model to experimental observables. We also compare our simulations with experiments on expanding colonies of MDCK cells, noting that these cells are motile, a feature absent in this version of our model. We find broad agreement for the area distribution and division rates, supporting the use of sizing and division mechanisms under the control of the cell-cycle.
In summary, we hope this work has helped to establish that the cell cycle can play a non-trivial role in the physics of dividing cells and that this might inform future model development. We hope the current paper motivates experimental work on growing cell colonies in 2D to isolate and further evaluate the role of the cell cycle in collective cell dynamics.
Cells are modelled as soft spheres (disks in 2D) with the interaction between cells described by an elastic repulsive force and constant adhesive force per unit contact area. The repulsive force experienced by cell i due to interactions with cell j is assumed to follow Hertzian contact mechanics [93] according to This force act along the unit vector n ij pointing from the center of cell i to the center of cell j. All variables are dimensionless (see Table I) and all lengths are measured in the quiescent cell size, √ V q in 2D. The overlap of two cells h ij is defined in terms of their center-to-center distance r ij , see Fig 8. The force per area f rep for solid elastic Hertzian spheres can be identified with f rep = (4/3)/ (1 − ν 2 i )/E i + (1 − ν 2 j )/E j with R i , E i and ν i the radius, elastic modulus and Poisson ratio of the i th cell, respectively. Alternatively f rep can be treated merely as a single (here, constant for all i, j) adjustable prefactor, controlling repulsion. The dimensionless overlap between cells i and j is This vanishes for cells that don't overlap: such cells naturally have no interactions. The intercellular adhesion force between cells i and j follows from the approximation that receptor-ligand interactions scale with the dimensionless contact area with f ad a constant that can be related to a more microscopic model for receptor-ligand binding, if desired. In our 2D simulations this contact becomes a dimensionless length, defined as This gives a total force between cells i and j of F ij = F rep ij + F ad ij . Summing all interactions over j gives the force on the i th cell Smoothly varying interactions of the kind chosen here are numerically convenient but the precise form of the interactions is unlikely to be important provided there is some weak attraction and strong enough repulsion to suppress excessively large cell indentations. τ r = τ v = 5τ div Figure 9. Time series of the instantaneous estimate of the dimensionless half-decay length for τr=τv=5τ div . The mean values are stable after ∼ 20τ div and all results in the main text are reported in such a steady-state regime using improved averaging over all data in that regime.

Pressure
During the simulation the instantaneous dimensionless pressure on each cell, measured in units of p r , is calculated by summing the scalar (inward-pointing) force and dividing by the 2D equivalent of the cell's surface area -its circumference, noting that the pressure can be negative if a cell is mainly experiencing attraction towards its neighbours.

Half-decay lengths
In order to estimate the half decay lengths described in the main text we process the simulation results using spatial binning along the channel direction and then use linear interpolation between the binned values. Before sampling, we ensure these have reached steady-state values, see Fig. 9.

Equations of motion
The over-damped equation of motion for the position r i of cell i is where the friction constant ζ is a friction per quiescent cell area, V q in 2D, and so the factor 1 πR 2 i scales this to the current contact area of the i th cell. Eq A7 is integrated forward in time for each particle using a time step ∆t chosen to be small enough for good numerical stability and computational accuracy, see Table I.

Model parameters and units
The standard values of all parameters, used in all simulations unless specified otherwise, are given in Table I. The boundary value problem given by eqs. (14) to (17) and boundary conditions eq. (18) for z > 0 are solved numerically as follows: At the outset, the value of the front speed s is unknown because it self-consistently depends on the full solution. However, we assert that the speed of the cells far away from the front must eventually decay to zero, u(z) → 0 for z large enough. The standard approach for such a system is the shooting method: By varying s systematically, a front speed can be found for which the cell speed far away from the front correctly decays to 0. For any given s, the system becomes an initial value problem. This problem can then be solved by numerically integrating from the growing front, z = 0, to a distance z max that is large enough for all quantities to have approximately decayed to their bulk values. This distance is typically of the order of many hundreds to thousands of cell diameters. We confirm afterwards that z max is much larger than all the decay lengths of the system. Since s and u(z max ) are scalar, the front speed for which the bulk speed decays to zero can be found with a scalar root finder.

Asymptotic expansion around the bulk state
We study how the variables asymptotically approach the constant values characterising the tissue deep in the bulk, u(z) → 0, p(z) → p bulk , r(z) → 0, and v(z) → 1, for z → ∞, with p bulk some unknown pressure, see insets to Fig 4. Expanding about the bulk state as u(z) = δu(z), p(z) = p bulk + δp(z), r(z) = δr(z), v(z) = 1 + δv(z) the equations eqs. (14) to (17) reduce to The system of equations for the perturbations can be written as a matrix equation with two amplitudes A 1 and A 2 which have to be matched to the data. Due to the repeated eigenvalue 0, the solution would normally also contain a constant term and a term linear in z, but we can exclude that because of the condition that the solution must tend to zero in the bulk, z → ∞. We observe that exponential decay with exponent e 2 = −(p bulk − 1)/(sτ r ) tends to describe the data well, see insets of fig. 4. For all our data, we find that e 1 > e 2 , which implies that the asymptotic solution always becomes dominated by the latter exponent eventually. In practice, the bulk pres-sure has to be extracted from the numerical solution; we use the pressure for the largest simulated z for each data set. The fact that a single eigenvalue dominates deep in the bulk means that only a single parameter value, or combination thereof, could be inferred by fitting to corresponding data. This supports the focus on the neighbourhood of the growing front adopted in the main text: The front region provides for better model discrimination and parameter inference.