Criticality in Cell Adhesion

We illuminate the many-body effects underlying the structure, formation, and dissolution of cellular adhesion domains in the presence and absence of forces. We consider mixed Glauber-Kawasaki dynamics of a two-dimensional model of nearest-neighbor interacting adhesion bonds with intrinsic binding-affinity under the action of a shared pulling or pushing force. We consider adhesion bonds that are immobile due to being anchored to the underlying cytoskeleton as well as adhesion molecules that are transiently diffusing. Highly accurate analytical results are obtained on the pair-correlation level of the Bethe-Guggenheim approximation for the complete thermodynamics and kinetics of adhesion clusters of any size, including the thermodynamic limit. A new kind of dynamical phase transition is uncovered -- the mean formation and dissolution times per adhesion bond change discontinuously with respect to the bond-coupling parameter. At the respective critical points cluster formation and dissolution are fastest, while the statistically dominant transition path undergoes a qualitative change -- the entropic barrier to complete binding/unbinding is rate-limiting below, and the phase transition between dense and dilute phases above the dynamical critical point. In the context of the Ising model the dynamical phase transition reflects a first-order discontinuity in the magnetization-reversal time. Our results provide a potential explanation for the mechanical regulation of cell adhesion, and suggest that the quasi-static and kinetic response to changes in the membrane stiffness or applied forces is largest near the statical and dynamical critical point, respectively.


I. INTRODUCTION
Cell adhesion refers to the specific binding of cells to neighboring cells or the extracellular matrix. It plays a major role in cell regulation [1], intercellular communication [2], immune response [3], wound healing [4], morphogenesis [5], cellular function [6], and tumorigenesis [7,8]. Cellular adhesion domains form as a result of the association of transmembrane cellular adhesion molecules (CAMs) that interact with the actin cytoskeleton [9] and can translocate over the membrane [10]. There are four major superfamilies of CAMs; the immunoglobulins, integrins, cadherins, and selectins, and throughout we generically refer to them as CAMs. Biological adhesion bonds are typically non-covalent with binding energies on the order of a few k B T corresponding to forces on the order of 4 pN · nm at T 300 K [11,12]. As a result of thermal fluctuations these bonds have finite lifetimes -they can break and re-associate depending on the receptorligand distance, their respective conformations and local concentrations, and depending on internal and external mechanical forces [12,13]. While it was originally thought that the strength of adhesion is determined by the biochemistry of CAMs alone, more recently, cellular mechanics [14] and adhesion bond interactions induced by thermal undulations of the membrane [15][16][17][18][19] emerged as essential physical regulators of cellular adhesion.
These observations imply many-body physics to be at play, i.e. an interplay between the coupling of nearby adhesion bonds through deformations of the fluctuating membrane and mechanical forces acting on the membrane [3, 15-19, 41, 44-50]. Supporting the idea are experimental observations of cells changing the membrane flexibility and/or membrane fluctuations through ATPdriven activity [51][52][53][54], decoupling the F-actin network [55] or remodelling the actomyosin cytoskeleton [54], and through acidosis [45], in order to alter adhesion binding rates and strength [41,45,[56][57][58][59] or to become motile [60]. There is also a striking correspondence between membrane stiffness and the metastatic potential of cancer cells -the stiffness of cancer cells was found to determine their migration and invasion potential [60]. The effect is not limited to cells; the elastic modulus was similarly found to significantly affect the specific adhesion of polymeric networks [61].
Most of our current understanding of the formation and stability of adhesion clusters derives from the analysis of individual [11] and non-interacting adhesion bonds [62][63][64], and studies of collective effects in biomimetic vesicular model systems with floppy membranes [48,65] and mobile CAMs [66]. These results therefore do not  1) and (2), depicting an adhesion domain on a cell-patch with 16 CAMs anchored to stiff substrate. Adhesion bonds are arranged on a 4 × 4 square lattice and can assume two states, σi = ±1, where +1 corresponds to an open (red) and −1 to a closed bond (green). Nearest-neighbor bonds experience an effective interaction J induced by undulations of the anchoring membrane. An external force h is pulling/pushing on the adhesion domain. Each adhesion bond has an intrinsic binding-affinity µ ≥ 0 that favors a bound state. A small number of bonds is depicted for convenience only. In the work we consider different system sizes including the thermodynamic limit. (b) Glauber and Kawasaki transition. A Glauber transition changes the binding state of a single adhesion bond to σi → −σi with transition rate wi({σj}) (see Eq. (3)). A Kawasaki transition interchanges two nearest-neighbor adhesion bonds σi ↔ σ k with transition rate w ik ({σj}) (see Eq. (5)), corresponding to lateral diffusion.
necessarily apply to cells, where membranes are stiffened by the presence of, and receptors are anchored to, the stiff actin cytoskeleton that can actively exert forces on the membrane [9].
While it is omnipresent in biological systems, cell adhesion displays subtle differences in the specific microsopic details. Here we aim to capture the essential general features of the physics of cell adhesion. In order to arrive at a deeper understanding of the mechanical regulation of cellular adhesion that would explain the collective dynamics of adhesion bonds on the level of individual (un)binding events we here consider mixed Glauber-Kawasaki dynamics of a generic, two-dimensional model of diffusing nearest-neighbor interacting adhesion bonds with intrinsic affinity µ under the action of a shared force h (see Fig. 1a).
Highly accurate analytical results on the Bethe-Guggenheim level reveal the many-body (that is, beyond "mean field") physics underlying biological adhesion. We consider in detail cluster-sizes ranging from a few CAMs to the thermodynamic limit. In the thermodynamic limit we determine the equation of state and complete phase behavior that displays a phase separation and co-existence of dense and dilute adhesion domains. The critical behavior is investigated in detail and striking differences are found between pulling-and pushingforces. Strikingly, we prove the existence of a seemingly new kind of dynamical phase transition -the mean first passage time to cluster formation/dissolution is proven to change discontinuously with respect to the coupling strength. This dynamical phase transition, and more generally the non-linear and non-monotonic dependence on the membrane flexibility, may explain the puzzling cooperative behavior of effective association and dissociation rates measured experimentally.

A. Outline of the work
The paper is structured as follows. In Sec. II we present an effective mesoscopic model of adhesion clusters and provide a practical roadmap to the diverse calculations and analyses. In Sec. III we present explicit analytical results for the thermodynamic equation of state and complete phase behavior of adhesion clusters, and in Sec. IV we present analytical results for the kinetics of cluster formation and dissolution both in the presence and absence of forces. In Sec. V we discuss the biological implications of our results and in particular the suggestive rôle of criticality in the context of equilibrium adhesion strength and the kinetic dissolution and formation rates, respectively. Finally, in Sec. VI we highlight the relevance of our results in the context of the Ising model. We conclude in Sec. VII with a summary and a perspective on the importance and limitations of our results, and mention possible extensions to be made in future studies. Details of calculations, explicit asymptotic results and further technical information is presented in a series of Appendices.

A. Equilibrium
We consider a two-dimensional patch of a cell surface with N adhesion molecules embedded in the cell membrane, their lateral positions forming a lattice with coordination number z (see Fig. 1). The results we derive hold for any lattice but we focus the discussion mainly on the square lattice with free boundary conditions. Opposing the patch is a stiff substrate or a neighboring cellpatch with complementary adhesion molecules occupying a commensurate lattice. The state of individual bonds is denoted by σ i , i = 1, 2, ..., N , where σ i = +1 if bond i is broken an σ i = −1 if it is closed.
In the presence of a timescale separation the opening/closing of nearest neighbor bonds is coupled via membrane fluctuations. Following closely the arguments of Ref. [17] we can integrate out the membrane degrees of freedom to obtain an effective Ising-like model for the bonds within the patch with effective Hamiltonian (1) where J ≥ 0 is the membrane-induced short-range coupling between the bonds, ij denotes all nearestneighbor pairs, µ is the effective chemical potential (i.e. intrinsic affinity) of individual bonds, and H h ({σ i }) is the Hamiltonian describing the effect of the mechanical force. The first term in Eq. (1) represents the effective coupling between nearest neighbor bonds, and is isomorphic to the interaction term in the Ising model [73]. It is an effective measure of bond-cooperativity, i.e. it reflects that the (free) energy penalty of closing/breaking a bond is smaller if neighboring bonds are closed/open, respectively [17]. Such an effective description in terms of bonds coupled via a short-range membrane-mediated interaction is feasible when bonds are flexible and/or the patch of the cell membrane is quite (but not completely) stiff and is thus rather pulled down as a whole instead of being locally strongly deformed by the binding of individual bonds [17]. In this limit the coupling strength is determined by the effective bending rigidity of the cell membrane, κ, via J ∝ 1/ √ κ (see [17] and Appendix A). That is, in this regime a relatively floppier cell membrane with lower bending rigidity induces a stronger cooperativity between neighboring bonds than a relatively stiff membrane. Notably, a detailed comparison between the full model of specific adhesion (i.e. reversible adhesion bonds explicitly coupled to a dynamic fluctuating membrane) and the lattice model captured by the first term of Eq. (1) revealed a quantitative agreement (see e.g. Fig. 5 in [17]) in the range 0 ≤ J 1.2 k B T that lies entirely within the rather stiff limit [17]. This is the range of J we are interested in and includes the values relevant for cell adhesion (see Sec. V below).
The second term in Eq. (1) reflects that each closed bond stabilizes the adhesion cluster by an amount −µ. Aside from the last term H h ({σ i }) the Hamiltonian (1) is isomorphic to the lattice gas model developed in [17], and a mapping between the two models is provided in Appendix A.
The third term in Eq. (1), H h ({σ i }), accounts for the mechanical force h acting on the membrane-embedded bonds that we assume to be equally shared between all N c closed bonds of a given configuration where δ ik is Kronecker's delta. More precisely, the force h destabilizes the bound state by introducing an elastic (free) energy penalty on all closed bonds whereby broken bonds remain unaffected. If all bonds are closed, N c = N , this penalty is set to be hx 0 , where x 0 is a microscopic length-scale specific for a given CAM that merely sets the energy scale associated with the elastic strain caused by h. Conversely, the penalty must vanish in a completely dissolved configuration with N c = 0, and is assumed to be a smooth and monotonic function of N c . A mathematically and physically consistent definition is A 'pulling force', h > 0, favors the dissociation of bonds while a 'pushing force', h < 0, favors their association. We are interested in strain energies on the order of the thermal energy per bond, i.e. |h|x 0 /N = O(k B T ). Note that the assumption of an equally shared force in Eq. (2) is valid if either of the following conditions is satisfied: the anchoring membrane has a large combined elastic modulus (i.e. stiff membranes or membrane/substrate pairs), individual bonds are flexible, the bond-density is low, or the membrane is prestressed by the actin cytoskeleton [43,74,75]. In the limit of a rather stiff membrane both, a spin representation with effective coupling J and a uniform force load are valid approximations to describe cell adhesion under force over a broad range of physically relevant parameters, as we detail below. The implications of a non-uniform force load are addressed in detail in Sec. VII and Appendix E 2.

B. Kinetics
The breaking/closure and lateral diffusion of adhesion bonds are assumed to evolve as a discrete time Markov chain with mixed single-bond-flip Glauber dynamics [76] and two-bond-exchange Kawasaki dynamics [77] (see Fig. 1b). For a single jump in the Markov chain we define the probability to attempt a Glauber transition as p k ∈ [0, 1] which controls the diffusion rate, and for the sake of generality is allowed to depend on the number of closed bonds k. Similarly, the probability to attempt a Kawasaki transition is given by 1 − p k ∈ [0, 1]. We consider two distinct scenarios, one in which adhesion bonds are immobile due to being anchored to the underlying cytoskeleton (i.e. p k = 1 ∀k), and the other in which adhesion molecules are allowed to transiently diffuse (i.e. 0 < p k < 1 ∀k; see e.g. [10]). Conversely, permanently associated/dissociated freely diffusing bonds (i.e. p k = 0 ∀k) will not be considered, since these are not relevant. Further details about the respective transition rates are given below.
Glauber transitions: Let {σ j } i denote the bond configuration obtained by flipping bond i while keeping the configuration of all other bonds fixed, i.e.
the energy difference associated with the transition. These rates can be specified uniquely by limiting interactions to nearest-neighbors, imposing isotropy in position space, and requiring that w i satisfies detailed balance, i.e.
, where α is an intrinsic attempt-frequency that sets the fastest timescale [76], and time will throughout be expressed in units of α −1 . Introducing furthermore the dimensionless quanti-tiesJ = βJ,μ = βµ andh = βhx 0 /N this leads to where we defined the auxiliary function Kawasaki transitions: Let {σ j } ik denote the bond configuration upon interchanging the state of the nearest neighbor bonds σ i and σ k while keeping the configuration of all other bonds fixed, i.e.
is the energy difference associated with the transition. Imposing the same symmetry constraints as for the Glauber rates as well as detailedbalance yields the general expression [77] where we have used that (σ i − σ k )/2 ∈ {−1, 0, 1}. As pointed out in [77], the transition is only meaningful when σ k = −σ i , otherwise the transition brings the system to an identical state, which is equivalent to no transition. Note that the Kawasaki rates given by Eq. (5) do not depend on the external forceh nor the bindingaffinityμ, since the Kawasaki transition conserves the total number of open and closed adhesion bonds. However, if in addition we introduce a position-dependent force/binding affinity the Kawasaki rates also depend oñ h andμ, which we analyze in Appendix E 2. Figure 2. Strategy roadmap. Small system sizes N ≤ 5 × 5 are solved for exactly. The thermodynamics of larger systems is treated on the level of the highly accurate Bethe-Guggenheim approximation and the kinetics by assuming local equilibrium. Within the Bethe-Guggenheim approximation we take the thermodynamic (TD) limit N → ∞ and determine the phase behavior, master scaling of dissolution/formation kinetics and analyze the statical and dynamical critical behavior.

C. Strategy roadmap
We focus in detail on both, the equilibrium properties as well as the kinetics of cluster formation and dissolution for all cluster sizes. A roadmap to our extensive analysis is presented in Fig. 2.
For small to moderate cluster sizes, i.e. up to 50 bonds for the equilibrium properties and up to 25 bonds in the case of formation/dissolution kinetics, we obtain exact solutions using standard algebraic methods [78]. To circumvent the explosion of combinatorial complexity for large system sizes we employ a variational approachthe so-called Bethe-Guggenheim approximation [79] -to derive closed-form expressions for the partition function, and finally carry out the thermodynamic limit to derive explicit closed-form results for large adhesion clusters. When considering the formation/dissolution kinetics of large clusters and in particular in the thermodynamic limit, we employ the local equilibrium approximation, where we assume that the growth and dissolution evolves like a birth-death process on the free energy landscape.
We systematically test the accuracy of all approximations by comparing them with exact results for system sizes that are amenable to exact solutions. The results reveal a remarkable accuracy that improves further with the size of the system (e.g. see Fig. B1).

A. Small and intermediate clusters
In order to quantify the equilibrium stability of adhesion clusters we first analyze the equation of state for the average fraction of closed bonds, ϕ ≡ N c ({σ i }) /N at givenμ,J andh. To this end we require Q k , the partition function constrained to the number of closed bonds N c ({σ i }) = k. We therefore write the total canonical partition function Q for a system of N adhesion bonds as ,k is the partition function of the Ising model at zero field conditioned to have a magnetization N/2−k. The free energy density (per bond) in units of thermal energy k B T constrained to a given fraction of closed bonds ϕ,f N (ϕ), and the equation of state, ϕ(μ,J,h) , are given bỹ We note that e −Nf N (ϕ) /Q = Prob(N c = N ϕ) in an equilibrium ensemble of N bonds. The sum over constrained configurations in Z k contains N k terms. Whereas it can be performed exactly for N 50 it explodes for larger system sizes. To overcome the computational complexity we employ a variational approach -the Bethe-Guggenheim approximation [79], yielding (see derivation in Appendix B 1) is the average coordination number in a cluster with local coordination z i that accounts for finite-size effects, X * k ≡ k(N − k)/N , we have defined and introduced the auxiliary function (10) where Γ(z) stands for the Gamma function. Note that by settingX k = X * k in Eq. (8) we recover the mean field result Z MF k (which happens automatically forJ = 0 or k = 0 ∨ N ) which is discussed in Appendix C 1. Fig. 3(a-c) shows a comparison of the free energy densityf N (ϕ) for a cluster of 40 bonds for various affinitiesμ and external forcesh, and confirms the high accuracy of the Bethe-Guggenheim approximation on the one hand, and the systematic failure of the mean field result on the other hand. This signifies that correlations between adhesion bonds decisively affect cluster properties. Moreover, pairwise correlations captured by the Bethe-Guggenheim approach are apparently dominant, whereas three-body and higher order correlations that were ignored are apparently insignificant.
Similarly, in Fig. 3(d-f) we depict the equation of state for a cluster of 40 bonds. The Bethe-Guggenheim approximation (blue lines) is very accurate for all values of J whereas the mean field approximation (red lines) fails for intermediate values of the coupling. We observe striking differences in the dependence of ϕ on the coupling J (and hence membrane rigidity) with respect to the intrinsic binding-affinityμ in the presence of a pulling force (see Fig. 3f). At strong coupling between adhesion bonds ϕ depends strongly onμ. In the presence of a pulling force adhesion bonds with a weak affinity are on average all broken, whereas they are all closed if the affinity is large. Notably, the dependence of ϕ on the couplingJ at zero force (see Fig. 3e) agrees qualitatively well with experimental observations [21,39,68] and hints at some form of critical behavior underneath, which we discuss in more detail in Sec. V.

B. Thermodynamic limit
To explore the phase diagram in detail and analyze the critical behavior we consider the thermodynamic limit of the Bethe-Guggenheim (BG) and mean field (MF) free energy density, i.e. the scaling limit which exists and is given bỹ where X * ϕ /N = ϕ(1 − ϕ), Ω ϕ ≡X ϕ /N , and we have introduced the auxiliary function Ξ(x) ≡ x ln x. The result forf MF (ϕ) is given in Appendix C 1. Somewhat surprisingly the free energy density of a finite system, f BG N (ϕ), converges to the thermodynamic limitf BG (ϕ) already for N 100. For convenience we henceforth drop the superscript BG when considering the Bethe-Guggenheim result, i.e.f BG (ϕ) →f(ϕ).
The equation of state in the thermodynamic limit is determined by means of the saddle-point method (for derivation see Appendix D), yielding a weighted sum over ϕ 0 i , the M global minima off(ϕ): wheref(ϕ 0 i ) =f min , ∀i, and stands for asymptotic equality in the thermodynamic limit. In practice M is either 1 (unique minimum) or 2 (two-fold degenerate minima). The minima have the universal form ϕ 0 m = ξ 4 µ,J,h /(1 + ξ 4 µ,J,h ) with the coefficients ξμ ,J,h and weights c i given explicitly in Appendix D. The equation of state ϕ for a finite cluster seems to converge to the saddlepoint asymptotic ϕ TD already for N 400 for any value of the forceh, bond affinityμ, and couplingJ (see Fig. 4(a-c)), and is qualitatively the same as for smaller clusters (compare Fig. 4(a-c) with Fig. 3(d-f)). However, important differences emerge in the thermodynamic limit -the system may undergo a phase transition and phaseseparate into dense ("liquid") and dilute ("gas") phases of closed bonds with composition ϕ l and ϕ g , respectively (see also [49]).

C. Phase diagram and critical behavior
To determine the phase diagram we require the bin-odalJ b (ϕ) and spinodalJ s (ϕ) line. The binodal linẽ J b (ϕ) denotes the onset of phase separation and is determined by the "common tangent" construction, i.e. from the solution of the coupled equations where the prime denotes the derivative with respect to ϕ at constantJ. The spinodal lineJ s (ϕ), also known as the stability boundary, denotes the boundary between the metastable and unstable regimes and is determined byf (ϕ) = 0. For a non-zero force,h = 0, we determinẽ J b (ϕ) numerically, whereas we obtain an exact result for a vanishing forceh = 0 that reads (see derivation in Appendix B 2) where we have introduced χ ϕ = ϕ/(1 − ϕ). The spinodal line for any forceh is in turn given exactly bỹ with the auxiliary function . Note that it follows from their respective definitions that neitherJ b (ϕ,h) norJ s (ϕ,h) depends onμ (for a proof see Appendix B 2). The phase diagram for a pushing, zero, and pulling forceh is shown in Fig. 4(d-f) and displays, above the critical coupling strengthJ >J s crit , a phase separation into a dense and dilute phase of closed bonds with compositions ϕ l and ϕ g , respectively. A pushing forceh < 0 lifts the critical coupling and "tilts" the phase diagram towards higher density, i.e. at a given couplingJ >J s crit the density of both phases increases. Conversely, a pulling forceh > 0 lowers the critical coupling and "tilts" the phase diagram towards lower density, i.e. at a given couplingJ >J s crit the density of both phases decreases. The biological implications of these results will be discussed in Sec. V. The binodal and spinodal line in the mean field approximation are given in Appendix C 2.
We now address in detail the behavior of the statical critical point (ϕ s crit ,J s crit ) -the point where the binodal and spinodal merge, The critical point denotes the onset of phase separation and is the solution off (ϕ) = 0, which in absence of the force yields (for derivation see Appendix B 2) (ϕ s,0 crit ,J s,0 crit ) ≡ 1 2 (1, lnz z−2 ). In the presence of a forceh = 0 we obtain the exact solution using a Newton's series approach [80][81][82]   Symbols depict the exact solution using a converged Newton's series and the gradient line depicts the two-term (so-called quadratic) approximation of the complete Newton's series, which is very accurate for any pulling-and up to a moderate pushing force, i.e.h ≥ −1. Explicit expressions are given in Appendix E. The black line corresponds to the prediction of second order perturbation theory from Eqs. (18) and (19) that is valid for small forces.
and correspondingly ϕ s The dependence of the statical critical point on the external force is depicted in Fig. 5. A pulling force (red) pulls the critical point towards lowerJ and lower ϕ, whereas a pushing force (blue) effects the opposite and shifts the critical point towards larger couplingJ and higher density ϕ. The mean field statical critical point can be derived exactly as a function of the forceh, and the result is given in Appendix C 2.

A. Small and intermediate clusters
We are interested in the kinetics of cluster formation from a completely unbound state, and cluster dissolution from a completely bound state. More general initial conditions are treated in Appendix E. We quantify the kinetics by means of the mean first passage time τ d,f , where the subscripts d and f stand for dissolution and formation, respectively, and τ d,f is the first passage time defined as where {σ i } t denotes the instantaneous state at time t.
A cluster with N adhesion bonds has 2 N possible states {σ i }. We enumerate them such that the first state corresponds to all bonds closed and the final state to all bonds broken. The transition matrix of the Markov chain describing mixed Glauber-Kawasaki dynamics on this statespace has dimension 2 N × 2 N , whereby we must impose absorbing boundary conditions on the fully dissolved and fully bound states, respectively. An exact algebraic result for τ d,f is given in Eq. (E2) in Appendix E 1 but requires the inversion of a (2 N − 1) × (2 N − 1) sparse matrix, followed by a sum over 2 N − 1 terms, which is feasible only for N 5 × 5. As a result of the non-systematic cluster formation and dissolution at zero couplingJ = 0, and motivated by the intuitive idea that the dynamics is dominated by low energy (i.e. minimum action) paths at large couplingJ 1, we make the so-called local equilibrium approximation to treat large clusters. Thereby we map the dynamics of the 2 N × 2 N state-space onto a one-dimensional birth-death process for the instantaneous number of closed bonds k (see Fig. 6) with effective transition rates where we have defined the re-weighted canonical partition functionQ k ≡ Q k /p k where p k is the Glauber attempt probability in state k, and we have introduced the exit rates from configuration Note that only the Glauber transitions, given by Eq. (3), enter in Eq. (22). The Kawasaki transitions given by Eq. (5), which conserve the total number of closed bonds, enter the dynamics through the diagonal of the transition matrix as the waiting ratesw k→k = 1 −w k→k+1 − w k→k−1 , where the right hand side follows from conservation of probability. Within the local equilibrium approximation the mean first passage time for cluster dissolution and formation become, respectively where one can further use the detailed balance rela-tionQ kwk→k−1 =Q k−1wk−1→k (which we prove in Appendix E 3) to interchange the backward and forward rate in the second line and change the summation according to N k=1w In Appendix E 4 we prove that Eq. (23) holds for any birth-death process where the transition rates obey detailed balance. A comparison of the exact result given by Eq. (E2) with the local Figure 6. Mapping the full dynamics onto a birthdeath process. For convenience, and without any loss of generality, we here show an example of a system composed of 3 adhesion bonds on a 1-dimensional lattice. The mapping holds for any lattice geometry. In the full dynamics each lattice configuration represents a different node, comprising a 2 N × 2 N transition matrix, whereas in the local equilibrium approximation we only need to distinguish between states with a different number of closed/open bonds, comprising a (N + 1) × (N + 1) transition matrix. equilibrium approximation in Eq. (23) shown in Fig. 7 demonstrates the remarkable accuracy of the approximation already for N ∼ 20 bonds, which increases further for larger N . The reason for the high accuracy can be found in the large entropic barrier to align bonds in an unbound/bound state, effecting a local equilibration prior to complete formation/dissolution. Moreover, the local equilibrium approximation is expected to become asymptotically exact even for small clusters in the ideal, non-interacting limitJ → 0 as well as forJ → ∞ that is dominated by the minimum-action, "instanton" path. A further discussion of the local equilibrium approximation and an approximate closed form expression for Eq. (23) for larger systems is given in Appendices E 5 and E 6.
The mean first passage times for cluster dissolution/formation shown in Fig. 7 both display a strong and non-monotonic dependence on the coupling parameterJ with a pronounced minimum, hinting at some form of critical dynamics. As we prove below this minimum in the thermodynamic limit indeed corresponds to a dynamical critical coupling.

B. Thermodynamic limit
We now consider dissolution and formation kinetics in very large clusters, i.e. in the limit N → ∞. Note that while the mean first passage time formally diverges, i.e. lim N →∞ τ d,f = ∞, it is expected to do so in a "mathematically nice", well-defined "bulk scaling". In  Fig E3). Colored lines correspond to exact results obtained from Eq. (E2) for various values of the Glauber attempt probability p, which we set to be constant p k → p, and symbols denote the local equilibrium approximation Eq. (23) evaluated with the exact Q k andw k→k±1 from Eqs. (6) and (21) respectively.
anticipation of an exponential scaling of relevant timescales with the system size N we define the mean formation/dissolution time per bond in the thermodynamic limit as t d,f ≡ lim N →∞ τ d,f 1/N . Using the local equilibrium approximation for the mean first passage time given by Eq. (23), and assuming that the Glauber attempt probabilities p k are strictly sub-exponential in N , we prove via a squeezing theorem in Appendix E 7 that the exact mean dissolution and formation time per bond in the thermodynamic limit reads where Eq. (24) shows that the mean first passage per bond in the thermodynamic limit is determined exactly by the largest left/right-approaching free energy barrier between the initial and final point, and is completely independent of the Glauber attempt probability p k . We obtain analytical results for Eqs. (24) and (25) for arbitraryJ,μ  6) and (8)) andw k→k+1 from Eq. (21). The discrepancy between the lines and symbols is due to finite-size effects. (c) In the thermodynamic limit and more generally for large clusters the mean dissolution/formation time t d,f depends only on the largest free energy barrier (see Eq. (24)). For small coupling (regime I) the latter corresponds to the difference between the free energy minimum and the fully dissolved or bound configuration, ∆f † = ∆f † 0,1 , respectively. At the statical critical coupling value, J s crit , (onset of regime II) a free energy barrier emerges separating the meta-stable from the stable phase, ∆f † H L , but the largest free energy barrier is still ∆f † = ∆f † 0,1 . At the dynamical critical coupling,J d crit , (onset of regime III) the free energy barrier separating the meta-stable from the stable phase becomes dominant, ∆f † = ∆f † H L . The depicted free energy landscapesf(ϕ) correspond to Eq. (12) withμ = 0.05 andh = 0.
andh. Since these results are somewhat complicated for µ > 0 andh = 0 we present them in Appendix E 8 and Fig E4. In the force-free case with zero intrinsic affinity, i.e.μ =h = 0, they turn out to be surprisingly compact and given by such that forJ = 0 andJ → ∞ we have t d,f = 2 being the maximum, and the minimum occurs atJ = ln (1 Fig. 8a,b shows a comparison of the prediction of Eq. (24) with the results for finite system given by Eqs. (23) and (21) rescaled according to τ d,f 1/N . Already for N = 900 a nearly complete collapse to the thermodynamic limit (24) is observed for both, cluster formation as well as dissolution. The mean field analogue of Eq. (26) is given by Eq. (E35) for a generalz and remarkably has a universal (i.e.z-independent) minimum value of t d,f MF ≈ 1.0785 at the dynamical critical couplingJ = 2 ln (2)/z (see Eq. (E37)). Moreover, t d,f MF displays an unphysical divergence in the limitJ → ∞ (see Fig. E5).

C. Dynamical phase transition and critical behavior
Strikingly, the mean dissolution and formation time in the thermodynamic limit (24) display a discontinuity as a function of the couplingJ (see jumps in ∂J t d,f depicted in the insets in Fig. 8a,b). In particular, for zero affinity and external force we find from Eq. (26) This implies the existence of a first order dynamical phase transition at the dynamical critical couplingJ d crit and hence a qualitative change in the dominant dissolution/formation pathway. Coincidentally, the Bethe-Guggenheim dynamical critical point forμ =h = 0 coincides with the exact (Onsager's) statical critical point for the two-dimensional zero-field Ising model [83]. Similarly, the mean field dynamical critical point forμ =h = 0 coincides with the Bethe-Guggenheim statical critical point (for a more detailed discussion see Appendices E 8 and E 9). Strikingly, the dynamic critical point always corresponds to the minimum of t d,f . The explanation of the physics underneath the dynamical phase transition and the meaning ofJ d crit is given in Fig. 8c. The qualitative behavior of t d,f has three distinct regimes. In regime I, where 0 ≤J <J s crit , the free energy landscapef(ϕ) has a single well and according to Eq. (24) t d,f is determined by ∆f † 0,1 -the free energy difference between the minimum and the absorbing point (i.e. ϕ = 0 for dissolution and ϕ = 1 for formation, respectively). ∆f † 0,1 is a decreasing function ofJ. At the statical critical couplingJ s crit , which marks the onset of regime II, a second free energy barrier emerges delimiting the phase-separated low (L) and a high (H) density phase. We denote this free energy barrier by ∆f † H L where → and ← stand for dissolution and formation, respectively. ∆f † H L is an increasing function of J. In regime II, that is whenJ s crit ≤J <J d crit , the dissolution and formation first evolve through a (thermodynamic) phase transition and, finally, must also surmount the second, predominantly entropic barrier to the complete dissolved/bound state. In regime II, as in regime I, the largest free energy barrier remains the free energy difference between the minimum and the absorbing point, i.e. ∆f † 0,1 > ∆f † H L . Exactly at the dynamical critical couplingJ d crit the two barriers become identical, ∆f † 0,1 = ∆f † H L and for J >J d crit we always have ∆f † 0,1 < ∆f † H L . Therefore, in regime III the rate-limiting event becomes the phase transition itself, whereas the fully dissolved/bound state is thereupon reached by typical density fluctuations. Since ∆f † 0,1 decreases withJ while ∆f † H L increases with J, the mean dissolution/formation time per bond at the dynamical critical couplingJ d crit must be minimal. This explains the dynamical phase transition completely.
Note that the dynamical phase transition is preserved under initial conditions that lie beyond the largest free energy barrier (from the final/absorbing state). For example, we may consider ϕ( L,H is the (meta)-stable minimum in the high and low density region for cluster dissolution and formation, respectively. In the thermodynamic limit the equilibration time from the initial condition ϕ({σ i }) = 0 ∨ 1 to the (meta)-stable minimum ϕ 0 L,H becomes exponentially faster than the total transition time, which renders t d,f unaffected.
The dependence ofJ d crit onμ andh is determined in the form of a Newton's series in Appendix E 8, and is depicted in Fig. 9. Depending on the intrinsic affinitỹ µ, the dependence ofJ d crit may be non-monotonic. Note that in contrast to the statical critical couplingJ s crit that Figure 9. Dynamical critical point. Dynamical critical couplingJ d crit as a function of the external forceh for several values of the intrinsic binding-affinityμ. Note thatJ d crit as a function ofh may be non-monotonic with a global minimum whose location depends onμ. Bending rigidity βκ 4 ∼ 400 [55,[88][89][90] is independent ofμ, the dynamical critical couplingJ d crit depends on the particular value ofμ.

V. MANY-BODY PHYSICS IN THE MECHANICAL REGULATION OF ADHESION
Our results tie the effective bending rigidity, κ, and in turn interactions between neighboring adhesion bonds, J ∝ κ −1/2 (see Appendix A), to the collective phase behavior of adhesion clusters at equilibrium, and to distinct dynamical phases of cluster dissolution and formation. Based on the quantitative relationship between the coupling strengthJ and bending ridigity κ given by Eq. (A3), and an order-of-magnitude estimation of the relevant parameters listed in Table I we find that the coupling strength in cellular systems lies within the range 0 J 2.5. Notably, both the statical and dynamical critical point at moderate values of the external force values and/or intrinsic binding-affinity lie within said range (see Fig. 5 and 9). Yet, it remains to be explained why a near-critical coupling may be beneficial for cells, and how it may be regulated.
Our results provoke the hypothesis that the membrane rigidity (and hence the coupling strength) may lie close to the statical critical value for quasi-static, and near the dynamical critical value for transient processes. Mechanical regulation of the bending rigidity can be achieved through hypotonic swelling [91], (de)polymerization of the F-actin network [92,93], by decoupling the Factin network from the plasma membrane [55], through changes of the membrane composition [88,89,94,95] or integral membrane proteins [95], membrane-protein activity [96], temperature modulation [28,35,88], and acidosis [45], to name but a few. Moreover, it has been shown experimentally that temperature modulations affects adhesion strength through changes in membrane fluidity [28], cell elasticity [35], or via a temperature cooperative process [26], albeit the denaturation of the binding proteins also provides a possible explanation [36].
Below we argue that the change in the response of a cell to a perturbation, defined as a change in the equilibrium binding strength or association/dissociation rates, is largest near criticality. This results in either a very small or very large response, depending on the change of the underlying parameter. Here we follow the same kind of reasoning as rooted in the criticality hypothesis, which states that systems undergoing an order-disorder phase transition achieve the highest trade-off between robustness an flexibility around criticality [97].

A. Criticality at equilibrium
In Fig. 10a we depict how oscillations in the coupling strength (arising through oscillations in the bending rigidity κ) around the statical critical point affect the average fraction of closed bonds. Similar oscillatory patterns and their effect on the adhesion strength have been observed in vascular smooth muscle cells, where changes in the bending rigidity were concerted by the remodeling of the actin cytoskeleton [34,57,58]. Minute changes in the amplitude, δJ, can drive the systems's behavior from oscillations within a dense phase with ϕ(t) > 0.5 to intermittent periods of nearly complete dissolution (compare full and dashed lines in Fig. 10a). Hence we find that the response (i.e. ϕ(t) ) is most sensitive to a change in the amplitude δJ whenJ lies close to the statical critical point.
Similarly, in Fig. 10b we show the response of ϕ(t) to a mechanical perturbation oscillating quasi-statically between a pulling and a pushing force,h(t) =h ref + δh sin(ωt) (for practical examples see e.g. [98,99]). Such mechanical perturbations can for example arise through changes in active stresses generated within the cytoskeleton [100]. Here as well, a small change in the force δh acting on the cluster, can lead to stark differences in the cluster stability ϕ(t) . The sensitivity to a change in the the force is most amplified near the statical critical cou-plingJ s crit (compare full and dashed lines in Fig. 10b), where a small change in the amplitude, δh, can cause intermittent periods of essentially complete cluster detachment.
Drastic changes in the average number of closed bonds have been observed experimentally in adhesion frequency assays and single-molecule microscopy [21,68]. There it was shown that binding affinities and binding dynamics for a T-cell receptor (TCR) interacting with the peptidemajor histocompatbility complex (pMHC) are more than an order of magnitude smaller in solution (i.e. in 3D) as compared to when they are anchored to a cell membrane (i.e. in 2D). One possible contribution to the discrepancy between the 3D and 2D binding kinetics is the difference in the reduction of the entropy upon binding, which is larger in 3D than in 2D [40]. However, it has been explicitly remarked that this contribution alone does not explain the measured difference in the binding affinities [40]. The authors of Ref. [21,68] rationalize these differences in binding in terms of a cooperativity between neighboring TCRs due to the anchoring membrane. In particular, Fig. 11a in the Supplementary Material of Ref. [68] shows the adhesion frequency P a (t c ) ∈ [0, 1], defined as the fraction of observed adhesion events between the TCR and pMHC as a function of the contact time t c between the anchoring membranes, derived from Monte-Carlo simulations. Upon introducing a heuristic neighbor-dependent amplification factor in the binding rates the authors observe an amplification of the adhesion frequency P a (compare squares with diamonds), indicating an increase in binding events in agreement with their experimental observations.
We may relate our results to the observations in Ref. [68] by recalling the relation between P a and N c = N ϕ , i.e. P ss a ≡ lim tc→∞ P a (t c ) = 1 − exp (−N ϕ ) (see [101] as well as Eqs. (1) and (2) in [68]). In our model the aforementioned amplification factor arises naturally from a nonzero coupling strengthJ due to the anchoring membrane. Indeed, in Fig. 3e an increase inJ leads to an increase in ϕ , which in turn causes an increase of the steady state adhesion frequency P ss a . Hence we find that the amplification factor in [68] and couplingJ in our model have the same effect on the adhesion frequency.
A similar observation was made in [39] on the basis of a detailed analysis of the binding affinities of the adhesion receptor CD16b placed in three distinct environments: red blood cells (RBCs), detached Chinese hamster ovary (CHO) cells, and K562 cells. Based on Fig. 4a,b in [39] the adhesion frequency for RBCs is around a 15fold larger than for CHO and K562 cells. In the discussion the authors point towards the modulation of surface smoothness as an explanation for the observed differences in adhesion frequency [39]. Since K562 cells are known to have a larger bending rigidity than RBCs [102, 103] (we were unfortunately not able to find the corresponding information for CHO cells in the existing literature), it is expected that the coupling strengthJ is generally higher in the latter (see Appendix A), which provides a potential explanation for the observed difference in adhesion frequencies between RBCs and K562 cells.

B. Criticality in kinetics
Many biological processes [104][105][106][107] and experiments [108][109][110] involve adhesion under transient, nonequilibrium conditions, where cells can become detached completely from a substrate (for a particular realization with a constant force see [108]). The duration of these transients may be quantified by the mean dissolution and formation time, t d,f (see Fig. 8). Imagine that the cell can change the bending rigidity by an amount ∆κ that in turn translates into a change in coupling, J =J + ∆J ∝ 1/ √ κ + ∆κ. If the mechanical regulation is to be efficient, a small change of ∆J should effect a large change of t d,f .
The efficiency of the regulation, expressed as the change of mean dissolution/formation time in response to a change ∆J, ∆ t d,f ≡ t d,f (J + ∆J) − t d,f (J) , is shown in Fig. 10(c-d). The results demonstrate that the regulation is most efficient, that is gives the largest change, whenJ is poised near the dynamical critical coupling,J J d crit , regardless of the magnitude of the change ∆J. Recall that the formation and dissolution rate, 1/ t f and 1/ t d , respectively, are highest at the dynamical critical coupling (see Fig. 8). Therefore not only do we find the largest response to a change inJ, but also the fastest formation and dissolution kinetics at the dynamical critical couplingJ J d crit . An example where fast kinetic (un)binding and a large sensitivity to the bending rigidity can be beneficial is found in tumor cells that undergo metastasis -the process through which tumor cells spread to secondary locations in the host's body. Recent studies suggest that cancer cells are mechanically more compliant than normal, healthy cells [111]. Moreover, experiments with magnetic-tweezers have shown that membrane stiffness of patient tumor cells and cancer cell-lines inversely correlates with their migration and invasion potential [60], and an increase of membrane rigidity alone is sufficient to inhibit invasiveness of cancer cells [90]. Cells with the highest invasive capacity were found to be five times less stiff than cells with the lowest migration and inva-sion potential, but the underlying mechanism behind this correlation remained elusive [60].
Based on our results a decrease in the bending rigidity, and hence the membrane stiffness, can alter both, the equilibrium strength of adhesion (see Fig. 3) as well as the kinetics of formation and dissolution of adhesion domains (see Fig. 8). This may provide a clue about the mechanical dysregulation of cell adhesion in metastasis in terms of a softening of the cell membrane.

VI. CRITICALITY IN THE ISING MODEL
By settingh = 0 and writing δ σj ,−1 = (1 − σ j )/2 we find that Eq. (1) is, up to a constant, identical to the Hamiltonian of the isotropic ferromagnetic Ising model in a uniform external magnetic field M ≡ µ/2. Therefore, our findings, and in particular the uncovered dynamical phase transition, also provide new insight into equilibrium and kinetic properties of the Ising model in the presence of a uniform external magnetic field.
The equilibrium properties of the two-dimensional Ising model in the absence of a magnetic field, such as the total free energy per spin, statical critical point, and binodal line were obtained in the seminal work by Onsager [83]. The effect of a uniform magnetic field has mostly been studied numerically [112,113], e.g. by Monte-Carlo simulations [114] and renormalization group theory [115], but hitherto no exact closed-form expression for the free energy per spin has been found. On the Bethe-Guggenheim level the free energy density, binodal line, spinodal line, and statical critical point were known [116], but to our knowledge we are the first to provide an exact closed-form expression for the equation of state in the presence of a uniform magnetic field (see Appendix D).
The kinetics of the two-dimensional Ising model have been studied in the context of magnetization-reversal times (i.e. the time required to reverse the magnetization) [117][118][119], nucleation times [120,121], and critical slowing down [122,123]. Here we report a new type of crit dynamical critical phenomenon related to a first-order discontinuity and a global minimum of the magnetization reversal time at the concurrent dynamical critical point (see Fig. 8), which is fundamentally different from the statical critical point. The dynamical phase transition reflects a qualitative change in the instanton path towards magnetization reversal, and has not been reported before.
In Table II we summarize the values of the statical and dynamical critical points obtained by the mean field and Bethe-Guggenheim approximation in the absence of a magnetic field, and for a general coordination numberz (for a derivation of the dynamical critical points see Appendices E 8 and E 9). We also state the exact statical critical point of the two-dimensional Ising model. Conversely, the exact dynamical critical point of the twodimensional Ising model remains unknown as it requires the exact free energy density as a function of the fraction of down spins (see Eq. (12) for the result within the Bethe-Guggenheim approximation). A lower bound on the dynamical critical point is set by the statical critical point, as the latter denotes the onset of an interior local maximum that is required for the dynamical critical point (see Fig. 8). The exact dynamical critical point may provide further insight into the nature of the dynamical phase transition. Moreover, it also sets a lower bound on the magnetization reversal times per spin in ferromagnetic systems in the absence of an external force.

VII. CONCLUDING REMARKS
The behavior of individual [11] and non-interacting [62][63][64]74] adhesion bonds under force, the effect of the elastic properties of the substrate and pre-stresses in the membrane [43,75], as well as the physical origin of the interaction between opening and closing of individual adhesion bonds due to the coupling with the fluctuating cell membrane [15-19, 44, 46, 124, 125] are by now theoretically well established. However, in order to understand the importance of these interactions and their manifestation for the mechanical regulation of cell adhesion in and out of equilibrium one must go deeper, and disentangle the response of adhesion clusters of all sizes to external forces and how it becomes altered by changes in membrane stiffness. This is paramount because interactions strongly change the physical behavior of adhesion clusters under force both, qualitatively as well as quantitatively.
Founded on firm background knowledge [11, 15-19, 43, 44, 46, 62-64, 74, 75, 124, 125] our explicit analytical results provide deeper insight into cooperative effects in cell-adhesion dynamics and integrate them into a comprehensive physical picture of cell adhesion under force. We considered the full range of CAM binding-affinities and forces, and established the phase behavior of twodimensional adhesion clusters at equilibrium as well as the kinetics of their formation and dissolution.
We have obtained, to the best of our knowledge, the first theoretical results on equilibrium behavior and dynamic stability of adhesion clusters in the thermodynamic limit beyond the mean field-level (existing studies, even those addressing non-interacting adhesion bonds [63,64,74], are limited to small clusters sizes [16-18, 43, 75]). We explained the complete thermodynamic phase behavior, including the co-existence of dense and dilute adhesion domains, and characterized in detail the corresponding critical behavior.
We demonstrated conclusively the existence of a seemingly new kind of dynamical phase transition in the kinetics of adhesion cluster formation and dissolution, which arises due to the interactions between the bonds and occurs at a critical couplingJ d crit , whose value depends on the external forceh and binding-affinityμ. At the dynamical critical couplingJ d crit , and in turn critical bending rigidity κ d crit ∝ (J d crit ) −2 , the dominant formation and dissolution pathways change qualitatively. BelowJ d crit the rate-determining step for cluster formation and dissolution is the surmounting of the (mostly) entropic barrier to completely bound and unbound states, respectively. Conversely, aboveJ d crit the thermodynamic phase transition between the dense and dilute phase for dissolution, and between the dilute and dense phase for cluster formation, becomes rate-limiting, whereas the completely bound and unbound states, respectively, are thereupon reached by typical density fluctuations.
We expect the non-monotonic dependence of the mean first passage time to cluster dissolution/formation on the coupling strengthJ that is asymmetric around the minimum to be experimentally observable (though the notion of a fully bound state during cluster formation may be experimentally ambiguous). According to our theory the existence of such a minimum and its asymmetric shape would immediately imply a dynamical phase transition in the thermodynamic limit.
Measuring the mean dissolution/formation time (in the absence or presence of an external force) for an ensemble of cells adhering to a stiff substrate seems to be experimentally feasible. The effective membrane rigidity (and thus the couplingJ) could in principle be controlled by varying the membrane composition (e.g. increasing the cholesterol concentration that in turn increases membrane rigidity [88]), by tuning the osmotic pressure of the medium [126], or by the depolymerization of F-actin [127]. Testing for signatures of the theoretically predicted dynamical phase transition thus seems to be experimentally (at least conceptually) possible, and we hope that our results will motivate such investigations.
We discussed the biological implications of our results in the context of mechanical regulation of the bending rigidity around criticality. Based on our results we have suggested that the response of a cell to a change in the bending rigidity may be largest near the statical critical point for quasi-static processes, and near the dynamical critical point for transient processes. This observation agrees with the criticality hypothesis, and might expand the list of biological processes hypothesized to be poised at criticality [128].
Finally, we discussed the implications of our result for the two-dimensional Ising model. The observed dynamical phase transition is related to a first-order discontinuity in the magnetization reversal time, and the exact dynamical critical point for the two-dimensional Ising model remains elusive (see Table II).
We now remark on the limitations of our results. The mapping onto a lattice gas/Ising model (i.e. Eq. (1) and Appendix A; see also [16,17]) may not apply to genuinely floppy membranes encountered in biomimetic vesicular systems [48,65]. Moreover, since we only allow for two possible states of the bonds, i.e. associated and dissociated, we neglect any internal degrees of freedom (e.g. orientations of the bonds) which may contribute to the entropy loss upon binding [40], thereby changing the free energy.
Likewise, the assumption of an equally shared force is generally good for stiff membranes (stiffened by the presence of, or anchoring to, the stiff actin cytoskeleton [9]) or stiff membrane/substrate pairs, flexible individual bonds, low bond-densities, or the presence of pre-stresses exerted by the actin cytoskeleton [43,74,75]. In Appendix E 2 we provide an analysis of the effect of a non-uniform force load. Based on this analysis we find that in the case of rather floppy membranes, corresponding to large values of the coupling strengthJ = O(1), the difference between a uniform and a non-uniform force load is negligible for a broad range of realizations of the non-uniform force distribution. Only under the extreme, non-physiological condition that the ratio of forces experienced by inner and outer bonds is larger than an order of magnitude, we observe significant differences. Therefore, the dependencies of the statical and dynamical critical points on the external force (see Figs. 5 and 8, respectively) are expected to remain valid for a non-uniform force distribution over a large range of force magnitudes.
In their present form our results may not apply to conditions when cells actively contract in response to a mechanical force on a timescale comparable to cluster assembly or dissolution [129], as well as situations in which cells actively counteract the effect of an external pulling force and make adhesion clusters grow (see results in terms of a change in membrane stiffness as well. Finally, throughout we have considered clusters consisting of so-called "slip-bonds", whereas cell adhesion may also involve "catch-bonds" that dissociate slower in the presence of sufficiently large pulling forces [130]. The reason lies in a second, alternative dissociation pathway that becomes dominant at large pulling forces [131][132][133][134]. Our results therefore do not apply to focal adhesions composed of catch-bonds and would require a generalization of the Hamiltonian (1-2) and rate (3). These open questions are beyond the scope of the present work and will be addressed in forthcoming publications.

VIII. DATA AVAILABILITY
The open source code for the evaluation of the equation of state and mean first passage times to cluster dissolution and formation for finite-size systems is available at [135].

IX. ACKNOWLEDGMENTS
We thank the anonymous referees for their valuable comments. The financial support from the "Deutsche Forschungsgemeinschaft" (DFG) through the Emmy Noether Program "GO 2762/1-1" (to AG), and an IM-PRS fellowship of the Max-Planck-Society (to KB) are gratefully acknowledged.

Appendix A: Relation between membrane rigidity and coupling strength
Here we provide a quantitative relation between the effective bending rigidity κ and the coupling strengthJ based on the Results of Ref. [17]. Consider a set of adhesion bonds at fixed positions {r i } coupled to a fluctuating membrane. The effective bending rigidity quantifies the amount of energy needed to change the membrane curvature, and is supposed to depend on the membrane composition [88,89], state of the actin network [55], and other intrinsic factors that determine the mechanical stiffness of the cell. Let {b i } describe the state of all bonds, where b i = 1 denotes a closed and b i = 0 an open bond. The bonds are represented by springs with constant k, resting length l 0 , and binding energy b . Non-specific interactions between the membrane and the opposing substrate are described by a harmonic potential with strength γ, which arises from a Taylor expansion around the optimal interaction distance h 0 between the membrane and the substrate. Assuming a timescale separation between the opening/closing of individual bonds and membrane fluctuations, the following partition function for the state of bonds {b i } can be derived [17] whereμ plays the role of an intrinsic binding-affinity and J ij is an effective interaction between the bonds given bỹ with m (r) = − 4 π kei 0 r(κ/γ) 1/4 , and kei 0 (x) is a Kelvin function defined as kei 0 (x) ≡ ImK 0 (xe 3πi/4 ) where K 0 (z) is the zero-order modified Bessel function of the second kind [17]. A systematic comparison of the equation of state ϕ of the full/explicit model (i.e. reversible adhesion bonds coupled to a dynamic, fluctuating membrane) and of the lattice gas governed by Eq. (A2) has been carried out in [17]. Using the following set of parameter values βκ = 80, βγ = 10 −5 nm −4 , βk = 2.25 × 10 −2 nm −2 , h 0 − l 0 = 45.9 − 50.3 nm, and m (1.5) = 0.42194(6), corresponding to a coupling strength ofJ ≈ 1.0 − 1.2, the authors found a quantitative agreement between the full and lattice gas models (see Fig. 5 in Ref. [17]). Note that the lattice gas model becomes exact in the limit κ → ∞ corresponding toJ = 0.
Our effective Hamiltonian, given by Eq. (1), is directly derived from Eq. (A1) by considering the following arguments: First we note that the effective interactioñ J ij decays exponentially fast as a function of the lattice distance between the bonds, and therefore it suffices to only take into account nearest neighbor interactions [17]. Moreover, since we place the adhesion bonds on a lattice with equidistant vertices, the position dependence in Eq. (A2) drops out and we get |r i − r j | = ∆r. Finally, upon introducing the variables σ i ≡ 1−2b i ∈ [−1, 1], and applying the transformationsμ →μ −zJ andJ → 4J, we arrive at our effective Hamiltonian Eq. (1).
Here we find the relationJ ∝ 1/ √ κ, as mentioned in the main text. (B3) where we have used Eqs. (B2). The core idea is to approximate Ψz k zN (N oc ) by the variational, Bethe-Guggenheim approximation [136][137][138] (for an excellent explanation of the method see [79,116,139]). ForJ = 0 the degeneracy factor must obey To implement this constraint it is convenient to normalize Eq. (B3) at zero coupling and write (B5) where we will determine S 1,2 variationally. Then we have forJ = 0 that Z k = N k . Let us now consider placing pairs of adhesion bonds randomly onto the lattice. The total number of unique lattice configurations for fixed N oo , N oc and N cc may be approximated by Notice that for a two-dimensional square lattice that is either infinite or has periodic boundary conditions, the number of open-closed pairs is always even, and therefore the term (N oc /2)! is well-defined. For a finite two-dimensional square lattice with free boundary conditions the number of open-closed adhesion pairs can be odd, which forces us to consider the generalized factorial (i.e. Gamma function). Replacing the factorial with the Gamma function we get (B7) where Γ (n) = (n − 1)! for n ∈ Z + . Substituting Eqs. (B2) for N cc and N oo leaves N oc as the only free parameter in Eq. (B7).
We approximate S 1,2 by an analytic continuation of the maximum term method [140,141] to real numbers. First we analytical continue the summands over positive real numbers using Eq. (B7). We now approximate both sums in Eq. (B5) by their respective largest term, i.e. by the solutions of the pair of optimization problems which yields . We used Stirling's approximation for the Gamma function to find the local maxima, i.e. Γ (z) = 2π/z(z/e) z [1 + O(1/z)], and therefore expect the accuracy of Eq. (B9) to increase with increasing N . Using Eq. (B9) and introducing ψz k zN (x) = Γ(1 +zN/2)/Ψz k zN (x) we obtain Eq. (8) in the main text. The mean field solution is in turn recovered by settingX k = X * k , which happens automatically forJ = 0 or k = 0 ∨ N . In Fig. B1 we compare the accuracy of the Bethe-Guggenheim (a) and mean field (b) approximations by means of the relative error BG,MF between the exact and approximate free energy density for a twodimensional lattice with ϕ = 1/2 closed bonds as a func-tion of the system size N . We find that the Bethe-Guggenheim approximation converges to the exact free energy density with increasing N , regardless of the coupling strengthJ. Conversely, the mean field approximation diverges with both, increasing system size N and coupling strengthJ.

Phase diagram
a. Independence of the binodal and spinodal line onμ Let us writef (ϕ) =μϕ +g (ϕ), whereg (ϕ) is the remainder of the free energy density after leaving out all linear terms in ϕ. Plugging this expression into the binodal line equations given by Eq. (14) gives Clearlyμ cancels and therefore does not affect the binodal line. The spinodal line -also known as the stability boundary -denotes the boundary between the metastable and unstable state, and is given byf (ϕ) = 0. The spinodal line is as well not affected byμ, since the second derivative of the linear term vanishes.

c. Spinodal line
To determine the spinodal line we calculate the second derivative of the Bethe-Guggenheim free energy densitỹ Eq. (B15) contains Ω ϕ , which we want to express in terms of Ω ϕ . Therefore we use Eq. (B11) and differentiate both sides with respect to ϕ for fixedJ, yielding Solving Eq. (B17) for Ω ϕ yields Φ(ϕ,h) in Eq. (16), which must be non-negative thus implying that Eq. (16) is valid for
For non-zero force we solve Eq. (B20) by means of a "quadratic" Newton series as explained in more detail in Appendix D. The main result for the statical critical fraction obtained by the quadratic Newton series reads ϕ s crit,BG ≈ where g (ϕ) = ∂ ϕ g(ϕ), g (ϕ) = ∂ 2 ϕ g(ϕ), and the auxiliary functions δz(h) and νz(h) are defined as respectively. The statical critical couplingJ s crit is obtained by inserting Eq. (B22) into Eq. (16), and the result is depicted by the gradient line Fig. 5 in the main text, where the black symbols represent the fully converged Newton's series (D17), as well as in Fig. C1 where the gradient line depicts the fully converged Newton's series.
Appendix C: Mean field approximation

Partition function
Within the mean field approximation the partition function reads with X * k ≡ k (N − k) /N , such that the corresponding free energy density in the thermodynamic limit attains the form

Phase diagram
We now evaluate the binodal and spinodal line, as well as the statical critical point within the mean field approximation. The corresponding exact solutions for the Bethe-Guggenheim approximation are given by Eqs. (15)-(19).

a. Binodal line for zero force
In the absence of an external field,h = 0, the binodal line is given bỹ This solution is well known and has been reported in the literature extensively [79,116,139,143]. Forh = 0 we solve for the binodal line numerically.

b. Spinodal line
The spinodal line is in turn given by the solution of f MF (ϕ) = 0 and reads J s,MF (ϕ) = 1 4z The statical critical point is given by the solution of f MF (ϕ) = 0, which after introducing the parameter αh ≡ 12h translates into solving the algebraic equation Notice thatz, the average coordination number, does not enter Eq. (C5). Whenh = 0 the solution is ϕ s crit,MF = 1/2, and the corresponding statical critical point is given byJ s crit,MF =z −1 . To solve Eq. (C5) for non-zero force we first note that 0 ≤ ϕ s crit,MF ≤ 1, and therefore we can divide Eq. (C5) by (1 + ϕ) 4 . Upon introducing the variable w = 2ϕ − 1 we get the equation with f (w) = (w + 1) 2 (w − 1) 2 (w + 3) −4 . Now we recall the Lagrange inversion theorem: Let f (w) be analytic in some neighborhood of the point w = 0 (of the complex plane) with f (0) = 0 and let it satisfy the equation Then ∃a, b ∈ R + such that for |α| < a Eq. (C7) has only a single solution in the domain |w| < b. According to the Lagrange-Bürmann formula this unique solution is an analytical function of α given by Notice that Eq. (C6) has the form of Eq. (C7), and therefore we can use Eq. (C8) to obtain ϕ s crit,MF . To evaluate the derivative inside Eq. (C8) it is convenient to write f (w) k = g(w)h(w), with g(w) ≡ (w + 1) 2k (w − 1) 2k = (w 2 − 1) 2k and h(w) ≡ (w + 3) −4k . Using the fact that Plugging Eq. (C12) into Eq. (C8) and using the relation ϕ = (1+w)/2 we find ϕ s crit,MF = ϕ s,0 (C13) The statical critical coupling for non-zero force, J s crit,MF (h), is obtained by plugging ϕ s crit,MF into Eq. (C4).
In Fig. C1 we plot the statical critical point for the Bethe-Guggenheim and mean field approximation as a function of the forceh. For large pulling force (h ≥ 2) the statical critical point is pushed towards lower values ofJ, and as a consequence the Bethe-Guggenheim and mean field solution start to coincide. For other values ofh, however, the Bethe-Guggenheim and mean field solutions disagree strongly, in particular for weak forces |h| → 0.
To compare the critical points for weak pulling/pushing forces in more detail we inspect the perturbation series given in Eqs. (18) and (19) and compare it with the first two terms of Eq. (C13).
Interestingly, whereas the Bethe-Guggenheim critical point depends onz (see Eq. (19)), the mean field result in Eq. (C15) does not. In the limitz → ∞ we find that the Bethe-Guggenheim statical critical point converges to the mean field solution, which is to be expected.

Appendix D: Equation of state in the thermodynamic limit
Here we derive the equation of state in the thermodynamic limit using the saddle-point technique, i.e.
and ϕ 0 1 , ϕ 0 2 , ..., ϕ 0 M denote the locations of the local minima of the Bethe-Guggenheim free energy densityf(ϕ) ≡ f BG (ϕ) in Eq. (12). The idea behind Eq. (D1) is that in the large N limit we expect the integral over ϕ to be dominated by the immediate neighborhood of the local minima off (ϕ). We may therefore approximate the exponent by its Taylor expansion around these extremal points. In general special care has to be taken when one of the global minima lies at the boundary of the integration interval [144], which turns out not to be the case here.
The locations of the local minima, maxima, and saddle points of the Bethe-Guggenheim free energy density are given by the solution off BG (ϕ) = 0. Notice that here we do not set the intrinsic binding-affinityμ to zero, since we are interested in the stationary points and not the binodal line. We solve the former equation for Ω ϕ and substitute the solution into Eq. (B11) which gives where χ ϕ ≡ ϕ/(1 − ϕ), α ≡ (z − 1)/z, and Rewriting Eq. (D3) gives For a two-dimensional square lattice in the thermodynamic limit we havez = 4 and so α = 3/4. Upon introducing the auxiliary variable ξ ≡ χ 1/4 ϕ we obtain the transcendental equation To solve for the roots of g(ξ) we first consider the forcefree scenarioh = 0 and afterwards solve for the general case.

Zero force
For zero force g(ξ) reduces to a quartic in ξ. The roots of a general quartic equation are known and are given by where Y ± = e 2J±μ/4 and S = 1 2 To determine the local minima we must analyze the properties of these roots starting with the sign of the discriminant, given by For ∆ < 0 there are two distinct real roots and two complex conjugate roots, whereas for ∆ > 0 there are either four real roots or four imaginary roots, where the former scenario applies here. The discriminant is zero at the critical coupling valuẽ Increasing the coupling strength aboveJ ∆=0 gives rise to a local maximum in the free energy landscape, its position being ξ 4 4 /(1 + ξ 4 4 ) (see purple line in Fig. D1).

a. Zero intrinsic binding-affinity
For zero intrinsic binding-affinity we note thatJ ∆=0 coincides with the statical critical pointJ s crit = ln (2)/2. In this limit the four solutions (D7a) and (D7b) simplify substantially, and the corresponding locations of the local minima -which are also global minima -are given by (see Fig. D1a) (D14) ForJ ≤J s crit there is a single global minimum in the free energy landscape, and therefore the weight c 1 given by Eq. (D2) becomes unity. ForJ >J s crit the global minima are two-fold degenerate and located equidistantly from the local maximum at ϕ = 1/2. Since both global minima have the same curvature, we find that the weights are given by c 1,2 = 1/2. Combining these results we obtain the average fraction of closed bonds in the thermodynamic limit for zero intrinsic binding-affinity in the form which was to be expected since the couplingJ does not favor bonds to be open nor closed.

b. Non-zero intrinsic bining-affinity
For non-zero intrinsic binding-affinity,μ = 0, the free energy landscape is tilted, resulting in a unique global minimum with corresponding weight c 1 = 1. As a result, the average fraction of closed bonds, which is dominated by this minimum ϕ 0 1 , is given by (see Fig. D1(b) and (c)) where ln s 0 ≡ ln 2 √ 2e −μ/4 denotes the coupling strength that solves for the root of S in Eq. (D8).

Non-Zero force
We determine the roots of g(ξ) for a non-zero force by means of a convergent Newton series yielding the exact result [80,82,142] where ξ 0 is an initial guess in a convex neighborhood around ξ, g (i) (ξ 0 ) denotes the ith derivative of g(ξ) at the point ξ 0 , and A k (ξ 0 ) are almost triangular matrices of size (k − 1) × (k − 1) with elements , (D18) where θ (x) denotes the Heaviside step function, i.e. θ (x) = 1 if x ≥ 0 and 0 otherwise, and we symbolically set det A 1 = 1. The determinant of almost triangular matrices, also known as upper/lower Hessenberg matrices, can be efficiently calculated using a recursion formula [145], for which a numerical implementation can be found in [146]. If we set g (3) = g (4) = ... = 0 in the almost triangular matrices, the resulting matrixÃ k becomes triangular, implying that its determinant is simply given by the product of its diagonal elements. Making the substitution A k →Ã k in Eq. (D17) yields the socalled "quadratic approximation" [82,142] (D19) that becomes exceedingly accurate when the root moves close to ξ 0 . For the initial point ξ 0 we use the ansatz which is derived by considering an adapted form of Eqs. (D7a) and (D7b) in combination with the implementation of the force term. The weight 2/9 is derived from the term (1+ξ 4 ) 2 /2(1+2ξ 4 ) 2 in Eq. (D4) evaluated at the point ξ = 1 (corresponding to ϕ = 1/2), and the weight 6/5 in front ofJ was selected empirically. This choice assures that Eq. (D6) satisfies the Lipschitz condition between ξ 0 and the root ξ and thus assures the convergence of the Newton's series. Plugging Eq. (D20) into Eq. (D19), and using the relation ϕ 0 1 = ξ 4 /(1+ξ 4 ), we obtain the location of the global minimum -and thus ϕ TD -for non-zero force. Notably, the ansatz given by Eq. (D20) also provides a numerically correct solution for a zero force and non-zero intrinsic binding-affinity. For completeness we write down explicitly all the terms which are used to evaluate Eq. (D19) (higher order terms entering the fully converged series in Eq. (D17) are omitted as they are lengthy).
Let g (ξ) be given by Eq. (D6). Introducing the auxiliary functions the first and second derivative can be written as where cμ ,h (ξ) is defined in Eq. (D4). Eqs. (D15), (D16), and (D19) form our main result for the equation of state in the thermodynamic limit. In Fig. 4 we show the results for various values of the force and intrinsic affinity.
Appendix E: Kinetics of cluster formation and dissolution

Exact algebraic result for small clusters
It is well known that the transition matrix for an absorbing discrete-time Markov chain with a set of recurrent states has the canonical form [78] where 1 is the identity matrix, T d,f is the submatrix of transient states in dissolution/formation, and R the submatrix of recurrent states. In the particular case of cluster dissolution the (2 N − 1) × (2 N − 1) matrix T d entering Eq. (E1) is obtained by removing the last column and row, and the (2 N − 1) × (2 N − 1) matrix T f entering Eq. (E1) by removing the first column and row. If we introduce the column vectorê k with components (ê k ) i = δ ki and the column vector e whose elements are all equal to 1, the mean first passage times for cluster formation and dissolution read exactly In applying Eq. (E2) one must invert a (2 N −1)×(2 N −1) sparse matrix and afterwards sum over 2 N − 1 terms, which is feasible for N 5×5. For a system of N = 4×5 the exact results are shown with the blue line in Fig. 7. Larger clusters are treated within the local equilibrium approximation.

Finite-size results for a non-uniform force distribution
Under the condition of a small combined elastic modulus, corresponding to large values of the coupling strength J 1, the assumption of an equally shared force load is no longer valid [43,74,75]. We therefore address how a non-uniform force distribution affects the equation of state and mean first passage time to cluster dissolution/formation for finite system sizes. Based on Eq. (7) in [75] and Eq. (4) in [147] we introduce a non-uniform force load by making the substitution h → is a normalization constant such that initially, i.e. when all bonds are closed, the total force load is h. The load on bond i, denoted as h i , is given by where ξ ≥ 1, and¯ i ≡ ( i − r)/(d − r) ∈ [0, 1] is a normalized distance of bond i to the center of the lattice, with i defined as the eccentricity of node i, which is the maximum number of edges between node i and any other node in the lattice. The radius r ≡ min i and diameter d ≡ max i of the lattice are defined as the minimum and maximum eccentricity, respectively. With the force distribution given by Eq. (E3), which is depicted in Fig. E1a, closed bonds located at the outer edge of the lattice ( i = 1) experience a larger external force than bonds located at the inner part of the lattice ( i = 0). The parameter ξ is an indicator for the spread in force load among the individual bonds. For ξ = 1, which holds when limJ → ∞ [75], the force distribution at the edge of the cluster is singular and nonphysical. On the contrary, for lim ξ → ∞, which is valid for limJ → 0, we recover the uniform force distribution. In Fig. E1(b-d) we depict the equation of state (b) and mean first passage time to cluster dissolution (c) and formation (d) for mixed Glauber-Kawasaki dynamics with a constant Glauber attempt probability p k → p = 0.5 and for various values of ξ under a pulling or pushing force (h = ±0.5). The results were obtained by exact summation/algebraic techniques. Interestingly, for ξ ≥ 1.1 the equation of state and mean first passage times are almost identical to the uniform force load solutions that correspond to ξ → ∞. Only for ξ < 1.1, which is valid for very large coupling values corresponding to extremely floppy membranes, we observe deviations from the uniform force results. The origin of the deviations is the extreme force load on the outer bonds, which is ξ/(ξ − 1) times larger than the force load on the inner bond. For ξ = 1.01 this leads approximately to a factor of ×10, and for ξ = 1.001 this leads approximately to a factor of ×32. Hence, for most physically meaningful realizations of a non-uniform force distribution (i.e. distributions based on Eq. (E3)) the results converge to the uniform force solutions. Only under the extreme conditions where the force load on the outer bonds becomes at least an order of magnitude larger compared to the inner bonds we find large deviations from the uniform force load.
Note that the relative fraction of edge bonds in the limit of larger system sizes (and specifically in the thermodynamic limit) vanishes. Therefore we expect a nonuniform force load, which mainly penalizes the edge bonds for ξ → 1, to have an even weaker effect on the equation of state and mean first passage times in large systems.

Proof of detailed-balance for local equilibrium rates
Before stating the explicit result for the mean first passage time to dissolution/formation in the local equilibrium approximation, we prove that the local equilibrium transition ratesw k→k±1 given by Eq. (21) obey detailed-balance w.r.t Q k defined in Eq. (6). The effective transition rates are obtained by mapping the full mixed Glauber-Kawasaki dynamics onto an effective birth-death process over the number of closed bonds (see Fig. 6), where we assume that the dynamics reaches a local equilibrium at any number of closed bonds before any transition. As a result, the birth-death process is a Markov chain on the free energy landscape for the fraction of closed bonds ϕ. Recall that the original Glauber rates in Eq. (3) obey detailed-balance w.r.t. the Hamiltonian H({σ j }), and therefore where we have explicitly incorporated the constraints arising from single-bond-flip dynamics by means of the Kronecker delta s. Upon summing the left hand side of Eq. (E4) over all initial configurations with N c ({σ j }) = k and over all rates that jump to a configuration with N c ({σ j } i ) = k ± 1, we reach all possible final configurations with N c ({σ j } i ) = k ± 1, with a backward rate given by the sum of all rates that jump to a configuration with N c ({σ j }) = k. Hence we find the equality Comparing Eq. (E5) with Eq. (21), we recognize the left and right hand side asQ kwk→k±1 andQ k±1wk±1→k , respectively, which proves the effective detailed-balance re-lationQ kwk→k±1 =Q k±1wk±1→k . (E6)

First passage time statistics within the local equilibrium approximation
The local equilibrium approximation maps the complete mixed Glauber-Kawasaki dynamics onto an effective birth-death process with a right-acting tri-diagonal transition matrix P le of size (N + 1) × (N + 1) with elements To obtain the mean first passage time we use the same algebraic technique as for small clusters. Upon removing the first/last row and column of P le we obtain the submatrix T le d,f for cluster dissolution and formation, respectively. We can invert the tridiagonal submatrix exactly, which leads to the following LU/UL decomposition where A d and B d are the lower and upper triangular matrix with elements and A f and B f are the upper and lower triangular matrix with elements A proof that Eq. (E8) is indeed the inverse of 1 − T le d,f is given in the SM. Let us denote with τ le d,f m the mean first passage time to cluster dissolution and formation, starting from the state with m closed bonds. Using Eq. (E8) we obtain an exact expression for the first moments whereê m is the column vector with dimension N with components (ê m ) i = δ mi , and e is the column vector with all components equal to 1. Notice that Eq. (E8), and therefore Eqs. (E11) and (E12), are applicable to any right-acting tri-diagonal transition matrix with rates obeying detailed-balance. Although we only present the mean first passage time here, we can easily obtain any higher order moments of the first passage time to cluster dissolution and formation using Eq. (E8) [78]. Notably, Eqs. (E11) and (E12) appear to have a similar structure as the largest eigenvalue of the transition matrix in classical nucleation theory [120,148].

A bound on the effective transition rates
Here we present a bound on the local equilibrium rates given by Eq. (21) which proves that the transition rates are strictly sub-exponential in N . First, we consider a bound for the exit rates w ± exit ({σ i }) defined in Eq. (22), which contain a sum over the original Glauber rates that are defined in Eq. (3). Since 1 − tanh (x) ≥ 0 ∀x ∈ R the Glauber rates are non-negative, and therefore the exit rates obey the bound with c + k = N − k and c − k = k denoting the number of terms inside the sum of Eq. (22), and w max k→k±1 denotes the largest transition rate to go from a state with k to k ± 1 closed bonds. The largest transition rate can be written as where ∆H min k→k±1 ≡ inf denotes the smallest possible energy change between two configurations {σ j } and {σ j } i with N c ({σ j }) = k and N c ({σ j } i ) = k ± 1, respectively. To obtain a closed-form expression for w max k→k±1 we first note that the contribution to ∆H min k→k±1 from the external force and intrinsic binding-affinity are fixed and given by the second and third term in Eq. (3). Therefore we are left to consider the smallest energy change due to the coupling strength, which we denote as ∆J H min k→k±1 . For a square lattice with free boundary conditions the minimal energy "forward transitions" with energy difference ∆J H min k→k+1 for various values of k are depicted in Fig. (E2). Similarly, the minimal energy "downward transitions" with energy difference ∆J H min k→k−1 are obtained by interchanging the open (red) and closed (green) adhesion pairs in Fig. (E2). Combining these two results yields ∆J H min k→k±1 = 2m ± kJ , with , (E16) and delivers the expression for w max k→k±1 . Finally, since w max k→k±1 is independent of the specific configuration {σ i } at fixed k, it drops out of the sum over {σ i } in Eq. (21) for the effective transition rates, and therefore the bound in Eq. (E13) can directly be applied to the effective transition rate upon multiplying both sides with the Glauber attempt probability which yields the bound on the effective transition rate. The lower and upper bound for the effective transition rate are used to determine an upper and lower bound for the mean first passage time to cluster dissolution and formation, respectively. The specific result for a rectangular lattice of size N = 6 × 7 for pure Glauber dynamics (i.e. p k = 1 ∀k) is shown in Fig. E3. For small values of the coupling strengthJ we find that the upper bound in Eq. (E17), corresponding to the lower bound in Fig. E3 (since t d,f ∝ 1/w k→k±1 ), is saturated by the exact effective transition rate. Conversely, for large values of the coupling strength it seems that the lower bound in Eq. (E17) is saturated.

Approximate effective transition rate
For systems larger than N ≈ 50 bonds the combinatorics involved in the computation of Q k defined in Eq. (6) andw k→k+1 in Eq. (21) become prohibitive, and thus forces us to make further approximations. To get Eq. (23) fully explicit we make an "instanton" approximation forw k→k±1 using the Bethe-Guggenheim approximation with the bound given by Eq. (E17), and reads with α k = δ k0 + δ kN + 2X k /N ∈ [0, 1], c + k = N − k, c − k = k, andX k given by Eq. (9). The prefactor α k c ± k is a measure for the number of "favorable" adhesion bonds that are most likely to flip in a configuration with k closed bonds. For k = 0∨N all bonds have an equal surrounding in the thermodynamic limit (or for a periodic lattice), and therefore all c ± k open/closed bonds are equally likely to attempt a flip. For 0 < k < N it becomes energetically more favorable to flip a bond which is part of an openclosed adhesion pair (see Fig. E2). To determine the number of bonds that constitute an open-closed pair, we recall thatzX k is a measure for the number of openclosed pairs in a lattice of size N with k closed bonds. Upon dividing by the total number of pairs in the system, given byzN/2, we obtain the probability to select an open-closed pair in the lattice that is given by To prove that the approximate effective rate given by Eq. (E18) obeys the bound given by Eq. (E17) we apply a chain of inequalities. First we note that 0 ≤ 2X k /N ≤ 1/2, where the upper bound follows from consideringJ = 0 and k = N/2 in Eq. (E19), and the lower bound is given for k = 0 ∨ N or the limitJ → ∞. From this it follows that 0 ≤ α k ≤ 1, and therefore 1 ≤ max(1, α k c ± k ) ≤ c ± k . Finally, since we use p k w max k→k±1 in Eq. (E18), it cancels on both sides of the inequality in Eq. (E17), which leaves to prove the inequality we have proven above and thereby completes the proof. Fig. E3 shows the mean first passage time to cluster dissolution and formation obtained with the approximate effective rates (E18) in combination with the Bethe-Guggenheim approximation for Q k for a lattice of size N = 6 × 7 for pure Glauber dynamics (i.e. p k = 1 ∀k). The results obtained with the approximate rates (black symbols) agree to a high degree with the results obtained by the exact effective rates (blue solid line).
, the inequality (E22) is preserved when exponentiated to 1/N . The thermodynamic limit of Eq. (E22) is a scaling limit, i.e. lim s ≡ lim N →∞ | l/N =ϕ l k/N =ϕ k , and thus the inequality (E22) becomes ] 1/N = 1, all four limits in Eq. (E23) exist and thus may be taken separately, implying the convergence of the upper bound in Eq. (E23) to the lower bound. Thereby t d,f becomes squeezed in-between rendering the inequality an equality. Since lim sfN (k † /N ) = f(ϕ max ) and lim sfN (l † d,f /N ) =f(ϕ d,f min ) we finally obtain Eqs. (24) and (25), thus completing the proof. 8. Evaluation of the Bethe-Guggenheim mean first passage time in the thermodynamic limit In the previous section we have proven that in the thermodynamic limit the mean first passage time to cluster dissolution/formation scales as τ d,f t d,f N = e N∆f † , where ∆f † denotes the largest left/right barrier in the Bethe-Guggenheim free energy density Eq. (12). In this section we determine ∆f † and thereby obtain a closedform expression for the mean first passage time per bond in the thermodynamic limit.
a. Case 1:J ≥ 0,h =μ = 0 We first consider the mean first passage time to cluster dissolution/formation in the absence of an external force and intrinsic binding-affinity. Due to the Z 2 symmetry of the coupling strength, we note that t d = t f . Our first task is to find the locations of the global maximum and minimum in the free energy landscape, denoted by and sets the dynamical critical coupling value for the zero-field Ising model under the Bethe-Guggenheim approximation. Surprisingly, for the two-dimensional square lattice withz = 4 we exactly recover the static critical point obtained by Onsager [83]. To check whether this is a mere coincidence we note that for the honeycomb lattice withz = 3 the exact statical critical point is given bỹ J s crit = 1 2 ln 2 +  Figure E4. Master scaling of mean dissolution and formation times per bond for finite clusters and in the thermodynamic limit. t d,f for cluster dissolution (a-b) and formation (c-d) as a function of the couplingJ for a pair of intrinsic affinitiesμ = 0 andμ = 0.5 and various cluster sizes (symbols) as well as the thermodynamic limit (lines) in the presence of an external pushing (a,c) and pulling (b,d) force; Symbols are evaluated with local equilibrium approximation Eqs. (23) using Q BG k (Eqs. (6) and (8)) andw k→k+1 from Eq. (E18) with p k = 1 ∀k (i.e. pure Glauber dynamics) . The discrepancy between the lines and symbols is due to finite-size effects.
givesJ d crit = −1 2 ln 2 1/3 − 1 = 0.67..., and so we find direct evidence that the dynamical critical point in the Bethe-Guggenheim approximation does not (at least not always) coincide with the exact statical critical point.
Combining our results for the locations of the global maximum and minimum we obtain the following result for the mean first passage time per adhesion bond in the thermodynamic limit for the zero field Ising model on a In the strong coupling limit we find from Eq. (E26) limJ →∞ t d,f = 2, which is identical to the result obtained for zero coupling. The physical intuition behind this result comes from considering the average number of steps required to change the state of a single independent adhesion bond. For zero force and intrinsic bindingaffinity the probability to associate/dissociate is a 1/2, and therefore the average dissolution/formation time is given by For an infinite coupling strength the interaction between the bonds is so strong that effectively the system behaves as one "super bond", and therefore the average dissolution/formation time is equal to that of a single independent adhesion bond.
b. Case 2:J ≥ 0,μ = 0,h = 0 Here we use the results obtained in Appendix D which leads to Using a quadratic Newton series (which is defined in Appendix D), Eqs. (E27) and (E28) are directly applicable to the non-zero force scenario upon applying the transformation ϕ i → ϕ * i , where ϕ * i = ξ * 4 i /(1 + ξ * 4 i ) and Figure E5. Mean dissolution/formation time in the thermodynamic limit: Bethe-Guggenheim versus mean field approximation. The Bethe-Guggenheim approximation is given by Eq. (E26), and the mean field approximation is given by Eq. (E35), where we solve for Eq. (E31) numerically.
Our analytical results for the mean first passage time to cluster dissolution and formation per adhesion bond are depicted in Fig. 8 and Fig E4 for zero and nonzero external force respectively; note the remarkable agreement between the black solid line depicting the thermodynamic limit and the results for finite system sizes on the order of N ≥ 10 × 10. Similarly to our previous analysis we must determine the global minimum and maximum, respectively, of the mean field free energy density given by Eq. (C2). For convenience we here consider only zero force and zero binding-affinity.
Below the statical critical couplingJ s crit,MF = 1/z there is a unique global minimum located at ϕ d,f min = 1/2. Above the statical critical coupling there exist two global minima, which are given by the non-zero solutions of the transcendental mean field equation In this Supplemental Material we prove the expression for the inverse of the matrix 1 − T le d,f used in Appendix E4, where T le d,f is the sub-matrix for cluster dissolution/formation, The elements of 1 − T le d,f are given by where The inverse can be expressed as where A d and B d are the lower and upper triangular matrix, respectively, with elements and A f and B f are the upper and lower triangular matrix, respectively, with elements andQ j ≡ Q j /p j where p j is the Glauber attempt probability and Q j the partition function constrained to j, the number of closed bonds defined in Eq. (6). To prove that Eq. (S3) is the inverse of 1 − T d,f we evaluate the corresponding matrix product. We first write the elements of the fundamental matrices as and For cluster dissolution, the matrix product is given by The difference between two consecutive elements of the fundamental matrix can be obtained from Eq. (S6), and reads and similarly