Spontaneous buckling of contractile poroelastic actomyosin sheets

Shape transitions in developing organisms can be driven by active stresses, notably, active contractility generated by myosin motors. The mechanisms generating tissue folding are typically studied in epithelia. There, the interaction between cells is also coupled to an elastic substrate, presenting a major difficulty for studying contraction induced folding. Here we study the contraction and buckling of active, initially homogeneous, thin elastic actomyosin networks isolated from bounding surfaces. The network behaves as a poroelastic material, where a flow of fluid is generated during contraction. Contraction starts at the system boundaries, proceeds into the bulk, and eventually leads to spontaneous buckling of the sheet at the periphery. The buckling instability resulted from system self-organization and from the spontaneous emergence of density gradients driven by the active contractility. The buckling wavelength increases linearly with sheet thickness. Our system offers a well-controlled way to study mechanically induced, spontaneous shape transitions in active matter.

C ontracting materials have many potential applications related to artificial muscles, sensing of mechanical stimuli, and active self-assembly of three-dimensional structures 1,2 . In particular, two-dimensional intrinsically contractile sheets may find applications as actively functional band-aids that speed up wound healing or in form of beating patches that could support heart function. Artificial contractile materials are often based on synthetic molecular motors 3 or use liquid crystals 1 that require external cues such as electrical stimuli, light, or temperature for activation and control. However, tremendous challenges must be overcome before synthetic active materials can be used for human health. An alternative approach is to employ reconstituted biological materials, notably viscous or elastic actomyosin networks similar to those found in living tissues. In addition to being obviously biocompatible, such systems demonstrate spontaneous contractility and self-organization [4][5][6][7] . The characteristics of contraction can be readily controlled by changing global parameters, such as the concentration of molecular motors or of passive cross-linkers [4][5][6][7] . In addition, contraction patterns and dynamics can be regulated by spatially patterned activation of myosin motors 8,9 . Contracting materials are part of essentially all biological organisms. They allow cells to sense mechanical properties of the environment 10 or enable deformations that can be used to deform the organism, which could be exploited for morphogenetic processes during development 11 . In these cases, buckling has been observed when a contractile epithelium was subjected to anisotropic contraction 12 and/or coupled to an elastic substrate 13 . How intrinsic contractions in absence of coupling to a substrate can intrinsically generate three-dimensional structures, however, has not yet been explored.
In fact, in developing biological systems, the formation of folds and other three-dimensional structures from two-dimensional objects is often attributed to in-plane growth [14][15][16][17][18][19][20] . However, there are systems, where structure formation occurs in absence of growth, for example in hornbean leaves 21 or during gastrulation of animal embryos 11,13,22,23 .
These examples suggest that spontaneously contractile materials could promote out-of-plane deformations of sheets. Currently, however, we lack good model systems for fundamental studies of such phenomena. A notable exception is provided by thin, elastic chemically crosslinked gel sheets that show largescale buckling upon non-uniform lateral contraction 24 . However, in this case, spatial variations in contraction result from preimposed gradients in the cross-linker density. Furthermore, in many cases-naturally occurring or artificial-contractile sheets are coupled to a substrate or other elastic structures 8,9,11,[20][21][22][23] , presenting an additional difficulty for studying contraction induced folding.
Here, we demonstrate the fabrication and dynamical properties of initially homogenous, thin, suspended, elastic acto-myosin sheets of controllable extent and elastic properties. These gel sheets spontaneously contract through the activity of myosin motors with no need for external stimuli save for the presence of ATP in solution. Even more remarkable, the gel sheets can selforganize into folded structures: buckles appear spontaneously in the absence of externally applied mechanical constraints and without imposed spatial variations of the gel composition or mechanical properties. The characteristic length scale of the folded pattern is proportional to the gel thickness and can thus be easily controlled. By analyzing the solvent flow during the gel contraction dynamics, we show that the acto-myosin gel behaves as a poroelastic active material. Using a theory of actively contractile, poroelastic gels, extended to include the catch-bond behavior of myosin motors, we explain the observations and the physical basis of the contraction dynamics. Overall, our results show that buckling can be spontaneously generated by motor activity and does not require mechanical coupling to the environment or pre-imposed gradients in sheet material properties.

Results
Formation of thin contractile actin sheets. Internally driven contractile gels were generated by polymerizing actin in the presence of the strong passive cross-linker fascin and clusters of myosin II motors 5,25 (Fig. 1a, Methods). Only if actin polymerization occurred in the presence of motors and cross-linkers, contractile elastic gels were produced where the motors were embedded in the network 5,7,25 . To generate thin sheets, we simultaneously introduced actin monomers, myosin, and fascin in chambers with a large aspect ratio ffiffiffiffiffi ffi A 0 p =h, where A 0 is the chamber area in the x-and y-directions and h the height in the z-direction (Fig. 1a). The chamber surfaces were treated to avoid sticking of the actin gel (Methods). After adding a drop of the protein mixture, the chamber was sealed with a cover glass. The chambers were only partially filled to avoid sticking of the gel to the lateral boundaries. In this way, a circular symmetric drop was obtained. Its size depended on the chamber height, which varied between 70 and 250 μm, and the drop volume, which was between 0.3 and 7 μL. Typically, the drops had initial radii R of 0.12-0.5 cm.
After mixing 5 M G-actin with 16.7 nM myosin motors and 280 nM fascin, defining the time t = 0, the system evolved in three phases (Fig. 1b, Supplementary Movie 1). In the first phase, which lasted 4 min, actin filaments nucleated and polymerized spontaneously. The protein distribution remained homogenous down to scales of several μm as one can infer directly from the fluorescence images and from the actin and myosin pair correlation functions G(r) (Fig. 1b, c, Supplementary Figure 1a-d, Methods). Note that concomitant with actin polymerization, the distribution of motors, while remaining homogenous, coarsens (Supplementary Figure 1e, f). In the second phase, which lasted 1 min, the system evolved into an interconnected network of actin bundles if both, fascin and myosin clusters, were present 5,25 (Fig. 1a, b). The network remained macroscopically homogenous and isotropic (Fig. 1b, c gray dots) with a mean mesh size ξ 0 = 60 ± 16 μm (mean ± SD, N meshes = 75; N gels = 1) at the end of this phase (Fig. 1d). This property is also reflected by the pair correlation function G(r) (Supplementary Figure 2). After this phase, the myosin motors dominantly localized at intersection points of filament bundles (Fig. 1a). In the third phase, the gel started to contract (Fig. 1b, c). The existence of this phase required prior formation of filament bundles that are mechanically sufficiently stable and a connection of these bundles by myosin motor clusters 5,7,25 . Consequently, the velocity of myosin clusters was the same as the local gel displacement velocity (Supplementary Figure 3, Supplementary Movie 2). Note that the three phases were probably not mutually exclusive as network reorganization could still have been accompanied by actin polymerization.
Contraction was isotropic throughout the process. This is in line with observations made on contracting actin gels that are attached at the periphery, where motors have been activated by light in specific areas and where it was shown that the final state presents the same aspect ratio as the initial state 9 . In the course of contraction, the radius r, which is determined from the projected gel area A gel , r ¼ The contraction phase can be divided into two sub-phases: in the beginning, the edge velocity v ¼ dr dt increased approximately linearly with time until it reached its maximal value v max at time t max (Fig. 1g). For times later than t max , the edge velocity decayed to zero exponentially with a characteristic time τ = 20 s. As particle image velocimetry (PIV) analysis shows, contraction always began at the gel boundaries (Fig. 1e, Supplementary Movie 3). The local gel contraction velocity decreased linearly from the gel boundary towards the gel center (Fig. 1h). During contraction, the network mesh size decreased (Fig. 1i), whereas the filament bundles remained straight, as confirmed by comparing the end-to-end distance l end-to-end and contour length l cont of filament bundles (Supplementary Figure 4c). This is different from contracting, dilute actin networks, where the weak cross-linker α-actinin was used instead of the strong cross-linker fascin and where individual actin filaments buckled 26,27 . Since our network started contracting from the boundary, the mesh size increased from the boundary to the center of the gel corresponding to decreasing gel volume fraction in that direction (Supplementary Figure 4). For an initially homogenous sheet, this contraction mode resulted in an increase of the gel density first at the periphery (Supplementary Figure 5f-i) and only later propagated into the gel interior (Supplementary Figure 5j). The density increase occurred gradually until it saturated at late times (Supplementary Figure 5j). In the stationary state, the contractile forces and the gel elasticity were equilibrated and the gel was homogenous again.
Contracting actin sheets can spontaneously buckle. We now turn to the dynamics along the short axis, the vertical z-direction. Similar to its lateral contraction, the gel contracted from an initial thickness e of 160 µm (the chamber height) to about 50 µm (Fig. 2a, Supplementary Movie 4). The vertical contraction began  Figure 6). Lateral contraction started only after the vertical contraction was essentially completed: during lateral contraction, the gel thickness e decreased slightly from 57 to 52 µm (Fig. 2b). This corresponds to an effective Poisson ratio ν ≈ 0. At advanced stages of planar contraction, the gel spontaneously started to buckle at about 580 s which corresponds to 1.5τ after t max (Supplementary Figure 6), as is evident from the developing bright and dark stripes that emerge perpendicularly to the boundary (Fig. 2c, Supplementary Figure 7, Supplementary Movie 5). The emerging folds were directed perpendicularly to the boundary (Fig. 2d, white arrows) and had a roughly sinusoidal shape (Fig. 2d, e, Supplementary Movie 6). We characterize this structure by means of "wavelengths", i.e., the distances between adjacent intensity maxima and adjacent intensity minima. The steady state thickness of that gel was b = 23.5 ± 3.4 μm (mean ± SD, N thicknesses = 69 for N gels = 1, Methods) and was uniform throughout the gel (Fig. 2e). We determined the average wavelength of the buckled state by measuring the distances between adjacent folds. To this end we used a confocal slice and analyzed the fluorescence intensity along a line tracing the border of the gel (Supplementary Figure 8, Methods). Analyzing the gel profile in the xz and yz planes close to the gel boundary gave very similar values (Fig. 2d, e). The distribution of wavelengths of the folds had an average value of λ = 196 ± 33 μm (mean ± SD, N wavelengths = 24 for N gels = 1. Supplementary Figure 8). In two-dimensional cuts across the xy plane, the folds appeared as areas of increased and decreased fluorescence intensity (Fig. 2e). The resulting folds were visible as straight bright strips at low magnification (Supplementary Figure 9). The contracting actomyosin sheet is a poroelastic gel. To obtain a better understanding of the system dynamics, we now characterize the material properties of the system taking into account the actin gel, the myosin motors, and the aqueous phase within the gel. It has been proposed that actin gels in living cells behave as poroelastic materials [28][29][30] . There are two distinguishing features of such materials 28 : First, the generation of a flow of the penetrating solvent caused by active myosin contractility. Second, diffusive stress relaxation characterized by an effective diffusion constant that depends on the drained elastic modulus, the solvent viscosity, and the gel porosity. We first checked whether myosin contractility generated a solvent flow. To measure the solvent flow, we added fluorescent beads to the solution. We used beads of 2.3 μm in diameter during the initial and intermediate phases of contraction, when the average pore size was larger than 15 μm, and beads of 200 nm diameter in the later contraction phase, when the pore size was smaller. The fluorescent beads, on average, moved in the outward radial direction as the gel contracted toward its center ( Fig. 3a, b, Supplementary Movie 7). The beads' radial velocity v r was up to 20 times larger than the gel edge velocity as long as they were inside the gel (Fig. 3c, d, Supplementary Figure 10, Methods). As the beads moved outwards from the gel center, their velocities initially increased dramatically by more than an order of magnitude. This was consistent with the gel contraction velocity profile, which also increased from the cell center to the periphery (Fig. 1h). The beads slowed down as they approached the boundary, where the gel densified. The filled circles in Fig. 3d indicate the time, when the beads left the gel. Outside the gel, the beads still moved for some time. Note that this movement was not due to inertia, because the Reynold number is less than 2 × 10 −4 even for the largest bead velocity of 100 μm s −1 and bead size of 2.3 μm. For later times, the average velocity of the beads decreased and fluctuated (Fig. 3d).
We attributed the fluctuations in the bead velocities to the porous structure of the actin network. To test this hypothesis, we investigated with higher resolution bead trajectories towards the end of the contraction process, when all movements had slowed down. Indeed, the trajectories of the beads were tortuous (Fig. 4a, Supplementary Figure 11) and their velocities fluctuated between 2 and 10 µm s −1 (Fig. 4b). The bead velocities increased with their distance from the actin bundles (Supplementary Figure 12).

Actin488
Nile red beads x (μm) Hence, the beads moved fastest in the center of a pore and slowest, when they were close to a bundle reflecting the local structure of the gel. In this advanced phase of gel contraction, the local solvent speed was on average faster than the local gel contraction velocity by a factor of 2.5 ± 0.7. The gel velocitysolvent velocity correlation function shows that locally, the fluid flow was always opposite to the gel flow (Fig. 4c, Methods). Again, the bead speed increased from the gel center towards the  Figure 13). All this indicates that as the gel contracted, the solvent was driven out as expected for a poroelastic material.
We now turn to the second property of poroelastic systems, namely that stress relaxation is determined by an effective poroelastic diffusion constant D. Explicitly, the relaxation time scales as r 2 max =D, where r max is the system radius at the beginning of the relaxation process, that is, at the time when the gel contraction velocity has reached its maximum. The value of D is proportional to an effective, concentration-dependent elastic modulus κ of the gel and inversely proportional to an effective friction constant γ that accounts for the permeation of water through the actin network, such that D $ κ=γ (Supplementary Notes). The friction coefficient γ scales as η=ξ 2 , where η is the solvent viscosity and ξ scales with the mesh size of the gel-i.e., the distance between crosslinks in the actin gel. The elastic modulus has units of energy per unit volume and is inversely proportional to the volume of a pore, whereas the friction coefficient γ depends on the pore facet perpendicular to the solvent motion. For the actin network, the pore size can be approximated by the mesh size. Initially, the gel is isotropic. Yet, at the time the lateral contraction velocity reached its maximum, the contraction in the z-direction was already essentially completed as mentioned above (Fig. 2b). Thus, the distance between two crosslinks of filament bundles is ξ k within the contraction plane and ξ ? perpendicular to it, consequently, κ $ 1=ξ 2 k ξ ? and γ $ η=ξ jj ξ ? 31,32 . Again, we take the values of these quantities at the time when the contraction velocity has reached its maximum, so that ξ k $ r max . Altogether, we thus have τ $ ηr 3 max . To test the dependence of the relaxation time on the solvent viscosity, we changed the latter by adding various amounts of glycerol. We conducted three experiments for each amount of glycerol. With increasing values of the viscosity, the duration of the polymerization and of the reorganization phases were extended (Supplementary Figure 14). In contrast, the linear acceleration and maximal contraction velocity did not change significantly (Supplementary Figure 14), suggesting that network organization during the polymerization and the reorganization phases adapts to changes in the solvent viscosity. The measured relaxation times obeyed the scaling relation derived above (Fig. 4d) further confirming the poroelastic nature of our system.
Having shown that the system is poroelastic, we can formulate a physical description of the contraction dynamics within the framework of active generalized hydrodynamics 33,34 . The system state is given in terms of the gel and solvent densities. Their time evolution is determined by the continuity equations that describe mass conservation of the two phases (Supplementary Notes). The respective gel and solvent velocity fields that appear in the continuity equations are obtained by balancing the forces due to gradients in the total mechanical stress in the gel with the friction forces resulting from motion of the actin gel relative to the surrounding fluid: Here, u denotes the displacement field of the actin gel and v is the solvent's velocity field. Furthermore, σ el is the drained elastic stress of the gel, which depends on the drained bulk and shear elastic constants K and μ, respectively (Supplementary Notes). The embedded motors generate an active contractile stress σ act . Within the framework of active generalized hydrodynamics it is expressed as −QζΔμII 35 . In this expression, Δμ is the difference between the chemical potentials of ATP and its hydrolysis products, which drives the contraction process. The factor ζ is a phenomenological constant linking the active processes to the causing thermodynamic force Δμ. In general, ζ depends on the filament and motor densities. We take it to be constant and account for the density of crosslinking motors through the fraction Q of bound motor molecules. Equivalently, the magnitude of the active stress is determined by the active force dipoles, generated by the molecular motors 36 , and is given by the force dipole density. We assume the active stress to be macroscopically isotropic and thus multiply its expression by the identity operator II. Equation (1) together with mass conservation (gel plus solvent) and the Stokes equation for the solvent describe the contraction of an active poroelastic gel (Supplementary Notes). This approach is simpler than those taken in refs. 8,37 , where a nonlinear dependence of the stress on the filament or myosin densities, respectively, were considered for contracting cytoskeletal networks. In ref. 9 , a static spring network without shear elasticity was used to describe contraction. Furthermore, all these works neglected the solvent dynamics, which we have shown above to be important at least in our system. We next hypothesized that the acceleration phase of the contraction dynamics was due to changes in the fraction of bound myosin motors that increases with the mechanical stresses in the gel. In terms of force dipole stresses, this means that the force dipole moment changes. Indeed, myosin forms catch bonds 38 that are strengthened in a macroscopically contracting network, because the actin filaments in-between the contractile myosin clusters (that act on filaments of opposite orientation) are stretched during network contraction. The fraction of bound motors Q is thus determined by The rates k on and k off , respectively, denote the myosin attachment and detachment rates. Whereas we take k on to be independent of the elastic stress, the detachment rate k off decreases with increasing elastic stress of the actin network σ el reflecting the motor's catch-bond characteristic according to Bell's law 39

(Supplementary Notes).
We solved the dynamic equations numerically using a finite difference scheme (Supplementary Notes). We use the experimental parameter values estimated from the relaxation dynamics and the buckling of the contractile sheets, see below. Notably, we use an effective elastic modulus of 0.05 Pa. For determining the friction coefficient between the solvent and the gel, we use a mesh size ξ || , ξ ⊥ = 1200 μm 2 taken from the experiment in Fig. 1b at t max , where we took ξ ⊥ = ξ || /3 (taken from experiments, see below), and an effective viscosity of 0.05 Pa·s. This value accounts for molecules dissolved in the penetrating solvent that increase the viscosity with respect to pure water. In our numerical calculations, we mimicked the polymerization and reorganization phases by increasing the fraction Q of bound motors from 0 to Q 0 linearly in time. During this phase, Q is kept spatially homogenous. We took Q 0 = 0.15 as estimated in ref. 5 . Afterwards, the value of Q evolved according to Eq. (2).
Our numerical solution of the dynamic equations for an initially homogenous, circular symmetric, elastic gel disc with an initial radius of 1.5 mm represented a radial contraction of the gel that starts from the boundaries (Fig. 4e). Contraction from the boundary was a consequence of the initially homogenous distribution of motor-induced force dipoles, which generated a homogenous active stress, such that active forces canceled each other everywhere but at the boundaries. The gel volume fraction was always highest at the boundary and decreased toward the gel center, until, in the final state, the gel density was again homogenous throughout the gel (Supplementary Figure 15). In the final state, the gel radius was 150 μm corresponding to a final lateral strain of 90%. Similar to the gel volume fraction, the elastic stress was also highest at the gel boundary but was uniform toward the end of the contraction process (Supplementary Figure 16). In the final state, mechanical equilibrium was reached as the effective stress diffusion current vanishes and the uniform elastic stress balances the uniform active stress. During the contraction process, the solvent was squeezed out of the gel at a velocity that increased with the distance from the gel center (Fig. 4e). The edge velocity of the contracting gel increases approximately linearly during the initial polymerization and reorganization phases, where the fraction of bound myosin Q was increased but remained approximately spatially homogenous (Fig. 4f). Subsequently, contraction accelerated dramatically until it reached a maximal velocity of about 20 μm s −1 and then relaxed exponentially in time toward 0 (Fig. 4f). The acceleration concurred with a strong increase in the fraction of bound motors, which is due to the catch-bond behavior (Fig. 4f). Throughout the contraction process, the fraction of bound motors remained essentially spatially homogenous (Supplementary Figure 17). This is due to the fact that the motor binding and unbinding rates are large compared to the time-scale of contraction and because we chose the motor unbinding rate to react very sensitively to stresses. The contraction velocity profile changed between the acceleration and the relaxation phases (Supplementary Figure 18); during the acceleration phase it was concave. It exhibited a convex part after the maximal contraction velocity was reached. The profile of the solvent velocity showed a similar change, however, it was delayed with respect to the change of the gel's contraction velocity profile. In the late stages of the contraction process, the fluid speed exceeded the gel contraction speed (Supplementary Figure 18) similarly to what was measured experimentally (Figs. 3d and 4b, Supplementary Figures 10, 12). The good quantitative agreement between the experiments and the numerical results points to the physics underlying the observed behavior: contraction starts from the periphery as the active forces are initially imbalanced only at the boundaries, the contraction velocity increases due to the increased stress in the gel and the motors' catch-bond characteristics, and poroelastic diffusion of the elastic stress yields eventually a uniform state of mechanical equilibrium.
Next, we looked at how the contraction of the gel propagated towards the gel center. To this end we chose to follow a value of the gel density as it propagates from the gel boundary inwards. The distance d traveled by the density in the reference frame of the laboratory shows different time dependence before and after t max (Fig. 4g). Before t max , it is concave showing that the velocity at which it moves inwards increases. After t max , the curve is convex corresponding to a slowing down of this velocity. For higher values of the density, the time to reach the center increases even though the distance that needs to be traversed decreases. In the numerical solutions, we find the same dependence of the time needed to reach the center on the initial distance to the center as in the experiments. However, the dependence of the distance travelled on time is always convex in the numerical solutions suggesting that nonlinear elastic properties that are not captured by our description affect the dynamics of contracting actin gels (Fig. 4h).
Estimation of the buckling wavelength. Having established the poroelastic nature of our reconstituted actin gels, we further investigated the buckling instability. In Fig. 5, we present the contraction (Fig. 5a-e) and the spatiotemporal evolution of the density (Fig. 5f-j) of a gel that buckles. Compared to gels that did not buckle (Supplementary Figure 5), the increase of the gel density at the periphery was relatively larger and persisted in steady state. These two situations differed furthermore in that the density gradients in buckling gels were steeper than in nonbuckling gels. In both cases, the density increased first at the periphery (Fig. 5a, b, f, g). In the case of buckling, at some point, the density at the periphery saturated, whereas the bulk density continued to increase (Fig. 5h-j, Supplementary Movie 9). In this case, contraction proceeded apparently in presence of a fixed perimeter, which is only possible if the sheet buckles. This is reminiscent of minimal surfaces under the constraint of a fixed perimeter that is longer than the perimeter of a circle of the same surface area. In very rare cases, we observed under the same conditions as in Fig. 5 that the gel ruptured instead of buckling (Supplementary Figure 19, Supplementary Movie 10). In these cases, the gel periphery became much denser than the bulk and was apparently much stiffer than in the buckling cases. The exact conditions responsible for this behavior remain to be explored.
How does buckling depend on the system's mechanical parameters? To answer this question, we estimated these parameters and expressed the buckling wavelength in terms of these. We can estimate the effective drained gel elastic modulus κ ¼ 4μ 3 þ K from the relaxation time τ by κ $ γr 2 max =τ, where r max is the radius at t max (Supplementary Notes). Using the above expression for γ in terms of the solvent viscosity η and the pore sizes ξ || and ξ ⊥ at t max , we obtain where R and ξ 0 are the initial radius and pore size, respectively, and b is the final gel thickness. The active stress can be estimated from force balance at steady state, where QζΔμ ¼ κε, whereε is the lateral strain at the end of the planar contraction phase relative to the extension of the sheet at t = t max (Methods). We considered the deformation after t max , because for that phase the exponential decay strongly suggests that κ is constant, which is implicit in the last scaling relation. We can now estimate the wavelength λ c at buckling onset. For a spatially homogenous activity, the contracting gel is analogous to a system under a compressional external force (Supplementary Notes). For a given external force there exists a minimal extension for which an elastic sheet would buckle (Supplementary Notes). This length was taken as an estimate for the buckling wavelength λ c . In this way, we get σ act = QζΔμ~κ(2πb/λ c ) 2 . With the above expression for the active stress and the effective elastic modulus, we finally obtain After having established a relation between the buckling wavelength and the final lateral strain, we checked its validity. We started by determining the final thicknesses of gels assembled in chambers of different heights and imaged with confocal microscopy (Fig. 6). Qualitatively, we observed an increase of the final sheet thickness and of the buckling wavelength with the chamber height h (Fig. 6). For each chamber height, we investigated 3-4 gels, and found a linear dependence of the final gel thickness b on the chamber height h (Fig. 7a). In a range of chamber heights between 80 to 160 μm, b ≈ h/3. As a function of the final sheet thickness b, the final strainε varied by less than 10% (Fig. 7b). Similarly, the active stress and the effective elastic modulus κ did not depend significantly on the gel thickness ( Fig. 7c) showing that in the range of system heights used these are intrinsic material properties that do not depend on the system geometry. Explicitly, we found κ ≅ 0.02 Pa. In contrast, an elastic modulus of 10 kPa has been measured for cells 40 . The difference between this value and the one found for our gel is due to the much larger distance between crosslinks in our system (50 μm at t max ) compared to 200 nm in cells 41 . The ratio between the cellular elastic modulus and that of our gel is given by the volume ratio ξ 3 gel =ξ 3 cell , where ξ cell and ξ gel are the distance between crosslinks in cells and gels, respectively 31 . This explains the factor of 10 6 between the two moduli. Finally, comparison of the experimental values of the buckling wavelength at steady state λ with the values obtained from the scaling relation (Eq. 4) is remarkably good (Fig. 7d). Sinceε is essentially constant, the wavelength is proportional to the final gel thickness. The value of the strain averaged over all heights,ε ¼ 0:8 ± 0:02 (mean ± SD, N gels = 13), yields 2π= ffif ε p ¼ 7:0 ± 0:1, which agrees very well with the measurements (Fig. 7d). Beyond a gel thickness of 80 μm we did not observe buckling of the gel edges.
As we have shown, initially homogenous elastic sheets contracting under the influence of molecular motors are wellcontrolled systems to study mechanically induced spontaneous shape transitions in active matter. We found that the network behaves as a poroelastic material where a flow of fluid is generated during contraction. The buckling instability resulted from system self-organization and relied on spontaneous formation of density gradients driven by active contractility.
The detailed conditions for observing buckling still need to be explored. In particular it remains to be seen if the fluid flow is essential for generating strong enough lateral gradients that lead to buckling. Moreover, further experiments are necessary to explore the range of thicknesses for which buckling occurs and how it depends on motor and passive cross-linker concentrations and the lateral extension of the gel.
With respect to applications, it will be interesting to combine our system with recently introduced methods to spatially and temporally structure the activation of myosin motors 8,9 . In this way, it should be possible to generate active origami, where the system contracts into complex, pre-designed three-dimensional structures. By tuning the thickness of the assembly chambers, one can generate growing and actively contracting sheets that can mimic the growth and folding of epithelia. In particular, we anticipate that our experimental system can be used to study the similarity and differences between systems that fold because of differential growth and those that spontaneously fold due to contraction in the absence of growth.

Methods
Protein purification. G-actin was purified from rabbit skeletal muscle acetone powder by gel filtration (HiPrep TM  skeletal muscle was done according to standard protocols 42 . Myosin II was labeled with Alexa-Fluor 568 at pairs of engineered cysteine residues (see details in ref. 5 ).
Network assembly-first, myosin II aggregates were prepared by bringing the stock motor solution (at 0.5 M KCl) to the final KCl concentration used in the experiment. Actomyosin network formation is initiated by transferring the preformed motor aggregates into the motility medium (see above). Depending on the experiment, we placed 0.3-7 μL of that solution between a glass slide and a glass coverslip (the exact drop volumes are given for each experiment in the corresponding figure captions). We varied the drop volume to generate circular gels of different diameters. The gel initial diameter was determined by the chamber height h and drop volume. To prevent protein adsorption, the glass coverslip and slide were cleaned using piranha solution and coated with an inert polymer (PEGmal M w = 5000 g mol −1 (Nanocs)). The effect of viscosity on contraction dynamics was determined by adding glycerol to the motility medium. Specifically we varied the percentage of glycerol from 0 to 34 to obtain an increase in solvent viscosity from η w to 2.76η w where η w is the viscosity of water at 25°C 43 .
Cell chamber preparation. Chambers were prepared using a glass slide and a glass coverslip or two glass coverslips. The chambers were sealed and delineated by Teflon, Parafilm, or tape spacers located between a glass slide and a coverslip. The choice of a given type of spacer sets the cell height h to a given range. We used various drop volumes to generate actin sheets of different lateral extensions. In our experiments, the chamber height h varied between 70 and 250 μm, and the drop radius ranged between 0.12 and 0.5 cm.
Microscopy techniques. Samples were imaged within 1-2 min after mixing, with an Olympus IX-71 or Leica DMI3000 inverted microscopes. The samples were excited at 488 and 561 nm and the images were recorded simultaneously in two channels using a Dual view Simultaneous Imaging System (Photometrics) with an Andor DV887 EM-CCD camera. Movies overlaying both channels' acquisitions were created using the Metamorph software (Molecular devices). The 3D structure of the sheet at steady state was determined by laser scanning confocal microscopy. Confocal micrographs of the gel at steady state were collected using Leica SP5 laser scanning confocal system on a DMI6000 microscope and analyzed with the LAS X software (Leica Microsystems). To follow the dynamics of contraction in 3D, we used spinning disk confocal microscopy. We used an Olympus IX-81 inverted microscope, supplemented with an Andor XD system with Andor DU-897 camera and acquisition with Andor IQ software. Analysis of the gel contraction and solvent flow dynamics was performed using Matlab (MathWorks) and Metamorph.
Edge velocity of contractile actomyosin sheets. The analysis of the lateral contraction velocity (i.e., edge velocity) was performed by measuring the changes in sheet radius (extracted from changes in sheet area), which is possible for gels with a radius of up to 1.5 mm. The sheet radius was determined by first thresholding the fluorescence images. The area of the super-threshold region was then determined. For each velocity plot we extracted the acceleration (a), the maximal velocity (v max ), the time at which the maximal velocity was reached (t max ), the lateral decay time τ, and final lateral strains ε andε. Let R denote the initial radius of the sheet, r max the radius of the sheet at t = t max , and r end its radius at the end of the contraction process (t = t end ). Then ε ¼ RÀr end R is the lateral strain at the end of the planar contraction phase relative to the initial extension of the sheet at t = 0 and ε ¼ r max Àr end r max is the lateral strain at the end of the planar contraction phase relative to the extension of the sheet at t = t max .
The contraction velocity of the gel edge in the z-direction was determined by measuring the changes in gel thickness e with time. The edge velocity corresponds to changes in the gel half-thickness. The decay time τ z was then evaluated by fitting the data with an exponential function.
Actin and myosin spatial velocity profiles. Particle image velocimetry (PIV) was used to determine the vector velocity fields and velocity profiles of contractile networks. PIV was also used to compare the vector velocity fields of actin (v actin ) and myosin motors (v myosin ). Analysis was performed using Matlab particle image velocimetry (PIV) statistical tool 44 . Velocity vectors calculation was based on a chosen interrogation square window with edge sizes of 40, 32, or 20 pixels with an overlap of 20, 16, or 10 pixels between each window, respectively. Determination of the average particle displacement field was accomplished by computing the spatial auto-correlation of the particle images. The same procedure has been repeated for all images of the contraction movie, providing a description of the network velocity dynamics during contraction.
Actin velocity-myosin velocity correlation. To obtain the local actin velocitymyosin velocity correlation (Supplementary Figure 3), we considered pairs of velocities of actin and myosin in the gel and determined the angle θ between them: . This process was performed for all velocity pairs across the gel surface.
Measuring fluid flow during contraction. The fluid flow was characterized by adding fluorescent beads (200 nm diameter Fluoresbrite YG Microspheres from Polysciences and 2300 nm diameter Nile red beads from Spherotech) to the motility medium. Large beads were used to analyze the fluid flow during the initial and intermediate contraction phases, whereas small beads were used in the final stage of contraction, where the mesh size would have hindered the motion of the large beads. For the experiments, the beads were first incubated with actin monomers and then Bovine Serum Albumin (BSA) to reduce their interaction with the actin network. Particle tracking (Metamorph and Matlab) was used to extract the center-of-mass position of the beads as a function of time.
Local bead velocity-gel velocity correlation. The bead velocity v bead was calculated by taking the bead's center-of-mass position at times t, x t ð Þ; y t ð Þ ð Þ bead and t þ Δt x t þ Δt ð Þ ; y t þ Δt ð Þ ð Þ bead , and dividing by Δt, that is To determine the local gel velocity v gel , we used defects in the gel or points of filament crossings as fiducial markers. The gel velocity was determined by taking the position at times t, x t ð Þ; y t ð Þ ð Þ gel , and t þ Δt x t þ Δt ð Þ ; y t þ Δt ð Þ ð Þ gel , and dividing by Δt. Thus, v gel ¼ Buckling wavelength λ, gel thickness b, chamber height h. Laser scanning confocal microscopy was used to visualize the 3D structure of the gel sheet at steady state and determine the colocalization of actin and myosin motors within the gel. From the micrographs, we extracted parameters value including: the chamber height h, the final sheet thickness b, the buckling wavelength λ at steady state (as a measure of λ c ), and final vertical strain ε z . The thickness b was evaluated from three confocal cross-sections in the xz direction-one passing through the gel center, the others at the two gel boundaries-and three similar cross-sections in the yz direction. For each cross section, we measured the gel thickness at equidistant positions along the gel. Typically, 30-90 values were thus extracted for each gel. The error bars for b indicate standard deviation of these values. The mean final thickness was then used to calculate the final vertical strain ε z . Let h denote the initial thickness of the sheet and b the final thickness of the gel sheet. Then ε z ¼ hÀb h is the vertical strain at the end of the vertical contraction phase relative to the initial thickness of the sheet. The buckling wavelength λ was estimated from the fluorescence intensity along a line tracing the gel periphery in xy confocal cross-sections, where we measured the distance between adjacent intensity maxima and between adjacent intensity minima (see Supplementary Figure 8 for example). The measurements of the buckling wavelength were done on gels with an initial radius of at least 3 mm to have several maxima in steady state; typically, 6-12 maxima and 6-12 minima, depending on the gel size and thickness. The positions of the maxima have been obtained by visual inspection. In addition, we looked at cuts through the gel profile in the xz and the yz directions. We determined the position of the maxima of the profile by visual inspection. Both methods yielded very similar distributions of the buckling wavelength.
Estimating the gel mesh size. The initial mesh size ξ 0 was estimated directly from the actin images at contraction onset. In a first step, we identified the pores in the network by visual inspection. This was possible, because they were delimited by actin bundles, which were clearly detectable in the fluorescence micrographs. Then, we determined the distances between two pairs of opposing bundles of a mesh and obtained the corresponding mesh size from their geometric mean. After contraction was initiated, the mesh size ξ was always determined at the center of the sheet. For each gel, we have measured the size of 55-90 pores.
Actin bundles contour length l cont , end-to-end distance l end-to-end . We considered actin bundles were limited by branching or by crosslinking points (see Supplementary Figure 4). To determine the contour length of a bundle (l cont ), we traced the latter manually and used the 'Region Measurements' function of Metamorph. To determine the end-to-end distance (l end-to-end ), we measured the distance between the two end points of a bundle.
Estimating the values of κ, σ act , and λ c . The values of the gel effective elastic modulus are κ ¼ ηRr max h bξ 2 0 τ (Eq. 3), active stress σ act ¼ κε, and wavelength at buckling onset λ c ¼ 2π ffi ffiε p b (Eq. 4) (see main text). For a single gel, the estimation of κ required the values of b, ξ 0 , R, r max , h, and τ. The value of b used in the expression for κ was the average of the 30-90 values taken at different positions across the gel surface, see above in Buckling wavelength λ, gel thickness b, chamber height h'. Also, we used the average values of ξ 0 for each gel, see above. Finally, the values of R, r max , and τ were evaluated as explained above in 'Edge velocity of contractile actomyosin sheets'. The value of the initial sheet thickness h was extracted from laser confocal micrographs. For a single gel, the value of the active stress σ act was calculated from κ andε, whereε is the lateral strain at the end of the planar contraction phase relative to the extension of the sheet at t ¼ t max (see definition in 'Edge velocity of contractile actomyosin sheets'). In Fig. 7, the values of b, λ, κ,ε, and σ act were averaged over 3-4 gels that were analyzed for each chamber height h; error bars indicate the standard deviations for these 3-4 gels. Note that for the calculation of κ,ε, and σ act , we used small gels (R < 1.5 mm), whereas b and λ were measured on large gels (R > 3 mm). The coefficient linking λ c and b (Fig. 7d) was obtained from λ c ¼ 2π ffi ffiε p b, where we took the value ofε averaged over the gels for all initial heights h (N gels ¼ 13 gels).
Data availability. The data that support the findings of this study are available from the corresponding author upon reasonable request.