Disentangling geometrical, viscoelastic and hyperelastic effects in force-displacement relationships of folded biological tissues

Drosophila wing discs show a remarkable variability when subject to mechanical perturbation. We developed a stretching bench that allows accurate measurements of instantaneous and time-dependent material behaviour of the disc as a whole, while determining the exact three-dimensional structure of the disc during stretching. Our experiments reveal force relaxation dynamics on timescales that are significant for development, along with a surprisingly nonlinear force-displacement relationship. Concurrently our imaging indicates that the disc is a highly heterogeneous tissue with a complex geometry. Using image-based 3D finite element modelling we are able to identify the contributions of size, shape and materials parameters to the measured force-displacement relations. In particular, we find that simulating the stretching of a disc with stiffness patterns in the extra-cellular matrix (ECM) recapitulates the experimentally found stretched geometries. In our simulations, linear hyperelasticity explains the measured nonlinearity to a surprising extent. To fully match the experimental force-displacement curves, we use an exponentially elastic material, which, when coupled to material relaxation also explains time-dependent experiments. Our simulations predict that as the disc develops, two counteracting effects, namely the discs foldedness and the hardening of the ECM lead to force-relative displacement curves that are nearly conserved during development.


Introduction
The role of mechanical forces in developing tissues has been attracting increased attention in recent years. This is due to technical advancements both in experimental techniques and computational methods. In this work, we focus on unravelling the interplay of material properties and complex shapes in living tissues.
Several approaches have been developed to measure the material parameters of living tissues [1]. Those include local perturbation experiments, such as atomic force microscopy [2], laser microsurgery [3][4][5][6], optical and magnetic tweezers [7,8], but also force inference methods [9][10][11] and FRET sensors [12]. Those approaches are very powerful to estimate local properties, but inherently require a certain previous knowledge of the tissue. Inferring large-scale material properties from local measurements is only feasible when the material is somewhat homogeneous. A different strategy is to measure the mechanical properties at the tissue scale, for example by stretching the tissue (see [13][14][15][16][17]). These measurements are still limited in the case of complex geometrical structures, such as those of a e-mail: aegerter@physik.uzh.ch the wing disc, but, combined with computational modelling, can be used to infer local mechanical properties.
The wing disc has been a paradigm to study the role of mechanical feedback on growth. Especially 3rd instar discs have been studied extensively [18][19][20][21]. The discs grow [22] and change shape [23] to a surprising extent during 3rd instar. Furthermore, 3rd instar discs can be dissected out of the larvae, which has permitted a number of ex vivo studies investigating the role of mechanical feedback [16,[24][25][26][27][28][29][30]. Yet, it is still unclear how the disc grows into the folded geometry of late 3rd instar. A crucial role in tissue shaping has been assigned to the extra-cellular matrix (ECM), the structural scaffold onto which cells attach, both in wing disc [31,32,23] and in other Drosophila tissues [33]. Recently, though, whether the ECM regulates growth through mechanical, rather than biochemical signals has been questioned [34], which calls for a direct experimental investigation of the mechanical properties of the ECM, since genetic studies alone cannot completely disentangle the roles of biochemistry, shape and mechanical properties.
Here, we focus on understanding how tissue scale mechanical properties and wing disc shape affect the disc response to large-scale mechanical perturbation. To do so, we used the same device as in [17], similar to that of [24], which allows the ex vivo stretching of live 3rd instar wing discs. Subsequently, we developed computer simulations to predict what the disc composition would have to be to reproduce the experimental force-displacement relationships. In particular, our focus is on disentangling the role of geometry and material properties. In a recent publication, it was shown that the wing disc exhibits an extreme strain-stiffening behaviour when stretched that could be explained by combining linear materials and the folded geometry of the disc in a 2D simulation of a cut of the disc [17]. Here, we develop 3D finite element simulations on an image-based wing disc geometry and show that linear materials can indeed account for a considerable, yet not sufficient amount of the nonlinearity when taking the full three-dimensional geometry into account. Our simulations, when exponential elasticity is assumed recapitulate the instantaneous material response, and can also account for time-dependent dynamics if force relaxation is included. Finally, we use simulations to predict how the tissue scale dynamics vary as the disc grows and folds up during 3rd instar. Our results call for a careful characterization of wing disc geometries over time and for mechanical perturbation experiments that address local material properties.

Stretching the wing disc highlights nonlinearity and relaxation
The wing disc is an extremely soft living tissue. To measure its mechanical properties, it must be kept intact and alive. The discs were dissected from late 3rd instar larvae and inserted in a stretching bench (see [17] for a drawing of the setup). The discs were submerged in wing medium 1 (WM1) [35]. The two dorso-ventral extremities of the disc were attached to two glass slides (see fig. 1(g) for nomenclature). The left was kept fixed, while the other one was attached to a cantilever controlled through a purposemade software. The software was developed to stretch the disc imposing a force and measuring a displacement, or vice versa. In particular, force or displacement trajectories were imposed that increased linearly, harmonically oscillated, or were constant in time. Furthermore, the disc is accessible to confocal imaging from the bottom, which allowed the recording of 3D image stacks. Figure 1(a) shows discs that are stretched following different modes of operation of the device. The blue and red dots are obtained controlling displacement of the disc, while the yellow dots controlling the force the disc is subject to. The blue dots show a disc that is very slowly stretched through imposing a linearly increasing displacement. Slow stretching of the disc allowed imaging with minimised disc vibration and we chose to extract the geometry of this disc (see below). The red and yellow dots correspond to discs that were harmonically stretched: the first elongation of a disc is shown by the leftmost dots -as the discs are repeatedly stretched the displacement at a given force value increases. The disc in yellow had already been stretched before the experiment plot in the figure.
Therefore, it seems that the discs stress-stiffen as displacement increases, and soften over time when stretched. We will use computer simulations to disentangle those two behaviour. Force relaxation is better addressed by fast stretching and then holding of the disc at constant displacement while looking at the force dampening over time ( fig. 1(b)).
Hence, discs show extremely nonlinear material behaviour and, on top of that a considerable amount of variability, best visible when discs are stretched at the same rate (2 μm/s, fig. 1(c)). Furthermore, the discs also have a complex geometry ( fig. 1(g)), seemingly made of non-homogeneous materials. We used flies expressing Trol::GFP and Lac::YFP, red and green, respectively, in fig. 1(g), with the former being expressed in the ECM and the latter at cell membranes through the disc [36,37]. An interesting aspect to note on the views from the top (left) is that the sides of the hinge do not straighten up. This suggests that they are not bearing a large part of the stress. When visualising the central disc section (right column), the disc is easily divided in three parts, from left to right: pouch, hinge and notum. The depth of the hinge folds has been shown to have a large impact on force-displacement curves in numerical simulations [17].
Hence, to carefully study the impact of the folds, we decided to create a very accurate initial geometry. To do so, we segment the disc when attached but not yet stretched, and derive a starting geometry to use in simulations ( fig. 1(h)). We define four regions of the disc whose foldedness we quantify, from left to right: pouch, ventral fold (V-fold), dorsal fold (D-fold) and notum. Foldedness is calculated in the central slice of the disc, as arclength divided by length of each region as bounded by the points marked by crosses in fig. 1(h). Geometries were created with MATLAB (MathWorks, Natwick, MA, United States), using custom made routines and the toolbox iso2mesh [38]. The custom made routines were used to manually segment the attachment slides and the disc on ≈ 15 anterior-posterior slices (one of those slices is plot on the right side of fig. 1(g)). The disc structure was reconstructed using an adapted version of the interpolation routine in the MathWorks File Exchange entry interpmask -interpolate (tween) logical masks [39]. The interpolated structure was used to create a voxelised image, whose voxels were assigned to distinct disc regions as follows ( fig. 1(h)). The inside and outside of the disc were defined using erosion. The inner part of the disc is taken up by a single region. The outer geometry is assumed to have a constant thickness of ≈ 3 μm and is subdivided in three regions: ECM on the basal side, apical side and lateral sides. In this work we focus on varying the material properties of the ECM with the respect to the rest of the disc. iso2mesh was used to create initial geometries, consisting of 1 to 3 million tetrahedral elements. The number of elements depends on the discretisation required to properly mesh each region. Simulations were carried out with the finite element framework FEBio [40]. Simulations were run on the CPU cluster Hydra of the University of Zurich. (a) Force-displacement curves for 3 different discs: blue dots correspond to a disc that is slowly stretched to improve quality of imaging and extracted geometry; red and yellow dots correspond to discs that are stretched in an oscillating fashion; red dots correspond to a disc that is extracted from an earlier 3rd instar larva compared to the other two. (b) Holding a disc at constant displacement shows force relaxation. (c) Force-displacement curves for 5 discs, all stretched at a rate of about 2 μm/s. (d)-(f) Force and displacement over time for the three discs of panel (a). (g) Images of a stretched disc for increasing force-displacement (F -d) values (re-drawn from [17]). (h) Initial simulation geometry derived from the top image of panel (g).
Red crosses indicate the reference points used for the foldedness calculations for the four regions of the disc named herein (see sect. 2.2).  fig. 2(d)).

Simulating wing disc stretching with linear materials explains some, but not all of the nonlinearity
We start by modelling all of the wing disc regions as Neo-Hookean materials, which reduces to linear elasticity for small strains, but supports large deformations (see [41,42]). In FEBio the material is parametrised by E, Young's modulus, and ν, Poisson's ratio [43]. We assume that the entire structure is at rest when it is not stretched. We define E ECM and E DISC to be Young's moduli of the ECM and of the rest of the structure, respectively ( fig. 2(a)). Poisson's ratio is fixed as ν = 0.3 throughout.
Our first attempt consists of extending the 2D model of [17] to 3D. We use material parameters corresponding to the optimal parameter values found in [17] by numerical simulations: E DISC = 5 kPa, ν = 0.3, E ECM /E DISC = 100. With 3D simulations, the force-displacement curve corresponding to that scenario is straight, and reaches much lower displacements than the disc ( fig. 2(c)). Reducing E DISC the displacement values become comparable but still the curves remain essentially straight ( fig. 2(c)). Hereafter, we choose E DISC = 0.5 kPa in order to match the nearly horizontal force-displacement curve at low displacements. We choose a smaller E DISC compared to the 2D model, but care should be taken in comparing the values directly (see below). Similar results can be obtained with E DISC = 1 kPa. Examining the top views of the disc under stretching it seems that the lateral disc sides do not always visibly straighten out, especially the lower side between d = 0 μm and d = 100 μm ( fig. 1(g)). This suggests that the lateral disc sides should not bear much stress, and we wondered whether that is also the case in our simulations. Figure 2(f) shows that stress accumulates on the sides of the pouch, which ultimately prevents its elongation.
This suggests that the ECM on the sides of the pouch should be softer than the rest of the ECM. The same ought to be true for the ECM on the sides of the hinge -because it would otherwise prevent the flattening of the folds observed experimentally. This would indicate that [17] correctly modelled the central part of the disc and also correctly raised the warning that 3D models would be needed to properly take into account a geometry that is more complicated than just replicating the central slice over the width of the disc. Limiting oneself to the twodimensional case and thus purely planar stresses apparently does not capture the entire picture in the complicated three-dimensional structure of the wing disc.
We therefore designed a geometry where only the central region of the ECM is harder than the rest of the disc ( fig. 2(b)). We define E SIDES to be Young's modulus of the lateral parts of the ECM (yellow in fig. 2(b)). Here, force-displacement curves show a surprising amount of nonlinearity ( fig. 2(d)) considering that all of the involved materials are actually linear. It is interesting to compare this scenario with that above, in which the sides of the disc, already straightened out, were bearing the force and hence producing a linear force-displacement curve. Here, even though the materials are still linear, a nonlinear force-displacement curve is obtained, which may therefore be a consequence of the folded geometry, similar to the results of [17], as well as the inhomogeneity of the material. Indeed, geometry snapshots as force is increased show that the ECM straightens out (purple in fig. 2(g)).
Quantifying the foldedness as a function of displacement shows that the pouch flattens out at much lower values of displacement when the ECM sides are soft ( fig. 2(e)), similar to what is found experimentally. Comparing the most stretched geometry of fig. 2(g) with the experimental data of fig. 1(g), though, it seems that neither folds flattening, nor notum elongation are as pronounced as in the experiment.

ECM stiffness patterns predict wing disc deformation
To test whether combining harder and softer ECM regions would result in better matching geometries, we designed additional geometries where either the ECM is softened under the notum (light blue region in fig. 3(a)) or the ECM is softened in the deepest part of the folds (light blue region in fig. 3(b)) or both. Notably, we assumed that the softer ECM regions are only marginally softer than the rest of the ECM and still considerably harder than the disc proper. When different from E ECM , Young's moduli of the ECM under the notum and inside the folds are denoted as E NOT UM , E F OLDS , respectively. For this section we fix E DISC = 0.5 kPa and E ECM = 250 kPa, as that value of E ECM gave the closest force-displacement relationship at high displacements among those that we tested ( fig. 2(d)). Figure 3(c) shows that when part of the ECM is softened, the discs reach higher displacements for any given force, in a way that leads to reduced nonlinearity. We notice, though, that the softening factors affect in one case the relative elongation of the notum with the respect to the rest of the disc ( fig. 3(a)), in the other the extent to which the folds flatten out ( fig. 3(b)). If the foldedness is measured as a function of displacement, softening the notum also affects the foldedness of the D-fold, since the corresponding increased notum elongation shifts the entire curve to the right ( fig. 3(e)). Plotting foldedness against force in contrast, this effect is removed and the two parameters can be seen to affect the folding and the notum elongation independently ( fig. 3(b)). The effect is also visible when the V-fold is considered instead of the D-fold (data not shown). The same is true for the effect on the notum of softening the folds, but it is much less visible due to the much smaller size of the softened regions in this scenario.
Hence, we have shown that combining ECM regions of different stiffness may serve to find the best matching geometries, as indicated by comparing the foldedness of the simulation to the experimental values. This softening does however reduce the nonlinearity of the curves because of the decreased materials inhomogeneity overall. In practice, this means that the time evolution of the geometry in simulations can be fine tuned to match both the force-displacement curves and observed geometrical deformations by a specific stiffness map of different ECM regions.

Nonlinear material behaviour with relaxation fits experiments
Since combining linear materials does not explain all of the nonlinearity, we implement exponential elasticity to match the force-displacement curves. More specifically, we use a Holmes-Mow material [44,45]. In FEBio the Holmes-Mow materials is parametrised by E, Young's modulus, ν, Poisson's ratio and β, the exponential stiffening coefficient [43].
Interestingly, fig. 4(a) shows that varying E ECM and the exponential stiffening parameter, β, the nonlinearity of the force-displacement curve of the segmented disc can be achieved. For this specific disc, β = 5 provides a reasonable amount of nonlinearity. For some more nonlinear discs (e.g., the disc in purple in fig. 1(a)) a higher β = 10 could be used to increase the nonlinearity of the curves. Varying E ECM affects not only the slope of the curve at high force, but also the displacement value at which the slope of the curve increases most. In absence of relaxation, good agreement with experimental data is obtained by  fig. 4(a)). Linear material parameters: E ECM = 250 kPa (cf. fig. 2(d)). (c) Computer-generated geometry with flattened ECMside. (d) Computer-generated geometry with increased basal folding.
E ECM = 100 kPa. Also the geometry of stretched discs resembles that of experimental data ( fig. 4(b)), even though, better matching could be obtained by softening specific areas. As in the linear case, softening the ECM under the notum and folds leads to reduced nonlinearity of the curves ( fig. 4(c)) and to increased unfolding specifically in those regions ( fig. 4(d)). Indeed, a rough quantification of foldedness in experimental data gives a value of about 1.03 at F = 5 μN, which is in line with the scenarios with softener ECM inside the folds, however given the experimental uncertainty, the differences of these scenarios cannot be properly resolved. So far, we have disregarded the dynamics in time of stretching and material behaviour. To figure out the rate of material relaxation we consider a disc that is elongated by 150 μm in 10 seconds time and then held fixed ( fig. 1(b)). The force at time 10 s is ≈ 12 μN, while at time ≈ 1000 s it drops to ≈ 4 μN. To implement material relaxation we used the compressible viscoelastic implementation of FEBio [43], where a discrete relaxation spectrum is used, so that the relaxation function G reads [46]. Throughout this study we assume N = 1 and normalise g i so that Young's modulus of a purely elastic material and its viscoelastic equivalent with τ i = ∞ are equal. Hence, the relaxation function reads G(t) = 1 1+g (1 + g exp(−t/τ )). We assume that only the ECM is viscoelastic and choose E ECM = 100 kPa to match the near-instantaneous force and set the force to decrease by at most 2/3 of its maximal value (i.e. g = 2). Simulations predict the material relaxation of the disc ECM to be close to 100 s ( fig. 4(e)). It also seems that adding a second, slower relaxation timescale would provide better fitting. To avoid adding more parameters and in the interest of simplicity, we use g = 2 and τ = 100 s in subsequent simulations.
The slowly stretched disc was stretched at a rate of 0.05 μm/s. We implemented that stretching rate and the relaxation timescale of ≈ 100 s found above. Figure 4(f) shows that the force-displacement curve is very well matched by a nonlinear material with E ECM = 300 kPa. Note that this is a different value to the best matching one of fig. 4(a). This indicates that relaxation shifts displacement value at which the slope of the curve increases most. Furthermore, the strain-stiffening characteristic of Holmes-Mow materials leads to force-displacement curves that do not decrease in slope for high displacements when relaxation is introduced (i.e. τ = ∞ vs. τ = 100 s). The closest a linear material can get is well exemplified by E ECM = 750 kPa: the final displacement can be matched but not enough nonlinearity is preserved with material relaxation.

Predicting the source of intra-disc variability: size, foldedness, material parameters
Having shown we can accurately reproduce specific experiments we next focus on explaining the variability between different discs. In this context there are different aspects to keep in mind. First of all, the discs more than double in size during third instar [22]. In order to be able to compare their force-displacement curves, displacements should be normalised. Since the stress is largely borne by the central part of the ECM, we choose the minimal distance between attached disc extremities in the central slice as normalising factor. We do not normalise force since it is not known how the width of the hard ECM scales with disc size. When force is plot against the relative displacement, the curves are much closer to each other ( fig. 5(a)). Yet, they do not perfectly overlap. We explore two more sources of variation.
Firstly, the discs are known to increase their foldedness during third instar [23]. Since it is experimentally not feasible to vary the foldedness of a disc while leaving its size and material parameters unchanged, we designed two geometries that evaluate the impact of an ECM with more or less pronounced folds ( fig. 5(c), (d)). Interestingly, our simulations predict that the foldedness severely affects the force-displacement curves. Indeed, as foldedness increases, the amount of displacement for a given force value increases, and it does so increasing the difference between the initial and final slope of the curve ( fig. 5(b)).
Secondly, it is reasonable to assume that the ECM increases in thickness and thus stiffness during 3rd instar, since the collagen that makes up the ECM is continuously secreted by the fat body adjacent the disc and incorporated in the ECM [32]. As already shown in fig. 4, increasing the ECM stiffness leads to decreasing the displacement at a given force value. This means that as the discs grow, two effects counter-act each other. On the one hand, stiffening the ECM shifts the nonlinearity in the force-displacement curves to smaller displacements, since in that case considerable stresses are already incurred in the ECM at small displacements. On the other hand, increased folding of the ECM shifts the nonlinearity to larger displacements, corresponding to larger stretching because of the unfolding of the deeper fold-depths.

Discussion
The wing disc during 3rd instar larval development has attracted attention for the opportunity to study the role of mechanical forces during its development. To study the role of mechanics in shaping tissues, though, it is crucial to understand the material properties and the composition of the disc. To do this, we designed a stretching bench that we use to extend the disc. The disc as a whole exhibits nonlinear material behaviour at the tissue scale. Whether this nonlinearity is caused by nonlinear material behaviour or its folded geometry is unknown. The ECM on the basal side of the disc is thought to be the main mechanical contributor [17].
We set off to separate the contributions of material parameters and shape by designing an initial geometry from an imaged disc. Firstly, our simulations predict that only the ECM in the central region of the disc is hard. Otherwise, the pouch and the folds would not be straightening out. When the ECM is modelled as a linear hyperelastic material, a considerable amount of nonlinearity emerges, yet not enough to match the experimental curves. With linear materials, though, it is already possible to closely reproduce the unfolding of the ECM by assuming specific stiffness patterns of the ECM.
To explain the rest of the nonlinearity, we turn to a nonlinear material. In vitro experiments with ECM-like collagen have indeed shown that the material exhibits considerable strain-stiffening [47]. As a nonlinear material, we chose an exponentially stiffening material, Holmes-Mow. With this we can match the segmented disc. Furthermore, when the material is modelled as exponentially elastic, also time-dependent force relaxation dynamics can be closely reproduced. We modelled a material behaviour that is independent of biochemical pathways. In view of recent evidence showing what pathways respond to stretching [30], it would be very interesting to couple the local material behaviour to biochemical dynamics in the wing disc.
Finally, we use our simulations to predict how the features of the disc change during 3rd instar development. The discs are known to increase in size [22] and foldedness [23] during 3rd instar, as well as in ECM stiffness [32]. Normalising force-displacement curves according to the disc size, we are able to make predictions on the role of foldedness and ECM stiffness during development. In particular, our simulations predict that foldedness and ECM stiffness have opposite effects that balance out on forcedisplacement curves.
Financial support was provided by SystemsX.ch (Morpho-genetiX initiative) and the Swiss National Science Foundation (SNSF). The authors appreciate the support of the Service and Support for Science IT unit (S3IT) of UZH.

Author contribution statement
FL carried out experimental work. FA extracted initial geometries and performed simulations. CMA designed and supervised research. All authors wrote the manuscript.
Publisher's Note The EPJ Publishers remain neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access
This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.