Formation of semi-dilute adhesion domains driven by weak elasticity-mediated interactions.

Cell-cell adhesion is established by specific binding of receptor and ligand proteins anchored in the cell membranes. The adhesion bonds attract each other and often aggregate into large clusters that are central to many biological processes. One possible origin of attractive interactions between adhesion bonds is the elastic response of the membranes to their deformation by the bonds. Here, we analyze these elasticity-mediated interactions using a novel mean-field approach. Our analysis of systems at different densities of bonds, ϕ, reveals that the phase diagram, i.e., the binodal and spinodal lines, exhibit a nearly universal behavior when the temperature T is plotted against the scaled density x = ϕξ(2), where ξ is the linear size of the membrane's region affected by the presence of a single isolated bond. The critical point (ϕc , Tc) is located at very low densities, and slightly below Tc we identify phase coexistence between two low-density phases. Dense adhesion domains are observed only when the height by which the bonds deform the membranes, h0, is much larger than their thermal roughness, Δ, which occurs at very low temperatures T≪Tc. We, thus, conclude that the elasticity-mediated interactions are weak and cannot be regarded as responsible for the formation of dense adhesion domains. The weakness of the elasticity-mediated effect and its relevance to dilute systems only can be attributed to the fact that the membrane's elastic energy saturates in the semi-dilute regime, when the typical spacing between the bonds r≳ξ, i.e., for x≲ 1. Therefore, at higher densities, only the mixing entropy of the bonds (which always favors uniform distributions) is thermodynamically relevant. We discuss the implications of our results for the question of immunological synapse formation, and demonstrate that the elasticity-mediated interactions may be involved in the aggregation of these semi-dilute membrane domains.


I. Introduction
The cellular membrane has the ability to adhere to different biological elements, including the extracellular matrix (ECM), the cytoskeleton and other cells. Cellular adhesion is mediated by several adhesion proteins (e.g., cadherins, integrins, and proteins from the immunoglobulin superfamily) that form specific bonds with receptors embedded in the adhesive element. 1 These adhesion bonds often aggregate into macroscopically large adhesion clusters, such as focal adhesions, adherens junctions and gap junction plaques. [2][3][4] In addition to providing mechanical stability to cells, these adhesion domains are essential for numerous biological processes, including signal transduction, 5 T-cell activation, 6 and tissue formation. 7 Therefore, it is paramount to gain a comprehensive understanding of the biophysical principles that govern the formation of adhesion clusters.
Over the past two decades, many studies have been conducted in order to better understand the biophysical interactions playing a role in the formation of adhesion clusters. 8,9 Special attention has been directed to the effective interactions that are induced by membrane elasticity and thermal undulations. These non-specific interactions have also been studied in relation to condensation of trans-membrane proteins (membrane ''inclusions''), [10][11][12] and in the broader context of ''Casimir-like'' interactions in condensed matter. 13 Specifically to the problem of adhesion domains, membrane mediated interactions between adhesion bonds originate from two interrelated mechanisms operating in concert. The first mechanism is related to the suppression of membrane thermal fluctuations by the adhesion bonds, which locally fix the membrane's height. 14 The resulting loss in the membrane's fluctuation entropy can be partially mitigated if the adhesion bonds aggregate into a single domain. The second mechanism stems from local membrane deformations imposed by the pinning points, which can trigger a redistribution of the adhesion bonds in order to minimize the elastic curvature energy. 15 Thus, membrane elasticity and thermal fluctuations induce a potential of mean force (PMF) between the adhesion bonds. The main challenge in deriving expressions for membrane mediated interactions arises from their many-body character, 16 i.e., their non-trivial dependence on the spatial distribution of the adhesion bonds.
Theoretical studies of membrane mediated interactions are often based on Helfrich's elasticity theory. 17 Within the framework of this model, the membrane is considered as a two dimensional sheet fluctuating over a flat adhesive surface. Using the Monge gauge representation and assuming a small membrane curvature, the elastic energy can be expressed by the effective Hamiltonian where k is the membrane's bending modulus, h = h(r) is the membrane's height (relative to an arbitrary reference plain) at position r = (x, y), and the integration is taken over the membrane's projected area. The first term in eqn (1) stands for the bending energy of the membrane, while the second one denotes a non-specific confining potential due to the interactions between the membrane and its surroundings, specifically an underlying adhesive surface. The attachment between the latter and the membrane by N bonds can be incorporated by a set of height constraints satisfying h( and h 0 is the height of the adhesive surface. The free energy corresponding to Hamiltonian (1) under these constraints constitutes the PMF between the adhesion bonds.
A commonly used practice in membrane elasticity studies is to assume that the membrane's free energy has the same form as eqn (1), with a renormalized bending modulus and with V(h) representing an effective potential between the surface and the membrane's mean height profile. 18 Using this approach, Bruinsma, Goulian and Pincus studied the thermodynamics of domains of gap junctions. 19 Two regimes with distinct expressions for V(h) have been proposed, corresponding to different membrane-surface interactions. In the first regime, coined the Helfrich regime, the bending modulus k is small and, therefore, thermal fluctuations of the membranes are significant. The membrane interacts with the surface via thermal collisions, creating an effective repulsive potential 20 The resulting free energy has been analyzed within a mean-field picture assuming a lattice of equally spaced gap junctions, and was found to grow logarithmically with the lattice spacing. This result has received support from coarsegrained membrane simulations. 21 Another prediction of ref. 19 was that due to the fluctuation-induced attraction between the gap junctions, the temperature is renormalized downward. This prediction was later examined in several computational studies, which demonstrated that, indeed, the renormalized temperature is about third to half of the thermodynamic temperature. [22][23][24][25] These findings highlight the important role of thermal fluctuations in facilitating the conditions required for adhesion cluster formation. The fact that the renormalized temperature remains positive implies that in order to achieve aggregation of adhesion bonds, other attractive interactions must also be present.
The second regime examined in ref. 19, termed the van der Waals regime, is characterized by small thermal fluctuations, which allows one to consider a Lennard-Jones type potential between the membrane and the surface. For small deviations from the potential's minimum, a quadratic approximation for VðhÞ ¼ 1 2 gh 2 can be assumed. In contrast to the Helfrich regime where the range of the fluctuation-induced interactions diverges, the elasticity-mediated interactions in the van der Waals regime span over a characteristic healing length beyond which the membrane sets back to the minimum of the confining potential. As in the Helfrich regime described above, the mean-field free energy was calculated in ref. 19 for a lattice distribution of gap junctions. Coupling the effective interactions with the mixing entropy of the adhesion bonds yields the full free energy of the system, from which conditions for the condensation of adhesion bonds have been derived.
In the past few years, attempts to develop a more rigorous statistical mechanical treatment of the van der Waals regime have been made. Considering the elastic energy given by eqn (1) with a harmonic confining potential VðhÞ ¼ 1 2 gh 2 , the partition function of the system is given by where the N pinning points are accounted for through a series of Dirac-delta functions, and the integration is performed over all possible height profiles of the membrane. The partition function Z N can be evaluated by (i) taking the Fourier representations of the height function and the Dirac-delta functions, (ii) applying N Hubbard-Stratonovich transformations, and (iii) evaluating the resulting Gaussian integrals. [25][26][27][28] This leads to the following expression: where Z 0 is the partition function corresponding to Hamiltonian (1), with VðhÞ ¼ 1 2 gh 2 and without (N = 0) adhesion bonds. The coupling matrix M appearing in eqn (4) is given by where the sum runs over all independent Fourier modes q, A p is the projected area of the membrane, kei(x) is the Kelvin function, 29 and denotes the mean square of height fluctuations (thermal roughness) of the membrane in the absence of adhesion bonds.
For a given distribution of adhesion bonds, the PMF is given by the free energy The first term on the r.h.s. of eqn (7) gives the energy of the height function that minimizes Hamiltonian (1) with the harmonic potential, subject to the height constraints imposed by the bonds. The second term is the entropic contribution due to the thermal undulations around this profile. 27 Notice that the energetic and the entropic components in the free energy decouple in this model, which follows from the quadratic nature of the Hamiltonian in q-space. Also notice that both terms in eqn (7) depend on the elements of the matrix M ij (5) in a non-linear manner, which is a mathematical manifestation of the many body nature of the PMF. An interesting observation was made by Speck, Reister and Seifert, who argued that for small thermal roughness (D { h 0 ) the model depicted by eqn (7) belongs to the two dimensional Ising universality class. 28 Furthermore, if the healing length x is smaller than the typical distance between the bonds, the model can be mapped onto a lattice-gas with nearest neighbor interactions. By estimating the effective interaction parameter between adhesion bonds occupying neighboring sites, the authors of ref. 28 were able to draw the phase diagram of the system and estimate the critical temperature below which clusters appear.
Despite the insights gained from previous studies, a satisfactory description of the thermodynamic behavior of the model described by eqn (7) is still lacking. Here, we take another look at this problem and derive a more accurate picture of the phase diagram, for a wide range of healing lengths, x, and adhesion bond densities, f. Our investigation relies on a novel mean-field treatment of the system's free energy. We obtain the spinodal and binodal curves and locate the critical temperature of the system, T c , above which adhesion domains do not form. Results for different systems exhibit data collapse when (D/h 0 ) 2 B T/T c is plotted as a function of the rescaled density x 2 f. Interestingly, we find that the critical point is located at extremely low densities, which is linked to the many-body membrane mediated PMF. Therefore, close to critically, phase coexistence is found between two extremely dilute phases, while dense domains form only for T { T c , i.e., when each bond deforms the membrane considerably.
The paper is organized as follows: in Section II we introduce our mean-field theoretical treatment. This approach involves calculations of the elastic energy of systems with randomly distributed adhesion bonds at various densities. These calculations, which are described in Section II A, yield the expression for the mean-field energy of the system. The free energy is then obtained by combining the energy with the mean-field mixing entropy. In Section II B, we analyze the dependence of the free energy on the density of the bonds, and draw the phase diagram of the system, i.e., the binodal and spinodal lines. We discuss and summarize our findings in Section III.

II. Mean-field theory
The PMF, F N , given by eqn (7) corresponds to a system with a given spatial distribution of N fixed adhesion bonds. The thermodynamics of a system with N mobile bonds is characterized by the free energy F, which depends on the bond density f = aN/A p , where a is the microscopic unit area for which 0 r f r 1. The free energy F can be derived from the corresponding partition function F = Àk B T ln Z, where Þ=kBT is obtained by integrating out the translational degrees of freedom of the bonds. Since the exact calculation of the partition function is out of reach, we invoke a simpler mean-field approach.
Within the mean field approximation, the free energy can be written as where the first term accounts for the mixing entropy of the bonds and the second term represents the mean-field estimation of F N . We are interested in the so called van der Waals regime (see Section I), which is characterized by small thermal roughness D. Following previous studies, 19,28 we will also make the assumption that each adhesion bond causes a deformation h 0 significantly larger than D. This allows us to drop the second term on the r.h.s. of eqn (7) accounting for the entropy of the thermal fluctuations, which leaves only the first term representing the elastic energy of the ground state. The latter can be estimated by considering a lattice of adhesion bonds with spacing r $ ffiffi ffi a p f À0:5 , which gives an energy landscape that depends on the ratio r/x. This approach yields good analytical expressions for the elastic energy only in the limits r/x c 1 and r/x { 1; 19 however, it fails to capture the correct thermodynamic behavior in the intermediate regime r/x B 1 where the lattice distribution does not necessarily represent the energy of a typical random distribution of adhesion bonds. Here, we take a different approach and derive an empirical expression for the dependency of the elastic energy on the bonds' density. We computationally obtain this expression by (i) generating membranes with random, rather than ordered, distributions of adhesion bonds, (ii) finding the membrane profile that minimizes the Helfrich elastic energy of each realization, and (iii) describing the computational data for the elastic energy by a fitting function, which applies to the entire range of densities.

A. Energy calculations
The ground state Helfrich energy corresponding to a random distribution of adhesion bonds is given by the first term on the r.h.s. of eqn (7) and can, in principle, be computed by inverting the coupling matrix (5). In practice, this involves a computationally expensive process and, thus, we adopt a different strategy based on the direct minimization of the Helfrich Hamiltonian. This is done by considering a triangular lattice with lattice spacing l. Each site, i, represents a small membrane segment of area a ¼ ffiffi ffi 3 p l 2 =2, and is characterized by a local height variable h i . On the lattice, N sites are randomly chosen for the locations of the adhesion bonds, at which we set h i = h 0 . The discrete analogue of the Helfrich Hamiltonian (1) is where the discrete Laplacian at site i is given by (i.e., at zero temperature), m € h i ¼ Àa _ h i À @H=@h i , which quickly brings the system to the ground state profile. We measure all lengths in units of the lattice spacing l = 1, and the energy scale is set to k B T = 1. The density of bonds is given by f = N/N s , where N s is the number of lattice sites. Most of the calculations were performed on a triangular lattice of 104 Â 120 sites (with periodic boundary conditions) that has an aspect ratio close to 1. We calculate the elastic energy of numerous random realizations at various densities f r 0.1 and for several values of x varying from x = 5 to x = 10. These values for the correlation length are chosen such that (i) x is sufficiently larger than the lattice spacing l = 1, which reduces the numerical errors associated with the discrete nature of eqn (10) to less than a few percent, and (ii) x is much smaller than the system linear size, to avoid finite size effects.
From eqn (5)-(7) (omitting the second term on the r.h.s.), we infer that for a given set of model parameters (k, h 0 , x, f), the average elastic energy has the form where f (x) is a scaling function of the renormalized density x = x 2 f. Notice that the values of k and h 0 can be fixed arbitrarily since the energy scales like kh 0 2 [see eqn (11)], and this scaling behavior is automatically satisfied by Hamiltonian (10) which is linear in k and quadratic in h i p h 0 . The lowdensity (x -0) asymptotic limit of f (x) is found by considering a system with a single bond, which gives the energy per bond in dilute systems where the typical spacing between the bonds is much larger than the correlation length x. From eqn (7) for N = 1, we read that in this limit, f (x) -1. In the high density limit, i.e., when the spacing between bonds is much smaller than x, the membrane assumes a nearly flat configuration at height h 0 . Setting h i = h 0 in eqn (10) and normalizing the energy by the number of bonds, we obtain the following asymptotic expression akh 0 2 /2fx 4 for F N /N. Using eqn (6) and a ¼ ffiffi ffi this yields the decaying form f ðxÞ ¼ ffiffi ffi 3 p =ð16xÞ for x c 1.
Taking these considerations into account, we propose the following expression for the scaling function: This form ensures the correct asymptotic behavior at low and high densities, and involves two fitting parameters, B 1 and B 2 , to be determined by comparison with the numerical data over the entire range of densities. In Fig. 1 we plot the computational results (triangles) for the elastic energy per bond, normalized by 4k(h 0 /x) 2 , which defines f (x) in eqn (11). The data, which are plotted against the scaled density x = x 2 f, exhibit an excellent data collapse over the entire range x r 10. The solid curve represents the fitting of the data to the form f 1 (x) given by eqn (12), with the parameters B 1 C 5.08 and B 2 C 9.87 that give the best fit. The scatter of the computational data is due to the randomness of the simulated configurations. As expected, the scatter is larger for small values x { 1, where the interaction between the closer pairs of adhesion bonds dominates the energy of the configuration. In fact, for some configurations in this regime, we find f (x) to be slightly larger than unity. This feature is to be expected, and follows from the non-monotonicity of Kelvin's function defining the elements of the coupling matrix M [see eqn (5)]. For x { 1, the PMF between the bonds can be approximated by the sum of pair potentials, as was assumed in ref. 26. By setting N = 2 in eqn (5) and (7), it is easy to confirm that the pair PMF is slightly repulsive at large bond separations. We, therefore, conclude that the scaling function f (x) should be non-monotonic: it first increases for very small values of x, before dropping to zero at larger values. Furthermore, from the fact that Kelvin's function converges exponentially to zero for large arguments, one can also conclude that the derivative of the scaling function df/dx = 0 at x = 0. These features of f (x) in the x -0 limit are Fig. 1 The scaling function for the elastic energy f (x) [see eqn (11)] as a function of the scaled density x. The numerical results are represented by triangles. The solid and dashed curves depict, respectively, the fitting functions f 1 (x) [eqn (12)] and f 2 (x) [eqn (13)] to the data. The inset shows an enlarged view of the data and the fitting functions for x { 1.
not accounted for by the scaling form f 1 (x) proposed by eqn (12). Therefore, we also consider the three fitting parameter scaling function which, in contrast to f 1 (x), correctly captures the behavior of f (x) near x = 0. The fit of the scaling function f 2 (x) to the computational data is also plotted in Fig. 1 (dashed line) with C 1 C 74.8, C 2 C 2174 and C 3 C 1836 which produce the best fit. The difference between f 1 (x) and f 2 (x) is visible only for x C 0, as seen in the inset of Fig. 1. Interestingly, even though f 2 (x) is better suited to represent the scaling function close to the origin than f 1 (x), the latter seems to provide a better fit to the numerical data. In any case, we expect these two functions to yield similar binodal and spinodal curves, except for x C 0. This will turn out to be in the vicinity of the critical point, which is where the validity of the mean-field picture is questionable anyhow.

B. Phase diagram
Plugging eqn (11) into eqn (9), the mean-field free energy, F, of a system with adhesion bond concentration f and correlation length x reads where g(x) = xf(x). With this expression for F, we analytically obtain the spinodal curve, enclosing the region of thermodynamic instability, by solving q 2 F/qf 2 = 0, which yields The binodal curve, which defines the thermodynamic coexistence line, is obtained numerically using a common tangent construction for F. Fig. 2(A) and (B) show the phase diagrams calculated using the scaling functions f 1 (x) and f 2 (x), respectively. In each of these figures, we plot the spinodal curve for x = 5 (solid line) and x = 10 (dotted line), which turn out to be practically indistinguishable. The binodal curves for x = 5 and x = 10 are plotted using squares and circles, respectively. As for the spinodal lines, the binodals for different values of x also overlap with each other. Comparing the phase diagrams presented in Fig. 2(A) [for f (x) = f 1 (x)] and Fig. 2(B) [for f (x) = f 2 (x)], we conclude that the phase diagrams appear to be similar, except for x t 0.6. This is to be expected because only in this regime, the scaling functions are essentially different (see the inset in Fig. 1). Fig. 2

III. Discussion and summary
Looking at the phase diagram depicted in Fig. 2, the one feature that stands out is that the critical point is found at very low densities. The precise value of the critical scaled density x c is, of course, unknown since it depends on the form of the scaling function f (x) [see Fig. 2(C)], and because the mean-field picture is not adequate in the vicinity of the critical point. Nevertheless, it is fair to conclude from the data in Fig. 2 that x c o 0.1, which implies that f c = x c /x 2 { 10 À2 (unless the correlation length is microscopically small, i.e., x B 1). The critical temperature T c can be related to the elastic deformation energy due to a single bond F 1 /k B T = 0.5(h 0 /D) 2 . From Fig. 2 we read that the critical temperature satisfies F 1 C 2 À 3k B T c . Another noticeable feature in Fig. 2 is the fact that the spinodal and binodal curves of membranes with different values of x overlap with each other when plotted against the scaled density x. This does not a priori follow from the data collapse exhibited in Fig. 1, because of the mixing entropy contribution to the free energy. The latter depends on the density f rather than the scaled density x. At low densities, however, we can use the approximation (1 À f) ln(1 À f) C Àf in eqn (9), and then it can be easily shown that the spinodal line [r.h.s. of eqn (15)] becomes only a function of x. Thus, the observation in Fig. 2 that the phase diagram depends on the scaled density is related to the fact that our investigation focuses on membranes with low densities of bonds. The fact that the critical point is located at very low densities means that, slightly below T c , we expect phase coexistence between two low-density phases. From Fig. 2 we also notice that for x \ 1, phase separation occurs only when the temperature drops significantly to roughly T t 0.2T c . This implies that low density systems with large x will not phase separate unless the bonds strongly deform the membrane (h 0 c D). In the two phase region of such a system, the scaled density of the condensed phase x \ 1 which, depending on the value of x, could mean that the density f is quite low. We term low-density (f { 1) regions with scaled density x B 1 as semi-dilute, and conclude that the elasticity-mediated interactions may indeed lead to the formation of such semi-dilute domains.
The ''weakness'' of the elasticity-mediated effect and its inability to induce formation of dense adhesion domains can be understood by looking at the variation of the total elastic deformation energy [second term on the r.h.s. of eqn (14)] with the density of the bonds f. The elastic energy E, normalized per unit area, is plotted in Fig. 3 for membranes with (h 0 /D) 2 = 20 (corresponding to T B 0.2T c ) and x = 10. Also shown in Fig. 3 is the free energy of mixing ÀTS (S denotes the mixing entropy), per unit area, given by the first term on the r.h.s. of eqn (14). Both contributions to the free energy are given in units of the thermal energy k B T. We observe that the total elastic deformation energy increases with f but, somewhat surprisingly, saturates at extremely low densities. The dashed-dotted vertical line in Fig. 3 at f = 0.01 corresponds to x = x 2 f = 1, and one can read from the data that the elastic energy of the membrane barely increases for x \ 0.5. The interpretation of this finding is that one needs a semi-dilute distribution of about one bond per area x 2 to cause the membrane to adopt nearly flat configurations with h B h 0 . Above the scaled density x B 0.5, the membrane elastic energy becomes thermodynamically irrelevant, leaving us with only the mixing entropy term which always favors uniform distributions. This explains why phase separation into regions with distinct concentrations of bonds is possible only at densities below f B 0.5x À2 . To state the last conclusion somewhat differently -the elasticity-mediated PMF induces an attraction between the bonds only if their separation is larger than x. This is an interesting collective (many-body) effect, exhibiting an ''opposite'' trend compared to the pair PMF, which is attractive at separations smaller than x and is screened off at larger distances. The pair PMF may play an attractive role only between two relatively isolated bonds in inhomogeneous distributions, but such configurations fall outside the framework of the mean-field picture presented in this work.
To put our findings in a biological context, we look at the example of the immunological synapse (IS), which forms the contact area between the T-cell lymphocyte and a target cell. Specifically, the cell-cell adhesion is mediated via binding between T-cell receptors (TCR) and MHC-peptide (MHCp) complexes, and between integrin LFA1 and its ligand ICAM-1. 31 These two types of adhesion bonds form a unique structure, in which TCR-MCHp bonds are clustered in its center, while the LFA1-ICAM1 bonds aggregate in the periphery of synapse. It is believed that the central domain, i.e., the TCR-MHCp rich area, plays a pivotal role in T-cell activation. 32 Typically, the bond density within the synapse is around 100 bonds per square micrometer, and the bond lengths are 14 nm and 41 nm for TCR-MHCp and LFA1-ICAM1 bonds, respectively. 33 We recall that in the model presented here, h 0 represents the local membrane deformation imposed by a bond relative to the resting height of the membrane. Thus, if we consider the resting separation between the two membranes in the IS to be dictated by the longer bonds, we can estimate the deformation to simply be the difference between the two bond lengths, h 0 C 27 nm. Taking the membrane bending rigidity to be k C 15k B T and the harmonic potential strength as g C 6 Â 10 5 k B T mm À4 , 34 we arrive at the values x C 0.5 and (D/h 0 ) 2 C 0.057 for the coordinates of this point in the phase diagram displayed in Fig. 2. Remarkably, the point lies in the two-phase region of the phase diagram, close to the binodal line. This raises the possibility that the TCR-MHCp rich domain may be the semi-dilute phase coexisting with a dilute phase of vanishingly small density. Thus, we speculate that the elasticity-mediated interactions may play an important role in the condensation of the TCR-MHCp signaling domain. They provide attraction which enables the TCR-MHCp bonds to spontaneously aggregate into domains with density comparable to that existing in the IS central zone. This finding is in line with several recent studies suggesting that passive thermodynamic processes can describe the short-time condensation of adhesion clusters of the IS, without evoking any active processes in the cytoskeleton (see, e.g., ref. 35 and references therein). Forces stemming from cytoskeletal activity may be essential during the later stages of IS pattern formation and stabilization. 36,37 Introducing such active processes into the equilibrium thermodynamic framework presented here is a task for future studies.