Mechanical coupling of supracellular stress amplification and tissue fluidization during exit from quiescence

Significance Most cells in the human body reside in a dormant state characterized by slow growth and minimal motility. During episodes such as wound healing, stem cell activation, and cancer growth, cells adapt to a more dynamic behavior characterized by proliferation and migration. However, little is known about the mechanical forces controlling the transition from static to motile following exit from dormancy. We demonstrate that keratinocyte monolayers install a mechanical system during dormancy that produces a coordinated burst of intercellular mechanical tension only minutes after dormancy exit. The activated forces are essential for large-scale displacements of otherwise motility-restricted cell sheets. Thus, cells sustain a mechanical system during dormancy that idles in anticipation of cell cycle entry and prompt activation of motion.


Mathematical model
Here we give a detailed description of the model presented in the main text. We assume that the cellular monolayer is a circular layer of radius R and average thickness heq. In previous models of epithelial monolayers both cell based vertex or Voronoi models (1)(2)(3)(4) and continuum models (1,5) have been considered for which its constitutive equations described the monolayer as an elastic (6)(7)(8)(9) or a viscous (10)(11)(12) material. Here, we assume that the monolayer can be described as a two-dimensional active gel (13,14) with a displacement field u(x, t) with x being the spatial coordinate, the bold style indicates vector notation and t is time. The displacement at a specific position is coupled to an activation concentration c(x, t), which regulates the contractile stresses within the cellular monolayer e.g. actinomyosin, which we adopt into the model to model the experimental measured behavior. Additional contributions to the model such as polarization of the cells are discussed in a following section. With negligible inertial forces and assuming no shear stress on the monolayer free surface, the in-plane force balance averaged over the two-dimensional monolayer's height is given as where T(x, t) is the traction forces generated by the friction between the monolayer and the substrate and ( , ) is the cellular stress tensor. The traction force is given as T(x, t) = Γ∂tu(x, t) with Γ being the cell-substrate friction coefficient. ∂ is the partial derivative with respect to time and ∇ is the two-dimensional gradient operator. The cellular stress consists of the sum of the passive elastic stress and the active stress due to the activation concentration (15) ( , ) = ( , ) + ( , ). The stress is modeled as a sum of a linearly elastic material σ = ϵ, where E is the Young's modulus, ϵ = 1/2 (∇ ( , ) + (∇ ( , )) ) is the strain tensor and a viscous component σ ( , ) = η ∂ ϵ( , ) with η the cellular viscosity. The active contribution is modeled with a linear dependence to the amount of active contractile units mimicked by the concentration field ( , ) giving ( , ) = ( , ) , where α sets the magnitude of the isotropic contractile force and is the identity matrix. The dynamics of the cellular monolayer is described by the force balance between the intracellular stresses generated by the amount of activated actinomyosin ( ( , )) and the traction forces acting on the monolayer from the underlying substrate. Considering the monolayer to be an active gel, the force balance in the continuum limit is expressed as The dynamics of the concentration field is given as a convection-reaction-diffusion equation τc is the relaxation time scale towards an equilibrium concentration ceq that eventually will lead to a flat homogenous stationary monolayer, β is the rate of production of concentration due to cellular contraction and D is the diffusion coefficient. The terms on the left-hand side of equation [III] describes the rate change and convection of ( , ), respectively. On the righthand side, the first term gives the kinetic growth/decay around an equilibrium concentration. The second term represents the increase or decay in activation during stretching/compression of the cells and the last term describes the diffusion of the signal inside the layer.

Biophysical parameters
All model parameters are cell-type dependent and are chosen accordingly when rescaling the data obtained from the numerical simulations. The monolayer thickness, heq, the monolayer radius, R, and the Young's modulus of the monolayer, E, are measured experimentally and are listed in Supplementary Table 1. The remaining values are chosen within the range reported in the literature (6,9), and their exact value are determined by comparing our numerical results with the experimentally measured time and length scales when mapping out the phase space in cell layer dynamics. . We assume the gradient ∇ to scale as −1 , with being a characteristic length scale of the system and ∂ ∼ −1 and scale ( , ) ∼ c 0 . By using these scaling arguments we can deduce the cellular displacement scales as

Supplementary
[IV] The velocity in the monolayer,  u/t, can then be estimated as At early times and small displacement we see that the monolayer velocity scales linearly with the concentration gradient, as is recovered with a small offset in the numerical simulations and shown in Supplementary Figure 5, where the maximal velocity is plotted as a function of the initial concentration ratio c0. Moreover, as the monolayer contracts the concentration in the compressed regions will be deactivated and thus reducing the monolayer velocity. From equation [V] we see that the time scale for the velocity is (heqE/(ΓL 2 ) ) -1 , which for the values listed in Supplementary Table 1 yields t ≈ 28h in correspondence with the time scale observed in the experiments when using the monolayer radius R as the characteristic length.

Numerical simulations
We non-dimensionalize equations [II] -[III] by introducing the scaling parameters x = x̂R, t = t̂ ΓR 2 /(Eheq), u = ûR and c = ĉceq with the hat notation indicating non-dimensional variables. with the dimensionless parameters η = ηℎ /(Γ 2 ) which is the ratio of viscous damping over friction forces, α̂ = αceq/E being the ratio between the active contractile strength and the elastic stiffness, Γĉ = Eheq/(DΓuceq) controls the strength between the elastic response and the friction due to diffusion of concentration, τĉ = R 2 /(Dτc) is the kinetic time scale towards an equilibrium concentration, and β̂ = βR 2 /(Dceq) is the ratio between the deactivation of concentration due to the monolayer compression and concentration diffusion. We solve the set of equations [VI] -[VII] coupled using an implicit Newton solver from the finite element library FEniCS on a circular mesh. The numerical simulations are initiated with a zero displacement condition but with a randomized concentration to mimic the stress that is generated in the monolayer during quiescence. The randomization of the initial concentration is performed such that each spatial coordinate has a 60% chance of being seeded with a value drawn from a Gaussian distribution. This initial concentration is then multiplied with a pre-factor that corresponds to the starvation period in the experiments, i.e. large pre-factor equals long starvation period. The initial concentration is defined as c0 =∑ĉ(x, t̂ = 0) dA / ceq , with dA the area of the triangulated mesh. We perform 10 simulations from which we compute the average of the system variables. The boundary conditions are set to reflect what is observed experimentally and we therefore use a homogenous Dirichlet condition on the displacement field u(x̂ = ∂Ω, t) = 0 and a homogenous Neumann condition on the concentration, ∇ĉ(x̂ = ∂Ω, t) · n = 0 with n being the normal vector to the domain boundary ∂Ω. We multiply equations with the test functions ϕ(̂), ψ(̂) and integrate by parts over our numerical domain Ω to obtain the variational formulation of our system equations as As our equations are reduced to first order equations, we discretize them using linear elements with the implicit discrete time derivatives as with being the current time step we evaluate, Δt̂ is the time step spacing and , represents the remaining terms in equations [VIII] -[VIIII], respectively.

Parameter sensitivity -contraction center formation
To determine the effect from the non-dimensional system parameters on the dynamics we performed a parameter sensitivity study, where we systematically vary two parameters while keeping all others fixed. All simulations are initialized with the same initial condition depicted in Supplementary Figure 7 with a ratio c0 = 1.2.
The results from the study are shown in Supplementary Figures 8, 9, 10, 13 with the red dotted lines highlighting the parameter range that yields collective motion and formation of a contraction center, and is summarized in the following points:  α̂ determines the magnitude of the monolayer displacement for a given value of c0. A sufficiently large value of α̂ is crucial to obtain the large-scale migration needed to form one big central contraction center.  Large values of τĉ increase the deactivation rate of the concentration to such an extent that collective migration can be halted but very large values are needed to significantly affect the dynamics.
 Increasing values of η changes the time scale of the cellular contraction. We also see that the magnitude of the contraction centers is fairly insensitive to the value of η.
As η does not affect the characteristics of the large-scale contraction as shown in Supplementary Figure 13 we neglect the viscous contribution in the following where the viscous term can also be absorbed into the uncertainty associated with the cell-substrate friction coefficient. Moreover, from our non-dimensional number η = ηℎ /(Γ 2 ) we see, using the parameter values in Supplementary Table 1, that obtaining a value of η = 1 requires in this model a fairly large cellular viscosity. A much smaller value of η is likely relevant in our system that would not significantly affect the dynamical behavior. Although the viscosity does not appear to influence the dynamics of the active contractile motion, it can be predominant in the long time relaxation of the monolayer or in the presence of cell division and apoptosis, which has been observed in other studies (17,18).
To extend our analysis in the absence of viscous forces we plotted the maximum monolayer velocity magnitude, max(|∂̂̂|), from the simulation data that provided the displacement fields in Supplementary Figures 8-10, shown in Supplementary Figure 11. The previously summarized conclusions are further enhanced and we see that the velocity magnitude is approximately unaffected by changes in both τĉ and β̂. Moreover, we see that the monolayer velocity magnitude increases linearly with α̂ and the square root of Γ̂c, making them both important to achieve large-scale collective migration.
We can thus separate the four dimensionless parameters into two categories, with α̂ and Γ̂c being the parameters important for initiating the collective migration, and τĉ and β̂ are important to slow down the migration at late times, preventing possible diverging compression in the monolayer and to enforce the inevitable relaxation to a flat homogenous equilibrium monolayer at t̂ → ∞.

Numerical solution in a square geometry
In order to check if the collective migration towards the center of the circular mesh is a feature induced by the geometry, we solve the same system on a square mesh. In Supplementary Figure  12A the displacement field in the square geometry shows that the formation of a large central contraction center occurs also for non-circular meshes. Moreover, the displacement field magnitude only differs by  10% to that in the circular mesh (Supplementary Figure 12B) given the same non-dimensional parameters α̂ = β̂ = τĉ = Γ̂c = 1.

Model and experiments -comparison
We performed numerical simulations using a wide range of initial concentration ratios and verified that an initial concentration can cause a collective cellular migration behavior towards the center of the monolayer. The initial response to the random concentration is to form small local contraction centers of high cell density, where the peaks in the initial concentration are located, accompanied by a large traction force pointing in towards these centers. After these initial density centers were formed the monolayer displacement was directed inwards toward the center of the monolayer, where one large cell density center was formed. Below a threshold ratio value c0, there was no collective motion towards the center but the monolayer relaxed to equilibrium after the initial formation of the small density centers (Supplementary Video 7). Around the threshold ratio value stronger local density centers were formed with a subsequent collective motion towards the center of the mesh. However, with small gradients in the concentration the dynamics stagnated before a single central contraction center was formed (Supplementary Video 9). Above the threshold value a grand collective migration occurred towards the center of the mesh, after the initial formation of small local contraction centers. At these conditions the dynamics did not stagnate until a high density central contraction center was formed (Supplementary Video 8). Thus, the dynamics produced in our numerical simulations using a sufficiently high value of the initial concentration ratio c0 recapitulates the dynamic behavior observed in serum-stimulated quiescent cell monolayers (Figure 1, Supplementary Video 2). The traction force was shown to peak at the early formation of the small local contraction centers and decay as the cells began to migrate collectively, which is in agreement with the experimental observations (see Figure 2A and F). The results from the numerical simulations are shown in Figure 4J. These data also highlights that the traction force magnitude increases with the initial concentration ratio c 0 . The overall traction force profile is qualitatively similar to that measured in the experiments, shown in Figure 2F-G. In addition to the increased traction force magnitude, a larger value of c 0 also indicated an increased migration velocity. The averaged radial migration velocity calculated from the numerical data is shown in Figure 4K. The data from the numerical simulations exhibit an increase in the radial velocity for initial concentrations ( 0 ) larger than 1, followed by saturation in migration speed for values above 1.3. From the numerical data we could fit a sigmoid function that describes the migration velocity as a function of initial concentration ( Figure 4K). To bridge the gap between the initial concentration ( 0 ) in the model and the starvation period in the experiments we plotted the extracted sigmoid function together with the experimental observations representing cell migration velocity versus quiescence depth (starvation period) in Figure 5B. This demonstrates a good agreement between the experimental results and the phenomenological model. We also perturbed our system to mimic the monolayer dynamics in the presence of actinomyosin inhibitors. To achieve this, we reduced the value of the non-dimensional parameter α, which controls the contractile strength in the model. During the parameter study (see Supplementary Figure 8 and Figure 4G-I) we show that our model exhibits no large-scale migration for small values of α, and in Supplementary Figure 4 we demonstrate that reducing the value of α leads to decreased traction force magnitude, which is consistent with the experimental observations presented in ( Figure 3A and B and Figure 5A).

Polarization effects
From the experiments there is no indication that the collective motion and formation of a contraction center is driven by cellular polarization. However, we know that cellular polarization can have an effect on the monolayer dynamics (3, 6-9, 12, 19).
To investigate the effects that cellular polarization has on the dynamics we introduced a polarization vector p(x, t). Following (9), the dynamics of the polarization field can be described by the following equation with a controlling the rate of relaxation towards a homogeneously polarized monolayer, κ is the strength of the nearest neighbor alignment, and ω is the coupling coefficient to the concentration gradient, such that local cell polarization points towards regions of high concentration. Further, we coupled the polarization field to the monolayer displacement through the traction force, which now reads T(x, t) = Γ∂tu(x, t) -p(x, t), with f being the coupling coefficient between the cell polarization and the monolayer displacement. Inserted into equation [I] we get Γ ( , t) = ℎ ∇ ⋅ ϵ( , t) + ℎ ∇ ( , t) + ( , t).

[XIII]
We with ̂= fR/(heqE) being the strength of the polarization-displacement coupling, â = aτ the polarization relaxation rate, ̂ = κτ/R 2 the local neighbor alignment coefficient, and ω̂ = ωτceq/R is the strength of the polarization alignment with the concentration gradient. We now have a non-dimensional coupled

three equation system consisting of equations [XIV], [VII] and
[XV] that we solve as described in the Numerical methods section. Using parameter values found in the literature, see Supplementary Table 1 (Supplementary Table 2 for non-dimensional units), we see the results from a numerical simulation in Supplementary Figure 14. At three different times we compared the displacement field with a polarization field to the displacement field without a polarization field using the same initial condition displayed in Supplementary Figure 7. There were no significant differences in displacement magnitude, length scale or time scale between polarized and unpolarized results. This indicates that cellular polarization has no significant impact on the large-scale collective migration observed in this study, which is consistent with the experiments and we use this as a justification to neglect this effect in our theoretical model.

Supplementary Figures
Supplementary Figure 1. Monolayer thickness measurements. Manual measurements of cell sheet thickness were performed on confluent keratinocyte monolayers subjected to serum depletion for two days followed by serum re-stimulation. Graph shows average values (blue dots) of cell sheet thickness at different cell densities. The red line shows the best fit curve based on linear regression.

Supplementary Video 7
Time lapse of the cellular density in the monolayer obtained from a numerical simulation using a normalized initial concentration c0 = 1.002. We observe that there is a small initial rearrangement in the monolayer density due to the inhomogeneous concentration. Due to the small gradients in the concentration there is no collective motion and the monolayer quickly adopts a quasi-static homogenous profile.

Supplementary Video 8
Time lapse of the cellular density in the monolayer obtained from a numerical simulation using a normalized initial concentration c0 = 1.2. After the formation of many strong local contraction centers there is a collective cell migration response towards the center of the monolayer. This collective migration results in the formation of one large global contraction center.

Supplementary Video 9
Time lapse of the cellular density in the monolayer obtained from a numerical simulation using a normalized initial concentration c0 = 1.09. At higher concentration levels we observe coordinated cellular motion in the monolayer after the initial formation of small local contraction centers. However, the concentration gradients are not large enough to stimulate the formation of a global contraction center.