The growth of a tissue, which depends on cell–cell interactions and biologically relevant processes such as cell division and apoptosis, is regulated by a mechanical feedback mechanism. We account for these effects in a minimal two-dimensional model in order to investigate the consequences of mechanical feedback, which is controlled by a critical pressure, *p*_{c}. A cell can only grow and divide if its pressure, due to interaction with its neighbors, is less than *p*_{c}. Because temperature is not a relevant variable, the cell dynamics is driven by self-generated active forces (SGAFs) that arise due to cell division. We show that even in the absence of intercellular interactions, cells undergo diffusive behavior. The SGAF-driven diffusion is indistinguishable from the well-known dynamics of a free Brownian particle at a fixed finite temperature. When intercellular interactions are taken into account, we find persistent temporal correlations in the force–force autocorrelation function (FAF) that extends over a timescale of several cell division times. The time-dependence of the FAF reveals memory effects, which increases as *p*_{c} increases. The observed non-Markovian effects emerge due to the interplay of cell division and mechanical feedback and are inherently a non-equilibrium phenomenon.

## I. INTRODUCTION

Life around us, spanning a bewildering array of length and time scales, is sustained through multicellular processes that are driven by non-equilibrium events such as cell growth and cell division.^{1,2} Although known for a long time,^{3} several recent experimental studies have emphasized that growth and division in cell collectives are governed by local stresses that the cells experience.^{4–7} A manifestation of the coupling of growth and division to local stress is the deviation of the growth law of the cell collective from exponential law.^{6} These experiments suggest that there must exist mechanical feedback between the local stress and cell division.

In a series of papers, we showed that single cells in a collective exhibit anomalous dynamics due to local stress-dependent cell growth and division.^{8–14} In these studies, we have uncovered the mechanism leading to spatially heterogeneous dynamics of cells in cell collectives^{15,16} from jamming to super-diffusion, which is a consequence of dynamical phase separation.^{8,10,11} In the present study, we explore the dependence of mechanical feedback, mediated by stress threshold *p*_{c}, on the dynamics of single cells. A recent interesting study^{17} has shown that mechanical feedback regulates the physical properties of jammed cell collectives, which supports experimental findings.^{4} However, how the mechanical feedback regulates individual cell migration in a collective is unknown and is the problem that we address in this study.

Using a two-dimensional off-lattice agent-based simulation model, we explore the role of mechanical feedback (*p*_{c}) on single-cell dynamics. The central results of the present study are as follows: (i) In the absence of cell growth and division, the dynamics is solely governed by short-ranged two body interactions. In this limit, the cells over the long-time behave like a glass-like solid. (ii) In the presence of cell division, when systematic interactions are absent, the cells exhibit Markov dynamics resulting in diffusive motion at long times. This finding is surprising because there is no thermal motion (temperature is an irrelevant variable). This is different from a free Brownian particle where temperature randomizes the particle motion. (iii) When both cell division and systematic interactions control the collective movement in tandem, the dynamics is regulated by the mechanical feedback that is parameterized using *p*_{c}. To quantify the dynamics, we calculated the force auto-correlation function (FAF) inspired by works in the theory of chemical reactions in liquids.^{18–20} We show that the FAF increases when *p*_{c} is increased. The emergence of long-time correlation in the FAF shows a departure from Markov dynamics and is suggestive of memory effects in growing cell collectives. (iv) The persistence of trajectories of individual cells increases as *p*_{c} is increased. The trajectories are strikingly different from simple Brownian motion. The enhanced persistence in cell dynamics, as *p*_{c} increases, is the origin of memory in active systems. Taken together, the present study establishes how mechanical feedback coupled with cell growth and division leads to non-Markovian cell dynamics whose importance has not been appreciated before.

## II. METHODS

We simulated the spatial and temporal dynamics of a two-dimensional (2D) growing tissue using the agent-based off-lattice model in which the cells are represented as interacting deformable disks. This simplified assumption of representing cells as deformable disks was also used in previous studies,^{21,22} although the details differ. In the model, the cells grow stochastically in time and divide upon reaching a critical size (*R*_{m}), the mitotic radius. The cell-to-cell interaction is characterized by elastic and adhesive forces. We also consider cell-to-substrate damping as a way of accounting for the effects of friction experienced by a moving cell by the substrate.

### A. Physical interactions

Each cell is modeled as a deformable disk with a time dependent radius. A cell is characterized by physical properties such as the radius, elastic modulus, membrane receptor, and ligand concentration. In addition, the cells attract each other through E-Cadherin mediated adhesive interactions. This model is inspired by previous studies on 3D off-lattice multicellular tumor growth models.^{8–13,23,24} The elastic (repulsive) force between two disks with radii *R*_{i} and *R*_{j} is given by

where *E*_{i} and *ν*_{i}, respectively, are the elastic modulus and Poisson ratio of cell *i*. The overlap between the disks, if they interpenetrate without deformation, is *h*_{ij}, which is given by $max[0,Ri+Rj\u2212|r\u20d7i\u2212r\u20d7j|]$ with $|r\u20d7i\u2212r\u20d7j|$ being the center-to-center distance between the two disks.

Cell adhesion, mediated by receptors on the cell membrane, is the process by which cells can attach to one another. For simplicity, we assume that the receptor and ligand molecules are evenly distributed on the cell surface. Consequently, the magnitude of the adhesive force, $Fijad$, between two cells *i* and *j* is expected to scale as a function of their contact line-segment, *L*_{ij}. Keeping the 3D model as a guide,^{8} we calculate $Fijad$ using

where $cirec$$(cilig)$ is the receptor (ligand) concentration (assumed to be normalized with respect to the maximum receptor or ligand concentration so that $0\u2264cirec,cilig\u22641$). The coupling constant *f*^{ad} allows us to rescale the adhesion force to account for the variabilities in the maximum densities of the receptor and ligand concentrations. We calculate the contact length, *L*_{ij}, using the length of contact between two intersecting circles, $Lij=(|4rij2Ri2\u2212(rij2\u2212Rj2+Ri2)2|)rij$. Here, *r*_{ij} is the distance between cells *i* and *j*. As before, *R*_{i} and *R*_{j} denote the radius of cell *i* and *j*.

Repulsive and adhesive forces considered in Eqs. (1) and (2) act along the unit vector **n**_{ij} pointing from the center of cell *j* to the center of cell *i*. The total force on the *i*th cell is given by the sum over its nearest neighbors (*NN*(*i*)),

The nearest neighbors satisfy the condition *R*_{i} + *R*_{j} − |**r**_{i} − **r**_{j}| > 0. Figure 1(a) shows the plot of the total force as a function of inter-cellular distance.

### B. Three scenarios

In order to elucidate the dramatically different dynamical behavior, we consider three limits. (I) The collective movement arises solely from the systematic forces, given in Eq. (3). (II) Cell movement with **F**_{i} = 0 (no inter-cellular interactions) but allowing for cell division and growth. Note that since the inter-cellular interactions are absent, mechanical feedback (*p*_{c}) does not play a role. In this limit, we show that the dynamics can only arise due to active forces generated upon cell division. The limits (I) and (II) are not relevant in describing collective movements in Multicellular Spheroids (MCSs)^{15} or evolving cell monolayers.^{6} (III) In this limit, we not only include interactions between cells [Eq. (3)] but also allow for cell growth, division, and apoptosis. Most importantly, the time-dependent growth of the tissue colony is limited by mechanical feedback, which prohibits the biologically important process of cell growth and division if the local stress on a cell exceeds a critical non-zero value, *p*_{c}.

### C. Equation of motion

The damped dynamics of the *i*th cell, assumed to have two degrees of freedom (*x* and *y* coordinates), is computed based on the equation of motion,

Here, *γ*_{i} is the friction coefficient and **r**_{i}(*t*) is the position of the *i*th cell at time *t*. We assume *γ*_{i} to be equal to *γ*_{o}*R*_{i}(*t*), where *γ*_{o} is a constant. This form of *γ*_{i} is inspired by the simulations of three-dimensional models for solid tumors where *γ*_{i} = 6*πηR*_{i} with *η* being the viscosity. Note that we do not consider the effect of temperature (set to zero in the simulations) as we assume the friction coefficient, that in reality arises from the extracellular matrix in 3D or substrate in 2D, to be so high^{21} that thermal motion is irrelevant. The equation of motion in Eq. (4) is similar to the case for soft granular materials where the role of temperature is neglected.^{25} However, it is crucial to note that in the growth of the tissue colony, scenarios II and III in our case, there is a self-generated active force (SGAF) that arises due to the biologically important processes of cell growth and division.^{10}

### D. Cell growth, division, and apoptosis

In our model, cells can be either in the dormant (*D*) or in the growth (*G*) phase depending on the local pressure associated with a cell [Fig. 1(b)]. Using Irving–Kirkwood definition, we track the pressure (*p*_{i}) experienced by the *i*th cell due to contact with its neighbors.^{26} The expression for *p*_{i} is given by

where $Ai=\pi Ri2$ is the area of the cell. If the local pressure, *p*_{i}, exceeds a critical limit (*p*_{c}) the cell stops growing and enters the dormant phase. Note that the cell can switch back to a growing phase if $pipc<1$ as the tissue evolves. The critical pressure *p*_{c} serves as a mechanical feedback, which is known to regulate the growth of tissues.^{7}

For growing cells $(pipc<1)$, their area increases at a constant rate *r*_{A}. The cell radius is updated from a Gaussian distribution with the mean rate $R\u0307=(2\pi R)\u22121rA$. Over the cell cycle time *τ*,

where *R*_{m} is the mitotic radius. The cell cycle time (*τ*) is related to the growth rate (*k*_{b}) by $\tau =ln2kb$. A cell divides once it grows to the fixed mitotic radius (*R*_{m}). To ensure area conservation, upon cell division, we use *R*_{d} = *R*_{m}2^{−1/2} as the radius of the daughter cells. The two resulting cells are placed at a center-to-center distance *d* = 2*R*_{m}(1 − 2^{−1/2}). The direction of the new cell location is chosen randomly from a uniform distribution on the unit circle. One source of stochasticity in the cell movement in our model is due to random choice for the mitotic direction. In our simulations, the cells may undergo apoptosis at the rate *k*_{a}. Throughout this work, the apoptosis rate was fixed to 10^{−6} s^{−1}. Table I depicts the parameters used in the simulations.

Parameters . | Values . | References . |
---|---|---|

Timestep (Δt) | 10 s | This paper |

Critical radius for division (R_{m}) | 5 µm | 8 and 24 |

Friction coefficient (γ_{o}) | 0.1 kg/(μm s) | This paper |

Cell cycle time (τ) | 54 000 s | 8 and 27–29 |

Adhesive coefficient (f^{ad}) | 10^{−4} µN/μm | This paper |

Mean cell elastic modulus (E_{i}) | 10^{−3} MPa | 8 and 30 |

Mean cell Poisson ratio (ν_{i}) | 0.5 | 8 and 24 |

Apoptosis rate (k_{a}) | 10^{−6} s^{−1} | 8 |

Mean receptor concentration (c^{rec}) | 1.0 | 8 |

Mean ligand concentration (c^{lig}) | 1.0 | 8 |

Parameters . | Values . | References . |
---|---|---|

Timestep (Δt) | 10 s | This paper |

Critical radius for division (R_{m}) | 5 µm | 8 and 24 |

Friction coefficient (γ_{o}) | 0.1 kg/(μm s) | This paper |

Cell cycle time (τ) | 54 000 s | 8 and 27–29 |

Adhesive coefficient (f^{ad}) | 10^{−4} µN/μm | This paper |

Mean cell elastic modulus (E_{i}) | 10^{−3} MPa | 8 and 30 |

Mean cell Poisson ratio (ν_{i}) | 0.5 | 8 and 24 |

Apoptosis rate (k_{a}) | 10^{−6} s^{−1} | 8 |

Mean receptor concentration (c^{rec}) | 1.0 | 8 |

Mean ligand concentration (c^{lig}) | 1.0 | 8 |

### E. Initial conditions

We begin the simulations by placing 100 cells on a 2D plane whose coordinates are chosen from a normal distribution with a mean zero and standard deviation of 25 *µ*m. All the parameters apart from critical pressure (*p*_{c}) are fixed.

## III. RESULTS

### A. Markov dynamics in the presence of delta-correlated random force

For comparison, we briefly summarize the well-known result for a stochastic process (over-damped Langevin equation) in one dimension for a free Brownian particle. The equation of motion is

where the random force obeys ⟨*η*(*t*)⟩ = 0 and ⟨*η*(*t*)*η*(*t*′)⟩ = *δ*(*t* − *t*′). The position of the particle is *x*. The solution *P*(*x*, *t*), for probability density for finding the particle at *x* at time *t*, is given by

Here, we have assumed that the particle was at the origin at *t* = 0, *P*(*x*, 0) = *δ*(*x*). The moments of *P*(*x*, *t*), which serve as the physical observables in cell tracking experiments,^{15} are readily calculated. For instance, the first moment ⟨*x*⟩ is given as

The second moment ⟨*x*^{2}⟩, also called the mean-squared displacement (MSD), is non-zero and is given as

The dynamics of a free Brownian particle is an example of a Markov process. There is no memory because $P(\delta x(t))=N(0,2D\delta t)$ is independent of *x*(*t*), where *δx*(*t*) = *x*(*t* + *δt*) − *x*(*t*) is the displacement of the particle from time *t* to *t* + *δt* and $N(0,2D\delta t)$ is the normal distribution with zero mean and variance 2*Dδt*. This example sets the stage for exploring the emergent force auto-correlation (FAF), with long temporal correlation, in an expanding tissue.

### B. Role of physical interactions and cell division

We first investigate the consequences of physical interactions and cell division when they are not coupled to one another.

#### 1. Physical interactions (limiting case I)

When the interactions between cells are based only on the systematic interactions, as given in Eqs. (1) and (2) without cell division and apoptosis, the dynamics of the interacting cells is governed by the elastic timescale $\gamma ERm$. In the absence of cell growth, division, and apoptosis, the number of cells, *N*(*t*), is a constant, which is confirmed in Fig. 2(a). Figure 2(b) shows the plot of mean-square displacement, Δ(*t*), defined as

where ⟨⋯⟩ denotes the ensemble average over 20 simulation runs and *N* is the initial 100 cells. As in Eq. (4), **r**_{i}(*t*) is the position of the *i*th cell at time *t*. Also, **r**_{i}(0) is the position of the *i*th cell at time *t* = 0. Note that the definition of Δ(*t*) is unchanged when cells undergo growth, division, and apoptosis. Figure 2(b) shows that Δ(*t*) relaxes rapidly to a plateau on a time scale $\u2248\gamma ERm$. Usually, Δ(*t*) ∼ *t*^{α}, with *α* = 0, is indicative of solid-like behavior. Figure 2(b) shows that in the long-time limit, $t\u226b\gamma ERm$, *α* = 0, and hence, the cell collective behaves as a solid in the sense there is absence of diffusion. Note that the dynamics is performed under athermal (temperature is not relevant) open boundary conditions and the scale of systematic interactions are short-ranged $(\u2248Rm)$. Hence, the cells cannot move after the initial relaxation process.

#### 2. Effect of cell division (limiting case II)

In this limit, during each cell division event, a cell is displaced by the distance $\u2248Rm$, randomly in space. Hence, when the time evolution of a cell colony is governed solely by cell division (absence of apoptosis and systematic interactions), we expect that with successive cell divisions, a cell would undergo a random walk that is uncorrelated in time and space. In other words, it would behave as a Brownian particle due to the SGAF induced by cell division. This is purely a non-equilibrium dynamical process. As a result, we expect that Δ(*t*) = *D*_{eff}*t* at long times. Because in this limiting case, *τ* is the only time scale and *R*_{m} is the only length scale, we obtain $Deff=Rm2\tau $. It is interesting to note that the origin of $Deff=Rm2\tau $ is reminiscent of non-equilibrium effects in active systems, and one can define a pseudo-temperature (*T*_{a}) like quantity based on analogies to the Stokes–Einstein relation, $Ta=Deff\gamma KB$. In the absence of systematic interactions, the pressure [Eq. (5)] on an individual cell is zero, and all the cells grow and divide independently. Hence, the number of cells, *N*(*t*), increases exponentially [see Fig. 2(c)]. In this limit,

and, therefore, $N(t)=Noekbt$. Interestingly, in accord with the arguments given above, the dynamics of individual cells is diffusive in this limit and the mean squared displacement Δ(*t*) is given by

as shown in Fig. 2(d). Although there is no thermal motion, the scaling behavior of Δ(*t*) is similar to the standard Brownian dynamics given by Eq. (10). The dynamics is diffusive because during every cell division, the cell is displaced by distance *R*_{m} randomly and, hence, mimics a Brownian motion in two dimensions. In this limit, the growth of the cell does not play a role in determining Δ(*t*) as there are no systematic interactions. Hence, if the size of the cell increases, the position of the cell does not change. Note that when only cell apoptosis is present (*k*_{b} = 0), the cells do not move, and *N*(*t*) decreases exponentially to zero at the rate *k*_{a}.

The dynamics is non-trivial when all the components (systematic forces, cell division and apoptosis, and mechanical feedback) are included. For this case, the growing tissue develops a core where the cells are jammed and exhibit glass-like dynamics.^{11} In contrast, the cells in the periphery are predominantly in the growth (*G*) phase. As a result, the cells exhibit anomalous spatially heterogeneous dynamics with a super-diffusive (sub-diffusive) periphery (core).^{11} Our previous work has shown that this dynamical phase separation arises due to SGAFs that arise due to local stress-regulated cell growth and division.^{10} In the model, cell division and growth are regulated by the mechanical feedback parameter *p*_{c}, the consequences of which are explored in Sec. III C.

### C. Highly correlated force correlations in an expanding tissue (limiting case III)

For a simple Brownian motion, the random force is delta-correlated and, hence, the dynamics is Markovian [see Eq. (7)]. Therefore, a signature of such dynamics is the fast decay of force auto-correlation function (FAF), in comparison to the smallest time-scale of the problem, which in the present case is $\gamma ERm$. To explore the nature of the dynamics, when both cell division and systematic interactions are present in tandem, we calculated FAF (*t**) given as

Here, **F**(*t*) is the force on the cell at time *t*, ⟨⋯⟩_{t} is the time average, and *t** is the delay (or waiting) time. The time averaging is performed over $\u22482000$ cells. Figure 3 shows the plot of FAF for *p*_{c} = 10^{−3}, 10^{−4}, and 10^{−5} Nm^{−1}. Figure 3 shows that the FAF decays on two time scales: long $(\u223c1kb\u2212ka)$ and short $(\gamma ERm)$. In order to extract the two time scales, we fit the FAF with $Ae\u2212t*\tau c+B$ in both the regimes.

In the short time regime (see the inset of Fig. 3), for *p*_{c} = 10^{−3} Nm^{−1}, $A=0.5,\tau c=1.2\gamma ERm$ and *B* = 0.41. For *p*_{c} = 10^{−4} Nm^{−1}, $A=0.75,\tau c=0.97\gamma ERm$ and *B* = 0.16. Finally, for *p*_{c} = 10^{−5} Nm^{−1}, $A=0.81,\tau c=0.95\gamma ERm$ and *B* = 0.11. As anticipated, we find that in the short time regime, the relaxation time is approximately close to the elastic time scale $\gamma ERm$, which is negligible compared to $1kb\u2212ka$. However, in the long time limit, the FAF shows memory effects, especially for *p*_{c} = 10^{−3} Nm^{−1}. For *p*_{c} = 10^{−3} Nm^{−1}, $A=0.41,\tau c=2.2kb\u2212ka$ and *C* = −0.06. For *p*_{c} = 10^{−4}, $A=0.12,\tau c=2.3kb\u2212ka$ and *B* = −0.02. Finally, for *p*_{c} = 10^{−5} Nm^{−1}, $A=0.04,\tau c=0.2kb\u2212ka$ and *B* ≈ 0. For *p*_{c} = 10^{−5} Nm^{−1}, A is negligible, which indicates the absence of memory. In the other two cases, A for *p*_{c} = 10^{−3} Nm^{−1} is four times larger than for *p*_{c} = 10^{−4} Nm^{−1}. Larger magnitude of FAF in the long time regime leads to higher degree of migration for *p*_{c} = 10^{−3} Nm^{−1}. The emergence of highly correlated forces elucidates the departure from Markovian dynamics in a system comprising many interacting cells whose time evolution occurs under non-equilibrium conditions and in the absence of fluctuation-dissipation theorem.^{14}

### D. Persistence of trajectories increases with increasing *p*_{c} (continuation of limiting case III)

To better understand, the significance of memory effects that are embodied in FAF, we probed the trajectories of individual cells in the growing tissue colony. We investigated the three cases: *p*_{c} = 10^{−5}, 10^{−4}, and 10^{−3} N/m in Fig. 4. Note that the trajectories are for the cells, which were present between 1.8 ≤ *t** = (*k*_{b} − *k*_{a})*t* ≤ 3.7. Hence, the trajectories correspond to both the initial parent cells and the subsequent daughter cells. Thus, the differences in the trajectories emerge due to the interplay of systematic interactions, cell growth, and division. Figure 4 shows that the nature of the trajectories is strikingly different when *p*_{c} values are changed. For *p*_{c} = 10^{−5} N/m in Fig. 4(c), the cells move the smallest compared to *p*_{c} = 10^{−4} N/m and *p*_{c} = 10^{−3} N/m and are not persistent in nature. This is due to the fact that the memory effect in the FAF is negligible, as illustrated in Fig. 3. The cells exhibit persistent trajectories for *p*_{c} = 10^{−4} N/m and *p*_{c} = 10^{−3} N/m, with the degree of persistence being higher for the latter. The emergence of persistence in trajectories implies the presence of memory or (non-Markovian) effects in the dynamics, which has also been observed in recent studies of growing cell colonies.^{22} We also calculated the number of cells, *N*, as a function of time for three values of *p*_{c} in Fig. 4(d). For *p*_{c} = 10^{−3} N/m, *N* increases exponentially (linear in a log-linear plot), for *p*_{c} = 10^{−4} N/m, *N* increase exponentially initially but transitions to sub-exponential growth at later times, and for *p*_{c} = 10^{−5} N/m, *N* increases sub-exponentially throughout. To quantify the degree of persistence as a function of *p*_{c}, we also calculated Δ(*t*) ∼ *t*^{α} in Fig. 4(e). In accordance with the results in Figs. 4(a)–4(c), *α* decreases as *p*_{c} is decreased. The dynamics is hyper-diffusive (*α* > 2) for *p*_{c} = 10^{−3} N/m, super-diffusive (1 < *α* < 2) for *p*_{c} = 10^{−4} N/m, and sub-diffusive for *p*_{c} = 10^{−3} N/m. The rich diversity in dynamics arising solely by tuning the value of *p*_{c} is of biological origin and has no obvious analogy to passive abiotic systems.

## IV. CONCLUSION

We have shown that in an evolving tissue colony, the dynamics of an individual cell may be approximated as a stochastic process where the force is correlated over many cell division times. The emergence of correlation in force is a manifestation of memory effects, which is a hallmark of non-Markovian dynamics. Memory effects in a growing tissue arise due to cell division and mechanical feedback. There is no potential or energy function, which solely governs the dynamics in over-damped athermal systems, whose derivative results in these processes. This immediately implies that there is no equilibrium in the growing tissue, and hence, the system is always out of equilibrium. The present study also provides a mechanism for persistent motion observed in many out of equilibrium active matter systems.^{31–33}

It is tempting to describe the simulated time-dependent force autocorrelation using a reduced description, something like the Generalized Langevin equation (GLE), to describe the effective dynamics of a cell in the evolving tissue. It is natural to use such an approach in thermally controlled barrier crossing problems, as was done decades ago in the most insightful studies^{18,19} and also in contemporary studies.^{34–37} For reasons stated above, construction of a similar set of equations, if it exists at all, in any reduced variable (the analog of the reaction coordinate in barrier crossing problems), which must also include large spatial heterogeneity, is likely to be difficult. Cell division and apoptosis require energy input and dissipation, which can only be described using a physical picture and a framework that goes beyond the usual description based on Hamiltonians or energy functions. The work here may provide an impetus to develop a general theoretical framework for describing feedback controlled dynamics in active systems.

## ACKNOWLEDGMENTS

We would like to thank Abdul N. Malmi-Kakkada and Himadri S. Samanta for valuable comments on the manuscript. This work was supported by grants from the National Science Foundation (Grant Nos. PHY 17-08128 and PHY-1522550). Additional support was provided by the Collie-Welch Reagents Chair (Grant No. F-0019).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding authors upon reasonable request.