Adhesion dynamics of confined membranes

Adhesion dynamics of confined membranes

Tung B.T. To, Thomas Le Goff, and Olivier Pierre-Louis
Abstract

We report on the modeling of the dynamics of confined lipid membranes. We derive a thin film model in the lubrication limit which describes an inextensible liquid membrane with bending rigidity confined between two adhesive walls. The resulting equations share similarities with the Swift-Hohenberg model. However, inextensibility is enforced by a time-dependent nonlocal tension. Depending on the excess membrane area available in the system, three different dynamical regimes, denoted as A, B and C, are found from the numerical solution of the model. In regime A, membranes with small excess area form flat adhesion domains and freeze. Such freezing is interpreted by means of an effective model for curvature-driven domain wall motion. The nonlocal membrane tension tends to a negative value corresponding to the linear stability threshold of flat domain walls in the Swift-Hohenberg equation. In regime B, membranes with intermediate excess areas exhibit endless coarsening with coexistence of flat adhesion domains and wrinkle domains. The tension tends to the nonlinear stability threshold of flat domain walls in the Swift-Hohenberg equation. The fraction of the system covered by the wrinkle phase increases linearly with the excess area in regime B. In regime C, membranes with large excess area are completely covered by a frozen labyrinthine pattern of wrinkles. As the excess area is increased, the tension increases and the wavelength of the wrinkles decreases. For large membrane area, there is a crossover to a regime where the extrema of the wrinkles are in contact with the walls. In all regimes after an initial transient, robust localised structures form, leading to an exact conservation of the number of adhesion domains.

\setstretch

1.125 0.8cm

footnotetext: Institut Lumière Matière, UMR5306 Université Lyon 1-CNRS, Université de Lyon, 69622 Villeurbanne, Francefootnotetext: Instituto de Física, Universidade Federal Fluminense, Avenida Litorânea s/n, 24210-340 Niterói, RJ, Brazilfootnotetext: Aix Marseille Univ, CNRS, IBDM-UMR7288, 13009 Marseille, Francefootnotetext: Olivier.Pierre-Louis@univ-lyon1.fr

1 Introduction

Bilayer lipid membranes are abundant in biological systems 1. They are found in cell membranes, skin, eyes, articulations and pulmonary organs 2, 3, 4. Since their elasticity is dominated by bending rigidity, lipid membranes exhibit specific morphologies and dynamics, classifying their shape in a different class as compared to surface tension dominated phenomena, which govern the physics of capillarity and wetting. Helfrich has first proposed a model energy functional within which the behavior of lipid membranes could be explored 5. In the past two decades, much work has been devoted to the analysis of the consequences of the Helfrich energy on the morphology and dynamics of lipid membranes. Various successes include the shape of vesicles at equilibrium 6 or under hydrodynamic flow 7, 8, or the behavior of membrane stacks, and of supported membranes on substrates 9, 10, 11.

Here we wish to investigate the consequences of confinement on the dynamics of membranes. We study the dynamics of a fluid membrane with bending rigidity and area conservation confined between two attractive walls. This geometry is firstly motivated by the suggestion that membranes can experience a double-well adhesion potential in cell adhesion processes 12, or in biomimetic experiments 13. In these experiments, a short-range potential well located in the vicinity of the substrate results from the molecular binding of ligand-receptor pairs. In addition, a free energy barrier at intermediate ranges is provided by the entropic repulsion of a polymer brush grafted to the substrate, which mimics the glycocalyx 14. Finally, a long range attraction, resulting either from Van der Waals forces 12, or from gravity 13 enforces a second potential well for the membrane. In parallel with these works, other studies in the literature allow one to formulate a different picture for membrane adhesion, which is also based on the concept of a double-well potential. Indeed, blebbing 15, 16, 17 —the local detachment of cell membranes from the cytoskeleton, suggests that the link between membranes and the cytoskeleton can be described to some extent by simple adhesion concepts. Cell adhesion could therefore be mimicked by a competition between this adhesion of the membrane to the cytoskeleton, and adhesion to a substrate. Such a scenario is reminiscent of that proposed in Ref. 18, where membranes experience competing adhesion between a substrate and the cytoskeleton. Moreover, a third and somewhat different situation involving a double-well potential arises when two types ligand-receptor pairs with different lengths compete for adhesion, leading to two possible equilibrium separation distances, as proposed in Ref. 19.

Our study is inspired by these diverse pictures of membrane adhesion in two-state, or double-well potentials. Our aim here is to capture some generic features of these systems by exploring the simple case of a membrane confined between two flat adhesive walls. Adhesion in biological cells involving e.g. signaling, the remodeling of the cytoskeleton, or the diffusion and clustering of ligands and receptors, is certainly more complex than our minimal modeling approach. However, we hope that our results will provide hints to understand systems involving more physical ingredients.

On a more fundamental and theoretical level, our model for membrane dynamics defines a novel universality class for phase separation in two dimensions with unique features. Indeed, standard models for phase separation were developed to study spinodal decomposition in alloys 20, binary fluids 21, reaction-diffusion, magnetism and wetting phenomena. They are generically described by the time-dependent Ginzburg-Landau (TDGL) equation –also called the Cahn-Allen equation 22, or its conserved version the Cahn-Hilliard equation 23, and give rise to power-law coarsening (i.e. perpetual increase of the domain size) via curvature-driven motion of domain walls. Here, we show that the dynamical equations governing confined membranes share similarities with the Swift-Hohenberg equation, however with a time-dependent tension that enforces membrane area conservation. Within this model, membrane adhesion domains exhibit a transition to coarsening controlled by the total excess area of the membrane.

The results reported here build on our previous study of membrane adhesion dynamics based on a one-dimensional (1D) model with bending rigidity, but without area conservation 24, 25, 26. In 1D, we found that bending rigidity induces oscillatory interactions between domain walls. As a consequence of these oscillations, the dynamics freezes into a disordered or ordered profile depending on the permeability of the walls. In contrast to this behavior, we show here that in two-dimensions (2D) and without area conservation these oscillatory interactions between domain walls do not affect the coarsening behavior. Indeed, we recover standard coarsening with the same exponents as that of usual phase separation models (TDGL or Cahn-Hilliard). The study of such a model without area conservation can be motivated by the investigation of the coarsening dynamics of 2D systems at the Lifshitz-point, defined as the point where the prefactor of the gradient-squared term in the Landau free energy density vanishes 27, and where higher-order bending-like squared-Laplacian terms come into play.

However, the present study focuses on the case of 2D membranes where local area conservation should be imposed. As discussed above, this leads to the presence of a time-dependent tension in the dynamical equations. The numerical solution of these equations reveals three regimes, hereafter denoted as A, B and C, depending on the excess area of the membrane.

For small excess area (regime A), the membrane freezes into a state with flat adhesion patches of finite size. This freezing can be understood as follows. Domain wall motion is driven by an effective positive wall tension, and acts so as to reduces the total domain wall length. Due to area conservation, domain wall length decrease implies an increase of the membrane excess area per unit length in domain walls. This increase induces a cancellation of the domain wall tension, leading to the arrest of domain wall motion. The dynamics in regime A can be analyzed within a simple model for the coupled dynamics of domain wall motion and of the nonlocal membrane tension. Since it corresponds to the cancellation of wall tension, the asymptotic value of the nonlocal tension is equal to the threshold tension for linear instability of domain walls in the SH equation.

For intermediate excess area (regime B), the membrane exhibits coarsening with a coexistence of flat and wrinkled domains. The wrinkle phase forms spontaneously when domain walls collide. The fraction of the system occupied by the wrinkle phase reaches a constant value at long times, which increases linearly with the excess area. Concurrently, the size of the flat and wrinkled domains increases indefinitely with time. Since the system reaches coexistence of wrinkles and flat domains, the asymptotic value of the nonlocal tension corresponds to the threshold for nonlinear stability of domain walls in the SH equation.

For larger excess area (regime C), the wrinkle phase invades the whole system and the membrane freezes into a labyrinthine wrinkle pattern. The amplitude and wavelength of the wrinkles are analyzed within a simple sinusoidal ansatz, that reveals the presence of two different regimes within regime C. For large excess area the amplitude of the wrinkles is fixed by the contact with the two walls, while wrinkles with smaller excess area exhibit a free amplitude smaller than the distance between the two walls.

We focus on the case with permeable walls, but also briefly report on the behavior of a simplified model for impermeable walls which exhibits similar dynamics. We also find unexpectedly that the number of adhesion domains is strictly conserved in all cases during the late stages of the dynamics due to the formation of robust localised structures which forbid the complete disappearance of the adhesion domains.

2 Model Equations

Fig.  1: (Color online) Membrane between two adhesive substrates. (a) 3D view; (b) 2D schematics: left black curve: double-well potential, right red curve: a slice of the membrane profile (profiles are obtained from simulation with at , corresponding to the last figure in regime B2 of Fig. 2).

A schematic representation of a membrane of height along confined between two parallel flat walls located at is shown in Fig. 1(a,b). The interaction between the membrane and the walls is modeled via a double-well potential . Adding this interaction energy with the Helfrich bending energy, we obtain the total energy of the membrane as

(1)

where denotes the infinitesimal area element of the membrane, denotes the membrane local curvature and denotes the bending rigidity.

In addition, the local conservation of the membrane area can be expressed as a conservation of the membrane area density :

(2)

where is the membrane velocity parallel to the walls, is the gradient in the (x,y) plane. In order to enforce this constraint, we make use of a local space-dependent and time-dependent Lagrange multiplier (which can be interpreted as a local membrane tension 28, 29). The membrane energy is thus generalised as

(3)

The local tension leads to an additional contribution to membrane forces, which is constrained to obey Eq. (2) at all times. Such a constraint allows one to determine .

Since all the phenomena described here occur at small scales, we consider the Stokes regime where inertia is negligible. Hence, membrane forces resulting from energy variations and inextensibility have to balance viscous forces exerted by the surrounding liquid on the membrane:

(4)

where is the stress tensor in the liquid, denotes the liquid above or below the membrane, and is the normal to the membrane. The liquid with velocity obeys the incompressible Stokes equations, with and , where where is the fluid viscosity and its pressure. In addition, we assume no-slip 30, 31 and impermeability 32, 33, 34 at the membrane, leading to

(5)

In order to account for the permeability of the walls, we impose

(6)

where is a reference constant pressure, and is the wall permeability. Finally, we assume a no-slip condition at the walls .

3 Lubrication limit: permeable walls

We apply the lubrication limit where is small. In this limit, the hydrodynamic flow above and below the membrane are simple Poiseuille flows parallel to the walls. Using the boundary conditions presented in the previous section, the hydrodynamic flows can be obtained explicitly. We then obtain dynamical equations for the membrane profile using Eq. (5).

The derivation of the evolution equation for in this limit actually follows the same steps as the simpler case of a one-dimensional membrane without area conservation reported in Ref. 24. The main difference comes from the constraint of membrane incompressibility. The details of the calculations and the general result are reported in Appendix A. We will mainly focus on the limit of large normalised permeability, , with

(7)

where is the amplitude of the interaction potential.

In physical systems, the value of depends crucially on the permeability of the substrate. One possibility to evaluate is to use Darcy’s law for porous media, using 35 , where is the scale of the pores and the thickness of the wall. This suggests . In the case of scaffolded actin cytoskeleton, the typical pore size is m and the thickness m. Moreover, assuming that the energy-scale of the potential is dictated by the cytoskeleton-lipid membrane adhesion energy, we find 36 . Using 37, 38 and , we find . For collagen, a very common extracellular matrix, the typical pore size is m and thickness is m 39 which gives same order of magnitude for .

However, if the substrate is covered by another lipid membrane, the permeability of the substrate becomes very small  32, 33, 34. An increase of about 1 atm of osmotic pressure induces a speed of about m.s for the water across a lipid membrane 32, 33, 34. This means that from Eq. (6) m.s.kg. Then using water viscosity, and the same values as above , m and J.m, we find . We will briefly discuss the case of small permeabilities in the end of the paper.

In the limit of strongly permeable walls , the dynamical equation reads

(8)
(9)

Here, denotes the Laplacian operator in planar coordinates . One important and non-trivial result emerging from the lubrication limit is that, to leading order, the local Lagrange multiplier is constant in space, with a value , thereby leading to a nonlocal area conservation constraint.

The evolution equations Eqs. (8,9) decrease the energy, here written in the small slope limit as

(10)

Indeed, as shown in Appendix B using the Schwarz inequality, we always have .

4 Normalization and numerical methods

In order to discuss the dynamical behavior of Eqs. (8,9) and to solve the equations numerically, we normalise space and time, leading to

(11)
(12)

where , , , , , with , and . Using the numerical values of the previous section, we obtain the order of magnitude of tensions J.m, and the typical lengthscale parallel to the membrane nm. The small slope limit amounts to considering is small as compared to 1. Although these slopes are not very small, the lubrication limit is known to provide a qualitatively good and robust description of the physical behavior for moderate slopes 40. In addition, smaller adhesion energies, e.g., using Jm as suggested by Ref.16, or Jm in Ref.13, lead to even smaller slopes. In general, the adhesion energy, on which the validity of this limit depends crucially, is system-dependent and varies strongly with the type and the density of binders. Considering the different case of physical adhesion with a porous solid substrate, a crude approximation consists in multiplying the physical adhesion potential of Ref.41 —based on Van der Waals and hydration interactions, with the solid fraction . Assuming for example , we then find a value Jm similar to those reported above. Summarizing this discussion, we expect in most cases to nm for nano-confinement with nm.

We choose a specific form of the double-well adhesion potential

(13)

where

(14)

is a short-range repulsion, the aim of which is to avoid collision of the membrane with the walls. In the simulations, we have chosen , . Assuming that nm, the position of the minimum of the potential corresponds to a distance to the substrate nm, which is similar to those reported in the literature for physical interactions41, or binders 13. During the simulations, we usually assumed that . However, we used to accelerate long simulations in the regimes where the membrane did not approach the walls.

The dynamics is integrated by means of a first-order exponential-time-differencing method (ETD1) in Fourier space 42 with space and time discretization bins , and . We have modified the integration scheme in order to conserve exactly the excess area

(15)

Our scheme amounts to choosing a specific discretization of the expression of the tension in Eq. (12) in order to enforce exact area conservation. The details of this scheme are presented in Appendix C.1.

The simulation box sizes where usually , or corresponding to physical sizes from to m.

5 Area-preserving vs tensionless membranes

We start with noisy initial conditions (details on these conditions are described in Appendix C.2). Three different regimes are obtained depending on the excess area density

(16)

where is the system size in the plane (if the simulations are performed in a rectangular box of size , we have ). For increasing excess area density, we first find a regime with frozen flat domains, then a regime with coarsening and with coexistence between flat domains and wrinkles, and finally a regime with a disordered pattern of frozen wrinkles. These three regimes will be denoted as regimes A, B and C in the following. Some snapshots of these evolutions are shown in Fig. 2.

Fig.  2: (Color online) Membrane adhesion dynamics for various excess areas. Yellow: adhesion patches on the upper wall. Dark blue: adhesion patches on the lower wall. Regime A with freezing of flat domains for small excess area : . Regime B with coexistence of the flat-domain phase and wrinkle phase with coarsening for intermediate excess area: B1 with , and B2 with . Regime C with frozen wrinkles for larger excess area . System size .
(a)
(b)
Fig.  3: (Color online) (a) TDGL4 (Eq. (8) with ), , . (b) CH4 (constant mobility) (Eq. (35a) with and ), , .

Since frozen states were also observed in a one-dimensional model without area conservation in Ref. 24, we wish to investigate the behavior of the model without the constraint of area conservation as a preamble to the analysis of the full model. This is done by simulating Eq. (11) with . The resulting equation was called the fourth order Time-dependent Ginzburg-Landau equation (TDGL4) in Ref. 24. Such an equation corresponds to a system that would exhibit bending rigidity, but for which the extension of the area of the interface would occur at no cost (physically corresponding to dynamics at a Lifshitz point).

Frozen disordered states obtained in simulations of TDGL4 in 1D 24 can be seen as a consequence of trapping of domain-walls (called kinks in 1D) into their mutual oscillatory interactions. These oscillations, which can be traced back to the oscillatory tails of the domain wall profiles, are still present in 2D. The first oscillation can indeed be observed in the vicinity of all domain walls in Fig. 1. However, the TDGL4 dynamics in 2D, shown in Fig. 3(a), actually leads to a simple coarsening behavior with flat adhesion patches, the size of which grows like . This is identical to the 2D behavior of the TDGL equation. Such a coarsening behavior is usually interpreted as a consequence of motion of domain walls driven by their curvature. Thus, the coarsening behavior of TDGL4 suggests that, in absence of area-conservation constraint, motion by curvature of domain walls dominates over oscillatory interactions.

As mentioned above, membranes with area conservation exhibit strikingly different dynamics. Hereafter, we discuss the three dynamical regimes arising in membranes with area conservation.

6 Small excess area

Fig.  4: (Color online) Time-evolution of the lengthscale and the tension . Size of the simulation box . (a) and (c): results in the limit of large wall permeability. (b) and (d): results for impermeable walls. In (a) and (b), thick lines correspond to , and thin lines to . Excess area: ( ) ; ( ) ; () ; ( ) ; ( ) .

Regime A corresponds to simulations of Eqs. (11,12) at small excess area. In this regime, flat adhesion domains expand initially (Fig. 2A, ), and the tension decreases to negative values. Later, these domains freeze (2A, ), and reaches a constant negative value independent of initial conditions and of the imposed excess area. The evolution of the domain size and of the tension are shown in Fig. 4(a,c).

When is constant in time, our dynamical equation (11) is identical to the much-studied Swift-Hohenberg (SH) equation 43, 44, 45, 46, 47 (neglecting the contribution to the potential). Hence, the steady-states of our model are steady-states of the SH equation. However, the stability of these steady-states can be different. We observe that the value towards which converges in our simulations is close to the limit of linear stability of flat domain walls with respect to transverse perturbation in the SH model at as reported by Hagberg et al.47. Such a limit of stability is associated with the cancellation of the energy of domain walls in the SH equation.

In order to explore the consequence of this cancellation, we need to relate the energy in the two models. In normalised form, they read:

(17)
(18)

The energies and per unit length of a straight and isolated domain wall in our model and in the SH model therefore read

(19)
(20)

where is the coordinate orthogonal to the domain wall. Assuming that the whole excess area is stored in domain walls, we have

(21)

where is the length of domain walls in the system, and

(22)

is the excess area per unit domain wall length. Combining Eqs. (19,20,21), we find

(23)

Since when , we obtain an expression for the size of the domains in the frozen state of our model

(24)

where is the value of for .

Using independent 1D simulations of the SH equation in a periodic box with two opposite domain walls, we have determined as a function of in steady-state, as seen in Fig. 5. In particular, we find . Inserting this value in Eq. (24), we obtain a prediction of in good agreement with numerical results, as shown by the left red curve on Fig. 6.

Fig.  5: (Color online) Domain wall energy and excess area. Continuous curve: . Dashed curve: . Dotted-dashed curve: . The vertical lines indicate (left red), and (right blue).
Fig.  6: (Color online) Typical length scale at long times as a function of for . Symbols: for , and for . The critical values of the normalised excess area and are indicated with dashed black lines (these lines are doubled to indicated the accuracy of our measurement of the transition). Thin red curve and thin magenta curve: analytical results. Inset: linear dependence of on . Symbols for simulation results . Magenta line: .

We now discuss why the tension converges to the special value . In order to investigate this point, we make use of an effective model for domain wall motion. The derivations are inspired from that of Ref. 25, and some details are reported in Appendix D. As expected, we find that domain wall motion can be driven either by wall-wall interactions or by curvature, and the local normal velocity of the wall obeys

(25)

where is the local domain wall curvature, and is an interaction term (see Appendix D for its detailed expression). As discussed above, interactions between two domain walls are known to be oscillatory 25, and to decay exponentially. They become negligibly small when wall-wall distances exceed a few domain wall thicknesses.

In regime A, domain walls are seen to be far apart, and as a consequence, domain wall motion should be mainly driven by curvature. Such curvature-driven domain wall motion leads to changes in the system configuration that reduce the total length of the domain walls when , i.e. when from Fig. 5. Indeed, using the geometric relation , and inserting the expression of from Eq. (25) neglecting the interaction term, we find . Since the total excess area is conserved, this decrease of leads to an increase of the excess area stored per unit length in domain walls. Such a change in is intuitively associated to a decrease of the tension toward more negative values. Indeed, we expect , i.e., more positive tensions correspond to pulling the membrane out from the domain walls which decrease the excess area inside the domain wall, while more negative tensions correspond to pushing towards the domain walls leading to an increase of excess area. The inequality is confirmed by 1D simulations of the SH equation in steady-state reported in Fig. 5.

As a summary, curvature-driven wall motion leads to an increase of accompanied by a decrease of the tension . This decrease is governed by the equation

(26)

where is the arclength along the domain walls, and the integration runs over all domain walls. This relation is a consequence of area conservation, and its derivation is reported in Appendix D. Once again, we neglect the interaction terms proportional to . As seen in Fig. 5, is an increasing function of . Since is proportional to , and recalling that , Eq. (26) shows that will decrease up to the point where , where it reaches the constant . Since , the motion of domain walls freezes from Eq. (25).

As already mentioned above, this discussion is based on the assumption that interactions between domain walls are weak. Since interactions decrease exponentially with the distance, such an assumption should be valid in the limit where the distance between domain walls is large enough. However, as increases, this distance decreases, and interactions between walls should become relevant. Indeed, a different regime, discussed in the next section and hereafter denoted as regime B appears for larger .

Fig.  7: (Color online) Comparison of the energies of an isolated wrinkle, and of two domain walls. Left red vertical line: . Right blue vertical line: . Solid line: . Dashed line: . Dashed-dotted line: .

7 Intermediate excess area

As announced in the previous sections, a different regime is observed for larger excess area. This regime, denoted as regime B is found for , with

(27)
(28)

The errors on and are based on the difference between the closest upper and lower bounds for the transition obtained by long-time simulations in a system of size .

In regime B, both flat and wrinkled domains coexist and expand perpetually (see regimes B1 and B2 in Figs. 2). Numerical results for the dynamics of the tension and of the size of the domains are shown in Figs. 4(a) and 4(c), with the symbols ( ), ( ), (), ( ) and ( ) in the order of increasing . Regime B corresponds to the symbols (), ( ) and ( ). The typical sizes of flat domains and the size of wrinkle domains both increase with time. The method for extracting these lengths from simulation data is discussed in Appendix C.4. The growth rates of and decrease with increasing . However, no universal exponent related to this coarsening process could be found within our simulations. In addition, we observe that the fraction of the system covered with the wrinkle phase reaches a constant at long times, as shown in Fig. 8(a).

As seen in Fig. 4(c), the tension in regime B first decreases to a minimum value close to or larger than , and then increases to another constant asymptotic value. This value, hereafter denoted as , corresponds to a point of coexistence of flat and wrinkled states in our system. Following the same procedure as for the definition of generalised thermodynamic potentials, coexistence between two states correspond to the equality of Legendre-transformed energy with respect to the Lagrange multiplier conjugate to the fixed quantity . Coexistence therefore corresponds to the point where the -energy densities of the flat and wrinkled states are equal: , where the indexes "FD" and "wr" respectively correspond to the flat domains regions and wrinkles regions.

Note that the energy a priori contains not only the contribution of flat domains with , but also that of domain walls inside this region. However, the density of domain walls is low, and we shall only consider the contribution of flat domains, for which (such a cancellation relies on the assumption that the minimum is far enough form the walls for the short-range potential to be negligible at , i.e., ). Following the same lines, the energy within the wrinkle phase contains a contribution due to wrinkle bending and wrinkle defects, which is also neglected. We therefore end up with a simpler condition of coexistence , calculated for a periodic phase of parallel rolls. This is exactly the condition of nonlinear stability of domain walls with respect to the formation of a wrinkle phase in SH as discussed in Ref. 47, leading to in agreement with our simulation results.

Using this concept of coexistence, together with the observation of the variation of the tension with time, we propose a scenario in three stages for regime B. First, the dynamics looks like that of regime A, and the evolution is dominated by domain wall motion, leading to a decrease of toward . Second, once , the flat states become unstable with respect to the formation of the wrinkle state. In simulations at low excess area, i.e., close to regime A, this instability usually occurs as follows: two domain walls collide and form an isolated wrinkle, which grows by zipping more domain walls. This is consistent with the inequality , where is the -energy per unit length of a straight and isolated wrinkle, which can be checked on Fig. 7. For larger excess area, the separation between the two initial stages is less clear since many domain walls are already close to each other initially. Finally, in the third stage, coarsening occurs with coexistence of wrinkle domains and flat domains. It is tempting to speculate that the coarsening process is driven by the decrease of the length of frontiers between the coexisting wrinkled state and the flat state. However, annihilation of defects within the wrinkled state and motion of simple domain walls between flat domains could also play a role.

When the excess area is increased, the fraction of the system covered by wrinkles increases. When the excess area exceeds a threshold value, the full system is covered by wrinkles, and a different regime, denoted as regime C is found.

(a)
(b)
Fig.  8: (Color online) Time-evolution of (a) the area fraction covered by wrinkle domains for (top three curves) and (bottom three curves). (

) . () . () . ( ) . ( ) . () . (b) Area fraction covered by wrinkle domains at long simulation times. () asymptotic value extracted from simulations. Blue line: Theoretical prediction from Eq. (29).

Some analysis of the fraction of the system occupied by wrinkles in regime B is possible assuming once again that the contribution of defects such as domain walls in the flat phase and defects in the wrinkle phase are negligible. The wrinkle state is then composed of parallel rolls of wavelength where the superscript indicates that this quantity is evaluated for . We define wrinkle length , i.e., the total wrinkle length summed over all wrinkles (formally, this can be defined as, e.g., the total length of all the lines of local maximums of the wrinkles in the whole system). Assuming that all the excess area is stored in the wrinkle phase, one has . The area covered by wrinkles then reads . Combining these relations, we find that the fraction of the system covered by the wrinkle state is proportional to the normalised excess area :

(29)

In order to determine the quantities appearing in the right hand side of Eq. (29), we simulated the SH equation with in 1D in a box of size . We found that a stable wrinkle profile is reached with 66 or 67 wavelengths. Taking the average of these values, we obtain , , and . These numerical results in 1D confirm that in the roll phase vanishes at . Indeed, . In addition, we find . Using this result in the prediction Eq. (29) provides good agreement with simulations, as seen in Fig. 8(b).

Assuming that the transition to the frozen wrinkle state at large excess area corresponds to the filling of the system with the wrinkle state , we obtain a prediction for the critical excess area corresponding to the transition to regime C:

(30)

Using the numerical values reported above, we find which is very close to the value of observed in the simulations, Eq. (28).

8 Large excess area

Regime C is obtained for . In this regime, after some transient dynamics, the adhesion pattern is frozen into labyrinths as shown in Fig. 2C. Such labyrinthine patterns have also been observed in the case of the SH equation 46.

In our simulations of regime C, the tension again converges to a stationary value. However, as opposed to the previous cases in regimes A and B, the asymptotic tension depends on . Furthermore, we observe in the simulations that, as the excess area increases in regime C, the width of the wrinkles decreases, whereas the maximum height and the variance of the membrane profile increase.

In order to analyze this behavior, we perform the changes of variables , and . The excess area is then written as

(31)

where the integral in the right hand side only depends on the shape of the wrinkle profile. An inspection of Eq. (31) shows that, if the shape of the wrinkle does not vary much, the quantity should be a constant. Indeed, the relation is in good agreement with numerical results. Using the value of measured in simulations, this expression provides the dashed magenta line on the main plot of Fig. 6.

Assuming that the profile of the wrinkles is sinusoidal , we find , and , suggesting that

(32)

The sine ansatz therefore leads to an overestimation of the constant, showing that the profile of the rolls is different from a simple sinusoidal profile.

Furthermore, the asymptotic value of the tension is controlled by the balance between the negative tension which tends to store excess area by forming the wrinkles, and the bending rigidity which tends to flatten the membrane. A simple balance between the two terms in the energy Eq. (18) suggests that , or . Simulation results reported in the inset of Fig. 6 indicate that this relation is in fair agreement with the numerical results with .

If we assume again a sine profile for the membrane rolls, the potential energy density is not affected by the wavelength. The contributions that depend on the wavelength are the bending energy and the tension contribution. Minimizing the energy with respect to the wavelength, we find

(33)

This relation between the tension and the wavelength appears in better agreement with the simulations as compared to Eq. (32), but is still not very accurate.

In the following, we proceed further with the sine ansatz to predict the wavelength, amplitude and tension from a direct minimization of the energy density. The details of the derivations are reported in Appendix E. Depending on the excess area, we find two regimes. In both regimes Eqs. (32,33) are valid. First, when the amplitude of the sine profile is small enough, the amplitude and the wavelength are both varied to minimise the energy density. However, for large enough excess area, the extrema of the sine profiles touch the walls. In this case, the short range repulsion potential prevents the membrane from crossing the wall. We therefore fix the amplitude and minimise the wavelength only. The crossover to this wall-contact regime occurs for with

(34)

The details of the calculations are reported in Appendix E.

The results of the sine-profile ansatz, shown on Fig. 9, reproduce the trends obtained with the full simulations in regime C.

Fig.  9: (Color online) (a) Tension, (b) amplitude and (c) wavelength of the wrinkles in regime C. Red dots represent the simulation results. The solid blue lines are the result of the sine-profile ansatz. The dashed and dotted vertical lines respectively indicate the start of regime C at , and the crossover to the wall-contact regime at .

Note that the above value of corresponds to a crossover rather than a sharp transition in the simulations. In addition, our convenient decomposition of the potential into a smooth double-well potential and a sharp short-range repulsion is valid only when the minimum of the potential is far enough from the walls. In contrast, when is close to , the short range repulsion affects the membrane profile in the potential well. One consequence of this is the inconsistency of Eq. (34) when . Our choice in simulations however corresponds to the regime where Eq. (34) is valid.

9 Impermeable walls

Fig.  10: (Color online) Membrane adhesion dynamics for impermeable wall for various excess areas. Yellow: adhesion patches on the upper wall. Dark blue: adhesion patches on the lower wall. Regime A with freezing of flat domains for small excess area : . Regime B with coexistence of the flat-domain phase and wrinkle phase with coarsening for intermediate excess area: B1 with , and B2 with . Regime C with frozen wrinkles for larger excess area . System size .

In this section, we discuss the limit of very impermeable walls with . The full lubrication equations are reported in Appendix A. These equations include not only a term accounting for the conservation of the total flow, as already found in a one-dimensional model 24, but also a space-dependent tension that drives tangential forces along the membrane due to local area conservation. For the sake of simplicity, we neglect these two terms, and discuss a simplified model which is a straightforward transposition of Eq. (8,9) to the case of conserved dynamics:

(35a)
(35b)

The nonlinear mobility 24