Localization of proteins to a membrane is an essential step in a broad range of biological processes such as signaling, virion formation, and clathrin-mediated endocytosis. The strength and specificity of proteins binding to a membrane depend on the lipid composition. Single-particle reaction-diffusion methods offer a powerful tool for capturing lipid-specific binding to membrane surfaces by treating lipids explicitly as individual diffusible binding sites. However, modeling lipid particle populations is expensive. Here, we present an algorithm for reversible binding of proteins to continuum surfaces with implicit lipids, providing dramatic speed-ups to many body simulations. Our algorithm can be readily integrated into most reaction-diffusion software packages. We characterize changes to kinetics that emerge from explicit vs implicit lipids as well as surface adsorption models, showing excellent agreement between our method and the full explicit lipid model. Compared to models of surface adsorption, which couple together binding affinity and lipid concentration, our implicit lipid model decouples them to provide more flexibility for controlling surface binding properties and lipid inhomogeneity, thus reproducing binding kinetics and equilibria. Crucially, we demonstrate our method’s application to membranes of arbitrary curvature and topology, modeled via a subdivision limit surface, again showing excellent agreement with explicit lipid simulations. Unlike adsorption models, our method retains the ability to bind lipids after proteins are localized to the surface (through, e.g., a protein-protein interaction), which can greatly increase the stability of multiprotein complexes on the surface. Our method will enable efficient cell-scale simulations involving proteins localizing to realistic membrane models, which is a critical step for predictive modeling and quantification of in vitro and in vivo dynamics.

Proteins localize to membranes via multiple different binding modes, including recognition and binding to highly specific lipid head-groups [e.g., PI(4,5)P2],1 electrostatically driven adherence to negatively charged membranes2 (as performed by many BAR domain proteins),3,4 and binding to membrane-inserted proteins (e.g., Ras).5,6 Once proteins have localized to the surface, they can be further stabilized by interactions with additional lipids or transmembrane proteins, and these subsequent binding events are effectively two-dimensional (2D) in their search.6,7 Capturing these various forms of membrane binding is critical for effective kinetic and spatial modeling of cell-scale systems with quantitative comparison to the experiment, both in equilibrium and nonequilibrium systems. Because molecular modeling approaches that could capture the fundamental electrostatic or hydrophobic nature of these interactions are too expensive at this scale,8,9 rate-based approaches such as reaction-diffusion provide a powerful alternative.10–14 These modeling approaches can be used not only to study systems with known experimental binding constants15 but also as a means of fitting and extracting information from new experiments on protein-membrane interactions. We present a new algorithm for efficient rate-based and reversible binding of proteins to flat or curved continuum membranes. We then characterize how treating lipids as explicit particles vs an adsorptive surface can alter reaction kinetics and equilibria, and illustrate our implicit lipid (IL) model’s ability to recover the behavior of the full explicit lipid (EL) method.

Our model offers a hybrid explicit-particle scheme that combines the efficiency of surface adsorption models with the enhanced spatial and molecular details of single-particle binding sites (Fig. 1). In a simple surface adsorption model, the surface is a boundary where adsorption occurs upon collision using an adsorption rate κ with units of length/time (see, e.g., Ref. 16). The more expensive but ultimately more detailed and flexible approach is to give the surface a density of reactive particles, and solution species can react according to a bimolecular reaction rate kon with units of volume/time. This second model, generally known as the Langmuir model,17 assumes a 1:1 binding stoichiometry between the surface species and the solution species. Both surface adsorption16,18–23 and Langmuir7,11,23,24 models have been used for single-particle reaction diffusion (RD) and partial differential equation (PDE)-based modeling.25,26 The Langmuir model has a clear advantage in that it treats surface binding sites as a variable, accounting for heterogeneity and exclusion of other solution particles from binding occupied sites. This second model is also the only one that can naturally account for subsequent binding to surface sites as a 2D reaction, that is, after a solution particle has first been localized via a separate binding site.7 However, for single-particle RD, propagating individual membrane binding sites (such as lipids or patches of lipids) is quite costly. Lipids are quite abundant. For example, at about 1 mol. % at the plasma membrane, the number density of PI(4,5)P2 is 25 000 µm−2,27 translating to about 50–80 µM given a yeast cell surface area and volume.

FIG. 1.

Schematic diagrams of the implicit lipid model. Binding sites or lipids on the surface, indicated by the NL black points (a), are replaced by a uniform field with density ρL = nL/S [(b) and (c)]. For a solution particle (green points), instead of evaluating binding to each of the lipids within a separation r0 < Rmax, we can instead evaluate the binding to the surface by integrating over all positions on the surface that are accessible (r0 < Rmax) given a height z0 (b). The separation r0=r2+z02 cannot fall below σ, hence a donut-shaped patch on (c). (d) Effects of the reflecting-surface RS3D. A diffusing particle cannot pass through the reflecting-surface, unless a reaction binding occurs in which case it is put on the membrane surface (denoted as S) directly. The dissociated particle is moved from the membrane surface and put on the reflecting surface.

FIG. 1.

Schematic diagrams of the implicit lipid model. Binding sites or lipids on the surface, indicated by the NL black points (a), are replaced by a uniform field with density ρL = nL/S [(b) and (c)]. For a solution particle (green points), instead of evaluating binding to each of the lipids within a separation r0 < Rmax, we can instead evaluate the binding to the surface by integrating over all positions on the surface that are accessible (r0 < Rmax) given a height z0 (b). The separation r0=r2+z02 cannot fall below σ, hence a donut-shaped patch on (c). (d) Effects of the reflecting-surface RS3D. A diffusing particle cannot pass through the reflecting-surface, unless a reaction binding occurs in which case it is put on the membrane surface (denoted as S) directly. The dissociated particle is moved from the membrane surface and put on the reflecting surface.

Close modal

While more expensive than continuum PDE25,26 or stochastic reaction-diffusion master equation (RDME) approaches,28–31 single-particle RD simulations are unique in capturing excluded volume,13,14,32,33 capturing fluctuations driven by small copy numbers and diffusional collisions,11,13,14,34 and enabling simulations of macromolecular self-assembly,35 thus providing a mechanism to bridge cell-scale modeling with molecular modeling approaches.36 For single-particle binding to surfaces, our method provides critical advantages over existing adsorption methods. Several algorithms for surface adsorption based on discretizing the diffusion equation require short time steps to provide accurate solutions.12,20,21,37 Methods using analytical solutions to the diffusion equation with a reactive boundary provide accurate solutions even for larger steps16,18,19,23,38 but have restricted parameter regimes,16 are model specific,38 and in all cases, are derived only for planar surfaces.12,16,18–21,23,38,37 Surfaces with mixed boundaries (inhomogeneity) have been considered for receptor arrays on planar surfaces,38,39 but these receptor disks are immobile. Perhaps most importantly, unlike these existing surface adsorption models, our implicit lipid model treats the density of the target lipid as a variable, thus naturally retaining the ability of the explicit lipid model to capture binding occupancy, nonhomogeneous distributions, and 2D lipid binding. Critically, it applies on flat and curved surfaces, retains accuracy with large time steps, and applies to reversible binding. This approach also provides a natural way to incorporate single-particle RD with hybrid continuum approaches,40–42 where the surface density is naturally a continuous variable. We achieve this versatility in our model by deriving reaction probabilities for a particle in solution to bind to a surface distribution of lipids that replace individual binding sites with an approximation by a continuous field (Fig. 1). By defining reaction probabilities based on the Green’s function (GF) solution to the Smoluchowski model, we can take relatively large time steps (on the order of microseconds) while correctly accounting for all possible collisions with the surface. Importantly, we propagate solution particles using simple Brownian updates such that our method can be integrated with most single-particle simulation tools.11,12,14,24 The desorption probabilities we define enforce detailed balance to reproduce the proper equilibrium distribution for reversible reactions. Our approach applies to all types of surface models, with any reaction parameters, and can be widely adopted by existing single-particle methods.

By comparing adsorption and explicit lipid models of surface binding, we show how the choice of model can lead to significant differences in the kinetics of membrane association and occupancy at the membrane. Surface adsorption models can be designed to recover the same equilibrium occupancy of the Langmuir model simply by making the adsorption rate κ a function of lipid density, κ(ρL) = kon · ρL. This useful formula highlights the relationship between the two models and guides the comparison between them. As suggested by this formula and shown below, for parameter regimes with a high density of surface binding sites, both models produce nearly identical results. With limited binding sites, explicit lipid models produce substantially slower kinetics. Our implicit lipid method, on the other hand, shows excellent agreement with the kinetics of explicit lipid models. We further compare the microscopic explicit lipid model with macroscopic binding models using spatially unresolved ordinary differential equations (ODEs). In small and rate-limited systems, these macroscopic rate equations produce quite comparable kinetics to the microscopic kinetics of the explicit lipid model, with results diverging for diffusion limited reactions.

By establishing kinetic data for a range of reaction parameters and surface topologies, from simple (planar) to complex using spatially resolved methods, we provide a standard for the comparison of other rate-based models across simulation software implementations. Curved surfaces present several challenges relative to flat surfaces, particularly, if they are not analytically defined (e.g., a sphere or a cylinder). We achieve binding to surfaces of arbitrary topology by using a sublimit surface mesh model,43 and by defining accurate methods for detection of the surface and reflection off the surface,44,45 and diffusion on the surface.46 In this paper, we first introduce the single-particle model for binding to explicit particles embedded in surfaces, before deriving our new implicit lipid model and the associated algorithm. We define relationships between rates of surface adsorption models, the explicit lipid model and our implicit lipid model to establish accurate and comparable definitions of macroscopic rates. We show how the models produce distinct kinetics in low lipid regimes, with converging kinetics at high surface coverage. Our method can be incorporated with both existing single-particle methods and hybrid PDE based models, as the solution particles can bind to a spatially dependent density of surface particles. We further demonstrate how 2D binding to surface sites (i.e., via localization by separate interfaces) can be done with our implicit lipid model, which is not possible with surface adsorption models. We apply this method to heterogeneous densities of surface lipids and to curved mesh surface models as an important step in integrating RD methods with continuum deformable membrane models.43 Finally, we note that while none of these rate-based models capture the molecular properties of lipid bilayers (e.g., electrostatic attraction and repulsion), some molecular features are implicitly represented by the lipid density and binding rate. This simpler model is a necessary starting point for ultimately capturing how binding responds to collective changes of the surface, such as curvature.3,47,48

1. Binding to explicit lipids on a flat surface

Binding between a particle in solution at r1 to a lipid particle restricted to a reflective 2D surface at r2, thus transitioning the 3D solution particle to the 2D surface, follows the same general theoretical approach as between two particles in solution (Fig. 2). The Smoluchowski model for binding between reactive pairs has been quite effectively used for single-particle RD with large time steps;13,14,24 we define the model here as there is a modification required relative to standard treatments of two solution particles49,50 due to the reflective surface. The diffusion equation for the particle pair can be written as

pr,t|r0t=Drpr,t|r0,
(1a)

with their separation r = r1r2 and initial separation r0, and p(r, t|r0) is referred to as the Green’s function (GF). The diffusion tensor D(r) = D1 + D2 is a generally nonisotropic matrix because the diffusion coefficients on the surface are nonzero in only 2 dimensions (for 2D diffusion on a horizontal plane, D2z = 0). The reaction is modeled via the reactive boundary at r = σ,

Jr^|r=σ=κpr=σ,t|r0,
(1b)

where J = −D(r)p(r12, t|r120) is the flux at r and r^ is the direction vector of r. The reactivity κ occurs over the accessible surface of a sphere of radius σ such that κ = ka/B(σ, d), where ka is the microscopic (or intrinsic) binding rate constant. Unlike in a full three-dimensional (3D) solution, where B(σ, d = 3) = 4πσ2, here only a hemisphere is accessible, 2πσ2, due to the additional boundary of the reflective surface,

Dr1pr1,t|r10r1nsurfΩ=0,
(1c)

where Ω indicates the boundary surface and nsurf is the direction vector of the surface. This affects the macroscopic rate as described in Sec. II C. The other boundary condition is

p(r=,t|r0)=0,
(1d)

and the initial condition is a delta function, p(r, t = 0|r0) = δ(rr0)/B(σ, d).

FIG. 2.

Model for binding to explicit lipids on a flat surface. (a) Geometric schematic in the spherical coordinate system. The green point represents a protein particle in solution, and the black point represents the lipid particle on the membrane. The green particles reflect off the membrane surface (blue color). This constrains the polar angle θ to be between 0 and π/2, while the radius r and azimuth angle φ are not influenced. (b) Reaction probability between the particles depicted in (a) as a function of time, given an initial separation r0 = 1 nm. Results from Brownian dynamics simulations of (a) are shown in red and blue, at different initial polar angles, with the average over all orientations in green. The theoretical solution provided by Eq. (2) is shown in the black dashed line, agreeing extremely well with the BD average in green, using the weighted D¯ = 1.667 nm2/μs. For comparison, we show the theory with D = D1x + D2x = 2 nm2/μs, which is slightly lower than the average. Results for ka = 347.18 nm3/μs and σ = 1 nm.

FIG. 2.

Model for binding to explicit lipids on a flat surface. (a) Geometric schematic in the spherical coordinate system. The green point represents a protein particle in solution, and the black point represents the lipid particle on the membrane. The green particles reflect off the membrane surface (blue color). This constrains the polar angle θ to be between 0 and π/2, while the radius r and azimuth angle φ are not influenced. (b) Reaction probability between the particles depicted in (a) as a function of time, given an initial separation r0 = 1 nm. Results from Brownian dynamics simulations of (a) are shown in red and blue, at different initial polar angles, with the average over all orientations in green. The theoretical solution provided by Eq. (2) is shown in the black dashed line, agreeing extremely well with the BD average in green, using the weighted D¯ = 1.667 nm2/μs. For comparison, we show the theory with D = D1x + D2x = 2 nm2/μs, which is slightly lower than the average. Results for ka = 347.18 nm3/μs and σ = 1 nm.

Close modal

For a lipid particle diffusing in 2D in a planar surface, the problem is perfectly symmetric around the plane, e.g., at z = 0. Reflection simply flips the polar angle θ to remain between 0 and π/2 without affecting r or φ [Fig. 2(a)]. The Green’s function (GF), p(r, t|r0), can thus be defined from the full 3D problem, with the only complication being the anisotropy of the diffusion matrix D (Dx = DyDz). Unfortunately, although the diffusion anisotropy can be removed by a change in variables, this then introduces an anisotropy on the radiative boundary condition (the z-axis is rescaled),23 and thus the GF can only be solved exactly as a function of initial and final particle separation and angular orientation.

Instead, we will approximate the dynamics as occurring under an isotropic D¯, weighted by the x, y, and z components so that our algorithms can use the analytically known GF solution in 3D dependent only on particle separation r, p(r, t|r0).14 As usual, the reaction probabilities for a pair of particles at a separation of r0 are defined based on the survival probability such that

preactΔtr0=1SΔtr0=1σ4πr2pr,Δtr0dr,
(2)

producing the same functions used in 3D. To define D¯, we note that the anisotropic diffusion in x, y, and z means that the diffusion in spherical coordinates along the relative separation r is a function of the angular orientation, D(θ, φ) = Dx cos2(φ)sin2(θ) + Dy sin2(φ)sin2(θ) + Dz cos2(θ). Since Dx = Dy, then we have D = Dx sin2(θ) + Dz cos2(θ). If we weight D over all polar angles, then D¯=0π/2D sin θdθ=2/3Dx+1/3Dz. We verify that using the known 3D reaction probabilities with this D provides excellent agreement with the exact reaction probabilities determined using Brownian dynamics simulations [Fig. 2(b)]. Brownian dynamics simulations were run using the method of Zhou,51,52 similar to recent work,35,53 with the addition of the reflecting surface around the lipid particle.

2. Binding to explicit lipids on a curved surface

For binding to a lipid diffusing on a curved surface, the system is no longer perfectly symmetric across the reflecting surface, and thus the reaction probabilities are not exactly described via the GFs discussed above. Here, we will simply approximate the reaction probabilities based on the GFs described above, thus ignoring the effect of the local curvature on the reactive flux at r = σ. The curvature of the reflecting surface would have two main effects on the GF solutions. First, reflection off a concave surface, for instance, tends to redirect particles closer to the surface compared to reflections off flat or convex surfaces, resulting in a small enhancement to reactive collisions (see the supplementary material for further details). Second, for a concave surface, the accessible surface of a sphere of radius σ around the lipid particle would be slightly less than that of the hemisphere for a flat surface, resulting in a small reduction in net reactive flux. These effects are thus somewhat compensatory. Additionally, the effect of the curvature is reduced by propagating a smaller time step, rendering the flat-surface GF solution above a more and more accurate approximation. In our results below, we find that evaluating binding probabilities to a particle on a curved surface assuming a flat reflection model of Sec. II A 1 does not introduce statistically measurable deviations in the known equilibria, suggesting it is a good approximation.

3. Dissociation

As performed in previous work,35 dissociation between explicit particles is treated as a Poisson process,

Pdissoc,SP(Δt)=1exp(kbΔt),
(3)

where kb is the microscopic dissociation constant. To recover the detailed balance and thus the proper equilibrium, dissociation must return particles to the same separation where the reaction occurs, which for the Smoluchowski model is at contact, r = σ.

The explicit lipid model for binding to surfaces described above achieves high resolution of the spatial distribution and its dependence on diffusion. However, it is expensive to propagate all of the lipid particles; as small, densely packed constituents of membranes, they are often much more numerous than solution molecules. Our new implicit lipid model replaces the lipid particles by a field or density distribution (see Fig. 1). This approach reduces the spatial details of the lipids but eliminates the need to propagate the lipids on the membrane surface. Diffusion of solution particles can still be propagated using simple Brownian updates, and binding probabilities are closed-form analytic functions for planes and spheres. An important aspect of this approach is that the depletion of the lipid particles due to binding of proteins can be captured in the lipid density, ρL = nL/S (nL, the number of unbound lipids and S, the area of membrane surface). Once a binding event occurs, nL decreases, which will reduce the density and change the reaction probability for subsequent events. Dissociation events return lipids to the pool and increase the density. By integrating over a surface patch of the membrane accessible to the solution particle, thus considering all possible collisions with the surface, this model can be applied to surfaces of varying curvature and lipid density, and to particles restricted to 2D. We enforce detailed balance by specifying an additional reflection surface parameter, RS.

1. Derivation of the implicit lipid model for a planar surface

For a solution particle that interacts with a field of lipid particles located on a planar surface at z = 0, where first we assume a fixed density of ρL = NL/S, we can define the implicit lipid survival probability for the solution particle in a time Δt, given a height above the membrane of z0, by integrating over all collisions with the surface (Fig. 1). This is given by the product of the pairwise survival probability between the solution particle and each position on the surface, weighted by the lipid density,

SILΔt|z0=r=r1r2SΔt|r02πrρLΔr=expr1r22πrρLlnSΔt|r0dr,
(4)

where from the geometry of the problem, we define r0=r2+z02 and r1 = 0 if z0σ, and r1=σ2z02 if z0 < σ (see Fig. 1). The well-known reaction probability preactt|r0) = 1 − St|r0) [Eq. (2)] can only be evaluated for r0σ14 [Fig. 1(c)]. The upper integration limit r2 can be set to infinity, although it can be convenient to truncate the integral. Because the integrand drops to zero rapidly, we define a cutoff based on the 3D distance where the particles are too far apart to collide in a time step, Rmax=σ+36DΔt, which gives r2=Rmax2z02, where D = D¯ (Sec. II A 1).

Then, the reaction probability for the implicit lipid model to a planar surface is written as

PILΔt|z0=1SILΔt|z0=1expr1r22πrρLln1preactΔt|r0dr.
(5)

This integral must be solved numerically. However, if the time step Δt is relatively small (order of microseconds) and the lipid density ρL on the membrane is relatively low, then the reaction probability is very well approximated analytically as

PILΔt|z0r1r22πrρLpreactΔt|r2+z02dr=2πρLσ2kakd(ka+kd)2expα2+2αZerfcα+Z2αZ+1erfcZ+2απexpZ2,
(6)

where α=DΔtka+kd/σkd, Z=z0σ/4DΔt, kd = 4πσD, and ka is the intrinsic binding rate. When z0 < σ, the probability is independent of the specific value of z0 and is in fact always given by PILt|z0 = σ). A typical curve of PILt|z0) is shown in Fig. 3, where Eq. (6) is indistinguishable from the numerical solution to Eq. (5) and is faster to evaluate. This approximation works because one can show that Eq. (5) is given by Eq. (6) minus corrections to prevent multiple events per step. These corrections are negligible for small Δt.

FIG. 3.

The reaction probability for the implicit lipid model, PILt|z0) [red line—Eq. (5)], can be extremely well approximated by a closed form analytical function [blue dashed line—Eq. (6)]. Parameters: σ = 1 nm, D¯ = 12.67 nm2/μs, Δt = 0.1 µs, ρL = 0.025 nm−2, ka = 347.18 nm3/μs, and kb = 2.09 s−1.

FIG. 3.

The reaction probability for the implicit lipid model, PILt|z0) [red line—Eq. (5)], can be extremely well approximated by a closed form analytical function [blue dashed line—Eq. (6)]. Parameters: σ = 1 nm, D¯ = 12.67 nm2/μs, Δt = 0.1 µs, ρL = 0.025 nm−2, ka = 347.18 nm3/μs, and kb = 2.09 s−1.

Close modal

The dissociation of particles from lipids follows a Poisson process,

PdissocILΔt=1expkoffΔt,
(7)

where unlike the explicit lipid model, here, we must use koff, the macroscopic dissociation rate that satisfies koff = konKD = kon · 2kb/ka, and kon, a simple function of ka, D, and σ as described in Sec. II C. This is necessary because the explicit lipid method returns particles exactly to contact, where the possibility of geminate rebinding creates a distinction between kb, the intrinsic off rate, and koff, the macroscopic off rate. For rate-limited reactions, kb and koff are the same, but for diffusion-influenced reactions, where ka is large, kb is faster than koff (kbkoff). Here, there are no explicit lipids to drive geminate rebinding, and thus the dissociation follows the macroscopic rate.

The explicit solution particles are propagated using free diffusion. Unlike in the FPR algorithm,14 here, we do not reweight the reaction probabilities as the position updates for particles that did not bind result from an effectively many-body interaction with the surface, rather than a pairwise reaction. To guarantee the correct equilibrium in the implicit lipid model, one final parameter is needed, the position of the reflecting surface, which excludes particles from getting within a distance RS3D of the surface. Detailed balance requires that at the equilibrium state, during a time step Δt, the number of solution particles binding to the surface should be equal to the number of particles dissociating from the surface. This is expressed as

NAnAPdissocILΔt=RS3DPILΔt|z0nAVSdz,
(8)

where V is the volume of the system, S is the area of the membrane surface, NA is the total number of solution particles, nA is the number of unbound solution particles at the equilibrium state, and RS3D is the position of the reflecting surface above the membrane. The above function has an analytical solution

RS3D=σ1expkbkdΔtka+kdka+kd24πσ2kbkd4DΔt12αexp(α2)erfc(α)1+2απ+α2exp(α2)erfc(α)1+2α/π,
(9)

where α and kd are the same as in Eq. (6). If both kb and Δt are relatively small, then RS3D has the simple expression,

RS3Dσkaka+kd,
(10)

where kd = 4πσD. We can see that RS3D has a value of 0 < RS3D < σ. The reflecting surface depends only on σ, ka, and diffusion constant D. It does not depend on the system size or particle numbers, as we would expect for this reaction model.

In sum of the implicit lipid model, solution particles move according to free diffusion and measure distances to the membrane surface properly, but their position updates must exclude them from occupying the layer RS3D above the membrane, via reflection. When dissociation occurs, particles are returned to a distance of RS3D above the membrane surface. Then, using the reaction probability of Eq. (6) and the dissociation probability of Eq. (7), we recover the correct equilibrium and kinetics very similar to the explicit lipid model, as shown below.

2. Implicit lipid model: Modification for a curved membrane

If the membrane is not planar, the binding probability must account for the curvature of the membrane. This alters the distance relative to the rule used above, where previously r0=r2+z02. Now, we use the more general expression for Eq. (6),

PILΔt|r0=SpreactΔt|rρLds,
(11)

where S represents a patch on the surface accessible to the diffusing solution particle within a step ∆t, r0 represents the spatial position of the particle in solution, and r is the distance from the particle to each point within the surface patch. Here again, S need not be the entire surface, as the reaction probability will drop to zero for r > Rmax. In addition, the surface integral must be evaluated over distances rσ. This integral must be computed numerically to account for the changes in distances over the surface patch, but can be thought of as the average of the reaction probabilities over all separations, multiplied by the lipid density and the surface area of the patch. For a spherical surface, the reaction probability in Eq. (11) is analytically solvable, providing faster evaluations of binding probabilities (supplementary material II).

The dissociation probability is not dependent on curvature. However, the reflecting surface RS3D depends on the reaction probability, Eq. (11), which changes slightly with curvature. Following the conservative detailed balance, we need to combine Eqs. (7), (8), and (11) to calculate the value of RS3D. In the supplementary material, we show the calculation for binding to both the interior and exterior of a sphere (positive and negative curvatures), illustrating that RS3D can be reproducibly simplified to the same result of Eq. (10). One can thus generally assume that when the membrane curvature is not extreme over the length of the time step, Eq. (10) is a very good approximation of RS3D.

3. Implicit lipid model for 2D reactions

An important feature that cannot be captured by surface adsorption models is binding to lipids by proteins that are already localized to the surface. This can happen either because proteins are localized to the surface through a protein-protein interaction or because a protein has multiple lipid binding sites. For these 2D reactions, the relative displacement z0 above the membrane is always zero such that, using the simpler formulation of Eq. (6), the integral is given by

P2DILΔt=RS2Dr22πrρLpreact,2DΔt|rdr.
(12)

Here again, the upper limit of integration can be chosen to a point where the reaction probability drops to zero at Rmax=σ+34DΔt. RS2D is the reflecting distance (we will introduce it below), which similar to 3D, is a parameter used to enforce detailed balance. The reaction probability requires using the known solution to the 2D association problem, preact,2Dt|r) (see Ref. 53), which is itself an integral over Bessel functions. By taking preact,2Dt|r) into Eq. (12), the binding probability is expressed as

P2DILΔt=4kaρL0expDu2Δt1u3γ2+ω2γaωbdu,
(13)

where γ = huY1() + kaY0(), ω = huJ1() + kaJ0(), a = uRmaxJ1(uRmax) − uR2DJ1(uRS2D), and b = uRmaxY1(uRmax) − uR2DY1(uRS2D), h = 2πσD. Jn denotes the order n Bessel function of the first kind, and Yn denotes the order n Bessel function of the second kind. The integral must be solved numerically. We use the GNU scientific library (GSL) quadrature routines for numerical evaluation.

The dissociation of particles from lipids still follows a Poisson process,

Pdissoc,2DILΔt=1expkoffΔt,
(14)

where koff is the macroscopic dissociation rate that satisfies koff = konKD = konkb/ka and a macroscopic association rate kon for 2D is determined based on our previous work,53 

kon=1ka+18πD4 ln(bρ/σ)(1σ2/bρ2)221σ2/bρ211,
(15)

where bρ=2S/πmax(nA,nL)+σ2, where S is the membrane area, nA is the number of unbound particles, and nL is the number of unbound lipids.

Similar to the 3D case, we need to introduce a reflecting distance, RS2D that will enforce detailed balance between binding and unbinding reactions, which will be ≥σ. At the equilibrium state, during a time step Δt, the number of particles binding to lipids must be equal to the number of particles dissociating from lipids, which can be expressed as

NAnAPdissoc,2DILΔt=nAP2DILΔt.
(16)

Taking Eqs. (13)–(15) into Eq. (16) and considering that kb is usually small, we have

0exp(Du2Δt)1u3(γ2+ω2)(γaωb)du=konΔt4ka,
(17)

where γ, ω, a, and b are defined as before in Eq. (13) (a and b are functions of RS2D). By numerically solving Eq. (17), we obtain the value of RS2D. It is usually larger than σ (Fig. S1).

On a curved surface, the approach is similar to that of the 3D case, where we explicitly integrate over all points along the surface in terms of their geodesic distance from the protein particle,

P2DILΔt|r0=Spreact,2DΔt|rρLds.
(18)

Note that the value of r will start from RS2D. The dissociation is still defined by Eq. (14), and the reflecting distance RS2D will be calculated by Eqs. (14), (16), and (18).

4. Implicit lipid model for nonuniform lipid densities

In the above examples, the lipid density was used as a fixed variable ρL. Nonuniform lipid densities could arise due to phase separation54 or some other confining mechanism. Our method can straightforwardly accommodate such inhomogeneous lipid densities: the same integrals of Eqs. (5), (6), and (11), or (18) apply but the integration over reaction probabilities is weighted by a spatially varying lipid density, ρL(r, φ), where φ is the azimuthal angle. While this requires numerical integration even for Eq. (6), this enables application of this model to a hybrid PDE-based method with explicit solution particles, where the lipid density on the surface could be defined as a continuous and diffusive distribution. To facilitate direct comparison with the explicit lipid simulations, we demonstrate this approach where the lipid densities on the surface are restricted to zero in one region and a fixed value ρL in the other region.

For both the explicit lipid model and the implicit lipid model that is derived from it, the relationship between the microscopic or intrinsic binding rate ka and the macroscopic rate kon is different from the known relationship for both particles in solution due to the reflective surface. The total reactive flux across the radiation boundary at r = σ is reduced, and for a planar surface, the reduction is by exactly 1/2. As a result, the macroscopic rate differs from the rate observed in 3D solution by a factor of 1/2 for the same value of ka. Specifically, the macroscopic rate is known to be related to the survival probability through51 

kt=Γκr0Str0ds0,
(19)

where r0 are all the points on the reactive surface Γ. For a spherically symmetric interaction, the reactive surface is at r0 = σ and φ = [0, 2π], with the reflective surface restricting the polar range to θ = [0, arccos(σc/2)]; thus, the total reactive surface area is 2πσ2(1 − σc/2). Here, c represents the local curvature of the surface, which can be positive or negative (c = 1/R for a sphere of radius R) or for a planar surface, c = 0. The survival probability is provided by the full 3D problem, where κ = ka/4πσ2. As a result,

k(t)=λkaS(t|σ),
(20)

where the reduction factor λ = 1/2 (1 − σc/2). Given the well-known survival probability in 3D,14 the macroscopic rate is

kon=k(t)=λ1ka+14πσD1
(21)

or λ times the rate that is produced in solution (1/2 for a planar surface). This means that for a single-particle RD simulation to produce a specific macroscopic rate, Eq. (21) must be used to extract the ka value, which will be up to a factor of λ−1 higher than for a solution reaction.

Critically, this results in a similar modification to define the equilibrium constant from the intrinsic rates, which applies for any reaction where one reactant is localized to a reflective surface. In 3D, we have that KD=koff/kon=kb/ka3D. For transitions between solution and the surface, we have that

KD=koffkon=λ1kbka3D2D.
(22)

This requires care when reactant pairs both start off in solution with a specified ka3D and one becomes localized to the surface. Subsequent binding reactions using ka3D2D = ka3D will follow the weaker KD of Eq. (22), despite involving the same two reactants. We noted in a recent publication that this correction is important in models of membrane localization to producing a proper equilibrium, rather than nonequilibrium, steady-state.7 To account for this, we use a simple approach in our software. All ka values are defined based on a 3D reaction, and if one reactant is localized to the membrane, we multiply ka by the zero-curvature value of λ−1 = 2 to preserve the same value of KD (ka3D2D = 2ka3D), although this results in slightly different kon values. In Table I, we list the relationships between the binding rates, unbinding rates, and equilibrium constant, for 3D, 2D, and the 3D → 2D transition. In Fig. S2, we illustrate the application of these relationships by comparing the microscopic explicit lipid simulations with an ODE solution using the macroscopic rates. Finally, we note that the macroscopic rates for a curved surface will not be exactly given by these formulas, as the survival probability and the size of the reactive surface change slightly with curvature (Sec. II A 2).

TABLE I.

Relationships between microscopic and macroscopic reaction rates.

Systems3D3D → 2D2D
Relation to 3D microscopic binding rate ka3D ka3D2D=2ka3D ka2D=ka3D/2σ 
Microscopic dissociation rate kb kb kb 
Macroscopic binding rate kon=1ka3D+14πσD1 kon=121ka3D2D+14πσD1 kon=1ka2D+18πD4 ln(bρ/σ)(1σ2/bρ2)221σ2/bρ211 
Macroscopic dissociation rate koff=konkb/ka3D koff=2konkb/ka3D2D koff=konkb/ka2D 
Equilibrium state KD=kb/ka3D KD=2kb/ka3D2D KD=kb/ka2D 
Systems3D3D → 2D2D
Relation to 3D microscopic binding rate ka3D ka3D2D=2ka3D ka2D=ka3D/2σ 
Microscopic dissociation rate kb kb kb 
Macroscopic binding rate kon=1ka3D+14πσD1 kon=121ka3D2D+14πσD1 kon=1ka2D+18πD4 ln(bρ/σ)(1σ2/bρ2)221σ2/bρ211 
Macroscopic dissociation rate koff=konkb/ka3D koff=2konkb/ka3D2D koff=konkb/ka2D 
Equilibrium state KD=kb/ka3D KD=2kb/ka3D2D KD=kb/ka2D 

In adsorption models, collisions with the surface result in “binding” with an adsorption rate κ with units of length/time. To compare against models of bimolecular association to surface sites with rate kon (the explicit or implicit lipid models), we can use κ = ρLkon, with ρL being the density of lipids per area. We show results here, both where κ is fixed and where κ(ρL(t)) varies as the lipid density changes in time. Desorption is treated like a standard dissociation reaction. An efficient and spatially dependent model for adsorption/desorption in single-particle reaction-diffusion was introduced by Andrews in Ref. 16. Because this method cannot be used to simulate reversible adsorption when κ is low (requires κ2 ⩾ 4koffD), it cannot be used for arbitrary binding strengths, as we do in some simulations here. We thus briefly describe the application here of the Smoluchowski model in 1D to surface adsorption, which has been already used in published algorithms.18,23,38

1. Adsorption using the Smoluchowski model in one dimension

For adsorption of a diffusing particle to a planar surface, the height above the surface is the only variable influencing the reaction. We can, therefore, describe the particle position based on the product of the free diffusion propagator in x and y, and diffusion in z subject to a radiation boundary at the surface (e.g., z = σ). Specifically,

pirrx,y,z,t|x0,y0,z0=pfreex,t|x0pfreey,t|y0pirrz,t|z0,
(23)

where pfree(x,t|x0)=exp[(xx0)2/(4Dt)]/4πDt and the same for y. The solution to the 1D Smoluchowski model pirr(z, t|z0), parameterized by the diffusion constant D = Dz and the rate κ (with adsorption units of length/time), is known to be19,49

pirrz,t|z0=14πDtexp(zz0)24Dt+exp(z+z02σ)24DtκDexpκ2tD+κDz+z02σ× erfcz+z02σ4Dt+κtD,
(24)

where σ can be used to move the surface to any arbitrary height. The reaction probability is then given by

PAD-GFΔt|z0=1σpirrz,Δt|z0dz=erfcz0σ4DΔtexpκDz0σ+κ2ΔtD× erfcz0σ4DΔt+κΔtD.
(25)

This probability depends on z0, which can then be compared with the result for PILt|z0) derived above from the full 3D GF. The probability distributions are not identical (Fig. S3), which is not surprising, considering that one has 3D diffusion and the other does not. If we set σ = 0 for this model, for a small intrinsic binding rate ka, they are comparable (Fig. S3A), but for a large ka, they are more distinct (Fig. S3B). Therefore, it is perhaps not surprising that the many-body simulations we show below produce better agreement for rate-limited (small ka) reactions. Importantly, for a curved surface, the height above the surface is not independent of x and y diffusion, and thus this model would have to be adjusted to correct for this. We thus apply it here only to planar surfaces. Finally, for this single-particle method, desorption is again treated as a Poisson process, and the equilibrium for a solution population of A particles will be defined by κ/kdesorb = [A]mem/[A]sol.

To evaluate our method, we set up four systems to represent different membrane surfaces or different characteristics of lipid distributions. Figure 4(a) shows a planar surface in a cube, while in Fig. 4(b), the lipids are distributed nonuniformly, occupying only half of the bottom plane. Figures 4(c) and 4(d) present systems with curved membrane surfaces, both a sphere and a complex surface with four joined spheres that present a range of curvature over convex and concave regions.

FIG. 4.

Schematics of four systems. A flat membrane surface (a), half-planar surface (b), spherical surface (c), and irregular curved surface (d) are shown. All the membrane surfaces are in blue color.

FIG. 4.

Schematics of four systems. A flat membrane surface (a), half-planar surface (b), spherical surface (c), and irregular curved surface (d) are shown. All the membrane surfaces are in blue color.

Close modal

The implicit lipid model always recovers the proper equilibrium, thanks to the use of the detailed balance expression in Eq. (10) to define the reflection surface RS3D (Fig. 5). For simulations with 200 solution particles and a range of lipid densities and reaction-rate parameters, we further show that the kinetics are in excellent agreement with the results of explicit lipid simulations (Fig. 5). In Fig. 5(b), we simulate a diffusion limited reaction (ka3D = 173.59 nm3/μs) at increasing lipid densities, illustrating how even with a low concentration of lipids (100/2002 nm−2), the kinetics are nearly the same. In Fig. 5(c), we simulate a rate-limited reaction (ka3D = 0.166 nm3/μs). Because rate-limited reactions are relatively insensitive to diffusion and the spatial distributions of reactants, the close agreement between the models here is less surprising. Rate-limited reactions also produce similar kinetics in a nonspatial, ODE simulation of the same system (Fig. S2A). Diffusion to the surface does start to have an impact when the separation from the surface can reach 2 μm, but the change is relatively small compared to the substantial difference observed in diffusion-limited reactions (Fig. S2B). For diffusion-limited reactions, most collisions are reactive, and thus the assumption of the ODE that the lipids are well-mixed, compared to localized to the surface, produces much faster kinetics.

FIG. 5.

The implicit lipid (IL) and explicit lipid (EL) models with the planar surface shown in Fig. 4(a) are in very close agreement for many-body simulations. The cubic system has volume V = 2003 nm3 and surface area S = 2002 nm2. All simulations use the following parameters: σ = 1 nm, Δt = 0.1 µs, D¯ = 12.67 nm2/μs, and particle number NA = 200 (all in solution at the beginning of simulations). (a) Kinetics of IL. Simulation parameters: ka = 347.18 nm3/μs and kb = 2.09 × 103 s−1; the reflecting surface RS3D = 0.69 nm; and lipid number NL = 500 (initially all unbound at the beginning of simulations). The equilibrium state is KD = 2kb/ka = 20 µM, which is also defined as KD ≡ [A][L]/[AL] = nA[NL − (NAnA)]/[V(NAnA)]; thus, we have the number of free solution particles at equilibrium nA = 0.5[(NANLKDV) + (NANLKDV)2+4KDVNA] = 43.77. (b) Kinetics of IL and EL in diffusion-limited situations are also in excellent agreement: ka = 347.18 nm3/μs, kb = 2.09 s−1, and RS3D = 0.69 nm. (c) Kinetics of IL and EL in rate-limited situations: ka = 0.332 nm3/μs, kb = 1.00 s−1, and RS3D = 0.002 nm. In this figure, EL curves are shown as mean calculated by 50 trajectories, while IL curves are shown as mean ± Standard Deviation (SD) calculated also by 50 trajectories.

FIG. 5.

The implicit lipid (IL) and explicit lipid (EL) models with the planar surface shown in Fig. 4(a) are in very close agreement for many-body simulations. The cubic system has volume V = 2003 nm3 and surface area S = 2002 nm2. All simulations use the following parameters: σ = 1 nm, Δt = 0.1 µs, D¯ = 12.67 nm2/μs, and particle number NA = 200 (all in solution at the beginning of simulations). (a) Kinetics of IL. Simulation parameters: ka = 347.18 nm3/μs and kb = 2.09 × 103 s−1; the reflecting surface RS3D = 0.69 nm; and lipid number NL = 500 (initially all unbound at the beginning of simulations). The equilibrium state is KD = 2kb/ka = 20 µM, which is also defined as KD ≡ [A][L]/[AL] = nA[NL − (NAnA)]/[V(NAnA)]; thus, we have the number of free solution particles at equilibrium nA = 0.5[(NANLKDV) + (NANLKDV)2+4KDVNA] = 43.77. (b) Kinetics of IL and EL in diffusion-limited situations are also in excellent agreement: ka = 347.18 nm3/μs, kb = 2.09 s−1, and RS3D = 0.69 nm. (c) Kinetics of IL and EL in rate-limited situations: ka = 0.332 nm3/μs, kb = 1.00 s−1, and RS3D = 0.002 nm. In this figure, EL curves are shown as mean calculated by 50 trajectories, while IL curves are shown as mean ± Standard Deviation (SD) calculated also by 50 trajectories.

Close modal

The CPU time of these simulations is independent of the number of lipids, and thus as the density of lipids or the size of the surface area increases, the implicit lipid model provides huge speed-ups relative to the explicit lipid model (Fig. 6). The implicit model is actually slightly slower at the lowest lipid density; this is because most of the solution particles remain in solution and must continue evaluating binding to the surface throughout the simulation, which requires more evaluations than if the proteins are bound to the surface. This is also why the explicit lipid simulations do not show linear slowdowns with increasing lipid copies—as more particles are bound to the surface, the speed of a time step increases.

FIG. 6.

CPU time comparison for explicit lipid vs implicit lipid models. The implicit lipid model is independent of the number of lipids, whereas the explicit lipid model becomes slower with increasing lipids. Each simulation was run with 600 solution particles for 105 steps in a cubic system (V = 2003 nm3 and S = 2002 nm2), with the average CPU time per step being plotted vs lipid copies.

FIG. 6.

CPU time comparison for explicit lipid vs implicit lipid models. The implicit lipid model is independent of the number of lipids, whereas the explicit lipid model becomes slower with increasing lipids. Each simulation was run with 600 solution particles for 105 steps in a cubic system (V = 2003 nm3 and S = 2002 nm2), with the average CPU time per step being plotted vs lipid copies.

Close modal

For a standard surface adsorption model, the value of the adsorption coefficient κ does not change throughout the simulation.16 This neglects occupancy of surface sites by other solution particles. For high lipid density, this does not matter and adsorption models produce very similar kinetics to the explicit lipid simulations due to the excess of surface binding sites. However, when lipids are not in excess, the models diverge in both equilibrium and kinetics (Fig. 7). One adaptive approach to the adsorption model is to impose lipid density changes on it by dynamically setting κ = kon · ρL(t) (the macroscopic binding rate kon stays the same as in the explicit lipid method) to account for occupancy of lipid binding sites.

FIG. 7.

Explicit lipid (EL) method compared to the surface adsorption (AD) model using Ref. 16, where the adsorption rate is fixed. For limited lipids, the fixed adsorption model produces a distinct equilibrium (red curve) because it does not account for occupancy of lipid binding sites. With high lipid density, the models converge (blue curve). All simulations use the following parameters: σ = 1 nm and Δt = 0.1 µs; the cubic system with volume V = 2003 nm3 and surface area S = 2002 nm2; and particle number NA = 200 (all in solution at the beginning of simulations). For the explicit lipid model, we use D¯ = 12.67 nm2/μs (diffusion constant of a solution particle with D = 12 nm2/μs relative to a lipid on the membrane with D = 1 nm2/μs), ka = 347.18 nm3/μs, kb = 2.09 s−1, and thus the equilibrium KD = 2kb/ka = 0.02 µM. However for the adsorption model, we use D = 12 nm2/μs (diffusion constant of a single particle in solution not considering explicit lipids on the membrane). For the limited lipids case, ρL = 0.0025 nm−2 and for the sufficient lipids case, ρL = 0.05 nm−2.

FIG. 7.

Explicit lipid (EL) method compared to the surface adsorption (AD) model using Ref. 16, where the adsorption rate is fixed. For limited lipids, the fixed adsorption model produces a distinct equilibrium (red curve) because it does not account for occupancy of lipid binding sites. With high lipid density, the models converge (blue curve). All simulations use the following parameters: σ = 1 nm and Δt = 0.1 µs; the cubic system with volume V = 2003 nm3 and surface area S = 2002 nm2; and particle number NA = 200 (all in solution at the beginning of simulations). For the explicit lipid model, we use D¯ = 12.67 nm2/μs (diffusion constant of a solution particle with D = 12 nm2/μs relative to a lipid on the membrane with D = 1 nm2/μs), ka = 347.18 nm3/μs, kb = 2.09 s−1, and thus the equilibrium KD = 2kb/ka = 0.02 µM. However for the adsorption model, we use D = 12 nm2/μs (diffusion constant of a single particle in solution not considering explicit lipids on the membrane). For the limited lipids case, ρL = 0.0025 nm−2 and for the sufficient lipids case, ρL = 0.05 nm−2.

Close modal

Using the surface adsorption model described in Sec. II D 1 (the 1D Smoluchowski model) with adaptive κ produces the same equilibria as the explicit lipid model for all lipid densities (Fig. 8). The kinetics of the surface adsorption model are faster when the reaction is diffusion-limited (large kon), and lipids are of low density [Fig. 8(a)]. As lipid density (and thus κ) is increased, the models converge. For rate-limited reactions (small kon), both models have similar kinetics even for low lipid densities, which is again not too surprising, considering that these reactions are not highly sensitive to the specific particle distributions [Fig. 8(b)]. We note that κ is not unique to a single lipid density, and a lower density ρL with a stronger binding kon can produce the same κ. The models in Fig. 8(c) have similar kinetics but distinct equilibria due to the use of adaptive κ; otherwise they would be identical. A major shortcoming of surface adsorption models is that they can only describe binding from a particle diffusing above the membrane and thus colliding with it. Proteins that are localized to the membrane surface are not freely diffusing in the direction normal to the surface, and for these molecules to bind to additional lipids, we can apply our implicit lipid model, below.

FIG. 8.

Surface adsorption model using the 1D Smoluchowski model (AD-GF) compared to the explicit lipid (EL) method for the planar surface system [Fig. 11(a)]. Adsorption rate is adaptive with lipid density. All simulations use the following parameters: σ = 1 nm, Δt = 0.1 µs, and the cubic system with volume V = 2003 nm3 and surface area S = 2002 nm2, and particle number NA = 200 (all in solution at the beginning of simulations). For AD-GF, we use D = 12 nm2/μs (diffusion constant of the single particle in solution, no explicit lipids on the membrane), and for the EL model, D¯ = 12.67 nm2/μs. (a) Diffusion-limited situations. Parameters for AD-GF: ka = 173.59 nm3/μs, kb = 2.09 s−1, and the equilibrium KD = kb/ka = 0.02 µM. Parameters for EL, ka = 347.18 nm3/μs, kb = 2.09 s−1, and KD = 2kb/ka = 0.02 µM. (b) Rate-limited situations. Parameters for AD-GF: ka = 0.166 nm3/μs, kb = 1 s−1, and the equilibrium KD = kb/ka = 10 µM. Parameters for EL, ka = 0.332 nm3/μs, kb = 1 s−1, and KD = 2kb/ka = 10 µM. (c) Microscopic split of the adsorption coefficient κ(t) in the AD-GF model. Blue curve used ka = 173.59 nm3/μs, kb = 2.09 s−1, NL = 1000, ρL = NL/S = 0.025 nm−2, and thus kon = 40.35 nm3/μs [evaluated by Eq. (21)], which gives the adsorption coefficient κ = kon · ρL = 1.01 nm/μs. The red curve used ka = 8.53 nm3/μs, kb = 2.09 s−1, NL = 10 000, ρL = NL/S = 0.25 nm−2, and thus kon = 4.035 nm3/μs, which also gives the adsorption coefficient κ = kon · ρL = 1.01 nm/μs at the beginning of the simulations. However, their kinetics and equilibrium states are different due to adapting to lipid occupancy. In this figure, each EL curve is shown as mean calculated by 50 trajectories, while AD-GF curves are shown as mean ± SD calculated also by 50 trajectories, and SD is shown in the shaded region.

FIG. 8.

Surface adsorption model using the 1D Smoluchowski model (AD-GF) compared to the explicit lipid (EL) method for the planar surface system [Fig. 11(a)]. Adsorption rate is adaptive with lipid density. All simulations use the following parameters: σ = 1 nm, Δt = 0.1 µs, and the cubic system with volume V = 2003 nm3 and surface area S = 2002 nm2, and particle number NA = 200 (all in solution at the beginning of simulations). For AD-GF, we use D = 12 nm2/μs (diffusion constant of the single particle in solution, no explicit lipids on the membrane), and for the EL model, D¯ = 12.67 nm2/μs. (a) Diffusion-limited situations. Parameters for AD-GF: ka = 173.59 nm3/μs, kb = 2.09 s−1, and the equilibrium KD = kb/ka = 0.02 µM. Parameters for EL, ka = 347.18 nm3/μs, kb = 2.09 s−1, and KD = 2kb/ka = 0.02 µM. (b) Rate-limited situations. Parameters for AD-GF: ka = 0.166 nm3/μs, kb = 1 s−1, and the equilibrium KD = kb/ka = 10 µM. Parameters for EL, ka = 0.332 nm3/μs, kb = 1 s−1, and KD = 2kb/ka = 10 µM. (c) Microscopic split of the adsorption coefficient κ(t) in the AD-GF model. Blue curve used ka = 173.59 nm3/μs, kb = 2.09 s−1, NL = 1000, ρL = NL/S = 0.025 nm−2, and thus kon = 40.35 nm3/μs [evaluated by Eq. (21)], which gives the adsorption coefficient κ = kon · ρL = 1.01 nm/μs. The red curve used ka = 8.53 nm3/μs, kb = 2.09 s−1, NL = 10 000, ρL = NL/S = 0.25 nm−2, and thus kon = 4.035 nm3/μs, which also gives the adsorption coefficient κ = kon · ρL = 1.01 nm/μs at the beginning of the simulations. However, their kinetics and equilibrium states are different due to adapting to lipid occupancy. In this figure, each EL curve is shown as mean calculated by 50 trajectories, while AD-GF curves are shown as mean ± SD calculated also by 50 trajectories, and SD is shown in the shaded region.

Close modal

In addition to binding from solution to a surface, we demonstrate that the implicit lipid model can also accurately reproduce the ability of particles, previously localized to the surface through one interface, to perform 2D binding reactions with the surface using another interface. In Fig. 9, we compare the results of implicit lipid model 2D with the explicit lipid method. For all the systems, whether diffusion-limited (ka = 86.80 nm2/μs), rate-limited (ka = 0.083 nm2/μs), or diffusion influenced (ka = 8.68 nm2/μs), the implicit lipid model produces the proper equilibrium states due to enforcing detailed balance through the RS2D parameter. For these 2D reactions, we find that the kinetics of the implicit lipid model are slower than the explicit lipids, particularly for diffusion-limited reactions. This is likely due to the fact that in 2D, reaction dynamics are always sensitive to the spatial distribution of particles and are not characterized by a single macroscopic kon.53 By replacing individual particles with a field, we eliminate rapid reactions due to pairs of particles that are close together, the effect of which is exaggerated in 2D relative to 3D, and by assigning RS2Dσ to ensure equilibrium, we also slow the reaction kinetics. RS2D here increases from σ at the smallest ka value, up to 1.46σ for the largest ka (see Fig. S1).

FIG. 9.

The 2D implicit lipid (IL) model reproduces the proper equilibrium as the explicit lipid (EL) model, but the kinetics are slower for diffusion-limited systems. The system used the flat surface as shown in Fig. 4(a). All simulations used the following parameters: σ = 1 nm, D = 1 nm2/μs, Δt = 0.1 µs, particle number NA = 100, and lipids number NL = 1000 (all are free at the beginning of simulations). Three different simulations are carried out. The strongly diffusion influenced black curve used ka = 86.80 nm2/μs, kb = 2.09 s−1, the equilibrium KD = kb/ka = 2.41 × 10−2µm−2, surface area S = 2002 nm2, and the reflecting distance RS2D = 1.46 nm. The diffusion-influenced blue curve used ka = 8.68 nm2/μs, kb = 2.09 × 103 s−1, KD = kb/ka = 2.41 × 102µm−2, S = 7002 nm2, and RS2D = 1.35 nm. The rate-limited red curve used ka = 0.083 nm2/μs, kb = 1 s−1, KD = kb/ka = 12.05 µm−2, S = 2002 nm2, and RS2D = 1.00 nm. In this figure, each EL curve is shown as mean calculated by 50 trajectories, while IL curves are shown as mean ± SD calculated also by 50 trajectories, and SD is shown in the shaded region.

FIG. 9.

The 2D implicit lipid (IL) model reproduces the proper equilibrium as the explicit lipid (EL) model, but the kinetics are slower for diffusion-limited systems. The system used the flat surface as shown in Fig. 4(a). All simulations used the following parameters: σ = 1 nm, D = 1 nm2/μs, Δt = 0.1 µs, particle number NA = 100, and lipids number NL = 1000 (all are free at the beginning of simulations). Three different simulations are carried out. The strongly diffusion influenced black curve used ka = 86.80 nm2/μs, kb = 2.09 s−1, the equilibrium KD = kb/ka = 2.41 × 10−2µm−2, surface area S = 2002 nm2, and the reflecting distance RS2D = 1.46 nm. The diffusion-influenced blue curve used ka = 8.68 nm2/μs, kb = 2.09 × 103 s−1, KD = kb/ka = 2.41 × 102µm−2, S = 7002 nm2, and RS2D = 1.35 nm. The rate-limited red curve used ka = 0.083 nm2/μs, kb = 1 s−1, KD = kb/ka = 12.05 µm−2, S = 2002 nm2, and RS2D = 1.00 nm. In this figure, each EL curve is shown as mean calculated by 50 trajectories, while IL curves are shown as mean ± SD calculated also by 50 trajectories, and SD is shown in the shaded region.

Close modal

In Fig. 10, we again find excellent agreement between the implicit and explicit lipid simulations, but here with binding to curved surfaces [Figs. 4(c) and 4(d)]. As we noted above, the implicit lipid reaction probabilities are analytical for a sphere but must be numerically integrated over the separation between the solution molecule and the surface for the triangulated mesh surface of Fig. 4(d). For the explicit lipid simulations, the lipids must also perform diffusion on a curved surface. The implicit lipid model produces exact equilibrium populations and excellent agreement with the kinetics of the explicit lipid method, whether in spherical systems [Figs. 10(a)–10(d)], in a cylindrical tube (Fig. S4), or with an irregular curved surface [Fig. 10(e)].

FIG. 10.

For binding to curved surfaces, the implicit lipid (IL) model shows very close agreement with the explicit lipid (EL) model results, with curved surfaces shown in Figs. 4(c) and 4(d). All simulations used the following parameters: σ = 1 nm, D¯ = 12.67 nm2/μs, Δt = 0.1 µs, ka = 347.18 nm3/μs, kb = 2.09 × 103 s−1, and thus the equilibria KD = 2kb/ka = 20 µM. (a) Spherical surface with the radius R = 75 nm. The solution particle number NA = 21, lipids on the surface NL = 112. Particles diffuse inside the sphere or diffuse while bound onto the sphere surface but cannot move outside the sphere. At the beginning of the simulations, all solution particles and lipids are unbound. (b) Spherical surface, R = 100 nm, NA = 50, and NL = 199. (c) Spherical surface, R = 150 nm, NA = 168, and NL = 448. (d) Spherical surface, R = 250 nm, NA = 776, and NL = 1244. (e) Irregular surface. Particles are closed inside the surface and diffuse within the enclosed volume or on while bound to the surface. The total volume accessible for solution particles, is V = 2.19 × 107 nm3 and total surface area S = 6.02 × 105 nm2. Solution particle number NA = 252 and lipids on the surface NL = 934. Each EL curve is the mean calculated by 100 trajectories, and each IL curve is the mean from 100 trajectories.

FIG. 10.

For binding to curved surfaces, the implicit lipid (IL) model shows very close agreement with the explicit lipid (EL) model results, with curved surfaces shown in Figs. 4(c) and 4(d). All simulations used the following parameters: σ = 1 nm, D¯ = 12.67 nm2/μs, Δt = 0.1 µs, ka = 347.18 nm3/μs, kb = 2.09 × 103 s−1, and thus the equilibria KD = 2kb/ka = 20 µM. (a) Spherical surface with the radius R = 75 nm. The solution particle number NA = 21, lipids on the surface NL = 112. Particles diffuse inside the sphere or diffuse while bound onto the sphere surface but cannot move outside the sphere. At the beginning of the simulations, all solution particles and lipids are unbound. (b) Spherical surface, R = 100 nm, NA = 50, and NL = 199. (c) Spherical surface, R = 150 nm, NA = 168, and NL = 448. (d) Spherical surface, R = 250 nm, NA = 776, and NL = 1244. (e) Irregular surface. Particles are closed inside the surface and diffuse within the enclosed volume or on while bound to the surface. The total volume accessible for solution particles, is V = 2.19 × 107 nm3 and total surface area S = 6.02 × 105 nm2. Solution particle number NA = 252 and lipids on the surface NL = 934. Each EL curve is the mean calculated by 100 trajectories, and each IL curve is the mean from 100 trajectories.

Close modal

Finally, in Fig. 11, we illustrate how the implicit lipid approximation also applies to surfaces with a nonuniform distribution of lipids, which could occur due to lipid phase separation or lipid microdomains. Here, the explicit lipid model had individual lipids restricted to one patch of the surface, and the implicit lipid model thus had the density also restricted to a patch of the surface [Fig. 4(b)]. This model is also an important demonstration of how this method can be coupled to a hybrid single-particle/PDE based approach, with lipid densities that can vary along the surface. The implicit lipid model shows excellent agreement with the explicit lipid model for both rate-limited and diffusion limited reactions.

FIG. 11.

The implicit lipid (IL) model agrees very well with the explicit lipid (EL) model for an inhomogeneous lipid distribution generated by restricting lipids to half the surface, as shown in Fig. 4(b). The box system has a volume V = 2003 nm3 and surface area S = 2002/2 nm2. All simulations use the following parameters: σ = 1 nm, D¯ = 12.67 nm2/μs, Δt = 0.1 µs, solution particle number NA = 100, and lipids on the surface NL = 1000. The blue and pink curves used ka = 347.18 nm3/μs, kb = 2.09 s−1, KD = 2kb/ka = 0.02 µM, and RS3D = 0.69 nm. The red and green curves used ka = 0.332 nm3/μs, kb = 1 s−1, KD = 2kb/ka = 10 µM, and RS3D = 0.002 nm. Each line is mean ± SD from 50 trajectories.

FIG. 11.

The implicit lipid (IL) model agrees very well with the explicit lipid (EL) model for an inhomogeneous lipid distribution generated by restricting lipids to half the surface, as shown in Fig. 4(b). The box system has a volume V = 2003 nm3 and surface area S = 2002/2 nm2. All simulations use the following parameters: σ = 1 nm, D¯ = 12.67 nm2/μs, Δt = 0.1 µs, solution particle number NA = 100, and lipids on the surface NL = 1000. The blue and pink curves used ka = 347.18 nm3/μs, kb = 2.09 s−1, KD = 2kb/ka = 0.02 µM, and RS3D = 0.69 nm. The red and green curves used ka = 0.332 nm3/μs, kb = 1 s−1, KD = 2kb/ka = 10 µM, and RS3D = 0.002 nm. Each line is mean ± SD from 50 trajectories.

Close modal

Single-particle RD methods are a powerful tool for studying cell-scale dynamics, and by integrating them with increasingly advanced membrane models, we push the development of cell-scale membrane remodeling simulations. We here introduced an implicit lipid RD model for reversible binding of particles to membranes of varying lipid density and curvature, demonstrating that we could quite accurately reproduce the equilibrium and kinetics of explicit lipid simulations at speed increases of orders of magnitude. This method enables simulations of proteins interacting with large and arbitrarily shaped membranes without the considerable expense of propagating lipid binding sites. Unlike models of surface adsorption, we developed our method to integrate over all possible collisions with the surface even for large time steps, thus ensuring and demonstrating its ability to apply to curved or inhomogeneous membranes, with excellent agreement to explicit lipid binding kinetics. Uniquely, our method also allows for binding between proteins on the surface and additional lipids (a 2D reaction), which represents an important stabilizing contact for proteins that are initially localized to the membrane via a protein binding partner.7 Our RD code is open source, incorporated into a new software tool to facilitate user-defined models, and immediately downloadable from Github at https://github.com/mjohn218/fpr_implicit, including both the implicit lipid model and the surface adsorption model for comparison.

To ensure that our implicit lipid model here is transferrable and comparable to other rate-based tools, whether other single-particle,11,23,24,33,35 RDME,29–31 or ODE/PDE based tools,25,26 and quantitatively comparable to experimental studies of binding,55 we have taken care to quantify how our model parameters translate to macroscopic rates and equilibrium constants. We additionally characterized how surface adsorption rates, which have distinct units relative to bimolecular association rates, can be adapted to, if necessary, better mimic the binding and kinetics of explicit lipid simulations. A limitation of our approach is that the implicit lipid model, unlike binding models between explicit particle pairs,56,57 is not a formally exact solution to a well-defined dynamical system. The particles are propagated according to free diffusion, which is algorithmically much simpler and transferrable to other RD tools, but neglects the effect of the reactive lipids on particle positions. We correct for this approximation by introducing the RS distance in both 3D and 2D to ensure detailed balance between binding and unbinding is maintained. Remarkably, this model nonetheless produces kinetics that are in very close agreement with explicit lipid models, indicating the robustness of replacing individual particle interactions with their spatially integrated average. An exception is for purely 2D reactions, where the implicit lipid model has slower kinetics for diffusion-limited reactions, likely due to the sensitivity of reaction dynamics in 2D to spatial distributions.53,58

A limitation of our theoretical approach is that for binding to curved surfaces, we approximate one aspect of the model, the pairwise reaction probability [Eq. (2)], as unaffected by curvature. We do account for the effect of curvature on determining separation from the membrane and thus the integrated reaction probability [Eq. (11)], on placement of particles that displace beyond the membrane, and on the macroscopic rate [Eq. (21)]. We show in the supplementary material that for a single step, membrane curvature has only a small influence on the pairwise reaction probability (supplementary material). In our simulations of binding on curved surfaces, our results for both implicit and explicit-lipid models produced excellent agreement with the theoretically predicted equilibrium, further indicating that these curvature effects are quite small, even for a sphere of R = 75 nm. Reducing the time step reduces the scale of the local curvature, thus improving the assumption of a planar reflecting surface on pairwise reaction probabilities and providing a simple mechanism for evaluating the impact of curvature on binding. We found that using the reaction probabilities from the model of Eq. (1) worked even for our largest time step of 1 μs, and the probabilities were independent of time step, as expected for a convergent method (Fig. S5). We acknowledge that our algorithm for arbitrary curved surfaces does require numerical integration for determining reaction probabilities, which is more expensive than simple adsorption models, but is nonetheless straightforward to manage using widely existing libraries such as GSL. We also note that one advantage of surface adsorption methods is that they can more easily introduce forces between planar surfaces and solution particles18 because they use solutions to the 1D GF rather than the 3D version we use here.

Our implicit lipid model can be combined with complex networks of protein-protein interactions to efficiently simulate a range of systems involving reversible localization to the membrane, such as models of clathrin-cage assembly7,35,59,60 and cell polarization,37 or experiments involving binding and oligomerization on membranes.61 Integration with continuum models of surfaces, as we have done here, is especially critical for developing quantitative and dynamical models of membrane dynamics as they are driven by proteins.62–65 These models are key for understanding processes from endocytosis, exocytosis, and virion formation at the smaller scale64–66 to cell division at the larger scale. These current models lack single-particle resolution of the proteins localizing from solution to the membrane, as we do here, which is a necessary step in capturing more molecular-level dynamics over cellular length scales. We focused here on protein binding to a single lipid-type population, and an important future work would be to generalize the implicit lipid model to capture binding to membranes with multiple lipid populations with varying affinity. In addition, although we have focused here on biological problems involving membranes, our rate-based approach applies equally well to study any type of reaction-dynamics involving surfaces, such as heterogeneous catalysis.67,68 Single-particle RD is more expensive than PDE or RDME simulations of reaction dynamics, but it is unique in providing the structure to species that can be used to capture short-range spatial correlations and structure-resolved self-assembly.35 While RD methods in general lack the physicochemical details of the molecular dynamics approaches,8,9 single-particle RD at least has the capacity to build in deterministic forces between species,14,36 and RD dramatically exceeds molecular dynamics in reaching cellular length and time scales, with applications to equilibrium and nonequilibrium systems. Our method here provides an essential component for performing tractable simulations involving large networks of proteins that localize dynamically to membranes, to ultimately build predictive models for cell biology.

Simulations of the explicit lipid model in 3D and 2D were implemented here with the FPR algorithm.14,35,53 Initial configurations were produced by a randomly sampled uniform distribution of particles in solution or on the surface. The binding probability of a solution particle to a lipid particle at a distance r < Rmax is determined by the pairwise association probability preactt|r) [Eq. (2)] multiplied with time-dependent reweighting ratios determined from the 3D GF.14 In 2D, the same procedure applies, with the 2D GFs used. A uniformly distributed random number between 0 and 1 is generated and compared to the binding probability. If the random number is smaller than the binding probability, then the binding reaction occurs, otherwise the particle survives from this reaction and needs to check whether to bind to other lipid particles. If the binding takes place, the particle is put adjacent to the lipid particle on the membrane surface at a separation of σ. If a particle does not bind, then it is propagated according to the free diffusion, where its new position is rejected and resampled if the particle overlaps (r < σ) any reaction partners. Particles that cross the membrane or boundaries are reflected back inside the box. For bound species, the dissociation probability is evaluated by Eq. (3). If dissociation occurs, the particle becomes free and is left at contact with its former partner, at r = σ. Only in the next time step, it is allowed to attempt a new reaction or diffuse.

For implicit lipid simulations with a flat membrane surface, if a solution particle is within Rmax of the surface, based on the shortest distance to the surface, then the probability of binding will be evaluated by Eq. (6). In addition, a uniformly distributed random number between 0 and 1 is generated and compared to the binding probability. If the random number is smaller than the binding probability, then the binding reaction occurs, otherwise the particle survives. If the reaction occurs, then the solution particle is placed on the membrane surface at the shortest distance. If the particle survives from binding, then it is freely diffused, where now the reflection occurs across the plane located at RS3D above the membrane. For bound species, the probability of dissociation is calculated by Eq. (7). If the dissociation occurs, then the particle will be put at the closest point on the reflecting surface and will be treated as a free particle in the next time step. Once a binding or unbinding event occurs, the number of free lipids is updated such that the density ρL(t) of free lipids is correct at all times.

For 3D cases with the curved membrane surface, the simulation details are similar, except that the binding probability needs to be integrated over the accessible area on the surface [for spheres we can use the closed-form expression Eq. (S9)]. We use a triangular mesh to build up the curved surface,43 and thus the accessible area is composed of many triangular elements. Each element contributes to the binding probability by pairwise reaction probability preactt|ri) multiplying element area si and the local density of lipids ρi. Then, the total binding probability is accumulated over all these elements based on Eq. (11), ri<Rmaxpreact(Δt|ri)siρi. For curved membranes with analytical shapes such as the sphere, the reaction probability can be calculated analytically for each height z0 (supplementary material).

For simulations of explicit protein particles restricted to the 2D surface with implicit lipids, simulations are similar to 3D cases, but Dz = 0 and distances are based only on the 2D separation between points. When protein particles bind to the lipids, their position does not need to change, but their diffusion coefficient is updated due to the formation of a complex and they are no-longer free to bind. Dissociation then reverses these changes, without requiring specific placement of the protein.

For simulations of surface adsorption using the 1D GF of Sec. II D, the binding probability of a solution particle is evaluated by Eq. (25). If the binding occurs, the particle will be put at the nearest point on the surface z = σ (i.e., radiation boundary condition). The position for this 1D approach after a survival would require a free diffusion move in x and y and a movement in z sampled from the pirr solution shown in Eq. (24), for a new z value by solving RAND=σzpirr(z,Δt|z0)dz, where RAND is a uniformly distributed random number between 0 and the value of the survival probability 1 − PAD-GFt|z0). This sampling ensures an exact solution to the model defined in Sec. II D 1 and is the same approach used in Green’s function reaction dynamics (GFRD).23 For the bound particle, the dissociation probability is evaluated by Eq. (3). If the dissociation occurs, the particle becomes free but remains in contact with the surface at z = σ. In the next time step, it can then either react or diffuse.

A subdivision limit surface (SLS) mesh algorithm described in Refs. 69 and 70 was used for those complex surfaces not representable by a closed form. Briefly, the SLS is determined by subdividing each triangle of a mesh into four smaller triangles repeatedly. Each iteration of the subdivision creates a finer and finer mesh. A rule is applied to generate the vertices of the subdivided triangles from the previous iteration. That rule is codified as a matrix transformation of the previous iteration’s vector of vertex coordinates that yields the coordinates of the new vertices. Taking the limit of infinite applications of the rule yields a set of points constituting a surface that is tangent plane continuous. In practice, the limit surface itself is written as a polynomial spline that is proportional to the original vertices. The spline can be computed from the mesh that results following any amount of subdivision.

At any iteration of the subdivision, the spline of the limiting surface is contained completely within the convex hull defined by the spline control points (that is, the vertex coordinates of that subdivision). Therefore, if the convex hull of a second object does not intersect the convex hull of the control points, it does not intersect the limit surface (Fig. S6A). If an object does intersect the limit surface, it will intersect the convex hull of the mesh at any level of subdivision (Fig. S6B). Therefore, collisions of any convex object with the mesh can be excluded to arbitrary thresholds by checking for intersections between convex hulls at a particular depth of subdivision. Here, the Gilbert-Johnson-Keerthi algorithm44 was used to detect convex hull collisions. A recursive function was defined to detect collisions with the target object and the box spline of a mesh element. If a collision was detected at a particular subdivision level, the recursive function was called on each of the four subdivisions of the triangle. If collision is excluded, the recursive branch terminates (Fig. S7). This recursion proceeds until a particular depth level (ten) at which point a collision was declared. This basic framework allows for computing the point of the closest approach of a particle to the surface (by running bisection to determine the largest sphere that does not contact the surface) as well as detecting the point of intersection of a linear path through the surface (using that line segment for collision detection) (Fig. S8).

Convenient surface coordinates for a face of the SLS are its (nonorthogonal) spline parameters, here denoted as u and v. These take the place of x and y that might be used on a planar surface, where the Laplacian is appropriate for describing the time dependence of diffusion for x and y coordinates on a Cartesian plane, diffusion on the SLS must be propagated using the Laplace-Beltrami operator as done in Ref. 46. This corrects for both the units and nonorthogonality of the spline coordinate.

See the supplementary material for two parts of theoretical proofs—I: pairwise reaction probability influenced by curved surfaces and II: implicit lipid model on curved surfaces—and 8 supplementary figures.

Research reported in this publication was supported by the National Science Foundation CAREER Award to M.E.J. (Award No. 1753174). The research used the NSF XSEDE computational resources of SuperMic under Grant No. XRAC MCB150059 and the MARCC supercomputer at JHU. Funding to A.J.S., R.K., and G.T. was provided by the intramural research program of the Eunice Kennedy Shriver National Institute of Child Health and Human Development of the NIH.

1.
S.
McLaughlin
 et al, “
PIP(2) and proteins: Interactions, organization, and information flow
,”
Annu. Rev. Biophys. Biomol. Struct.
31
,
151
175
(
2002
).
2.
J. W.
Yu
 et al, “
Genome-wide analysis of membrane targeting by S. cerevisiae pleckstrin homology domains
,”
Mol. Cell
13
(
5
),
677
688
(
2004
).
3.
A.
Frost
,
V. M.
Unger
, and
P.
De Camilli
, “
The BAR domain superfamily: Membrane-molding macromolecules
,”
Cell
137
(
2
),
191
196
(
2009
).
4.
Y.
Yoon
,
X.
Zhang
, and
W.
Cho
, “
Phosphatidylinositol 4,5-bisphosphate (PtdIns(4,5)P2) specifically induces membrane penetration and deformation by Bin/amphiphysin/Rvs (BAR) domains
,”
J. Biol. Chem.
287
(
41
),
34078
34090
(
2012
).
5.
S. M.
Christensen
 et al, “
One-way membrane trafficking of SOS in receptor-triggered Ras activation
,”
Nat. Struct. Mol. Biol.
23
(
9
),
838
846
(
2016
).
6.
B. N.
Kholodenko
,
J. B.
Hoek
, and
H. V.
Westerhoff
, “
Why cytoplasmic signalling proteins should be recruited to cell membranes
,”
Trends Cell Biol.
10
(
5
),
173
178
(
2000
).
7.
O. N.
Yogurtcu
and
M. E.
Johnson
, “
Cytoplasmic proteins can exploit membrane localization to trigger functional assembly
,”
PLoS Comput. Biol.
14
(
3
),
e1006031
(
2018
).
8.
E.
Lindahl
and
M. S.
Sansom
, “
Membrane proteins: Molecular dynamics simulations
,”
Curr. Opin. Struct. Biol.
18
(
4
),
425
431
(
2008
).
9.
I.
Yu
,
M.
Feig
, and
Y.
Sugita
, “
High-performance data analysis on the big trajectory data of cellular scale all-atom molecular dynamics simulations
,”
J. Phys.: Conf. Ser.
1036
,
012009
(
2018
).
10.
D. T.
Gillespie
,
E.
Seitaridou
, and
C. A.
Gillespie
, “
The small-voxel tracking algorithm for simulating chemical reactions among diffusing molecules
,”
J. Chem. Phys.
141
(
23
),
234115
(
2014
).
11.
J.
Schoneberg
and
F.
Noe
, “
ReaDDy—A software for particle-based reaction-diffusion dynamics in crowded cellular environments
,”
PLoS One
8
(
9
),
e74261
(
2013
).
12.
R. A.
Kerr
 et al, “
Fast Monte Carlo simulation methods for biological reaction-diffusion systems in solution and on surfaces
,”
SIAM J. Sci. Comput.
30
(
6
),
3126
3149
(
2008
).
13.
J. S.
van Zon
and
P. R.
ten Wolde
, “
Green’s-function reaction dynamics: A particle-based approach for simulating biochemical networks in space and time
,”
J. Chem. Phys.
123
,
234910
(
2005
).
14.
M. E.
Johnson
and
G.
Hummer
, “
Free-propagator reweighting integrator for single-particle dynamics in reaction-diffusion models of heterogeneous protein-protein interaction systems
,”
Phys. Rev. X
4
(
3
),
031037
(
2014
).
15.
R. V.
Stahelin
 et al, “
Contrasting membrane interaction mechanisms of AP180 N-terminal homology (ANTH) and epsin N-terminal homology (ENTH) domains
,”
J. Biol. Chem.
278
(
31
),
28993
28999
(
2003
).
16.
S. S.
Andrews
, “
Accurate particle-based simulation of adsorption, desorption and partial transmission
,”
Phys. Biol.
6
(
4
),
046015
(
2009
).
17.
I.
Langmuir
, “
The adsorption of gases on plane surfaces of glass, mica and platinum
,”
J. Am. Chem. Soc.
40
(
9
),
1361
1403
(
1918
).
18.
G.
Lamm
and
K.
Schulten
, “
Extended Brownian dynamics. II. Reactive, non-linear diffusion
,”
J. Chem. Phys.
78
(
5
),
2713
2734
(
1983
).
19.
G.
Lamm
and
K.
Schulten
, “
Extended Brownian dynamics approach to diffusion-controlled processes
,”
J. Chem. Phys.
75
(
1
),
365
371
(
1981
).
20.
A.
Singer
 et al, “
Partially reflected diffusion
,”
SIAM J. Appl. Math.
68
(
3
),
844
868
(
2008
).
21.
R.
Erban
and
S. J.
Chapman
, “
Reactive boundary conditions for stochastic simulations of reaction-diffusion processes
,”
Phys. Biol.
4
(
1
),
16
28
(
2007
).
22.
T. M.
Bartol
 et al, “
Monte-Carlo simulation of miniature end-plate current generation in the vertebrate neuromuscular-junction
,”
Biophys. J.
59
(
6
),
1290
1307
(
1991
).
23.
T. R.
Sokolowski
 et al, “
eGFRD in all dimensions
,”
J. Chem. Phys.
150
(
5
),
054108
(
2019
).
24.
S. S.
Andrews
, “
Smoldyn: Particle-based simulation with rule-based modeling, improved molecular interaction and a library interface
,”
Bioinformatics
33
(
5
),
710
717
(
2017
).
25.
I. I.
Moraru
 et al, “
Virtual cell modelling and simulation software environment
,”
IET Syst. Biol.
2
(
5
),
352
362
(
2008
).
26.
B. R.
Angermann
 et al, “
Computational modeling of cellular signaling processes embedded into dynamic spatial contexts
,”
Nat. Methods
9
(
3
),
283
289
(
2012
).
27.
Y.
Yoon
 et al, “
In situ quantitative imaging of cellular lipids using molecular sensors
,”
Nat. Chem.
3
(
11
),
868
874
(
2011
).
28.
M.
Malek-Mansour
and
J.
Houard
, “
A new approximation scheme for the study of fluctuations in nonuniform nonequilibrium systems
,”
Phys. Lett. A
70
,
366
368
(
1979
).
29.
B.
Drawert
,
S.
Engblom
, and
A.
Hellander
, “
URDME: A modular framework for stochastic simulation of reaction-transport processes in complex geometries
,”
BMC Syst. Biol.
6
,
76
(
2012
).
30.
E.
Roberts
,
J. E.
Stone
, and
Z.
Luthey-Schulten
, “
Lattice microbes: High-performance stochastic simulation method for the reaction-diffusion master equation
,”
J. Comput. Chem.
34
(
3
),
245
255
(
2013
).
31.
B.
Drawert
 et al, “
A framework for discrete stochastic simulation on 3D moving boundary domains
,”
J. Chem. Phys.
145
(
18
),
184113
(
2016
).
32.
J.
Schoneberg
,
A.
Ullrich
, and
F.
Noe
, “
Simulation tools for particle-based reaction-diffusion dynamics in continuous space
,”
BMC Biophys.
7
,
11
(
2014
).
33.
G.
Gruenert
 et al, “
Rule-based spatial modeling with diffusing, geometrically constrained molecules
,”
BMC Bioinf.
11
,
307
(
2010
).
34.
S. S.
Andrews
 et al, “
Detailed simulations of cell biology with Smoldyn 2.1
,”
PLoS Comput. Biol.
6
,
e1000705
(
2010
).
35.
M. E.
Johnson
, “
Modeling the self-assembly of protein complexes through a rigid-body rotational reaction-diffusion algorithm
,”
J. Phys. Chem. B
122
(
49
),
11771
11783
(
2018
).
36.
A.
Vijaykumar
 et al, “
Multiscale simulations of anisotropic particles combining molecular dynamics and Green’s function reaction dynamics
,”
J. Chem. Phys.
146
(
11
),
114106
(
2017
).
37.
W. X.
Chew
 et al, “
Surface reaction-diffusion kinetics on lattice at the microscopic scale
,”
Phys. Rev. E
99
(
4
),
042411
(
2019
).
38.
L.
Batsilas
,
A. M.
Berezhkovskii
, and
S. Y.
Shvartsman
, “
Stochastic model of autocrine and paracrine signals in cell culture assays
,”
Biophys. J.
85
(
6
),
3659
3665
(
2003
).
39.
A. M.
Berezhkovskii
 et al, “
Boundary homogenization for trapping by patchy surfaces
,”
J. Chem. Phys.
121
(
22
),
11390
11394
(
2004
).
40.
M. B.
Flegg
,
S. J.
Chapman
, and
R.
Erban
, “
The two-regime method for optimizing stochastic reaction-diffusion simulations
,”
J. R. Soc., Interface
9
(
70
),
859
868
(
2012
).
41.
C. A.
Smith
and
C. A.
Yates
, “
The auxiliary region method: A hybrid method for coupling PDE- and Brownian-based dynamics for reaction-diffusion systems
,”
R. Soc. Open Sci.
5
(
8
),
180920
(
2018
).
42.
C. A.
Smith
and
C. A.
Yates
, “
Spatially extended hybrid methods: A review
,”
J. R. Soc., Interface
15
(
139
),
20170931
(
2018
).
43.
F.
Feng
and
W. S.
Klug
, “
Finite element modeling of lipid bilayer membranes
,”
J. Comput. Phys.
220
(
1
),
394
408
(
2006
).
44.
E. G.
Gilbert
,
D. W.
Johnson
, and
S. S.
Keerthi
, “
A fast procedure for computing the distance between complex objects in three-dimensional space
,”
IEEE J. Rob. Autom.
4
(
2
),
193
203
(
1988
).
45.
M.
Lin
and
S.
Gottschalk
, “
Collision detection between geometric models: A survey
,” in
Proceedings of IMA Conference on Mathematics of Surfaces
(
IMA
,
1998
), Vol. 8, pp.
602
608
.
46.
E.
Reister-Gottfried
,
S. M.
Leitenberger
, and
U.
Seifert
, “
Hybrid simulations of lateral diffusion in fluctuating membranes
,”
Phys. Rev. E
75
(
1
),
011908
(
2007
).
47.
W. T.
Snead
 et al, “
BAR scaffolds drive membrane fission by crowding disordered domains
,”
J. Cell Biol.
218
(
2
),
664
682
(
2019
).
48.
D. J.
Busch
 et al, “
Intrinsically disordered proteins drive membrane curvature
,”
Nat. Commun.
6
,
7875
(
2015
).
49.
H. S.
Carslaw
and
J. C.
Jaeger
,
Conduction of Heat in Solids
(
Clarendon Press
,
Oxford
,
1959
).
50.
N.
Agmon
and
A.
Szabo
, “
Theory of reversible diffusion-influenced reactions
,”
J. Chem. Phys.
92
(
9
),
5270
5284
(
1990
).
51.
H. X.
Zhou
, “
Kinetics of diffusion-influenced reactions studied by Brownian dynamics
,”
J. Phys. Chem.
94
(
25
),
8794
8800
(
1990
).
52.
H. X.
Zhou
, “
Brownian dynamics study of the influences of electrostatic interaction and diffusion on protein-protein association kinetics
,”
Biophys. J.
64
(
6
),
1711
1726
(
1993
).
53.
O. N.
Yogurtcu
and
M. E.
Johnson
, “
Theory of bi-molecular association dynamics in 2D for accurate model and experimental parameterization of binding rates
,”
J. Chem. Phys.
143
,
084117
(
2015
).
54.
T.
Baumgart
 et al, “
Large-scale fluid/fluid phase separation of proteins and lipids in giant plasma membrane vesicles
,”
Proc. Natl. Acad. Sci. U. S. A.
104
(
9
),
3165
(
2007
).
55.
W. W.
Cho
,
L.
Bittova
, and
R. V.
Stahelin
, “
Membrane binding assays for peripheral proteins
,”
Anal. Biochem.
296
(
2
),
153
161
(
2001
).
56.
M.
von Smoluchowski
, “
Attempt to derive a mathematical theory of coagulation kinetics in colloidal solutions
,”
Z. Phys. Chem.
92
,
129
(
1917
).
57.
F. C.
Collins
and
G. E.
Kimball
, “
Diffusion-controlled reaction rates
,”
J. Colloid Sci.
4
(
4
),
425
437
(
1949
).
58.
U. C.
Tauber
,
M.
Howard
, and
B. P.
Vollmayr-Lee
, “
Applications of field-theoretic renormalization group methods to reaction-diffusion problems
,”
J. Phys. A: Math. Gen.
38
(
17
),
R79
R131
(
2005
).
59.
A. P.
Schoen
 et al, “
Dynamic remodelling of disordered protein aggregates is an alternative pathway to achieve robust self-assembly of nanostructures
,”
Soft Matter
9
(
38
),
9137
9145
(
2013
).
60.
J. J.
VanDersarl
 et al, “
Rheology and simulation of 2-dimensional clathrin protein network assembly
,”
Soft Matter
10
(
33
),
6219
6227
(
2014
).
61.
W. F.
Zeno
 et al, “
Molecular mechanisms of membrane curvature sensing by a disordered protein
,”
J. Am. Chem. Soc.
141
,
10361
(
2019
).
62.
K.
Sapp
and
L.
Maibaum
, “
Suppressing membrane height fluctuations leads to a membrane-mediated interaction among proteins
,”
Phys. Rev. E
94
(
5
),
052414
(
2016
).
63.
K.
Sapp
,
A. J.
Sodt
, and
L.
Maibaum
, “
Modeling relaxation timescales of coupled membrane/protein systems
,”
Biophys. J.
116
(
3
),
363a
(
2019
).
64.
J. E.
Hassinger
 et al, “
Design principles of robust vesiculation in clathrin-mediated endocytosis
,”
Proc. Nat. Acad. Sci. U.S.A.
114
(
7
),
E1118
E1127
(
2017
).
65.
N.
Walani
,
J.
Torres
, and
A.
Agrawal
, “
Endocytic proteins drive vesicle growth via instability in high membrane tension environment
,”
Proc. Natl. Acad. Sci. U. S. A.
112
(
12
),
E1423
E1432
(
2015
).
66.
N.
Cordella
 et al, “
Membrane indentation triggers clathrin lattice reorganization and fluidization
,”
Soft Matter
11
(
3
),
439
448
(
2015
).
67.
K.
Reuter
and
M.
Scheffler
, “
First-principles kinetic Monte Carlo simulations for heterogeneous catalysis: Application to the CO oxidation at RuO2(110)
,”
Phys. Rev. B
73
(
4
),
045433
(
2006
).
68.
A.
Chatterjee
and
D. G.
Vlachos
, “
An overview of spatial microscopic and accelerated kinetic Monte Carlo methods
,”
J. Comput.-Aided Mater. Des.
14
(
2
),
253
308
(
2007
).
69.
F.
Cirak
,
M.
Ortiz
, and
P.
Schröder
, “
Subdivision surfaces: A new paradigm for thin-shell finite-element analysis
,”
Int. J. Numer. Methods Eng.
47
(
12
),
2039
2072
(
2000
).
70.
C.
Loop
, “
Smooth subdivision surfaces based on triangles
,” M.Sc. thesis,
Department of Mathematics, University of Utah
,
1987
.

Supplementary Material