Stochastic elastohydrodynamics of contact and coarsening during membrane adhesion

Abstract

Contact between fluctuating, fluid-lubricated soft surfaces is prevalent in engineering and biological systems, a process starting with adhesive contact, which can give rise to complex coarsening dynamics. One representation of such a system, which is relevant to biological membrane adhesion, is a fluctuating elastic interface covered by adhesive molecules that bind and unbind to a solid substrate across a narrow gap filled with a viscous fluid. This flow is described by the stochastic elastohydrodynamic thin film equation, which incorporates thermal fluctuations into the description of viscous nanometric thin-film flow coupled to elastic membrane deformation. The average time it takes the fluctuating elastic membrane to adhere is predicted by the rare event theory, increasing exponentially with the square of the initial gap height. When the forces arising from spring-like adhesive molecules are included in the simulations, thermal fluctuations initiate phase separation of domains of bound and unbound molecules. The coarsening process of these unbound pockets displays close similarities to classical Ostwald ripening; however, the inclusion of hydrodynamics affects power-law growth. In particular, we identify a new bending-dominated coarsening regime, which is slower than the well-known tension-dominated case.

1. Introduction

Adhesion between soft fluctuating surfaces is found in a range of engineering applications and in biological systems. Complex lifeforms rely on cells adhering to each other, facilitated by the binding of membrane-anchored adhesive molecules across the gap between the membranes of the cells. The dynamics of adhesion involves a rich interplay of membrane deformation, chemical kinetics of the aforementioned molecular binding, as well as fluid flow in the narrow space between the membranes, which can give rise to intricate dynamical processes such as the coarsening of adhesion patches. A number of essential physiological processes depend on adhesion dynamics, such as cadherin-mediated adhesion (van Roy & Berx 2008), lumen formation between cells in a growing embryo (Dumortier et al. 2019) and the immune synapse which facilitates an immunological response (Grakoui et al. 1999; Dustin & Cooper 2000). Experiments and numerical results have shown that the features of these phenomena can be passive processes only relying on the physical forces involved (Qi, Groves & Chakraborty 2001; Dustin 2010; Le Verge-Serandour & Turlier 2021).

One way to describe the adhesion of membranes is to adopt the Helfrich model with an adhesion potential (Seifert & Lipowsky 1990; Smith et al. 2003, 2008; Agrawal 2011; Hill & Al-Amodi 2024). This approach successfully predicts the equilibrium shape of the membrane, but neglects the motion of the extracellular fluid in the cleft between the two membranes. In such nanometrically thin but micrometrically wide channels, however, the forces required to squeeze the viscous extracellular fluid are not negligible, motivating a viscous thin film description of the flow (Leong & Chiam 2010). The lubrication theory allows for a relatively simple inclusion of forces due to membrane deformation, molecule/protein binding, as well as thermal fluctuations (Carlson & Mahadevan 2015a , b ). Membranes resist lateral deformation due to membrane tension $\gamma$ (Evans, Waugh & Melnik 1976), as it increases their free energy in a manner analogous to interfacial surface tension, whose influence on the dynamics of thin liquid films has been studied extensively (Oron, Davis & Bankoff 1997; Craster & Matar 2009), particularly in the case of the dewetting of nanoscale thin liquid films (Zhang & Lister 1999; Mecke & Rauscher 2005; Grün et al. 2006; Nguyen et al. 2014; Zhang, Sprittles & Lockerby 2019; Zhao et al. 2022). In addition, membranes of finite width $d$ resist bending (Evans 1974), as one side is compressed and the other side is stretched, which is characterised by a bending modulus $B=Ed^3/12(1-\nu ^2)$ , where $E$ is the Young’s modulus and $\nu$ is the Poisson ratio. Both tension and bending resist the deformation of a flat membrane, but there are subtle differences in how they act. A bendocapillary length $l_{BC}=\sqrt {B/\gamma }$ can be obtained by balancing the forces from tension and bending, which is approximately $100$ nm for the properties of most cell membranes, and bending dominates for length scales smaller than the critical length scale (Roman & Bico 2010; Deserno 2015).

For adhesion to occur, the membranes must first come close enough to each other to allow the adhesive molecules to start to form bonds. In the absence of directed motion due to active cytoskeletal forces or protein–membrane interactions (Liese & Carlson 2021; Yuan et al. 2021), the forces required to push out the extracellular fluid in the channel can be attributed to thermal fluctuations (Aarts 2004; Delgado-Buscalioni, Chacon & Tarazona 2008) as well as the other fluctuations inherent to living matter (Guo et al. 2014; Gupta & Guo 2017). The fluctuations in the width of the channel, although small in amplitude and random in direction, can still bring membranes close enough to initiate adhesive molecule/protein binding if given sufficient time. This process is similar to the spontaneous thermal dewetting of a linearly stable thin viscous film coated on a solid substrate, where thermal fluctuations are needed to bring the liquid free surface close enough to the solid substrate for the disjoining pressure to rupture the film, for which the average waiting time for rupture can be predicted by rare-event theory (Sprittles et al. 2023; Liu, Sprittles & Grafke 2024).

After the initial binding of adhesive molecules, the adhesion patches grow in size. If more than one type of adhesive molecule is involved, the adhesion patches can separate into distinct regions (phases) where one specific type of molecule is bound, as is seen in the immune synapse, receptor tyrosine kinases (Janes, Nievergall & Lackmann 2012; Case, Ditlev & Rosen 2019) and simulation of membrane adhered to heterogenous substrate (Shelby et al. 2023). Excess fluid from the adhesion patches is squeezed into the regions where molecules are unbound, which further increases the distance between membranes, creating pockets of fluids known as lumens, which can be observed in mammalian embryo development (Dumortier et al. 2019). These lumens can also be reproduced in reconstituted systems by applying an osmotic shock to a giant unilamellar vesicle (GUV) that is adhered to a supported lipid bilayer (Dinet et al. 2023). In these systems, the separated phases often undergo a passive, physically driven coarsening process where small patches diminish, while larger patches grow. These dynamics are reminiscent of other phase separation processes occurring in the cooling of metal alloys (Allen & Cahn 1979; Komura et al. 1985; Bray 1994; Livet et al. 2001), liquid–liquid phase separation in biological systems (Su et al. 2024), as well as in droplet aggregation when a thin liquid or polymer film is adhered to a solid substrate by an attractive potential due to intermolecular forces (Derrida, Godrèche & Yekutieli 1991; Otto, Rump & Slepcev 2006; Gratton & Witelski 2008; Lal et al. 2020). In those systems, an effective interfacial tension drives the coarsening process, providing a well-established $t^{1/3}$ power law (Bray 1994) for the growth the characteristic length scale $L_c$ . For the coarsening observed in biological membranes, while membrane tension could give rise to the same $1/3$ power law, the effect of membrane bending has not been studied and will potentially introduce different dynamics.

In this article, numerical simulations of lubrication flow including thermal fluctuations are used to study membrane adhesion, with particular emphasis on the effects of membrane bending as opposed to analogous systems involving interfacial tension. The average waiting time for adhesion to occur due to fluctuations is predicted using rare-event theory. We then show that membrane bending gives rise to power-law coarsening behaviour that is slower than for tension-driven coarsening. We also demonstrate that the strength of thermal fluctuations does not significantly alter the coarsening rate, but that the initial height of the membranes does so in the hydrodynamic regime.

2. Mathematical model and numerical methods

To simultaneously study the fluid flow, elastic membrane bending, protein-like binding and stochastic fluctuations during membrane adhesion, we turn to the elastohydrodynamic thin film (Williams & Davis 1982; Carlson & Mahadevan 2016), describing the physical scenario of figure 1 $(a)$ , i.e. a thin layer of viscous liquid confined by an elastic membrane. For simplicity as well as reflecting in vitro experimental conditions in which a GUV interacts with a supported lipid bilayer (Fenz et al. 2017; Dinet et al. 2023), the substrate is static while the upper membrane ‘rests’ on the viscous fluid film at height $\hat {h}(\hat {x},\hat {y},\hat {t})$ (where the hats indicate dimensional variables) and can move vertically. As in the biological and synthetic systems described previously (Dumortier et al. 2019; Dinet et al. 2023), the film height $\hat {h}$ (which is typically in the range of tens of nanometres) is much smaller than the length of the domain in the lateral directions, $L$ (typically several microns). Across the gap, membrane-bound proteins may bind or unbind to the solid surface.

Figure 1. $(a)$ A sketch of an elastic membrane with thickness $d$ in close proximity to a rigid wall separated by a thin layer of viscous fluid of height $\hat {h}(\hat {x},\hat {y},\hat {t})$ . Membrane molecules may bind across the channel only if their distance is below a critical value $h^*$ . $(b)$ Left: a contour map of the non-dimensional height profile $h(x,y,t)$ at the time $t^*$ , which is the onset of adhesion between the membrane and solid surface. Data are shown for a membrane with non-dimensional initial height $h_0=1.6$ and fluctuation intensity $Q_B=( {1}/{l})\sqrt {{2 k_B T }/{B^{1/2}(\kappa c_0)^{3/2}}}=0.01$ , with $l$ the equilibrium length of the adhesive molecules, $k_BT$ the thermal energy, $B$ the bending stiffness, $c_0$ the equilibrium concentration and $\kappa$ the molecule spring stiffness coefficient. Right: a cross-section of the profile along the red line in the contour map; the film height at the point of contact drops to a non-dimensional value of $h=\hat {h}/l\approx 1$ , matching the equilibrium length of the adhesive molecules. $(c)$ Left: a contour map at a later time when most of the membrane is bound, but liquid is collected in unbound patches. Right: a cross-section of the profile along the red line in the contour map, illustrating the formation of blisters (regions where the adhesive molecules are unbound and $h=\hat {h}/l\gtrapprox 1$ ) during the coarsening process.

2.1. Stochastic thin film equation

By assuming a small aspect ratio of the viscous channel, $\hat {h}/L\ll 1$ , one can apply the lubrication approximation to describe the flow of a Newtonian fluid with viscosity $\mu$ in the channel, leading to a parabolic velocity profile (Batchelor 2000). A random stress tensor is introduced in the momentum equations to account for the thermal fluctuations in the fluid, which are significant on the relevant nanometric scale (Landau & Lifshitz 1987). By imposing no-slip and kinematic boundary conditions at the membrane, one arrives at the following stochastic thin film (Davidovitch, Moro & Stone 2005; Mecke & Rauscher 2005; Grün et al. 2006; Durán-Olivencia et al. 2019) that can be used to describe $\hat {h}(\hat {x},\hat {y},\hat {t})$ :

(2.1) \begin{equation} \frac {\partial \hat {h}(\hat {x},\hat {y},\hat {t})}{\partial \hat {t}} = \hat {\boldsymbol{\nabla }}\boldsymbol{\cdot }\left (\frac {\hat {h}^3(\hat {x},\hat {y},\hat {t})}{12\mu }\hat {\boldsymbol{\nabla }} \hat {p}(\hat {x},\hat {y},\hat {t}) \right) + \sqrt {\frac {k_B T}{6\mu }}\hat {\boldsymbol{\nabla }} \boldsymbol{\cdot }\left (\hat {h}^{3/2}(\hat {x},\hat {y},\hat {t})\boldsymbol{\hat {\eta }}(\hat {x},\hat {y},\hat {t})\right) \end{equation}

where $\hat {\boldsymbol{\nabla }}$ represents the two-dimensional (2-D) gradient operator $(\partial /\partial \hat {x},\partial /\partial \hat {y})$ and the dependent variable $\hat h$ represents height in the third spatial direction.

The first term on the right-hand side of (2.1) represents the change in film height due to a lateral flux driven by the horizontal pressure gradient $\hat {\boldsymbol{\nabla }} \hat {p}$ . The second term on the right-hand side incorporates the effect of thermal fluctuations at temperature $T$ , by introducing a random flux in accordance with the fluctuation–dissipation theorem and averaged across the channel (Davidovitch et al. 2005; Mecke & Rauscher 2005; Grün et al. 2006). Here, $\boldsymbol{\hat {\eta }} (\hat {x},\hat {y},\hat {t})$ , is a random vector with Gaussian white noise uncorrelated in both time and space, i.e. $\langle \hat {\eta }_i(\hat {x},\hat {y},\hat {t})\rangle =0$ and $\langle \hat {\eta }_i(\hat {x},\hat {y},\hat {t})\hat {\eta }_j(\hat {x}^{\prime},\hat {y}^{\prime},\hat {t}^{\prime})\rangle =\delta _{ij}\delta (\hat {x}-\hat {x}^{\prime})\delta (\hat {y}-\hat {y}^{\prime})\delta (\hat {t}-\hat {t}^{\prime})$ , where $\langle \; \rangle$ is the ensemble average, $\delta _{ij}$ is the Kronecker symbol, $\delta$ is the Dirac distribution and $k_B$ is the Boltzmann constant.

Here, we note that the nonlinear prefactor $\hat {h}^3/(12\mu)$ arises from viscous resistance to the flow, describing the mobility of the fluid. In fact, (2.1) can be re-written as a gradient flow (Durán-Olivencia et al. 2019; Cates 2022; Sprittles et al. 2023) of the form

(2.2) \begin{equation} \frac {\partial \hat {h}(\hat {x},\hat {y},\hat {t})}{\partial \hat {t}} = \hat {\boldsymbol{\nabla }}\boldsymbol{\cdot }\left (M(\hat {h})\hat {\boldsymbol{\nabla }} \frac {\delta \hat {F}}{\delta \hat {h}} + \sqrt {2 k_B TM(\hat {h})}\boldsymbol{\hat {\eta }}\right)\!, \end{equation}

where $M(\hat {h})$ is the mobility and $\hat {F}[\hat {h}(\hat {x},\hat {y})]$ is an energy functional that gives rise to the pressure. Depending on the problem and assumptions, the mobility can take other forms (Glasner 2008). For example, in a crowded, narrow section of cytoplasm, one could use Darcy’s law to describe the flow through a porous medium, which reduces the mobility to ${\sim} \hat {h}$ (Gnann & Petrache 2018). A slip boundary, however, would enhance mobility by introducing an additional ${\sim} \hat {h}^2$ term (Zhang, Sprittles & Lockerby 2020). The ${\sim} M^{1/2}$ prefactor of the fluctuation term ensures that detailed balance is satisfied (Durán-Olivencia et al. 2019; Liu et al. 2024). The free energy $\hat {F}$ depends on what forces drive the flux in the channel, as described in the following sections.

2.2. Membrane dynamics

Deformation of the membrane, as shown in figure 1 $(a)$ , results in tension and bending forces, which lead to a change in the pressure $\hat p(\hat {x},\hat {y},\hat {t})$ of the fluid in the channel. As the membrane changes shape, so does the fluid flux in accordance with (2.1). The membrane is idealised as an isotropic linearly elastic solid with Young’s modulus $E$ , Poisson ratio $\nu$ and thickness $d$ (Landau & Lifshitz 1986). In the limit of small deflections, i.e. small spatial gradients in $\hat {h}(\hat {x},\hat {y},\hat {t})$ , which is natural from the scale separation, and small strains, the coupling between the pressure in the channel and the deflection of the membrane can be represented by the Föppl–von Kármán (FvK) equations (Audoly & Pomeau 2008; Howell, Kozyreff & Ockendon 2008). The FvK equations incorporate the effects of both bending and tension, the latter of which is highly coupled to the deformations in the membrane. We consider a membrane subject to constant tension $\gamma$ , which could arise from some external stretching of the membrane at its edges (or a pressure inside an adhered cell), which is much greater than the tension caused by local membrane deflections. Under this assumption, the relationship between membrane deflections and channel pressure becomes (Evans 1985; Peng & Lister 2019)

(2.3) \begin{equation} \hat {p}_{\textit{elastic}}(\hat {x},\hat {y},\hat {t}) =B \hat {{\nabla }} ^4\hat {h}(\hat {x},\hat {y},\hat {t}) - \gamma \hat {{\nabla }} ^2\hat {h}(\hat {x},\hat {y},\hat {t}), \end{equation}

where $B=Ed^3/12(1-\nu ^2)$ is membrane bending rigidity and $\gamma$ is the membrane tension (which is analogous to the surface tension at a fluid–fluid interface). Both of these components are generally present, but their relative strength will depend on the ratio of the characteristic horizontal length scale of the system, $L$ , to the bendocapillary length $l_{BC}= \sqrt {B/\gamma }$ (Deserno 2015). For $L/l_{BC}\ll 1$ , bending should be the prominent driving force, whereas for $L/l_{BC}\gg 1$ , we expect membrane tension to dominate. Note that the tension in a membrane can also be described by an integral constraint for the membrane length (Kodio, Griffiths & Vella 2017), or by solving the full FvK to give a more complete description of membrane deformation (Pihler-Puzović et al. 2013).

The pressure terms in (2.3) can be formulated as a gradient flow of (2.2) as the two contributions can be rewritten as functional derivatives of free energies. The bending energy is $\hat {F}_{{bend}}[\hat {h}(\hat {x},\hat {y})] = \int ({B}/{2}) |\hat \nabla^2 \hat h |^2\, {\rm d}\hat {A}$ and the surface energy is $\hat {F}_{\gamma }[\hat {h}(\hat {x},\hat {y})] = \int ( {\gamma }/{2}) \lvert \hat {\boldsymbol{\nabla }} \hat {h}\rvert ^2 \,{\rm d}\hat {A}$ .

2.3. Dynamics of adhesive molecules

When proteins or other adhesive molecules bind across the gap between the two membranes, additional forces are generated that contribute to the pressure $\hat p(\hat {x},\hat {y},\hat {t})$ (Carlson & Mahadevan 2015a , b ). Here, we use a Hookean elastic spring model to describe how bound proteins resist stretching and compression. Given a concentration $\hat {c}(\hat {x},\hat {y},\hat {t})$ of the bound adhesive molecules, an additional contribution to the pressure $\hat p_{\textit{adh}}$ takes the form

(2.4) \begin{equation} \hat p_{\textit{adh}}(\hat {x},\hat {y},\hat {t})=\kappa (\hat {h}(\hat {x},\hat {y},\hat {t})-l)\hat {c}(\hat {x},\hat {y},\hat {t}), \end{equation}

where $l$ is the equilibrium bond length and $\kappa$ is the spring coefficient.

To determine $\hat {c}(\hat {x},\hat {y},\hat {t})$ , we use a minimal description of the chemical kinetics of a single protein species, a full description of this model is given in Appendix A (Bell, Dembo & Bongrand 1984; Carlson & Mahadevan 2015a ). The adhesive molecules are assumed to be uniformly distributed on both surfaces at a constant concentration $c_0$ . The concentration of pairs of unbound proteins is thus $c_0-\hat {c}(\hat {x},\hat {y},\hat {t})$ . The molecules can bind or unbind at a rate of $(c_0-\hat c)K_{{on}}$ or $\hat cK_{\textit{off}}$ , respectively. The rate constant for binding, $K_{{on}}$ , is a Gaussian distribution about the equilibrium length of the bond, $l$ , with a standard deviation $\sigma$ ; whereas the rate constant for unbinding, $K_{\textit{off}}$ , is constant (Bell et al. 1984; Qi et al. 2001). Furthermore, we assume that the molecular dynamics of binding/unbinding are much faster than the time scale of viscosity-mediated membrane deformation, so that the adhesive molecule concentrations immediately adjust to the equilibrium values for a specific membrane height. Balancing the binding and unbinding rates then gives us the following expression for the concentration of a single species of bound molecules (Carlson & Mahadevan 2015b ):

(2.5) \begin{equation} \hat {c}(\hat {x},\hat {y},\hat {t})=c_0\frac {\exp \left (- ((\hat {h}(\hat {x},\hat {y},\hat {t})-l)/\sigma )^2\right)}{\exp \left (-((\hat {h}(\hat {x},\hat {y},\hat {t})-l)/\sigma)^2\right)+\tau _{{on}}/\tau _{\textit{off}}}, \end{equation}

where a value of $\sigma /l=0.2$ is used in the results of this study and the ratio of kinetic times for binding to unbinding is set to $\tau _{on}/\tau _{\textit{off}}=1/3$ Carlson & Mahadevan (2015a , b ). We note that with $\hat {c}(\hat {x},\hat {y},\hat {t})$ given by (2.5), the molecule kinetics can be reduced to a single function of $\hat {h}(\hat {x},\hat {y},\hat {t})$ , meaning that we do not need to solve for $\hat {c}(\hat {x},\hat {y},\hat {t})$ as a variable in our system, which would be the case if the time scale for molecules to bind or unbind were of the order of the time scale for viscous fluid transport in the channel. The pressure given by (2.4) and (2.5) can conveniently be written as a free energy $\hat {F}_{\textit{adh}}[h]$ , as is described in Appendix A.

2.4. Non-dimensional analysis

We work with non-dimensional versions of (2.1) in the remainder of this article. Different scalings of the lengths and time are used in accordance with the dominant physical effects for the work presented in the subsequent sections. In each case, (2.1) is left with one non-dimensional parameter, $Q$ , that is directly proportional to $\langle |\delta \hat {h}|\rangle$ , the average amplitude of thermal fluctuations in the film (Dhaliwal et al. 2024). Here, we outline the non-dimensionalisation used for each case; further details can be found in Appendix B.

In § 3, we study (2.1) on a one-dimensional (1-D) domain (which represents a 2-D film since $h$ is the height in the $z$ -direction), with $\gamma =0$ in (2.3) and no adhesive molecules/proteins. In this case, the horizontal length $\hat {x}$ is scaled by the domain size $L$ , as it is the only relevant horizontal length scale, the film height $\hat {h}(\hat {x},\hat {t})$ is scaled by the initial film height $\hat {h}_0$ and time $\hat {t}$ is scaled by $\tau _\mu = {12\mu L^6}/{B\hat {h}_0^3}$ , which is the time scale of viscous relaxation of an elastohydrodynamic thin film and can be obtained by scaling the left-hand side of (2.1) with the bending term on the right (Carlson 2018). With these scalings, the dimensionless version of (2.1) is

(2.6) \begin{equation} \frac {\partial h}{\partial t} = \frac {\partial }{\partial x} \left (h^3\frac {\partial }{\partial x} \left (\frac {\partial ^4}{\partial x^4} h\right)\right) + Q_{1\text{-}\textrm{D}}\frac {\partial }{\partial x} \left (h^{3/2}\eta \right)\!, \end{equation}

where $Q_{\textrm{1-D}}= ({L}/{\hat {h}_0})\sqrt { {2k_B T L}/{Bw}}$ . Here, $w$ is the width of the quasi-one-dimensional (quasi-1-D) film in the y-direction.

In § 4, we study (2.1) on a 2-D domain (i.e. a three-dimensional (3-D) film), with either  $\gamma$ or $B$ set to zero in (2.3) and the protein pressure contribution described in § 2.3. In this case, a physical horizontal length scale can be obtained by balancing the relevant membrane pressure with the protein spring pressure (Carlson & Mahadevan 2015a ). When bending dominates, this gives $L_B= (B/c_0\kappa)^{1/4}$ and while tension dominates, $L_\gamma =(\gamma /c_0\kappa)^{1/2}$ . We scale $\hat {x}$ and $\hat {y}$ by this length scale, $\hat {h}$ by the protein equilibrium length $l$ , the protein concentration $\hat {c}$ by $c_0$ , and $\hat {t}$ by a time scale constructed as for the 1-D case, but with the gradients scaled by the new horizontal length scale $L_\gamma$ or $L_B$ . For the bending-driven case, inserting this scaling, i.e. $\hat {x}=L_Bx$ , $\hat {y}=L_By$ and $\hat {t}=t( {12\mu B^{1/2}}/{l^3(c_0\kappa)^{3/2}})$ , gives the following dimensionless equation:

(2.7) \begin{equation} \text{Bending} \,(\gamma =0): \quad \frac {\partial h}{\partial t} = \boldsymbol{\nabla }\boldsymbol{\cdot }\left (h^3\boldsymbol{\nabla }\left ({\nabla} ^4h+ (h-1)c\right) \right) + Q_{B}\boldsymbol{\nabla }\boldsymbol{\cdot }\left (h^{3/2}\boldsymbol{\eta }\right)\!, \end{equation}

where $Q_{B}= ({1}/{lL_B})\sqrt { {2k_B T }/{c_0\kappa }}$ . For tension-driven films, the equivalent scalings, i.e. $\hat {x}=L_\gamma x$ , $\hat {y}=L_\gamma y$ and $\hat {t}=t({12\mu \gamma }/{l^3(c_0\kappa)^{2}})$ , give us the appropriate form of the non-dimensional thin film equation:

(2.8) \begin{equation} \text{Tension}\, (B=0): \quad \frac {\partial h}{\partial t} = \boldsymbol{\nabla }\boldsymbol{\cdot }\left (h^3\boldsymbol{\nabla }\left (-{\nabla} ^2h+ (h-1)c\right) \right) + Q_{\gamma }\boldsymbol{\nabla }\boldsymbol{\cdot }\left (h^{3/2}\boldsymbol{\eta }\right)\!, \end{equation}

where $Q_{\gamma }=({1}/{lL_{\gamma }})\sqrt { {2k_B T }/{c_0\kappa }}$ .

2.5. Finite element solver

In this article, (2.6), (2.7) and (2.8) are solved using the finite element method on a rectangular domain. In all simulations, periodic boundary conditions are imposed in the horizontal directions at the four boundaries. Since the film height $h$ is a dependent variable, simulations are performed in both 1-D and 2-D domains, which represent 2-D and 3-D films, respectively. The source code for the simulations in this article is available on GitHub (Dhaliwal 2025).

To solve (2.1), we set up a separate equation for the pressure $p$ , with contributions coming from (2.3) and (2.4). Equation (2.1) is divided into a system of two (plus an additional one for the film curvature ${\nabla} ^2 h(x,y,t)$ when there is a bending pressure) coupled second-order partial differential equations. The system is reformulated into weak form where boundary terms can be neglected due to the periodic boundary conditions. The scalar fields ${\nabla} ^2 h(x,y,t)$ , $p(x,y,t)$ and $h(x,y,t)$ are then discretised with linear elements, and the system of equations is solved using a Newton’s method solver from the FEniCS finite element library (Logg et al. 2012). Time integration is performed using an implicit first-order finite difference scheme. The domain and its discretisation in space and time are chosen according to the relevant non-dimensional form presented in § 2.4. For the 1-D simulations presented in § 3.2, a domain of length $L=1$ is used with $\Delta x=0.01$ and $\Delta t =1\times 10^{-8}$ . For the 2-D simulations presented in § 4.2, a square $100\times 100$ -cell grid is used with both $\Delta x$ and $\Delta t$ varying in the range $1{-}3$ (which means that $L$ is in the range $100{-}300$ ) for the various results presented.

The random vector $\boldsymbol{\eta }(x,y,t)$ is implemented by choosing random numbers using the ‘normal’ function in the ‘random’ class of NUMPY (Harris et al. 2020). For the numerical solution at each time step, the two components of $\boldsymbol{\eta }(x,y,t)$ are each assigned a new value at every point in the mesh. Values are drawn from a Gaussian distribution with zero mean and a variance of $1/(\Delta x \Delta t)$ in 1-D and $1/(\Delta x^2 \Delta t)$ in 2-D. Due to the stochastic nature of the problem, each individual run is not to be considered as deterministic. For each set of parameters, we run multiple independent realisations and report the ensemble averages of the extracted/predicted parameters.

3. Waiting time for the initiation of adhesion

Membrane adhesion is facilitated by the binding of membrane-anchored adhesive molecules, which requires them to be in close range. It is then natural to ask: how long does it take for fluctuations to deform the membranes sufficiently so that they can bind to initiate adhesion? In this section, we use the model described in § 2 to study the average waiting time for initial contact to occur. We limit the study to a 1-D periodic domain, to avoid excessive computational costs, but expect that our results should be qualitatively similar for a 2-D domain. The physical picture of our simplified model is shown in figure 1 $(a)$ ; the two membranes are separated by a channel of viscous fluid with initial uniform thickness $h_0$ . Although the channel height is initially only fluctuating about $h_0$ , some part of the domain eventually reaches a critical thickness $h^*$ at which the two adhesive molecules come in contact.

3.1. Rare-event theory

It is well known that thermal fluctuations can generate waves across the membrane, where the average wave amplitudes $\delta h_q$ of frequency $q$ , or ‘roughness’ of the membrane, depends on bending moduli and tension coefficient (Helfrich & Servuss 1984; Rautu et al. 2017). If the $h_0-h^*$ is of the same order as $\delta h_q$ , then thermal fluctuations can easily initiate binding. If $h_0-h^*$ is sufficiently larger than $\delta h_q$ , attachment of the molecules then requires the membrane profile $h(x,t)$ to attain a rather unlikely shape, which may take a long time. The process can thus be thought of as $h(x,y,t)$ fluctuating in an energy landscape that does not favour large deviations from $h_0$ , until it eventually gets ‘lucky’ and reaches $h^*$ at some point in the domain. Although the waiting time for protein binding is a random variable, the rare event theory can provide a prediction for the ensemble-averaged waiting time, known as the Eyring–Kramers law, which states that $\langle t_B \rangle \sim C \exp ({2(F[h_B(x)]-F[h_0])/Q_{\text{1-D}}^2})$ , where $\langle t_B \rangle$ is the ensemble-averaged binding time, $C$ is a prefactor, $F[h(x,t)]$ is the energy of a profile $h(x,t)$ and $h_B(x)$ is the final binding profile (Eyring 1935; Kramers 1940; Liu et al. 2024).

If the membrane primarily exhibits tension rather than bending, i.e. $L\gt l_{BC}$ , the energy is determined by the second term in (2.3), and the problem is analogous to the rupture of a thin viscous liquid film with a free surface due to van der Waals forces. This problem was recently studied by Sprittles et al. (2023), where the rare-event description based on the Eyring–Kramers formula was validated by both numerical simulations of the stochastic thin film equations and in molecular dynamics simulations. In this paper, these methods are adapted to determining the binding time between two membranes when rather the bending dominates the pressure term in (2.3).

When bending dominates over tension, i.e. $L\lt l_{BC}$ , the non-dimensional energy functional for the 1-D membrane becomes

(3.1) \begin{equation} F[h(x)] = \int ^1_0 \frac {1}{2} \left (\frac {\partial ^2 h}{\partial x ^2}\right)^2 \,{\rm d}x. \end{equation}

The pressure term in (2.3) can be recovered by taking a functional derivative of $F$ with respect to $h(x)$ . Finding the average waiting time for adhesion requires finding the profile $h_B(x)$ that minimises $F[h]$ while also maintaining conservation of mass, which enters as the constraint $\int _0^1(h-h_0)\,{\rm d}x = 0$ , and satisfying the periodic boundary condition. This constrained optimisation problem can be solved using the Euler–Lagrange equation, which gives the fourth-order polynomial profile shown by the dashed blue line in figure 2 $(a)$ . The energy corresponding to this profile is $F_{B}= 360 (h_0-h^*)^2$ . With a sharp asymptotics analysis (see Appendix C for details), the predicted average attachment time is given by

(3.2) \begin{equation} \langle t_B \rangle = \frac {1}{\beta (h_B)}\sqrt {\frac {Q_{\text{1-D}}^2}{(2\pi)^7}} \exp \left (720\left (\frac {h_0-h^*}{Q_{\textrm{1-D}}}\right)^2\right)\!, \end{equation}

where the parameter $\beta$ is given by

(3.3) \begin{equation} \beta (h_B) = (h_0-h^*)(2\pi)^6\left (\frac {1}{2}h_0^3 +\frac {3}{8}h_0(h_0-h^*)^2\right)\!. \end{equation}

We note that the expressions in (3.2) and (3.3) predict the attachment time across the channel without any fitting parameters. They are valid in the limit of small noise, i.e. either large $h_0-h^*$ or small $Q_{\text{1-D}}$ , such that the trajectory of initial attachment must cross $h_B(x)$ in the energy landscape.

Figure 2. $(a)$ Film profiles at time of attachment obtained from 15 independent solutions (centred around the point of ‘contact’) with parameters $Q_{\text{1-D}}=5$ and $h^*=0.3$ . The dotted black line represents the average of the individual simulations and the dashed blue line represents the theoretical prediction from the Euler–Lagrange equation. $(b)$ Average waiting time for adhesion $\langle t^* \rangle$ as a function of $(h_0-h^*)^2$ for different values of the noise amplitude $Q_{\text{1-D}}$ . The lines represent the predicted value of $\langle t_B \rangle$ from (3.2) with no free parameter. Error bars represent the standard deviation for a set of $N=15$ simulations for each data point. The shaded blue colour is intended as a guide to the eye to emphasise the region where rare-event theory is valid, i.e. attachment events are sufficiently unlikely.

3.2. Numerical predictions

Numerical simulations of (2.6) are used to test the prediction in (3.2) considering only the effect of membrane bending. The simulations are initiated with a flat membrane of height $h_0=1$ . Since there are no adhesive molecules, i.e. no force pulling/pushing on the membrane, the minimum film height fluctuates randomly with time, occasionally reaching a lower value than it has before. The simulations are run until the lowest point in the film reaches the cutoff height $h^*$ , indicating that the molecules are now within binding range for initial attachment, and we denote this time as $t^*$ . The point in the periodic domain at which attachment occurs is random. If the deformation required for film rupture, $h_0-h^*$ , is small, attachment occurs rapidly and the profile shape at $t^*$ varies significantly from simulation to simulation. When $h_0-h^*$ is larger, however, the minimum film height fluctuates for a long time before reaching attachment (the number of time steps being many orders of magnitude), often coming close many times before finally reaching $h^*$ . When this is the case, the individual profiles at $t^*$ are consistent, as shown in figure 2 $(a)$ , which suggests that attachment does indeed tend to occur at a specific minimum energy profile.

When the value of $h_0-h^*$ is large, as is the case for the profiles shown in figure 2 $(a)$ , one would expect that the final profile is similar to the fourth-order polynomial predicted theoretically in § 3.1. The dashed blue line in figure 2 $(a)$ shows that this is indeed the case if the profiles are centred about their minimum and then averaged. This indicates that it is truly rare for fluctuations to make the film reach $h^*$ and that this almost always occurs very close to the minimum energy profile that allows attachment.

The average attachment time, $\langle t^* \rangle$ , is plotted for various values of $Q_{\text{1-D}}$ and $h^*$ in figure 2 $(b)$ . It is clear that the time increases rapidly when $h^*$ is reduced and eventually grows exponentially with $(h_0-h^*)^2$ for constant $Q_{\text{1-D}}$ . The theoretical prediction for the attachment time across a membrane channel from (3.2) provides an excellent quantitative prediction for the average rupture time when $h_0-h^*$ is large. For smaller values of $h_0-h^*$ , the rare-event theory is not valid, as the membrane profile does not necessarily cross the saddle point on its way to attachment. We note that the performance of the theoretical prediction also depends on $Q_{\text{1-D}}$ , as fluctuations are more likely to cause large deformations to the film profile, making the event less rare as well as introducing larger errors in the Laplace method used to approximate the mean first passage time (see Appendix C). Generally, the rare-event prediction is valid when the attachment event is sufficiently rare, meaning that the average attachment time is sufficiently large. This can be achieved either by having a small value of $Q_{\text{1-D}}$ or a large height difference $h_0-h^*$ , as in the blue shaded region of figure 2 $(b)$ .

4. Coarsening of adhesion patches

We now turn our attention to the next stage in the dynamics, when the adhesive molecules have formed bonds ( $t \gt \langle t_B\rangle$ ) that lead to a rich coarsening dynamics. To do this, we study (2.2) with the protein binding model of (2.4) and (2.5) on a 2-D domain, where we vary the membrane pressure term and the mobility factor.

Once adhesive molecules start to bind across the thin film, the membrane is pulled towards the substrate at the binding site, leading to further binding of adjacent proteins (Bell et al. 1984). This process brings the surfaces closer to each other, squeezing liquid out laterally through the channel. Conservation of the liquid, however, may lead to distinct regions where the membrane separation is small and proteins are bound or where membrane separation is large and proteins unbound, as is depicted in figure 1 $(c)$ (Fenz et al. 2017; Dinet et al. 2023). The network of unbound domains containing excess liquid may then coarsen to reduce the deformation of the membrane, i.e. reducing the overall energy of the system. Smaller pockets/lumen-like structures shrink as liquid is transported towards larger pockets/lumens, which subsequently grow in size. Such coarsening is reminiscent of when liquid droplets form in a dewetting liquid thin film (Otto et al. 2006). It also falls into the broader category of coarsening dynamics in physical systems, which has been extensively studied both theoretically and experimentally (Hohenberg & Halperin 1977; Furukawa 1985; Bray 1994; Camley & Brown 2011; Cates 2022). We will now describe how our model fits into this context, demonstrating that elastohydrodynamic thin films display distinct coarsening behaviour due to the combined effects of elasticity and viscosity.

4.1. Domain coarsening of adhered patches

Equation (2.1) can be seen as the governing equation of a dynamical phase field, in which the film height $h(x,y,t)$ is the sole order parameter (McLeish et al. 2003; Cates 2022). When coupled with the pressure term given by (2.3), it is in fact quite similar to the widely studied dynamical model B, which is used to describe the coarsening of a conserved order parameter (Hohenberg & Halperin 1977). Nevertheless, a number of features distinguish our system. First, a nonlinear ${\sim} h^3$ mobility (as described in § 2.1) arises from the lubrication flow in the narrow channel. Second, rather than the typical double-well potential, our system imposes an energy landscape arising from the adhesive molecules/proteins described in § 2.3 and Appendix A. Finally, and perhaps most significantly, the classical Cahn–Hilliard Laplacian free energy term (which, in this case, represents isotropic membrane tension) is supplemented by the fourth-order bending term as described in (2.3).

A quantitative understanding of phase separating systems is achieved through the dynamical scaling hypothesis, which states that in the later stage of a coarsening process, the domain structure of a system (as quantified by the structure function $S(\boldsymbol {k},t)=\langle h_{\boldsymbol {k}}(t)h_{\boldsymbol {-k}}(t)\rangle$ ) is constant in time if one rescales the lengths by a single characteristic length scale $L_c(t)$ (Bray 1994; Camley & Brown 2011). This concept has been shown to be valid both numerically and experimentally in many systems, with power-law growth displayed when $L_c(t)$ is computed as a moment of $S(\boldsymbol {k},t)$ or the radius of individual particles (Furukawa 1985; Shinozaki & Oono 1993; Sung et al. 1996; Tateno & Tanaka 2021; Saiseau et al. 2024; Su et al. 2024). For coarsening of the aforementioned model B system, $L_c(t)$ is known to grow according to a ${\sim} t^{1/3}$ power law, as has been validated by a number of numerical studies (Vladimirova, Malagoli & Mauri 1998; Tiribocchi et al. 2015), and can also be predicted theoretically using scaling arguments, renormalisation group theory or Lifshitz–Slyozov–Wagner theory (in the limit where the minority phase occupies a small volume fraction) (Lifshitz & Slyozov 1961; Wagner 1961; Bray 1994; Camley & Brown 2011).

For the case where the bending energy drives coarsening rather than surface tension, it is unclear if the coarsening rate will follow the ${\sim} t^{1/3}$ power law. This essentially corresponds to replacing the surface tension contribution to the classical free energy functional, $F_{\gamma }(h(x,y)) = \int ({1}/{2}) \lvert \boldsymbol{\nabla }h\rvert ^2 \,{\rm d}A$ , by the corresponding bending term $F_{\textit{bend}}(h(x,y)) = \int ({1}/{2}) ({\nabla} ^2 h)^2 \,{\rm d}A$ . We now follow the derivation of Bray (1994) to make a scaling prediction for how bending changes the coarsening dynamics (details are provided in Appendix D.2). We ignore the effects of the nonlinear mobility as well as the specifics of the potential function. In analogy with the coarsening literature, we note that the pressure in our system takes a form that can be interpreted as a chemical potential $\varPhi$ (Hohenberg & Halperin 1977; Bray 1994; Cates 2022). We thus adopt the notation $ p\equiv \varPhi =\delta F/\delta h$ , where $F$ consists of the aforementioned bending term and a symmetric double-well potential function $V(h)$ . The chemical potential is thus

(4.1) \begin{equation} \varPhi = \frac {\partial V}{\partial h} +{\nabla} ^4h \end{equation}

and the film height then evolves according to a continuity equation with flux $\boldsymbol{j}=-\boldsymbol{\nabla }\varPhi$ .

In late-stage coarsening when motion of the interface between two domains is slow, diffusion of the order parameter $h$ and chemical potential $\varPhi$ in the bulk is fast, so both should satisfy Laplace’s equation ${\nabla} ^2h={\nabla} ^2\varPhi =0$ in the bulk regions. The flux through the interface, and thus by continuity the velocity of the interface, can then be determined by finding $\varPhi$ at the interface.

At an interface between the two phases with radius of curvature $R$ , $\varPhi$ can be re-expressed in terms of the coordinate $g$ representing distance along the unit vector $\boldsymbol{\hat {g}}$ in the direction perpendicular to the interface ( $g=\pm \infty$ in the bulk and $g=0$ at the centre of the interface). Noting that $\boldsymbol{\nabla }h = (\partial h/\partial g)\boldsymbol{\hat {g}}$ and $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\hat {g}}=1/R$ near the interface, we find that

(4.2) \begin{equation} \varPhi = \frac {\partial V}{\partial h} +\frac {\partial ^4 h}{\partial g^4}+\frac {2}{R}\frac {\partial ^3 h}{\partial g^3}-\frac {1}{R^2}\frac {\partial ^2 h}{\partial g^2}+\frac {1}{R^3}\frac {\partial h}{\partial g}. \end{equation}

The value of $\varPhi$ at the interface is then determined by multiplying (4.2) by $\partial h/\partial g$ and integrating across the interface from $g=-\infty$ to $g=\infty$ . By imposing the boundary condition that $h$ is constant in the bulk and assuming that the chemical potential is a double well with equal potential on both sides, integration gives

(4.3) \begin{equation} \varPhi \Delta h = -\frac {2}{R}\int ^\infty _{-\infty }\left (\frac {\partial ^2 h}{\partial g^2}\right)^2\,{\rm d}g+\frac {1}{R^3}\int ^\infty _{-\infty }\left (\frac {\partial h}{\partial g}\right)^2\,{\rm d}g, \end{equation}

where $\Delta h$ is the height difference between the domains on the inside and outside of the interface, and the integrals $-2\int ^\infty _{-\infty} ({\partial ^2 h}/{\partial g^2})^2\,{\rm d}g$ and $\int ^\infty _{-\infty} ({\partial h}/{\partial g})^2\,{\rm d}g$ can effectively be interpreted as 2-D line tension-like coefficients $\varGamma _1$ and $\varGamma _2$ that represent an energy per unit length associated with the existence of the interface between the two phases in 2-D space. Equation (4.3) thus represents a Gibbs–Thomson-like boundary condition that determines the value of $\varPhi$ on interfaces with radius of curvature $R$ . Since the flux of fluid is proportional to the gradient of the chemical potential (assuming a constant mobility), the velocity of the interface scales as $-\partial \varPhi /\partial g$ . Assuming that the domains are circular with radius $R$ corresponding to the macroscopic length scale $L_c$ , this then gives the following growth law for bending-driven coarsening:

(4.4) \begin{equation} \frac {{\rm d}L_c}{{\rm d}t}\sim \frac {\varGamma _1}{L_c^2}+\frac {\varGamma _2}{L_c^4}. \end{equation}

Equation (4.4) stands in contrast to the model B case which has only a single contribution to the line tension, giving rise to an ${\sim} L_c^{-2}$ growth rate. This suggests that bending-driven coarsening should have a coarsening rate with an exponent somewhere between $1/3$ and $1/5$ with time, with the value being determined by the relative sizes of $\varGamma _1$ and $\varGamma _2$ , which in turn depend on the shape of the potential function. We note that this simple scaling prediction only describes the difference between surface tension-driven and bending-driven coarsening. It does not take into account the other differences between our model and the classical model B system described by Bray (1994) such as the ${\sim} h^3$ mobility and the protein binding potential energy. Since the potential energy function $V(h)$ is not specified, we expect these results to be valid for other adhesive potentials such as van der Waals interactions. Equation (4.4) suggests a fundamental distinction between surface tension-driven and bending-driven coarsening, which will be investigated numerically for our specific system in § 4.2.

4.2. Numerical investigation of bending-driven coarsening dynamics

We now present the results of numerical simulations of the phase separation of adhesive molecules/proteins during membrane adhesion when dominated by bending of the membrane, as described in (2.7). The membrane is always initialised as a flat profile with non-dimensional initial height $h(x,y,0)= \hat {h}_0/l\gt 1$ and the non-dimensional width of the protein binding rate distribution, $\sigma _{{on}}/l=0.2$ . Figure 3( $a$ ) shows the contour plots (see supplementary movie 1 available at https://doi.org/10.1017/jfm.2025.10872) of the film profile at three different times when the initial film height is just at the edge of the range in which the proteins can bind. Fluctuations are clearly present for the first time steps, where the film is pulled to the substrate by the adhesive molecules. Most of the film profile moves towards the equilibrium length of the adhesive molecules ( $h=1$ in non-dimensional terms) to reduce the energy of the system. Due to conservation of liquid mass, excess liquid must flow somewhere else, but periodic boundary conditions prevent fluid from leaving the domain. Thus, in some regions, the film is pushed up, forming circular ‘pockets’ or ‘lumens’ in which proteins are unbound and the film height is significantly larger than the initial height, as can be seen in the first panel to the left in figure 3( $a$ ).

Figure 3. Contour plots illustrating the height $h(x,y,t)$ of thin films in a 2-D domain under the influence of protein binding and unbinding for different times $t$ , with $h_0=1.4$ . In panel $(a)$ , the fluctuation parameter is $Q_{B}=0.005$ , whereas in panel $(b)$ , it is set to $Q_{B}=0.5$ . Although similar coarsening dynamics occur at late times regardless of $Q_{B}$ , the size of the domains at $t=350$ is somewhat larger for $Q_{B}=0.5$ due to early-time coalescence.

With time, the detached domains coarsen in an Ostwald ripening-like process, resulting in fewer, larger domains as can be seen in the second and third panels of figure 3( $a$ ). Once the circular pockets form, their centres remain essentially fixed as fluid is slowly transported from smaller to larger domains in a diffusion-like manner. When the size of a single pocket drops below a critical threshold, it suddenly ‘implodes’, which causes a minor horizontal rearrangement of nearby domains (see supplementary movies). Figure 4( $a$ ) shows the time evolution of the characteristic size $L_c$ for these domains, as computed using the method of Shinozaki & Oono (1993), averaging the structure factor $S(\boldsymbol {k},t)$ across 15 independent realisations (details of how we calculte $L_c$ are provided in Appendix D.1). Although the distinct ‘implosion’ events lead to discrete jumps in the growth for an individual simulation, the ensemble-averaged $L_c$ grows smoothly with time following a power law. The growth is significantly slower than the ${\sim} t^{1/3}$ power law observed for a model B system (Bray 1994; Tiribocchi et al. 2015), instead, the exponent is closer to $1/5$ , which can be expected based on our analysis in § 4.1.

Figure 4. Length scale $L_c$ computed using (D5) from the average of $N=15$ individual simulations under the same conditions as the data in figure 3. In panel $(a)$ , the fluctuation parameter is $Q_{B}=0.005$ , whereas in panel $(b)$ , it is set to $Q_{B}=0.5$ . Power law coarsening with an exponent well below $1/3$ is observed in both cases, but starts off with a larger domain size when $Q_{B}=0.5$ .

The contour plots in figure 3( $b$ ) show the evolution of height profiles of a coarsening film when the amplitude of the thermal fluctuations, represented by the parameter $Q_{B}$ , is increased by two orders of magnitude. With higher noise amplitude, the lumens are irregular in shape and constantly deforming (see supplementary movie 2). The fluctuations also cause significant lateral motion of the domains, which leads to some instances of coalescence at very early times. Despite this, at the later time, large droplets seem to repel each other and domain growth is primarily caused by fluid flow through adhered patches, as was observed in the low noise amplitude case. Interestingly, the power law for domain growth is relatively unaffected by the strength of the fluctuations, as shown in figure 4( $b$ ). The size of the domains before the late-stage power law growth is somewhat higher when the fluctuations are larger, perhaps due to coalescence events at early times, as can be seen in the first panels of figure 3(a,b). The primary effect of fluctuations is thus to provide the perturbations that trigger the instabilities giving rise to the initial configuration, rather than driving coarsening.

To better understand how and why the coarsening rate of our system differs from the more familiar model B system, we perform additional simulations in which we vary the initial membrane height $h_0$ , the fluctuation strength $Q$ , the mobility prefactor in the thin film equation and which term in (2.3) we use. In general, we find that coarsening behaviour is a preserved feature of the system even when these parameters are changed. To ensure an initial configuration with many small pockets in the domain, $h_0$ should be close to the edge of the binding range of the adhesive molecules. If $h_0$ is too small, the entire domain is already in a bound state and few pockets form. If, however, $h_0$ is too large, it takes a very long time for contact to occur, which is typically at only one point in the domain, leading to a small number of initial pockets when coarsening starts.

To begin with, we simulate coarsening using the tension term of (2.3) (corresponding to the well-studied Cahn–Hilliard free energy) instead of the bending term, as in (2.8). As expected, coarsening behaviour is observed, as shown in figure 5. Comparing the first snapshots of figures 5 $(a)$ and 5 $(b)$ demonstrates that larger initial heights lead to fewer domains at early times. At later times, we observe that the coarsening actually appears to be slower when $h_0$ is increased.

Figure 5. Contour plots illustrating the height $h(x,y,t)$ of thin films in a 2-D domain under the influence of protein binding and unbinding for different times $t$ , with $Q_{\gamma }=0.01$ . In panel $(a)$ , the initial height is $h_0=1.25$ , whereas in panel $(b)$ , it is set to $h_0=1.45$ . Increasing $h_0$ leads to a larger initial domain size, but slower coarsening at late times.

Figure 6 $(a)$ shows how the characteristic length scale $L_c$ of a viscous film profile grows with time when the membrane exhibits only interfacial tension. When $h_0$ is small, the growth exponent is close to ${\sim} t^{1/3}$ , in accordance with what one expects for the well-known model B system. For higher values of $h_0$ , the growth rate for late-stage coarsening decreases. This seems to be caused by the nonlinear ${\sim} h^3$ term arising from the viscous resistance to flow in the governing equation. To confirm this, we perform additional simulations in which we replace the nonlinear mobility by a constant mobility with value  $1$ . The inset in figure 6 $(a)$ shows that $L_c\sim t^{1/3}$ growth is observed regardless of $h_0$ when the mobility is constant, as expected.

Figure 6. $(a)$ Scaling length $L_c$ as a function of time for tension-driven coarsening in a viscous film with varying $h_0$ when $Q_{\gamma }=0.01$ . $L_c$ is calculated from the average $S(\boldsymbol {k},t)$ for $N=22$ independent simulations. The inset shows the results when a constant mobility is used instead. $(b)$ Scaling length $L_c$ as a function of time for bending-driven coarsening in a viscous film with varying $h_0$ when $Q_B=0.01$ . $L_c$ is calculated from the average $S(\boldsymbol {k},t)$ for $N=20$ independent simulations. The inset shows the results when a constant mobility is used instead.

In figure 6 $(b)$ , similar results are presented, but this time, the tension term in (2.3) is replaced by the bending term, i.e. we set $\gamma \gt 0$ and $B=0$ , and solve 2.8. As suggested by the theory in § 4.1, the power law for bending-driven coarsening is reduced from $L_c\sim t^{1/3}$ . The inset in figure 6 $(b)$ shows that this is indeed the case when the mobility is kept constant, with a growth exponent of approximately $1/4$ consistently observed. For the viscous case, there is a decrease in the growth exponent when the value of $h_0$ is increased, but this effect is not as pronounced as it is for the tension-driven films.

The decreased coarsening rate as $h_0$ is increased may seem counterintuitive at first glance, as a thicker film should lead to less viscous resistance to fluid flow and thus a higher mobility. In the coarsening process, however, the fluid flow between pockets must pass through the adhered region, where the film height is fixed at the equilibrium height of the adhesive molecules ( $h=1$ ), thus restricting the mobility of the flow through this region. Increasing $h_0$ thus only increases the depth of the pockets, which decreases the flux through the interface between the pocket and the attached region.

Figure 7. $(a)$ Scaling length $L_c$ as a function of time for tension-driven coarsening in a viscous film with varying $Q_{\gamma }$ . These results are for $h_0=1.25$ and the $L_c$ is calculated from the average $S(\boldsymbol {k},t)$ for $N=22$ independent simulations. $(b)$ Scaling length $L_c$ as a function of time for bending-driven coarsening in a viscous film with varying $Q_B$ . These results are for $h_0=1.4$ and the $L_c$ is calculated from the average $S(\boldsymbol {k},t)$ for $N=15$ independent simulations.

Figure 7 shows the coarsening behaviour for both tension-driven and bending-driven viscous films as $Q$ is varied by two orders of magnitude. In both cases, the value of $L_c$ at the beginning of coarsening is larger when the fluctuations are increased to $Q=0.5$ , due to increased coalescence at very early times. In the bending case, large fluctuations also lead to an earlier onset of the power-law coarsening regime, as can be seen in figure 7 $(b)$ . During the late coarsening stage, however, fluctuations do not play a significant role. In figure 7 $(a)$ , we observe that increasing the magnitude of fluctuations leads to a slight decrease in the power law for tension-dominated films, while figure 7 $(b)$ shows that changing $Q_B$ does not seem to significantly affect the power law of bending-dominated films. This seems to suggest that the main role of the fluctuations is simply to provide perturbations at early times which initiate the adhesion process.

5. Discussion

We have demonstrated that the waiting time for membranes to come into contact across a viscous channel can be effectively predicted by the rare-event theory. Although computational concerns limited the present study to 1-D films, the methods we present can be extended to two dimensions. In that case, the rare-event prediction will have the same form, but a modified prefactor. Our results demonstrate the utility of rare-event theory in predicting the average time for important but unlikely events to occur in stochastic systems. Such statistical techniques are relevant in many biological systems where rare events caused by stochastic fluctuations may be necessary for important biological processes (Chou & D’Orsogna 2014). Although the adhesion of live cells often involves active cellular machinery such as actomyosin and filopodia to help cells bind, our results may be useful for estimating the likelihood that cells may come into contact when such active mechanisms are absent and may help to distinguish between actively and passively driven dynamics.

The membrane coarsening results of § 4 present a minimal mathematical description of coarsening behaviour similar to the recent experiments of Dinet et al. (2023). In their experiments, an osmotic shock is applied to GUVs initially attached through biotin–neutravidin bonds to a supported lipid bilayer, pushing liquid from inside the GUV into the membrane–membrane channel, forming small unattached lumens. As the system evolves, the smaller lumens shrink and eventually disappear, while larger lumens grow. Our work suggests a physical model to explain the intrinsic coarsening process in those experiments. Based on the results of § 4, one would expect to observe power-law coarsening of the domains in such experiments. Unfortunately, the data in the experiments of Dinet et al. (2023) are not quite amenable to study the coarsening law because the separate processes of fluid flux through the membrane into the channel, flux of fluid from pocket to pocket and flux of fluid from the pockets to the exterior region are all occurring simultaneously. Nevertheless, it points to experimental realisations allowing for the observation of the predicted power-law coarsening dynamics by isolating the interluminal coarsening process, e.g. by using a larger GUV so the centre of the adhesion patch is not as affected by flux to the outside.

Our work provides a framework for understanding coarsening exponents in thin films. The results of § 4.2 provide a scaling-based prediction that the power-law exponent is smaller for bending-driven than for tension-driven coarsening. A more rigorous prediction of the growth rate for bending-driven coarsening might be feasible applying renormalisation group theory, as is described by Bray (1994). Whether a physical system falls under the bending- or tension-dominated regime depends on the ratio of the characteristic length scale of the domains to the bendocapillary length scale, $L_{bc} = \sqrt {B/\gamma }$ . As such, there may be a cross-over between the two power-law regimes as the domain size grows past $L_{bc}$ during coarsening. Observing this, however, would require simulations of growth across a larger range of length scales than that shown in this study. An estimate of $L_{bc}$ based on typical values for biological membranes is smaller than the size of lumens in the experiments of Dinet et al. (2023), which suggests that tension should dominate for that system. Even so, we would expect a growth exponent with time that is smaller than ${\sim} t^{1/3}$ due to the nonlinear mobility coming from the fluid viscosity. It would also be interesting to see if the coarsening rate would be altered if one were to change the membrane properties to have a larger value of $L_{bc}$ .

The membrane deformation containing just the bending and isotropic tension terms (see § 2.2) is a simplification that ignores the fact that the tension will be affected by deformation (Landau & Lifshitz 1986; Peng & Lister 2019). This allows us to characterise some aspects of how bending and isotropic tension act differently, but prevents us from getting a full understanding of how the deformation of a thin elastic sheet affects the dynamics of adhesion and coarsening. To complete the picture, the full Föppl–von Kármán equations could be implemented to properly describe stretching in the membrane (Lister, Peng & Neufeld 2013; Pihler-Puzović et al. 2013).

Another observation to highlight is that the noise amplitude seems to have no effect on the power law for the coarsening process in both tension- and bending-dominated regimes. This is not true, however, for the spreading of droplets on a flat substrate, which can be seen as the coarsening of a single lumen. In the tension-dominated regime, as the droplet spreads, the growth of its radius follows the classical Tanner’s law (Tanner 1979) if the system is deterministic and follows a fluctuation-enhanced Tanner’s law (Davidovitch et al. 2005) in the presence of noise. In the bending-dominated regime, a different power law for growth was also found for deterministic (Lister et al. 2013) and stochastic (Carlson 2018; Pedersen et al. 2019) systems. One key difference between our set-up and droplet spreading is the presence of bound adhesive molecules, which keeps the membrane fixed at a constant height in the adhered region. Another difference is that instead of a single droplet, we have a network of pockets connected by adhesion patches.

Our numerical results predict that coarsening is slower for larger initial heights in viscous, tension-driven coarsening, which is directly relevant to coarsening in dewetting liquid films. Experiments on dewetting nanometric polymer films by Limary & Green (2003) indeed showed that the coarsening exponent was smaller when the initial film height was increased. Although the driving potential in such a system comes from the surface energies of the materials rather than protein-like binding, the governing equations should otherwise be the same and we would expect similar dynamics. Further work on this problem could lead to a more complete understanding of how and why an increased initial film height delays coalescence.

Finally, a future interesting avenue to explore would be looking into experimental systems of elastohydrodynamic coarsening beyond a biological context. The predicted coarsening exponent (§ 4.1) is not specific to the molecular binding-based adhesive potential in our model, so similar coarsening dynamics should be observed even if the adhesive force in (2.4) is replaced by some other physical mechanism. For example, if liquid is trapped between two thin elastic sheets grafted with polymer brushes that are attracted to each other (Mani, Gopinath & Mahadevan 2012), one might observe coarsening of non-adhered pockets quite similar to the observations in our simulations.