Mechanical control of cell proliferation patterns in growing tissues

Cell proliferation plays a crucial role in regulating tissue homeostasis and development. However, our understanding of how cell proliferation is controlled in densely packed tissues is limited. Here we develop a computational framework to predict the patterns of cell proliferation in growing tissues, connecting single-cell behaviors and cell-cell interactions to tissue-level growth. Our model incorporates probabilistic rules governing cell growth, division, and elimination, while also taking into account their feedback with tissue mechanics. In particular, cell growth is suppressed and apoptosis is enhanced in regions of high cell density. With these rules and model parameters calibrated using experimental data, we predict how tissue confinement influences cell size and proliferation dynamics, and how single-cell physical properties influence the spatiotemporal patterns of tissue growth. Our findings indicate that mechanical feedback between tissue confinement and cell growth leads to enhanced cell proliferation at tissue boundaries, whereas cell growth in the bulk is arrested. By tuning cellular elasticity and contact inhibition of proliferation we can regulate the emergent patterns of cell proliferation, ranging from uniform growth at low contact inhibition to localized growth at higher contact inhibition. Furthermore, mechanical state of the tissue governs the dynamics of tissue growth, with cellular parameters affecting tissue pressure playing a significant role in determining the overall growth rate. Our computational study thus underscores the impact of cell mechanical properties on the spatiotemporal patterns of cell proliferation in growing tissues.

A common choice of the Hamiltonian [1] includes terms that represent contact energies and an area elastic energy, as in The contact energy, J i, j , is given by the difference between the cellular surface tension and intercellular adhesion energy. This energy is represented in the Hamiltonian by a component that is proportional to the number of pixels shared between cells or between a cell and the medium. The sum goes over all adjacent lattice sites i and j, ignoring self-adhesion. The inner term, τ, gives the cell type at the lattice site i and j respectively. A value for contact energy is defined between each cell type, which for our simulations consists of interactions between cells and medium J cell-cell = J cc , J cell-medium = J cm , and J medium-medium = J mm . The second term in the Hamiltonian is the area elastic energy. This is given by the product of the area elastic modulus, λ , with the sum of squared differences between the actual area of the cells, A α , and their target areas, A T α . These two mechanical energy components regulate individual cell shape and size, in the absence of other rules. The CPM simulations were run using CompuCell3D [4], a software package that handles all of the computations necessary for the CPM and the Metropolis-Hastings algorithm.

Rules for active cell behaviors
The decision-making layer of the model encodes probabilistic rules for active cell behaviors including changes in cell target area, cell cycle regulation, and cell elimination.
Dynamics of cell growth -Cell growth dynamics are implemented by increasing the cell's target area as where G is the growth rate in the absence of crowding, k quantifies cellular sensitivity to crowding, and A α is the current area of cell α. In the presence of crowding, a cell's area will deviate from its target area, and this deviation is directly proportional to the strength of contact inhibition. Consequently, cellular growth exhibits an exponential decline as crowding intensifies, ultimately leading the tissue to enter an arrested state.
Cell cycle regulation -We implemented a two-phase cell cycle model, comprising a G1 growth phase operating through a sizer mechanism, followed by an S/G2/M timer phase. When a cell is born, it is assigned a G1 area threshold A S and an S/G2/M phase timer τ. The cell remains in the G1 phase until its area exceeds A S . Afterward, it enters the S/G2/M phase, which lasts for a duration of τ. Upon completion of the S/G2/M phase, the cell immediately undergoes division, perpendicular to the semi-major axis, resulting in two daughter cells, both of which enter the G1 phase.
Mechanisms for cell elimination -TIn the simulation, cell elimination involves two mechanisms: programmed cell death (apoptosis) and live-cell extrusion. Experimental data have shown that increasing cell density leads to a higher apoptosis rate [5][6][7]. We defined the local cellular density as the sum of inverse cell areas within a local neighborhood of a cell. The density of cell α is given by: To model the probability of apoptosis in wild-type MDCK cells as a function of local cell density [7], we fit the data to a sigmoid curve: In the simulations, the probability of apoptosis is calculated every 10 Monte Carlo steps, which, on the time scale of the timer threshold, reproduces the experimentally observed probabilities. When a cell is eliminated through apoptosis, its target area is set to zero, causing it to shrink rapidly. During this process, surrounding cells grow to occupy the space left by the apoptotic cell.
In addition to apoptosis, cells can also be eliminated via live cell extrusions. This occurs when a cell's area significantly decreases compared to the average population area. In experiments [5], extrusion happens when tissue experiences compressive stress, causing a cell to delaminate from the substrate and eject from the monolayer. In our simulations, a cell whose area falls below A /4 (where A is the average area of the colony) is immediately removed to capture the short time scale over which extrusion occurs [5,8].

Choice of Default Model Parameters
Converting lattice units in the CPM framework to physical units is essential for comparing our simulations with experimental observations. There are two key conversions: translating pixels and Monte Carlo Steps (mcs) into SI units for distance and time. To achieve this, we compared the experimentally reported values for average cell volume during subconfluence for MDCK cells [9][10][11]  The single-cell growth rate G is sampled from a Gaussian distribution to determine the change in the target area at each Monte Carlo step (Eq. (S.2)). The mean value of G is set such that isolated cells maintain the average subconfluent areaĀ SC , G 0 = 1 2Ā SC τ 0 , and the standard deviation is set to 0.05 G 0 . While G sets cellular growth rate (dA T /dt) in an isolated state, tissue crowding reduces the available space for growth to bulk cells, causing their actual area to deviate from the target area. This results in the domination of the exponential term of Eq. (S.2), e −k(A−A T ) 2 , and drives the growth dA T /dt to zero until the cell's area increases. The sensitivity to crowding, represented by the heuristic parameter k, quantifies the allowable deviation between a cell's area and its target area The default value, k = 0.1 pixels −2 , was determined empirically to ensure the tissue growth rate captures the expected growth dynamics over several cell cycles.
In our implementation of the CPM, there are four parameters within the Hamiltonian that set the energy scale; elastic energy has λ and contact energy has J cc , J cm , and J mm . The parameter λ represents the area elastic modulus for the cell, and we choose a default value λ = 1. The parameter J represents the disparity between cell surface tension and intercellular adhesion, implying that higher J values indicate lower cellcell adhesion. Consequently, we chose J cc = 6 and J cm = J mm = 10, which allowed cells to adhere more strongly to each other than the medium, resulting in a circular morphology for the growing colony.

Video Legends
Video S1. Simulation of a growing epithelial colony with apoptosis turned off. The simulation uses default parameters as listed in Table 1. Individual cells are color coded by their growth rates.
Video S2. Simulation of a growing epithelial colony with apoptosis turned on. The simulation uses default parameters as listed in Table 1. Individual cells are color coded by their growth rates.
Video S3. Simulation of a growing epithelial colony exhibiting boundary growth due to strong contact inhibition. The simulation uses default parameters as listed in Table 1, except elastic modulus λ = 1 and sensitivity to crowding k = 1. Left: Cells are color coded by their growth rates, Right: Cells are color-coded by the cell cycle phase (blue: G1 phase, fuchsia: S/G2/M phase).
Video S4. Simulation of a growing epithelial colony exhibiting uniform bulk growth due to weak contact inhibition. The simulation uses default parameters as listed in Table 1