Nucleation of cadherin clusters on cell-cell interfaces

Cadherins mediate cell-cell adhesion and help the cell determine its shape and function. Here we study collective cadherin organization and interactions within cell-cell contact areas, and find the cadherin density at which a ‘gas-liquid’ phase transition occurs, when cadherin monomers begin to aggregate into dense clusters. We use a 2D lattice model of a cell-cell contact area, and coarse-grain to the continuous number density of cadherin to map the model onto the Cahn-Hilliard coarsening theory. This predicts the density required for nucleation, the characteristic length scale of the process, and the number density of clusters. The analytical predictions of the model are in good agreement with experimental observations of cadherin clustering in epithelial tissues.

Fenz et al. 1 quantify the binding and unbinding rates of cadherin monomers to their trans binding partners in their Supplementary materials. They adjust expressions from Dembo 2 and Bell 3 to express these rates as depending on the relative distances between the two membranes h(x i , t): and where k 0 is the intrinsic reaction rate, α the binding range between the outer cadherin domains (taken to be 1nm in their text), b the binding enthalpy of the bond, λ c the elastic constant of the cadherin and ∆h xi ≡ h xi 0 − 2l 0 the instantaneous separation between the outer cadherin domains.
For a cadherin domain to stop growing indefinitely, local membrane fluctuations need to be sufficiently large that k off = k on . This occurs when: Here, the two important parameters are the binding enthalpy b and the cadherin stiffness λ c . Fenz et al. 1 give a value b = 7k B T , which is on the lower end of reported stiffnesses (wild-type cadherin bonds strengthen with time from 9 to 13k B T 4 ). They also use a value of λ c = 5 × 10 −2 k B T nm −2 for cadherin stiffness. This value is on the upper end of stiffness values extrapolated from single molecule cadherin AFM studies 5 , which show that force-extension curves are best-fitted by a persistence length l p = 0.5 ± 0.1nm and a domain length L 0 ≈ 60nm. In the linear regime, the WLC model yields a stiffness λ c ≈ k B T lpL0 ≈ 3 × 10 −2 k B T nm −2 for cadherin. Substituting these values into Eqn. A.3 above, we find that the size of the membrane fluctuations required for two cadherin clusters to dissociate is ∆h xi 0 ≈ 30 − 40nm. This value must now be compared with the observed sizes of membrane fluctuations in cells -which have a well-defined cytoskeleton that modifies the amplitude of membrane fluctuations as opposed to GUVs used in experiments and simulations by Fenz et al. 1 . Biswas et al. 6 find that spatial fluctuations in cells adhering to flat surfaces fluctuate spatially by 7.2 ± 1.5 nm and temporally by 4.9 ± 0.7 nm in adhesive regions (First Branch Regions, which cover most of the cell-surface contact). This value matches with similar results reported by Dos Santos et al. 7 , who also find in their Figures 5 and 6 that fluctuations are lower in focal adhesions -regions where the membrane is closer to the surface.
Cadherins adhere cells to their neighbours, rather than to a smooth matrix, so local fluctuations in the inter-membrane distance will result from active cortexmediated fluctuations from both cells. Whilst it is very difficult to exactly quantify this effect, we expect that it should increase the amplitude of local fluctuations in the separation between the two membranes (assuming two independent Gaussian distributions of fluctuations) by ≈ √ 2, leading to local 10 − 15 nm fluctuations in the inter-membrane separation within adhesion regions. This is much less than the 30 − 40 nm figure required for the formation of cadherin trans-bonds. In other words, most cadherin trans bonds are quickly formed within the larger adhesion area, and thousands of cadherins are brought together by diffusion until a sharp boundary forms around the entire zonula adherens.
However, within the large adhesion region, most cadherin trans dimers have now formed, and simulations by Fenz et al. 1 show that sufficiently small height differences (< 25 nm) between cadherins do not lead to any structure formation, regardless of the size of membrane fluctuations (see Figure SI 4 in their text). Because the membrane fluctuations, and hence the height differences between cadherins are much smaller than 25 nm in larger adhesive areas, cadherins within these regions do not form dense micro-structures without another mechanism, namely direct cis bonds. These computational predictions are corroborated by experiments which show that punctate adhesions do not develop when cis-abolishing V81D/V175D mutations are introduced.

B. CALCULATION OF THE PARTITION FUNCTION
The full partition function is the product of A partition functions for each site, subject to the constraint N = Σ i η i of the constant total number of individual sensors. This condition can be expressed with a delta function, and rearranged into matrix form as: with the non-dimensional coupling constantJ = βJ. Using the N -dimensional Gaussian integral identity, we can write the second exponential term in Eqn.(B.1) as: where x is the Hubbard-Stratonovich vector field, which we need to later obtain the order parameter description. Note the ambiguity in the sign of the term with the exponent linear in x in Eqn. (B.2); both formulations will later be shown to be identical when expressed with respect to the physical mechanosensor density. Treating x as independent of η, we can factorise the partition function (B.1) as follows: 3) The integral of the linear part of Eqn.(B.3) can be simplified: An action functionalS ± [x,μ] can now be defined following the steps of Landau-Ginzburg-Wilson theory: and found to be in the most general form: (B.6) The next step is to calculate the expectation value of the Hubbard-Stratonovich field: After introducing an auxiliary N -component column vector y = (y 1 , ..., y N ) T , we can re-express the Hubbard-Stratonovich field components as: We need a variable whose expectation value can be identified with the average occupation of a site, which is the order parameter: This identifies the discrete field ρ i as the mechanosensor concentration at site i. We will later transform this to a continuous density ρ which depends on a the position along the adhesion contact ring.
An immediate consequence of defining the density field ρ in Eqn. (B.9) is that the previously noted sign ambiguity is lifted: both formulations lead to identical expressions when expressed in terms of ρ. We may therefore without loss of generality only consider the case where the exponent of the linear term e x T η in (B.2) is positive.
It is convenient at this stage to define the fluctuation of the Hubbard-Stratonovich field over all A lattice sites: x . At equilibrium, the action is maximum with respect tõ µ, or equivalently with respect to the modified variable m = iμ (according to the condition of stationary phase). Carrying out the µ-differentiation of Eqn. (B.6) allows us to express m in terms of the Hubbard-Stratonovich field x and the system variables A and N in equilibrium: (B.10)

Model I: Low cadherin concentration
The low-occupancy approximation gives the following condition for the Boltzmann factors: When N A, we can simplify the constraint (B.10) in order to analytically solve it:

This can be rearranged to:
As expected, this is a large quantity if the number of sites is large and the number of sensors is relatively small.

Model II: High cadherin concentration
On the other hand, we find an opposite situation in the high-occupancy approximation: A, we simplify the constraint (B.10): e −x i and rearrange it to: This is a large quantity if the number of sites and the number of sensors are both large.

C. SERIES EXPANSION OF THE ACTION: MODEL I, LOW CADHERIN CONCENTRATION
By explicitly substituting m(x i ) as obtained in Eqn. (B.10) into Eqn. (B.6), we can obtain a useful expression for the action S(x i , m) in terms of x i only. Doing so mathematically incorporates the N = const. constraint on the number of mechanosensing complexes into Ginzburg-Landau action and will allow us to subsequently formulate it in terms of the physical sensor concentration fluctuation φ. We begin by writing the exponential e x i as a series: Summing over all of the states i gives: To keep track of the order of fields, we introduce the parameter a through: We now expand the above in powers of a: A similar trick can be used to track the order of the terms in the expansion of −Σ A i=1 ln [1 + e (m+h+xa+x i ) ], substituting in the expression for the Boltzmann factor. We are interested in the expansion of: Substituting in the expansion of the Boltzmann factor e c , we find the series expansion of this term up to 4 th order in a: Setting the series parameter a = 1, we forget for the time being about the explicit dependency of the terms on A and N , and write: where we make the identifications: The action (or the effective free energy) can be reexpressed in terms of the average value of the Hubbard-Stratonovich field in order to substitute the above: We will also need to check that the total contribution of the sixth order term is positive in order for the Hamiltonian to be bounded from below. Using the connection, Eqn. (B.9), between the Hubbard-Stratonovich field x and the density ρ, we finbally obtain the action S[ρ].

D. SERIES EXPANSION OF THE ACTION: MODEL II, HIGH CADHERIN CONCENTRATION
Here, we write the exponential e −x i as a series: Summing over all of the states i gives: Introducing the order parameter a as above i , we expand the exponential in powers of a: A similar trick can be used to track the order of the terms in the expansion of −Σ A i=1 ln [1 + e (h+m+xa+x i ) ], substituting in the expression for the Boltzmann factor. We are interested in the expansion of: Substituting in the expansion of the Boltzmann factor e c , we find the series expansion of this term up to 4 th order in a: Setting the series parameter a = 1, we forget for the time being about the explicit dependency of the terms on A and N , and write: where we make the identifications: The action (or the effective free energy) can be reexpressed in terms of the average value of the Hubbard-Stratonovich field in order to substitute the above: We will also need to check that the total contribution of the sixth order term is positive in order for the Hamiltonian to be bounded from below. Using the connection, Eqn. (B.9), between the Hubbard-Stratonovich field x and the density φ, we finbally obtain the action S[φ].

E. CONTINUOUS SENSOR CONCENTRATION
The partition function is given by the path integral of the action, given by Eqns. (C.6),(D.6) above, which plays the role of effective free energy in our thermallyequilibrated problem: where the two parts of the action are given by Eqn. (??) after substituting Eqn. (B.9): where the couplings J ij of the original theory have been restored, and the enumeration variable ι (hereafter not explicitly written except for clarity) indicates that the expression can refer to either Model I or II. The order parameter can be transformed into wave vector space, imposing periodic boundary conditions: where s i is the position of the sensor i along the contact ring, and the wave vectors are quantised as k µ = 2πnµ Nµ with n µ = 0, 1, ..., N µ − 1 where N µ is the number of lattice sites in the direction µ (this naturally allows the problem to be extended...) s.t. Π D µ=1 N µ = A is the total number of lattice sites. Note that both real-space and reciprocal space variables are non-dimensional at this stage (as is the lattice size A), with the underlying natural length scale a defining the sensor size. We use the identity: to obtain the Fourier transform of the terms in the truncated effective action: where J k is the Fourier transform of the exchange couplings J ij ≡ J(s i − s j ): φ i and J ij are real, and the couplings are reciprocal: J(−s) = J(s), so φ −k = φ * −k and J −k = J k . Thus in Fourier space, the truncated effective action may be written: .
Sufficiently close to the transition point, only long-range fluctuations (so k small) contribute significantly. So J k can be expanded in powers of k. If the coordination number is z = 2D, or explicitly z = 4 on our contact plane, then: So the quadratic term in the field ρ has the following coefficient: where |T − T c | T c is assumed (so the coefficient of the quadratic term in k simplifies) and the important constants r 0 and c 0 are defined as: The difference in the constants r 0 and c 0 between the two models arises from the difference between g 2,I and g 2,II .
In the limit when the discrete set of allowed wave vectors merges into a continuum, we replace the momentum sums by integrations according to: Normalising the fields by the continuum fields φ(k) = a √ V φ k instead, we define the coupling constants: In the main text, we treat the evolution of a sensor concentration on the adhesion contact plane between a cell and a surface as a 2D problem. It is clear that D = 2, if we set V = A, and set the inter-site spacing to be 1 (scaling all lengths by a), then the above constants simplify considerably. Wavevectors should also be rewritten without any confusion in scalar form.
To check the dimensionality, remember that in this situation, there are in effect not four but three integrals, as one is eliminated in the evaluation of the delta function. The effective action becomes (note, there is an ultraviolet cut-off to the integration, hence the subscript Λ 0 ): Transforming into real space, we obtain the effective free energy for the continuous variable φ(r) and its gradients:

F. BEHAVIOUR AWAY FROM THE TRANSITION POINT
The action in Eqn.(E.11) is hard to analyse without making some simplifying assumptions. We saw in the previous section that fluctuations in the sensor concentration are well-described by a set of sinusoidal perturbations with a given wavenumber k near the transition point. We can examine each of these contributions separately while fluctuations are still sinusoidal (i.e. provided that the cubic, quartic and higher order terms in the action do not dominate over the quadratic term), considering: where we continue using the non-dimensional lengths scaled by the sensor soze a. In the long-time limit, a fastest growing mode will dominate, and so we will approximately find that: where φ max is the amplitude of fluctuations. This enormously simplifies the problem, as the integral in the fourth order term can be approximated as: The Ginzburg-Landau action can then be written in this limit as: Now cos 2 θ ≥ cos 4 θ, ∀θ, so if we show that u 2 > 0 and that the fourth order terms in the action above are positive for k max s = 2mπ, for integer m, then they will be positive at every point s along the boundary of the cell.
Using the results from part C and D, we can show that when k max s = 2mπ, the fourth order terms in both Models combine to: cell boundary, for all sensor concentrations for which our Boltzmann factor approximation is appropriate (so A − N A) and when the modes have grown for a sufficiently long time that the mode closest to the fastest growing spatial frequency k max dominates.
In this limit, a fourth order expansion of the Ginzburg-Landau action is therefore sufficient to describe the destabilization and initial growth of spatial modes (by analysing the quadratic terms) as well as the development of two slightly asymmetrical minima in the action profile (due to the presence of non-zero third order terms) relative to the cell-wide average sensor concentration. This behaviour is examined in Fig. S1. text, we find: Next we solve for the transition point where r 0 = 0. This gives us the value b = b 0 , around which we perform our series expansion: where we use the shorthand: Inserting j = βJ = 7, we find b 0 = N/A = 0.01. The series expansion is asymmetrical: the square root in Eqn. H.1 is only defined when b > b 0 . We convert the wavevector into mode numbers m x /L x ≈ m y /L y ≈ (a k max )/(2 √ 2π) and find a series expansion for the total number of clusters M = m x m y around the transition point as: where the coefficient Ψ I (j) takes the messy form: b0 itself depends strongly on j, so it is not immediately obvious how ΨI (j) will depend on j, but we find in Fig. S2 that it is quite constant above j ≈ 10. In Fig. S3, we show a comparison between the seriesapproximated expression and the exact expression for the number density of clusters at the transition, for A = 20000 (1µm 2 lattice) and βJ = 5. This means that the series expression is good for small changes in contact area (up to 5 − 10%) away from the transition point. . For small changes in (A0 − A)/A, of the order of up to 5 − 10%, the series approximation is good. We find in the main text that cadherin clusters aggregate at a number density M * ≈ 30, so the series approximation for m is good for our purposes.
Finally, we need to re-express Eqn. (H.2) into a more useful form for the reasoning in the main text, which leads to the estimate of the time to destabilize the uniform density field ρ = N/A. Instead of using the ratio b = ρ = N/A, we express m in terms of changes in the size of the lattice ∆A = A0 − A = A(b/b0 − 1). The ratio of mechanosensors to lattice sites only depends on the scaled interaction energy j = βJ, so we write: The experimentally useful value is not the change in the dimensionless lattice size, but rather of the area of the cell. To convert to correctly dimensional quantities, we see that ∆A = t1Σ shrink /a 2 , where a is the distance between lattice sites (7nm). We find that the number of adhesions corresponding to the fastest growing mode number M (t1) depends on the time elapsed since contact area passed the transition point t1: Next we solve for the transition point where r0 = 0. This gives us the value b = b0, around which we perform our series expansion: b0(j) = 1 12j ζ(j) + (16j − 1) + 1 + 16j − 32j 2 ζ(j) with the same shorthand as above: ζ(j) = 3 24 √ 3 48j 6 − 16j 5 + 20j 4 + j 3 − 224j 3 − 48j 2 − 24j − 1