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.

## I. INTRODUCTION

Proteins localize to membranes via multiple different binding modes, including recognition and binding to highly specific lipid head-groups [e.g., PI(4,5)P_{2}],^{1} electrostatically driven adherence to negatively charged membranes^{2} (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 constants^{15} 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 *k*_{on} 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 adsorption^{16,18–23} and Langmuir^{7,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)P_{2} is 25 000 *µ*m^{−2},^{27} translating to about 50–80 *µ*M given a yeast cell surface area and volume.

While more expensive than continuum PDE^{25,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 steps^{16,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}) = *k*_{on} · *ρ*_{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}

## II. THEORY

### A. Single-particle reaction-diffusion: Solution particles binding to explicit lipids on surfaces

#### 1. Binding to explicit lipids on a flat surface

Binding between a particle in solution at *r*_{1} to a lipid particle restricted to a reflective 2D surface at *r*_{2}, 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 particles^{49,50} due to the reflective surface. The diffusion equation for the particle pair can be written as

with their separation ** r** =

*r*_{1}−

*r*_{2}and initial separation

*r*_{0}, and

*p*(

**,**

*r**t*|

*r*_{0}) is referred to as the Green’s function (GF). The diffusion tensor

**(**

*D***) =**

*r*

*D*_{1}+

*D*_{2}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,

*D*

_{2z}= 0). The reaction is modeled via the reactive boundary at

*r*=

*σ*,

where ** J** = −

**(**

*D***)**

*r**p*(

*r*_{12}, t|

*r*_{120}) is the flux at

**and $r^$ is the direction vector of**

*r***. The reactivity**

*r**κ*occurs over the accessible surface of a sphere of radius

*σ*such that

*κ*=

*k*

_{a}/

*B*(

*σ*,

*d*), where

*k*

_{a}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,

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

and the initial condition is a delta function, *p*(*r*, *t* = 0|*r*_{0}) = *δ*(*r* − *r*_{0})/*B*(*σ*, *d*).

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*|

*r*_{0}), can thus be defined from the full 3D problem, with the only complication being the anisotropy of the diffusion matrix

**(**

*D**D*

_{x}=

*D*

_{y}≠

*D*

_{z}). 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\xaf$, 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*|*r*_{0}).^{14} As usual, the reaction probabilities for a pair of particles at a separation of *r*_{0} are defined based on the survival probability such that

producing the same functions used in 3D. To define $D\xaf$, 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*(*θ*, *φ*) = *D*_{x} cos^{2}(*φ*)sin^{2}(*θ*) + *D*_{y} sin^{2}(*φ*)sin^{2}(*θ*) + *D*_{z} cos^{2}(*θ*). Since *D*_{x} = *D*_{y}, then we have *D* = *D*_{x} sin^{2}(*θ*) + *D*_{z} cos^{2}(*θ*). If we weight *D* over all polar angles, then $D\xaf=\u222b0\pi /2D\u2009sin\u2009\theta d\theta =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,

where *k*_{b} 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* = *σ*.

### B. Implicit lipid model: Solution particles binding to implicit lipids on surfaces

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} = *n*_{L}/*S* (*n*_{L}, the number of unbound lipids and *S*, the area of membrane surface). Once a binding event occurs, *n*_{L} 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} = *N*_{L}/*S*, we can define the implicit lipid survival probability for the solution particle in a time Δ*t*, given a height above the membrane of *z*_{0}, 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,

where from the geometry of the problem, we define $r0=r2+z02$ and *r*_{1} = 0 if *z*_{0} ≥ *σ*, and $r1=\sigma 2\u2212z02$ if *z*_{0} < *σ* (see Fig. 1). The well-known reaction probability *p*_{react}(Δ*t*|*r*_{0}) = 1 − *S*(Δ*t*|*r*_{0}) [Eq. (2)] can only be evaluated for *r*_{0} ≥ *σ*^{14} [Fig. 1(c)]. The upper integration limit *r*_{2} 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=\sigma +36D\Delta t$, which gives $r2=Rmax2\u2212z02$, where *D* = $D\xaf$ (Sec. II A 1).

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

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

where $\alpha =D\Delta tka+kd/\sigma kd$, $Z=z0\u2212\sigma /4D\Delta t$, *k*_{d} = 4π*σD*, and *k*_{a} is the intrinsic binding rate. When *z*_{0} < σ, the probability is independent of the specific value of *z*_{0} and is in fact always given by *P*^{IL}(Δ*t*|*z*_{0} = *σ*). A typical curve of *P*^{IL}(Δ*t*|*z*_{0}) 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*.

The dissociation of particles from lipids follows a Poisson process,

where unlike the explicit lipid model, here, we must use *k*_{off}, the macroscopic dissociation rate that satisfies *k*_{off} = *k*_{on}*K*_{D} = *k*_{on} · 2*k*_{b}/*k*_{a}, and *k*_{on}, a simple function of *k*_{a}, *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 *k*_{b}, the intrinsic off rate, and *k*_{off}, the macroscopic off rate. For rate-limited reactions, *k*_{b} and *k*_{off} are the same, but for diffusion-influenced reactions, where *k*_{a} is large, *k*_{b} is faster than *k*_{off} (*k*_{b} ≥ *k*_{off}). 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 *RS*_{3D} 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

where *V* is the volume of the system, *S* is the area of the membrane surface, *N*_{A} is the total number of solution particles, *n*_{A} is the number of unbound solution particles at the equilibrium state, and *RS*_{3D} is the position of the reflecting surface above the membrane. The above function has an analytical solution

where *α* and *k*_{d} are the same as in Eq. (6). If both *k*_{b} and Δ*t* are relatively small, then *RS*_{3D} has the simple expression,

where *k*_{d} = 4*πσD*. We can see that *RS*_{3D} has a value of 0 < *RS*_{3D} < *σ*. The reflecting surface depends only on *σ*, *k*_{a}, 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 *RS*_{3D} above the membrane, via reflection. When dissociation occurs, particles are returned to a distance of *RS*_{3D} 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),

where *S* represents a patch on the surface accessible to the diffusing solution particle within a step ∆t, *r*_{0} 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* > *R*_{max}. 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 *RS*_{3D} 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 *RS*_{3D}. 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 *RS*_{3D} 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 *RS*_{3D}.

#### 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 *z*_{0} above the membrane is always zero such that, using the simpler formulation of Eq. (6), the integral is given by

Here again, the upper limit of integration can be chosen to a point where the reaction probability drops to zero at $Rmax=\sigma +34D\Delta t$. *RS*_{2D} 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, *p*_{react,2D}(Δ*t*|*r*) (see Ref. 53), which is itself an integral over Bessel functions. By taking *p*_{react,2D}(Δ*t*|*r*) into Eq. (12), the binding probability is expressed as

where *γ* = *huY*_{1}(*uσ*) + *k*_{a}*Y*_{0}(*uσ*), *ω* = *huJ*_{1}(*uσ*) + *k*_{a}*J*_{0}(*uσ*), *a* = *uR*_{max}*J*_{1}(*uR*_{max}) − *uR*_{2D}*J*_{1}(*uRS*_{2D}), and *b* = *uR*_{max}*Y*_{1}(*uR*_{max}) − *uR*_{2D}*Y*_{1}(*uRS*_{2D}), *h* = 2*πσD*. *J*_{n} denotes the order *n* Bessel function of the first kind, and *Y*_{n} 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,

where *k*_{off} is the macroscopic dissociation rate that satisfies *k*_{off} = *k*_{on}*K*_{D} = *k*_{on}*k*_{b}/*k*_{a} and a macroscopic association rate *k*_{on} for 2D is determined based on our previous work,^{53}

where $b\rho =2S/\pi \u22c5max(nA,nL)+\sigma 2$, where *S* is the membrane area, *n*_{A} is the number of unbound particles, and *n*_{L} is the number of unbound lipids.

Similar to the 3D case, we need to introduce a reflecting distance, *RS*_{2D} 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

where *γ*, *ω*, *a*, and *b* are defined as before in Eq. (13) (*a* and *b* are functions of *RS*_{2D}). By numerically solving Eq. (17), we obtain the value of *RS*_{2D}. 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,

#### 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 separation^{54} 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.

### C. Macroscopic rates from the microscopic Smoluchowski model

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 *k*_{a} and the macroscopic rate *k*_{on} 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 *k*_{a}. Specifically, the macroscopic rate is known to be related to the survival probability through^{51}

where *r*_{0} are all the points on the reactive surface Γ. For a spherically symmetric interaction, the reactive surface is at *r*_{0} = *σ* 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 *κ* = *k*_{a}/4*πσ*^{2}. As a result,

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

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 *k*_{a} 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

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 $ka3D\u21922D$ = $ka3D$ will follow the weaker *K*_{D} 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 *k*_{a} values are defined based on a 3D reaction, and if one reactant is localized to the membrane, we multiply *k*_{a} by the zero-curvature value of *λ*^{−1} = 2 to preserve the same value of *K*_{D} ($ka3D\u21922D$ = 2$ka3D$), although this results in slightly different *k*_{on} 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).

Systems . | 3D
. | 3D → 2D
. | 2D
. |
---|---|---|---|

Relation to 3D microscopic binding rate | $ka3D$ | $ka3D\u21922D=2ka3D$ | $ka2D=ka3D/2\sigma $ |

Microscopic dissociation rate | k_{b} | k_{b} | k_{b} |

Macroscopic binding rate | $kon=1ka3D+14\pi \sigma D\u22121$ | $kon=121ka3D\u21922D+14\pi \sigma D\u22121$ | $kon=1ka2D+18\pi D4\u2009ln(b\rho /\sigma )(1\u2212\sigma 2/b\rho 2)2\u221221\u2212\sigma 2/b\rho 2\u22121\u22121$ |

Macroscopic dissociation rate | $koff=konkb/ka3D$ | $koff=2konkb/ka3D\u21922D$ | $koff=konkb/ka2D$ |

Equilibrium state | $KD=kb/ka3D$ | $KD=2kb/ka3D\u21922D$ | $KD=kb/ka2D$ |

Systems . | 3D
. | 3D → 2D
. | 2D
. |
---|---|---|---|

Relation to 3D microscopic binding rate | $ka3D$ | $ka3D\u21922D=2ka3D$ | $ka2D=ka3D/2\sigma $ |

Microscopic dissociation rate | k_{b} | k_{b} | k_{b} |

Macroscopic binding rate | $kon=1ka3D+14\pi \sigma D\u22121$ | $kon=121ka3D\u21922D+14\pi \sigma D\u22121$ | $kon=1ka2D+18\pi D4\u2009ln(b\rho /\sigma )(1\u2212\sigma 2/b\rho 2)2\u221221\u2212\sigma 2/b\rho 2\u22121\u22121$ |

Macroscopic dissociation rate | $koff=konkb/ka3D$ | $koff=2konkb/ka3D\u21922D$ | $koff=konkb/ka2D$ |

Equilibrium state | $KD=kb/ka3D$ | $KD=2kb/ka3D\u21922D$ | $KD=kb/ka2D$ |

### D. Background on surface adsorption models

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 *k*_{on} (the explicit or implicit lipid models), we can use *κ* = *ρ*_{L}*k*_{on}, 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} ⩾ 4*k*_{off}*D*), 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,

where $pfree(x,t|x0)=exp[\u2212(x\u2212x0)2/(4Dt)]/4\pi Dt$ and the same for *y*. The solution to the 1D Smoluchowski model *p*_{irr}(*z*, *t*|*z*_{0}), parameterized by the diffusion constant *D* = *D*_{z} and the rate *κ* (with adsorption units of length/time), is known to be^{19,49}

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

This probability depends on *z*_{0}, which can then be compared with the result for *P*^{IL}(Δ*t*|*z*_{0}) 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 *k*_{a}, they are comparable (Fig. S3A), but for a large *k*_{a}, 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 *k*_{a}) 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 *κ*/*k*_{desorb} = [A]_{mem}/[A]_{sol}.

## III. RESULTS

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.

### A. Implicit lipid vs explicit lipid models for planar systems

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 *RS*_{3D} (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 nm^{3}/*μ*s) at increasing lipid densities, illustrating how even with a low concentration of lipids (100/200^{2} nm^{−2}), the kinetics are nearly the same. In Fig. 5(c), we simulate a rate-limited reaction ($ka3D$ = 0.166 nm^{3}/*μ*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.

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.

### B. Adsorption method vs explicit lipid method for planar systems

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 *κ* = *k*_{on} · *ρ*_{L}(*t*) (the macroscopic binding rate *k*_{on} stays the same as in the explicit lipid method) to account for occupancy of lipid binding sites.

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 *k*_{on}), and lipids are of low density [Fig. 8(a)]. As lipid density (and thus *κ*) is increased, the models converge. For rate-limited reactions (small *k*_{on}), 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 *k*_{on} 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.

### C. Binding to lipids in 2D: Implicit vs explicit lipid model

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 (*k*_{a} = 86.80 nm^{2}/*μ*s), rate-limited (*k*_{a} = 0.083 nm^{2}/*μ*s), or diffusion influenced (*k*_{a} = 8.68 nm^{2}/*μ*s), the implicit lipid model produces the proper equilibrium states due to enforcing detailed balance through the *RS*_{2D} 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 *k*_{on}.^{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 *RS*_{2D} ≥ *σ* to ensure equilibrium, we also slow the reaction kinetics. *RS*_{2D} here increases from *σ* at the smallest *k*_{a} value, up to 1.46*σ* for the largest *k*_{a} (see Fig. S1).

### D. Implicit lipid vs explicit lipid model on curved surfaces

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)].

### E. Implicit lipid vs explicit lipid with nonuniform lipid distribution

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.

## IV. DISCUSSION/CONCLUSIONS

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 particles^{18} 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 assembly^{7,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 scale^{64–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.

## V. METHODS

### A. Simulation details of the explicit lipid method

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* < *R*_{max} is determined by the pairwise association probability *p*_{react}(Δ*t*|*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.

### B. Simulation details of the implicit lipid model

For implicit lipid simulations with a flat membrane surface, if a solution particle is within *R*_{max} 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 *RS*_{3D} 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 *p*_{react}(Δ*t*|*r*_{i}) multiplying element area *s*_{i} and the local density of lipids *ρ*_{i}. Then, the total binding probability is accumulated over all these elements based on Eq. (11), $\u2211ri<Rmaxpreact(\Delta t|ri)\u22c5si\rho i$. For curved membranes with analytical shapes such as the sphere, the reaction probability can be calculated analytically for each height *z*_{0} (supplementary material).

For simulations of explicit protein particles restricted to the 2D surface with implicit lipids, simulations are similar to 3D cases, but *D*_{z} = 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.

### C. Simulation details of the surface adsorption model

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 *p*_{irr} solution shown in Eq. (24), for a new *z* value by solving $RAND=\u222b\sigma zpirr(z\u2032,\Delta t|z0)dz\u2032$, where *RAND* is a uniformly distributed random number between 0 and the value of the survival probability 1 − *P*_{AD-GF}(Δ*t*|*z*_{0}). 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.

### D. Triangulated mesh model for complex surfaces

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.

### E. Detecting collisions with the triangulated mesh surface

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 algorithm^{44} 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.

## SUPPLEMENTARY MATERIAL

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.

## ACKNOWLEDGMENTS

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.

## REFERENCES

_{2}(110)