## Abstract

Morphogenesis often yields organs with robust size and shapes, whereas cell growth and deformation feature significant spatio-temporal variability. Here, we investigate whether tissue responses to mechanical signals contribute to resolve this apparent paradox. We built a model of growing tissues made of fiber-like material, which may account for the cytoskeleton, polar cell-cell adhesion, or the extracellular matrix in animals, and for the cell wall in plants. We considered the synthesis and remodeling of this material, as well as the modulation of synthesis by isotropic and anisotropic response to mechanical stress. Formally, our model describes an expanding, mechanoresponsive, nematic, and active fluid. We show that mechanical responses buffer localized perturbations, with two possible regimes - hypo-responsive and hyper-responsive, and the transition between the two corresponds to a minimum value of the relaxation time. Whereas robustness of shapes suggests that growth fluctuations are confined to small scales, our model yields growth fluctuations that have long-range correlations. This indicates that growth fluctuations are a significant source of heterogeneity in development. Nevertheless, we find that mechanical responses may dampen such fluctuations, with a specific magnitude of anisotropic response that minimizes heterogeneity of tissue contours. We finally discuss how our predictions might apply to the development of plants and animals. Altogether, our results call for the systematic quantification of fluctuations in growing tissues.

Variability has emerged as an inherent feature of many biological systems (1, 2), spanning molecular scales — such as in cytoskeletal dynamics (3) — to tissular scales — such as in organ expansion (4). For instance, cell growth was found to be spatially heterogeneous (5-9), cell cycle length may appear random (10), and there is extensive evidence of stochastic gene expression (11, 12). Such variability has been hypothesised to be required for the emergence of complex shapes since it favors symmetry breaking (13) and self-organisation (14) during development. Nevertheless, growth variability would need to be restrained to ensure robust morphogenesis. In plant tissues, an increase in the spatial correlations of growth fluctuations was shown to reduce the robustness of floral organ size and shape (15). In animal tissues, work on the wing imaginal disc of the fruit fly indicates that robust wing development involves cell competition and requires the modulation of cell division and apoptosis (16, 17).

Mechanical signals are natural candidates for the regulation of growth variability because spatial differences in growth or in deformation rates induce mechanical stress (18-20). In animals, a mechanical feedback affecting the rate of cell divisions was hypothesized (21) and then supported by experiments in Drosophila and in zebrafish (22-25). Actomyosin cables are reinforced by mechanical tension in the wing imaginal disk of Drosophila (26). In plants, mechanical sensing is required to prevent growth fluctuations in roots (27). The deposition of cellulose fibers, the main load-bearing component of the cell wall, depends on wall tension (28, 29), which stiffens the cell wall in the direction of maximal tensile stress (30).

Previous theoretical studies have modelled how mechanical feedback regulates proliferation (21) and how transitions in tissue rheology are induced by proliferation and apoptosis (31, 32). Here, we build upon such studies; in addition, we account for small sources of stochasticity and investigate the consequences on large scale tissue growth. We focus on generic aspects of tissue growth, so that our results may be broadly applicable to active matter (33).

## GROWING TISSUES AS MECHANORESPONSIVE ACTIVE FLUIDS

We built a continuous two-dimensional model of tissue growth. The tissue is assumed to be made of a material with a preferred orientation (i.e. fiber-like), accounting for its main mechanical elements: cytoskeleton, polar cell-cell adhesive junctions, extra-cellular matrix (ECM) in animals; cellulose within the cell wall in plants. Hence, the state of the tissue is locally described by two order parameters, the density of fibers and the nematic field describing the orientation of fibers and their degree of alignment, which confer isotropic and anisotropic mechanical properties to the material, respectively. In our continuous description, we only account for variations in density, orientation, and degree of alignment at supra-cellular length scales, though the material may be patterned at smaller scales (sub-cellular or cellular). We account for fiber synthesis and remodeling, which may be modulated by responses to mechanical stress: reinforcement of actin stress fibers or of the ECM, enhancement of myosin activity, or fluidisation by cell division, in animals; increase in cell wall synthesis, cellulose synthesis, or cell division, in plants. Synthesis has a small random component, considered as a stochastic, uncorrelated source. Stochasticity in synthesis induces growth heterogeneity, which results in mechanical stress and feeds back on synthesis. We use a viscous description of long-term tissue remodeling, so that we cannot account for short-term elastic tissue behavior, which would be better captured by an elastic model as performed in a parallel study (34). Formally, the model describes an expanding, mechanoresponsive, nematic, and active fluid.

### A model of nematic viscous fluid

We describe the fibers with a density and a nematic 2 × 2 tensor , that vary with the position vector . The nematic tensor may be defined as an average over a small region around position , , where the unit vector defines the polarization of fiber monomers and is the unit tensor. Regions of the material move at velocity , which may also vary spatially. In the following, we use the gradient of the velocity field, decomposed into strain rate, , and vorticity , where stands for the partial derivative with respect to position and * ^{T}* for the transpose of the preceding tensor.

We neglect diffusion of fibers in the tissue. The equations of continuity for density and nematic tensor are then
where *t* is time, *κ _{ρ}* is the rate of synthesis of material, and is a nematic tensor that describes the orientation of synthesis and its degree of alignment. The second terms in the left-hand sides of [1-2] account for the dilution of material density and degree of alignment due to tissue expansion; the third term of [2] accounts for the rotation of fibers due to the flow.

Expansion of the tissue is assumed to be driven by a uniform and isotropic tension, *p*, which may correspond to turgor pressure in plants, or to a pressure induced by cell divisions in animals (31); this tension is one of the active components of our model. The mechanical stress, , then follows the force balance equation , supplemented with standard boundary conditions: no shear stress and normal stress equal to tension, *p*, at system boundaries. We consider time scales long enough for tissue remodeling to occur, so that we neglect elastic behavior, assuming that depends on the strain rate tensor, on the density, and on the nematic tensor. This dependence is the constitutive law that characterizes the rheology of the tissue.

In the following, we consider small fluctuations around an average state. The statistical averages of variables are denoted by brackets. For convenience, tensorial fields are decomposed into hydrostatic and deviatoric components, being traceless. On average, the tissue has uniform density, , and is isotropic, ; mechanical stress has only a hydrostatic component ; areal growth rate is on average uniform; through an appropriate change of reference frame, the averaged velocity may be written as . Assuming small fluctuations of all fields, we linearize the constitutive equation as a function of the hydrostatic strain rate, γ, the deviatoric strain rate, , the density, ρ, and the nematic tensor ,
where *η* is an effective viscosity coefficient, *c _{ρ}* and

*c*are effective compressibilities, and

_{s}*v*is analog to Poisson’s ratio. In the right-hand side of [3-4], the first term accounts for viscous-like remodelling and the second for tissue compressibility (changes in density and alignment under stress).

### Activity: mechanical responses and fluctuations

On the one hand, mechanical stress orients cell divisions (22, 23, 35) and plant cell wall reinforcement (30). On the other hand, synthesis of ECM or of cell wall and cytoskeleton polymerization are not uniform in space, having some level of randomness (3, 36, 37). The two classes of phenomena are incorporated in the other active component of our model, namely synthesis. Without loss of generality, synthesis may be written at linear order in fluctuations as
where the first terms in the right-hand side of [5-6] describe the mechanical feedback on synthesis. *τ _{ρ}* and

*τ*are the response times of the mechanical feedbacks,

_{s}*σ*and

_{ρ}*σ*determine the amplitudes of the mechanical feedbacks, and

_{s}*ξ*and are the hydrostatic and deviatoric part of the noise, respectively. Noise is assumed to be white and Gaussian with zero mean, and has extended space correlation characterized by noise strengths

_{ρ}*K*and

_{ρρ}*K*and by a correlation length

_{ss}*ℓ*, which is typically sub-cellular or cellular. The correlations functions of noise take the form , and .

*δ*is the delta distribution and

*g*(

*x*) is a positive function decaying quickly to zero as

*x*→ +∞. (We use

*g*(

*x*) =

*e*

^{−}

*in calculations.) Cartesian coordinates of the 4-tensor*

^{x}**K**are constrained by the traceless nature of to be of the form

_{ss}**K**

_{SS}*=*

_{abcd}*K*{

_{ss}*δ*+

_{ad}δ_{bc}*δ*−

_{ac}δ_{bd}*δ*}, where

_{db}δ_{cd}*δ*is the Kronecker delta. We do not take the limit

_{ij}*ℓ*→ 0 for spatial correlations because otherwise the problem would have no characteristic lengthscale. Accordingly, we hereafter use

*ℓ*as a unit of length.

### Dimensionless parameters

We rescale all fields and variables as follows.

P and are the dimensionless density fluctuation and nematic tensor. They are associated to the dimensionless random components of synthesis and is the dimensionless fluctuation of strain rate tensor, is the dimensionless stress fluctuation. T, , and are, respectively, the dimensionless time, position vector, and velocity fluctuation. The dimensionless versions of Eqs. [1-6] are given in **[SI]**.

This rescaling shows that the model has 8 dimensionless parameters. *ω _{ρ}* = 1 + 1/(2

*τ*〈

_{ρ}*γ*〉) and

*ω*= 1 + 1/(2

_{s}*τ*〈

_{s}*γ*〉) characterize the relaxation of the tissue in absence of mechanical feedback.

*β*

_{0}=

*c*〈

_{ρ}*ρ*〉/(2

*η*〈γ〉) compares the contributions of density variations and growth to mechanical stress.

*v*is the dimensionless difference between effective dilatational and shear viscosities. K

*=*

_{ρρ}*Kρρ*/(16

*η*

^{2}〈

*γ*〉

^{4}) and K

*ss*=

*Kss*/(16

*η*

^{2}〈

*γ*〉

^{4}) are the rescaled magnitudes of random synthesis.

*β*=

_{ρ}*c*〈

_{ρ}*ρ*〉/(2

*τ*〈

_{ρ}*γ*〉

*σ*) and

_{ρ}*β*=

_{s}*c*/(2

_{s}*τ*〈

_{s}*γ*〉

*σ*) are the measures of isotropic and anisotropic responses to stress.

_{s}## RESPONSE TO PERTURBATIONS IN SYNTHESIS

The general formulation is given in the Appendix. Here, we discuss tissue response to an isotropic perturbation that is localised in space – a disk of initial radius *ℓ* – and in time – a duration that is small with respect to all other time scales. Formally, the perturbation to density synthesis is , with H the Heaviside function, while the perturbation to synthesis of nematic order vanishes, . The fields have self-similar forms, , where represents the amplitude of the perturbation and its spatial patterns. The dynamics of the amplitude is specific to each field, whereas the pattern always expands with a characteristic lengthscale *ℓ* exp(〈*γ*〉*t*) (in dimensional units). and are represented in Fig. 1a-e and are explicitly given in **[SI]**. An immediate consequence of the perturbation is to stiffen the tissue, which reduces expansion and increases stress levels ; then the anisotropic mechanical response gradually induces radial fibers and reinforcement in the direction of the main stress. The behavior at longer times depends on the level of mechanical response. For low anisotropic response, i.e. for *β _{s}* smaller than a threshold computed in

**[SI]**, the tissue is hypo-responsive and all amplitudes evolve monotonously as a function of time and vanish at times that are long with respect to the correlation time

*τ*; tissue nematic orientation, strain rate, and mechanical stress are all mainly radial. For high anisotropic response, i.e. for

_{c}*β*above this threshold, the tissue is hyper-responsive, and amplitudes show an underdamped-like dynamics: they change sign before decaying to 0; after well-defined times, density becomes slightly smaller than average density, and all of nematic order, strain rate, and mechanical stress become circumferential. This hyper-responsive regime can be understood as follows. An initially high mechanical anisotropy of the tissue reduces the radial strain until strain becomes circumferential, leading to circumferential stress and then circumferential nematic order.

_{s}These dynamics occur on a time scale *τ _{c}*, which is the maximal relaxation time scale in response to a perturbation

**[SI]**. The time scale

*τ*depends on the magnitudes of mechanical responses, as shown in Fig. 1f (see Fig. S2 for the effect of other parameters). Isotropic mechanical feedback makes perturbations more persistent in time, because

_{c}*τ*increases with

_{c}*β*. The effect of the anisotropic feedback on relaxation is more complex:

_{ρ}*τ*first decrease and then increase as

_{c}*β*is increased; the minimum of

_{s}*τ*corresponds to the transition between hypo-response and hyper-response. This characteristic time

_{c}*τ*will also appear to be important for the effect of noise.

_{c}## GROWTH FLUCTUATIONS

We find that growth (areal strain rate in 2D) has long-range correlations, with a correlation function that spatially decays with an exponent −4/*τ _{c}*

**[SI]**. In practice, growth is measured at the scale of the spatial resolution of experimental measurements, which depends on the landmarks used and is often at cell scale. We therefore define a coarse-grained growth rate and we consider the time correlation function

*G*(R, T) of the growth of a disk of radius R, where R is the coarse-graining size i.e. the resolution size. It is simply related to velocity fluctuations (see Appendix) by . . This correlation function is plotted in Fig. 2. Panel

**a**shows time correlations for high and low anisotropic mechanical feedback, respectively corresponding to the hypo-response and hyper-response. The negative correlations for high feedback are related to underdamped relaxation of hyper-responsive tissues. The correlation function decays quickly to 0 with a characteristic time scale that is exactly the relaxation time,

*τ*, shown in Fig. 2e. Areal growth mean square deviation appears roughly scale-invariant, see Fig. 2b; it is exactly scale-invariant for hypo-response and oscillates around a scale-invariant for hyper-response. Two regimes characterise the decay. In a weakly-correlated regime, when

_{c}*τ*< 2, growth mean square deviation scales with the inverse of the coarse-graining area,

_{c}*G*(R,0) ~ R

^{−2}, an exponent due to the central limit theorem. In a strongly-correlated regime, when

*τ*> 2, growth mean square deviation decays more slowly , see

_{c}**[SI]**for a rationale.

## FLUCTUATIONS OF ORGAN SHAPE

We are now interested in the effects of noise in synthesis on tissue contours or on organ shape. In a homogeneous and isotropic tissue, quantifying the fluctuation of contours is equivalent to determining the fluctuation of a vector joining two landmarks followed throughout growth of the tissue. Hence, we use a Lagrangian description and consider the position, , at time *T* of a landmark initially at position , which is determined by the dimensionless velocity field through . The fluctuations of are computed in the Appendix. Heterogeneity of contours is assessed using the coefficient of variation, , which is plotted in Fig. 3a. Its asymptotic trend for long times depends on the value of the correlation time *τ _{c}*. If

*τ*< 2, then , and for

_{c}*τ*> 2, CV scales for low feedback and oscillate around this trend for high feedback.

_{c}We represent in Fig. 3b the maximal value of the coefficient of variation, normalized by . This enables us to quantify contour fluctuations in the tissue or of organ shapes. In absence of anisotropic feedback, *β _{s}* = 0, we find that heterogeneity increases with isotropic feedback. Accordingly, a positive isotropic feedback maintains perturbations and induces long-range correlations as seen in previous sections. Conversely, negative isotropic feedback dampens perturbations. Whatever the level of isotropic feedback,

*β*, we find that heterogeneity, as a function of anisotropic feedback, has a single minimum, which corresponds to the transition between hypo- and hyper-response. In the hypo-responsive regime, increasing anisotropic feedback dampens perturbations, whereas in the hyper-responsive regime, increasing anisotropic feedback enhances perturbations due to the oscillatory overshoot. Finally, we note that the behavior of heterogeneity in Fig. 3b is qualitatively similar to the behavior of correlation time (

_{ρ}*τ*) in Fig. 1f, indicating that the correlation time is a major determinant of heterogeneity because the correlation time sets how the tissue keeps the memory of its previous state.

_{c}## DISCUSSION

We built a continuous viscous model of tissue growth, describing density and nematic order of the tissue, and modelled material synthesis and responses to mechanical stress. Here, the responses are characterized by two parameters, *β _{ρ}* and

*β*, corresponding to isotropic response - increase in density due to increase in stress when

_{s}*β*> 0 - and anisotropic response - increase in tissue anisotropy due to increase in stress anisotropy when

_{ρ}*β*> 0, and conversely when these parameters are negative. In plants, it is believed that cell wall synthesis is enhanced when tension increases (38), which corresponds to

_{s}*β*> 0. The alignment of cortical microtubules with maximal stress orientation leads to the anisotropic stiffening of the cell wall in this direction (30, 39), while cell divisions are associated with new cell walls oriented in the direction of maximal stress (35); both processes yield

_{ρ}*β*> 0. In animal tissues, experiments indicate that tissues are fluidised by cell divisions (22-25): proliferation is enhanced by tensile stress and daughter cells tend to separate along the direction of highest mechanical stress, which corresponds to

_{s}*β*< 0 and

_{ρ}*β*< 0, respectively. At shorter time scales, actomyosin cables are reinforced in the direction of applied stress (26), which yields

_{s}*β*> 0. At intermediate time scales, ECM would also be reinforced against mechanical stress (

_{ρ}*β*> 0), though its role in morphogenesis has not received attention until recently (40-43).

_{ρ}In this study, we determined tissue response to a localized perturbation, depending on the mechanical feedback parameters *β _{ρ}* and

*β*. We generalized predictions that anisotropic mechanical feedback buffers such a perturbation (29), in agreement with observations on trichomes – a cell type with transient faster growth -–in Arabidopsis sepals (29). Here, we unravelled two possible regimes: hypo-response at low anisotropic feedback – perturbations decay monotonously and hyper-response at high anisotropic feedback – perturbations oscillate before decaying, with a characteristic time that is minimal at the transition between the two regimes. This case study provides an assay of mechanical responses in both plant and animal systems, for instance by inducing clones with altered growth rate and quantifying the relaxation timescales in backgrounds with different levels of mechanical response.

_{s}We then investigated the statistical properties of tissue growth, unravelling long-range correlations, with slowly decaying correlation functions. To test this, it would be crucial to examine correlation functions in live imaging data of growing organs, e.g. (15, 44-47). Given that a larger correlation time (*τ _{c}*) corresponds to longer-range correlations, higher levels of anisotropic feedback (independently of its sign) or higher levels of positive isotropic feedback yield more slowly decaying correlations functions. Nevertheless, long-range correlations could be mediated by both chemical and mechanical signals. Further experiments would be required to test whether mechanical signals are involved in such correlations.

Finally, we found that heterogeneity of contours and shapes is minimal for a well-determined level of anisotropic mechanical response. This generalizes a similar conclusion reached for local heterogeneity using a cell-based toy model (48). Here we also accounted for isotropic mechanical responses and considered heterogeneity at all scales. We identified the correlation time as a key parameter determining the extent of spatial correlations and the level of heterogeneity of organ shape. Based on our results, we make the following predictions. In plants (*β _{ρ}* > 0 and

*β*> 0), heterogeneity in development can be significantly high unless anisotropic feedback is close to the value that minimizes heterogeneity. In animals, if we discard possible contributions of the ECM,

_{s}*β*< 0 and

_{ρ}*β*< 0 at long time scales, heterogeneity in development is minimal when anisotropic feedback is negligible. Such predictions can be tested by measuring correlations of growth in space and time, as well as the strength of isotropic and anisotropic variations in growth due to an external force or to the induction of clones. More generally, characterizing the fluctuations of cell properties appears as a promising avenue to shed light on how signals orchestrate organismal development.

_{s}## APPENDIX

### Response to perturbations in synthesis

We assume the tissue to be infinite, the perturbations to not induce rotation , and the reference frame to satisfy . We investigate the response to generic perturbations. Because the average strain rate profile stretches patterns, we consider a modified Fourier transform defined as
with the position rescaled by the average growth factor e^{T/2}. Defining as the direction of the wavevector, the Fourier transform of the nematic tensor can be written as . The linear response for material density and nematic tensor is then given by
where is an exponential involving the relaxation matrix , see **[SI]**, and the fields and are the components of the noise Fourier transform. Finally, the Fourier transform of the strain rate tensor is decomposed as , so that the linear response of strain rate is given by
where the integrand in the r.h.s. of [7] is a sum over the values {+, −} of the index *ε* and *ω _{±}* are the two eigenvalues of [

*ω*]. The expression of

*ω*and the coefficients can be found in

_{±}**[SI]**. We thus obtain the full response of the tissue to any perturbation of synthesis, in terms of the modified Fourier transform of the sources of density and of nematic order.

#### Growth fluctuations

Using the linear response of flow velocity to synthesis perturbations [7-8], we derived the velocity fluctuations, as detailed in **[SI]**. The correlation tensor of velocity fluctuations is proportional to the unit tensor,
for ΔT > 0, with the same coefficients as in [7-8] and
where . This explicitly defines the correlation functions.

#### Fluctuations of organ shape

We look for the probability that a material point follows a path . For small fluctuations, this probability can be simplified as
where the velocity correlation function is given by [9]. We determined the asymptotic statistics of the Lagrangian flow by applying the saddle point method (49) to is maximized by the average trajectory and the correlation tensor of the position is given by [SI]
with T_{1} ≤ T_{2}. This correlation tensor does not depend on T_{2} − T_{1} because of the lack of time translation invariance.