Peeling dynamics of fluid membranes bridged by molecular bonds: moving or breaking

Biological adhesion is a critical mechanical function of complex organisms. At the scale of cell–cell contacts, adhesion is remarkably tunable to enable both cohesion and malleability during development, homeostasis and disease. It is physically supported by transient and laterally mobile molecular bonds embedded in fluid membranes. Thus, unlike specific adhesion at solid–solid or solid–fluid interfaces, peeling at fluid–fluid interfaces can proceed by breaking bonds, by moving bonds or by a combination of both. How the additional degree of freedom provided by bond mobility changes the mechanics of peeling is not understood. To address this, we develop a theoretical model coupling diffusion, reactions and mechanics. Mobility and reaction rates determine distinct peeling regimes. In a diffusion-dominated Stefan-like regime, bond motion establishes self-stabilizing dynamics that increase the effective fracture energy. In a reaction-dominated regime, peeling proceeds by travelling fronts where marginal diffusion and unbinding control peeling speed. In a mixed reaction–diffusion regime, strengthening by bond motion competes with weakening by bond breaking in a force-dependent manner, defining the strength of the adhesion patch. In turn, patch strength depends on molecular properties such as bond stiffness, force sensitivity or crowding. We thus establish the physical rules enabling tunable cohesion in cellular tissues and in engineered biomimetic systems.


Introduction
Cell-cell adhesion is an essential mechanical function that is required to maintain tissue integrity under mechanical stress [1][2][3], disrupted during cancer [4] and finely tuned during development [5][6][7] or wound healing [8]. Cell-cell adhesion needs to manage a contradiction between stability and malleability. Such versatile and adaptable interfaces avoid unspecific adhesion [9], and instead rely on the collective effect of weak transmembrane bonds, notably of the cadherin family [10]. Cell-cell adhesion is a multi-scale and highly regulated function that involves bond clustering, coupling to the cytoskeleton through mechanosensitive adapter proteins, turnover through endocytosis [11,12] and Ca 2+ -mediated control of the molecular properties of the binders and bonds, including diffusivity [10,13], stiffness [13] and force sensitivity [14]. The distinguishing physical feature of cell-cell adhesion as compared with cell-matrix adhesion, and in general compared with conventional specific adhesion at solid-solid or solid-fluid interfaces, is that both free binders and bonds are embedded in fluid membranes and hence are laterally mobile. As a result, the dynamics of adhesion between cells, and more generally between fluid membranes bridged by transient bonds, depend on binding/unbinding reactions between partner molecules and on the lateral motion of bonds and free binders. Despite this fact being long acknowledged [9,[15][16][17], its consequences for the dynamics of peeling are not understood and a mapping of the different dynamical scenarios of decohesion is lacking.
To isolate the physical aspects of cell-cell adhesion, previous studies have focused on minimal artificial models based on lipid membranes decorated with adhesion molecules [18][19][20][21][22][23][24]. While the theoretical understanding of equilibrium in such systems is established [9,25], adhesion dynamics under force have been barely studied theoretically even though the mechanical environment of cell-cell adhesions is fundamentally dynamical. Here, we focus on the dynamics of unbinding of two vesicles held together by an adhesive patch made of mobile adhesion molecules forming transient trans bonds. Our model is also pertinent to the forced unbinding of adjacent cells, since adhesion molecules attached to the cortex are still mobile owing to the turnover of cortical components.
During peeling, an adhesion patch shrinks, possibly until complete separation. In a tear-out limiting scenario, shrinkage of the patch may proceed by sequential bond breaking [15,16,[26][27][28][29][30]. When a vesicle with mobile binders adheres to a substrate with fixed receptors, spreading critically depends on diffusion of free binders on the vesicle [31][32][33]; however, for immobile bonds, peeling necessarily proceeds by tear-out. Peeling by bond breaking has been extensively studied theoretically [15,[34][35][36]. In a competing limiting scenario, the patch may shrink by lateral motion of bonds leading to an increasingly crowded patch [16,17], a situation observed during cell-cell separation in vitro and in developing embryos [7,[37][38][39]. In general, the dynamics of decohesion may depend on a combination of bond breaking and bond motion (or reaction and diffusion), but this interplay has not been systematically examined [15,16,27,40] despite the fact that bond mobility has been shown to strongly influence adhesions in hybrid cell-supported bilayer studies [41] and in purely artificial systems [21].
To understand the physical principles governing peeling of an adhesive interface bridged by mobile bonds, we develop a self-consistent continuum dynamical model capturing the reaction kinetics of bond formation and dissociation, the lateral diffusion of adhesion molecules and the mechanics of the adhesion patch and of the adhering vesicles. We then identify and characterize distinct regimes pertinent to (i) long-lived mobile bonds, (ii) short-lived bonds with reduced mobility such as cadherin molecules linked to the cell cytoskeleton, and (iii) short-lived mobile bonds such as cadherins in a lipid bilayer. Each of these regimes exhibits fundamentally different dynamics (self-similar, travelling or mulitphasic) with multi-scale features in space and time. We further examine the relation between molecular properties of bonds and the effective behaviour of the adhesive patch.

Theoretical and computational model
The state of the system is defined by the shape of the adhering vesicles and the number concentration of bonds on the adhesion patch, c 1 , and that of free binders on the entire vesicle, c 2 . We assume that the adhering vesicles are identical and are made of a fluid membrane where bonds and free binders are mobile. Focusing first on a dilute limit and non-compliant bonds, the chemical potentials of bonds and free binders take the form where m 0 i is the standard chemical potential, k B T is the thermal energy scale and c 0 is an arbitrary reference concentration; see [42] for a detailed equilibrium statistical mechanics treatment. In equilibrium, shape and concentrations obey chemical and mechanical equilibrium conditions [25]. Chemical equilibrium requires that μ 1 (and hence c 1 ) is uniform on the adhesion patch, that μ 2 (and hence c 2 ) is uniform on the entire vesicle including the adhesion patch and that μ 1 = 2μ 2 over the patch since two free binders react to form a bond. This last condition implies that with K the non-dimensional equilibrium constant of the binders in the two-dimensional environment of the adhesion patch [43], further discussed below. At high membrane tension γ expected for a mature adhesion patch rich in bonds, the length scale , where κ ∼ 10 −19 J is the bending stiffness, is much smaller than the vesicle size. In this capillary limit, mechanical equilibrium at the edge of the adhesion patch is formally a Young-Dupré equation, k B T c 1 = 2γ(1 − cosθ), where the left-hand side is the dilute approximation of the osmotic tension of the bonds in the adhesion patch and θ is the contact angle between the free part of the vesicle and the symmetry plane [17,25] (figure 1a). Being commensurate with the osmotic tension of bonds, we estimate membrane tension for a typical bond concentration of c 1 4 Â 10 3 molecules=mm 2 to be γ ∼ 1.6 × 10 −5 N m −1 and hence ℓ 1 ∼ 70 nm. With the precise parameters considered in our examples below, ℓ 1 is close to 20 nm. Besides justifying the capillary limit, high tension suppresses thermal fluctuations of the membrane [44,45], which are not considered here but play a prominent role during low-tension adhesion spreading [21,24]. Mechanical equilibrium in the non-adhered part of the vesicle is expressed by Laplace's law. These conditions determine the equilibrium shape of the vesicles, the size of the adhesion and the concentration of free binders and bonds.
Describing the out-of-equilibrium dynamics under force requires accounting for diffusion and chemical kinetics, which in turn depend on mechanics in different ways. Diffusion of bonds is biased by their tendency to leave regions where they are highly stretched. Chemical kinetics are influenced by mechanics since unbinding rates are force sensitive and rebinding rates depend on the distance between potential partners [34,36]. We thus need to resolve the force borne by bonds and the separation between membranes required to compute these rates. However, a single mechanical modelling resolving simultaneously the vesicle-scale capillary mechanics at length scales of tens of micrometres and the separation profile near the edge of the adhesion patch with sub-nanometre resolution is very challenging from a computational point of view. For this reason, exploiting scale separation, we combine a vesicle-scale capillary model in terms of macroscopic quantities such as patch size, contact angle or vesicle pressure with a model for the micro-mechanics of the adhesion patch, taking the macroscopic quantities as parameters and resolving the force on bonds by accounting for the bending rigidity of the membrane κ and the compliance of the molecular bonds (figure 1c,d). In this model, the length scale over which the tension of the freestanding membrane is transmitted to the adhesion patch can be estimated by balancing the bending and the bond normal pressures as ' 2 ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ½4k=ðkc 1 Þ p [46], where k is the bond stiffness and c 1 is a typical bond concentration. Considering k = 2.5 × 10 −4 N m −1 and the values given above for κ and c 1 , we find that ℓ 2 ∼ 20 nm, much smaller than the typical size of an adhesion patch.
To focus on the mechano-chemistry of forced decohesion and to simplify all other aspects of the model, we restrict ourselves to a two-dimensional geometry where the membrane royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20220183 becomes a line whose arc-length coordinate is denoted by s. We summarize here the theoretical model and provide a detailed derivation based on Onsager's variational formalism in electronic supplementary material, note 1. The vesicle is connected to a loading device, which controls membrane tension γ by drawing or supplying membrane length (area in three dimensions) and applies a vertical force F (figure 1c). The area enclosed by the curve defining the two-dimensional vesicle is kept constant by a pressure difference P. The mechanical ensemble thus simulates in two dimensions a vesicle whose enclosed volume is constant and which is aspirated by a micropipette that controls membrane tension [37][38][39]. Mechanical relaxation is much faster than chemical relaxation, and thus we treat the mechanics quasi-statically. Given the prescribed γ and F and the current size of the adhesion patchŝðtÞ, the capillary model provides the shape of the vesicle (figure 1b), in particular the contact angle θ and the pressure P. With this information and the current concentration of bonds c 1 (s, t), the micromechanical model provides the membrane separation profile h(x, t) (figure 1d), by minimizing the functional capturing the mechanics of the adhesive patch Ph dx À g sin u hðaŝÞ: ð2:2Þ The first term is a capillary energy and results from ds ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ h 02 p dx. The second term is then bending energy, since the curvature of the curve can be written as h 00 /(1 + h 0 ) 3/2 . The third term is the energy stored in the elastic foundation, whose stiffness is given by the number concentration of bound binders c 1 multiplied by the stiffness to stretching of one bond, k. The fourth term accounts for the pressure inside the vesicle, which presses against the molecular bonds. The last term is the potential energy of the tension force at the boundary of the domain. α > 1 is a parameter such that ða À 1Þŝ . ' 1 , and hence the boundary of the micro-mechanical model in the top-right of figure 1d is not affected by bending elasticity and the results are insensitive to this parameter. Minimization of equation (2.2) allows us to compute the out-of-plane force per molecule k h(s, t) within the patch region. In turn, this information allows us to evolve the bond and free-binder concentrations c 1 (s, t) and c 2 (s, t) and the position of the interfaceŝðtÞ as discussed next.
The reaction-diffusion dynamics for c 1 and c 2 are given by # 0 þk on c 2 2 À k off c 1 in ð0,ŝðtÞÞ; ð2:3Þ where dots and primes denote time and space derivatives, D 1,2 are diffusion constants of bonds/free binders, k on is the binding rate, k off is the unbinding rate, L 0 is the total membrane length and is the scale of thermal fluctuations of binders. These partial differential equations are defined on a time-dependent domain. The transport term in equation (2.3) includes a diffusive term and a bias, according to which bonds try to reduce the mechanical contribution of their chemical potential [9,40], by moving away from regions where they are highly stretched. An additional consequence of the stretching term in the chemical potential of bonds identified by our model is the dependence of the binding rate on separation, in agreement with previous treatments [9,34,36]   where k on is the reference binding rate for h = 0. This expression captures the fact that the probability of bond formation depends on the proximity of potential partners as they thermally fluctuate. By equating μ 1 (c 1 , h) in equation (2.6) with twice the chemical potential of free binders, m 2 ðc 2 Þ ¼ m 0 2 þ k B T log c 2 =c 0 , we can identify a separation-dependent two-dimensional equilibrium constant KðhÞ ¼ K exp½Àðh=x g Þ 2 , where K is given by equation (2.1). Interestingly, previous theories accounting for orientational entropy and membrane thermal fluctuations but not for bond extensional compliance have related three-dimensional to two-dimensional equilibrium constants to find the same dependence of K(h) on separation with a different interpretation for x g [43,47]. In the present setting of high tension and forced peeling, membrane fluctuations may not be as important as bond deformability, but in general one can expect that the effective compliance of the binding domain depends on both the intrinsic compliance of the adhesion molecule and that of the membrane supporting the bond, which for adhesion molecules attached to the actin cortex includes adaptor molecules and the actin network. This suggests that x g , or alternatively k, should be regarded as an effective parameter rather than one strictly describing the stiffness of an adhesion molecule.
Although cadherin bonds are thought to shift between ideal, slip or catch bonds depending on environmental conditions and conformation [14,48], here we only consider the slip-bond behaviour as described by Bell's model [49] where f b is the force sensitivity and where we introduce a separation sensitivity The governing equations (2.3)-(2.5) need to be supplemented by initial, boundary and interface conditions at s ¼ŝðtÞ. Since free binders can cross the interface, their concentration and flux are continuous. By contrast, the interface is by definition a barrier for bonds. Consequently, the diffusive flux of bonds at the interface must be compensated by bond transport owing to interface motion, where now c 1 depends on space and time and is governed by equation (2.3). Comparison of this equation with Rivlin's classical theory of peeling [50] shows that osmotic tension of bonds at the interface, k B Tc 1 ðŝ, tÞ, plays the role of the adhesion fracture energy. However, rather than a material property of the interface as in classical peeling and in the case of tear-out of a vesicle against a substrate with immobile receptors [17], here this quantity is a dynamical variable. Equations (2.3), (2.4), (2.9) and (2.10) supplemented by the initial and boundary conditions at s = 0 and s = L 0 allow us to solve for c 1 (s, t), c 2 (s, t) andŝðtÞ. As these variables evolve, we need to update the mechanical variables θ and h(s, t), which in turn affect the reaction-diffusion interface dynamics. The self-consistent finite-element numerical solution of the model is described in detail in electronic supplementary material, note 1.

System preparation and parameters
Before driving the system out of equilibrium, we prepared the system at an equilibrium state for non-compliant and ideal bonds. In all calculations, we set c 0 ¼ 2:5 Â 10 3 molecules=mm 2 Â ' lat , and the length of the half vesicle to L 0 = 35 μm. Thus, the total number of molecules is N tot-= c 0 × L 0 . Here, ℓ lat is an arbitrary depth of our one-dimensional membrane to make it a ribbon, allowing us to map twodimensional to one-dimensional number densities. Without loss of generality, m 0 1 ¼ 0 and, from equation (2.1), we conclude that m 0 2 ¼ ðk B T log KÞ=2. In all figures, we non-dimensionalize chemical potentials by m 0 2 . With these data, conservation of the total number of molecules, N tot ¼ c 1ŝ0 þ c 2 L 0 , and the law of mass action, c 0 c 1 =c 2 2 ¼ K, provide two equations to solve for the equilibrium concentrations, obtaining c 1 = 1.58c 0 and c 2 = 0.89c 0 . Since R 0 andŝ 0 determine the contact angle as sin u 0 ¼ŝ 0 =R 0 , the force balance at the interface provides an equation for the membrane tension, for which, taking k B T = 4.11 × 10 −21 J, we find γ = 2.55 × 10 −4 N m −1 × ℓ lat . This equilibrium state is illustrated in figure 1a.
Starting from this state and driving it out of equilibrium by suddenly increasing the applied force F, we tracked the mechano-chemical dynamics of forced decohesion, which tend to uniformize μ 1 and μ 2 over the patch and vesicle, and tend to equilibrate them as μ 1 = 2μ 2 over the patch. The time scale of bond diffusion in the adhesion patch is t diff,1 ¼ŝ 2 0 =D 1 , whereas that of free binders on the entire vesicle is t diff,2 ¼ L 2 0 =D 2 . Since a bond connects two binders, in the simplest situation its mobility is half that of a free binder, hence D 2 = 2 D 1 . With our choice of parameters, τ diff , 2 ≈ 100 τ diff , 1 . Regarding reactions, the natural time scale is t reac ¼ 1= k off . Once k off is fixed, we determine the binding rate from K=c 0 ¼ k on = k off , which results from equations (2.1) and (2.3). The ratios between the reaction time scale and the diffusive time scales introduce Damkohler numbers weighing the relative importance of reactions and diffusion. For compliant bonds/binders, we consider a reference stiffness of k = 2.5 × 10 −4 N m −1 .

Diffusion-dominated regime
We first focused on the situation in which reactions are extremely slow compared with diffusion by setting k on and k off to zero and thus τ reac = +∞. In this limit, the dynamics of free binders become uncoupled from the rest of the system. Thus, we are only left with one time scale associated with the diffusion of bonds. For a typical diffusion coefficient of free binders on lipid membranes, D 2 = 2D 1 = 0.5 μm 2 s −1 , this time scale is τ diff,1 ≈ 25 s. We first considered non-compliant bonds by setting h(s, t) = 0; see electronic supplementary material, note 1. For this model, the state in figure 1a is an exact equilibrium state at F = 0. We then suddenly applied a separation force F. We reasoned that, in response to force application, the contact angle should rapidly increase θ > θ 0 , requiring a concomitant increase in bond concentration at the interface according to equation (2.10). In turn, this should create a sharp positive gradient in bond concentration and chemical potential, and hence a motion of the interface leading to patch shrinkage (equation (2.9)). As the patch becomes smaller, bond chemical potentials should equilibrate faster since they diffuse over a smaller distance and, since the number of bonds is constant, the patch should become increasingly concentrated with bonds. Our simulations confirmed this physical picture (figure 2a,b; electronic supplementary material, figure S1a,c and movie S1), and the system reaches a new equilibrium state where μ 1 is high and uniform and the higher bond osmotic tension balances the larger out-of-equilibrium membrane force at the interface royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20220183 (equation (2.10)). Hence, the system is able self-adjust the effective adhesion fracture energy, k B Tc 1 ðŝ, tÞ, to balance the higher peeling driving force. Although the final chemical potential of bonds is much larger than that of free binders, equilibration between species is not possible in the absence of reactions ( figure 2a,b).
As F increases, θ also increases and a larger osmotic tension, and hence a larger bond concentration, is required to balance the mechanical force at the interface. Since the number of bonds is constant, this in turn requires a smaller equilibrium patch (electronic supplementary material, figure S1b, figure 2c). Interestingly, a similar process of shrinkage and concentration of adhesive patches has been observed in cell doublets under force in vitro, and during cell-cell hydraulic fracture in developing embryos [7,25].
Increasing the force also leads to faster dynamics (figure 2c; electronic supplementary material, movie S1). To further understand peeling dynamics, we sought an analytical solution. Starting from the equilibrium state in figure 1a with contact angle θ 0 , patch sizeŝ 0 and uniform bond concentration c 1 (s, 0) = C 0 = 2γ(1 − cosθ 0 )/(k B T ), we suddenly increased the force and hence θ > θ 0 . In the actual system, asŝðtÞ decreases so does θ(t) owing to vesicle capillarity (figure 2d, inset). To develop the analytical solution, we simplified the problem by assuming that θ remains constant during peeling. Introducing non-dimensional space x ¼ s=ŝ 0 À 1, time τ = t/τ diff, 1 , driving parameter U = (1 − cosθ)/(1 − cosθ 0 ), position of the interface XðtÞ ¼ŝðt diff,1 tÞ=ŝ 0 À 1 and concentration of bonds uðx, tÞ ¼ c 1 ðŝ 0 ðx þ 1Þ, t diff,1 tÞ=C 0 À U, the governing equations (2.3), (2.9) and (2.10) take the form of a classical Stefan problem; see electronic supplementary material, note 2, which describes a myriad of phase transformation problems [51]. This problem admits an analytical solution valid at short times or for large domains in terms of the similarity variable x= ffiffi ffi t p , uðx, tÞ ¼ ffiffiffi ffi p p l e l 2 U erfðlÞ À erf x 2 ffiffi ffi t p ð3:1Þ and where erf is the error function and λ is a constant depending on the driving parameter U and is implicitly determined by ffiffiffi ffi p p le l 2 U½erfðlÞ þ 1 ¼ 1 À U. To linear order in U − 1, lðUÞ % ÀðU À 1Þ= ffiffiffi ffi p p / ðcos u À cos u 0 Þ, showing how the interface motion depends on the horizontal force imbalance. Interestingly, the diffusion-controlled spreading of a membrane with mobile molecules binding to fixed receptors exhibits analogous self-similar dynamics [31,32].
We then compared these analytical predictions, valid at short times/large domains and assuming that θ remains constant, with our numerical calculations, which did not make any of these assumptions. At short times, equation (3.2) predicts very well the motion of the interface, particularly at high forces where the interface moves significantly before the finite-size effects leading to self-stabilization of the interface start to play a role (figure 2c). To further test the theory, we plotted the rescaled bond concentration at different time instants against the similarity variable, finding a remarkable  Next, we examined numerically the more realistic situation of deformable bonds while keeping all other model parameters equal, including the total number of bonds. In this case, the state in figure 1a is not in equilibrium for F = 0 since bonds are stretched in a boundary layer near the interface, which modifies their chemical potential (equation (2.6) and figure 2e). Despite the small size of the perturbed region, the new equilibrium state is significantly smaller and more concentrated (figure 2f; electronic supplementary material, figure S1c,d). The width of the boundary layer where bonds are stretched and depleted is commensurate with ℓ 2 and thus decreases with increasing bond stiffness, but the new position of the interface is quite insensitive to k (figure 2f ). Interestingly, the results for noncompliant bonds are not the limit as k → +∞ of those for compliant bonds. We note, in this respect, that in this limit both h and x g tend to zero, and hence the last term in equation (2.6) does not necessarily tend to zero to recover the non-compliant case. Upon force application following equilibration, the dynamics proceed similarly to the case of non-compliant bonds, albeit at a faster rate (figure 2g; electronic supplementary material, figure S1c,d and movie S2). To understand the faster dynamics, we note that, close to the edge of the adhesion patch, the entropic and stretching components of the chemical potential of bonds (equation (2.6)) compete and become similar (electronic supplementary material, figure S1e), with bonds tending to move towards the depleted edge because of entropy and tending to move away from it to avoid stretching. During peeling, bond transport due to stretching dominates at the edge, driving inward interface motion (equation (2.9)). Thus, the effect of bond stretching is limited to a very small region close to the edge but has macroscopic effects in setting the equilibrium patch size and bond concentration and the kinetics of the peeling process. Although the similarity solution does not account for bond compliance, bond concentration at different times remarkably collapses when rescaled and plotted in the similarity variable (figure 2h), showing that the Stefan problem captures the physics of this diffusion-dominated regime at short times.
We finally examined the force distribution in the adhesion patch. The out-of-plane traction k c 1 (s, t)h(s, t) needs to balance F and hence increases with it. However, increasing F does not lead to further separation. Instead, the system adapts to higher F (higher θ) by increasing c 1 (s, t) near the interface (equation (2.10)), while keeping a rather constant profile of the force per molecule, k h(s, t), with a force scale emerging from the competition of mixing entropy and bond stretching (equation (2.6)) and given by f g ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi k Á k B T p % 1:0 pN (electronic supplementary material, figure S2).

Reaction-dominated case
We then studied a different extreme scenario characterized by fast reaction rates, k off ¼ 10 s À1 , representative of weak bonds such as cadherins, and very low diffusivity, which should result in a reaction-dominated regime similar to the tear-out of an adhesive vesicle from a solid substrate with immobile receptors [28][29][30]. We decreased D 1,2 by a factor of about 1000 (D 1 = D 2 /2 = D 0 = 0.25 × 10 −3 μm 2 s −1 ), which is comparable to the reduction of diffusivity of adhesion molecules from artificial lipid bilayers to cell membranes [52], and initially considered rather insensitive slip bonds with f b =f g ¼ 4. Upon force application, we observed that, in contrast to the diffusion-dominated case whereŝðtÞ Àŝ 0 / ffiffi t p , now both the size of the patch and the number of bonds decrease nearly linearly, with the net unbinding reaction rate localized in close vicinity to the interface (figure 3a; electronic supplementary material, movie S3). Complete decohesion is reached before significant diffusion of bond or free binders can take place at the patch or vesicle scales.
In view of these results, we hypothesized the existence of travelling solutions of the form c i ðs, tÞ ¼ f i ðs þvtÞ (figure 3b). To systematically examine this point, we again made the approximation of constant driving force (constant θ). Propagating fronts in reaction-diffusion systems requires non-generic nonlinearity, as in the prototypical Fisher-Kolmogorov-Petrovskii-Piscunov (FKPP) equation [53,54] or in the FitzHugh-Nagumo system [55]. For non-compliant ideal bonds, our model in the moving frame reduces to a twospecies advection-reaction-diffusion equation whose only  nonlinearity is the term k on c 2 2 , and our simulations did not develop travelling solutions. Instead, for compliant and force-sensitive bonds the model couples to mechanics, which predicts a separation profile h localized near the adhesion edge that in turn biases bond motion away from the edge, locally increases off rates and decrease on rates. In this case, our simulations readily developed travelling solutions with constant interface velocity, localized unbinding, a sharp transition of bonds to accommodate the concentration at s ¼ŝðtÞ (equation (2.10)) and a sharp transition of free binders to a higher plateau in the wake of the interface due to broken bonds (figure 3c; electronic supplementary material, figure S3).
To understand the physics controlling the peeling speedv and the local profiles around s ¼ŝðtÞ, we reasoned that, even though the time of bond diffusion in the whole patch is very long, τ diff,1 ≫ τ reac , there should be a small length scale ℓ 3 over which diffusion and reactions compete near the interface, For our parameter set, ℓ 3 is smaller than ℓ 2 , over which bonds are loaded. Dimensional analysis and the analogy with the FKPP model suggest that this small-scale diffusion may control the overall decohesion process by setting the front speed This situation is also analogous to fracture mechanics of solids, where the physics within a small process zone in the vicinity of the crack tip determines the effective fracture properties [56]. Our simulations confirmed that v 0 and ℓ 3 provide order-of-magnitude estimates of peeling velocity and size of disturbed region ( figure 3c,d). However, the system exhibits a more complex behaviour that depends on a delicate interplay not only of bond transport and reactions but also of mechanics through h (equations (2.3), (2.7) and (2.8)). This interplay controls the dynamical organization in the moving process zone (electronic supplementary material, figure S3) and ultimately front velocity. As in the diffusion-dominated case, we found that dynamics are faster for a larger mechanical driving force θ. However, herev is nearly proportional to the vertical component of force (sin θ) (figure 3d, inset). We finally examined the influence of the slip bond effect, finding a very strong increase ofv as bonds became more sensitive to force, i.e. when f b becomes smaller than the force scale close to the edge given by f g (figure 3e).
In summary, unlike the non-local and self-stabilizing peeling dynamics of the diffusion-dominated regime, here the front moves at a constant speed that depends only on the driving force and material parametersvðu; D i , k, f b , k off , k on Þ akin to the kinetic law for the motion of a material interface [57]. The unconventional tear-out described here also differs from the classical tear-out for immobile bonds [15] in that it fundamentally depends not only on marginal unbinding but also on small-scale diffusion near the front.

Reaction-diffusion regime
To examine an intermediate regime, we kept the off-rate k off ¼ 10 s À1 representative of weak bonds such as cadherins and considered diffusion constants typical of adhesion molecules on lipid membranes, D 2 = 2D 1 = 0.5 μm 2 s −1 . In this regime, depending on the magnitude of the applied force, the adhesive patch can either fail or reach a stable configuration with full uniformization and equilibration of the chemical potentials of bonds and free binders (electronic supplementary material, movie S4). With our choice of k off , reactions take place much faster than diffusion, τ reac ≈ 0.1 s ≪ τ diff,1 ≈ 25 s, and thus chemical potentials between bonds and free binders locally equilibrate very quickly in the adhesion patch (electronic supplementary material, figure S4a). For longer lived bonds, e.g. k off ¼ 0:1 s À1 , reactions and diffusion dynamically compete (electronic supplementary material, figure S4b), leading to stronger 100 Going back to the case of fast reaction rates ( k off ¼ 10 s À1 ) and first focusing on compliant ideal bonds (f b ¼ þ1), we studied in detail the peeling dynamics, which exhibit multiple scales in space and time, as shown the kymographs (figure 4b-e). At short times commensurate with the diffusion time scale in the patch τ diff,1 , the dynamics proceed similarly to a self-similar diffusion-dominated regime (figure 4d,f, inset), but now, as in the reaction-dominated regime, fast unbinding localizes at the edge of the patch to reduce the chemical potential of stretched and concentrated bonds (figure 4e), leading to a fast initial decrease of the total number of bonds (figure 4g, inset). Examining early times, t = t 1 , t 2 , even if μ 1 and μ 2 equilibrate and uniformize within the adhesion patch, the local excess chemical potential of free binders μ 2 resulting from fast unbinding has not had time to equilibrate in the rest of the vesicle, driving diffusion of free binders away from the patch ( figure 4a,b), thereby reducing μ 2 in the patch and further driving unbinding reactions (figure 4g). This process is much slower since it is controlled by the diffusive time scale over the vesicle, τ diff,2 = 100τ diff,1 . Since τ diff,2 ≫ τ diff,1 ≫ τ reac , now reactions take place nearly uniformly in a quasi-equilibrated adhesion patch. With fewer bonds, mechanical equilibrium at the interface requires reducing the size of the patch, which decreases the contact angle and increases bond concentration (equation (2.10)). In turn, the higher bond concentration favours further unbinding. Thus, the much slower dynamics during this second phase are complex, multi-phasic and depend on the diffusion of free binders over the entire vesicle with time scale τ diff,2 (figure 4f,g).
The speed and outcome of these dynamics depend on the magnitude of F. To characterize adhesion strength, we simulated the dynamics for different forces and tracked the time to complete failure t fail (figure 4h), finding that lifetime very abruptly increases as force is reduced [36]. This suggests the existence of a threshold force below which the patch is long-lived and above which decohesion occurs rapidly. Consistent with this, lifetime closely follows a power law t / ðF À F c Þ a with a ≈ −2.2 and the critical force F c a fitting parameter ( figure 4h, inset). Thus, F c can be interpreted as the strength of the adhesion patch. This mesoscopic notion of strength should depend on the microscopic strength of individual bonds given by f b . As in the FKPP regime, we found that force sensitivity of slip bonds only plays a significant role when f b , f g , in which case lifetime at fixed F and strength dramatically reduce (figure 4h).
Up to now, our model assumes a dilute limit of molecules on the membrane. However, force application leads to increasing molecular crowding, which should affect the  royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20220183 dynamics of decohesion by changing the interplay between reactions and diffusion. Crowding has been shown to facilitate passive endocytosis of particles with mobile binders [58] and to modify the mechanisms and morphology of nascent adhesions [59]. We thus developed a model accounting for crowding by considering a more general expression for the mixing entropy based on the Flory-Huggins theory that included cross-diffusion effects (electronic supplementary material, note 1). According to this model, close to a maximum concentration of molecules, c max , the chemical potential of bonds rapidly increases, which accelerates unbinding reactions and increases the effective diffusion coefficient. Close to the crowding limit, the osmotic tension of bonds takes the form −k B T(2c maxlog [(c max − c 1 )/c max ] + c 1 ), which blows up as c 1 → c max and reduces to the van't Hoff form k B T c 1 used before in the dilute limit c 1 ≪ c max . For high crowding (c max /c 0 = 5), our simulations show that concentrations uniformize much more rapidly than in the dilute limit and to a lower value as saturation is reached in the patch (figure 5a). The patch rapidly shrinks until failure because of unbinding reactions taking place throughout the adhesion driven by the higher chemical potential of bonds in a crowded situation (figure 5b-e). Systematic simulations at various forces show that crowding accelerates failure and reduces F c analogously to the effect of slip bond sensitivity (figure 5f ). Hence, crowding weakens the adhesion patch by favouring a tear-out mechanism with diffuse unbinding (figure 5c), different from the FKPP regime exhibiting highly localized marginal unbinding (figure 3a). For cadherins on lipid vesicles, c max /c 0 ≈ 20 [60], our model still predicts significant weakening caused by crowding relative to the dilute limit. Our analysis does not include the effect of other membrane molecules, which can act as crowders but not contribute to adhesion.

Summary and discussion
In summary, we have developed an out-of-equilibrium model self-consistently coupling diffusion, binding and unbinding reactions and mechanics to understand the dynamics of peeling between fluid membranes bridged by mobile adhesion molecules forming transient bonds. We have used this model to map various distinct and biologically relevant scenarios of forced decohesion amenable to experimental examination. In all of these scenarios, the macroscopic peeling behaviour depends on the physics occurring in a very small process zone close to the edge of the adhesion patch.
(i) For long-lived mobile bonds, adhesion patches shrink and become concentrated in a self-stabilizing process controlled by diffusion. At short times, the system evolves according to the self-similar dynamics of a classical Stefan problem with the interface moving as ðŝðtÞ Àŝ 0 Þ / ffiffi t p . (ii) For short-lived bonds with low diffusivity, such as cadherins partially immobilized by the cytoskeleton, we have identified a new unconventional tear-out regime characterized by FKPP-like travelling solutions with ðŝðtÞ Àŝ 0 Þ / t, which are localized reactions in the vicinity of the interface, but also by small-scale diffusion in a process zone of size ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi D 1 = k off q . The interplay between diffusion and reactions sets the order of magnitude of the front speed ffiffiffiffiffiffiffiffiffiffiffiffi ffi D 1 k off p , but this speed is strongly influenced by the applied force and by the ratio between force sensitivity f b and the characteristic force born by bonds close to the interface f g ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi k Á k B T p . (iii) For mobile short-lived bonds such as cadherins on a lipid membrane, the system exhibits a hierarchy of reaction and diffusion time scales resulting in multi-phasic dynamics. The reinforcing effect of bond motion and the weakening effect of bond breaking compete in a force-dependent manner, defining the strength of the patch below which peeling arrests and above which peeling rapidly leads to complete failure. Strength strongly decreases for sensitive bonds (f b , ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi k Á k B T p ) and with molecular crowding. Although our minimal model ignores important aspects of cell-cell adhesion, the physical rules identified here should bear biological relevance. We have shown how the ability of bonds to laterally move in fluid-fluid adhesive interfaces leads to a very rich repertoire of peeling scenarios that cells can use to stabilize cell-cell junctions during physiological stretch, or to selectively detach during morphogenesis. For instance, cells can effectively tune adhesive strength, and hence their ability to stay adhered or to disengage, by controlling molecular properties of bonds such as stiffness k and force sensitivity f b , e.g. through extracellular Ca 2+ , by controlling the number of transmembrane crowding molecules or by controlling the actively generated surface tension. Beyond cells, our study also provides a conceptual framework for artificial biomimetic systems with a comparable degree of adhesive tunability [45].
Data accessibility. The Matlab computer codes used to generate the results in the manuscript are available from: https://github.com/ pradeep927/Peeling_dynamics_matlab.git.
The data are provided in the electronic supplementary material [61].