Multiscale approaches to protein-mediated interactions between membranes—relating microscopic and macroscopic dynamics in radially growing adhesions

The membrane is described as a thin sheet with an energy given by the Helfrich–Hamiltionian [41]

$$\begin{eqnarray}&&{{\mathcal{H}}}_{M}[h({\bf{r}})]={\displaystyle \int }_{A}^{}{{\rm{d}}}^{2}{\bf{r}}\;\left(\displaystyle \frac{\kappa }{2}{(\Delta h({\bf{r}}))}^{2}+\displaystyle \frac{\gamma }{2}{\left[h({\bf{r}})-{h}_{0}\right]}^{2}\right).\end{eqnarray} \tag{ 1 }$$

Specifically, the lipid bilayer is parametrized in the Monge gauge, where the height $h({\bf{r}})$ is given as a function of the position ${\bf{r}}$ of the membrane above the substrate (figure 1). A list of the variables and parameters in equation (1) can be found in tables 1 and 2. Specifically, the first term in equation (1) is the deformation energy of the membrane, that is itself a product of the bending rigidity κ and the local mean curvature of the membrane.

Table 1.  Variables and parameters of our Helfrich–Hamiltionian (see equations (1) and (2)).

Quantity	Meaning	Unit
a	Lattice constant	$10\;\mathrm{nm}$
$h({\bf{r}})$	Membrane profile	a
${l}_{0}$	Rest length of the bond/receptor	a
α	Width of interaction potential	a
${{\bf{r}}}_{{\bf{i}}}$	Position of bond i	a
${h}_{0}$	Minimum of the interaction potential	a
A	Area of the simulation box	a2
${k}_{B}T$	Thermal energy at $300\;{\rm{K}}$	$4.14\times {10}^{-21}\;{\rm{J}}$
$\kappa$	Bending rigidity	${k}_{B}T$
${\epsilon }_{b}$	Binding enthalpy	${k}_{B}T$
$\gamma$	Curvature of the interaction potential	${k}_{B}T/{a}^{4}$
λ	Stiffness of the bond/receptor	${k}_{B}T/{a}^{2}$
${N}_{b}(t)$	Number of bonds	—

Table 2.  Parameters used in the simulations.

Parameter	Values in simulation units	Values in SI units
a	−	$10\;\mathrm{nm}$
${h}_{0}$	$8\;a$	$80\;\mathrm{nm}$
α	$1\;a$	$10\;\mathrm{nm}$
${l}_{0}$	$4\;a$	$40\;\mathrm{nm}$
d	$8\;a$	$80\;\mathrm{nm}$
${\rho }_{l}$	$1/64\;{a}^{-2}$	$1.5625\times {10}^{-4}\;{\mathrm{nm}}^{-2}$
A	$4096\;{a}^{2}$	$0.4096\;\mu {{\rm{m}}}^{2}$
${\epsilon }_{b}$	$5\;\;\mathrm{to}\;\;10\;{k}_{B}T$	$2.12\times {10}^{20}\;\;\mathrm{to}\;\;4.14\times {10}^{20}\;{\rm{J}}$
$\kappa$	$10\;{k}_{B}T$	$4.14\times {10}^{20}\;{\rm{J}}$
$\gamma$	$3.125\times {10}^{-3}\;{k}_{B}T/{a}^{4}$	$1.3\times {10}^{-27}\;{\rm{J}}\;{\mathrm{nm}}^{-4}$
λ	$0.75\;\;\mathrm{to}\;\;5\;{k}_{B}T/{a}^{2}$	$3.1\times {10}^{-23}\;\;\mathrm{to}\;\;2.1\times {10}^{-22}\;{\rm{J}}\;{\mathrm{nm}}^{-2}$
η	$2.4\times {10}^{-7}\;{k}_{B}T{\rm{s}}{a}^{-3}$	$1\;\mathrm{mPa}\;{\rm{s}}$
${k}_{0}$	$2\times {10}^{4}\;\;\mathrm{to}\;\;5\times {10}^{5}\;{{\rm{s}}}^{-1}$	$2\times {10}^{4}\;\;\mathrm{to}\;\;5\times {10}^{5}\;{{\rm{s}}}^{-1}$
D	$5\times {10}^{3}\;\;\mathrm{to}\;\;5\times {10}^{5}\;{{\rm{a}}}^{2}\;{{\rm{s}}}^{-1}$	$5\times {10}^{5}\;\;\mathrm{to}\;\;5\times {10}^{7}\;{\mathrm{nm}}^{2}\;{{\rm{s}}}^{-1}$

While the specific protein molecules embedded in the cell wall (or membrane) are usually considered to be responsible for cell adhesion, over the past two decades a realization emerged that the cell membrane itself, being a floppy sheet, adds another unavoidable, yet not fully understood, interaction with the opposing surface it binds to. Although this interaction does not depend at all on any specific proteins, it can have a major impact on the protein-mediated adhesion and can be viewed as a mechanism that controls the binding affinity to the cell-adhesion molecules [42]. Such steric interactions [43] typically maintain the two membranes at relatively large separations h₀, which can be modeled by introducing a nonspecific harmonic potential of a strength γ with the minimum at h₀ [33, 44, 45]. The strength of this potential depends directly on the average intensity of membrane fluctuations that are themselves regulated by the tension in a membrane but also by numerous other factors such as the thickness and the composition of the glycocalyx. In the mimetic systems, this contribution is dominated by continuous interactions between the membrane and the substrate, such as gravity, polymer-repulsion, or Van-der-Waals forces [44, 46]. The strength of this potential can be obtained experimentally by the analysis of membrane fluctuations [47, 48].

We assume that the receptors are thermalized springs with stiffness λ and rest length l₀. This leads to the following expression for the energy of the ${N}_{b}$ bonds in the membrane

$$\begin{eqnarray}&&{{\mathcal{H}}}_{B}[h({\bf{r}})]=\displaystyle \sum _{{\rm{i}}=1}^{{N}_{b}}\delta ({\bf{r}}-{{\bf{r}}}_{{\bf{i}}})\left[\displaystyle \frac{\lambda }{2}{(h({\bf{r}})-{l}_{0})}^{2}-{\epsilon }_{b}\right].\end{eqnarray} \tag{ 2 }$$

Here, ${\epsilon }_{b}$ accounts for the bond enthalpy gain for forming a bond and $\delta ({\bf{r}}-{{\bf{r}}}_{{\bf{i}}})$ is the Dirac-Delta function for the positions ${{\bf{r}}}_{{\bf{i}}}$ of the bonds.

For a given instantaneous membrane profile $h({\bf{r}},t)$, each bond fulfills a local detailed balance condition for the transitions between the bound and unbound state, determined by the rates ${k}^{\mathrm{off}}(h({\bf{r}},t))$ and ${k}^{\mathrm{on}}(h({\bf{r}},t))$ as

$$\begin{eqnarray}&&\displaystyle \frac{{k}^{\mathrm{off}}(h({\bf{r}},t))}{{k}^{\mathrm{on}}(h({\bf{r}},t))}=\mathrm{exp}\left[\displaystyle \frac{\lambda }{2}{\left(h({\bf{r}},t)-{l}_{0}\right)}^{2}-{\epsilon }_{b}-\displaystyle \frac{1}{2}\mathrm{ln}\left(\displaystyle \frac{\lambda {\alpha }^{2}}{2\pi }\right)\right].\end{eqnarray} \tag{ 3 }$$

Here, α is the range of the interaction potential of the ligand–receptor bond and for simplicity, we set $\beta ={({k}_{B}T)}^{-1}\equiv 1$. Equation (3) depends on the stretching energy of the bond (first term in the exponent), the binding affinity (second term) and an entropic contribution (last term) which describes the suppression of fluctuations if a receptor is bound to a ligand. This entropic contribution lowers the effective binding affinity [46], and has to be incorporated if the fluctuations of the receptor are not considered explicitly. The exact form of this term (see appendix B) emerges from the assumption that the structural fluctuations of free receptors occur on a faster time scale than the time scale of the membrane dynamics. This assumption is justified as the dynamics of the receptor is determined by thermal excitations of its secondary structure.

This local detailed balance condition (equation (3)) can be associated with the free energy gain for forming a bond [39], and will govern the stochastic binding and unbinding even if the thermodynamic equilibrium with respect to ligand–receptor binding has not been reached yet.

Inspired by [8, 49, 50], we set the rate ${k}^{\mathrm{on}}(h({\bf{r}},t))$ to form a bond proportional to the probability of the binders being in the binding range α and the intrinsic reaction rate k₀. This naturally depends on the fluctuations, and the barrier for the forward process in the exponent

$$\begin{eqnarray}&&{k}^{\mathrm{on}}(h({\bf{r}},t))={k}_{0}\sqrt{\displaystyle \frac{\lambda {\alpha }^{2}}{2\pi }}\mathrm{exp}\left[-\displaystyle \frac{\lambda }{2}{\left\{\left(h({\bf{r}},t)-{l}_{0}\right)-\alpha \right\}}^{2}\right].\end{eqnarray} \tag{ 4 }$$

From this local on-rate and the detailed balance condition the local, off-rate can be determined readily

$$\begin{eqnarray}&&{k}^{\mathrm{off}}(h({\bf{r}},t))={k}_{0}\mathrm{exp}\left[-{\epsilon }_{b}\right]\mathrm{exp}\left[\lambda \left(h({\bf{r}},t)-{l}_{0}\right)\alpha -\displaystyle \frac{\lambda {\alpha }^{2}}{2}\right].\end{eqnarray} \tag{ 5 }$$

We show the reaction rates equations (4) and (5) in figure 2. The association rate adopts the form of a Gaussian with a width inversely proportional to the stiffness of the receptor. The maximal association rate is obtained at ${l}_{0}+\alpha$, which is at the outer edge of the potential well associated with an unperturbed receptor. The off-rate increases exponentially with the distance between the receptor and the ligand. Interestingly, if the ligand is in the middle of the binding region (${l}_{0}+\alpha /2$), the stiffness of the receptor does not affect the breaking of the bond (the bond is not stressed).

Zoom In Zoom Out Reset image size

Due to the membrane fluidity, the molecules within the bilayer diffuse on its surface [51]. Even though there may be an influence of the membrane elasticity on the diffusion of embedded proteins (for example by the curvature that a protein induces in the membrane [32, 52–54]), these effects seem to be small for experimental relevant parameters [53]. Therefore, we simulate the mobility of binders by a random walk, whereby two proteins interact laterally by a hard-core potential. The time step of the random walk is given by

$$\begin{eqnarray}&&{\tau }_{d}=\displaystyle \frac{{a}^{2}}{4D},\end{eqnarray} \tag{ 6 }$$

with the diffusion constant D. In the current work, only the ligands embedded in the membrane of the vesicle are allowed to diffuse. However, it is straightforward to extend the simulation scheme to situations in which both binders retain lateral mobility and explore the surface of the membrane. The latter may be finite as in the case of vesicles and cells. These situations are simulated using periodic boundary conditions on the level of the system, with a selected area in the center of the simulation box representing the area of contact between two cells or the cell/vesicle and the substrate. Consequently, the formation of bonds can take place only within this region, and the remainder of the system will be depleted from the binders due to the accumulation in the zone of contact.

Binders can be also embedded in bilayers, which provides a constant chemical potential. For simulations of interactions with vesicles, a contact zone is defined, and the periodic boundary conditions are imposed for the bilayer grid. However, to maintain the constant chemical potential (constant concentration of binders in the bulk), entering and exiting of a binder from the contact zone is associated with placing or removing a binder from a random position outside the contact zone.

In this scheme, the membrane shape (see figure 3) is determined explicitly in every time step. Thereby, the system is propagated in time by means of the Langevin equation in Fourier space (see e.g. [4, 32]) derived from the equations (1) and (2)

$$\begin{eqnarray}\displaystyle \frac{\partial h({\bf{k}},t)}{\partial t} & = & -\Lambda ({\bf{k}})\left\{\left[\kappa {k}^{4}+\gamma \right](h({\bf{k}},t)-{\delta }_{{\bf{k}},0}{{Ah}}_{0})\right.\\ & & \left.+\displaystyle \sum _{i=1}^{{N}_{b}(t)}\lambda (h({{\bf{r}}}_{i},t)-{l}_{0})\mathrm{exp}\left(-{\rm{i}}{\bf{k}}\cdot {{\bf{r}}}_{i}\right)\right\}+\xi ({\bf{k}}).\end{eqnarray} \tag{ 7 }$$

Here, $\Lambda ({\bf{k}})$ is the Oseen tensor, describing the hydrodynamic interaction between membrane and surrounding fluid, N_b(t) is the number of bonds of the current bond configuration and A is the area of the membrane in the simulation. The stochastic force $\xi ({\bf{k}})$ in the Langevin equation above is set by the temperature of the surrounding fluid. Thereby, the Oseen tensor is connected to the stochastic force by the fluctuation-dissipation theorem

$$\begin{eqnarray}&&\langle \xi ({\bf{k}})\xi ({\bf{k}}^{\prime} )\rangle =2{k}_{B}T\Lambda ({\bf{k}})\delta ({\bf{k}}+{\bf{k}}^{\prime} ).\end{eqnarray} \tag{ 8 }$$

Zoom In Zoom Out Reset image size

The definition of the Fourier transformation of the Langevin equation is given by

$$\begin{eqnarray}&&h({\bf{k}})={\displaystyle \int }_{A}^{}{{\rm{d}}}^{2}{\bf{r}}\;\mathrm{exp}(-{\rm{i}}{\bf{k}}\cdot {\bf{r}})h({\bf{r}});\quad \quad h({\bf{r}})=\displaystyle \frac{1}{A}\displaystyle \sum _{{\bf{k}}}\mathrm{exp}({\rm{i}}{\bf{k}}\cdot {\bf{r}})h({\bf{k}}).\end{eqnarray} \tag{ 9 }$$

In general, it was shown that the Oseen tensor depends on the geometry of the membrane [55]. However, this dependence is very weak for membranes far away from the substrate and only relevant for the four largest modes of the membrane for the parameters used in the simulations. Such a choice neglects the Helfrich interaction of the membrane with the planar surface [43, 56]. The latter effect should be small for binding of proteins of a size of 20–40 nm when the membrane fluctuates with amplitudes of less than 15 nm, which seems to be the relevant range [9, 57, 58].

Thus, we use the Oseen tensor for a free membrane

$$\begin{eqnarray}&&\Lambda ({\bf{k}})=\displaystyle \frac{1}{4\eta k},\end{eqnarray} \tag{ 10 }$$

where η is the viscosity of the surrounding fluid. The Oseen tensor for the ${\bf{k}}=0$ mode diverges. Following [4], the Oseen tensor for this mode is set to

$$\begin{eqnarray}&&\Lambda ({\bf{k}})=\displaystyle \frac{3\sqrt{A}}{8\pi \eta }.\end{eqnarray} \tag{ 11 }$$

The Langevin-equation (7) is solved numerically with the help of the Euler–Maruyama scheme (see for example [59]). The time step in this scheme has to be set below the smallest time scale of the membrane

$$\begin{eqnarray}&&\tau ({k}_{\mathrm{max}})=\displaystyle \frac{4\eta {k}_{\mathrm{max}}}{\kappa {k}_{\mathrm{max}}^{4}+\gamma },\end{eqnarray} \tag{ 12 }$$

which is the typical relaxation time of the mode with the largest ${\bf{k}}$ in the simulation (${k}_{\mathrm{max}}=\sqrt{2}\pi /a\approx 4\;{\mathrm{nm}}^{-1}$) and is on the order of 10⁻⁹ s.

The simulation is performed following the algorithm shown on the left in figure 4. The first step initializes the system. This involves the thermal equilibration of a free membrane obtained by executing 10⁶ steps in the time loop explained below without the reactions and binder diffusion. After that, the ligands are placed randomly on their lattice, and the receptors are put on a grid of the second lattice.

Zoom In Zoom Out Reset image size

The second step is the initialization of the time loop, where the step accounts for the shortest characteristic membrane time scale ($\Delta t\equiv \tau ({k}_{\mathrm{max}})$). Every time step involves (i) the calculation of the force on the membrane induced by the formed bonds in real space; (ii) the transformation of this force to the Fourier space; and (iii) the determination of bending and unspecific forces in Fourier space (first term in equation (7)). The sum of this forces is input to the Euler–Maruyama step, within which the membrane profile is updated in Fourier space and transformed back to real space. This back-transformation is a prerequisite for the execution of the association and dissociation step. Here, the binding probabilities are obtained from the equations (4) and (5), in which the height of the membrane and the time step of the simulation are required. As the time scale of the reactions is much larger than the typical time scale of the membrane, these probabilities are rather small.

Finally, the mobile ligands need to be displaced to one of the neighboring unoccupied sites. In principle, the diffusion of binders is characterized by the time step given by equation (6), in which case the probability to jump in any direction would be $1/4$. However, as the time scale of the diffusion is typically several orders of magnitude larger than the step of the time loop (i.e. ${\tau }_{D}\gg \tau ({k}_{\mathrm{max}})$), the probability of a jump is rescaled to

$$\begin{eqnarray}&&p=\displaystyle \frac{\tau ({k}_{\mathrm{max}})}{4{\tau }_{D}}.\end{eqnarray} \tag{ 13 }$$

This new probability guarantees the correct diffusive behaviour of the ligands. Here we note that this implementation decouples the simulation of diffusion from the simulation of the membrane. As a result, the ligand shape and the coupling of diffusion to the membrane curvature is not considered, in contrast to earlier works [52, 53, 60]. However, in the context of adhesion, the curvature coupling is a small effect, justifying this approximation.

In the current simulations, the ligands are immobilized after they form a bond with a receptor, which means that only free ligands diffuse. This restriction is motivated by the experimental observation that the bonds change position only if they are subject to a significant lateral force [61]. After the diffusion has been resolved, a new iteration in the time loop is started, or the simulation is terminated.

Computationally, most time in this simulation scheme is consumed by fast Fourier transformations of the membrane profile and the forces, which scale like $N\mathrm{log}(N)$, (N is the number of considered lattice points), and not linearly like other operations (diffusion and reaction kinetics). Furthermore, the time step has to be chosen very small to accurately describe the time evolution of the membrane, and a large number of replicas must be produced to obtain a statistically sound representation of the system. These are the main reasons which make this simulation scheme computationally very expensive allowing only for length scales of up to $1\;\mu {{\rm{m}}}^{2}$ to be simulated for up to ${10}^{-1}\;{\rm{s}}$.

The difficulties that arise with Langevin simulations could be circumvented if the explicit treatment of the membrane could be avoided. We achieve this goal in an effective Monte Carlo scheme which is based on the recently acquired understanding of the effects of the membrane on the formation of bonds [39, 46, 48]. This scheme relies on the fact that the typical time scale of the membrane fluctuations depends on the viscosity η of the surrounding fluid

$$\begin{eqnarray}&&{\tau }_{\mathrm{mem}}=\displaystyle \frac{4\eta {q}_{0}}{\kappa {q}_{0}^{4}+\gamma }=\displaystyle \frac{2\eta {q}_{0}}{\gamma }\simeq 2\times {10}^{-5}\;{\rm{s}}.\end{eqnarray} \tag{ 14 }$$

Here, ${q}_{0}={(\gamma /\kappa )}^{1/4}$ is the inverse lateral correlation length for a membrane without tension.

Importantly, even the slowest modes are significantly faster than the reaction kinetics for ligand–receptor binding (the fastest avidin-biotin in membranes was reported to take place at $\sim {10}^{3}\;{{\rm{s}}}^{-1}$ [40, 42]), while other pairs are found at $\sim {10}^{2}\;{{\rm{s}}}^{-1}$ [8, 42]. Consequently, the membrane fluctuations can be regarded as equilibrated with fixed mean shape as long as the configuration of bonds interacting with the membrane remains unchanged. During this time, the fluctuating membrane, and with it the ligands, sample the entire probability distribution of distances between ligands and receptors. In the following, we denote the height distribution at the considered binding site ${\bf{r}}$ before the bond has formed by $p({h}^{r})$, and the height distribution after a bound ligand–receptor pair is formed by $p({h}^{b})$. Here h^r is the height at a binding site before a bond is established and h^b is the height after a bond is established. The distribution of the latter is nontrivial if the receptor or the bond itself maintains some flexibility.

The first and the second moment of these typically Gaussian distributions (equations (15) and (16)) can be calculated explicitly [46] yielding the mean height h^s and the fluctuation amplitude ${\sigma }^{s}$ of the membrane at the binding site for an arbitrary bond configuration. Here the superscript s denotes the state at the binding site, which can be either a free receptor (s = r) or an existing bond (s = b). Specifically, we calculate a functional integral over all membrane profiles weighted by their Boltzmann factor (see appendix A for technical details). As result, we obtain the mean height

$$\begin{eqnarray}&&\langle h({\bf{r}})\rangle \equiv {\bar{h}}^{s}({\bf{r}})=\displaystyle \frac{-\displaystyle \frac{4}{\pi }{\sum }_{{ij}}{l}_{0}{G}_{{ij}}{({\bf{r}})}^{-1}\mathrm{kei}\left({q}_{0}\ | {\bf{r}}-{{\bf{r}}}_{{\bf{i}}}| \right)}{8\sqrt{\kappa \gamma }+\displaystyle \frac{16}{{\pi }^{2}}{\sum }_{{ij}}\mathrm{kei}\left({q}_{0}\ | {\bf{r}}-{{\bf{r}}}_{{\bf{i}}}| \right){G}_{{ij}}{({\bf{r}})}^{-1}\mathrm{kei}\left({q}_{0}\ | {\bf{r}}-{{\bf{r}}}_{{\bf{j}}}| \right)}\end{eqnarray} \tag{ 15 }$$

and the fluctuation amplitude

$$\begin{eqnarray} & \langle {h}^{2}({\bf{r}})\rangle -{\langle h({\bf{r}})\rangle }^{2}\equiv {\left({\sigma }^{s}\right)}^{2}({\bf{r}})\\ & ={\left(8\sqrt{\kappa \gamma }+\displaystyle \sum _{{ij}}\displaystyle \frac{16\mathrm{kei}\left({q}_{0}\ | {\bf{r}}-{{\bf{r}}}_{{\bf{i}}}| \right){G}_{{ij}}{({\bf{r}})}^{-1}\mathrm{kei}\left({q}_{0}\ | {\bf{r}}-{{\bf{r}}}_{{\bf{j}}}| \right)}{{\pi }^{2}}\right)}^{-1}.\end{eqnarray} \tag{ 16 }$$

The sum runs over all pairs of bonds in the membrane at the positions ${{\bf{r}}}_{{\bf{i}}}$ and ${{\bf{r}}}_{{\bf{j}}}$, while $\mathrm{kei}(x)$ is the Kelvin function [62]. The elements of the coupling matrix ${G}_{{ij}}({\bf{r}})$ are the effects of the existing bonds on the shape and fluctuations at the arbitrary position ${\bf{r}}$, whereby the membrane mediated interaction between the bonds are comprised in the off-diagonal elements (see equation (A.9) for the explicit form of the matrix).

The average shape and the fluctuation amplitude of a membrane containing a small cluster of bonds are shown in the top panels of figure 5. At large distances from the cluster, the membrane is on average flat since it resides and fluctuates in the minimum of the nonspecific potential. Because of a relatively high concentration of bonds within the cluster the membrane is likewise flat on average, but much closer to the substrate. At the same time, its fluctuations are strongly suppressed. However, the shape and fluctuations of the membrane are significantly different in the vicinity of the bonds at the edge and in the center of the cluster.

Zoom In Zoom Out Reset image size

We use the two height distribution functions to average the Bell–Dembo rates (equation (4) and (5)) at the position of a free or a bound receptor giving rise to effective binding and unbinding rates

$$\begin{eqnarray}&&{K}^{\mathrm{on}}\equiv \displaystyle \int \;{\rm{d}}{h}^{r}p({h}^{r}){k}^{\mathrm{on}}({h}^{r})\;\\ &&{K}^{\mathrm{off}}\equiv \displaystyle \int \;{\rm{d}}{h}^{b}p({h}^{b}){k}^{\mathrm{off}}({h}^{b}).\end{eqnarray} \tag{ 17 }$$

Appropriately inserting equations (4), (5), (15), and (16) into the the above expression, and evaluating the integrals yields

$$\begin{eqnarray}&&{K}^{\mathrm{on}}={k}_{0}\displaystyle \frac{\sqrt{\lambda {\alpha }^{2}}}{\sqrt{2\pi (1+\lambda {\left({\sigma }^{r}\right)}^{2})}}\mathrm{exp}\;\left[-\displaystyle \frac{\lambda {\left[{\bar{h}}^{r}-(\alpha +{l}_{0})\right]}^{2}}{2(1+\lambda {\left({\sigma }^{r}\right)}^{2})}\right],\\ &&{K}^{\mathrm{off}}={k}_{0}\mathrm{exp}\left[\left(\displaystyle \frac{\lambda \alpha }{2}\{2\left({\bar{h}}^{b}-{l}_{0}\right)+\alpha [\lambda {\left({\sigma }^{b}\right)}^{2}-1]\}-{\epsilon }_{b}\right)\right].\end{eqnarray} \tag{ 18 }$$

In this manner, the effective rates describing the association and the dissociation of a bond depends on the exact position of a bond, and the time dependent configuration of all bonds in the system. Examples of such rates for one bond configuration can be seen in bottom panels of figure 5. Obviously, the rates reflect the average shape and fluctuations within the membrane (top panels), which are the result of the bond configuration around the respective binding site. The dissociation rate of a bond at the rim is up to two orders of magnitude larger than for a bond deep within the domain (see figure 5). This is due to the stabilization effects of the neighboring bonds, which share the deformation load and cooperatively suppress the fluctuations. On the other hand, the association rate is the largest near the bond domain and exponentially decreases on the length scale of the lateral correlation length with increasing distance to the domain.

To summarize, cooperative effects between bonds and binding sites are induced by the deformation of the membrane in the vicinity of the bond and the changes in membrane fluctuations. In the absence of cooperative effects, the rates are constant at every binding site (see figure 3, [39] for additional details), because the deformation is either not reaching the first neighbors (very soft membranes) or the membrane does not deform at all (receptors significantly more flexible than the membrane). However, if the formation of a bond induces changes in the membrane, the effect is relatively short range.

We can now construct a Monte Carlo simulation of the adhesion process in which only the reaction kinetics and the diffusion of binders must be treated explicitly. Thereby, the effective rates for breaking or forming a bond at the given site must be determined for each site in every time step. This in turns requires inverting the coupling matrix containing all bonds, for every site in every step. In order to make the simulation fast, we assume that only the first two sets of neighboring bonds affect the rates on a particular (un)binding site (figure 6). Consequently, only the configuration of bonds in the immediate environment is taken into account in the calculation of the effective rates. Since this environment consists only of 9 sites, all possible configurations can be explored a priori, and their respective rates used to create a lookup table. This restriction to the next-nearest neighbours is justified because the binding rates decay very fast with increasing distance between the bonds.

Zoom In Zoom Out Reset image size

A flow chart of the Monte Carlo scheme is shown on the right panel of figure 4. To initialize the system, all ligands and receptors are positioned on their respective grids as in the Langevin simulations (random and ordered distributions are possible). Furthermore, the characteristic time steps are determined. The time step of the simulation is given by the characteristic diffusion–time $\Delta {\tau }_{D}$ (equation (6)). The time step for the reaction kinetics is set to be $\Delta {t}_{B}=\Delta {\tau }_{D}/n$, where n is the smallest integer satisfying the inequality ${K}^{\mathrm{on}/\mathrm{off}}\Delta {\tau }_{D}/n\lt 1$. From here the probabilities for binding and unbinding are calculated as ${K}^{\mathrm{on}/\mathrm{off}}\Delta {t}_{B}$, and stored in a lookup table.

The simulation step starts with the reaction loop which consists of n iterations. In each iteration, for every binding site (i) the bond configuration is determined, (ii) the appropriate rate is retrieved from the lookup table, (iii) association or dissociation is attempted, and (iv) the bond configuration is updated. Following the reaction loop, each binder attempts to move to a neighboring site in a same manner as in the Langevin scheme. This completes the simulation step and the system is propagated in time until the program is terminated. While the program allows for the diffusion of both binder types, the following discussion will be restricted to the case when the receptors are immobilized.

The advantage of the Monte Carlo scheme is that it allows for a larger time step and avoids fast Fourier transformations limiting the Langevin code. This allows us to simulate length scales of several tens of micrometers and time scales of several seconds with the resolution of about 100 nm and 10⁻⁵ s, which is necessary to understand biological processes.

We first focus on the simulation of rare events such as is the nucleation of adhesion domains. The number of bonds in such a domain can be calculated explicitly within the capillary approximation [39]. Once this number is estimated we perform about 2000 simulations with each method to generate the distribution of nucleation times (figure 7). Specifically, each simulation is set to start from an equilibrated box with zero bonds. When a cluster of bonds of critical size is formed anywhere in the system, the simulation is interrupted, and the time necessary to achieve this domain size is recorded. As shown in figure 7, very good correspondence of the coarse-grained and the higher-level simulation approach is obtained for the entire distribution of nucleation times. This agreement could have been anticipated from the successful comparison of Langevin simulations with the analytic model for the nucleation dynamics of a single seed, based on a simplified version of the here used effective rates [39]. The current, more accurate approach fully validates the concept of the effective rates and enables studies of the early stages of the adhesion process in the regimes that are either not accessible to analytic modeling or are extremely demanding from the computational point of view. Examples of such regimes, which can be now addressed with ease, are fast nucleation, competitive growth of multiple seeds, or diffusion limited nucleation.

Zoom In Zoom Out Reset image size

Encouraged by our results reproducing the nucleation dynamics, we validate the Monte Carlo scheme by reproducing the results of the higher level scheme for the full dynamic adhesion process, i.e. nucleation, growth and saturation to equilibrium. More specifically, for each set of parameters we perform 200 runs over which we average the dynamic process. This level of accuracy was found previously to produce converged results for the Langevin scheme in thermal equilibrium [32, 33].

We first explore the correspondence of the two schemes when the diffusion of ligands is fast ($D=5\times {10}^{7}\;{\mathrm{nm}}^{2}\;{{\rm{s}}}^{-1}$), for soft, moderately stiff, and stiff receptors (figure 8). In each graph, the number of bonds as a function of time is presented for three different binding affinities (the smallest being at the phase transition to the unstable adhesion dominated by unbinding, two intermediate affinities, and one large affinity where the unbinding is negligible). In general, we find that the two approaches show very similar dynamics, irrespective of the choice of parameters.

Zoom In Zoom Out Reset image size

Somewhat larger deviations (up to 15%) are observed at the phase boundary (${\epsilon }_{b}=6.79\;{k}_{B}T$) between the 'stable' and 'unstable' adhesion. Here, stable adhesion denotes a state with an equilibrium number of bonds that is always larger than zero, whereas in the 'unstable' state the membrane occasionally unbinds and the number of bonds spontaneously drops to zero [33]. In this regime, the simulations are more sensitive to the size of the micro-environment, unlike in the rest of the phase space where considering the first two sets of nearest neighbors produces results which are very similar, and more accurate than considering larger neighborhoods with up to 20 sites.

Very similar results are obtained for slow diffusion of ligands (figure 9), where the adhesion dynamics is shown for three receptor rigidities, at an intermediate binding affinity. Equally good, quantitative agreement is obtained for all affinities above the transition energy (data not shown). These results validate convincingly the concept of effective rates, and establish the Monte Carlo approach as a reliable and versatile method for the simulation of protein mediated membrane interactions.

Zoom In Zoom Out Reset image size

It should be noted that different effective Monte Carlo schemes, based on the integration of membrane fluctuations in the Hamiltonian were successful in comparison with the Langevin simulation [37], with the time scale of reactions being a free fitting parameter. However, the accuracy of that scheme relied on the magnitude of the effective cooperativity parameter χ to be much smaller than one. This dimensionless parameter evaluates the fluctuations of the unbound membrane with respect to the fluctuations of free receptors

$$\begin{eqnarray}&&\chi \equiv \displaystyle \frac{\lambda }{8\sqrt{\kappa \gamma }}.\end{eqnarray} \tag{ 19 }$$

The accuracy of the current scheme does not depend on the effective cooperativity parameter, which for the systems shown in figure 9 range from 0.53, for the softest receptors, to 3.54, for the stiff receptors. Actually, the regime of large effective cooperativity parameters seems to be very important in the context of experiments with cells or vesicles [63].

We perform a series of Monte Carlo simulations, where we use two opposing square grids of a size of 40.96 × 40.96 μm² in the diffusion limited case and of 10.24 × 10.24 μm² in the reaction limited case (typical sizes of a giant unilamellar vesicle). The first grid carries $2.5\times {10}^{5}$ receptors (soft or stiff), immobilized on a lattice. To simulate diffusion or reaction limited growth, the second grid is decorated by randomly placing 5 × 10⁴ diffusing ligands or placing immobile ligands above the receptors, respectively. These concentrations, as well as the other parameters are strongly inspired by the analogous experimental realizations of the system [42]. Specifically, the height of the membrane (${h}_{0}-{l}_{0}=55\;\mathrm{nm}$), curvature of the nonspecific potential ($\gamma =3.125\times {10}^{-3}\;{k}_{B}T/{a}^{4}$), bending rigidity of the membrane ($\kappa =10\;{k}_{B}T$), binding affinity (${\epsilon }_{b}=10\;{k}_{B}T$), intrinsic reaction attempt frequency (${k}_{0}={10}^{5}\;{{\rm{s}}}^{-1}$) and the diffusion constant ($D\;=\;$ 5 μm²s⁻¹) is chosen such that the nucleation of domains and the unbinding of bonds are rare. Furthermore, we investigate the reaction and diffusion limited growth regimes for stiff ($\lambda =5\;{k}_{B}T/{a}^{2}$) and soft receptors ($\lambda =2\;{k}_{B}T/{a}^{2}$) mimicking bulky cell adhesion receptors and glycoprotein receptors, respectively.

For ligand densities larger than the receptor density, we expect a quadratic growth of the domain area [19, 42, 65] containing N_b uniformly distributed bonds

$$\begin{eqnarray}&&{N}_{b}(t)=\pi {({K}_{R}^{\mathrm{on}})}^{2}{a}^{4}{\rho }_{l}^{2}{t}^{2}.\end{eqnarray} \tag{ 20 }$$

Here, ${\rho }_{l}$ is the initial density of ligands and ${K}_{R}^{\mathrm{on}}$ is the effective rate at the rim. This expression explicitly takes into account the two-dimensional nature of the growth process.

The results of our Monte Carlo approach (blue full lines in figures 10(c) and (d) confirm that the growth of the domain is, quadratic as expected. This is confirmed by the very good agreement of the data which we fit by equation (20), shown in figure 10 with dashed orange lines. The observed processes show that growth is faster for stiff (${K}_{R}^{\mathrm{on}}$ = 3.7 × 10⁴ s⁻¹) than for soft (${K}_{R}^{\mathrm{on}}$ = 2.0 × 10⁴ s⁻¹) receptors, presumably because of stronger correlations between bonds. Clear deviations from the quadratic behavior take place when the finite size effects start to play a role and the domain begins to cover the whole simulation box.

Zoom In Zoom Out Reset image size

In order to relate the rates extracted from the fit to the microscopic rates which were actually used to grow the domains, we construct bubble charts for binding and unbinding rates (figures 10(e)–(h)), which are classified by the number of neighbors. A fixed number of neighbors can be organized in several different configurations around the receptor of interest, which results in the multiple bubbles for each number of neighbors. In the bubble charts, the area of the bubble is associated with the occurrence of a particular rate in the simulation.

Interestingly, the effective rate ${K}_{R}^{\mathrm{on}}$ corresponds very well to the average rate recorded in the simulation. Actually, for the stiff bonds the average rate at the rim, obtained by averaging all rates forming with up to five neighbors ${\bar{K}}^{\mathrm{on}}$, is ${\bar{K}}^{\mathrm{on}}$ = 3.7 × 10⁴ s⁻¹, and for soft bonds ${\bar{K}}^{\mathrm{on}}$ = 2.5 × 10⁴ s⁻¹. Rates for the formation of bonds with three to five bonds in the neighborhood are most commonly observed (largest bubbles), which is consistent with the formation of new bonds at the edge of the domain. The rates for the formation of bonds with six or more neighbors are considered to be the results of events from rebinding within the domain, in agreement with the large number of dissociation events with seven and eight neighbors (figures 10(g) and (h)).

The analysis of microscopic rates in the bubble plots shows that the binding rates have a tendency to increase up to five neighbors. This happens because the formation of additional bonds, in principle, reduces the distance between the receptor and the ligand at the position of the binding site. The rates for forming the bond with 3–5 neighbors are significantly larger for stiff receptors, which is the source of the difference in the speed of the overall growth process of the domain. The reason for this difference is that for stiff receptors, the membrane approaches closer to the substrate than for soft receptors, which themselves deform while forming a bond, leaving the membrane at a larger height. Rates for forming a bond with 6–8 neighbors decrease with increasing the number of adjacent bonds. This effect is more significant for stiff receptors, because the fluctuations in the membrane are suppressed to a larger extent, and while the distance from the receptor is relatively small, stronger thermal membrane excitation is necessary to bring the ligand into the reaction zone of the receptor. The unbinding rates occur less frequently. The most common unbinding rate is the one with eight adjacent bonds, which is clearly associated with unbinding within the domain. The unbinding rates decrease exponentially as a function of the number of neighbors, for both stiff and soft receptors, showing the stabilization effects that binding in the surrounding has on the respective bonds.

For ligand densities lower than the receptor density, the growth is diffusion limited (figure 11), and depends only implicitly on the effective reaction rates through the density of ligands and bonds at the edge of the domain. In other words, the growth explicitly depends only on the diffusion constant, and the area of the growing domain A(t) is given by [19, 21]

$$\begin{eqnarray}&&A(t)=4\pi {\alpha }^{2}{Dt}.\end{eqnarray} \tag{ 21 }$$

In this equation, α is a dimensionless speed factor (in contrast to earlier sections where it was the width of the interaction potential between the binders), which is, in two dimensions, determined from the implicit equation [21]

$$\begin{eqnarray}&&\displaystyle \frac{{\rho }_{l}-{\rho }_{e}}{{\rho }_{b}}={\alpha }^{2}\mathrm{exp}\left({\alpha }^{2}\right)\mathrm{Ei}\left({\alpha }^{2}\right).\end{eqnarray} \tag{ 22 }$$

Here, ${\rho }_{e}$ is the density of ligands at the edge of the domain and ${Ei}(x)$ is the so-called exponential integral [62]. This relation is obtained from the binder conservation at the rim of the domain and the respective solution of the diffusion equation (see Shenoy and Freund [21] for details).

Zoom In Zoom Out Reset image size

We numerically solve equation (22) using the densities of bonds and ligands at the edge of the domain evaluated from the radially averaged density profiles. We obtain $\alpha =0.34$ for stiff receptors and $\alpha =0.27$ for soft receptors. The difference in the speed factors emerges from the somewhat larger density of ligands at the rim of a domain with soft receptors in the simulation. This difference is due to the smaller binding rates at the rim of the domain with soft receptors. Using these speed factors, we can calculate the expected diffusion constant from the linear fit (orange dashed lines in figure 11). Specifically, we obtain a diffusion constant of 4.8 ± 0.6 μm² s⁻¹ for the large bond stiffness, and 5.2 ± 0.7 μm² s⁻¹ for the low bond stiffness (right column of figure 11). The error is due to the error of the bond density ${\rho }_{b}$ and the density at the rim ${\rho }_{e}$. The relative error of the values of the diffusion constant compared with the diffusion constant of the simulation (5.0 μm² s⁻¹) is 4% which is within the accuracy of the fits. This confirms our initial hypothesis that in this density regime the diffusion limited growth is well accounted by the mean field approach.

The obtained results from the simulations of the growth of ligand–receptor domains show that it is, in principle, possible to relate macroscopic measurements with the underlying microscopic processes. From diffusion limited processes we can extract the diffusion constant with excellent accuracy, which is also accessible from experimental data. However, as noted before [42], issues may arise if the crossover to the saturation of the growth curve due to the finite size of the vesicle or cell occurs relatively quickly and the vesicle runs out of free binders. Furthermore, it is possible to relate the mean reaction rate to the microscopic events.

The above presented analysis allows us to validate the Monte Carlo simulations, and in addition to show how the parameters characteristic for the mean field modelling relate to the microscopic stochastic processes of binding and diffusion. However, the behaviour described in this section is characteristic only in particular limits, whereas different growth laws are valid in different parts of the phase space. For example, if ligand–receptor bonds already form during the formation of the initial contact zone between the vesicle and the substrate, the dynamics has a different scaling behaviour [66]. The methods presented herein are a necessary prerequisite for the systematic study of the growth laws and will enable deeper understanding of the domain development in the regimes which are not accessible to analytic modelling.