Wrinkle patterns in active viscoelastic thin sheets

We show that a viscoelastic thin sheet driven out of equilibrium by active structural remodeling, such as during fast growth, 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 remodeling are unable to relax to a conﬁguration that minimizes the elastic energy and the sheet is inherently out of equilibrium. We argue that this nonequilibrium 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 conﬁgurations that minimize 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.


I. INTRODUCTION
Thompson set the mathematical foundation for describing and classifying the astonishing diversity of shapes and form in the living world [1]. A century later, our understanding of biological processes at the molecular level has been vastly improved [2], yet it is still largely unknown how the formation of large, functional structures such as tissues and organs arises from these molecular processes [3]. 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 [3]. The formation of large structures such as tissues and organs is a result of a complex set of guided collective mechanochemical 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 at equilibrium.
Out-of-equilibrium biological processes are naturally described within the framework of active matter physics, where the system is driven out of equilibrium by a constant input of energy at the microscopic scale [4]. Despite great progress in understanding the behavior of active fluids, much less is known about how activity affects the behavior of solid and viscoelastic materials, such as tissues [5][6][7]. 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 [8]. Such dynamical heterogeneity is a hallmark of an active glassy state [9], with epithelial cell monolayers being prime examples of such behavior [10][11][12][13]. The biological significance of dynamical heterogeneity is only starting to emerge. When it comes to describing bending deformations in active systems, only recently has a theoretical description been proposed [14,15] and applied to the systems with axisymmetric geometry [16]. Similarly, new in vitro experiments on active buckling in poroelastic contractile actin sheets [17] and kinesin-driven twodimensional (2D) microtubule sheets [18] have been reported.
In this paper, we study thin elastic and viscoelastic sheets with activity introduced as a dynamical change of the reference shape, with the simplest example of such process being growth. Recently, it has been argued that the imbalance of surface and volume growth and the rate of relaxation can account for a morphological diversity of prebiotic organisms [19]. Physically, activity provides structural remodeling that acts as a local time-dependent source of strain. The timedependent reference shape can be either stress free (embeddable metric) or contain residual stress (nonembeddable metric) [20]. While the distinction between the two cases has important consequences for the elastic ground state [20,21], 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 remodeling, depending on whether or not the particular reference state is embeddable in R 3 . It has been recently argued [22] that viscoelastic relaxation can stabilize cell shapes during morphogenesis. Such viscoelastic effects remove all stresses over a sufficiently long time. Here, we focus on the regime where active remodeling 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.  1. (a) Insertion of an elastic disk into an aperture smaller than the disk's size (top) induces residual stress in the disk due to compression. This residual stress can be released by buckling out of plane (bottom). (b) In a discrete picture where an elastic sheet is represented as a triangulation of a surface, the geometric incompatibility leading to the residual stress (i.e., nonembeddable metric) can be understood as two triangles that share an edge having mutually incompatible preferred shapes, e.g., the red (lower) and blue (upper) triangles are the preferred shapes for the corresponding gray mesh triangle shown underneath it.

A. Thin active viscoelastic sheet
We study a thin sheet of size L and uniform thickness h L with linear elastic response [23]. We assume that the surrounding fluid provides damping, but ignore all other hydrodynamic effects. The sheet is represented by the twodimensional midsurface, initially in the xy plane. The deformed midsurface with no overhangs can be parametrized as r = r(x, y) = [x, y, w(x, y)], where w(x, y) is a sufficiently smooth height function. One defines the metric, g αβ = ∂ α r · ∂ β r, and curvature, c β α = g βγ b αγ , tensors, where b αβ = −∂ α r · ∂ β n (α, β ∈ {x, y}) is the second fundamental form and n = (∂ x r × ∂ y r)/|∂ x r × ∂ y r| is the unit normal vector [24] [ Fig. 2(a)]. The elastic energy of the midsurface is [20,25] where u αβ = 1 2 (g αβ − g αβ ) is the strain tensor, g αβ is a reference metric tensor, dA = √ det gdxdy is the area element, A αβγ δ 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 R 3 space, while Greek indices are used to label intrinsic curvilinear coordinates. For an isotropic material, where Y is the Young's modulus, ν is the Poisson ratio, and g αγ g γ β = δ α β . 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 and with u β α = g βγ u αγ and c β α = g βγ b αγ [20]. With the mean curvature H = 1 2 c α α ≡ 1 2 Tr(ĉ) and the Gaussian curvature The sheet is represented as a two-dimensional midsurface, parametrized by (x, y) coordinates, with two tangent vectors (e 1 ≡ ∂ x r, e 2 ≡ ∂ y r) and a unit-length normal [n = (e 1 × e 2 )/|e 1 × e 2 |] assigned to each point of the surface. For numerical implementation, the surface is discretized in terms of triangles. (b) Red (outer) vectors form the reference metric tensor,ḡ αβ =ā α ·ā β , and blue (inner) vectors form the realized metric tensor, g αβ = a α · a β . The strain tensor is defined as u αβ = 1 2 (g αβ − g αβ ). (c) Viscoelasticity is modeled as a relaxation of the reference metric towards the realised metric, with a characteristic timescale τ v . K = det (c β α ), the bending energy given by Eq. (4) becomes where κ = h 3 Y/12(1 − ν 2 ) is the bending stiffness (B). The timescale of bending elastic deformations is τ el ∼ 1/κ. 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 remodeling does not affect elastic parameters. In reality, material properties are affected by the structural remodeling. However, imposing spatial and time dependence on the elastic parameters did not qualitatively change our findings and, for simplicity, in the following we assume them to be constant. Active remodeling is introduced by imposing dynamical changes of the reference metric. The precise functional form of active remodeling is not important, as long as one can associate a typical timescale τ a to it. Active remodeling can be thought of as a generalization of growth, with the quasistatic differential growth being described asḡ αβ (r, t ) = a(r)tḡ αβ (r, t = 0), where a(r) > 0 and τ growth a ≡ a −1 τ el . The advantage of expressing deformation with respect to the reference metric [26] is that the formalism can be directly generalized to include active remodeling and viscoelastic relaxation, without making only assumptions about the existence of a stress-free reference state.
Finally, the metric of the sheet is subject to viscous relaxation. This process acts to dissipate the energy introduced by activity. We model viscoelastic dissipation via internal rearrangement processes [27] leading to the relaxation of the reference metric towards the realized (i.e., current) metric, where τ v is the viscous remodeling timescale [ Fig. 2(c)].

B. Active sheet dynamics
We assume overdamped dynamics and solve the set of firstorder equations for each vertex i and discrete metric of each triangle, Here, r i ∈ R 3 is the position vector of vertex i, and η i (t ) ∈ R 3 is a weak random noise, obeying η i = 0 and η m with m, n ∈ {x, y, z}. γ is the friction coefficient modeling 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 of 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 γ δ αβ (t ) is a tensor that sets the rate of viscous relaxation. While, in general, V γ δ αβ is a function of time, here we assume it to be constant, V γ δ αβ = 1 τ v δ γ α δ δ β . We note that a similar metric-based approach to modeling viscoelastic behavior has recently been discussed in Refs. [28,29]. R αβ is a tensor function that prescribes the active remodeling rate. Here, R αβ models metric expansion that explicitly depends on time and, thus, describes dynamical changes of the active remodeling rate. Finally, discrete versions of the realized and reference metric tensors are defined in Fig. 2(b). Equations (7a) and (7b) are integrated numerically using the standard first-order Euler-Maruyama discretization scheme keeping connectivity of the triangulation fixed (see the numerical implementation in Appendix E). Expressions for ∇ r i E and R αβ are given in Appendices C and D, respectively.
A description based on the time-evolving reference metric is also suitable for direct discretization (Fig. 2) and efficient parallel implementation on GPUs. This allows us to simulate systems containing up to 2 × 10 6 triangles, removing the need to implement complex remeshing procedures to avoid reduction in accuracy in the vicinity of high-curvature folds. Note that in the current implementation, we do not include steric effects and the sheet can take unphysical self-intersecting configurations. The inclusion of self-avoidance is possible, but technically challenging to efficiently implement on GPUs [30]. The steric effect would indeed affect the folding patterns, but would not change our main conclusions. Values of the parameters used in simulations are given in Table I, Appendix E. Moreover, length is measured in units of h, time in units of t * = γ /Y h, and energy in units of κ.

A. Timescales
As noted above, there are three relevant timescales in this problem: the elastic relaxation time τ el , the viscous dissipation timescale τ v , and the active relaxation timescale τ a . Here we estimate these timescales for a typical biological system such as an early embryo.
We start by estimating τ el . We assume that a flat sheet is suspended in a fluid that only provides drag (i.e., the fluid acts as a simple sink for the sheet's momentum). Furthermore, considering only out-of-plane motion, the overdamped dynamics of the sheet in the continuum limit reads [31] ∂ where w(x, y) measures the vertical displacement from the flat configuration, 2 is the bi-Laplacian operator, and is the friction per unit area. The magnitude of the drag force on a disk of radius R moving perpendicular to its plane with velocity v in a fluid of viscosity η is F d = 16ηRv [32]. This gives = 16 π η R . Noting that for scaling purposes ∂ t ∼ 1/τ el and ∼ 1/R 2 [31], we estimate For an epithelial cell sheet in water, τ el ∼ 10 1 -10 2 s, consistent with [33]. Clearly, the timescale 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 [34], as well as cell growth and division. Processes related to the cytoskeleton typically occur at timescales of τ a ∼ 10 1 -10 2 s [35,36], while cell growth and division are slower and can span several hours [2]. Dissipation in tissues results from multicellular rearrangements (i.e., plastic events such as intercalations, ingressions, and extrusions) and subcellular cytoskeleton remodeling (i.e., cell-shape relaxation). We note that dissipation is accompanied by entropy production and, in general, an entropy production equation would be required [14]. Here, we are not concerned with the details of the dissipative processes (rendering the entropy production equation unnecessary) and assume that they occur on a timescale τ v .
We note, however, that cell rearrangements are typically slower (occurring on the scale ∼10 min) than the subcellular remodeling (occurring on the seconds to minutes scale). While it is not always the case, the out-of-equilibrium situation with τ a < τ el , τ v is, therefore, biologically plausible and, we argue, beneficial to access the diversity of shapes needed to form complex structures. In the following, we explore the range of possible dynamical shape patterns formed in the nonequilibrium regime.

B. Nonequilibrium wrinkling patterns
Here we explore out-of-equilibrium dynamics of flat disks of radius R subject to active remodeling and viscous dissipation (Fig. 3). The choice of the disk geometry is inspired by extensive work on wrinkling patterns due to tension [37,38] or resulting from a quasiequilibrium growth, e.g., during biofilm formation [39][40][41]. We note, however, that this regime corresponds to τ el τ a . We assume that a ring of radius r i < R is kept fixed, but can transmit stress. Active remodeling is assumed to occur only in the outer annulus, for r i < r < R. With no viscoelastic relaxation and slow active remodeling (lower-left corner in Fig. 3), the system is in the extensively 013165-3

FIG. 3.
A snapshot of the out-of equilibrium shapes obtained by numerical integration of Eqs. (7a) and (7b) starting from a flat disk configuration. The snapshots are taken at t = 10 4 t * . The vertical axis represents the rate of viscous (dissipative) relaxation with increasing values designating faster residual stress relaxation. On the horizontal axis, we plot the active structural remodeling rate, with larger values corresponding to faster changes of the local reference metric. The usual slow, quasiequilibrium elastic growth corresponds to the lower-left corner in this graph. Colors represent the height function, w(x, y). In these simulations, R αβ is time independent. Different geometries are shown in Appendix A. studied quasiequilibrium 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 remodeling 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 remodeling is dissipated by a fast relaxation of the reference metric of the sheet. As one increases the remodeling 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 remodeling is in part dissipated, which prevents wrinkles from growing rapidly (upper-right region in Fig. 3). As shown in Fig. 4, the ratio between active relaxation and viscous dissipation then determines the steady-state wrinkling patterns. If viscoelastic relaxation is introduced, the stress generated by active remodeling is in part dissipated, which prevents wrinkles from growing rapidly. The ratio between active relaxation and viscous dissipation determines the wrinkling pattern morphology. These patterns, however, do not correspond to the minima of elastic energy and thus exhibit far richer morphologies compared to the equilibrium states. Note that in our model, the system would typically not reach a steady state. Therefore, one would need to introduce an additional mechanism that suppresses active remodeling in order to stabilize the system. Biological systems indeed have such homeostatic mechanisms.
Furthermore, if the system is able to dynamically tune the active remodeling rate, it can reach conformations that would otherwise require overcoming large energy barriers. For example, for a fixed high value of τ −1 a , one needs to inject substantial energy in order to initiate wrinkling (Fig. 5, circles). On the other hand, if the initial value of τ −1 a is reduced, the wrinkling energy barrier is significantly lowered (Fig. 5,  triangles). This is not surprising as elastic relaxation is not fast enough to accommodate structural changes due to fast active remodeling. If τ −1 a is increased once the wrinkles are formed, however, it is easy to reach different wrinkling patterns (Fig. 5, pentagons) without the high initial energy cost.
If we note that wrinkles are predominantly in the radial direction, in order to analyze the effects of dynamical tuning 013165-4 FIG. 5. Energy per vertex as a function of simulation time. To reach a wrinkled configuration with a remodeling rate τ −1 a = 5 × 10 −5 would require the sheet to overcome a large energy barrier (peak of the orange (circles) curve). If the initial remodeling rate, however, is set to τ −1 a = 2 × 10 −5 , the system requires less energy to reach the wrinkling instability (peak of the blue (triangles) curve). Upon switching to τ −1 a = 5 × 10 −5 at t = 2 × 10 3 t * , the evolution continues along the green (pentagons) curve and the system reaches a wrinkling pattern, which is very similar to the one obtained by following the orange (circles) curve, as shown by the two snapshots on the right.
of the system's parameters on wrinkle patterns, we compute the power spectrum of the Monge parametrization w(r, θ ) represented in polar coordinates of an annular ring w R (θ ) = 1/N R r∈(R,R+ R) w(r, θ ), where N R is the number of mesh vertices in the annulus with the inner radius R and the outer radius R + R. The power spectrum S ww ( f θ ) = |w R ( f θ )| 2 , where w R ( f θ ) = θ w R (θ )e i f θ θ is the Fourier transform of w R (θ ). The position of the peaks of S ww ( f θ ) corresponds to the position of the dominant wave number f θ . The power spectra for the three systems discussed in Fig. 5 are shown in Fig. 6. The left column shows S ww ( f θ ) at t = 2 × 10 3 t * at the point when the system marked with green symbols (pentagons) had its parameters switched from those of the system marked with blue symbols (triangles) to those of the system marked with orange symbols (circles). As expected, the "blue" (triangles) and "green" (pentagons) power spectra are nearly identical along with the entire system. On the other hand, the system with a higher active remodeling rate (orange/circles) has distinctly different features for intermediate radii. At t = 4 × 10 3 t * , "orange" (circles) and "green" (pentagons) systems have the same energy (Fig. 5), yet their power spectra are significantly different. This is not surprising as the "green" (pentagons) system retains a memory of its past, which is an expected feature in a system far from equilibrium and also suggests that the total energy is not a good quantity to describe the actual wrinkling pattern. Finally, at t = 8 × 10 3 t * , the power spectra of the "orange" (circles) and "green" (pentagons) systems are similar to each other and distinct from the "blue" (triangles) system. The only exception is Fig. 6(c), where the "green" (pentagons) system  (30,40), and (d) (40,50). The × marks peak values of each power spectrum. The symbols and colors have the same meaning as in Fig. 5. All power spectra were computed using Welch's method [43]. still retains some similarity to the "blue" (triangles) system. This is, however, not surprising since in this region, wrinkles are irregular and would take the longest to reshape.
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, e.g., via adenosine triphosphate-dependent cellular processes [42].

IV. SUMMARY AND CONCLUSIONS
By applying an active solid model to viscoelastic thin sheets subject to active structural remodeling, 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 remodeling allowing local perturbations to grow. As a consequence, the shape patterns depend on 013165-5 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 [44]. Most theoretical approaches are based on continuum mechanics augmented to encode the effects of growth into Föppl-von Kármán equations [45][46][47]. The salient point in such treatments is that elastic relaxation occurs at the timescales that are short compared to growth and thus describe the regime where the system is always in quasistatic mechanical equilibrium [45,48]. We argue that the out-of-equilibrium regime studied here is of particular interest in developing a physical understanding of morphogenesis, where the embryo has to undergo a series of carefully orchestrated shape changes to develop the functioning organism. In birds, for example, this process also has to be energy efficient given the finite amount of nutrients available in the egg.
Our observations are in line with and extend the recent findings of Ruiz-Herrero et al. [19]. 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 remodeling 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 remodeling, relaxation, and mechanical parameters. This would be much easier to encode in the space available in the genome.

APPENDIX A: ALTERNATIVE SHEET GEOMETRIES AND MODELS OF ACTIVE REMODELING
Here we show two examples of geometries and structural remodeling. The first example (Fig. 7) shows the case of a strip of size L x = 100 and L y = 20 under uniform structural and viscous remodeling. As in the case of Fig. 3, if one increases the rate of viscoelastic relaxation while keeping the active remodeling, wrinkles are less pronounced or, in the case of very fast dissipative relaxation, do not form at all.
The second example (Fig. 8) mimics active compression of a flat disk of size R = 50. The compression is introduced by imposing a rapid strain through an instantaneous change on the reference metric in an external annulus 0.8R < r < R. Viscous remodeling is assumed to occur only in the inner annulus, for r < 0.8R.

APPENDIX B: BRIEF SUMMARY OF THIN SHEET ELASTICITY
Starting from the expression for the elastic energy of a thin three-dimensional solid of size L and uniform thickness h L with linear elastic response, we derive the expression for energy of its two-dimensional neutral surface [23]. We consider a growing body with homogenous and isotropic elastic properties, which configuration is described by a strain tensor u αβ = 1 2 (g αβ − g αβ ), where both the reference metric g αβ and the realized metric g αβ are assumed to be continuous and nonsingular [20,49]. Using Einstein summation convention, the elastic energy can be written as where dV = √ det gdxdydz is the volume element, and the elastic tensor is A αβγ δ = λg αβ g γ δ + μ(g αγ g βδ + g αδ g βγ ),  , w(x, y).

013165-6
where λ and μ are two elastic constants related to the Young's modulus E and Poisson's ratio ν via . (B3)

Two-dimensional plate energy density
Expression of the elastic energy density of the neutral surface can be derived under the Kirchhoff-Love assumptions [20,50]: (1) Body is in the state of plane stress, i.e., stress normal to the surfaces parallel to the neutral surface can be neglected.
(2) Points which lie on a normal to the neutral surface in the reference configuration remain on the same normal in the deformed configuration.
Using these assumptions, it follows that The elastic energy density in Eq. (B1) can now be rewritten as where the two-dimensional elastic tensor is Using the Kirchhoff-Love assumptions effectively decouples different sheets parallel to the neutral surface [25]. Therefore, one can obtain the expression for the total elastic energy of the neutral surface by integrating along the sheet's thickness (chosen to be the z direction), In the small strain approximation, we can neglect all terms that are cubic or higher power in u αβ to obtain where b i j = e i · ∂ j n is the second fundamental form, related to the curvature tensor c j i = g ik b k j [51]. Equation (B8) is the expression for the elastic energy of a thin shell expressed in terms of its neutral surface. The first term in the twodimensional 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. Therefore, one can write E 2D = E s + E b [20], where and Y = Eh. Similarly,

Stretching energy
Following Ref. [52], coordinates of a given point P inside a triangle can be written in terms of the two vectors that span the triangle, r 1 and r 2 , as where 0 ξ 1 and 0 η 1 are coordinates of vector r in the basis {r 1 , r 2 }. For convenience, we introduce a third vector Then, r P =ĥs P with s P = (ξ η ω ≡ 0) T . Further, the squared distance between P and any point Q [r Q =ĥs Q , and s Q = (φ ψ 0) T ] in the triangle is whereĝ =ĥ Tĥ is the (discrete) metric of the reference triangle, which can be easily computed in a simulation. Consider an affine deformation of the triangle that contains point P defined as then the displacement vector (P → P ) can be written as whereÎ is the identity matrix. In order to derive the expression for the nonlinear strain tensor in term of matricesĥ andĥ, we recall that in the mixed tensor form, where we have dropped the P subscript for convenience. We then find Inserting Eq. (C7) into Eq. (C6) leads tô We can now use Eq. (B9) to write where A T is the triangle area and Finally, the force associated with E s is then given by {(F 11 + νF 22 ) 2α 11 r 12 + α 12 r 13 α 12 r 12 + (F 22 + νF 11 ) α 12 r 13 2α 22 r 13 + α 12 r 12 + (1 − ν)[F 21 α 11 r 13 2α 12 r 13 + α 11 r 12 + F 12 2α 12 r 12 + α 22 r 13 α 22 r 12 ]} + e s 4A T g 22 r 12 − g 12 r 13 g 11 r 13 − g 12 r 12 −1 1 0 where α μν = (ĝ −1 ) μν and e s = E s /A T .

Bending energy
We start from Eq. (B10) and use κ G = −κ (1 − ν) to obtain From the definition of the second fundamental, it immediately follows that Following Ref. [53], the 1 2 κ∂ γ n · ∂ γ n term is a continuum version of the expression where T i is the triangle i, T j are three of its neighboring triangles,κ is the discrete value of the bending rigidity, and the subscript "SN" stands for "Seung-Nelson." The sum in Eq. (C15) can be written as where b is the distance between the centers of two neighboring triangles, and vectors v are v ( j) = cos ( 2π 3 j)e x + sin ( 2π 3 j)e y (Fig. 9). For b → 0, we expand n(r i + bv ( j) ) to the linear FIG. 9. A triangle T i and three of its nearest neighbors T j . In the continuum limit, b → 0. The shaded area is the vertex area element. order in b, Thus, We can now calculate the j-sum explicitly for each component φ = x, y and ψ = x, y, and we assume that ∂ φ n is calculated at point r i . We write . Thus, the total discrete energy is where H T i and K T i are the mean and Gaussian curvature of the triangle T i , and the sum goes over all triangles. From Fig. 9, we see that the area element A i = b 2 √ 3 2 , which leads to Comparing Eq. (C12) with the last expression, we obtain κ = 1 √ 3κ . Note that a different prefactor was obtained in Refs. [53] and [54], since in Eq. (C17) we truncated the expansion too early. The exact value of the constant of proportionality is, however, of the same order of magnitude, κ = √ 3 2κ [53]. In practice, this difference is inconsequential since we useκ as a tuning parameter. Next, the bending force matrix f p,e is [(r 02 · r 03 )r 02 −(r 02 · r 02 )r 03 (r 02 · r 03 )r 01 + (r 01 · r 02 )r 03 − 2(r 01 · r 03 )r 02 (r 01 · r 02 )r 02 − (r 02 · r 02 )r 01 ] + (A k · A l ) 1 |A k | 2 ((r 02 · r 02 )r 01 − (r 01 · r 02 )r 02 (r 01 · r 01 )r 02 − (r 01 · r 02 )r 01 0) + 1 |A l | 2 (0(r 03 · r 03 )r 02 − (r 02 · r 03 )r 03 (r 02 · r 02 )r 03 − (r 02 · r 03 )r 02 ) , with r i j defined in Fig. 10.

APPENDIX D: EXPRESSION FOR ACTIVE REMODELING
Remodeling is introduced as a change in the local reference metricḡ. Here we choose a circular geometry for which we have the following natural metric (Fig. 11 with β 11 and β 22 being the remodeling rates in the R 12 and R 13 direction of the reference configuration, and φ is the angle between these two vectors. Equations (D1) and (D2) can be easily discretized. FIG. 11. Active remodeling is introduced as a change of the reference metric of each triangle.

APPENDIX E: NUMERICAL IMPLEMENTATION AND SIMULATION PARAMETERS
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. We use the Brownian dynamics approach [55]. All computation-heavy tasks are fully implemented on the GPU, so that there are no transfers between device and 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 [56] as an external visualization software for testing and presentation purposes.
The coarse-grained triangular meshes used in the simulation were created using a public domain package GMSH [57], 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 a 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 I. It is important to note that all material parameters are assumed to be time independent and uniform across the entire mesh.
The active remodeling processes are assumed not to be uniform on the mesh. In particular, we have chosen to restrict active and viscous remodeling to an external annulus of 20 < r < 50. The remodeling and viscous remodeling rates 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, where μ is the inverse friction coefficient and F i is the total force acting on the vertex i due to the mesh deformation, and F R is a uniform random force whose magnitude fulfills the fluctuation-dissipation theorem for the given inverse friction coefficient and temperature, T ; in our simulation, we set μ = 1.0 and T = 10 −6 . In addition, the integration is set to be 10 −3 for remodeling rates equal to or smaller than 10 −3 and 10 −5 otherwise.