Minimal cellular automaton model with heterogeneous cell sizes predicts epithelial colony growth

,

Understanding the principles of dynamic tissue organization requires an interdisciplinary effort where quantitative biological observations must be challenged by mathematical models that integrate existent insight on a theoretical level [1,2].Cell-based mathematical models essentially consider a tissue as 'a society of cells' [3] where cells behave and interact with each other and the extracellular environment according to individual cellular characteristics summarizing the effects of sub-cellular processes, such giving rise to complex emergent processes at tissue scale.Existing cell-based model frameworks can be classified into three broad categories, (i) continuum descriptions of tissue movement which rarely describe individual cells explicitly [4][5][6][7][8][9], (ii) models that include an explicit description of the cell membrane, such as the cellular Potts model [10][11][12][13], vertex models [14][15][16][17][18][19] or disk models [18,20], and finally (iii) particle basedmodels or cellular automata where cells are represented as point-like particles [12,[21][22][23][24][25][26][27][28].Each class has its own merits and drawbacks: (i) Tissue scale descriptions allow to capture biophysical principles which govern large scale phenomena like tissue fluidity and deformations (elastic, plastic, viscoelastic) or jamming transitions [29,30].However, relating biophysical parameters of cell colonies to measurable single cell properties remains challenging and simulation of large collectives with single-cell resolution can cause high compu-tational cost [4,5].(ii) Models that describe cells as spatially extended objects, such as cellular Potts and vertex models, can incorporate heterogeneous cell sizes and shapes within the tissue and provide a related contact network for specifying intercellular interaction.However, the cellular behavior and interaction is often described indirectly via membrane adjustments controlled by free energy terms, the latter also involving purely technical parameters which have an impact on the emergent dynamics [12,31,32].This fact together with the considerable computational effort for full exploration of the parameter space limits the calibration to concrete experimental data.(iii) Cellular automata and interacting particle systems are computationally more efficient while being mathematically easier to analyze [33].They allow to describe cellular behavior and intercellular interaction in a rule-based manner based on clear mechanistic hypotheses and with biologically meaningful parameters [34].A major drawback, however, is the idealization of point-like cells hampering the modeling of tissue dynamics where cell sizes and shapes are supposed to play an essential role [35].
Here, we consider monolayers of growing epithelial tissue viewed as a dynamic 2D arrangement of individual cells, which grow in size, divide and migrate in response to mechanical constraints and, potentially, interaction on cell-cell contact.In-vitro experiments on the development of epithe-lial tissue of Madin-Darby canine kidney (MDCK) cells from Puliafito et al. (2012) [36] show two distinct growth regimes: At early times, in the confluent phase, cell density is low and spatially homogeneous and the colony grows exponentially in time.At longer times, in the postconfluent phase, cell density increases inside the colony and subexponential colony growth driven by more motile cells near the edge of the expanding colony is observed.The inhibition of cellular proliferation and migration in crowded environments as observed in the center of the growing monolayers is called contact inhibition of proliferation response and migration.It was found to be a not only a consequence of cell-cell contact but of mechanical constraints that cause successive cell divisions to reduce cell area.To further differentiate between the impact of biomechanical constraints and that of intercellular interaction on cell-cell contact, cell-based mathematical modeling was exploited.Puliafito et al. (2012) [36] used a 1D vertex model to support their findings qualitatively.However, the model parameters and the observables of the emergent behavior were not systemically related to experimental measurements.Aland et al. (2015) [4] could reproduce the observed cell culture topology and radial distribution function with a hybrid continuum-discrete model, however, the spatial scale of the culture had to be extrapolated to realistic culture sizes due to computational cost of the model.Jain et al. (2022) [37] incorporated proliferation in the energetic description of a multiphase field model to predict the coordination numbers of cells in growing epithelial colonies corresponding to different hypotheses on contact inhibition and to confirm the typical, experimentally observed linear boundary growth of the colony radius.The role of cell polarity and cellular cycling between motile and dividing states for contact inhibition of proliferation was explored using a disk model [20].However, the question whether mechanical constraints alone, without further intercellular interaction, can explain the observed tissue-scale dynamics remains still largely open.
We study this question by proposing a minimal model, an extension of probabilistic cellular automata (CA) that allows describing heterogeneous cell sizes and growth of individual cells.In our model, cells located on the nodes of a 2D regular grid can grow, divide and migrate, see Fig. 1 for an illustration.They only interact by the exclusion principle, which means that a cellular action is prevented if there is insufficient space.We calibrate the model parameters, which quantify single-cell dynamics, to previous in-vitro observations based on single-cell measurements from Puliafito et al. (2012) [36], and subsequently compare the emergent features of the calibrated CA model, such as culture size and cell density, to independent experiments from the same study.These features of CA model reproduce quantitatively the corresponding experimental results, implying that the exclusion principle, the only intercellular interaction incorporated in the model, paired with size dependent proliferation rates is sufficient to generate the observed contact inhibition.We thus demonstrate the effectiveness of the CA model for where each square of the grid can be occupied by several cells as long as their total size does not exceed the size of the square i Ai ≤ Amax.Each cell i can proliferate with rate γ(Ai), migrate with rate µ, and grow logistically with rate α to a maximal size of Amax (each action exemplarily sketched for the green highlighted cell in the central square).Cells interact via exclusion principle, i.e., migration and growth of individual cells is inhibited by lack of sufficient free space.Dynamics are numerically integrated using Markov-chain-Monte-Carlo method with time increment ∆t.Note that the position of a cell is not resolved beyond its grid square.
efficient quantitative predictions of tissue organization on the basis of rule-based cellular behavior with heterogeneous cell sizes.We discuss the potential to incorporate effective biomechanical interactions between cells in our CA model for future studies of contact inhibition.

I. DATA-DRIVEN APPROACH
Our data-driven approach is based on the mentioned experiments on the development of epithelial tissue of MDCK cells [36].In this study, the transition from a confluent phase to a postconfluent phase is measured during the expansion of the epithelial colony in terms of colony area, cell density, average cell cycle time τ2 , and average velocity of cell motion over time.Furthermore, it is observed that the cell behavior in the center of the expanding colony at the end of the confluent phase is similar to confluent, but homogeneously seeded cell cultures, except that the transition to a stationary, postconfluent state occurs soon after initiation of the culture.Thus, homogeneously seeded cultures are exploited to further quantify the long-term behavior at and after the transition from the confluent phase in terms of the asymmetry G(f ) of daughter cell sizes after division, distribution of cell sizes ρ(A), cell size-dependent division time τ 2 (A), and At the same time, a radial density gradient is visible, implying that cells at the boundary of the culture are able to grow larger and drive further outwards expansion.Model parameters calibrated as written in the main text.Note that the local average cell size is approximated by the displayed cell density as each pixel represents a lattice square of size Amax = 990 µm 2 , e.g., a grid square with density of 10 1 • 10 3 cells/mm 2 hosts 9.9 cells with average size ≤ 1/density = 100 µm 2 .
cell median size over time.
In the following, we develop an according CA model (II.CA model ), whose parameters are calibrated based on the measured asymmetry G(f ), distribution ρ(A), average cell cycle time τ2 , cell size-dependent proliferation rate γ(A) = ln(2)/τ 2 (A), and time point of the transition t T (III.Parameter calibration).The emergent dynamics of the calibrated CA model are then successfully compared to the experiment in terms of colony area, cell density, and average velocity of cell motion over time for the setting of an expanding colony as well as cell median size over time for the setting of homogeneously seeded cells (IV.Results).Thus, the parameter calibration (based on single-cell dynamics) and subsequent model assessment (based on emergent features of the tissue) rely on separate datasets from independent experiments.

II. CA MODEL
We present an efficient, agent-based dynamical model to simulate the development of two-dimensional tissue.The model is based on a probabilistic cellular automaton (CA) but additionally allows for size changes of individual cells during cell growth and division, see Fig. 1 for an illustration of the model and Fig. 2 for an exemplary simulation of epithelial tissue.In the CA, the two-dimensional space is discretized into a regular grid of M 2 squares each with size A max , which can be occupied by cells.Several cells i with size A i can occupy the same square as long as their cumulative size does not exceed its size i A i ≤ A max .While the position of a cell is not resolved beyond its grid square, the restriction of the total size ensures a monolayer without overlap of cells as in the experiment.Cell dynamics are defined by actions with corresponding rates, i.e., a cell proliferates with rate γ, migrates from one grid square to an adjacent one with rate µ, and grows with rate α.Further interactions could be incorporated, like cell-cell adhesion [23,27] or cellsubstrate adhesion, but are neglected to focus on a minimal dynamical model.
The proliferation rate γ(A) and growth rate α(A) are both assumed to depend on the size of the cell A as suggested by single-cell measurements [36], while the migration rate µ is held constant.Cells interact with each other solely via exclusion principle, that is a cell can only grow or migrate if sufficient free space is available in its current grid square or the adjacent ones.This means that the values α(A) and µ represent the maximal possible rates for growth and migration of a single cell of size A, whereas typically these processes are mitigated or inhibited by neighboring cells competing for the same free space.
The proposed model is a continuous-time cellular automaton, where updates are asynchronous with time to next reaction being an exponential random variable and its dynamics are as usually numerically approximated by a Markovchain-Monte-Carlo simulation with Monte-Carlo-time-step ∆t.Within a time step ∆t, if there are N cells present, 3N -times a random cell and a random action -proliferation, migration or growth -with corresponding rate ω are chosen.The cell then performs the chosen action with probability ω∆t.This procedure emulates numerically that on average each cell performs each action with corresponding rate, i.e., the probability that the action has not been performed decreases exponentially ∼ exp(−ωt), given a sufficiently small time step ∆t ≪ max{γ −1 , µ −1 }.

A. Proliferation
Proliferation is defined as a potentially asymmetric division of a mother cell into two daughter cells keeping the overall area preserved.Measurements of single-cell division times [36,Fig. 4E] suggest a sigmoidal function for the dependence of the proliferation rate γ(A) on the size A of the mother cell This rate increases with cell size and becomes constant at large γ(A ≫ A 0 ) ≈ γ max and small sizes γ(A ≪ A 0 ) ≈ 0, while the steepness of the transition at A 0 is set by the exponent m.Further, measurements of the cell areas during single-cell division [36, Inset of Fig. 4C] suggest an asymmetry between the areas of the two daughter cells, which is taken into account: If a mother cell of size A divides, one of the daughter cells is assigned a fraction f of the mother cell and while the other daughter is assigned the remaining part, meaning the daughter cells have size f A and (1 − f )A.Naturally, the probability density G(f ) of the fraction f ∈ [0, 1] is symmetric around 1/2.Note that the total cell area is not always conserved in the experiment [36,Fig. 4c, lower inset], which is ignored in the model for simplicity.

B. Migration
For migration cells move with rate µ to a randomly selected neighboring grid square.For grid square size A max this implies a maximal migration velocity v max = µ √ A max across the grid or Cell migration is only performed if the target grid square has sufficient free space to accommodate the migrating cell.

C. Growth
We assume a generic, logistic cell growth with maximal cell size A max and growth rate α.Note that the growth rate α indirectly affects cell proliferation, as it affects the position and width of the cell size distribution, which in turn determines the average proliferation rate according to Eq. ( 1).For instance, increasing the growth rate α leads to accelerated growth of the cell culture.While the logistic growth Eq.( 3) is time-continuous and could be integrated directly via forward-Euler-scheme or the corresponding analytical solution, we employ for simplicity the same Markov-chain-Monte-Carlo scheme as for proliferation and migration.We use α = 1/∆t as the rate of the growth action, meaning it is performed on average once per time step and cell, and define the action as increasing the cell size by ∆A = αA(1 − A/A max )∆t.Note that the time step is chosen sufficiently small such that ∆A/A ≪ 1.Furthermore, the increment ∆A is limited to the available free space in the grid square.It should be pointed out that in the previously proposed 1D vertex model, cell growth was assumed to be a consequence of the surrounding tissue stretching the cells [36].In contrast, we explicitly choose a growth function that does neither rely on cells actively pushing or pulling on each other, to devise a model with minimal cell-cell interactions.

III. PARAMETER CALIBRATION
The biologically motivated parameters γ(A), v max , α of the CA model are calibrated based on in-vitro observations, mostly of single-cell dynamics, during the growth of normally differentiating epithelial cell cultures of Madin-Darby canine kidney cells [36]: (III A) The division asymmetry G(f ) is directly fitted to data from tracking of single cell divisions.Then, this asymmetry G(f ) is incorporated into a partial differential equation model to optimize the proliferation rate γ(A) based on cell size distributions measured in cell tracking experiments.The obtained rate γ(A) is validated against independent data of cell division times in a stationary, postconfluent state.(III B) The maximal migration velocity v max , and consequently the migration rate µ, is estimated from the time point of the transition to a stationary, postconfluent state and the average cell cycle time.(III C) Finally, the growth rate α is optimized using simple Monte-Carlo-Simulations of independent cells based on the average cell cycle time using the already calibrated parameters γ(A), G(f ).Note that the two mathematical models supporting the calibration are much simpler than the CA model, without spatial resolution or inter-cellular interaction.

A. Proliferation
Firstly, the probability density G(f ) of the division asymmetry can be directly fitted to corresponding measurements from tracking single-cell divisions, see Fig. 3.The data is well described by a Gaussian distribution with expectation value f = 1 2 , standard deviation σ = 1 8 , boundaries of the fraction f ∈ [ 1 4 , 3  4 ], and corresponding normalization  5), dashed lines), which only incorporates cell size-dependent proliferation γ(A), Eq. ( 1), and asymmetric cell division G(f ), Eq. ( 4).The distribution at day 0 serves as initial condition for the PDE.Deviations from the experimental distributions originate mainly from overestimated tails at small cell sizes, which could hint to modification of the asymmetry of the cell division at cell sizes < 100 µm 2 , see main text.
The parameters γ max , A 0 , m defining the proliferation rate γ(A) are optimized based on the observed dynamics of the cell size distribution ρ(A, t), see solid lines in Fig. 4.These distributions are measured in a regime of the cell culture growth, where cells no longer migrate or increase their size due to their dense packing.Instead, these cells only divide into smaller and smaller cells whose proliferation rates decrease with their size.We model this loss of larger cells and rise of smaller cells in the distribution ρ(A) using a partial integro-differential equation (PDE) (5) The first term reflects the loss of cells with size A due to division with proliferation rate γ(A).The second and third term reflect the gain of cells with size A due to division of larger cells with sizes A/f or A/(1 − f ), where one of the resulting daughter cells obtains a fraction f or 1 − f , respectively, of the mother's size.Note that Eq. ( 5) simplifies to the equation proposed in Ref. [36], if a delta distribution G(f ) = δ(f − 1/2) is assumed.We solve Eq. ( 5) numerically by discretizing the cell size A → A i∈N and integrating the resulting set of ordinary equations with odeint [38].We optimize the parameters γ max , A 0 , m of the proliferation rate γ(A) by minimizing the distance D between the distributions ρ exp (A i ), ρ PDE (A i ) obtained from experiment and PDE, respectively, at selected times t using the python package scipy.optimize[38].We select experimental distributions at three characteristic time points from the experiment, where the distribution at the first time point ρ exp (A i , t ≈ 0) serves as initial condition for the PDE while the other two time points t ≈ 0.5 d and t ≈ 6 d are the reference points for the optimization in Eq. ( 6).While distributions are reported at more intermediate time points in Ref. [36], they do not significantly contribute to the optimization as these distributions transition smoothly into each other and thus can be neglected for clarity.
We obtain parameters γ max ≈ 1.3 d −1 , A 0 ≈ 80 µm 2 , m ≈ 3, for which we observe good agreement both with the dynamics of the cell size distributions, see Fig. 4, and the division times of single cells, see Fig. 5. Deviations from the experimental distributions in Fig. 4 originate mainly from overestimated tails at small cell sizes, left of the maximum of each distribution.These deviations seem to be diminished, when the domain of the size fraction p, Eq. ( 4), is further narrowed around the symmetric case p = 1/2.Thus, the deviations in Fig. 4 could hint at a narrowing of the asymmetry of the cell division at small at cell sizes < 100 µm 2 , e.g., due to a minimal absolute size a cell must have.In fact, the distribution G(p), which the PDE relies on, has been obtained for mother cell sizes A ∈ [100, 1000] µm 2 , due to the difficulty of tracking cells below 70 µm 2 , while the cell size distributions ρ(A) were mostly measured at A < 100µm 2 [36].Furthermore, the PDE assumes conservation of total cell area for simplicity, which is not always the case in the experiment [36,

B. Migration
According to Eq. ( 2), the migration rate µ is directly defined by the maximal migration velocity v max .This velocity can be estimated from the transition point t T from exponential to subexponential growth of the cell culture (transition to a stationary, postconfluent state [36]): As cells proliferate and grow in size, they inevitably migrate outwards causing the cell culture to expand.As long as this expansion is fast enough, each cell has sufficient space to grow and divide at maximal capacity and the cell culture grows exponentially, while the density in the culture remains constant, see black points in Figs. 6 and 7 before t < 5 d, respectively.However, at some point the outward migration of the cells cannot compensate for the increase in total cell size and consequentially the growth of invidiual cells in the center of the culture is inhibited by the lack of free space.This leads to subexponential growth of the culture and an increase in cell density, see black points in Figs. 6 and 7 after t > 6 d, respectively.For an exponential growth of a circular cell culture with the average proliferation rate γ = ln(2)/τ 2 , given in terms of the reported average cell cycle time τ2 = 0.75 ± 0.14 d [36].Note that in the context of the model, the rate γ should reflect the average proliferation rate over an ensemble of independently proliferating, Eq. ( 1), and growing, Eq. (3), cells, which is latter exploited for the calibration of (III C) cell growth.Along with the initial culture size A cult (0) ≈ 1.8 • 10 4 µm 2 , the maximal migration velocity v max can then be estimated from the condition The range of the cell cycle time τ2 and the transition point t T = 6 ± 1 d imply a range for the maximal migration velocity v max ∈ [9, 100] µm/h.In Ref. [36] a coarse consistency check between the observed average cell velocity, transition time t T , average cell cycle time τ2 , and colony area A cult (t T ) was performed based on the same assumption as Eq. ( 8).However, the ranges for t T and τ2 were neglected and the maximal and average cell velocity were assumed to be the same, while the former is a single-cell parameter and the latter an emergent feature from the interaction with other cells.Note that the condition in Eq. ( 8) is a simplification and rather estimates a lower limit for v max .The transition from exponential to subexponential growth is expected before the rim of the cell culture actually reaches an outward velocity v max .For instance, at t = 5 d the culture already exhibits deviations from exponential growth, while according to Eq. ( 8) the velocity of the boundary at this point is only 8 − 30 µm.In order to reflect the estimated range for v max , we use its upper limit v max = 100 µm/h in the following, but also report the results for v max = 50 µm/h in SI Figs.4-7.

C. Cell growth
The average proliferation rate γ derived from the measured average cell cycle time is also used to estimate the cell growth rate α: Firstly, the maximal cell size is set to A max = 990 µm 2 based on tracking of sizes of single cells [36, Fig. 4A] and the observed minimal cell density.For the phase of exponential growth of the culture, it is reasonable to assume that cells proliferate and grow independently of each other according to Eqs. (1), (3), (4).We compute these simple dynamics numerically for ensembles of cells using Monte-Carlo simulations.In particular, we use the parameters of proliferation γ(A), G(p) as calibrated above and additionally optimize the growth rate α of logistic cell growth Eq.( 3) such that the resulting exponential slope of the total size of the cell ensemble matches with the slope γ = ln(2)/τ 2 corresponding to the average cell cycle time τ2 = 0.75 ± 0.14 d measured independently by single-cell tracking.From this optimization a growth rate α = 1.3 ± 0.3 d −1 is obtained.

D. Initial and boundary conditions
For expanding culture experiments, the size of the spatial grid is chosen sufficiently large, such that cells do not reach the boundary of the grid during simulation.We obtain the initial condition of the CA model by starting from a single  1B]).The colony size exhibits initially an exponential growth (low-density, confluent regime) and slows around time tT = 6 ± 1 d to a subexponential growth (transition to a stationary, postconfluent state).Results from CA model are averaged over 15 independent runs and standard deviation is displayed as red shaded region.Model parameters calibrated as written in the main text.
cell of maximal size in the center of the grid and then simulating the dynamics until the size of the culture matches the initial one in the experiment.For experiments where cells are seeded homogeneously, a 12×12-grid with reflecting boundary conditions is used and as initial condition a single cell of maximal size is placed in every second grid point, reflecting the experimental seeding density of 600 cells/mm 2 .Finally, the time increment in the CA is set sufficiently small ∆t = 6 min such that ∆t ≪ max{γ −1 max , µ −1 } and ∆A/A ≪ 1.

IV. RESULTS
We find quantitative agreement between the calibrated CA model and independent experiments with respect to several emerging features of the epithelial cell culture.This validates our minimal model and showcases its power to predict emergent tissue wide dynamics from single cell behavior.The agreement also suggests that the exclusion principle as sole inter-cellular interaction paired with a size-dependent proliferation rate of each cell is sufficient to explain the emerging dynamics of contact inhibition: We implement the previously calibrated parameters into the CA model and compare the emerging dynamics of culture size and cell density with the experiment of epithelial colony growth, see Figs. 6 and 7.In the CA, the culture size is computed as the total area of all grid points containing at least a single cell and the cell density is computed as ratio  1C], corresponding data to areas displayed in Fig. 6).The cell density is initially constant (low-density, confluent regime) and increases rapidly around time tT = 6 ± 1 d (transition to a stationary, postconfluent state) implying dense packing of cells, which become smaller and smaller due to slowing but continued proliferation without cell growth.Note that the small variations (±30%) until day 6 appear larger due to the logarithmic y-scale and were interpreted as constant density by the original experimenter [36].Results from CA model are averaged over 15 independent runs and standard deviation is displayed as red shaded region.Model parameters calibrated as written in the main text.
between this area and the total number of cells.The resulting dynamics of culture size and cell density match with the experimental ones.In particular, the CA model does not only quantitatively reproduce the exponential growth and constant density in the initial low-density confluent regime and the time point of the transition to a stationary, postconfluent state, but also the curve at and after this transition.The predictions of the CA model are also consistent for the case of initially homogeneously seeded cells, see Fig. 8, which reflects the regime directly at and long after the transition to a stationary, postconfluent state.Furthermore, the transition is accompanied by a rapid decrease in the emergent cell velocities, which are comparable between CA model and experiment, see Fig. 9.In particular, there is a velocity gradient from the inside to the outside of the colony and an outwardly biased flow, see Fig. 10 as well as SI Figs. 1 and  2.
Note that the calibration of the model parameters is only based on observations from experiments independent from the one we compare the model predictions to.The only exception is the time point of transition t T used for the calibration of the maximal velocity v max , as the maximal velocity at the boundary of the culture has not been directly measured.Apart from that, the majority of the calibration is based on tracking of the dynamics of single cells and the model parameters themselves describe only single-cell behavior.
The match of the emergent features, i.e., dynamics of culture size and cell density, between experiment and CA model suggests that they are determined by the single-cell dynamics.In particular, the CA model does not incorporate any inter-cellular interactions, besides the exclusion principle.In the CA model, the transition of the dynamics at t T originates solely from the limitation that the boundary of the culture cannot expand faster outwards than the maximal migration velocity v max of the cells.Beyond the corresponding critical size of the culture, cells in the center of the culture no longer get access to the necessary space to grow, i.e. dA/dt → 0, and consequently become smaller due to continued proliferation.The latter gets slower due to the sizedependent, sigmoidal proliferation rate γ(A).The capacity of cells to grow separates their proliferation behavior into a pre-(A ≫ A 0 → γ(A) = γ max ) and posttransition regime (A ∼ A 0 → γ(A) ≪ γ max ) according to Eq. ( 1).Note that based on the assumption that for times t ≫ t T effectively only cells in a boundary layer with thickness ∆R of the culture can proliferate with rate γ, the asymptotic, posttransition growth may be estimated by dA/dt ≈ 2 , analogous to recent estimates for spherical cultures [39], while ∆R(v max ) may depend on the maximal migration velocity.We obtain an effective thickness of the proliferating rim ∆R = 370 ± 10 µm (t 0 = 3.3 ± 0.2 d) from the CA simulation, consistent with rescaled results from a previous theoretical model of epithelial colonies [4].Note that while the cell size-dependent division times t 2 exhibited by the calibrated CA model closely resemble the experimental ones, see SI Fig. 3, the distribution of these times is broader than in the experiment, see inset in Fig. 3.In fact, the division times in the CA model follow an exponential distribution by design for mathematical and numerical convenience, whereas biological cells divide at more regular intervals leading typically to a narrower Erlang-like distribution.However, our results imply that this difference in the cell cycle distributions does not affect the emergent dynamics of the colony and the implementation of cell division according to regular cell cycles into the CA model is straightforward.
In the CA model, the observed rapid decrease in cell velocity at the transition results from cells in the inner region of the colony slowing down due to insufficient space to migrate, which generates a velocity gradient from inside to outside, see Fig. 10.At the same time cell motion transitions from undirected to a collective outwards bias, see also Fig. 1.This outwards directed motion is an emergent feature of the interplay between the migration rate, the exclusion principle, and the fact that more free space is available outwards.It should be emphasized that for the CA model the maximal migration velocity v max is significantly bigger than the velocities that are measured for individual cells inside the culture, see Fig. 9, due to the exclusion principle.velocity by a factor of 3/4 to 1/4.Moreover, the migration velocity in the CA model represents the speed of instantaneous, undirected movement of single cells, whereas in invitro experiments typically collective, directed movement of groups of cells is measured, e.g., directed outwards in case of an expanding culture.Consequently, this directed motion exhibits much smaller velocities, consistent with previous experimental measurements of cell velocity in epithelial tissue ∼ 25µm/h depending on the setting (wound-cut technique, expanding colony, ...) [36,[40][41][42][43].Note that the exact value of the average velocity depends to some degree on the chosen definition, i.e., the chosen time increment between two cell position measurements and whether the average is taken over the whole colony or a specific location.In Ref. [36] the cell velocity was extracted from four local 450 × 336 µm 2 fields of view, while for the CA model it is possible to average over the entire cell culture.This explains the large range of velocities observed in the CA model, which reflects the radial velocity gradient, see Fig. 10, while the standard error for the experimental velocities is not reported.Nevertheless, we can infer from the match of simulation and experiment in Fig. 9 that the course of the cell velocities, initially constant with a rapid decrease around the transition point, is reproduced by the CA model.

V. DISCUSSION
We developed a mechanistic, data-driven model to simulate the dynamics of planar cell cultures by extending a corresponding data to Figs. 6 and 7).The average cell velocity is initially constant (low-density, confluent regime) and decreases rapidly around time tT = 6 ± 1 d (transition to a stationary, postconfluent state).Presented are the root-mean-square velocities, where the model velocities are computed from a single simulation run over time increments of 580 min to neglect small-scale oscillations.Due to the presence of a radial velocity gradient, see Fig. 10, the corresponding standard deviation of the velocities (red shaded region) is quite large.
probabilistic cellular automaton to incorporate size changes of individual cells during growth and cell division.We successfully apply this CA model to previous in-vitro experiments [36] of epithelial tissue composed of Madin-Darby canine kidney cells.In particular, the model parameters are calibrated using measurements based on single-cell tracking and subsequently the resulting CA model is validated against independent experiments on the development of culture size, cell density, and cell velocity.The agreement between model predictions and these experimental observations suggests that the exclusion principle as sole inter-cellular interaction paired with a size-dependent proliferation rate of each cell is sufficient to explain the emerging dynamics of contact inhibition.In particular, incorporation of additional mechanisms like cell-cell signaling or cell adhesion are not necessary.In the CA model, the reproduced transition from exponential to quadratic growth of the colony originates from the increasing lack of space in the center of the colony, inhibiting cell growth and division, when the expansion of the colony can no longer be compensated for by the migration of the exterior cells.This origin of the transition is consistent with previous interpretations and the observation that contact inhibition can take place several days after cells have been in contact [4,36].Moreover, the collective, outwards motion of the exterior cells and the rapidly increasing cell density in the center of the culture are a consequence of these simple mechanical constraints.
As mentioned there is a wide variety of theoretical models to describe and simulate the development of healthy and cancerous cellular tissue and even contact inhibition [44][45][46][47][48][49][50][51][52][53].Most of these models explicitly consider intra-and intercellular forces, allowing to study the delicate interplay between biomechanical quantities like tension, pressure, adhesion, as well as cellular stiffness and compressibility, which govern cell topology, e.g., the distribution of the number of direct cell neighbors, and large scale phenomena like elastic, plastic, viscoelastic deformations.For example, previously a mechanistic hybrid continuum-discrete model has been applied to the same experiment of epithelial cell colonies as our study [4].Besides culture size, cell density and cell area median, this model was able to reproduce the observed cell topology and radial distribution function.However, the spatial scale of the culture had to be fitted empirically to the culture size, as the model underpredicts the cluster size where the transition to a stationary, postconfluent state occurs, and the resulting cell velocities are reported to be ten times smaller than in the experiment.This is partially due to computational cost, which prevented the use of a larger parameter of mobility, which also reflects effects of cell-substrate adhesion.This example highlights two more general issues of the above-mentioned types of models, see also Ref. [12]: Firstly, they often rely on empirical or effective parameters, like coefficients of free energy functions or elastic properties of tissue, which are either biologically not known or at least very difficult to measure, in contrast to the properties of single cells.Secondly, they often entail high computational effort, which limits exploration of the parameter space and consequently the calibration to concrete experimental data.
While CA models have the advantage, that their parameters are usually based on biologically motivated, single-cell behavior and that they are computationally relatively cheap, the incorporation of biomechanical forces has been a challenge with few exceptions like cell adhesion [23,27].The type of CA we propose here has the potential to partially overcome this challenge: The fact that cells of different sizes can be located at the same grid point naturally defines local cell densities and density gradients.Together with the possibility to define cellular, size-dependent mechanisms and forces, beyond proliferation and growth, this allows to effectively introduce biomechanical processes from other types of models, e.g., homeostatic pressure [7,53] or cell bulk stress in terms of cell compression [4] from the ratio of the actual size of a cell and its (dynamic) target size.Note that the concept of local densities has already been exploited in previous spatial stochastic models, which allow multiple particles in each lattice site [54,55], e.g., to define migration probabilities based on the density differences between lattice sites, which is similar to our CA model if cell size and cell growth is omitted.Mean-field approximations of these models lead to partial differential equations of the reaction-advection-diffusion type exhibiting matching dynamics [55,56].However, these models neglect cell growth, a crucial effect which complicates the derivation of a corresponding continuum model.An example is the already mentioned hybrid continuum-discrete model, derived from a stochastic particle model using the framework of dynamic density functional theory [4].Furthermore, in our CA model cells are capable of moving past each other, which is typically not inherent to both common CA models and vertex models.This feature can be exploited to study effects like collective motion or fingering [57,58] by adding direction vectors of persistent motion to each cell, which are affected by neighboring cell motion or cell densities.This could help to reproduce the finger-like protrusions observed for the expanding epithelial colony [36], which neither our CA model nor the corresponding vertex model [4] exhibit.Moreover, the concept of local cell densities allows to couple the CA model to additional PDEs to incorporate the interaction of cells with external fields like oxygen, nutrients, and metabolic waste.Finally, the CA model can be directly extended to three-dimensional space.However, features of small spatial scales like cell shape [4] and correlation lengths [4,36] can not be captured by our model.While none of these concepts were in the end necessary for the application of the CA model to the considered experiments of epithelial growth, they may guide efficient, data-driven modeling of future scenarios of tissue dynamics.

VI. ACKNOWLEDGEMENT
We are grateful for discussions with S. Rühle, M. Luft and F. Freier.The authors acknowledge that this re-search has been co-financed by the EU, the European Social Fund (ESF), and by tax funds on the basis of the budget passed by the Saxon state parliament (project SAB-Nr.100382145), the Bundesministerium für Bildung und Forschung (BMBF 16dkwn001a/b), and the Sächsisches Staatsministerium für Wissenschaft und Kunst (SMWK) (project FORZUG II TP 3).The publication of this article is funded by the Open Access Publication Fund of Hochschule für Technik und Wirtschaft Dresden -University of Applied Sciences.The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

VII. CONFLICTS OF INTEREST
The authors declare no conflict of interest.

VIII. DATA AVAILABILITY
All relevant data are within the manuscript and its Supporting information files.The code is publicly available under https://github.com/j-schmied/cellularautomata.1B]).The colony size exhibits initially an exponential growth (low-density, confluent regime) and slows around time tT = 6 ± 1 d to a subexponential growth (transition to a stationary, postconfluent state).Results from CA model are averaged over 15 independent runs and standard deviation is displayed as red shaded region.Model parameters as Fig. 6 in the main text, but for a smaller maximal migration velocity vmax = 50 µm/h.1C], corresponding data to areas displayed in Fig. 4).The cell density is initially constant (low-density, confluent regime) and increases rapidly around time tT = 6 ± 1 d (transition to a stationary, postconfluent state) implying dense packing of cells, which become smaller and smaller due to slowing but continued proliferation without cell growth.Results from CA model are averaged over 15 independent runs and standard deviation is displayed as red shaded region.Model parameters as Fig. 7 in the main text, but for a smaller maximal migration velocity vmax = 50 µm/h.

FIG. 1 .
FIG. 1. Illustration of the cellular automaton (CA) model.Cells (red and green) with individual sizes Ai are positioned on a regular grid (gray 3 × 3 squares represent a section of this grid), where each square of the grid can be occupied by several cells as long as their total size does not exceed the size of the square

FIG. 2 .
FIG. 2. Exemplary simulation from the CA model of the growth of epithelial tissue in the confluent regime (1 d, 3.3 d, 5 d), as well as at (6.2 d, 7.3 d) and post-transition (8.3 d, 8.8 d, 9.8 d) to a stationary postconfluent state.For each time point, the local cell density is encoded in color.The growing culture expands outwards over time.Around day 6 the density increases drastically in the center of the culture, as the lack of free space inhibits migration and cell growth, reflecting the transition to a postconfluent state.At the same time, a radial density gradient is visible, implying that cells at the boundary of the culture are able to grow larger and drive further outwards expansion.Model parameters calibrated as written in the main text.Note that the local average cell size is approximated by the displayed cell density as each pixel represents a lattice square of size Amax = 990 µm 2 , e.g., a grid square with density of 10 1 • 10 3 cells/mm 2 hosts 9.9 cells with average size ≤ 1/density = 100 µm 2 .

FIG. 3 .
FIG. 3. Asymmetry G(f ) of the daughter cell sizes after cell division is well described by a Gaussian distribution with expectation value f = 1/2, standard deviation σ = 1/8, and a cut-off at large and small cell sizes f ∈ [1/4, 3/4], reflecting the absence of smaller or larger daughter cells (black dots represent data from Ref. [36, Fig. 4C (inset)] obtained from tracking single-cell divisions, Gaussian distribution shown as red line).It is assumed that during division of a cell, one daughter cell has a fraction f of the size of the mother cell, while the other cell has a fraction (1 − f ).

FIG. 4 .
FIG.4.Distribution of cell sizes over time after transition to a stationary, postconfluent state (solid lines selected from experimental data [36, Fig.4D]) is well described by partial differential equation (PDE Eq. (5), dashed lines), which only incorporates cell size-dependent proliferation γ(A), Eq. (1), and asymmetric cell division G(f ), Eq. (4).The distribution at day 0 serves as initial condition for the PDE.Deviations from the experimental distributions originate mainly from overestimated tails at small cell sizes, which could hint to modification of the asymmetry of the cell division at cell sizes < 100 µm 2 , see main text.

FIG. 6 .
FIG.6.Emergent dynamics of calibrated CA model (red line) reproduces experimentally observed size of expanding epithelial colony (black dots from Ref.[36, Fig. 1B]).The colony size exhibits initially an exponential growth (low-density, confluent regime) and slows around time tT = 6 ± 1 d to a subexponential growth (transition to a stationary, postconfluent state).Results from CA model are averaged over 15 independent runs and standard deviation is displayed as red shaded region.Model parameters calibrated as written in the main text.

FIG. 7 .
FIG.7.Emergent dynamics of calibrated CA model (red line) reproduces experimentally observed cell densities of expanding epithelial colony (black dots from Ref.[36, Fig. 1C], corresponding data to areas displayed in Fig.6).The cell density is initially constant (low-density, confluent regime) and increases rapidly around time tT = 6 ± 1 d (transition to a stationary, postconfluent state) implying dense packing of cells, which become smaller and smaller due to slowing but continued proliferation without cell growth.Note that the small variations (±30%) until day 6 appear larger due to the logarithmic y-scale and were interpreted as constant density by the original experimenter[36].Results from CA model are averaged over 15 independent runs and standard deviation is displayed as red shaded region.Model parameters calibrated as written in the main text.

FIG. 8 .
FIG. 8. Emergent dynamics of calibrated CA model (red line) consistent with experimentally observed cell size median for initially homogeneously seeded cells (black dots from Ref. [36, Fig. 3C], transition tT to a stationary, postconfluent state corresponds roughly to 1.5 d).Cell size median initially drops rapidly and then saturates at small cell sizes A ∼ 40 µm 2 .Results from CA model are averaged over 15 independent runs and standard deviation is displayed as red shaded region.

FIG. 9 .
FIG. 9. Emergent dynamics of calibrated CA model (red line) reproduces qualitatively the experimentally observed cell velocities of expanding epithelial colony (black dots from Ref.[36, Fig. 2C], corresponding data to Figs.6 and 7).The average cell velocity is initially constant (low-density, confluent regime) and decreases rapidly around time tT = 6 ± 1 d (transition to a stationary, postconfluent state).Presented are the root-mean-square velocities, where the model velocities are computed from a single simulation run over time increments of 580 min to neglect small-scale oscillations.Due to the presence of a radial velocity gradient, see Fig.10, the corresponding standard deviation of the velocities (red shaded region) is quite large.

v
FIG. 10.CA model displays an outwards biased cell motion (right panel) accompanied by a velocity gradient from inner region of colony to its periphery (left and right panel encoded by color).Displayed are cell velocities v of the last 580 min of the simulation (posttransition) spatially averaged (left panel: 4 × 4 grid squares, right panel: higher resolution 2 × 2 grid squares) where arrows represent direction of motion (left panel) and color magnitude of velocity (linear scale in left panel and logarithmic scale in right panel for visibility, missing arrows or white squares mean no measurable cell motion).For a radial profile of the velocities see SI Fig. 2.

FIG. 1 .FIG. 2 .
FIG. 1. CA model displays transition from undirected, uniform, fast cell motion to outwards biased motion with radial velocity gradient (from top to bottom, note shrinking scale bar).Panels analogous to Fig. 10 in the main text, but for different times around the transition point tT = 6 ± 1 d.For a radial profile of the velocities see Fig. 2.

FIG. 3 .
FIG.3.Emergent dynamics of the calibrated CA model (red triangles) reproduce the observed division times τ2 for an initially homogeneously seeded cells (black dots from Ref.[36, Fig.  4E],blue line corresponds to the fitted proliferation rate from Fig.5).The division times of the CA model are averaged over intervals of the cell size A and corresponding standard deviation, resulting from the underlying exponential distribution of the division times, is shown as red shaded region.In the inset the according distributions of the division times of CA model (red) and experiment (black [36, Fig.4Einset]) are compared.While division times in the CA model follow an exponential distribution, biological cells divide at more regular intervals leading typically to a narrower distribution.

FIG. 4 .
FIG. 4. Emergent dynamics of calibrated CA model (red line) reproduces experimentally observed size of expanding epithelial colony (black dots from Ref. [36, Fig.1B]).The colony size exhibits initially an exponential growth (low-density, confluent regime) and slows around time tT = 6 ± 1 d to a subexponential growth (transition to a stationary, postconfluent state).Results from CA model are averaged over 15 independent runs and standard deviation is displayed as red shaded region.Model parameters as Fig.6in the main text, but for a smaller maximal migration velocity vmax = 50 µm/h.

FIG. 5 .
FIG.5.Emergent dynamics of calibrated CA model (red line) reproduces experimentally observed cell densities of expanding epithelial colony (black dots from Ref.[36, Fig. 1C], corresponding data to areas displayed in Fig.4).The cell density is initially constant (low-density, confluent regime) and increases rapidly around time tT = 6 ± 1 d (transition to a stationary, postconfluent state) implying dense packing of cells, which become smaller and smaller due to slowing but continued proliferation without cell growth.Results from CA model are averaged over 15 independent runs and standard deviation is displayed as red shaded region.Model parameters as Fig.7in the main text, but for a smaller maximal migration velocity vmax = 50 µm/h.

FIG. 6 .FIG. 7 .
FIG.6.Emergent dynamics of calibrated CA model (red line) reproduces qualitatively the experimentally observed cell velocities of expanding epithelial colony (black dots from Ref.[36, Fig. 2C], corresponding data to Figs.4 and 5).The average cell velocity is initially constant (low-density, confluent regime) and decreases rapidly around time tT = 6 ± 1 d (transition to a stationary, postconfluent state).Presented are the root-mean-square velocities, where for the model velocities are evaluated from a single simulation run over time increments of 580 min to neglect small-scale oscillations and standard deviation is displayed as red shaded region.Model parameters as Fig.9in the main text, but for a smaller maximal migration velocity vmax = 50 µm/h.