Mathieu Leocmach

Mathieu Nespoulous

Sébastien Manneville

Thomas Gibaud
email: thomas.gibaud@ens-lyon.fr

Hierarchical wrinkling in a confined permeable biogel

Confined thin surfaces may wrinkle due the growth of an excess of material. Usually the wavelength is set by elasticity or gravity, but here we explore new selection mechanisms based on hydrodynamics. First, inspired by yoghurt-making processes, we use caseins, a family of milk proteins, as pH-responsive building blocks and the acidulent glucono-δ-lactone, to design a porous biogel film immersed in a confined buoyancy-matched viscous medium. Provided specific boundary conditions, but without any external stimulus, the biogel film spontaneously wrinkles in cascade. Second, using a combination of titration, rheology, light microscopy and confocal microscopy, we demonstrate that, during continuous acidification, the gel first shrinks and then swells, inducing wrinkling. Third, taking into account both Darcy flow through the gel and Poiseuille flow in the surrounding solvent, we develop a model that predicts correctly the wrinkling wavelength. Our results should be universal for acid-induced protein gels as they are based on pH-induced charge stabilisation/destabilisation, and therefore could set a benchmark to gain fundamental insights into wrinkled biological tissues, to texture food or to design surfaces for optical purposes.

INTRODUCTION

There are many ways and reasons for a film to wrinkle. An elastic film buckles due to excess area with respect to its boundaries and wrinkles when buckling is hindered by a substrate. The selection of the preferred wrinkling wavelength is peculiar to the situation and is related to the mismatch of elastic properties between the surface and the substrate [1][2][3][4][5][6]. On the biology side wrinklingcontrolled morphogenesis is ubiquitous. Ageing and the loss of elastic fibres makes our skin wrinkle [6,7]. Difference in growth rates between the gut tube and its dorsal anchoring is responsible for the vilification of guts [8][9][10]. Localised cell death in biofilms focuses mechanical forces and initiates 3D labyrinth pattern [START_REF] Trejo | Proc. Natl. Acad. Sci. U. S. A[END_REF][START_REF] Asally | Proc. Natl. Acad. Sci. U. S. A[END_REF]. On the physics side the last two decades have seen a bloom of methods to obtain well-controlled patterns via linear [4,6,[START_REF] Hu | [END_REF][14][15][16] or nonlinear wrinkling [17][18][19][20]. Such patterns can be triggered by temperature dilation [4], swelling [START_REF] Hu | [END_REF]14] or the removal of pre-strain [6].

Yet, benchmark experiments that explore the possibility of wrinkling in confined porous soft materials immersed in a buoyancy-matched viscous medium are in line, not only from a physics point of view, but also to get fundamental insights into biologically relevant situations. Indeed, stability analysis of a film lying on a thin viscous substrate was performed only recently [21] and remains theoretical because, on a free interface, gravity dominates over viscous substrate. In a biological context, however, thin tissues, e.g. epithelimum or endothelium, are naturally confined and immersed in a nearly buoyant medium, e.g. lymph, blood or mucus, a practical situation where a thin viscous substrate would set the wrinkling wavelength [22]. Moreover, most biological films are porous, and little is known on the interplay between permeability and wrinkling [23,24]. The main goal of the present study is therefore threefold: to create a buoyancy-matched material that wrinkles, to unveil the physical origin of the observed patterns and to identify the selection mechanism of the wavelength.

Here, we introduce model permeable biogels that produce hierarchical wrinkling as a result of their interplay between their gelation dynamics and the confinement conditions. As a biogel, we use a simplified yoghurt obtained by the slow acidification of caseins, a family of milk proteins [START_REF] Roefs | Structure of acid casein gels: a study of gels formed after acidification in the cold[END_REF][START_REF] Bremer | [END_REF][27][28]. While still liquid, the casein dispersion is injected in a home-made optical cell that is then sealed. This cell is designed such that the casein adhesion to both top and bottom walls is turned off and so that the gel is only anchored to the remaining four vertical sides of the cell. Although no external stimulus is exerted on the system, we observe the spontaneous formation of a wrinkling pattern, as shown in Figure 1 and Supplementary Videos 1 and 2. This pattern shows unique properties: a primary pattern appears simultaneously throughout the system with a wavelength λ that is much smaller than the cell width L and does not evolve with time. Furthermore, within this primary pattern, secondary and ternary patterns appear in cascade. Using a combination of titration, rheology, light microscopy and confocal microscopy, we demonstrate that the wrinkles at the millimetre scale result from the spontaneous swelling of the casein network at the micron level upon acidification. Finally we show that the wavelength is selected by a purely viscous mechanism combining porous Darcy flow and viscous Poiseuille flow. We systematically checked that this combined model pins down the dynamical origin of the constraints exerted on the gel and nicely predicts the wrinkling wavelength under a wide range of conditions.

RESULTS AND DISCUSSION

Dynamics of the wrinkling process

To form a biogel we start from an aqueous dispersion of caseins and glucono-δ-lactone (GDL). GDL slowly acidifies the medium and triggers gelation of the caseins [START_REF] Roefs | Structure of acid casein gels: a study of gels formed after acidification in the cold[END_REF][START_REF] Bremer | [END_REF][27][28]. While still liquid, this yoghurt-forming suspension is injected into a home-made optical cell, as sketched in Figure 1a. The cell is sealed and no external stimulus is applied during gelation or pattern formation. Both top and bottom walls are coated with acrylamide brushes, effectively turning off the casein adhesion, so that the gel is only anchored to the remaining four vertical sides of the cell.

The dynamics of the pattern formation are followed by light transmission microscopy and confocal microscopy, see Materials and Methods. As shown in Figure 1b and Supplementary Video 1, we observe the spontaneous formation of a first pattern that appears with a characteristic wavelength λ L. The 2D shape and size of this primary pattern are set as soon as the pattern forms and do not evolve in time. However inside the primary pattern a secondary pattern emerges, inside which a ternary forms. At each step of this Russian dolls-like cascade, the wavelength is divided by two, as shown in Supplementary Figure 1.

3D confocal microscopy, as shown in Figure 1c and Supplementary Video 2, reveals that those patterns actually correspond to wrinkles of the biogel film. Figure 2a shows successive side views of the same confocal measurements as in Figure 1c, see also Supplementary Videos 6 and 7. Initially the dispersion is homogeneous and fills the entire slit. The black margins in Figure 2a correspond to the glass walls. Around 23 min, the gel forms and simultaneously becomes thinner and more concentrated -the fluorescence signal which is proportional to the casein concentration is brighter compared to the initial situation. This effect is the fingerprint of syneresis: as gelation proceeds, solvent is expelled through both the upper and lower sides of the gel film.

Extracting the top and bottom surface of the gel phase allows for quantitative measurements of the relative volume, Figure 2b, and excess surface area of the gel film, Figure 2c, see Supplementary Materials. The evolution of the gel volume is surprisingly non-monotonous. After reaching a minimum around 35% of its initial value at 23 min, the volume increases. Shortly after, at 36 min, the planar film destabilises and starts to wrinkle in the zdirection creating excess area. This confirms that the instability leading to the primary pattern is due to swelling. Around 38 min, the amplitude of the wrinkles reaches the thickness of the slit. As swelling continues, the gel flattens on both top and bottom walls (44 min) and finally buckles back (53 min) to form the secondary pattern (1 h21). When the excess area saturates, the dynamics stop (3 h). We note that at any given point in time the thickness h of the film of gel is spatially homogeneous.

The primary, secondary and tertiary patterns follow the same path: the birth of an instability of wavelength λ of growing amplitude A(t) until saturation to A max , see Figure 2a. Throughout the whole process, λ remains constant. This path is dynamic in essence. In Figure 2d, we thus report the velocity v(t) of the wrinkles as it de-velops in the z-direction, see Supplementary Materials for the definition. We observe that the velocity peaks at different times corresponding to syneresis and to the respective growths of the primary, secondary and tertiary patterns. We associate each growth with a characteristic time τ = A max /v max . For example, in the present case the τ corresponding to the primary pattern is 123 s.

Therefore contrary to observations in, e.g., Ref [29], we do not report any coarsening of the patterns. We believe that the wrinkling pattern is stable against coarsening because of the bending rigidity of the gel [30]. Indeed when the amplitude saturates the gel flattens against the walls which creates sharp edges. To accommodate those sharp edges with finite bending rigidity, the flat part of the gel actually ripples, see Supplementary Figure 2. The ripples created by opposite sides of a same flat patch interfere, creating barriers against coarsening [30].

Based on those observations, we examine possible scenarios upon shrinking and swelling of the gel. As the gel forms it expels the solvent and detaches from top and bottom wall, see Figure 3b andc. For the swelling back we consider three options: (i) In Figure 3d, the swelling back is reversible, only releasing tensile stresses accumulated during shrinking, and the layer remains flat, as expected for chemical gels [14]. Since this case is not observed it means that in this physical gel, syneresis and swelling processes foster irreversible stress dissipation and structural changes. In the following, we will therefore consider that the shrunk state is stress-free. (ii) In Figure 3e, we assume that nothing hinders the buckling of the gel film. This simple buckling situation leads to a primary pattern composed of a single bump of size L. As the bump amplitude grows larger and saturates on the top and bottom wall, buckling of smaller wavelength may appear, as observed by [31]. This is also not what we observe since our primary wavelength λ is much smaller that L. (iii) We are left with the case of Figure 3f where a wrinkling pattern emerges spontaneously with a wavelength λ L selected before saturation on the top and bottom walls. To obtain such an instability, one must balance the bending rigidity of the gel film with a transverse load σ ⊥ that hinders buckling on large wavelengths. The transverse load could be due to gravity [32][33][34][35][36], however we confirmed that gravity plays no significant role by successfully repeating the wrinkling pattern in a cell held vertically, i.e. with the weight acting longitudinally rather than transversely. Furthermore, the transverse load cannot be adhesion [37] since the film is not adhering to the top and bottom walls. The binding of a longitudinal edge to a rigid substrate can induce wrinkles upon compression [16,38]. However in that case, the wavelength gets larger away from the bound edge. Since the wavelength of our pattern is uniform throughout the sample we exclude the influence of the side boundaries. Finally as the four side boundaries remain fixed σ ⊥ cannot come from the resistance to uniaxial stretching [5].

Therefore at this stage, our observations still raise two unresolved issues which we address in the next sections: (i) what triggers the the spontaneous syneresis and swelling back processes? And (ii) what is the nature of σ ⊥ responsible for the selection of the primary wavelength?

Syneresis and swelling back mechanism

The syneresis and swelling back mechanism depends on (i) the pH sensitivity of caseins, (ii) the route to gelation, (iii) the gel properties and (iv) the boundary conditions. Around the isoelectric pH, the gel forms and immediately expels solvent leading to (c). If tensile stresses were not released we would observe reversible swelling back leading to (d) a flat swollen gel layer. If tensile stresses are released, swelling back leads either to buckling (e) or to wrinkling (f) depending on whether the bending is free or hindered by a transverse load.

The stability of casein dispersions strongly depends on the pH. Supplementary Video 3 shows the titration of a casein dispersion at 4%w by 1 molar HCl. At pH=7, the caseins are well dispersed then they flocculate between pH=5.5 to 4 and finally at pH=3 the solution becomes clearer as caseins partially redisperse. The state diagram of caseins as a function of pH is thus related to charge stabilisation. At neutral pH caseins are charged negatively and repel each other. As the pH decreases to the isoelectric point of the caseins (pH≈ 4.6), casein molecules lose their negative net charge and tend to stick to each other [START_REF] Roefs | Structure of acid casein gels: a study of gels formed after acidification in the cold[END_REF][START_REF] Bremer | [END_REF][27][28]. Beyond the isoelectric point caseins acquire a positive net charge and therefore partially re-dissolve. In the inset of Figure 4b, we quantify this effect and measure with UV absorption spectroscopy the amount of free caseins within the 4%w dispersion as a function of pH, see Supplementary Materials. Starting at pH=6.8 all the casein molecules are well dispersed and the free casein fraction x free is 1. Upon decreasing the pH to the isolectric point the free casein fraction drops to almost 0: almost all the casein molecules are aggregated. Further decrease of the pH of the solution leads to a partial redissolution of the casein aggregates. At pH=2 the fraction of free casein is 0.35.

When acidified by a strong acid-like HCl, caseins flocculate and do not form a space-spanning gel. Therefore, we use glucono-δ-lactone (GDL) a molecule that slowly hydrolyses into gluconic acid and thus lowers slowly, continuously and homogeneously the pH of the suspension, Figure 4a. This GDL-controlled acidification enables the formation of a bulk gel phase below pH=5.3 with x free ≈ 0. This process is kinetically driven by diffusion limited cluster aggregation [START_REF] Bremer | [END_REF]. As shown in Figure 4a and b, when adding enough GDL the pH decreases beyond the isoelectric point shifting the balance between attraction and charge repulsion at the casein level, so that caseins partially dissociate from the fractal network. This results in the coexistence between a dilute phase x free ≈ 5 -20% and a dense glassy network, a scenario reminiscent of arrested phase separation [39,40].

We monitor the GDL-induced gelation with a rheometer in cone-plate geometry. As casein strongly adsorbs to surfaces, the strain and stress are well transmitted to the casein dispersion and we can measure the elastic and viscous moduli of the casein dispersion as a function of time, Figure 4c. The mechanical properties of the gel are closely related to the fraction of free casein. First, the dispersion solidifies then reaches its maximum elasticity (≈ 500 Pa) at the isoelectric point and finally the elastic modulus decreases down to ≈ 100 Pa. This finite value indicates that the sample is still a gel.

To gain further insight into the mechanical properties of the gel we turn to microstructural observations under confocal microscopy. As shown in Supplementary Videos 4 and 5 and in Supplementary Figure 4 the gel consists in a space-spanning fractal network with a cutoff length scale ξ that corresponds to the maximum pore size [41]. Permeability α is of the order of 67 000 nm 2 , see Supplementary Figure 5. The fractal model fits the structure at all times but its parameters vary non monotonously during GDL acidification. As shown in Figure 4d and Supplementary Figure 4, when reaching the isolectric point ξ increases up to a maximum length of 7 µm due to the solvent expulsion out of the casein network [27] like in viscoelastic phase separation [42]. This corresponds to a tighter packing of the caseins within the network strands. Such stiffening of the strands explains the observed increase of elasticity. Beyond the isoelectric point caseins partially regain solubility, the packing within the strands loosens which results in weaker, thicker strands and smaller pores as ξ decreases down to 4 µm, hence the long-term decrease of the elastic modulus.

Boundary conditions is the final parameter that is essential to form wrinkled patterns. In order to turn off the adsorption of the caseins on a surface, we coat the surface with polymer brushes, see Materials and Methods, effectively creating non adhesive boundary conditions for the gel. When the casein adhesion to the upper boundary is turned off, the gel detaches from the top wall. However the gel does not wrinkle because it remains attached to the bottom wall. In this situation we can monitor ξ, Fig- ure 4d. ξ remains mostly constant and variations manifests themselves at the macroscopic level, on the volume occupied by the gel.

Therefore, the syneresis and swelling back mechanism is related to charge stabilisation of the casein proteins as function of the pH. In particular, beyond the isoelectric point, the casein are charged and repel each other, leading to a swelling at every scale: proteins are further from each other, the strands become larger and longer and the gel sheet swells. Such a mechanism is only effective to produce wrinkles when adhesion of the caseins to both top and bottom walls is turned off. Consistently with our explanation, the amount of GDL, and thus the final pH reached by the suspension, controls the amount of swelling back. Thus at low GDL concentrations we only observe the primary pattern, at high GDL concentrations we observe the formation of higher order patterns.

Wavelength selection mechanism

The determination of the characteristic length λ of a wrinkling pattern is an old problem. Large deformations of thin sheets are governed by Föppl-von Kármán equations [START_REF] Landau | Course of theoretical physics: theory of elasticity[END_REF], which are essentially impossible to solve in analytical form. However using scaling and asymptotic arguments one can obtain the dominant wavelength of the wrinkling pattern up to dimensionless prefactors and isolate the physically meaningful ingredients [5,[START_REF] Davidovitch | [END_REF]. An elastic film buckles due to excess area with respect to its boundaries. Wrinkling is a buckling hindered by a substrate: the mismatch of elastic properties between the film and the substrate selects the preferred wrinkling wavelength [1,2]. λ is then set by the competition between the flexural modulus of the film B = E/(1 -ν 2 ) (E Young modulus, ν Poisson ratio) and a transverse load σ ⊥ . For small deflections of amplitude A, λ scales as [33,34]:

λ ∼ Ah 3 B σ ⊥ 1/4 . (1)
This framework has been successfully applied to situations where σ ⊥ comes from the resistance of an elastic substrate [1,2]. Dealing only with (semi)infinite substrates, Biot generalised the elasticity of both the film and the substrate to viscoelasticity [3] or poroelasticity [45], leading to time-dependent kinetic wavelength selection. Furthermore, it was shown that equation (1) remains valid even when σ ⊥ comes from the resistance to uniaxial stretching [5], from gravity [32][33][34][35][36], from capillary forces [46], from boundaries [16,38] or from adhesion [37]. However, we have excluded these various 5. Sketch of the two limit scenarios for the wavelength selection. In both cases the gel film of thickness h is destabilised by the excess area which gives rise to wrinkles of amplitude A(t) and wavelength λ. A single wavelength is represented and the longitudinal dimension L of the system is much larger than λ. The value of λ is set by the interplay between the bending rigidity B and a transverse load σ ⊥ due to a pressure gradient p2 -p1 in the solvent of viscosity η. (a) Darcy scenario: the gel initially sits without sticking to the bottom wall. σ ⊥ is due to the flow of the solvent through the porous gel of permeability α to fill the growing blister. (b) Poiseuille scenario: the gel film of negligible porosity initially lies in the middle of the cell, separated by a distance H from each wall. σ ⊥ is due to the lubrication flow in the top and bottom solvent layers.

origins for σ ⊥ in the discussion around Figure 3. The only possible explanation left for σ ⊥ involve hydrodynamic stresses due to viscous flow of the solvent, either through the gel porous matrix itself or in the surrounding solvent layers. Accordingly we consider the two limit scenarios sketched in Figure 5 that lead to wrinkling patterns. In a first limit, the film of gel sits without sticking to the bottom wall of the cell, i.e. H 1 → 0, see Figure 5a and Supplementary Video 6. Creating a two dimensional blister requires flowing solvent of viscosity η through the gel of permeability α. Darcy law [START_REF] Darcy | Les fontaines publiques de la ville de Dijon[END_REF]

λ D ∼ h 1/2 α 1/4 Υ 1/4 , with Υ = Bτ η . (2)
Υ is a dimensionless viscoelastic factor quantifying the relative stiffness of the film and the solvent at a characteristic time τ . The opposite limiting case is when an impermeable (α → 0) gel film lies in the middle of the cell (H ≡ H 1 = H 2 , see Figure 5b and Supplementary Video 7).

Here destabilisation over a wavelength λ creates a lubrication (H λ) flow in the viscous layers. By symmetry the transverse load across the gel σ ⊥ = p 2 -p 1 is also the pressure difference over the wavelength. Using a Poiseuille profile for the flow, dA dt ∼ σ ⊥ H 3

ηλ 2 [START_REF] Poiseuille | [END_REF] and injecting this expression in equation ( 1) yields the wavelength of the Poiseuille mode as

λ P ∼ (hH) 1/2 Υ 1/6
(3)

where we recognise the scaling for the wrinkles of an elastic film on a thin (H λ) elastic substrate [5] where Υ would be the ratio between the flexural modulus of the film and that of the substrate. Indeed within the lubrication approximation or at low Reynolds numbers one can consider a viscous film as elastic with effective Young modulus 3η/τ [3,49].

A more complete derivation of the models can be found in the Supplementary Materials and yields the following prefactors:

λ D = 2πh 1/2 α 1/4 Υ 12 1/4 (4)
λ P = π(hH) 1/2 2 9 Υ 1/6 ( 5
)
where H is redefined by 1/H 3 = 1/H 3 1 + 1/H 3 2 . The crossover between the two models where λ D = λ P is found for

H = H * ≡ 2 2/3 3 1/6 α 1/2 Υ 1/6 .
The characteristic time τ is the hallmark of the kinetic nature of the wavelength selection mechanisms and is present in the expressions of both λ D and λ P .

To sum up, we expect two regimes of viscositydominated wrinkling: for H H * , i.e. slow destabilisation of a highly permeable film close to a wall, the main dissipation mechanism being the flow through the porous film; whereas for H H * , i.e. fast wrinkling of a poorly permeable film far from any wall, the flow in the viscous layers dominates and sets the wavelength.

Experimental tests

For a typical gel film, the buoyancy after syneresis is estimated from the initial sodium caseinate mass concentration C cas as ∆ρ = C cas e/h = 130 kg m -3 leading to a gravity stress σ G ⊥ ≈ 0.04 Pa. This is much lower than typical Darcy or Poiseuille stresses σ D ⊥ ≈ σ P ⊥ ≈ 0.15 Pa which confirms that gravity can be neglected as already inferred from test experiments where the cell is held vertically rather than horizontally.

Darcy and Poiseuille models yield wavelengths that are of the same order in most of the cases under study. Thus only systematic measurements of the experimental λ together with estimates for all the parameters involved in either model shall allow one to discriminate between them, see Supplementary Materials about such measurements. As summarised in Supplementary Table I we systematically vary (i) the thickness of the cell, impacting h and the maximal amplitude, (ii) the gel composition and thus both its stiffness and its permeability, (iii) the viscosity of the solvent by adding glycerol, which also impacts the gel stiffness, permeability and τ . Finally, the only parameter that we have no control upon is the initial altitude H of the gel before wrinkling that varies randomly.

For each experiment we measure independently λ and all the parameters that go into the expressions of λ D , Equation (4) and λ P , Equation (5). λ is measured by optical microscopy; e, h, H 1 , H 2 , τ are measured by confocal microscopy; α is measured via capillary experiment and B by rheology, see Materials and Methods and Supplementary Materials.

The beauty of this approach is that it is not only valid to model the primary wavelength, but that it describes also the growth of the secondary pattern. Although the wavelength λ 2 of the secondary pattern is fully determined by the primary pattern λ 2 = λ 1 /2, Figure 1 and Supplementary Figure 1, the secondary pattern has its own dynamics with its own characteristic time τ 2 = τ 1 . Since the saturation of the primary pattern flattens the gel on the top and bottom walls, the secondary blisters always appear in a Darcy situation (H H * ). Therefore, we are able to obtain more data at small H/H * ratios by considering not only the primary patterns but also the secondary ones.

Despite changing the gel properties and the cell geometry we always observe wrinkling proving the robustness of this approach to form wrinkling patterns. We are in the position to vary H/H from 0.1 to 10 so that we can test the crossover between Darcy and Poiseuille mechanism. In Figure 6a and 6b we have plotted the measured wavelength as a function of the model prediction λ D and λ P respectively. If the predictions are correct we expect λ exp = λ D when H < H * and λ exp = λ P when H > H * . Indeed, both model seems to apply in their domain of validity. The failure of the Poiseuille model is only visible at small wavelengths where an unexplained offset (dashed line in Figure 6b) is necessary to fit the data correctly. To go further, we derive a complete model summing up the influence of both Poiseuille and Darcy dissipation, see Supplementary Materials. The results of this mixed model are displayed in Figure 6c and agree with λ exp over the whole H/H * range. Despite a very good agreement of the scaling, the exact prefactors of all models disagree by a factor 0.63-0.69. We suggest that this discrepancy comes from the fact that our model is 2D, i.e. uniaxial compression and deflection, whereas our system is 3D, i.e. biaxial compression and deflection. Since our system is isotropic in the plane, this does not change the scaling but might account for a dimensional factor.

CONCLUSION

To conclude we have taken advantage of the properties of the casein protein to engineer and study the wrinkling of a porous confined biogel. In particular, we have fine-tuned the interplay between the attraction and the repulsion at the protein level to assemble dense glassy networks that shrink and swell as a result of continuous acidification. Together with carefully controlled adhesion of the gels to the boundaries, this route to form permeable biogels is the key to produce an original cascade of wrinkling patterns. We developed a model that captures the dynamical origin of the constraints exerted on the gel and predicts the wrinkling wavelength: Poiseuille flows of the solvent above and below the gel when the film is far from any wall, Darcy flow of the solvent through the gel when the film is close to a wall. Such a model experiment, because it relies on pH-induced charge stabilisation/destabilisation, should be generalisable to different kinds of proteins. As such it could set a benchmark to explore potential applications in micro-fingerprinting [50] or in soft optical devices such as diffraction gratings or Fresnel lenses [16], food texturing [51,52] and the possibility of wrinkling in confined porous soft materials immersed in a buoyancy-matched viscous medium such as biological tissues [22][23][24].

MATERIALS AND METHODS

Sodium caseinate (Firmenich) is dissolved in deionised water by stirring for 20 min at room temperature. The mixing of the caseine dispersion with glucono-δ-lactone (GDL) in powder (Firmenich) sets the beginning of the experiments. We typically have about 10 min to pipet the dispersion into the slit before gelation starts. Sodium caseinate (Firmenich) is labelled with Dylight 550 NHS ESTER (Thermo Scientific). Excess dye is removed by centrifugation. For confocal microscopy measurements, we add 10% of fluorescently labelled caseins to the casein/GDL dispersion while it is fluid.

Home made optical cells are composed of a slide (RS) and a cover slip (Menzel-Gläser) spaced by paraffin film (Parafilm). The optical cell is briefly heated on a hot plate so that the paraffin welds to the slide and the cover slip. Gently applying a pressure on the cover slip while being heated allows us to tune the thickness of the cell from 50 to 150 µm. The slit is filled with the casein dispersion, immediately sealed using ultraviolet-cured glue (Norland Optical) and placed under the microscope before gelation can take place. To prevent adsorption of caseins to surfaces, microscope slides and cover slips used to make microscopy cells are cleaned and subsequently coated with a silane agent (3-(trimethoxysilyl) propyl methacrylate, Sigma Aldrich) from which a polyacrylamide brush is polymerised from methacrylate groups (Euromedex).

3D data are collected on a Zeiss LSM510 confocal microscope using 532 nm laser excitation. We use either a 10x (air) lens to measure the geometry of wrinkling, see Supplementary Material, or a 63x (oil) lens to observe the microstructure of the gel. Larger scope pictures, obtained either by stitching fluorescent microscopy images (Nikon Eclipse Ti) or by transmitted/reflected light macroscope (Nikon SMZ745T/Leica DMS1000) are used to measure the wavelength of the primary pattern, Supplementary Figures 1 and3. The secondary pattern wavelength is taken as half of the primary wavelength, according to Supplementary Figure 1. Data shown in Figure 4a and c are obtained with a SevenCompact pH-meter (Mettler Toledo) and an MCR 301 rheometer (Anton Paar). For every gel composition we monitor the gelation by small amplitude oscillations (strain 0.1%, frequency 1 Hz) and estimate the flexural modulus B from the maximum value of the shear modulus G as B = 2G /(1 -ν) taking ν = 0.3 for the Poisson ratio, a value typical of a spongy network [53]. Permeability is measured using a protocol found in [41], see Supplementary Method and Supplementary Figure 5. The amount of free casein protein present in the gel during the GDL acidification process is measured by titration and UV-spectrometry (Ocean Optics) according to the protocol found in [START_REF] Roefs | Structure of acid casein gels: a study of gels formed after acidification in the cold[END_REF] and detailed in Supplementary Method.

SUPPLEMENTARY MATERIALS

Wavelength measurements

To measure the wrinkle wavelength λ, we use large field pictures of the final state. These picture are obtained either by stitching fluorescent microscopy images obtained at 10x magnification, Supplementary Figure 3 except c ande, or by transmitted or reflected light macroscope at magnifications ranging from 1x to 6x, Supplementary Figure 3c ande. The contrast in such pictures is a measure of the magnitude of the altitude gradient of the gel film. White lines in fluorescence and reflection (respectively black lines in transmission) denote a sharp transition between two extreme altitudes (on top wall or on bottom wall). It is thus possible to reconstruct a binary picture of the pattern, inset of Supplementary Figure 1. Furthermore since the shape and size of a given generation of pattern do not change after formation, it is possible to revert the final binary picture to the configuration of the primary pattern.

Figure 1 shows Fourier spectra of binarised altitudes at different times. When only the first generation pattern is present a peak clearly indicates the dominant mode. When the secondary pattern has appeared, the dominant wavelength is halved and the peak much weaker. This allows us to use the secondary pattern is our analysis by setting their wavelength to half of the primary wavelength. Further generations do not seem to appear simultaneously and the corresponding peaks are weak and broad.

Confocal images processing

At 10x magnification and depending on cell thickness, we acquire stacks of 20 to 50 images spaced by 6.69 µm (Nyquist sampling). Each confocal slice is (256 px) 2 with pixels of 4.97 µm. At such low resolution details of the gel microstructure are undistinguishable as ξ is smaller than the pixel size for all compositions and the gel appears as a continuous medium.

We first de-noise each plane by a short-range Gaussian blur (σ = 2 px λ). Then, for each (x, y), we focus on the z-dependent intensity profile I(z). We define the position of the top surface of the gel z top as the point where I(z) first crosses half of its maximum value over the given vertical line. Since this criterion is local, we do not need to correct for background intensity variation.

In this way we obtain z top (x, y) with subpixel resolution (±0.1 px = 0.67 µm), and z bottom (x, y) in a symmetrical way.

From these data, we obtain for each stack the volume and thickness of the gel phase, area and peak-to-peak amplitude of the mid surface. Velocity is obtained by time differentiation. Once the instability has saturated on the top and bottom walls of the cell, peaks in the histograms of z top and z bottom yield the measure of the positions of both walls, thus of the cell thickness e and retrospectively of H 1 and H 2 .

Characteristic time measurements

In most of our samples it is possible to measure directly the peak-to-peak amplitude A(t) from our confocal measurements. However when the wavelength is too large compared to the field of view, this approach breaks down since we cannot observe both the highest and lowest points of the surface. Even when valid, such a local measurement is rather noisy and we prefer to focus on the measure of the vertical velocity v of the gel that can be integrated over the whole field of view:

v(t) = dz dt dxdy dxdy , (6)
where z(x, y, t) is the altitude of the mid surface of the gel.

Fraction of free casein as a function of pH

The fraction of free casein proteins of a 100 ml, 4%w sodium caseinate dispersion is obtained either by titration by 1 molar HCl, see Supplementary Video 3, or during the action of 4%w GDL. In both cases, 0.4 ml of the 100 ml dispersion is sampled regularly as a function of pH, and centrifuged (centrifuge 5424, Eppendorf) in Eppendorf tubes at 14500 rpm for 3 min. The supernatant is diluted and analysed with a UV absorption spectrometer (USB4000, Ocean Optics). According to [START_REF] Roefs | Structure of acid casein gels: a study of gels formed after acidification in the cold[END_REF], we deduce the fraction of free casein proteins from the adsorption spectrum as:

x free = A(pH) A(pH = 7) , (7)
with

A(pH) = A 280 -1.706A 320 , (8)
where A 280 is the absorption peak at 280 nm due to aromatic amino acids in the casein molecules plus a turbidity background whereas the absorption peak at 320 nm is only due to turbidity. Above the isoelectric point, the homogeneous acidification by GDL continuously creates very small protein aggregates that are difficult to centrifuge. Such suspensions are too turbid and too rapidly evolving to reliably measure x free .

Microstructure measurements

To follow the evolution of the microstructure of the gel, we acquire slices of (2048 px) 2 with pixels of 98 nm, Supplementary Figure 4a and Supplementary Videos 4 and 5. After Hann windowing we perform Fourier transform and radially average the resulting spectrum to obtain I(q), see Supplementary Figure 4b that is fitted by the following fractal form [41]:

I(q) = 2χΓ(d) 1 + ((d + 1)ξq) 2 d/2 (9)
where the susceptibility χ is related to the fluorescence contrast between the network and the solvent, the cutoff length ξ corresponds to the largest pore size and d is the fractal dimension. Γ is the Gamma function. Equation ( 9) is a good approximation to the 2D Fourier transform of a fractal pair correlation function g(r) ∼ 1 + e -r/ξ r d-2 . Figure 4d of the main text and Supplementary Fig- ure 4c-e show the evolution of the fitting parameters in two different sets of boundary conditions: (i) complete adhesion of the gel on all walls, preventing syneresis and swelling; (ii) complete adhesion except on the top wall, allowing syneresis and swelling. The latter situation is a proxy to the situation where we observe wrinkling, since it is necessary to keep the gel into the thin focal plane of the 63x objective lens of the confocal microscope. To do so, we coat only the microscope slide with acrylamide brush. The cover slip is cleaned but not coated, thus keeping strong adhesion on the side closer to the objective lens. Neither buckling nor wrinkling is detected in such a configuration.

In situation (i), the gel keeps its topology throughout the gelation and over-acidification. Left of Supplementary Figure 4a shows that the arms of the gel first shrink laterally while keeping their length, and then swell again. Syneresis and swelling are thus local processes. By contrast in situation (ii), the gel as a whole compresses while forming. Since our focal plane is at a fixed depth, we cannot follow the topology of the same subset of the network, however the pore size is clearly smaller than in situation (ii). The isoelectric gel, at minimum gel volume, is probably isotropic or only marginally anisotropic in z. This explains that subsequent swelling is not occurring only along z but also in plane, causing buckling.

ξ is also systematically measured in the same way at the end of each wrinkling experiment and reported in Supplementary Table I.

Permeability measurements

To measure the permeability of our gels, we adapt a method from [41]. First, the gel of interest is formed at the lower extremity of a thin glass tube. The tube is held vertical in a water saturated atmosphere to prevent evaporation until the pH has reached the isoelectric point (≈ 1 h for our most common composition). Then, we dip the tube in a bath of acetic acid/acetate buffer at pH=4. The buffer has the same glycerol content as the gel solvent, so that the viscosity is constant throughout the experiment. Alongside the first tube, see Supplementary Figure 5 a, we dip a second identical tube with no gel in it. A webcam (Logitech Webcam Pro 9000) captures the rise of the liquid in each tube. Since the liquid height in tube 2 reaches its final value quickly, we measure directly z ∞ -z(t).

Darcy law through the gel of permeability α and height

H reads dz dt = α η ∆P H (10)
where ∆P = ρg(z ∞ -z(t)) is the hydrostatic pressure drop through the gel. We use tabulated values of waterglycerol mixtures for the density ρ and viscosity η of the buffer. We can thus write

d(z -z ∞ ) dt = - z -z ∞ τ , with τ = ηH αρg . (11)
By fitting z(t) -z ∞ to an exponential, see Supplementary Figure 5 b-c, we obtain τ and thus α. Note that even with H of the order of a few millimetres high glycerol contents yield τ ≈ 100 h due to both high viscosity and low permeability. Obtaining z ∞ without tube 2 would be time prohibitive.

Mixed Darcy-Poiseuille model

General framework for wrinkling

The mechanical equilibrium of a plate of thickness h, Young modulus E and Poisson ratio ν submitted to a compression load σ along x and a transverse load σ ⊥ along z writes

Eh 3 12(1 -ν 2 ) ∂ 4 w ∂x 4 + σ h ∂ 2 w ∂x 2 = σ ⊥ , (12)
where w(x, t) is the deflection of the plate along z [3]. Since we deal with slow rates of deformation, inertia is neglected.

We decompose the deflection in sinusoidal modes w(x, t) = A sin(qx) of wavevector q. Equation ( 12) then reads

1 12 Bh 3 q 4 -σ hq 2 w = σ ⊥ , (13)
where B = E/(1 -ν 2 ) is the modulus of the plate. Following [3], we consider the case where the compression load σ is increased quasistatically. The observed mode is the one that arises at the lowest possible σ , i.e. the mode minimising σ .

Elastic film in an infinite elastic medium

In this framework, we recall here the case of an elastic film lying over an infinite elastic medium [1][2][3]5]. The deformation of the substrate of modulus B s gives the transverse load as σ ⊥ = -qB s w, and thus the compression load writes Minimising σ as a function of q yields hq = (6B s /B) 1/3 [3,5] so that the dominant wavelength is

σ = B h 2 12 q 2 + B s hq . (14)
λ = 2πh B 6B s 1/3 . (15)

Elastic film on a viscous layer

Let us now turn to the case of an elastic, impermeable film on a single viscous layer. As in [21] we neglect gravity, although it would be difficult to do so in practice at a free interface. Assuming the layer thickness H to be much smaller than λ d , we can use the lubrication approximation. Further considering only small deflections, we may neglect the displacement of the plate along x. The flux along x is thus of the Poiseuille form

Q = - H 3 12η ∂p ∂x . (16)
where p(x) is the pressure in the layer of viscosity η.

Mass conservation (on the lower gray area of Supplementary Figure 6) writes

∂Q ∂x + ∂w ∂t = 0. ( 17
)
Assuming an exponentially growing amplitude A(t) = A 0 exp(t/τ ) and keeping only linear order terms, one gets

12η τ w H 3 q 2 = -σ ⊥ . (18)
Combining Equations ( 13) and ( 18), the compression load writes

σ = 1 12 Bh 2 q 2 + 12 η τ 1 H 3 hq 4 (19)
and is minimised at q P such that

q 6 P = 2 × 12 2 1 h 3 H 3 Υ , (20)
with Υ = Bτ /η the viscoelastic factor. The dominant wavelength for the Poiseuille model is thus

λ P = π √ hH 2 9 Υ 1/6 . (21)
We note that the mode minimising σ is also maximising the growth rate

1 τ = H 3 q 2 12 2 η 12σ -Bh 3 q 4 . (22)
When the elastic film is sandwiched between two viscous layers of respective thickness H 1 and H 2 , the linearity of equation [START_REF] Asally | Proc. Natl. Acad. Sci. U. S. A[END_REF] allows us to use an effective substrate of thickness

1 H 3 = 1 H 3 1 + 1 H 3 2 . ( 23
)
This corresponds to the experimental case where the gel film is far from both of the cell walls. Since we have an almost buoyant gel and no more free surface, gravity can indeed be neglected.

Porous elastic film with no possible Poiseuille flow

Another limiting case is that of a porous elastic film immersed in a viscous medium where longitudinal resistance to flow is infinite. It is the case in our experiments when the film is at contact with one of the walls. The only way to deform the film is to flow some liquid through the film of permeability α. The transverse load derives from the Darcy law:

σ ⊥ = p 1 -p 2 = hη α ∂w ∂t . (24)
Combining equations ( 13) and ( 24), the compression load writes

σ = 1 12 Bh 2 q 2 + η τ 1 αq 2 (25)
and is minimised at q D such that

q 4 D = 12 h 2 αΥ . (26)
The dominant wavelength for Darcy model is thus

λ D = 2πh 1/2 α 1/4 Υ 12 1/4 . (27)
Here also the dominant mode maximises the growth rate.

Porous elastic film between two viscous layers

Finally the most realistic scenario mixes both Poiseuille and Darcy mechanisms, i.e. considers a porous elastic film between two viscous layers of arbitrary thickness. In this case, see Supplementary Figure 6, the mass conservation Q 1 + Q 2 = 0 over the whole height of the cell yields

H 3 1 p 1 + H 3 2 p 2 = 0. (28)
The mass conservation in the lower viscous layer is the same as equation ( 17) with an added leak v due to the porosity

∂Q 1 ∂x + ∂w ∂t + v = 0. ( 29
)
v can be expressed by the Darcy law

v = α η p 1 -p 2 h . (30)
Combining equations ( 28), ( 29) and ( 30) we obtain the transverse load

σ ⊥ = p 1 -p 2 = - η τ w H 3 q 2 12 + α h , (31)
Note that we recover equation ( 18) when α → 0 and equation ( 24) when H → 0. We can recast this intuition in terms of the dimensionless number H/H * where H * is obtained by equating q P and q D : H * = 2 2/3 3 1/6 α 1/2 Υ 1/6 . (

Analytic minimisation of σ knowing σ ⊥ is possible but cumbersome. Instead we rewrite equation [START_REF] Hu | [END_REF] in terms of previously calculated λ P and λ D :

12σ q 2 D h 2 B = q q D 2 + 1 2 q 2 D q 6 P q 4 + q q D 2 (33) = λ D λ 2 + 1 2 λ 6 P λ 2 D λ 4 + λ D λ 2 (34)
and we minimise this expression numerically.

Samples characteristics

To test the theoretical predictions we have fully characterised 9 samples. Their characteristics are gathered in Supplementary Table I

FIG. 1 .

 1 FIG. 1. Dynamics of pattern formation in a confined film of casein gel (caseinate 4%w, GDL 4%w in water). (a) Sketch of the cell where the adhesion on both top and bottom walls is turned off. The cell is sealed and the gel is only anchored to the four sides. Typical dimensions are L > 1 cm, e ≈ 100 µm. (b) Light transmission microscopy. Successive generations of patterns are highlighted in colour to stress the absence of coarsening after formation. The successive wavelengths are λ = 1.5 mm (yellow), λ = 0.75 mm (orange) and λ = 0.32 mm (brown). (c) 3D reconstruction from fluorescent confocal microscopy which highlights that patterns observed in (b) correspond to wrinkles. The contrast in (b) is not due to thickness inhomogeneities but to altitude gradients as indicated by the yellow vertical lines in (a). Scale bars are 1 mm.

FIG. 2 .

 2 FIG. 2. 3D analysis of the wrinkling process. (a) Confocal (x, z) cuts showing syneresis, swelling, wrinkling and cascade buckling. The scale bar is 100 µm (real size ratio). (b-d) Confocal microscopy measurement of the evolution of the volume of the gel phase relative to cell volume, of the excess area and of the velocity along the z-direction. Crosses correspond to times in (a).

FIG. 3 .

 3 FIG. 3. The wrinkling experiment. (a) Initial configuration where the sealed cell contains a homogeneous protein solution. (b)Around the isoelectric pH, the gel forms and immediately expels solvent leading to (c). If tensile stresses were not released we would observe reversible swelling back leading to (d) a flat swollen gel layer. If tensile stresses are released, swelling back leads either to buckling (e) or to wrinkling (f) depending on whether the bending is free or hindered by a transverse load.

FIG. 4 .

 4 FIG. 4. Acid-induced protein gel properties behave non monotonically with pH. (a) pH decrease over time in a 4%w sodium caseinate solution acidified by 4%w GDL in water. The horizontal line indicates isoelectric pH of caseins. (b) Corresponding evolution of the fraction of free caseins. Symbols are direct measures upon addition of GDL. The continuous line is deduced from x free at the same pH obtained by addition of HCl 1 molar and waiting for equilibration (inset). (c) Evolution of the elastic modulus G measured in a rheometer with full adhesion to the cone-plane geometry. (d) Evolution of pore size measured by confocal microscopy in a slit geometry with either adhesion to all walls (black line) or no adhesion to top wall, i.e. allowed syneresis and swelling (orange dashed line).

  FIG. 5. Sketch of the two limit scenarios for the wavelength selection. In both cases the gel film of thickness h is destabilised by the excess area which gives rise to wrinkles of amplitude A(t) and wavelength λ. A single wavelength is represented and the longitudinal dimension L of the system is much larger than λ. The value of λ is set by the interplay between the bending rigidity B and a transverse load σ ⊥ due to a pressure gradient p2 -p1 in the solvent of viscosity η. (a) Darcy scenario: the gel initially sits without sticking to the bottom wall. σ ⊥ is due to the flow of the solvent through the porous gel of permeability α to fill the growing blister. (b) Poiseuille scenario: the gel film of negligible porosity initially lies in the middle of the cell, separated by a distance H from each wall. σ ⊥ is due to the lubrication flow in the top and bottom solvent layers.

  relates the pressure gradient, (p 2 -p 1 )/h, to the volume flux per unit area v = α η p2-p1 h = dA dt . Identifying σ ⊥ = p 2 -p 1 and injecting this expression in equation (1) yields the wavelength of the Darcy mode as

FIG. 6 .

 6 FIG. 6. Comparing model predictions λD, λP and λD+P with measured wavelengths λexp. Dots come from primary pattern, from secondary blisters. Lines are the best linear fits through the origin taking into account only the points that should be (a) in Darcy mode H < H * , (b) in Poiseuille mode H > H * (c) all points. Prefactors are 0.63, 0.69 and 0.67 respectively. The dashed line in (b) is the best affine fit (λexp = 0.52λP + 0.33 mm) to all data points.

- 1 )FIG. 1 .

 11 FIG. 1. Fourier spectra of binary altitudes of Main Figure 1a corresponding to primary, secondary and tertiary patterns. Inset: reconstructed binary altitudes at time 2 h. The scale bar is 1 mm.

FIG. 2 .

 2 FIG. 2. Evidence for the ripples generated by the junction to flat patches, see [30]. (a) Position of the bottom surface of the gel obtained by confocal images. (b) one dimensional cut of the region of interest highlighted in (a) using an orange box. Only the first ripple is visible. (c) Schematic view of the interference between ripples.

FIG. 3 .

 3 FIG. 3. Patterns corresponding to the samples of Supplementary Table I. (a-d) increasing cell thickness. (e-f) larger solid content. (g-i) increasing glycerol content. All pictures are stitching of fluorescent microscopy images except (c,e) which are details of reflected light macroscope images. The scale is common to all panels (scale bar 1 mm). Arrows show the measured primary wavelength.

FIG. 4 .

 4 FIG. 4. Microstructure evolution of a sodium caseinate solution acidified with 4%w GDL in water. (a) Details of confocal images in fully adhesive conditions (top row) and with no adhesion on top (bottom row). The scale bar is 10 µm. Arrows indicates the pore size. (b) Corresponding Fourier spectra when no adhesion on top, same times as (a). Lines are fits to equation (9). The inset shows a schematic association between real space structures and wave numbers. (c-e) Evolution of cut-off length ξ, susceptibility χ and fractal dimension d in full adhesive (black line) and no adhesion on top (orange dashed line) situations. Symbols indicate the times of the images in (a).

FIG. 5 .FIG. 6 .

 56 FIG. 5. measurements. (a) Schematic representation of the experimental set up. (b-c) Evolution of the height of the interface in tube 1 relative to the final height in tube 2. The black line is the best exponential fit Ae -t/τ . (b) Gel is 4%w casein, 4%w GDL in water, H = 2.3 mm and τ = 57 min. (c) Same as (b) for a 50%w glycerol-water mixture, H = 4 mm and τ = 100 h.

  . Transmitted light microscopy images of the pattern formation in a 4%w casein and 4%w GDL aqueous dispersion. Top and bottom slides of the optical cell are coated with acrylamide brushes. The total duration of the Video is 2 h 30. The scale bar is 1 mm. Video 2: 3D reconstruction of the pattern formation using confocal microscopy in a 4%w casein and 4%w GDL aqueous dispersion, Supplementary Table Ic. 10% of the proteins are fluorescently labelled. Top and bottom slides of the optical cell are coated with acrylamide brushes. The total duration of the Video is 2 h 30. The scale bar is 1 mm. Video 3: Titration of a dispersion with 4%w casein by 1 molar HCl. The total duration of the Video is 14 min. The Becher diameter is 55 mm. Video 4: Confocal microscopy (x, y) cut of the evolution of the local structure of a 4%w casein and 4%w GDL aqueous dispersion when sticky boundary are on. 10% of the proteins are fluorescently labelled. The total duration of the Video is 132 min. The scale bar is 10 µm. Video 5: Confocal microscopy (x, y) cut of the evolution of the local structure of a 4%w casein and .9 38.3 6.3 29.3 0.11 1.62 0.54 0.89 1.20 0.88 1.22 40.2 1.7 32.9 0.09 0.70 0.09 0.45 1.00 0.39 1..5 65.1 7.9 14.4 0.02 16.30 1.52 1.23 1.50 1.85 1.96 TABLE I. Characteristics of the samples used for Figure6. Lines where preparation and properties are left blank correspond to the average of the secondary blisters of the previous line. 'cas', 'GDL' and 'gly' indicate the weight fraction of sodium caseinate, GDL and glycerol respectively. ξ is the pore size at the end of each experiment. h is the thickness of the gel film just before buckling.4%w GDL aqueous dispersion when the top wall is coated with acrylamide brushes and the bottom wall remains sticky. 10% of the proteins are fluorescently labelled. The total duration of the Video is 147 min. The scale bar is 10 µm.Video 6: Confocal microscopy (x, z) cut of the evolution of the gel film deflection of a 4%w casein and 4%w GDL aqueous dispersion when both top and bottom walls are coated with acrylamide brushes, Supplementary TableIa. 10% of the proteins are fluorescently labelled. The scale bar is 100 µm (with a pixel size ratio of 1.3). Here the gel initially lies close to the bottom wall, H < H * . Darcy mechanism dominates the formation of the primary pattern. The total duration of the video is 145 min. Video 7: Confocal microscopy (x, z) cut of the evolution of the gel film deflection of a 4%w casein and 4%w GDL aqueous dispersion when both top and bottom walls are coated with acrylamide brushes, Supplementary Table Ic. 10% of the proteins are fluorescently labelled. The scale bar is 100 µm (with a pixel size ratio of 1.3). Here the gel initially lies far from either top or bottom wall, H > H * . Poiseuille mechanism dominates the formation of the primary pattern. The total duration of the video is 181 min.

	Supplementary Videos Results Υ H/H * λexp λD nm 2 µm µm µm µm µm/s 10 8 Properties Measurements G α ξ h H1 H2 vmax mPa s Pa η mm mm mm mm λP λD+P 1.0 487 67000 400 1.0 487 67000 4.9 27.8 13.5 50.8 0.20 1.57 1.17 0.85 1.01 1.10 1.20 1.0 487 67000 4.9 32.7 30.7 42.9 0.14 2.54 2.23 1.52 1.24 1.85 1.90 1.0 487 67000 4.9 48.2 32.0 63.9 0.11 4.05 2.28 1.62 1.69 2.56 2.62 57.4 2.1 80.9 0.05 1.25 0.09 0.81 1.37 0.6 1.37 1.0 1817 2150 3.1 49.8 10.0 26.8 0.10 6.50 3.79 1.00 0.82 1.60 1.61 60.5 1.1 24.0 0.11 1.94 0.42 0.50 0.77 0.55 0.82 1.0 1817 2150 3.1 23.4 30.2 47.7 0.10 12.05 9.70 1.50 0.66 2.04 2.04 30.7 1.6 55.0 0.02 15.60 0.32 0.75 0.82 0.59 0.83 1.6 1142 24600 2.7 39.1 1.2 10.9 0.21 0.86 0.20 0.49 0.80 0.36 0.80 2.7 1858 25300 1.8 58.3 0.9 12.6 0.02 10.98 0.13 0.61 1.06 0.56 1.06 Video 1: Preparation # cas GDL gly e %w %w %w µm a 4 4 0 74.5 b 4 4 0 92.4 c 4 4 0 106.1 d 4 4 0 138.0 e 8 8 0 85.2 f 8 8 0 91.0 g 4 4 20 63.7 h 4 4 35 82.3 i 4 4 50 95.3 5.4 4660 5700 1

Acknowledgements Confocal experiments were conducted at SFR BioSciences Gerland -Lyon Sud (US8 / UMS3444). The authors thank José Bico, Arezeki Boudaoud, Cyprien Gay, L. Mahadevan and Olivier Pierre-Louis for theoretical insights at various stages of the research. Precious experimental ideas came from Catherine Barentin, Denis Bartolo, Elisabeth Bouchaud and Thibaut Divoux. Alan Parker at Firmenich provided the casein and GDL. Special thanks to Madame Biot who made M.A. Biot work accessible through the poronet website. ML thanks the Region Rhône Alpes and the Programme d'Avenir Lyon -Saint Etienne (PALSE NoGELPo) for postoctoral grant. TG and MN acknowledge funding from the Agence Nationale de la Recherche (ANR-11-PDOC-027). SM and ML acknowledge funding from the European Research Council under the European Union's Seventh Framework Program (FP7/2007-2013) / ERC grant agreement No. 258803 and from Institut Universitaire de France.

Author Contributions ML and MN conducted preliminary experiments. ML, SM and TG designed the experiments and interpreted the results. ML and TG performed quantitative experiments. ML performed data analysis and modelling. ML, SM and TG wrote the manuscript.

Author Information The authors declare that they have no competing financial interests. Correspondence and requests for materials should be addressed to TG (thomas.gibaud@ens-lyon.fr).