Suppression of soft spots and excited modes in the shape deformation model with spatio-temporal growth noise

The shapes of plant organs such as sepals are precise and reproducible, although the cells composing these organs are highly heterogeneous. We investigate the collective behavior of a shape deformation model with spatio-temporal fluctuations in elasticity (Hong et al., Dev. Cell, 38, 15, 2016). It is shown that the spatio-temporal growth noise has two opposing effect: the temporal noise promotes the reproducibility of the organ shape, whereas the spatial noise disturbs it, leading to the organ undergoing an abnormal shape deformation. To understand how such abnormality can emerge, we study the eigenvectors of the correlation matrix of the growth velocity. It is revealed that the anomalous shape deformations can be created by two independent factors: (1) excited modes where the orientations of the eigenvectors with m ≥ 2 are spatially correlated, and/or (2) soft spots where spot structures with high fluctuations in the magnitudes of the eigenvectors appear non-uniformly. © 2019 The Author. Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
Most organisms can reproduce uniformly and precisely shaped organs during their developmental processes ( Hong et al., 2016 ;Lander, 2011 ), even though the cell sizes and growth rates within an organ are highly variable and heterogeneous ( Meyer and Roeder, 2014 ;Roeder et al., 2010 ;Elsner et al., 2012 ;Tsugawa et al., 2017 ). This has underpinned a longstanding biological mystery because non-uniform growth at the cellular level and uniform growth at the organ level appear to be inconsistent. To resolve this conflict between different scales, the mechanical properties of the relevant cells have recently been investigated ( Niklas, 20 0 0 ;Niklas and Spatz, 2012 ;Uyttewaal et al., 2012 ;Sampathkumar et al., 2014 ). For example, it has been shown that in addition to gene regulation, mechanical regulation via stress-microtubule feedback plays a significant role in cell shape deformation ( Uyttewaal et al., 2012 ;Sampathkumar et al., 2014 ). Likewise, recent studies have established mechanical-based models to explain organ shape anisotropy ( Uyttewaal et al., 2012 ;Hervieux et al., 2016 ). Furthermore, newly developed mathematical models incorporating real data on organ shapes have begun to uncover an underlying mechanism in complicated organ shape deformations based on gene and mechanical regulation ( Coen et al., 2004 ;Hayakawa et al., 2016 ;Rebocho et al., E-mail address: stsugawa@bs.naist.jp 2017 ; Coen et al., 2017 ). Therefore, to understand organ shape reproducibility, a mathematical model that accounts for the mechanical properties of cells is required.
In addition to the mechanical properties of the cells, collective behavior in growth processes involving various cells has recently been intensively studied ( Vicsek and Zafeiris, 2012 ). Mathematical models of collective cell motility can be classified into five approaches, according to the cellular objects that these models describe: (1) the center of mass ( Bretschneider et al., 1999 ;Palsson and Othmer, 20 0 0 ; Morishita and Iwasa, 2008 ), (2) domains ( Graner and Glazier, 1992 ;Glazier and Graner, 1993 ;Szabó et al., 2010 ), (3) vertices ( Nagai and Honda, 2001 ;Honda et al., 2004 ;Fozard et al., 2013 ), (4) elastic deformations ( Hong et al., 2016 ;Uyttewaal et al., 2012 ;Hervieux et al., 2016 ;Hamant et al., 2008 ;Bonazzi et al., 2014 ;Boudon et al., 2015 ;Bozorg et al., 2016 ) and (5) active processes ( Goriely et al., 2019 ) of the cells. Among these, the fourth approach has a physical advantage in its connection to the stress feedback, particularly in plant science ( Hervieux et al., 2016 ;Bonazzi et al., 2014 ). In our previous study concerning the fourth approach ( Hong et al., 2016 ), an organ deformation model with spatio-temporal growth noise explained that a reproducible organ shape can be achieved through the relaxation of temporal and spatial fluctuations in elasticity. That means that even in the presence of spatio-temporal growth noise at the cellular level, the growth noise can be canceled out spatially and temporally by averaging the different growth rate and/or different growth directions. We therefore proposed a developmental mechanism that the spatio-temporal randomness at the cellular level can lead to the organ size and shape reproducibility. However, there remains a need to clarify which of the temporal and spatial effects is essential for organ reproducibility. We note that the idea of the spatio-temporally randomized elasticity in the present paper is based on the experimental observations of the spatio-temporal fluctuation in elasticity in the actual sepal data in the previous work ( Hong et al., 2016 ).
In the context of collective behavior, a useful analytical framework for trajectories of elements has recently been discussed ( Bottinelli et al., 2016 ;Toner and Tu, 1998 ;Ramaswamy, 2010 ). To understand the spatio-temporal effects of elastic properties, the dynamical properties of the trajectories of elements are important. For example, the velocity correlation matrix enables the determination of the mutual dependence of the speeds of the two arbitrary elements ( Bottinelli et al., 2016 ). Furthermore, in the theory of disordered materials, the eigenmodes of the correlation matrix can fully characterize the structural stability of a state, including mode softness (fluctuations in a displacement field) and mode extension (spatially correlated localized motion) ( Brito et al., 2010 ;Ashton and Garrahan, 2009 ;Manning and Liu, 2011 ). In particular, the mode with the largest eigenvalue plays a more significant role because it represents the softest mode of the system. Therefore, it is possible to characterize shape deformation dynamics by specifying the mutual dependence and fluctuation of growth velocity.
This study uses methodologies such as the correlation matrix and eigenmode analysis, which have been previously employed in studies regardng the collective behavior of active matter ( Bottinelli et al., 2016 ;Henkes et al., 2012 ). The effects of temporal and spatial noises are investigated using a shape deformation model similar to that used in Hong et al. (2016) , and shape deformation processes are characterized.

Model
We reimplement the two-dimensional elastic model developed in Hong et al. (2016) using the finite element method. The organ shape is modeled using a plate with a leaf-like shape (initially this is a half disk with a radius r = 1 . 1 (a.u.) being consistent with Hong et al. (2016) ; see Fig. 1 A), and its interior is assumed to be filled with small triangles. We use triangular elements based on the observation of subcellular heterogeneity in terms of elasticity in the real sepal data ( Hong et al., 2016 ). We labeled the triangle vertices without overlap as k ( k = 1 , · · · , K). We assume a single cell as approximately six triangles but we do not explicitly explore the cell contours. Cell division is also not considered because recent observations show that the cell growth with or without division do not differ significantly ( Tsugawa et al., 2017 ). The straight edge of the half disk is referred to as ∂ 1 , and its arc as ∂ 2 ( Fig. 1 A). We take the plane of the half disk as the xy -plane, with the x -axis parallel to the straight edge ∂ 1 and the origin located at the center of the half disk.
The time evolution is defined as follows. An external force, corresponding to the turgor pressure p , is applied to the arc ∂ 2 of the half disk. More precisely, the external force vector P on ∂ 2 is given by P = p n , where the unit vector n is the outward normal to the arc, and the force per unit area p is constant. In other words, the plate expands as if the water-balloon swells due to the turgor pressure. On the other hand, the straight edge ∂ 1 corresponding to the connected part to the stem is assumed to be fixed, i.e., the displacement u is zero on ∂ 1 : u x = u y = 0 . Under these boundary conditions, the half disk swells into a different shape for which a mechanical equilibrium is satisfied. At equilibrium, the surface stress balances the external force P . This final shape is determined by the elastic properties of the half disk (see below). Next, the additional external force P = p n is applied to the arc ∂ 2 with a fixed boundary on ∂ 1 . This additional external force on the boundary can be thought of as the increase of the turgor pressure corresponding to an additional water supply from the stem, or the relative force change due to weakening of the mechanical properties of the cell walls with same turgor pressure. Because the surface stress balances the external force at the previous equilibration step ( p n ), we assume the deformed shape as the input without tension and the external force P = p n is applied, which imitates plastic deformation in plant ( Lockhart, 1965 ;Ortega, 1985 ). As a result, the plate swells into yet another shape. These procedures are repeated, with the number of swelling processes denoted by t ( t = 1 , 2 , · · · , T ) , constituting the time variable in the model ( Fig. 1 B). We note that Fig. 1 B demonstrates the time evolution of the mesh of the plate in the case of one of the parameters (ii) with non-uniform elastic properties as defined below.
The elastic properties of the plate are defined as follows. Because the anisotropic growth of an organ can be modeled by mechanical anisotropy, we use a generalized Hooke's law for the relation between the stress tensor σ i j ( i, j = x, y ) and strain tensors ε kl ( k, l = x, y ) , as in ( Landau and Lifshitz, 1986 ): where D ijkl is the elastic modulus tensor of the infinitesimally small region around the k th vertex, which can be defined by Thus, the elastic properties around the k th vertex can be determined by the Young's modulus E k , Poisson's ratio ν, and anisotropic parameters α and β. The anisotropic parameters and direction are assumed to be constant. To consider the elastic properties of each vertex, the Young's modulus of an arbitrary point inside the triangle is calculated by a linear interpolation of the vertices of the triangle, as described in the supporting material. Then, the equilibrium shape of the plate can be determined by div σ = 0 under the boundary conditions mentioned above. To solve this equation numerically, we employ the principle of virtual work, described as where is the whole inner region of the plate and v is a virtual displacement (see the details in the supporting material). To determine the deformation u that satisfies Eq. (3) , we used the FreeFem ++ software platform ( Hecht, 2019 ).
Then, we assume that the Young's modulus E k (t) around the k th vertex varies in time and in space. We note that the shape deformation without the spatio-temporal variability in elasticity does not generate the organ shape variability ( Hong et al., 2016 ), which we call normal shape deformaiton in Fig. 1 C. The initial value of the Young's modulus, E k (0) , is assumed to follow a Gaussian probability distribution with a fixed mean μ 0 ( = 3 . 27 MPa ) and standard deviation (SD) s 0 ( MPa ). Hereafter, the unit MPa will be omitted for ease of reading. We incorporate the spatio-temporal noise as where the random elasticity ξ k (t) is chosen from the Gaussian probability distribution with fixed mean μ 0 and SD s 0 . The temporal parameter q ∈ [0, 1] means the contribution of the random elasticity to the current elasticity so that the value 1 − q indicates the memory of the elasticity. In other words, the value 1 − q can be thought of as the parameter of the relaxation time of the elastic property. We note that the temporal parameter q ∈ [0, 1] is not introduced independently from the spatial noise. Here, q = 0 for the case with a full memory, whereas q = 1 for the case with no mem- with a stationary value V st q = qs 2 0 2 −q (in this calculation of the variance the threshold is not taken into account, because the threshold only has a minor effect on the parameter range used in this study). Thus, E k ,q (t) undergoes a non-stationary process. A minimum threshold of 0.1 is introduced for the Young's modulus E k ,q (t) , in order to ignore an extremely small and negative elasticity. More precisely, in Eq. (4) the noise ξ k (t) is regenerated until we obtain a Young's modulus E k ,q ( t + 1 ) that is larger than 0.1.
Changes in the threshold values with 0.001, 0.01 do not significantly affect the results.

Velocity correlation matrix
To understand the mechanism of the shape deformation, we consider the trajectories of the small triangles (pieces of cells). The mesh is not reconstructed at each equilibration step to quantitatively track the spatio-temporal growth velocity. As described in the introduction, the dynamical properties of the triangle trajectories have potential to characterize shape deformation dynamics. Using velocity correlation matrices, the collective growing trajectories can be decomposed into different types of growth modes with meaningful information. First, the matrix eigenvalues provide the contribution ratio of the growth modes, i.e., the largest eigenvalue corresponds to the softest mode of the system. Second, the eigenvector orientations at a specific growth mode includes spatial correlation in terms of the growth velocity within the mode. Third, the eigenvector magnitudes at a specific growth mode characterize the degree of the speed fluctuation. Therefore, it is possible to identify the mechanism of the shape deformation from those quantities.
The components of the velocity correlation matrices R x and R y can be defined as where For harmonic oscillations in terms of Newtonian dynamics, vibrational eigenmodes (i.e., normal modes) capture the intrinsic low-dimensional motions of a system ( Brito et al., 2010 ;Henkes et al., 2012 ). However, the system in this study is non-equilibrium although it reaches a mechanical equilibrium at each step, so we apply a principal component analysis without expecting to relate the eigenvalues of the correlation matrix to the vibrational properties, as further discussed in Bottinelli et al. (2016) .

Organ shape reproducibility
We focus mainly on two parameters: the spatial randomness s 0 of the elasticity and the temporal randomness q . The parameter ranges are selected as s 0 ∈ [1, 3] and q ∈ [0, 0.1], for the reason that a qualitative change can be observed in these ranges. The other parameters are fixed as p = 0 . 068 ( MPa ) , α = 0 . 15 , β = 0 . 9 , and ν = 0 . 45 if there are no further explanations. The effects of these parameters are summarized in the supporting material.
Three characteristic dynamics are illustrated in Fig. 1 D, E, and F.
First, for ( s 0 , q ) = ( 2 . 8 , 0 . 09 ) (case (i) in Fig. 2 ), the initial shape is conserved during the time evolution with high spatial and high temporal fluctuations in elasticity. Second, for ( s 0 , q ) = ( 2 . 8 , 0 . 01 ) (case (ii) in Fig. 2 ), the plate deforms into anomalous shapes with high spatial and low temporal fluctuations. Third, for ( s 0 , q ) = ( 1 . 2 , 0 . 01 ) (case (iii) in Fig. 2 ), the shape of the plate deforms to some extent with low spatial and low temporal fluctuations. The dependence of the shape on the two parameters is captured by the phase diagram at t = 100 in Fig. 2 . The effect of the mesh size on the final shape variability are summarized in the supporting material. Note that the shape variability increases as a function of the mesh size, except for parameter regions with small q , in which an abnormal deformation appears (see the supporting material). This For each parameter set ( s 0 , q ), shapes at t = 100 are displayed for seven different random seeds. The higher the spatial randomness and the lower the temporal randomness, the greater the shape variability.
indicates that the mesh size also constitutes a system parameter. However, as the final shape do not change significantly for different spatial correlation length in elasticity with same mesh size (#TRI 93) and for different mesh size with same spatial correlation length (Fig. SI. 3 and 4), we mainly focus only on the parameters s 0 and q with #TRI 93. Now, to determine how such an abnormal shape deformation in the parameter region (ii) can emerge, we analyze the collective cell motility through an eigenmode analysis of the correlation matrix of the cell growth velocity.

Effects of spatial and temporal noise: eigenvalues
For comparison, we consider the simplest case of "uniform growth" of the half disk. To describe the uniform growth of the shape, we simplified the boundary condition accordingly, with u y = 0 at ∂ 1 and only the origin of the half disk fixed, and with α = 0 , β = 1 , and s 0 = 0 , that means that this model does not include mechanical anisotropy and spatio-temporal variability.
The eigenvalues λ x m and λ y m averaged over 30 different random seeds are illustrated in Fig. 3 A and B as functions of the mode number m . The eigenmodes are arranged in descending order of the eigenvalues, and thus m = 1 corresponds to the mode with the largest eigenvalue. The eigenvalues λ x m and λ y m for uniform growth are considerably smaller than in these three cases, and only a few eigenmodes are dominant. The eigenvalues in case (i) are larger than in cases (ii) and (iii), indicating that the deformation velocity in (i) fluctuates more than in (ii) and (iii). These results are consistent with the fact that the random elasticity in case (i) fluctuates more rapidly than in cases (ii) and (iii). However, in contrast to these microscopic fluctuations, the macroscopic final shapes for case (i) are less diverse than those of cases (ii) and (iii). This apparent contradiction between the microscopic randomness and macroscopic uniformity is reminiscent of the biological mystery described in the introduction. Therefore, we explore the mechanism of such macroscopic uniformity from the microscopic randomness in the next subsections.

Excited mode: orientations of eigenvectors
To resolve the contradiction mentioned above, we focused on one specific random seed for the uniform growth and cases (i) and (ii), as shown in Fig. 4 A, E, and I, respectively. The results for case (iii) are qualitatively similar to those of case (i). In Fig. 4 D, H,  tor in each panel, which corresponds to a dominant deformation, can be thought of as a "collective expansion", where the triangles appear to move collectively as one object. On the other hand, the higher eigenmodes, which correspond to subordinate deformations, can be thought of as "local expansions".
Here, let us quantify the spatial correlations of the eigenvectors for each mode, in order to characterize the collective and local expansions. The mean polarization −→ m for a vector field − → o m i is defined by ( Bottinelli et al., 2016 ) Fluctuations around the mean polarization are measured with a polarization correlation function, defined by where d i, j is the initial distance between the centers of the triangles i and j . Color plots of the correlation function C m ( d ) are displayed in Fig. 4 B, F, and J. It is clear that the first eigenmode, which corresponds to the collective expansion, has a longer correlation length than the higher modes in all cases. For the uniform growth, only four eigenmodes contribute to the deformation, and the correlations for the other eigenmodes are extremely small (see Fig. 3. Eigenvalues λ x m and λ y m as a function of the mode m . The eigenvalues for case (i) are higher than those for cases (ii) and (iii), indicating that the velocities in case (i) fluctuate more than those in cases (ii) and (iii). The eigenvalues for the uniform growth are relatively small, reflecting the fact that only a few eigenmodes are dominant. also the eigenvalues smaller than double precision for m > 5 in Fig. 3 ). For case (i), the spatial correlation monotonically decays as a function of the mode m , indicating that the eigenmodes gradually become localized as m increases. In contrast, for case (ii) the spatial correlation does not decay monotonically. For example, as shown in Fig. 4 J, the spatial correlation of eigenmode 4 is longer than those of eigenmodes 2 and 3. We call such a spatially correlated mode with m ≥ 2 an excited mode. This excited mode should relate to an abnormal shape deformation because the excited mode corresponds to a "locally collective" growth mode.
To confirm the emergence of these excited modes for different random seeds, we define a spatial correlation length D ( m ) as the minimum distance d that satisfies C m (d) = 0 (i.e., transition points from red to white in Fig. 4 B, F, and J). Note that D ( m ) is a correlation length for the orientations of the eigenvectors − → e m i . A numerical estimate of D ( m ), obtained by averaging over 30 different random seeds, is presented in Fig. 5 . Furthermore, the SD for D ( m ) is displayed in Fig. 5 with error bars. As seen from these error bars, the correlation length D ( m ) fluctuates more widely for the system with less temporal randomness in the elasticity (i.e., the system with smaller q ), for which an abnormal shape deformation occurs. The large fluctuations in D ( m ) imply that the locally collective growth mode can be observed for different random seeds and the excited modes cannot be attributed to the random seed, but to the fixed parameters. For the uniform growth ( Fig. 4 C), highly fluctuating triangles are located uniformly around the origin, and for case (i) ( Fig. 4 G) highly fluctuating triangles are located uniformly along the boundary ∂ 2 . In contrast, for case (ii) ( Fig. 4 K) highly fluctuating triangles are located non-uniformly, and are concentrated at topcentral, left, and right regions. In these regions, several circles form clusters, which we call "spots". Thus, in Fig. 4 K there are three spots, which are mechanically weak. We emphasize that this non-uniformity for highly fluctuating triangles and the formation of spot structures seem to be consistent with the non-uniform growth of the plate, i.e., the abnormal shape deformation shown in Fig. 4 I. Namely, in Fig. 4 K the three spot structures formed by highly fluctuating triangles correspond to the three bumps observed in the abnormal shape shown in Fig. 4 I. In this sense, we refer to these spot structures as soft spots. It should be noted that the typical spatial scale of the soft spots ranged from 0.3 to 0.6, but their parameter dependence is beyond the scope of this study because the scale of the soft spot is variable even with fixed parameters.

"Soft spots": magnitudes of eigenvectors
To summarize, when we focus on the temporal evolution of the spatial randomness the results can be interpreted as follows. A system with weak spatial randomness may be considered to have almost uniform elasticity, and the highly fluctuating triangles appear uniformly on the plate, thereby not resulting in the emergence of any soft spots. On the other hand, for a system with strong spatial randomness the highly fluctuating triangles appear non-uniformly, and therefore soft spots emerge. Accordingly, the soft spots facilitate the formation of abnormal organ shapes, as shown in Fig. 2 . To suppress the emergence of soft spots, the nonuniformity in the elasticity should be effectively reduced by shuffling the elasticity using temporal noise.

Discussion
In this study, it was observed that the spatial and temporal randomness play opposing roles in shape deformation processes. As is evident from the phase diagram ( Fig. 2 ), the temporal randomness is essential for obtaining a precise organ shape, whereas the spatial randomness disturbs the reproducibility of the final shape. We expected that an underlying mechanism may be hidden in the trajectories of the triangles, and analyzed the velocity correlation matrix along with its eigenvalues and eigenvectors. As a result, we found that an abnormal shape deformation can be created by two independent factors: (1) excited modes where the orientations of the eigenvectors with m ≥ 2 are spatially correlated, and/or (2) soft spots where spot structures with high fluctuations in the magnitudes of the eigenvectors appear non-uniformly.
These findings may shed light on the biological mystery concerning macroscopic uniformity and microscopic non-uniformity that was described in the introduction. There are two possible scenarios for achieving macroscopic uniformity (i.e., reproducibility) with microscopic non-uniformity (i.e., randomness in elasticity). (1) For the case of a weak microscopic spatial non-uniformity, the highly fluctuating triangles are distributed almost uniformly in space, resulting in no emergence of soft spots, as exemplified by case (iii). In this case, macroscopic uniformity is easily achieved, owing to the weakness of the spatial randomness, and the temporal randomness is not so important for the reproducibility of the final shape.
(2) For the case of strong microscopic spatial nonuniformity, soft spots and excited modes are created, as illustrated by case (ii). However, through temporal noise the spatial randomness can be effectively averaged out, resulting in no emergence of soft spots or excited modes, as exemplified by case (i). As a result, the reproducibility of the final shape is recovered, and thus in this case the temporal randomness plays a significant role in the precise shape deformation.
We could detect how such macroscopic uniformity properties (absence of soft spots and/or excited modes) can emerge from microscopic non-uniformity through spatio-temporal growth noise. In addition to the concept of the spatio-temporal averaging of growth noise proposed in the previous study ( Hong et al., 2016 ), a supporting methodology of the "suppression of soft spots and excited modes" is proposed, and should provide a basic bridge to resolving the apparent contradiction between microscopic non-uniformity and macroscopic uniformity.

Funding
This work was supported by MEXT KAKENHI Grant-in-Aid for Scientific Research on Innovative Areas "Plant-Structure Optimization Strategy" Grant Numbers JP18H05484 (to Taku Demura), and by the Special Postdoctoral Researcher Program in RIKEN (Award number K1731017 ).