Joseph Ackermann

Paul-Qiuyang Qu

Loïc Legoff
email: loic.le-goff@univ-amu.fr

Martine Ben Amar

Contributed Loïc Legoff

Modeling the mechanics of growing epithelia with a bilayer plate theory

Epithelia, which consists of cell sheets lying on a substrate, are prevalent structures of multicellular organisms. The physical basis of epithelial morphogenesis has been intensely investigated in recent years. However, as 2D mechanics focused most attention, we still lack a rigorous description of how the mechanical interactions between the cell layer and its substrate can lead to 3D distortions. This work provides a complete description of epithelial mechanics using the most straightforward model of an epithelium: a thin elastic bilayer. We first provide experimental evidence in Drosophila tissues that localized alterations of the cell substrate (the extracellular matrix) can lead to profound 3D shape changes in epithelia. We then develop an analytical model modifying the Föppl-von Kármán equation with growth for bilayers. We provide a complete description of all contributions from biophysical characteristics of epithelia. We show how any localized inhomogeneity of stiffness or thickness drastically changes the bending process when the two layers grow differently. Comparison with finite element simulations and experiments performed on Drosophila wing imaginal discs validates this approach for thin epithelia

Introduction

Biological cell assemblies are often organized in laminar structures called epithelia. Epithelia cover most of our hollow organs, such as the esophagus [START_REF] Li | Growth and surface folding of esophageal mucosa: a biomechanical model[END_REF], intestine [START_REF] Ben Amar | Anisotropic growth shapes intestinal tissues during embryogenesis[END_REF][START_REF] Hannezo | Instabilities of monolayered epithelia: shape and structure of villi and crypts[END_REF], and stomach. The skin is also an epithelium constituted of several cell layers. The structure of epithelia is of one or several attached cell layers, which rest on a more disorganized polymeric substrate, called the extra-cellular matrix (ECM). Even the simplest epithelium, constituted of a single cell layer, must rest on an abutting ECM layer.

In the past decade, a lot of emphasis has been put on the role of the cytoskeletal cortex on the apical side of the cells (the side away from the ECM layer) in setting mechanical properties of epithelia [START_REF] Pinheiro | Mechanical force-driven adherens junction remodeling and epithelial dynamics[END_REF]. However, recent investigations have established that the mechanical properties of the ECM may be as important in setting the shape of a tissue [START_REF] Khalilgharibi | To form and function: on the role of basement membrane mechanics in tissue development, homeostasis and disease[END_REF]. Both the cell layer and the ECM must be taken into account in the mechanical description of epithelia.

From the biomechanical viewpoint, an epithelium may be viewed as a bilayer made of two different soft materials with different stiffnesses and different ways to grow: the cell layer grows by mass increase and proliferation, while the ECM grows by the addition of polymer chains and swelling. A bilayer can also serve as a suitable approximation for more complex geometries, such as the multilayered skin, subdivided into the dermis and the epidermis. In between layers of cells or in between cells and connective tissues, the ECM may develop into a stiff membrane enriched in collagen filaments called the basal membrane [START_REF] Ciarletta | Papillary networks in the dermal-epidermal junction of skin: a biomechanical model[END_REF]. Epithelial cells may also have a reinforced cortex at their basal side, increasing the interfacial stiffness [START_REF] Bateman | The receptor tyrosine phosphatase dlar and integrins organize actin filaments in the drosophila follicular epithelium[END_REF][START_REF] Delon | The integrin adhesion complex changes its composition and function during morphogenesis of an epithelium[END_REF]. A more realistic mechanical view of epithelia would thus be that of a bilayer with the addition of stiff interfaces originating from the ECM or the cytoskeleton.

When one layer is a very soft substrate, several studies have attempted to deal with such bilayered systems [START_REF] Biot | Folding of a layered viscoelastic medium derived from an exact stability theory of a continuum under initial stress[END_REF][START_REF] Kücken | A model for fingerprint formation[END_REF][START_REF] Kücken | Fingerprint formation[END_REF][START_REF] Audoly | Buckling of a stiff film bound to a compliant substrate-part iii: herringbone solutions at large buckling parameter[END_REF][START_REF] Reis | Localization through surface folding in solid foams under compression[END_REF][START_REF] Yin | Deterministic order in surface microtopologies through sequential wrinkling[END_REF] by treating the soft layer as an elastic foundation or as an ad hoc resistance potential [START_REF] Kücken | A model for fingerprint formation[END_REF][START_REF] Kücken | Fingerprint formation[END_REF]. These models become essentially single-layered and cannot account for the diversity of behavior of two-layered systems where the two layers play similar roles and have similar properties. We aim here to present a simple mechanical model of thin bilayers using the Föppl-von Kármán equations (FvK) with growth. In the context of a thin plate hypothesis, the FvK equations allow a 3D to 2D reduction, which considerably simplifies the analytical treatment of plate mechanics. The FvK equations are thus well adapted to incorporate many different biological features of layered tissues, which are difficult to handle analytically in a fully three-dimensional formalism. It has recently been shown that FvK equations with growth allow to explain the buckling of thin objects such as flowers and algae in a rather simple way [START_REF] Dervaux | Morphogenesis of growing soft tissues[END_REF][START_REF] Dervaux | Morphogenesis of thin hyperelastic plates: a constitutive theory of biological growth in the Föppl-von Kármán limit[END_REF], when compared to the full treatment of finite elasticity with growth [START_REF] Rodriguez | Stress-dependent finite growth in soft elastic tissues[END_REF]. In addition, a slight modification of these equations allows to treat initially weakly curved membranes. This prompted us to use the same approach for a bilayer with growth. When modeling a single-layered epithelium, the two layers of the model are the ECM and the cells, while the interfacial stiffness represents the basal cytoskeletal cortex of cells (left inset in Fig. 1). When modeling a more complex system, such as the skin, one layer of the model represents the cell layers of the epidermis, the second layer is the dermis (a connective tissue) and the interfacial stiffness represents the extra-cellular basal lamina (right inset in Fig. 1), see reference [START_REF] Ciarletta | Papillary networks in the dermal-epidermal junction of skin: a biomechanical model[END_REF]. Our model can account for local variations in the stiffness or in the thickness of the layers as well as interfacial stiffnesses. Such variations are often present in biological tissues and are thought to be shape generators in morphogenetic processes. These local variations in mechanical properties can also be induced in the context of perturbative experiments.

As much as elasto-mechanical processes may drive the shape of epithelia, this shape builds upon a long and intricate history of growth, stress distribution, and changes in mechanical properties. On the other hand, our formalism, which uses linear elasticity with moderate nonlinear elastic strains, can only account for slight shape variations, spanning time intervals typically smaller than the doubling time of cells within the tissue. To mitigate this problem, we encapsulate all the previous history of development and morphogenesis in an initial shape and a pre-stress. This allows confrontation with experiments where the tissue, at the onset of observation, is rarely in a flat and stressless configuration and where the full history of development is also inaccessible. We show that the model can treat many aspects of biological tissues, which are not so common in material sciences. Although the formalism was initially Fig. 1 Geometry of the bilayer model. w 0 (x, y) is the deviation of the initial configuration from a flat surface; ζ (x, y) is the displacement from this initial configuration as a result of elastic energy minimization under the action of growth and pre-stress (N 0 ). The schema at bottom represents biological context that can be modeled by the theory: the two layers can either be cell or ECM layers of different height, Young's modulus and growth rate; the stiff interface can be a cytoskeletal cortex or a thin ECM developed to treat biological growth, it could also serve for the opposite case of resorption, which could be of interest both in a biological and material science context.

In this paper, we first demonstrate experimentally how alteration of the ECM can substantially impact the shape of an epithelium. For this, we genetically degrade the ECM in a band of cells in the Drosophila wing imaginal disc-the precursor of the adult wing. Subsequent sections are devoted to developing an analytical and numerical treatment of a bilayer model. We aim to account for the experimentally observed distortion of the tissue and provide a general framework to address the mechanics of growing bilayered tissues at large. Section 3 is devoted to the geometry of the sample under study and a reminder of the necessary approximations to incorporate the formalism of finite elasticity with growth into the FvK approach. In Sect. 4, we give the main equations for a bilayer and an approach of the treatment of the interface. In Sect. 5, we derive the Euler-Lagrange equations for the elastic bilayer equivalent to the FvK equations, and we demonstrate the peculiar role of the neutral surface for the bending of the plate. Its position is derived for arbitrary stiffness and thickness, not necessarily constant, of the two layers. Section 6 focuses on uniaxial deformations, which simplify the FvK-modified equations for a bilayer with additive terms at the origin of buckling. In this section, numerical results illustrate the various cases for pre-stressed plates but also for slightly curved membranes, and the theory is confronted with experimental data. In Sect. 7, results obtained with finite element simulations achieved in the same context of growth, thickness, stiffness and defects are presented. Finally, we conclude in Sect. 8 by giving some perspectives.

The model can be adapted to many other biological systems, as epithelia are ubiquitous in most living species. It could also be used for other thin objects such as leaves or algae.

Experimental motivation: biomechanics of Drosophila wing imaginal discs

The Drosophila wing imaginal discs (the precursors of the adult wing) are epithelial tissues that became, over the years, one of the most studied and best characterized systems to study (2022) 137:8 growth [START_REF] Beira | The legacy of drosophila imaginal discs[END_REF]. The growing wing imaginal disc displays a highly patterned field of mechanical stresses [START_REF] Legoff | A global pattern of mechanical stress polarizes cell divisions and cell shape in the growing drosophila wing disc[END_REF][START_REF] Mao | Differential proliferation rates generate patterns of mechanical tension that orient tissue growth[END_REF][START_REF] Rauskolb | Cytoskeletal tension inhibits hippo signaling through an ajuba-warts complex[END_REF]. Cells at the periphery of the epithelium sustain a strong mechanical stretch. This pre-stress builds up as the tissue grows. In addition, the tissue becomes curved as development proceeds, gradually changing from a simple flat surface to a more complex curved surface.

Until recently, most attention was brought upon the role of the apical surface of the epithelium in setting its mechanical properties. It was, however, recently demonstrated that the extracellular matrix upon which the epithelium sits plays a major role in shaping the tissue. Alterations in the basement membrane composition, which includes collagen IV, laminin, nidogen, and heparan-sulfate proteoglycans, have profound effects on the shape of the wing imaginal disc [START_REF] Pastor-Pareja | Shaping cells and organs in drosophila by opposing roles of fat body-secreted collagen iv and perlecan[END_REF]. Alterations of the basement membrane have also been shown to play an active part in setting the 3D shape of the wing disc by promoting the formation of the folds that gradually arise in the tissue [START_REF] Sui | The dorsocross t-box transcription factors promote tissue morphogenesis in the drosophila wing imaginal disc[END_REF][START_REF] Sui | Differential lateral and basal tension drive folding of drosophila wing discs through two distinct mechanisms[END_REF]. From these experiments, it now stands that wing imaginal discs, and most epithelia in general, are two-layered composite structures. Both the cell layer and the ECM layer take an active part in shaping the tissue [START_REF] Ramos-Lewis | Basement membrane mechanics shape development: lessons from the fly[END_REF].

To assess the role of basement membranes in setting mechanical properties of the imaginal discs, we used Drosophila genetics to express the matrix-cleaving metalloprotease Mmp2 in a band of cells, located at the center of the wing imaginal disc along the Dpp-Gal4 genetic pattern (Fig. 2a). The genotype of the observed tissues is ubi-cad: GFP UAS-GFP tub-Gal80 ts /+; Dpp-Gal4/UAS-Mmp2. It combines the necessary transgenes to express the metalloprotease locally and to visualize cells with a cadherin: GFP fusion. The biochemical action of Mmp2 is to degrade structural components of the ECM [START_REF] Page-Mccaw | Matrix metalloproteinases and the regulation of tissue remodelling[END_REF], which, in mechanical terms, implies a possible change of ECM stiffness as well as the thickness. We cannot rule out also an indirect effect on the basal cytoskeletal cortex of cells, as an altered ECM in the vicinity of the cell layer implies a reduced activation of integrin receptors. Our elastic model will need to take into account these different contributions. The metalloprotease expression is controlled in time via the thermo-sensitive mutant of Gal80, by switching Drosophila from 18°C to 29°C 24h before observations (Fig. 2b). Such a timing corresponds to approximately 18 hours of Mmp2 expression since 6 hours are required to reach full Gal4 expression after the temperature switch [START_REF] Bergantiños | Cell death-induced regeneration in wing imaginal discs requires jnk signalling[END_REF]. To image the imaginal discs, we performed ex vivo cultures of wing imaginal discs as in [START_REF] Legoff | A global pattern of mechanical stress polarizes cell divisions and cell shape in the growing drosophila wing disc[END_REF]. The living tissues were then imaged with a spinning disc confocal microscope.

Figure 2c shows the cross section of an imaginal disc 18 hours after metalloprotease expression. The region of perturbation, identified by RFP expression, is shaded in red on the figure. Figure 2d shows the profile of the apical surface of the epithelium, which has a region of inflection near the protease expression. This profile (Fig. 2c,d) corresponds to a wing imaginal disc that was initially flat. Alternatively, Fig. 2d, e shows the same perturbation outcome on a wing imaginal disc that was initially curved-curvature normally builds up at late stages of wing disc development (around 80 hours after egg laying). The resulting tissue is very different from the previous case. This time, the profile shows a more complex configuration: the naturally occurring curvature is prevalent on the borders of the wing disc; at the same time, the metalloprotease action in the region of perturbation induces an inversion of the curvature there.

To conclude, we observed upon ECM degradation (which may impact the stiffness, the thickness, and the interfacial mechanics) that the epithelium curves in the perturbation region (Fig. 2h). We observed two archetypal outcomes depending on whether the tissue was initially flat or curved, which can lead to the inversion of the naturally occurring curvature in the perturbation region. These experimental observations will be subsequently confronted to our FvK bilayer model and our finite element simulations (see Figs.

Modeling

The geometry of the model is represented in Fig. 1: two layers of different thickness (h f and h s ), stiffness and growth rate, are organized in a bilayer, with an interface which may also bear a stiffness. The shape of the initial configuration (see explanations in Fig. 1) is given by w 0 (x, y), the deviation from the horizontal line; the vertical deflection arising from the elastic minimization is denoted ζ (x, y). There are no explicit rules about the difference in stiffness or thickness of the two layers, but they tend to be of the same order of magnitude for biological systems. While there are also no rules regarding the thin interface layer, it will only play a role if it is much stiffer. Contrary to inert materials, changes in the shape of the structure arise from inner processes such as volumetric growth rather than through external loading, although this latter is not eliminated a priori from the model. The bilayer can model, with different degrees of accuracy, a single-or a multilayered epithelium.

A decade ago, Dervaux et al. established the formalism of a growing plate using the theory of finite elasticity with growth [START_REF] Dervaux | Morphogenesis of growing soft tissues[END_REF][START_REF] Dervaux | Morphogenesis of thin hyperelastic plates: a constitutive theory of biological growth in the Föppl-von Kármán limit[END_REF]. Here, we extend their analysis of a single layer to a thin bilayer, both layers having different growth rates and different elastic properties but remaining thin. Indeed, the FvK equations rest on some limitations: the order of magnitude of the vertical deflection ζ due to the buckling and the initial deviation from the horizontal line w 0 must remain small compared to the horizontal length L. (They may be larger than the thickness h of the plate.)

In the following, we first present the model for the general 2D case with arbitrary buckling deformations and subsequently treat the simpler uniaxial folding.

The geometric and elastic strain

We consider two different layers with thickness h f (top) and h s (bottom). When volumetric growth or resorption occurs in an elastic sample, each point is displaced and the gradient of this displacement is a geometric tensor called the deformation tensor: F. According to the famous Kroner-Lee decomposition [START_REF] Goriely | The mathematics and mechanics of biological growth[END_REF], F results from both the elastic tensor A and the growth tensor G in a simple way:

F = AG (1)
Here, plasticity is neglected. The amount of growth per unit volume is given by Det G -1, which is negative for mass resorption and positive for growth. When it is small, the growth tensor G can be decomposed into: G i j = δ i j + ḡi j with | ḡi j | 1. In addition, in the case where the external load maintains a low level of strain, the geometric deformation tensor is then also small leading to the following displacement of the points originally located at r = (x, y, Z ):

u(r) = u 0 (x, y) + (Z -Z n )u 1 (x, y) + w 0 (x, y) = U 0 1 (x, y) + (Z -Z n )U 1 1 (x, y) e 1 + U 0 2 (x, y) + (Z -Z n )U 1 2 (x, y) e 2 + (ζ (x, y) + w 0 (x, y)) e 3 (2)
where w 0 (x, y) represents the initial position of the interface bilayer, and ζ (x, y) the deflection due to the buckling event. Each e i represents a Cartesian unit vector, and the superscripts "0" or "1" stand for the order of perturbation. Indeed, the ratio between the thickness of the bilayer and the horizontal size L is one small parameter 1 = h/L, but the magnitude of the displacement in the Z direction compared to L is a second independent small parameter 2 = ζ /L. We assume that both layers have comparable thickness and h means either h f or h s . In the classical FvK approach, 2 √ 1 a scaling which is justified below and which will constrain ḡi j and w 0 . Since we examine a composite material here, the choice of the origin of the Z axis is free at this stage. We define this origin at the physical interface between both layers. We also introduce Z n , a possible shift in the expansion of order h. This surface defines the surface of separation between parallel surfaces in extension compared to the ones in compression due only to the buckling process. In the single plate geometry, Z n corresponds to the position of the neutral surface and is located at h/2. In the following, an analytical expression for Z n will be derived. We can then define the geometric deformation gradient:

F i j = δ i j + ∂u i ∂ x j (3)
(2022) 137:8 Page 7 of 29 8 where, as in Landau and Lifshitz's book of elasticity [START_REF] Landau | Théorie de l'élasticité[END_REF], the index i indicates an index varying between 1 and 3, while Greek letters restrict to 1 or 2. The elastic strain that we deduce from Eq. (1) when A and G are close to unity reads:

A i j = δ ik G -1 k j + ∂u i ∂ x k G -1 k j = δ i j -ρ i j + ∂u i ∂ x j - ∂u i ∂ x k ρ k j ( 4
)
where we have defined the tensor ρ such that G -1 = Iρ, each component of ρ remaining a small quantity to validate the FvK approach. ρ is not rigorously equal to (G -I). They are equivalent only at linear order. In addition, we assume no change of the thickness of the sample, which means ρ 33 = 0; this hypothesis simplifies the equations and can be easily revisited. A spatially constant value of ρ 33 over the sample simply enters in the definition of h. We also assume that ρ αβ scales as 2 2 , while ρ α3 and ρ 3α scale as 2 . We then introduce a new tensor g whose components are:

g αβ = ρ αβ + ρ βα -ρ 3α ρ 3β and g α3 = ρ 3α + ρ α3 (5)
This tensor is symmetric and has the same scaling behavior as ρ. It will simplify equations later on. In the general case, the growth tensor g is a function of the coordinates, including Z . This is of interest in the context of epithelia, as the cells may not grow in the same way as the ECM substrate. A discontinuity may therefore develop at the interface. Similarly, the structural material coefficients such as the Young's modulus will be different in both layers.

For simplicity, we do not expand these coefficients, and we keep their full variations. We expand only the strain and stress tensors in power of Z -Z n . This assumption is indeed perfectly exact if the bilayer is made of two perfect homogeneous layers. To order 2 in 2 , the elastic strain reads:

       αβ ≈ 1 2 ∂u α ∂ x β + ∂u β ∂ x α -g αβ -ρ 3α ∂u 3 ∂ x β + ρ 3β ∂u 3 ∂ x α + ∂u 3 ∂ x α ∂u 3 ∂ x β -∂w 0 ∂ x α ∂w 0 ∂ x β α3 = 1 2 ∂u α ∂ x 3 + ∂u 3 ∂ x α -g α3 (6)
where we only keep the quadratic term in u 3 = ζ + w 0 (x, y) for the deformation tensor, all the other contributions being neglected. Such simplification is justified by the choice of the traditional scalings of the FvK equations which favors the bending contribution, of order 2 , compared to the in-plane terms of order 2 2 . It is clear that any change in the scaling approximation for the strains as done for the enhanced FvK or Koiter-Sanders model [START_REF] Ruocco | Critical behavior of flat and stiffened shell structures through different kinematical models: a comparative investigation[END_REF] will complicate a lot the equations in the presence of growth. We will not treat this in the present work. For a full analysis of the scalings in initially curved plates, see Chap. 6 of [START_REF] Audoly | Elasticity and Geometry[END_REF]. Our formulation is slightly different from the FvK-Donnell formulation [START_REF] Kücken | A model for fingerprint formation[END_REF][START_REF] Bock | Dynamic contact problem for a von kármán-donnell shell[END_REF], but gives similar results after the variational process of the elastic energy so after integration by part. It recovers perfectly the results of Ciarlet and Miara [START_REF] Ciarlet | Justification of the two-dimensional equations of a linearly elastic shallow shell[END_REF]. The last terms in Eq. ( 6) represent the geometric stretch of the deformation in the x 3 direction.

For a consistent expansion, all terms in this relation must have the same order of magnitude. This implies several scaling relations. We get for the horizontal deformation: u ζ 2 /L, for the growth element ρ αβ ζ 2 /L 2 and for the representation of the initial stress-free configuration, w 0 ∼ ζ .

The monolayer case

We aim here to establish the equilibrium equations of a growing bilayer under initial lateral loading represented by N 0 in Fig. 1, when it is initially planar or weakly curved. To this (2022) 137:8 end, as in Landau and Lifshitz book [START_REF] Landau | Théorie de l'élasticité[END_REF] and in Refs [START_REF] Kücken | Fingerprint formation[END_REF][START_REF] Dervaux | Morphogenesis of growing soft tissues[END_REF][START_REF] Dervaux | Morphogenesis of thin hyperelastic plates: a constitutive theory of biological growth in the Föppl-von Kármán limit[END_REF], our strategy consists in calculating the elastic energy and deriving the Euler-Lagrange equations through variations. Most importantly, we write the stress and strain tensors as a function of x and y only. First, we consider each layer independently since each of them has its own elastic coefficients and its own growth characteristics. Due to the weakness of the deformation for slender objects, we can apply the Hooke's law, which is the constitutive equation of linear elasticity, but we maintain the nonlinearities of the strains according to the strategy of the FvK equations. We also consider the usual membrane hypothesis: σ i3 = 0, which is equivalent to the plane stress approach. Then, it reads for the Cauchy stress components σ αβ where α and β are restricted to x, y:

σ αβ = E 1 -ν 2 (1 -ν) αβ + ν γ γ δ αβ (7)
with ν 1/2 for the incompressible case, an assumption commonly used for living tissues.

In order to determine the elastic energy ϒ, we then write the elastic strain tensor (also a function of x and y), which we decompose in orders of (Z -Z n ) with the help of expansion 2), so that

αβ = 0 αβ + (Z -Z n ) 1 αβ (8)
The zeroth order reads:

0 αβ = 1 2 ∂ α U 0 β + ∂ β U 0 α + ∂ α ζ ∂ β ζ + ∂ α ζ ∂ β w 0 + ∂ α w 0 ∂ β ζ 0 -g αβ -ρ 3α ∂(ζ + w 0 ) ∂ x β -ρ 3β ∂(ζ + w 0 ) ∂ x α (9)
The first order results from the cancellation of the stresses of the third dimension: σ 33 = σ 13 = σ 23 = 0, so 13 = 23 = 0 (see [START_REF] Landau | Théorie de l'élasticité[END_REF]), and using the definition Eq. ( 6), we derive:

U 1 1 (x, y) ≈ - ∂ ∂ x (ζ + w 0 ) + g 13 ; U 1 2 (x, y) ≈ - ∂ ∂ y (ζ + w 0 ) + g 23 (10)
Thus, the strains at first order become:

1 αβ = -∂ 2 αβ (ζ + w 0 ) + 1 2 ∂ α g β3 + ∂ β g α3 (11)
Following the equilibrium equation Eq. ( 7) and the decomposition of the elastic strain (Eqs. ( 9), ( 11)), we can also decompose σ αβ as

σ αβ = σ 0 αβ + (Z -Z n )σ 1 αβ .

FvK equations for the bilayer case

As in the presentation achieved in Landau and Lifshitz's book of linear elasticity [START_REF] Landau | Théorie de l'élasticité[END_REF], we average the elastic energy in the thickness of the bilayer.

Averaging the FvK equation

An important problem in the construction of the model is to establish the position of the neutral surface when the layer is not perfectly homogeneous or for multiple layers. According to Ref. [START_REF] Landau | Théorie de l'élasticité[END_REF], the neutral surface defines the separation between layers in compression from layers in tension when the structure is weakly bent. For an homogeneous sample, the neutral surface is naturally the middle surface of the sample. But for a bilayer, the position of the (2022) 137:8 Page 9 of 29 8 neutral surface is unknown. The interface between both layers seems to be the natural choice for the definition of the origin of the Z coordinate, but it is not the neutral surface. So we focus first on its position, then we proceed to the process of variation of elastic energy.

Position of the neutral surface

Despite the fact that attention was given to this question for functionally graded beams [START_REF] Ruocco | Critical behavior of flat and stiffened shell structures through different kinematical models: a comparative investigation[END_REF][START_REF] Zhang | A theoretical analysis of fgm thin plates based on physical neutral surface[END_REF][START_REF] Larbi | An efficient shear deformation beam theory based on neutral surface position for bending and free vibration of functionally graded beams#[END_REF][START_REF] Cutolo | Third-order nonlocal elasticity in buckling and vibration of functionally graded nanoplates on winkler-pasternak media[END_REF], we did not find the determination of Z n in the literature. Here, we propose a simple argument and we suggest decomposing the full energy into:

ϒ(u α , ζ, Z n ) = 1 2 S 0 d S h f -h s d Z σ 0 αβ 0 αβ + (Z -Z n ) σ 1 αβ 0 αβ + σ 0 αβ 1 αβ +(Z -Z n ) 2 σ 1 αβ 1 αβ (12)
Since ϒ(u α , ζ, Z n ) must remain positive independently of the growth coefficients and the strains, this happens automatically if the intermediate term cancels then imposing:

h f -h s E(Z , x, y)(Z -Z n (x, y))d Z = 0 equivalent to Z n (x, y) = h f -h s E(Z , x, y)Zd Z h f -h s E(Z , x, y)d Z ( 13
)
average position of the bilayer weighted by the Young's modulus. In the general case, there is no guarantee that the neutral surface Z n is positioned inside the bilayer. However, if in each layer the Young's modulus is a constant (but different from top to bottom), the neutral surface is located inside the bilayer. Finally, we derive:

ϒ(u α , ζ, Z n ) = 1 2 S 0 d S h f -h s d Z E(Z ) 1 -ν 2 (1 -ν)( 0 αβ ) 2 + ν( 0 γ γ ) 2 +(Z -Z n ) 2 (1 -ν)( 1 αβ ) 2 + ν( 1 γ γ ) 2 (14)
This method has been employed by Kozlov and Winterhalter [START_REF] Kozlov | Elastic moduli for strongly curved monoplayers. position of the neutral surface[END_REF] for the determination of the neutral surface of strongly curved lipidic membranes.

Variational determination of the Euler-Lagrange equation

The Euler-Lagrange equations result from variations of ϒ with respect to ζ and u α , the linearity between σ and leading to δϒ = drσ αβ δ αβ . These variations, taken one after the other, must vanish at linear order. Special attention must be given if the growth tensor or the Young's modulus E is dependent on Z as it is obviously the case for the bilayer, and sometimes of the other coordinates x and y. Coming back to the definition of ϒ as:

ϒ = S 0 d S h f -h s d Z 1 2 σ αβ αβ - S 0 d Sζ P(x, y) ( 15
)
where P is a vertical external loading. The variation of ϒ with respect to the horizontal deformation u α is easily derived and yields: which gives that ∂ β σ αβ = 0 everywhere in the sample. This corresponds to the second equation in the FvK formalism [START_REF] Landau | Théorie de l'élasticité[END_REF], sometimes re-written with the Airy potential, which is very useful for analytical or numerical solutions of a 2D problem. Variation of ϒ with respect to ζ gives:

δϒ(u α , ζ, Z n ; u α ) = - S 0 d S h f -h s ∂ β σ αβ δu α d Z = 0 ( 1
δϒ(u α , ζ, Z n ; ζ ) = S 0 d Sδζ h f -h s d Z -∂ 2 αβ (Z -Z n )σ αβ -∂ α σ αβ • (∂ β (ζ + w 0 ) -ρ 3β ) -P = 0 ( 1 7 )
In Eq. ( 17), the bracketed term must vanish. We do not mention explicitly the contributions coming from integration by part which fixes the boundary conditions for free boundaries. In the following, we consider only clamped boundary conditions.

The averaged Euler-Lagrange equation for the growing bilayer

It is now possible to transform Eq. ( 17) once we define the new bending coefficient and the new spontaneous curvature due to growth:

D = 1 1 -ν 2 h f -h s E(Z , x, y)(Z -Z n ) 2 d Z; C αβ = 1 2 h f -h s E(Z , x, y) 1 -ν 2 (Z -Z n ) 2 ∂g α3 ∂ x β + ∂g β3 ∂ x α d Z (18)
Doing again the separation of σ into σ 0 and σ 1 , an intermediate step consists in integrating σ αβ over the Z variable, which gives first:

αβ = - h f -h s (Z -Z n )σ αβ d Z = (1 -ν)(D∂ 2 αβ ζ -C αβ ) + ν(D ζ -C γ γ )δ αβ + M αβ ( 19
)
where we have dropped the dependence of D and C αβ in x, y to simplify notations, introduced the 2D Laplacian = ∂ αα , and a new tensor M αβ :

M αβ = h f -h s E(Z , x, y) 1-ν 2 (Z -Z n ) (1 -ν) g αβ - 1 2 ∂ β (ζ +w 0 )ρ 3α - 1 2 ∂ α (ζ +w 0 )ρ 3β +ν g γ γ -∂ γ (ζ + w 0 )ρ 3γ δ αβ d Z (20)
The simplification of M is due to the definition of the neutral surface and the fact that only the growth component g αβ may depend on Z . If it is not the case, M αβ vanishes. At last, Eq. ( 17) gives the equivalent of the first FvK equation:

∂ 2 αβ αβ - h f -h s σ αβ d Z • ∂ α ∂ β (ζ + w 0 ) + h f -h s σ αβ ∂ α ρ 3β d Z = P (21)
This last equation, Eq. ( 21), governs mostly the bending deformation. In the nonlinear regime of the FvK formalism, it is coupled to the horizontal stresses σ αβ whose equilibrium is given by Eq. ( 16). When not stated otherwise, σ αβ is the total stress. D appears to be an effective bending coefficient of the composite structure and C αα as a spontaneous curvature associated with growth. For a homogeneous layer with no growth Z n = h/2, M αβ vanishes and one recovers the traditional FvK equations without growth, see [START_REF] Landau | Théorie de l'élasticité[END_REF]. Notice that in Eq. (21), we also recover the spontaneous curvature κ = ∂ α g α3 found in [START_REF] Dervaux | Morphogenesis of thin hyperelastic plates: a constitutive theory of biological growth in the Föppl-von Kármán limit[END_REF] and the curvatures, inverse of the radii of curvature, of the initial shape given by the tensor ∂ αβ w 0 . Let us focus now on the uniaxial folding where all the equations simplify.

FvK equations governing uniaxial folding

For uniaxial folding, the only dependence is along the x-axis and all the tensors are reduced to one element, which simplifies a lot the analysis but reduces the diversity of observed buckling patterns, as explored in previous works [START_REF] Audoly | Buckling of a stiff film bound to a compliant substrate-part iii: herringbone solutions at large buckling parameter[END_REF][START_REF] Yin | Deterministic order in surface microtopologies through sequential wrinkling[END_REF]. The elastic plane stress equilibrium equation, Eq. ( 16), gives ∂ x σ 11 = 0 so it is only a function of Z and N defined by:

N = h f -h s σ 11 d Z = -N 0 ( 22
)
N 0 indicates a lateral integrated compression, a situation easily realized in material sciences, but in embryos it results from accumulated pre-stress originated from the imaginal disc and the connected tissues at the boundaries [START_REF] Legoff | A global pattern of mechanical stress polarizes cell divisions and cell shape in the growing drosophila wing disc[END_REF][START_REF] Mao | Differential proliferation rates generate patterns of mechanical tension that orient tissue growth[END_REF] (see Fig. 1). The bending equilibrium equation now reads:

d 2 dx 2 D d 2 ζ dx 2 Elastic term + N 0 d 2 (ζ + w 0 ) dx 2 -N 0 dρ 31 dx Pre-stress = d 2 dx 2 (C 11 -M 11 )
Growth contribution

+ P(x)

External stress [START_REF] Sui | The dorsocross t-box transcription factors promote tissue morphogenesis in the drosophila wing imaginal disc[END_REF] where we have neglected the Z dependence of ρ 31 involved in Eq. ( 21). We detail M 11 in the simplified case of homogeneous growth occurring only in each layer.

       M 11 = h f -h s E(Z -Z n ) 1-ν 2 g 11 d Z = h f 0 E(Z -Z n ) 1-ν 2 g f 11 d Z + 0 -h S E(Z -Z n ) 1-ν 2 g s 11 d Z C 11 = h f -h s E(Z -Z n ) 2 1-ν 2 dg 13 dx d Z = h f 0 E(Z -Z n ) 2 1-ν 2 dg f 13 dx d Z + 0 -h S E(Z -Z n ) 2 1-ν 2 dg s 13 dx d Z (24)
We have derived an effective equation for the bilayer, which takes into account that the growth may differ from the top to the bottom layer via C 11 and M 11 . Contrary to C 11 , M 11 appears only if the growth tensor component g 11 is different in both layers. If, in addition, the second derivative of these two quantities does not vanish, they appear equivalent to a vertical loading acting on the bilayer giving it a curvature in the absence of an initial one (w 0 = 0). This will destroy the perfect symmetry up and down, once averaged over the thickness of the bilayer and explain the bending in the case of growth. We now examine in more detail how the structure of the bilayer and eventually how defects affect the bending one-dimensional equation.

Boundary conditions for the bilayer

Boundary conditions concern the top and the bottom interfaces between the surrounding fluid and the two layers, as well as the interface between the two layers (Z = 0). In solid mechanics, the interface between layers is often considered as a line of discontinuity. More realistically, it is a thin zone of sharp variations for physical constants-in our case for growth coefficients or stiffness. Moreover, it may happen that the elasticity of the interface is not simply the mean value of both elastic coefficients. As discussed above, this happens when the interface of a cell layer and the ECM is enriched in intracellular cytoskeletal filaments, or when a thin basal lamina separates two cell layers, both circumstances leading to a stiffer (2022) 137:8

interface. If the interface is really smaller than each layer, we can consider it as a transition zone around Z = 0, then avoiding the application of boundary conditions to a system of 3 layers. A way to avoid writing boundary conditions is to represent the sample as a unique layer where the elastic parameters such as the Young's modulus vary continuously. If the interface is very thin, one needs only to know its thickness relative to the thickness of the bilayer, the details of the description are not really important. The fact that we use the socalled membrane hypothesis with cancellation of the stresses σ 13 = σ 23 = σ 33 = 0 ensures automatically the continuity of the normal and shear stress at the upper and lower boundaries and at the interface of the two layers (Z = 0). So only the continuity of the displacements is required at the interface. Choosing a "continuous approach" and aiming to transform the bilayer into a unique layer whose entire thickness will be h f +h s , we assume that the Young's modulus E and any growth number g αβ can be written as:

E(Z ) = E f + E s 2 + (E f -E s ) 2 tanh(Z /l 0 ) with l 0 <<< h f and l 0 <<< h s ( 25
)
and

g i j (Z ) = g f i j + g s i j 2 + (g f i j -g s i j ) 2 tanh(Z /l 0 ) with l 0 <<< h f and l 0 <<< h s ( 26
)
When a thin stiff membrane superposes to the two main layers of the sample at a position Z m , which can be on top Z m = h f , or at the bottom Z m = -h s , or at the surface of separation: Z m = 0, as justified above, then we can slightly modify Eq. ( 25) into

E(Z ) = E f + E s 2 + (E f -E s ) 2 tanh(Z /l 0 ) + E m e -α 2 m (Z -Zm) 2 /l 2 0
with l 0 <<< h f and l 0 <<< h s [START_REF] Bergantiños | Cell death-induced regeneration in wing imaginal discs requires jnk signalling[END_REF] where the Young's modulus of the membrane, E m , can be much larger than both E f and E s . α m is a numerical coefficient characterizing the thickness of the membrane. Another boundary layer may exist on top or bottom of the bilayer, in addition to the interfacial one. The formulation of the Young's modulus, E(Z ), must then also include this new contribution. Such a definition of E(Z ) can be implemented into Eq. ( 18) very easily to deliver the new FvK equations. We focus now on the neutral surface position and on the bending coefficients when either the elastic coefficients or the thicknesses present a sharp but small in size variation along the sample.

Position of the neutral surface Z n for uniaxial loading

We can guess from the definition of the coefficients (Eq. 24) of the bending FvK equation (Eq. 23) that the position of the neutral surface is crucial for the shape of the sample when growing and buckling. We consider first the position of the neutral surface for layers of constant thickness and constant Young's modulus, and we introduce a linear perturbation dependent on the x variable. If the linear approximation is not valid, a numerical solution is always possible as in Sect. 7.

(2022) 137:8 Page 13 of 29 8

5.1 Position of the neutral surface for a bilayer without structural defects

In the case of an ideal bilayer with no structural defects, Z

N is determined implicitly via Eq. ( 13) with a Young's modulus E(Z ) defined by: E 18) and it reads:

(Z ) = E f if 0 ≤ Z ≤ h f and E(Z ) = E s if -h s ≤ Z ≤ 0. D is deduced from Eq. (
Z (0) n = E f h 2 f -E s h 2 s 2(E f h f + E s h s ) D (0) = E 2 f h 4 f + 4E f E s h 3 f h s + 6E f E s h 2 f h 2 s + 4E f E s h f h 3 s + E 2 s h 4 s 12(1 -ν 2 )(E f h f + E s h s ) (28)
where the superscript " (0) " for Z n and D reminds us that these formulas are restricted to ideal cases. One can notice that if the two layers have the same elastic coefficient, one recovers the standard bending stiffness:

D = E(h s + h f ) 3 /(12(1 -ν 2 )
) and the neutral surface is located at the middle of the layer of thickness (h f + h s )/2. We can convince ourselves that the neutral layer is located inside the bilayer since -h s < Z (0) n < h f , which is a necessary condition. These two quantities Z (0) n and D (0) play a deep role in the bending equation, first equation of the FvK set of equation Eq. ( 21) and for the definition given by Eq. [START_REF] Beira | The legacy of drosophila imaginal discs[END_REF]. It is why we consider now the departure from the ideal situation:

Z n = Z (0) n + δ Z n ; and D = D (0) + δ D (29)
Hereafter, several causes are investigated at linear order in the perturbation amplitude.

Position of the neutral surface for diffuse or stiff interfaces

We allow a diffuse interface according to the representation given by Eq. [START_REF] Ramos-Lewis | Basement membrane mechanics shape development: lessons from the fly[END_REF]. The thickness of the interface l 0 is small compared to h f or h s . Although an exact calculation of Z n is doable, it involves unusual analytical functions so we give here only an asymptotic formula for the result up to a correction of order e -h/l 0 :

Z n = Z (0) n + b l ; b l - π 2 l 2 0 (E f -E s ) 24(E f h f + E s h s ) (30)
This is a tiny effect being of order l 2 0 / h 2 . In the presence of a stiff interface with E I larger than E f or E s , the correction to Z n is more important:

Z n E f h 2 f -E s h 2 s -π 2 l 2 0 (E f -E s )/12 2 E f h f + E s h s + √ π E I l 0 /α I ( 31
)
where the symbol indicates also a correction of order e -h/l 0 , h being either h f or h s . In the presence of an additional stiff layer, on top or on the bottom, with a different stiffness E b , we derive, in the limit of vanishing e -h/l 0 , the following value for the neutral surface position:

Z n E f h 2 f -E s h 2 s + τ b E bj l 0 (h s /α b -l 0 /α 2 b ) -π 2 l 2 0 (E f -E s )/12 2 E f h f + E s h s + √ πl 0 (E b /(2α b ) + E I /α I ) (32)
with τ b = 1 when the stiff layer is on top, and τ b = -1 when the stiff layer is at the bottom.

Knowing that the interface stiffness E I or E b may be an order of magnitude larger than E f 123 or E s , this correction may modify significantly the values of the D coefficient in the FvK equations. The addition of thin layers, as described here, to the initial bilayered system will change the position of the neutral surface and the resulting D coefficient, but it will not change the structure of the FvK equations.

Localized defects of the thickness

Let us now consider a weak and localized variation of the thickness of one of the two layers.

The neutral surface will then be distorted. We derive its new position by a simple linear expansion of Z n . This induces an equivalent expansion on the parameter D, which multiplies the bending term in equation Eq. [START_REF] Sui | The dorsocross t-box transcription factors promote tissue morphogenesis in the drosophila wing imaginal disc[END_REF]. Assuming that the thickness varies as

h = h i (1 + β i (x)) then Z n Z (0) n + ∂ Z (0) n ∂h i h i β i (x) and D D (0) + ∂ D (0) ∂h i h i β i (x). (33)
where D (0) and Z (0) n are given by Eq. ( 28). and means an expansion restricted to first order. It reads -On the top layer, we get:

Z n Z (0) n + E f h f -Z (0) n (E f h f + E s h s ) h f β f (x); and D D (0) + E f (1 -ν 2 ) h f -Z (0) n 2 h f β f (x) (34)
-On the lower layer:

Z n Z (0) n - E s h s + Z (0) n (E f h f + E s h s ) h s β s (x); and D D (0) + E s (1 -ν 2 ) h s + Z (0) n 2 h s β s (x) (35)
It is to be noted that a dip in either the upper or the lower layer means a negative value for the coefficients β i (x) so, as expected, a local decrease of the stiffness of the sample.

Localized inhomogeneity of the Young's modulus

We now consider a localized variation of the Young's modulus, which we treat through the expansion:

E E i (1 + δ i (x)) and Z n Z (0) n + ∂ Z n ∂ E i E i δ i (x) and D D (0) + ∂ D ∂ E i E i δ i (x) (36)
Leading to:

-Inhomogeneity in the top layer: -Inhomogeneity in the bottom layer:

Z n Z (0) n + E f E s h f h s (h f + h s ) (E f h f + E s h s ) 2 δ f (x); and D D (0) + E f 3(1 -ν 2 ) h 3 f -3h 2 f Z (0) n + 3h f (Z (0) n ) 2 δ f (x) (37)
Z n Z (0) n - E f E s h f h s (h f + h s ) (E f h f + E s h s ) 2 δ s (x); and D D (0) + E f 3(1 -ν 2 ) h 3 s + 3h 2 s Z (0) n + 3h s (Z (0) n ) 2 δ s (x) (38)

Results for uniaxial folding

As mentioned above, the uniaxial case simplifies a lot the analytical analysis but reduces the diversity of buckling patterns that can be accounted for compared to the 2D case. Nevertheless, the 1D approximation is still of great experimental relevance for some experimental systems. In the case of the experiments of Fig. 2, for example, the tissue and the genetic perturbation can be assumed to be spatially invariant along the y-coordinate, making it essentially a 1D problem. In this context, the present section provides an in-depth analysis of bending contributions of the bilayer in the uniaxial case, with and without structural defects.

We first analyze the bending equation of the uniaxial folding analytically. We then address the problem numerically by selecting possible 1D examples using the software of resolution Mathematica. Finally, simulations from the FEM (finite element method) software COMSOL Multiphysics are presented.

Source of bending in bilayers

We will demonstrate in the following that the source of bending for growing bilayers is numerous and diverse and may have different biological origins. We analyze some cases here as examples. First, we define dimensionless parameters which will govern our FvK equations for a bilayer made of two layers with constant elastic coefficients.

Dimensionless bending equation

Using as length unit the length L of the sample, Eq. ( 23) reads:

d 2 d X 2 (1 + δ D) d 2 ζ d X 2
Elastic term

+ N 0 d 2 (ζ + W 0 ) d X 2 - dg 31 d X Pre-stress = d 2 d X 2 (C -M)
Growth contribution [START_REF] Bayly | A cortical folding model incorporating stress-dependent growth explains gyral wavelengths and stress patterns in the developing brain[END_REF] where we keep the same notation for ζ (X ), which becomes

ζ (x)/L, X = x/L, N 0 = N 0 L 2 /D (0) , δ D = (D -D (0) )/D (0) , M = M 11 L/D (0) , C = C 11 L/D (0)
are dimensionless quantities and W 0 (X ) = w 0 (x)/L represents the initial position of the curved sample. External vertical loading is discarded. It is worth noting, however, that the equivalent of a loading can arise from the growth contribution (right-hand term of Eq. ( 39)) owing to the fact that it involves spatial second derivatives, like the pre-stress term. For such a growth-induced loading to arise, the growth of the two layers must differ to get a non-vanishing M value, and at least one of the terms (C or M) must vary spatially. This growth-induced loading accords with previous investigations that demonstrated the need for differential growth of apposed layers to induce buckling of brain cortical folds for example [START_REF] Bayly | A cortical folding model incorporating stress-dependent growth explains gyral wavelengths and stress patterns in the developing brain[END_REF][START_REF] Ben Amar | Mimicking cortex convolutions through the wrinkling of growing soft bilayers[END_REF][START_REF] Ambrosi | Growth and remodelling of living tissues: perspectives, challenges and opportunities[END_REF][START_REF] Weickenmeier | Dimensional, geometrical, and physical constraints in skull growth[END_REF]. Even with homogeneous growth, the right-hand term can lead to a non-vanishing contribution, for example, in the case of a structural defect that will change Z n , and induce high spatial frequency components in C or M. Nevertheless, it is difficult to discuss a priori the ζ profile, which is added to W 0 as a result of the growth term, since it depends also on the boundary conditions. We can intuitively predict that inhomogeneous growth is responsible for a "fictitious pressure" or on the contrary for a tension added to the sample. The simplest case may be the one of a homogeneous bilayer that we consider first.

Bilayer without structural defects

This case corresponds to Sect. 5.1. For simplification, we consider the case of a sharp interface and the bilayer is made of two perfect layers of constant thickness and Young's modulus.

The growth process occurs mostly in one layer. Assuming first g f 11 (X ) in the upper layer and g s 11 (X ) = 0 in the lower layer, we have for M f 0 and C f 0 according to Eq. ( 24):

M f 0 = g f 11 (X ) E f h f (h f -2Z (0) n ) 2D (0) (1 -ν 2 ) ; C f 0 = dg f 13 d X E f h f (h 2 f -3Z (0) n (h f -Z (0) n ) 3D (0) (1 -ν 2 ) (40)
where Z (0) and D (0) are given by Eq. ( 28). It is easy to check that both parameters M 0 and C 0 have the same sign as the growth elements g f 11 and ∂ X g f 13 . If growth occurs in the lower layer and not in the upper layer, it reads

M s 0 = -g s 11 (X ) E s h s (h s + 2Z (0 n ) 2D (0) (1 -ν 2 ) C s 0 = dg s 13 (X ) d X E s h s h 2 s + 3Z (0) n (h s + 3Z (0) n ) 3D (0) (1 -ν 2 ) ( 41
)
Here again,

h s + 2Z (0)
n is a positive quantity and the sign of M s 0 is opposite to the sign of g s 11 , while the sign of C s 0 is automatically given by the growth coefficient derivative. Without pre-stress, it is possible to observe a buckling of an initially flat plate because of a differential and inhomogeneous growth process. There exists a competition of the origin of this buckling between g 11 or g 13 , or between M 0 and C 0 . Differences in growth of the two layers induce a symmetry breaking between up and down. It is nontrivial to compute the sign of the growth term-whether it contributes to a vertical pressure or tension-owing to the fact that it stems from the competition of two terms (M 0 and C 0 ) through their second derivatives.

Notice, however, if the growth is not x dependent, these terms will disappear from the buckling equation; on the contrary they will become more efficient if they strongly depend on x. Localized defects will increase the efficiency of the buckling as shown hereafter.

Analysis of structural defect in growing bilayers a) Localized stiffness variation in one layer

We first consider that the stiffness can be locally modified as discussed in Fig. 2, inducing a small change in the position of the neutral surface. We focus on the linear variation of M f or M s , and C f or C s , simply deduced from Eqs. [START_REF] Ben Amar | Mimicking cortex convolutions through the wrinkling of growing soft bilayers[END_REF][START_REF] Ambrosi | Growth and remodelling of living tissues: perspectives, challenges and opportunities[END_REF]. The modification of these quantities is deduced taking into account the variation of δ Z n reported by Eqs. [START_REF] Kozlov | Elastic moduli for strongly curved monoplayers. position of the neutral surface[END_REF][START_REF] Boudaoud | Dynamics of singularities in a constrained elastic plate[END_REF] Sect. 5.4, which reads:

δ M f = ∂ M f 0 ∂ Z (0) n ∂ Z (0) n ∂ E i δ E i ; δC f = ∂C f 0 ∂ Z (0) n ∂ Z (0) n ∂ E i δ E i (42)
where E i is the Young's modulus of the layer affected by the defect: E f or E s and δ E i is represented by E i δ i (X ), in a similar way to Eqs. [START_REF] Kozlov | Elastic moduli for strongly curved monoplayers. position of the neutral surface[END_REF][START_REF] Boudaoud | Dynamics of singularities in a constrained elastic plate[END_REF]. Similar results apply for δ M s which corresponds to the growth in the lower layer and for δC s . We give here only the results but we define first the following positive quantities:

m f = E 2 f E s h 2 f h s (h s + h f ) 2D (0) (1 -ν 2 )(E f h f + E s h s ) 2 c f = E f E s h f h s (h f + h s ) 2 2D (0) (1 -ν 2 )(E f h f + E s h s ) 3 (43)
-When the growth occurs in the upper layer, which also exhibits a stiffness defect, then

     δ M f = -g f 11 (X )m f δ f (X ); δC f = - dg f 13 dx c f δ f (X ) defect on top δ M f = g f 11 (X )m f δ s (X ); δC f = dg f 13 dx c f δ s (X ) defect on the bottom (44)
-When the growth occurs in the substrate, perturbation occurs now on M s 0 and C s 0 and one needs to consider again a defect of the stiffness either on top or on the lower layer.

     δ M s = -g s 11 (X )m f E s h s E h h f δ f (X ); δC s = dg s 13 dx c f δ f (X ) defect on top δ M s = g s 11 (X )m f E s h s E h h f δ s (X ); δC s = - dg s 13 dx c f δ s (X ) defect on the bottom (45)
b) Localized variation of the thickness of one layer.

A defect in the thickness affects the position of the neutral surface, but also, when the thickness defect occurs in the same layer as the growth, the averaged value of the coefficients of Eq. ( 24) related to the first bending FvK equation, Eq. ( 23). We give here only the schema when a dip appears in the lower layer, the growth occurring either on the top or the lower layer. A similar analysis gives

δ M f = d M f 0 d Z (0) n δ Z n = g 11 (X )M f β s (X ); δC f = dC f 0 d Z (0) n δ Z n = dg f 13 d X C f β s (X ) (46)
whereas in Sect. 5.3, β s (X ) describes the shape of a notch and M f and C f are dimensionless positive quantities.

M f = E f E s h f h s E s h 2 s + E f h f (2h s + h f ) 2D (0) (1 -ν 2 )(E f h f + E s h s ) 2 ; C f = E s h s (h f + h s ) E h h f + E s h S M f ( 47
)
We examine now the case where growth occurs in the bottom layer where also the defect is localized: .

δ M s = -g s 11 (X )M f β s (X ); δC s = d s g 13 dx C s β s (X ) (48)
with

C s = E 2 s h 3 s + E f E s h f h s (-h f + h s ) + E 2 f h 2 f (h f + 2h s ) 2E f h f (E f h f + E s h S ) M f (49)
c) Shape of defects 123 Fig. 3 Cylindrical mylar shell under a localized forcing. Originally the plate has a length of λ = 35 cm and a width L = 17.5 cm. The Young's modulus of the mylar plate is 3.810 9 Pa.s, and the Poisson ratio is ν = 0.4. The initial bending is realized by imposing the boundary conditions. On the left, one notices that for a low value of the imposed deviation Z = 4 mm, the pattern of sinking is quite elliptical, while it becomes cylindrical for Z = 15 mm representing a folding in the opposite direction of the initial one. Image extracted from [START_REF] Boudaoud | Dynamics of singularities in a constrained elastic plate[END_REF] Localized defects can be positioned in any place of the sample. From Eq. ( 39), we have noticed that they may induce an equivalent forcing on the sample if they contribute to strong variation of the second derivatives of the coefficients δ M and/or δC. Choosing a sharp defect in the center of the sample, a good representation may be δ = e Y e -X 2 / 2 where e Y is the relative amplitude of the stiffness defect, which can be positive or negative, 2 depicts qualitatively its width. The second derivative is ∂ X X δ, which is of order -2e Y / 2 which for a defect of width 0.05 corresponds to an amplitude of 200e Y at the center. It is to be noted that the sign is opposite to e Y . So the sign of the second derivative of δ M or δC can be easily found from paragraph (a). The convention for thickness defect is the following:h i β i (X ) where β i (X ) is negative for a notch and positive for a protrusion. Here again, ∂ x x β will be opposite to β. In any case, it remains that the prediction of the global sign of the second member of Eq. ( 39) is hard to predict due to the competition between g 11 and g 13 . Let us remember that a positive sign of a localized contribution will be equivalent to a localized forcing, which will mimic the experiment [START_REF] Boudaoud | Dynamics of singularities in a constrained elastic plate[END_REF] shown in Fig. 3. 6.2 Numerical investigations of the FvK bilayer for a plate or a shell.

In this section, we aim to demonstrate the role of defects on plates and slightly curved shell portions under pre-stress. To simplify, we focus first on a homogeneous growth process: It means a process independent of X and identical in both layers or equivalently a non-growing sample. Then, we will give numerical examples to illustrate the analytical results of the previous section when differential growth occurs contingently to defects.

Local decrease of thickness on a non-growing pre-stressed plate or shell

In this section, we will evaluate the role of a sharp variation on the D coefficient of a prestressed sample initially flat or curved. For the numerical purpose, the defect has the shape of a penetrating peak and is represented by a Gaussian so δ D = -D 0 e -(X -X d ) 2 /l 2 d . D 0 is a positive coefficient which can be evaluated with either Eq. [START_REF] Zhang | A theoretical analysis of fgm thin plates based on physical neutral surface[END_REF] or Eq. ( 35), l d is an estimation of the scale of the defect compared to the scale of the sample L and X d is its relative position in the sample. For a plate, a curvature is initially induced by the pre-stress and applying a slope at both ends. In Fig. 4a, the imposed slope is ∂ X ζ ± 0.03 at X ± 0.5 giving a deflection below the horizontal (black line). This deviation increases with the defect since it weakens the stiffness more and more with the amplitude D 0 = 0.2 (blue curve), D 0 = 0.3 amplitude (red curve). Even for a sharp-pointed defect, with l d = 1/10, the deviation ζ remains small and not very selective. In addition, a displacement on the position of the defect does not affect ζ too much (magenta curve). For a shell which is initially curved in the positive Z direction with clamped conditions at both ends, as shown in Fig. 4b,c, the deviation ζ is amplified by the initial curvature and the pre-stress. Here again, any deviation on the right of the defect weakens its role on ζ since we cannot distinguish the deviation with lateral defect (magenta line) from the black line. When the whole profile = ζ + w 0 (X ) is considered, we clearly observe an increase of the total curvature when the defect is centered at the maximum of deflection.

Thickness defects on growing plates and shells

The thickness variation on δ D gives a very intuitive result. A much less intuitive result concerns the sign of ∂ X 2 δ M and ∂ X 2 δC of Eq. ( 39) since the results are dependent of the second derivative of these growth coefficients. These quantities have been studied in detail in Sect. 6.1.3. We must notice that they have an opposite sign a priori, and in practical situations, we will have little information except by comparison with genetically modified tissues. So we will join these two contributions into a unique term δS G = -(δ M -δC) = s g e -(X -X d ) 2 /l 2 d in the right-hand side of Eq. [START_REF] Bayly | A cortical folding model incorporating stress-dependent growth explains gyral wavelengths and stress patterns in the developing brain[END_REF]. Intuitively, as soon as this term is differentiated 2 times, it behaves as either a pressure if it is positive or a tension in the opposite case, while the initial curvature of the concave shell behaves like a pressure and a convex one like a tension. The same definition and color codes apply as in Fig. 5: black no defect, blue continuous s g = 0.5,blue discontinuous s g = -0.5, red continuous curve s g = 1 and l d = 0.1, red discontinuous s g = 1 for the 3 panels. Clamped boundary conditions are applied. In (a), the shell is initially concave and l d = 0.1. In (b), the shell is initially convex. In (c), the defect is localized on the right x d = 0.2 and the shell is concave. W 0 (X ) = (1/4 -X 2 )/5., for concave shape and W 0 (X ) = (X 2 -1/4)/5 for convex shape Knowing that a groove in the substrate is represented by a negative profile, the second derivative is then positive for X = X d being given by ∂

X 2 δS R ∼ -2/l 2 d s g e -(X -X d ) 2 /l 2 d .
Illustration is given in Fig. 5 for a pre-stress plate with defect localized at X = 0 with two different length scales: l d = 0.1 and l d = 0.3 and different amplitudes for s g = ±0.5 and ±1. A defect put on the right is also shown with a noticeable distortion in panel (c) of Fig. 5. In Fig. 6, the same set of perturbations act on a shell either concave (panel(a)) or convex (panel(b)), and the deviation ζ is plotted with boundary conditions ζ = ∂ X ζ = 0 for X ± 0.5. For the concave case (a), the deviation due to the defect is decreased for s g > 0 and increased otherwise. The opposite result is obtained for in the convex case panel (b). The distortion is shown in panel (c) only for the concave case. Once the deviation is superposed on the initial shell, depending on the respective sign of the deviation versus the shell geometry, it is possible to observe an inversion of the curvature. This is demonstrated in the three panels of Fig. 7. Here, we confront our bilayer FvK formalism in the uniaxial geometry to experimental profiles extracted from the images of Fig. 2c,e. We first focus on the tissue presented in Fig. 2c, that we rescale in the unit of the lateral size of the imaginal disc (Fig. 8a). The perturbation of the ECM was induced before the naturally occurring curvature could develop in the tissue. We therefore assume no initial curvature: a horizontal plate only under pre-stress. We fix the slope conditions at the border, making the hypothesis that the tissue periphery is not affected by the perturbation. This hypothesis is justified by the fact that the defect is well localized at the center of the tissue. The fit of parameters is achieved by a limited number of trials. Fitted parameters converge to the following values: a pre-stress of order 9 in unit of D and a parameter D 0 = 0.3; a defect characterized by an amplitude s g = 4, position x d = 0.1, and width l d = 0.25, when the defect is modeled by the relation s g e -(X -X d ) 2 /l 2 d , see Sect. 6.2.2. Finally in Fig. 8b, the shape of the disc at the boundaries and the developmental timing of the perturbation suggest that the tissue was already curved at the onset of the perturbation. We then used the shell treatment with a best guess for the unperturbed shape w 0 (X ) = 0.7(1/4 -X 2 + 4(1/16 -X 4 )), N 0 has little effect and is put to zero, which suggests that in old epithelia, the outer stresses have relaxed, which is not the case for young epithelia. The amplitude of the defect has increased: it is now s G = 10; localization is X d = -0.05, and width l d = 0.1. Thus, the bilayered FvK theory can account for the mechanics at play in the growing wing imaginal disc in the presence of a localized defect in the ECM, both when the tissue is initially flat and curved.

Numerical simulations with finite element method

We now consider the problem of the bilayer bending in the linear Hookean elasticity approach using COMSOL Multiphysics, which computes the deformations in the real geometry of the bilayer, two different connected layers, without averaging the elastic properties. The bilayer we study in the following is made of a thin (h f = 10µm) and relatively soft (E f = 50 kPa) upper layer and a thicker (h s = 12µm) and stiffer (E s = 75 kPa) bottom layer. The length Theoretical fit of an older wing imaginal disc, for which the initial profile before perturbation was curved (same as Fig. 2e) of the bilayer is L = 500 µm. The layers are submitted to growth and pre-stress, and the boundaries are either free or constrained. We investigate a defect, as a local removal of the basal layer in the middle of the bilayer with an extension of 10 µm and a thickness h d ≤ h s , and as a local change of stiffness and growth in the bilayer. The growth is anisotropic and is introduced in a tensorial form: G i j = 1 + g δ 1i δ 1 j , with g ∼ 0 -0.1. The growth may be introduced in the top layer only as in Sect. 7.2, the bottom layer only in Sect. 7.4, or in both layers as in Sect. 7.3. The pre-stress consists in an anisotropic compression or tension, so that σ P S i j = -σ P S δ 1i δ 1 j where σ P S ∼ 0 -500 Pa. We also study the case of the shell on the further deviation of the sample. We set an initial curvature to the shell:

W 0 = H × (L/2) -2 ((L/2) 2 -X 2 )
where H is the displacement at X = 0. In order to compare the different results, we plot the profile of the interface between the upper and lower layers (the red line in Fig. 9e).

Finally, we fit FEM simulations to the experiment described in Sect. 2. The result is compared to the one obtained from the FvK calculation in Sect. 6.3.

Pre-stress on an initially curved bilayer

In Fig. 9a, we investigate different initial shell shapes with different curvatures characterized by H = -20 µm,0 µm, 20 µm with the same pre-stress σ P S = 200 Pa in the bilayer. For a plate, with no initial curvature, a bending toward the stiffer and thicker basal side is obtained, while for a convex or a concave shell, the initial curvature of the shell is reinforced in absolute value. This is consistent with the previous FvK model for shells. In fact, the bending deviation follows the initial shape of the shell.

Thickness and growth defect on growing plates and shells

In biological systems, the two layers in epithelia are not independent. For instance, the basal membrane regulates the proliferation and the differentiation in the epithelium. In general, the basal layer structures the upper layer of the epithelium, and a defect in the lower layer can have consequences on the proliferation and metabolism of the upper layer. Therefore, we now investigate the shape change at the level of a defect, which alters the thickness of the basal layer and consequently the growth of the upper growing layer. Figure 9b displays the results when the basal layer is removed in proportion of 85% (dotted lines) in thickness and the upper layer is not growing at the defect level, as well as the case with no defect (continuous lines). We observe a buckling toward the basal membrane, in the neighborhood of the defect.

Place of a defect

Defects can appear anywhere along a bilayer. In Fig. 9c, we impose defects on a planar bilayer, located at X = 0, 0.08, 0.16 for the blue, green, and red lines, respectively. Both layers grow to the same extent g = 0.1. The defect consists in removing a part of the bottom layer. Interestingly, with the same size of the defect, the amplitude of the deformation is different, and the shift of the defect is not only a shift in the deformation.

Dependence on boundary conditions

The comparison between confined and free geometries is a recurrent topic in cellular growth studies [START_REF] Alessandri | Cellular capsules as a tool for multicellular spheroid production and for investigating the mechanics of tumor progression in vitro[END_REF][START_REF] Kiran | Distinct modes of tissue expansion in free versus earlier-confined boundaries for more physiological modeling of wound healing, cancer metastasis, and tissue formation[END_REF] from multiple perspectives, such as the folding of epithelial sheets [START_REF] Inoue | Epithelial tissue folding pattern in confined geometry[END_REF].

The boundary conditions of bilayer plates or shells are also necessary for the analytical and numerical study of the system. But unfortunately, they are not always perfectly known. When the setup is symmetric, one can use those symmetries to constrain the problem, and the system has free-boundary conditions. Sometimes, the experiment is such that the system is clamped, or loses its symmetry. Confined systems are strongly dependent on boundary conditions. In Fig. 9d, we compare the buckling in a constrained and free geometry. The buckling is very different, and its direction (up or down) even changes when starting with a concave geometry. It may lead to strong inaccuracies in FEM simulations of rectangular geometry and semi-analytical treatments become necessary, see reference [START_REF] Ruocco | An analytical model for the buckling of plates under mixed boundary conditions[END_REF]. Notice that in the previous sections, Sects. 6.2.1, 6.2.2, boundary conditions are always fixed symmetrically for plates and clamped boundary conditions for shell deviations. This point is important for comparing our numerical results to the experiments as shown in the next section. There can be a few types of constrained geometry. In Fig. 9d, it is obtained by imposing a zero displacement at two points on the lateral boundaries. In Fig. 9f, the whole lateral boundaries are imposed a zero displacement. In the second case, the Z-displacement is close to having a zero derivative as boundary conditions.

Confrontation of the finite element simulations with experiments

We use the finite element method to reproduce qualitatively the side to which the imaginal wing curves. We assume that the stress inside the cell layer is induced by growth and is relaxed when a defect is introduced. Since in classical mechanics only the stress-free configuration and the configuration obtained after the process are to be taken into account, we make a multilayer growth with and without defects. In our simulations, we observed that a 2-layer model could not reproduce the experimentally observed change of curvature upon ECM degradation. Only when an apical membrane is introduced in the model, can we reproduce the behavior-making it essentially a 3-layer model. Indeed, contrary to the resolution of the FvK equation where boundary conditions are imposed, we use free boundary conditions with the finite element method. Therefore, a third layer has to be added in order to constrain the system to bend toward the correct side.

More precisely, we simulated a thick (15 µm) cell layer of length L = 500 µm, growing between two non-growing thin layers. This process creates a compressive stress in the two non-growing layers. The top layer has a thickness of 5 µm, and the lower layer has a thickness of 2 µm. The apical membrane may represent the apical cortex, and the apical adherens layer dominates apical mechanics in epithelia [START_REF] Pinheiro | Mechanical force-driven adherens junction remodeling and epithelial dynamics[END_REF], whereas the basal membrane represents the ECM. The growth deformation gradient is written as: G i j = 1 + δ 1i δ 1 j g 11 . The boundary conditions are free. We fit the experiment (Fig. 11a,b), and we find a basal ECM stiffer (E b = 100 kPa or E b = 30 kPa) than the upper membrane and that the cell layer (E = 10 kPa). In the case of the early imaginal wing, g 11 = 0.2 and E b = 30 kPa, and for the old imaginal wing g 11 = 0.25 and E b = 100 kPa. We also introduce an asymmetry for the place of the defect at X = -0.025 for the old epithelium and X = 0.15 for the early epithelium. The apparent growth is thus higher and the ECM stiffer for the old imaginal wing when compared to the early imaginal wing. This setup results in a bending toward the upper layer in the case with no defect for the old imaginal wing and a slight bending toward the basal side for the early imaginal wing Fig. 10a,c. However, when a defect of length l d = 100 µm is introduced in It seems that for the early imaginal disc this model is sufficient to properly reproduce the experimental profile Fig. 11a. However, the agreement was not as good for the older initially curved imaginal disc in Fig. 11b. We then modified the FEM model with the same strategy as for the FvK model, by changing the stiffness at the level of the defect, making the cell layer softer. This refined model provided a better fit to the experiments in Fig. 10c,e. The new values of the fit for the old imaginal wing are: g 11 = 0.43 and E = 1.7 kPa in the cell layer at the defect level. The other parameters remain unchanged. Contrarily, changing the growth rate at the level of the defect does not improve the fit. We deduce that the contribution s g of the defect in Sect. 6.3 is mostly due to a change of the stiffness in the growing cell layer at the level of the defect, rather than a change of growth. This is consistent with the duration of Mmp2 expression during the experiment (18-24 h), see Sect. 2. No pre-stress was introduced in the FEM simulations, since this pre-stress is assumed to be caused by growth, and the growth can be introduced explicitly with an order of magnitude larger than with FvK. For the same reason, growth generates the curvature, which is not introduced explicitly. To conclude, FEM simulations validate the FvK calculations.

Discussion and conclusion

Motivated by experimental evidence that the ECM plays an important role in shaping epithelia in 3D, as shown in Fig. 2, this work provides a mechanical description of growing epithelia (2022) 137:8 Theoretical fit of an early wing imaginal disc, for which the initial profile before perturbation was flat (same as Fig. 2c). (b): Theoretical fit of an older wing imaginal disc, for which the initial profile before perturbation was curved (same as Fig. 2e).

The stiffness of the cell layer is not altered by the defect. (c): Comparison between the refined fit when the stiffness of the cell layer is altered by the defect (in green) and the case when it is not (red)

in the Föppl-von Kármán framework. In the context of embryo-genesis and organo-genesis, previous theoretical studies have modeled epithelia using the vertex model [START_REF] Farhadifar | The influence of cell mechanics, cell-cell interactions, and proliferation on epithelial packing[END_REF][START_REF] Staple | Mechanics and remodelling of cell packings in epithelia[END_REF][START_REF] Graner | Simulation of biological cell sorting using a two-dimensional extended potts model[END_REF] inspired from the physics of foams [50]. Indeed, epithelia are cellular pavings where the main unit has a polygonal shape in 2D or a shape of polyhedrons in 3D. Vertex models approximate cells as polygons, make hypotheses about the mechanics of individual cells (e.g., their growth, surface or line tension and stiffness) and simulate the behavior of cell assemblies from these elementary rules. These simulations are based on energy minimization, as in the present work. One limitation of cell-based models is the fact that simulating cell assemblies may become difficult at very large cell numbers, notably when the parameter space must be explored. In addition, it is rather hard to extract macroscopic physical quantities such as Young's modulus or surface tension from microscopic parameters describing the cells. Most of the aforementioned theoretical works focused on the 2D mechanics of the apical surface of epithelia, with only few studies addressing the 3D aspects of tissues [START_REF] Bielmeier | Interface contractility between differently fated cells drives cell elimination and cyst formation[END_REF][START_REF] Hannezo | Theory of epithelial sheet morphology in three dimensions[END_REF]. Epithelia are usually thin layers, and the elasticity of slender elastic objects such as plates, tubes or membranes can be good candidates to describe their mechanical properties, especially in the presence of bending. However, a single plate may be too simple a model when it comes to explaining the behavior of epithelia, which are intrinsically multilayered systems. In the present paper, we have considered the multicomponent structure of epithelia from the view point of elasticity. Local changes in elastic coefficients and differential growth are at the origin of stresses able to buckle the sample. Not only could our bilayered elastic model predict a broad range of morphogenetic behavior, it also yielded surface profiles that were in very good agreement with experimentally observed ones. Notably, it could account for the localized degradation of the ECM on two experimental observations that were limit cases: one initially flat and one initially curved wing imaginal discs (Figs. 2,8). The morphogenetic processes that we observed and modeled-fold formation upon ECM degradation-were artificially induced via an ectopic genetic expression of the metalloprotease Mmp2. These experiments, as such, fall in the realm of "synthetic morphogenesis". Nevertheless, similar processes are naturally occurring in tissues. For example, deep folds that develop at late stage of the wing imaginal disc are thought to arise through such a local degradation of the ECM [START_REF] Sui | The dorsocross t-box transcription factors promote tissue morphogenesis in the drosophila wing imaginal disc[END_REF][START_REF] Sui | Differential lateral and basal tension drive folding of drosophila wing discs through two distinct mechanisms[END_REF].

In this paper, we performed experimental observations of tissues exclusively after the action of the genetic perturbation. Ideally, one would need to image the tissues before and after the onset of the perturbation to completely disentangle the origin of the buckling. This is experimentally possible in Drosophila with chronic imaging [START_REF] Heemskerk | Dynamic clonal analysis based on chronic in vivo imaging allows multiscale quantification of growth in the drosophila wing disc[END_REF] and potentially through prolonged imaging (i.e not just a few snapshots) with an appropriate method to reduce the phototoxicity associated with the light dose [START_REF] Abouakil | An adaptive microscope for the imaging of biological surfaces[END_REF].

(2022) 137:8

Page 27 of 29 8 In the uniaxial geometry (with only one dimension for the sample), defects generate buckling distortions of the shape due to variation of thickness, stiffness but also of growth process. For a 1D bending process, the in-plane elastic equation is trivially solved and only the bending equation exhibits interesting degrees of distortions. It is not the case for 2D modes of deformations where spatial patterns happen on the outer periphery of the samples ( [START_REF] Dervaux | Morphogenesis of growing soft tissues[END_REF][START_REF] Dervaux | Morphogenesis of thin hyperelastic plates: a constitutive theory of biological growth in the Föppl-von Kármán limit[END_REF]). One can then expect interesting patterns in these geometries.

  8 and 10).

Fig. 2

 2 Fig. 2 Folding an epithelium by local ECM degradation. a Transmission light imaging of a Drosophila wing imaginal disc (left) and corresponding sketch with the region of localized expression of the protease Mmp2 shaded in red (right). b To control expression of the metalloprotease, the sample is switched from 18°C to 29 °C several hours before observation. c-f Experimental demonstration. c An xz cross section demonstrates the shape of an initially flat wing imaginal disc under the effect of a localized metalloprotease expression (in the red colored region). The schema in d demonstrates the change in shape of the tissue. e, f Same as c, d but with an initially curved wing imaginal disc. g A schematic representation of the shape changes observed in the wing imaginal disc

Fig. 4

 4 Fig.[START_REF] Pinheiro | Mechanical force-driven adherens junction remodeling and epithelial dynamics[END_REF] Buckling deviation due to a variable stiffness coefficient D(X ), in a pre-stressed non-growing plate or shell. In (a), the profile ζ of an initial elastic plate under pre-stress N 0 = 4 is plotted as a function of the X coordinate. The X -domain is reduced to the interval {-0.3, 0.3} for clarity. At both ends,X ± 0.5, ζ = 0 andX = ±0.5, ∂ X ζ ± 0.03. δ D = -D 0 e -(X -X d ) 2 /l 2d with D 0 = 0 for the black line and D 0 = 0.3 for the blue, red, magenta colored lines. Except for the magenta line, the dip is centered on X = 0; otherwise, it is located at X d = 0.3. The blue and magenta curves correspond to a narrow dip with l d = 0.1, the red curve to a larger one with l d = 0.3. Notice that a dip located far from the center has no effect and the deflection remains weak in any case (scale of the vertical line). In (b), ζ for a pre-stressed shell of equation W 0 (X ) = (1/4-x 2 )/5 of negative curvature and clamped boundary conditions: X ± 1/2,ζ = ∂ X ζ = 0. In (c), the full profile = ζ + W 0 with a dotted black line for W 0 . The same code of colors is applied in the three panels (a, b, c)

Fig. 5 Fig. 6

 56 Fig. 5 Buckling deviation ζ in a pre-stressed growing plate. The prescribed boundary conditions are: ζ = 0 and ∂ X ζ = ±0.03 at X ± 0.5, N 0 = 4 and the black line is the initial shape of the plate under pre-stress. The defect is as follows: for D , δ D (0) = -0.3e -(X -X d ) 2 /l 2 d , δS G = s g e -(X -X d ) 2 /l 2 d with s g = 0.5 (resp. -0.5) represented by a continuous blue line (resp. a dot-dashed blue line) and s g = 1 (resp. s g = -1) represented by a continuous red line (resp. dotted-dashed red line) for the 3 panels. In (a), X d = 0 and l d = 0.1. In (b),X d = 0 and l d = 0.3. In (c), the same as in (a) for an eccentric defect positioned at X d = 0.2. Notice the strong asymmetry of all curves in panel (c)

Fig. 7

 7 Fig.7 In (a, b, c) The shell profile (X ) = W 0 (X ) + ζ (X ). The colored lines correspond to the buckling conditions of Fig.6). Notice the possibility to observe an inversion of curvature in panel (a) and (b) for a suitable size of s g with respect to the curvature of the shell

Fig. 8

 8 Fig.[START_REF] Delon | The integrin adhesion complex changes its composition and function during morphogenesis of an epithelium[END_REF] Fitting experimental wing disc profiles with the uniaxial FvK bilayer theory. (a): Theoretical fit of an early wing imaginal disc, for which the initial profile before perturbation was flat (same as Fig.2c). (b): Theoretical fit of an older wing imaginal disc, for which the initial profile before perturbation was curved (same as Fig.2e)

8 Fig. 9 (

 89 Fig. 9 (a): Profile of the interface between the two layers for different curvatures and a pre-stress σ P S = 200 Pa, with no defect in the lower layer. The continuous, dashed, and dotted lines correspond, respectively, to W 0 + ζ , W 0 and ζ . The red, green, and blue lines correspond, respectively, to H = {-20, 0, 20}µm. (b): Profile of the interface between the two layers for a case with a defect of relative thickness h d / h s = 0.85 with h d the defect thickness and h s the bottom layer thickness. The growth parameter is g = 0.05. The green, blue, and red curves correspond, respectively, to cases with an initial zero, concave, and convex curvatures H = 0, 20, -20. (c): Different shifts for the defect, which is at X = 0, 0.08, 0.16, respectively, for the blue, green, and red lines. The amplitude of the deformation depends on the localization of the defect. The growth parameter takes the value g = 0.1. (d): Comparison between constrained (dashed lines) and free (continuous lines) boundary conditions for zero, convex, and concave initial curvatures (H = 0, -20, 20 for the green, red, and blue lines) and a growing bottom layer (g = 0.03). (e): Interface between the lower and upper layer in red. (f): Buckling for a constrained bilayer with a growing bottom layer and an initial concave curvature

Fig. 10

 10 Fig. 10 Cuts of the imaginal wing for the different FEM simulations with and without defects. The colors reflect the stiffness: red for the basal ECM, blue for the upper membrane and the cell layer, and cyan for the zone of the cell layer around the defect that softens in the refined model

Fig. 11

 11 Fig.11Fitting experimental wing disc profiles with the FEM model. (a): Theoretical fit of an early wing imaginal disc, for which the initial profile before perturbation was flat (same as Fig.2c). (b): Theoretical fit of an older wing imaginal disc, for which the initial profile before perturbation was curved (same as Fig.2e). The stiffness of the cell layer is not altered by the defect. (c): Comparison between the refined fit when the stiffness of the cell layer is altered by the defect (in green) and the case when it is not (red)

Acknowledgements

The authors acknowledge the support of ANR (Agence Nationale de la Recherche) under the contracts MecaTiss (ANR-17-CE30-0007) and GuideFusion (ANR-18-CE13-028); Excellence Initiative of Aix-Marseille University -A*Midex (capostromex), a French Investissements d'Avenir programme. PQQ acknowledges the support of the China Scholarship Council (CSC), file N 0 201706100182.

Author contributions

The two first authors contributed equally. Both LL and MBA conceived the project and wrote the manuscript.