Mechanical feedback controls the emergence of dynamical memory in growing tissue monolayers

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,³ 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.⁶ 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¹⁷ has shown that mechanical feedback regulates the physical properties of jammed cell collectives, which supports experimental findings.⁴ 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.

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.

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−|r⃗i−r⃗j|]$ with $|r⃗i−r⃗j|$ 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,⁸ 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≤cirec,cilig≤1$⁠). 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−(rij2−Rj2+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 ith 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.

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)¹⁵ or evolving cell monolayers.⁶ (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.

The damped dynamics of the ith 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 ith cell at time t. We assume γ_i to be equal to γ_oR_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²¹ 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.²⁵ 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.¹⁰ 

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 ith cell due to contact with its neighbors.²⁶ The expression for p_i is given by

where $Ai=π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.⁷ 

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 $Ṙ=(2πR)−1rA$⁠. 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 $τ=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_m2^−1/2 as the radius of the daughter cells. The two resulting cells are placed at a center-to-center distance d = 2R_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⁻⁶ s⁻¹. Table I depicts the parameters used in the simulations.

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.

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,¹⁵ are readily calculated. For instance, the first moment ⟨x⟩ is given as

The second moment ⟨x²⟩, 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(δx(t))=N(0,2Dδ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δt)$ is the normal distribution with zero mean and variance 2Dδt. This example sets the stage for exploring the emergent force auto-correlation (FAF), with long temporal correlation, in an expanding tissue.

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

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 $γ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 ith cell at time t. Also, r_i(0) is the position of the ith 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 $≈γERm$⁠. Usually, Δ(t) ∼ t^α, with α = 0, is indicative of solid-like behavior. Figure 2(b) shows that in the long-time limit, $t≫γ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 $(≈Rm)$⁠. Hence, the cells cannot move after the initial relaxation process.

In this limit, during each cell division event, a cell is displaced by the distance $≈Rm$⁠, 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_efft 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τ$⁠. It is interesting to note that the origin of $Deff=Rm2τ$ 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γ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.¹¹ 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).¹¹ 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.¹⁰ 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.

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 $γ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 $≈2000$ cells. Figure 3 shows the plot of FAF for p_c = 10⁻³, 10⁻⁴, and 10⁻⁵ Nm⁻¹. Figure 3 shows that the FAF decays on two time scales: long $(∼1kb−ka)$ and short $(γERm)$⁠. In order to extract the two time scales, we fit the FAF with $Ae−t*τc+B$ in both the regimes.

In the short time regime (see the inset of Fig. 3), for p_c = 10⁻³ Nm⁻¹, $A=0.5,τc=1.2γERm$ and B = 0.41. For p_c = 10⁻⁴ Nm⁻¹, $A=0.75,τc=0.97γERm$ and B = 0.16. Finally, for p_c = 10⁻⁵ Nm⁻¹, $A=0.81,τc=0.95γ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 $γERm$⁠, which is negligible compared to $1kb−ka$⁠. However, in the long time limit, the FAF shows memory effects, especially for p_c = 10⁻³ Nm⁻¹. For p_c = 10⁻³ Nm⁻¹, $A=0.41,τc=2.2kb−ka$ and C = −0.06. For p_c = 10⁻⁴, $A=0.12,τc=2.3kb−ka$ and B = −0.02. Finally, for p_c = 10⁻⁵ Nm⁻¹, $A=0.04,τc=0.2kb−ka$ and B ≈ 0. For p_c = 10⁻⁵ Nm⁻¹, A is negligible, which indicates the absence of memory. In the other two cases, A for p_c = 10⁻³ Nm⁻¹ is four times larger than for p_c = 10⁻⁴ Nm⁻¹. Larger magnitude of FAF in the long time regime leads to higher degree of migration for p_c = 10⁻³ Nm⁻¹. 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.¹⁴ 

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⁻⁵, 10⁻⁴, and 10⁻³ 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⁻⁵ N/m in Fig. 4(c), the cells move the smallest compared to p_c = 10⁻⁴ N/m and p_c = 10⁻³ 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⁻⁴ N/m and p_c = 10⁻³ 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.²² 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⁻³ N/m, N increases exponentially (linear in a log-linear plot), for p_c = 10⁻⁴ N/m, N increase exponentially initially but transitions to sub-exponential growth at later times, and for p_c = 10⁻⁵ 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⁻³ N/m, super-diffusive (1 < α < 2) for p_c = 10⁻⁴ N/m, and sub-diffusive for p_c = 10⁻³ 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.

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.

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).

The authors have no conflicts to disclose.

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