Wrinkle patterns in active viscoelastic thin sheets

TL;DR

This paper investigates how active viscoelastic thin sheets develop diverse wrinkle patterns due to the interplay of viscous relaxation and active remodeling, especially when activity outpaces relaxation, leading to non-equilibrium morphologies relevant in biological systems.

Contribution

It introduces a formalism of evolving target metric to model active remodeling in viscoelastic sheets, revealing new out-of-equilibrium wrinkle morphologies.

Findings

Diverse wrinkle patterns emerge from active viscoelastic sheets.

Out-of-equilibrium dynamics enable access to unique morphologies.

Active processes dominate relaxation, preventing energy-minimizing configurations.

Abstract

We show that a viscoelastic thin sheet driven out of equilibrium by active structural remodelling develops a rich variety of shapes as a result of a competition between viscous relaxation and activity. In the regime where active processes are faster than viscoelastic relaxation, wrinkles that are formed due to remodelling are unable to relax to a configuration that minimises the elastic energy and the sheet is inherently out of equilibrium. We argue that this non-equilibrium regime is of particular interest in biology as it allows the system to access morphologies that are unavailable if restricted to the adiabatic evolution between configurations that minimise the elastic energy alone. Here, we introduce activity using the formalism of evolving target metric and showcase the diversity of wrinkling morphologies arising from out of equilibrium dynamics.

Tables

Table 1. Table 1: Simulation parameters.

Streaching Potential	Value
Young’s modulus, $E$	${10}^2$
Plate thickness, $h$	${10}^0$
Poissson’s ratio, $\nu$	$1/3$
Seung-Nelson Bending Potential
Bending modulus, $\kappa$	$5\times {10}^{-2}\left(E{h}^3\right)$

Equations

$$E=\int dA{\mathcal{A}}^{\alpha \beta \gamma \delta }\left(\frac{h}{2}{u}_{\alpha \beta }{u}_{\gamma \delta }+\frac{h^3}{24}{b}_{\alpha \beta }{b}_{\gamma \delta }\right),$$

$$\gamma {\dot{\mathbf{r}}}_i$$

$${\dot{\bar{g}}}_{\alpha \beta }\left(\mathbf{r},t\right)$$

$${\mathbf{r}}^{\prime }=\mathbf{r}+\mathbf{u},$$

$$u_{ij}=\frac{1}{2}\left(g_{ij}-{\bar{g}}_{ij}\right).$$

$${\bar{g}}_{ij}={\partial }_i\mathbf{r}\cdot {\partial }_j\mathbf{r}={\mathbf{e}}_i\cdot {\mathbf{e}}_j,$$

$$\begin{array}{rl}{g}_{ij}& ={\bar{g}}_{ij}+{\mathbf{e}}_i\cdot {\partial }_j\mathbf{u}+{\mathbf{e}}_j\cdot {\partial }_i\mathbf{u}+{\partial }_i\mathbf{u}\cdot {\partial }_j\mathbf{u},\end{array}$$

$$u_{ij}=\frac{1}{2}\left({\nabla }_j{u}_i+{\nabla }_i{u}_j+{\nabla }_i{u}_s{\nabla }_j{u}^s\right),$$

$$u_{ij}=\frac{1}{2}\left({\partial }_i{u}_j+{\partial }_j{u}_i+{\partial }_i{u}_s{\partial }_j{u}^s\right).$$

$$\begin{array}{rl}{E}_{el}& \approx E\left({\bar{g}}_{ij}\right)+{\frac{\partial E}{\partial g_{ij}}|}_{u_{ij}=0}{u}_{ij}+\frac{1}{2}{\frac{{\partial }^{2}E}{\partial g_{ij}\partial g}|}_{u_{ij}=0}{u}_{ij}{u}_{kl}\\ & +o\left(u^3\right)\\ & =E\left({\bar{g}}_{ij}\right)+\frac{1}{2}{A}^{ijkl}{u}_{ij}{u}_{kl}+o\left(u^3\right),\end{array}$$

$$E_{el}=\frac{1}{2}{A}^{ijkl}{u}_{ij}{u}_{kl}.$$

$$E_{tot}=\frac{1}{2}{\int }_V\sqrt{\left|g\right|}{A}^{ijkl}{u}_{ij}{u}_{kl}.$$

$$A^{ijkl}=\lambda g^{ij}{g}^{kl}+\mu \left(g^{ik}{g}^{jl}+g^{il}{g}^{jk}\right),$$

$$\begin{array}{rl}E& =\frac{\mu \left(3\lambda +2\mu \right)}{\lambda +\mu }\\ \nu & =\frac{\lambda }{2\left(\lambda +\mu \right)}.\end{array}$$

$${\sigma }^{i3}=0i=x,y,z$$

$$g_{ij}=\left(\begin{array}[]{cc}g_{\alpha\beta}&0\\ 0&1\end{array}\right)\,\,\,\,\mathrm{or\,\,\,\,\epsilon_{\alpha 3}=0.}$$

$$u_3^3=u_{33}=-\frac{\lambda }{\lambda +2\mu }{u}_{\alpha }^{\alpha }.$$

$$\begin{array}{rl}{E}_{tot}^{2D}& =\frac{1}{2}{\mathcal{A}}^{\alpha \beta \gamma \delta }{u}_{\alpha \beta }{u}_{\gamma \delta },\end{array}$$

$${\mathcal{A}}^{\alpha \beta \gamma \delta }=2\mu \left(\frac{\lambda }{\lambda +2\mu }{g}^{\alpha \beta }{g}^{\gamma \delta }+g^{\alpha \gamma }{g}^{\beta \delta }\right).$$

$$\begin{array}{rl}{E}_{tot}^{2D}& =\frac{1}{2}{\int }_S{\int }_{-\frac{h}{2}}^{\frac{h}{2}}dz\sqrt{\left|g\left(z\right)\right|}{\mathcal{A}}^{\alpha \beta \gamma \delta }{u}_{\alpha \beta }\left(z\right)u_{\gamma \delta }\left(z\right).\end{array}$$

$$E_{tot}^{2D}={\int }_S\sqrt{\left|g\right|}{\mathcal{A}}^{\alpha \beta \gamma \delta }\left(\frac{h}{2}{u}_{\alpha \beta }{u}_{\gamma \delta }+\frac{h^3}{24}{b}_{\alpha \beta }{b}_{\gamma \delta }\right),$$

$$E_{stretch}=\frac{Y}{2\left(1+\nu \right)}\left(\frac{\nu }{1-\nu }{u}_{\alpha }^{\alpha }{u}_{\beta }^{\beta }+u_{\beta }^{\alpha }{u}_{\alpha }^{\beta }\right),$$

$$\begin{array}{rl}{E}_{bend}& =\frac{h^3}{24}\frac{E}{\left(1+\nu \right)}\left(\frac{\nu }{1-\nu }{b}_{\alpha }^{\alpha }{b}_{\gamma }^{\gamma }+b_{\alpha }^{\beta }{b}_{\beta }^{\alpha }\right)\\ & =2\kappa H^2+{\kappa }_{G}K,\end{array}$$

$$\mathbf{r}=\xi \mathbf{a}+\eta \mathbf{b},$$

$$\begin{array}{rl}\mathbf{a}& ={\mathbf{r}}_2-{\mathbf{r}}_1,\\ \mathbf{b}& ={\mathbf{r}}_3-{\mathbf{r}}_1,\end{array}$$

$$\mathbf{c}=\mathbf{a}\times \mathbf{b}.$$

$$\hat{h}=\left(\begin{array}[]{cc}\mathbf{a}&\mathbf{b}\end{array}\right)\equiv\left(\begin{array}[]{ccc}a_{x}&b_{x}&c_{x}\\ a_{y}&b_{y}&c_{y}\\ a_{z}&b_{z}&c_{z}\end{array}\right),$$

$$\mathbf{r}=\hat{h}\mathbf{s},$$

$$\mathbf{s}=\left(\begin{array}[]{ccc}\xi&\eta&\omega\end{array}\right)^{T},$$

$$r_i=h_i^{\alpha }{s}_{\alpha },$$

$${\mathbf{r}}_Q=\hat{h}{\mathbf{s}}_Q.$$

Full text

D’Arcy Thompson set the mathematical foundation for describing and classifying the astonishing diversity of shapes and form in the living world (Thompson, 1942). A century later, our understanding of biological processes at the molecular level has been vastly improved (Alberts et al., 2014), yet it is still largely unknown how the formation of large, functional structures such as tissues and organs arises from these molecular processes (Wolpert et al., 2015). A unifying feature of all higher organisms is that they start as a single cell, a zygote, and autonomously develop into an individual, without external input. The genome provides a template that steers development towards the desired body plan (Wolpert et al., 2015). The formation of large structures such as tissues and organs is a result of a complex set of guided collective mechano-chemical processes. To select a specific morphology, the phase space of possible shapes has to be large. Furthermore, transition between shapes should be possible at a reasonably low cost, which is hard to achieve in equilibrium.

Out of equilibrium biological processes are naturally described within the framework of the active matter physics, where the system is driven out of equilibrium by a constant input of energy at the microscopic scale (Marchetti et al., 2013). Despite great progress in understanding the behaviour of active fluids, much less is known about how activity affects the behaviour of solid and viscoelastic materials, such as tissues (Harris et al., 2012; Berthoumieux et al., 2014; Matoz-Fernandez et al., 2017). Numerical simulations of dense self-propelled elastic disks, for example, showed that part of the energy intake is diverted into local elastic deformations leading to prominent spatial and temporal heterogeneities in observed velocity fields (Henkes et al., 2011). Such dynamical heterogeneity is a hallmark of an active glassy state (Berthier and Kurchan, 2013), with epithelial cell monolayers being prime examples of such behaviour (Petitjean et al., 2010; Angelini et al., 2011; Chepizhko et al., 2018; Henkes et al., 2019). The biological significance of dynamical heterogeneity is only starting to emerge. When it comes to describing bending deformations in active systems, only recently a theoretical description has been proposed (Salbreux and Jülicher, 2017).

In this paper, we study thin elastic and viscoelastic sheets with activity introduced as a dynamical change of the reference shape. Physically, activity provides structural remodelling that acts as a local time-dependent source of strain. The time-dependent reference shape can be either stress free (embeddable metric) or contain residual stress (non-embeddable metric) (Efrati et al., 2009). While the distinction between the two cases has important consequences for the elastic ground state (Efrati et al., 2009; Kang et al., 2014), it is not essential for the present discussion. As shown in Fig. 1, bending out of plane can fully or partly remove the residual stresses due to remodelling, depending on whether the particular reference state is embeddable or not in $\mathbb{R}^{3}$. It has been recently argued (Clément et al., 2017) that viscoelastic relaxation can stabilise cell shapes during morphogenesis. Such viscoelastic effects remove all stresses over a sufficiently long time. Here, we focus on the regime where active remodelling is faster than both elastic and viscoelastic relaxation, leading to the system being inherently out of equilibrium. This regime is expected to be of particular importance to early embryonic development.

We study a thin sheet of size $L$ and uniform thickness $h\ll L$ with linear elastic response (Audoly and Pomeau, 2010). We assume that the surrounding fluid provides damping but ignore all other hydrodynamic effects. The sheet is represented by the two-dimensional mid-surface, initially in the $xy$ plane. The deformed mid-surface with no overhangs can be parametrised as $\mathbf{r}=\mathbf{r}\left(x,y\right)=\left(x,y,w\left(x,y\right)\right)$, where $w\left(x,y\right)$ is a sufficiently smooth height function. One defines the metric, $g_{\alpha\beta}=\partial_{\alpha}\mathbf{r}\cdot\partial_{\beta}\mathbf{r}$, and curvature, $c_{\alpha}^{\beta}=g^{\beta\gamma}b_{\alpha\gamma}$, tensors where $b_{\alpha\beta}=-\partial_{\alpha}\mathbf{r}\cdot\partial_{\beta}\mathbf{n}$ ($\alpha,\beta\in\left\{x,y\right\}$) is the second fundamental form and $\mathbf{n}=\left(\partial_{x}\mathbf{r}\times\partial_{y}\mathbf{r}\right)/\left|\partial_{x}\mathbf{r}\times\partial_{y}\mathbf{r}\right|$ is the unit normal vector (Do Carmo, 1976) (Fig. 2a). The elastic energy of the mid-surface is (Koiter, 1966; Efrati et al., 2009; SI, )

[TABLE]

where $u_{\alpha\beta}=\frac{1}{2}\left(g_{\alpha\beta}-\overline{g}_{\alpha\beta}\right)$ is the strain tensor, $\overline{g}_{\alpha\beta}$ is a reference metric tensor, $dA=\sqrt{\det g}dxdy$ is the area element, $\mathcal{A}^{\alpha\beta\gamma\delta}$ is the elastic tensor, and summation over pairs of repeated indices is assumed. Latin indices refer to the components of vectors in the embedding Euclidean $\mathbb{R}^{3}$ space, while Greek indices are used to label intrinsic curvilinear coordinates. For an isotropic material, $\mathcal{A}^{\alpha\beta\gamma\delta}=\frac{Y}{1+\nu}\left(\frac{\nu}{1-\nu}\overline{g}^{\alpha\beta}\overline{g}^{\gamma\delta}+\overline{g}^{\alpha\gamma}\overline{g}^{\beta\delta}\right)$, where $Y$ is the Young’s modulus and $\nu$ is the Poisson ratio and $\overline{g}^{\alpha\gamma}\overline{g}_{\gamma\beta}=\delta_{\beta}^{\alpha}$. The first term in Eq. (1) is the stretching energy and the second term accounts for bending. For an isotropic material, stretching and bending energies simplify to $E_{s}=\frac{h}{2}\int dA\frac{Y}{1+\nu}\left(\frac{\nu}{1-\nu}u_{\alpha}^{\alpha}u_{\beta}^{\beta}+u_{\alpha}^{\beta}u_{\beta}^{\alpha}\right)$ and $E_{b}=\frac{h^{3}}{24}\int dA\frac{Y}{1+\nu}\left(\frac{\nu}{1-\nu}c_{\alpha}^{\alpha}c_{\beta}^{\beta}+c_{\alpha}^{\beta}c_{\beta}^{\alpha}\right)$, with $u_{\alpha}^{\beta}=\overline{g}^{\beta\gamma}u_{\alpha\gamma}$ and $c_{\alpha}^{\beta}=\overline{g}^{\beta\gamma}b_{\alpha\gamma}$ (Efrati et al., 2009; SI, ). With the mean curvature $H=\frac{1}{2}c_{\alpha}^{\alpha}\equiv\frac{1}{2}\text{Tr}\left(\hat{c}\right)$ and the Gaussian curvature $K=\det\left(c_{\alpha}^{\beta}\right)$, the bending energy becomes $E_{b}=\int dA\kappa\left(2H^{2}-\left(1-\nu\right)K\right)$, where $\kappa=h^{3}Y/12\left(1-\nu^{2}\right)$ is the bending stiffness. In general, material properties and the reference metric can be position dependent and the sheet can have a spontaneous curvature, $H_{0}$. Here we assume that $H_{0}=0$ and the active remodelling does not affect elastic parameters. In reality, material properties are affected by the structural remodelling. However, imposing spatial and time dependence on the elastic parameters did not qualitatively change our findings and, for simplicity, in following we assume them be constant. Finally, we estimate that the relaxation time associated with bending, $\tau_{el}\sim\eta L^{3}/\kappa$, here $\eta$ is the dynamical viscosity of the surrounding fluid (SI, ). For an epithelial cell sheet in water, $\tau_{el}\sim 10^{1}-10^{2}$ s, consistent with (Marmottant et al., 2009). Clearly, the time scale of relaxation associated with stretching deformation is much shorter and consequently of no importance for the present discussion.

Active effects in a tissue result, for example, from myosin driven contractions and turnover of the actin cytoskeleton (Joanny and Prost, 2009) as well as cell growth and division. Processes related to the cytoskeleton typically occur at time scales of $\tau_{a}\sim 10^{1}-10^{2}$ s (Rauzi et al., 2008, 2010), while cell growth and division are slower and can span several hours (Alberts et al., 2014). Dissipation in tissues results from multi-cellular rearrangements (i.e., plastic events such as intercalations, ingressions and extrusions) and sub-cellular cytoskeleton remodelling (i.e., cell shape relaxation). We note that dissipation is accompanied by entropy production and, in general, an entropy production equation would be required (Salbreux and Jülicher, 2017). Here, we are not concerned by the details of the dissipative processes (rendering the entropy production equation unnecessary) and assume that they occur on a time scale, $\tau_{v}$. We note, however, that cell rearrangements are typically slower (occurring on the scale $\sim 10$ min) than the sub-cellular remodelling (occurring on the seconds to minutes scale). While it is not always the case, the out of equilibrium situation with $\tau_{a}<\tau_{el},\tau_{v}$ is, therefore, biologically plausible and, we argue, beneficial to access diversity of shapes needed to form complex structures. In the following, we explore the range of possible dynamical shape patterns formed in the non-equilibrium regime.

The advantage of expressing deformation with respect to the reference metric (Sknepnek and Olvera de la Cruz, 2012) is that the formalism can be directly generalised to include active remodelling and viscoelastic relaxation, without making only assumptions about the existence of a stress free reference state. Here, active remodelling is introduced by imposing dynamical changes of the reference metric. The precise functional form of active remodelling is not important, as long as one can associate a typical time scale, $\tau_{a}$, to it. Active remodelling can be thought of as a generalisation of growth, with the quasi-static differential growth being described as $\bar{g}_{\alpha\beta}\left(\mathbf{r},t\right)=a\left(\mathbf{r}\right)t\bar{g}_{\alpha\beta}\left(\mathbf{r},t=0\right)$, where $a\left(\mathbf{r}\right)>0$ and $\tau_{a}^{growth}\equiv a^{-1}\ll\tau_{el}$. We model viscoelastic effects as a relaxation of the reference metric towards the realised metric (Fig. 2c). Therefore, viscoelastic relaxation has the opposite effect of elasticity, for which the reference state conforms to the realised shape rather than the other way around. A description based on the time-evolving reference metric is also suitable for direct discretisation (Fig. 2) and efficient parallel implementation on GPUs (SI, ). This allows us to simulate systems containing up to $2\times 10^{6}$ triangles removing the need to implement complex remeshing procedures to avoid reduction in accuracy in the vicinity of high-curvature folds.

We assume overdamped dynamics and solve the set of first-order equations for each vertex $i$ and discrete metric of each triangle,

[TABLE]

Here $\mathbf{r}_{i}\in\mathbb{R}^{3}$ is the position vector of vertex $i$, $\boldsymbol{\eta}_{i}\left(t\right)\in\mathbb{R}^{3}$ is a weak random noise, obeying $\left\langle\boldsymbol{\eta}_{i}\right\rangle=0$ and $\left\langle\eta_{i}^{m}\left(t\right)\eta_{j}^{n}\left(t^{\prime}\right)\right\rangle=\sqrt{2\gamma k_{B}T}\delta_{ij}\delta_{mn}\delta\left(t-t^{\prime}\right)$ with $m,n\in\left\{x,y,z\right\}$. $\gamma$ is the friction coefficient modelling dissipation by the surrounding fluid and $T$ is the temperature kept very low and used only for numerical convenience to avoid being trapped in shallow local minima. All our simulations were effectively at $T=0$ as thermal fluctuations are not expected to play an appreciable role in biological systems, i.e., relevant energy scales far exceed $k_{B}T$. V_{\alpha\beta}^{\gamma\delta}$$\left(t\right) is a tensor that sets the rate of viscous relaxation. While in general $V_{\alpha\beta}^{\gamma\delta}$ is a function of time, here we assume it to be constant, $V_{\alpha\beta}^{\gamma\delta}=\frac{1}{\tau_{v}}\delta_{\alpha}^{\gamma}\delta_{\beta}^{\delta}$. $R_{\alpha\beta}$ is a tensor function that prescribes active remodelling rate. Here, $R_{\alpha\beta}$ models metric expansion and is given in Eq. (S34) in (SI, ). Furthermore, $R_{\alpha\beta}$ explicitly depends on time and, thus, models dynamical changes of the active remodelling rate. Finally, discrete versions of the realised and reference metric tensors are defined in Fig. 2b. Eqs. (2a) and (2b) are integrated numerically using standard first-order Euler-Maruyama discretisation scheme keeping connectivity of the triangulation fixed. Expressions for the gradient of energy in Eq. (2a) are straightforward but lengthy (SI, ). Note that in the current implementation, we do not include steric effects and the sheet can take unphysical self-intersecting configurations. Including self-avoidance is possible but technically challenging to efficiently implement on GPUs. Steric effect would indeed affect the folding patterns but would not change our main conclusions. Values of parameters used in simulations are given in (SI, ). Moreover, length is measured in units of $h$, time in units of $t^{*}=\gamma/Yh$ and energy in units of $\kappa$.

We explored out of equilibrium dynamics of flat disks of radius $R$ subject to active remodelling and viscous dissipation (Fig. 3). The choice of the disk geometry is inspired by extensive work on wrinkling patterns due to tension (Jagla, 2007; Davidovitch et al., 2011) or resulting from a quasi-equilibrium growth, e.g., during biofilm formation (Amar and Wu, 2014; Yan et al., 2019). This regime corresponds to $\tau_{el}\ll\tau_{a}$. We assume that a ring of radius $r_{i}<R$ is kept fixed but can transmit stress. Active remodelling is assumed to occur only in the outer annulus, for $r_{i}<r<R$. With no viscoelastic relaxation and slow active remodelling (lower left corner in Fig. 3), the system is in the extensively studied quasi-equilibrium differential growth regime. Free expansion of the outer boundary can relieve part of the stress produced by growth. There is, however, no such stress relief mechanism in the tangential direction and the sheet forms a regular pattern of radial wrinkles. The inner disk, on the other hand, is compressed in both directions leading to wrinkles with no preferred orientation. If one instead allows for viscoelastic relaxation while keeping the active remodelling slow (left column in Fig. 3), wrinkles are less pronounced or, in the case of very fast dissipative relaxation, do not form at all (top left in Fig. 3). This is easy to understand, as in this regime the stress generated by active remodelling is dissipated by a fast relaxation of the reference metric of the sheet. As one increases the remodelling rate (second and third columns in Fig. 3), wrinkling patterns become more pronounced and less regular, especially close to the inner ring, where stress accumulation is strong. Without viscous dissipation (bottom right in Fig. 3) the sheet continues to expand and quickly reaches unphysical self-intersecting configurations. In a real system, steric repulsion and intrinsic biological processes such as apoptosis due to hypoxia and nutrient deprivation would prevent this uncontrolled growth. If viscoelastic relaxation is introduced, the stress generated by active remodelling is in part dissipated, which prevents wrinkles from growing rapidly (upper right region in Fig. 3). The ratio between active relaxation and viscous dissipation then determines the steady state wrinkling patterns. These patterns, however, do not correspond to minima of elastic energy and thus exhibit far richer morphologies compared to the equilibrium states (Fig. S1 in (SI, )).

Furthermore, if the system is able to dynamically tune the active remodelling rate, it can reach conformations that would otherwise require ovecoming large energy barriers. For example, for a fixed high value of $\tau_{a}^{-1}$, one needs to inject substantial energy in order to initiate wrinkling (Fig. S1, circles). On the other hand, if the initial value of $\tau_{a}^{-1}$ is reduced, the wrinkling energy barrier is significantly lowered (Fig. S1, triangles). This is not surprising as elastic relaxation is not fast enough to accommodate structural changes due to fast active remodelling. If $\tau_{a}^{-1}$ is increased once the wrinkles are formed, however, it is easy to reach different wrinkling patterns (Fig. S1 pentagons) without the high initial energy cost. This simple example shows that an out of equilibrium system is not only able to develop a rich variety of morphologies but it also can avoid costly energy barriers between different patterns by dynamically tuning its parameters, which most biological systems are equipped to do.

By applying an active solid model to viscoelastic thin sheets subject to active structural remodelling, we showed that the interplay between activity and viscous relaxation leads to a diverse morphology of out of equilibrium wrinkling patterns. Of particular interest in this study is the regime where active processes are faster than elastic and viscoelastic relaxation. In this case, the system has no time to fully relax local stresses produced by active remodelling allowing local perturbations to grow. As a consequence, the shape patterns depend on the initial conditions and local fluctuations. This is in stark contrast to the mechanics of growth, in particular in plants, that has been extensively studied with great success (Goriely, 2017). Most theoretical approaches are based on continuum mechanics augmented to encode the effects of growth into Föppl-von Kármán equations (Amar and Goriely, 2005; Goriely and Amar, 2005, 2007). The salient point in such treatments is that elastic relaxation occurs at the time scales that are short compared to growth and thus describe the regime where the system is always in quasi-static mechanical equilibrium (Rodriguez et al., 1994; Amar and Goriely, 2005). We argue that the out of equilibrium regime studied here is of particular interest in developing physical understanding of morphogenesis.

We note that a similar observation has been recently made in a study of the dynamics of growth and form in prebiotic vesicles (Ruiz-Herrero et al., 2019) where the observed diversity of shapes was associated with the imbalance of surface and volume growth and the rate of relaxation. This suggests that keeping a growing system out of equilibrium significantly increases the range of available morphologies. The development of higher organisms is too complex to be captured by a simple mechanical model of actively remodelling sheets. Our observations, however, point to a mechanism by which a system that is kept out of equilibrium could be steered towards a desired shape by a careful regulation of remodelling, relaxation and mechanical parameters. This would be much easier to encode in the space available in the genome.

RS would like to thank C. J. Weijer for his valuable insights into developmental biology. FD, NSW and DMF were funded by the UK BBSRC (Award BB/P001335/1). RS acknowledges support by the UK BBSRC (Award BB/N009789/1).

I Energy

Here we show the elastic energy for different value of $\tau_{a}^{-1}$ and $\tau_{v}^{-1}$. Fig. S1 shows that the elastic energy is far from the global minimum with prominent regions of highly concentrated bending energy. If viscoelastic relaxation is introduced, the stress generated by active remodelling is in part dissipated, which prevents wrinkles from growing rapidly (upper right region in Fig. 3 main text). The ratio between active relaxation and viscous dissipation then determines the steady state wrinkling patterns. These patterns, however, do not correspond to minima of elastic energy and thus exhibit far richer morphologies compared to the equilibrium states. It is also easy to transition between different wrinkling patterns by tuning system parameters, which most biological systems are equipped to do. Note that while the precise morphology of wrinkling patterns depends on the geometry of the system, the mechanism that leads to such out of equilibrium structures does not.

II Other geometries and remodelling tensors

Here we show two examples of geometries and structural remodelling. The first example (Fig. S2) shows the case of a strip of size $L_{x}=100$ and $L_{y}=20$ under uniform structural and viscous remodelling. As in the case of Fig. 3 in the main text, if one increase the rate of viscoelastic relaxation while keeping the active remodelling wrinkles are less pronounced or, in the case of very fast dissipative relaxation, do not form at all.

The second example (Fig. S3), mimics active compression of a flat disk of size $R=50$. The compression is introduced by imposing a rapid strain trough an instantaneous change on the reference metric in an external annulus $0.8R<r<R$. Viscous remodelling is assumed to occur only in the inner annulus, for $r<0.8R$.

III Elasticity

Starting for the energy expression for a thin three-dimensional solid we derive energy expression for its two-dimensional neutral surface. The neutral surface is placed midway along the thin direction (see Fig. S4) and on it bending and stretching are decoupled.

III.1 Strain tensor

A point at position $\mathbf{r}$ under the deformation is displaced to (Green, )

[TABLE]

where vector $\mathbf{u}\left(\mathbf{r}\right)$ is the displacement vector. The elastic energy of the body cannot depend on actual displacement but it depends on the derivatives of $\mathbf{u}$. In other words, the elastic energy depends on the changes in the metric of the body. We define the strain tensor

[TABLE]

If we recall that

[TABLE]

and

[TABLE]

the strain tensor becomes

[TABLE]

where $\nabla_{j}$ is the covariant derivative. In the Euclidean space, $\nabla_{i}=\partial_{i}$, and the last expression reduces to the familiar definition of the strain tensor defined in standard text books on elasticity

[TABLE]

III.2 Three-dimensional elastic energy density

The elastic energy density depends on the strain tensor, i.e., on the metric, $E_{el}=E_{el}\left(g_{ij}\right)$. If the strain is small we can expand $E_{el}$ in powers of $u_{ij}$ around the reference configuration ($\bar{g}_{ij}$ and $u_{ij}=0$)

[TABLE]

where we assumed the linear term in the expansion vanishes and introduced a contravariant elastic tensor $A^{ijkl}=\left.\frac{\partial^{2}E}{\partial g_{ij}\partial g_{kl}}\right|_{u_{ij}=0}$. We can omit the unimportant constant term $E\left(\bar{g}_{ij}\right)$ to obtain the expression for elastic energy density in the small strain approximation

[TABLE]

We need to make a distinction between small strain and small displacement approximations, that is, one does not imply the other. If the strains are small and the elastic response of the body is directly proportional to the applied stress, the small strain approximation is applicable and Eq. (S5) is valid, i.e. constitutive laws are linear (Hookean elasticity). However, when studying thin shells, we can have a situation that although strains are small displacements are large. In this case Eq. (S5) is still valid, but we cannot omit nonlinear terms in $u_{ij}$ and we have to use Eq. (S3).

Total elastic energy is Koiter ; Efrati ,

[TABLE]

For an isotropic body there the elastic tensor has only two independent components and can be written as (e.g., Ref. (Efrati, ))

[TABLE]

where $\lambda$ and $\mu$ are two Lame coefficients. We can introduce Young’s modulus $E$ and Poisson’s ratio $\nu$ via

[TABLE]

III.3 Two-dimensional plate Energy density

Expression of the elastic energy density of the neutral surface can be derived (Koiter2, ) under the Kirchhoff-Love assumptions (Refs. (Koiter2, ) and (Efrati, )):

Body is in the state of plane-stress, i.e., stress normal to the surfaces parallel to the neutral surface can be neglected. 2. 2.

Points which lie on a normal to the neutral surface in the reference configuration remain on the same normal in the deformed configuration.

These assumptions translate into

[TABLE]

where $\sigma^{ij}$ is the contravariant stress tensor, and

[TABLE]

From Eqs. (S7) and (S8) we have

[TABLE]

The elastic energy density in Eq. (S5) can be now rewritten as

[TABLE]

where the two-dimensional elastic tensor is

[TABLE]

Using the Kirchhoff-Love assumptions effectively decouples different sheets parallel to the neutral surface (Koiter, ). Therefore, we can obtain the expression for the total elastic energy of the neutral surface by integrating along the sheet thickness (chosen to be the $z$ direction),

[TABLE]

In the small strain approximation, we can neglect all terms that are are cubic or higher power in $u_{\alpha\beta}$ to get

[TABLE]

where $b_{ij}=\mathbf{e_{i}}\cdot\partial_{j}\mathbf{n}$ is the second fundamental form, related to the curvature tensor $c_{i}^{j}=g^{ik}b_{kj}$ (doCarmo, ). Eq. (S11) is the expression for the elastic energy of a thin shell expressed in terms of its neutral surface. The first term in the two-dimensional energy expression is stretching energy and it describes energy penalty of stretching or compressing of the neutral surface. The second term is the bending energy, which describes energy penalty of flexing the sheet. We can write $E^{2D}=E_{stretch}+E_{bend}$ where

[TABLE]

and $Y=Eh$. Similarly, using $\mathrm{Tr}\left(b_{\alpha}^{\beta}\right)=2H$ and $\mathrm{Tr}\left(\left(b_{\mu}^{\gamma}\right)^{2}\right)=4H^{2}-2K$,

[TABLE]

with the bending modulus, $\kappa=Eh^{3}/12\left(1-\nu^{2}\right)$ and Gaussian modulus $\kappa_{G}=-Eh^{3}/12\left(1+\nu\right)$.

IV Discrete model

We discretise the surface using a triangular mesh (Seung, ), Fig. S5.

IV.1 Stretching energy

We start with the stretching energy term. Closely following Ref. (Parrinello, ), coordinates of a given point $P$ inside a triangle can be written in terms of the two vectors $\mathbf{a}$ and $\mathbf{b}$ (Fig. S5),

[TABLE]

where $0\leq\xi\leq 1$ and $0\leq\eta\leq 1$ are coordinates of vector $\mathbf{r}$ in the basis $\left\{\mathbf{a},\mathbf{b}\right\}$. Note that vectors $\mathbf{a}$ and $\mathbf{b}$ are themselves three dimensional vectors constructed as

[TABLE]

where $\mathbf{r}_{1}$, $\mathbf{r}_{2}$ and $\mathbf{r}_{3}$ are positions of the three corners of the triangle. For convenience, that is to be able to work with square matrices that are invertible, we will introduce a third vector

[TABLE]

We can now construct a $3\times 3$ matrix

[TABLE]

such that

[TABLE]

where

[TABLE]

where $\omega\equiv 0$. In term of coordinates Eq. (S16) can be written as

[TABLE]

and we have used Latin indices to count components in the embedding space and Greek indices to count components of vector $\mathbf{s}$. If we now introduce point $Q$ in the same triangle with coordinates $\phi$ and $\psi$, i.e., with \mathbf{s}_{Q}=\left(\begin{array}[]{ccc}\phi&\psi&0\end{array}\right)^{T} then

[TABLE]

Square of the distance $l^{2}$ between points $P$ and $Q$ is

[TABLE]

and we have used parentheses in the sub- and superscript to designate that $P$ and $Q$ are not indices and have freely changed the name of the repeated summation indices. Matrix

[TABLE]

is the (discrete) metric of our reference triangle. Explicitly,

[TABLE]

i.e., matrix $\hat{\bar{g}}$ is a $3\times 3$ matrix with a $2\times 2$ sub-matrix corresponding to the metric tensor of the triangle and the $\mathbf{c}\cdot\mathbf{c}$ term that is added for convenience. Matrix $\hat{\bar{g}}$ can be easily computed.

After the deformation edges of the triangle change and the basis vectors become $\mathbf{a}^{\prime}$ and $\mathbf{b}^{\prime}$, while the point $P$ has moved to the new position $P^{\prime}$ with coordinates

[TABLE]

where \hat{H}=\left(\begin{array}[]{ccc}\mathbf{a}^{\prime}&\mathbf{b}^{\prime}&\mathbf{c}^{\prime}\end{array}\right)^{T}\left(\begin{array}[]{ccc}\mathbf{a}^{\prime}&\mathbf{b}^{\prime}&\mathbf{c}^{\prime}\end{array}\right). Note that we assume that the deformation is linear and as such automatically affine. An important of this restrictions is that the point $P^{\prime}$ has the same coordinates $\left(\xi,\eta,0\right)$ in the deformed triangle $\left\{\mathbf{a}^{\prime},\mathbf{b}^{\prime}\right\}$ as the point $P$ had in the original undeformed state. Point $P^{\prime}$ is point $P$ after deformation. If we recall Eq. (S1) and the definition of the displacement vector

[TABLE]

where in the second line we have inverted Eq. (S16) to get $\mathbf{r}=\hat{h}^{-1}\mathbf{s}$. Finally,

[TABLE]

where $\hat{I}$ is the identity matrix. In terms of coordinates we have

[TABLE]

In order to derive the expression for the non-linear strain tensor in term of matrices $\hat{H}$ and $\hat{h}$ we use Eq. (S4) repeated for convenience (in a form of a mixed tensor),

[TABLE]

Now we compute

[TABLE]

If we plug the last expression into Eq. (S19) we obtain

[TABLE]

where we have used $\left(AB\right)^{T}=B^{T}A^{T}$ and $A^{-T}\equiv\left(A^{-1}\right)^{T}$. We can now use Eq. (S12) to write

[TABLE]

where $A_{T}$ is the triangle area and the tensor $\hat{F}$ is given by,

[TABLE]

IV.1.1 Vertex stretching force

From Eqs. (S17) and (S21) is easy to see that it is more convenient to express stretching energy in terms of the distance between vertices of one triangle, $\mathbf{r}_{ij}$ ($\mathbf{r}_{12},\mathbf{r}_{13}$) (Fig. S6). Thus, using the chain rule we can write the force over the triangle vertices as,

[TABLE]

After lengthy but straightforward algebra, the expression for the force on vertex $i$ due to stretching reads,

[TABLE]

where $\alpha_{ij}=\left(\bar{g}^{-1}\right)_{ij}$. This expression can be directly implemented in a simulation.

IV.2 Bending energy

In order to derive expressions of the discrete version of the bending energy we start from Eq. (S13) and use the fact that $\kappa_{G}=-\kappa\left(1-\nu\right)$ to obtain

[TABLE]

We can now use, $\mathrm{Tr}\left[\left(b_{\mu}^{\gamma}\right)^{2}\right]=4H^{2}-2K$ to write

[TABLE]

If we use the definition of the second fundamental form we obtain

[TABLE]

Following Ref. (Seung, ) the $\frac{1}{2}\kappa\partial_{\gamma}\mathbf{n}\cdot\partial^{\gamma}\mathbf{n}$ term is a continuum version of the expression

[TABLE]

where $T_{i}$ is the triangle $i$, $T_{j}$ are three of its neighbours and $\tilde{\kappa}$ is the discrete value of the bending rigidity and subscript “SN” stands for “Seung-Nelson”. $\tilde{\kappa}$ is related proportional to $\kappa$ and we’ll discuss the constant of proportionality below. The sum in Eq. (S26) can be written as

[TABLE]

where $b$ is the distance between centers of two neighbouring triangles and vectors $\mathbf{v}$ are (see Fig. S7)

[TABLE]

For $b\to 0$ we can use Taylor series to expand $\mathbf{n}\left(\mathbf{r}_{i}+b\mathbf{v}^{\left(j\right)}\right)$ to the linear order in $b$,

[TABLE]

Thus,

[TABLE]

We can now calculate the $j$-sum explicitly for each component $\phi=x,y$ and $\psi=x,y$,

[TABLE]

and

[TABLE]

and we assume that $\partial_{\phi}\mathbf{n}$ is calculated at point $\mathbf{r}_{i}$. We write $E_{SN}=-\frac{3}{2}\frac{b^{2}}{2}\tilde{\kappa}\partial^{\phi}\mathbf{n}b_{\phi}^{\mu}\mathbf{e}_{\mu}=\frac{3}{2}\frac{b^{2}}{2}\tilde{\kappa}\mathrm{Tr}\left(b_{\phi}^{\mu}\right)^{2}=\frac{3b^{2}}{4}\tilde{\kappa}\left(4H^{2}-2K\right)=\frac{3b^{2}}{2}\tilde{\kappa}\left(2H^{2}-K\right)$. Thus, the total discrete energy is

[TABLE]

where $H_{T_{i}}$ and $K_{T_{i}}$ are mean and Gaussian curvature of the triangle $T_{i}$ and the sum goes over all triangles. From Fig. S7 we see that the area element $A_{i}=\frac{b^{2}\sqrt{3}}{2}$, which leads to

[TABLE]

which in the limit $A_{i}\to 0$ becomes an integral

[TABLE]

Comparing Eq. (S25) with the last expression we obtain $\kappa=\frac{1}{\sqrt{3}}\tilde{\kappa}$. Note that this is a different prefactor than obtained in Refs. (Seung, ) and (Schmidt, ). The reason is that in Eq. (S27) we have truncated the expansion to early. The exact constant of proportionality is of the same order of magnitude and is given as $\kappa=\frac{\sqrt{3}}{2}\tilde{\kappa}$ Therefore, we have showed that Eq. (S26) is a good discrete approximation for the continuum elastic energy. The advantage of the last expression is that it can be easily computed in a simulation.

IV.2.1 Vertex bending force

The bending energy of a triangular mesh can be expressed as Eq. (S26),

[TABLE]

where $e$ is an edge of the mesh and $\mathbf{n}_{e,1(2)}$ are the unit-length normals of two neighbouring triangles. For the shake of simplicity we consider a sole edge and two triangles $k$ and $l$ as is shown in Fig. S8. Note that $\mathbf{n}_{k(l)}=\mathbf{A}_{k(l)}/A_{k(l)}$, where $\mathbf{A}_{k(l)}$ is a vector normal to the triangle having length equal to the triangle’s area.

Now the bending force over a vertex $p$ is the negative gradient $E_{B}$ calculated at $\mathbf{r}_{p}$,

[TABLE]

is easy to see that since we are dealing with a sum over all the edges then it suffice to calculate the gradient $\nabla_{\mathbf{r}_{p}}$ for one generic edge,

[TABLE]

After straightforward but lengthy algebra, the force matrix $\mathbf{f}_{p,e}$, Eq. (S30) is,

[TABLE]

with

[TABLE]

V Active Remodelling

Remodelling is introduced as a change in the local reference metric $\bar{g}$. Here we choose a circular geometry for which we have the following natural metric (Fig. S9),

[TABLE]

and for simplicity we impose a non-shear linear uniform remodelling, i.e.,

[TABLE]

with $\beta_{11}$ and $\beta_{22}$ being the remodelling rates in the $\mathbf{R}_{12}$ and $\mathbf{R}_{13}$ direction and $\phi$ is the angle between $\mathbf{R}_{12}$ and $\mathbf{R}_{13}$. Eqs. (S33) and (S34) can be easily discretised by expressing the metric tensor in the triangle laboratory coordinates.

VI Viscoelastic relaxation

As in Ref. Muoz we model viscoelastic dissipation via an internal rearrangement processes leading to the relaxation of the reference metric towards the realised (i.e., current) metric, i.e.,

[TABLE]

where ${g}_{ij}$, $\bar{g}_{ij}$ are the current and reference metric tensors at time $t$, respectively, and $\tau_{ve}$ is the remodelling time scale. The physical interpretation of Eq. (S35) is that the energy is dissipated in local rearrangement processes.

VII Numerical implementation

We have built our own parallel GPU-based (NVidia CUDA) implementation of the discrete model outlined in the previous sections. Our code is specifically designed to introduce different sources of activity into the system. The general workflow is shown in Fig. S10. All computation-heavy task are fully implemented on the GPU, so that there are no transfers between DEVICE-HOST during the execution. The only routines executed by the host are those required by the user in order to save data.

Our CUDA kernels are moderately optimized, trying to keep aligned and coalesced memory access avoiding threads divergence and only using atomic functions when absolutely necessary. Finally, we used ParaView paraview as an external visualisation software for testing and presentation purposes.

VII.1 Simulation parameters

The coarse-grained triangular meshes used in the simulation were created using a public domain package Gmsh (gmesh, ) setting the edge target length to $l=0.35$ and the plate radius equal to $R=50$, with all lengths measured in units of thickness, $h$. In order to obtain different initial configuration the vertices are moved randomly in $(x,y)$ around the initial configuration using a normal distribution with standard deviation equal to $10^{-3}$. After this procedure the Device Mesh is created and the reference metric is set to the mesh actual metric.

The potentials used in our simulation with its respective parameters are listed in Table 1. It is important to note that all material parameter are assumed to be time-independent and uniform across the entire mesh.

The active remodelling processes are assumed not to be uniform on the mesh. In particular, we have chosen to restrict remodelling and remodelling to an external annulus of $20<r<50$. The remodelling and viscous remodelling rate are set to be uniform inside of the annulus, for the respective values used in the simulation, see Fig. 3.

To integrate the vertex equation of motion, we have implemented a Brownian dynamics integrator,

[TABLE]

where $\mu$ is the inverse friction coefficient and $\mathbf{F}_{i}$ is the total force acting on the vertex $i$ due the mesh deformation and $\mathbf{F}_{R}$ is a uniform random force whose magnitude fullfil the fluctuation-dissipation theorem for the given inverse friction coefficient and temperature, $T$; in our simulation we set $\mu=1.0$ and $T=10^{-6}$. In addition, the integration is set to be $10^{-3}$ for remodelling rates equal or smaller than $10^{-3}$ and $10^{-5}$ otherwise.

VIII Elastic relaxation time

Here we make a rough estimate of the elastic relaxation time scale assuming that the a nearly flat sheet is suspended in a fluid. We assume that the fluid only provides drag and do not consider any effects of its flow, i.e., the fluid acts as a simple sink for the sheet’s momentum. Assuming only out of plane motion described by the high function $w(x,y)$, in the overdamped limit the equation of motion is,

[TABLE]

where $\Delta^{2}$ is the bi-laplacian operator and $\Gamma$ is the friction coefficient due to fluid. Note that $\Gamma$ has dimensions of $\frac{\text{mass}}{\text{length}^{2}\times\text{time}}$, and is thus interpreted as the friction per unit area. If we recall the well-known result in fluid dynamics (Lamb, ) that the drag coefficient on a disk of radius $R$ moving perpendicular to its plane in a fluid of viscoelastic $\eta$ is $\zeta=16\eta Rv$, where $v$ is the velocity, we obtain $\Gamma=\frac{16}{\pi}\frac{\eta v}{R}$. Therefore, we estimate

[TABLE]

or

[TABLE]

For an epithelial tissue of size $R\sim 1$ mm in water, assuming bending rigidity $\kappa\sim 10^{-12}$ J, we estimate $\tau_{el}\sim 10$ s.