|Home | About | Journals | Submit | Contact Us | Français|
Particle-based simulators represent molecules of interest with point-like particles that diffuse and react in continuous space. These simulators are often used to investigate spatial or stochastic aspects of biochemical systems. This article presents new particle-based simulation algorithms for modeling interactions between molecules and surfaces; they address irreversible and reversible molecular adsorption to, desorption from, and transmission through membranes. Their central elements are: (i) relationships between adsorption, desorption, and transmission coefficients on the one hand, and simulator interaction probabilities on the other, and (ii) probability densities for initial placements of desorbed molecules. These algorithms, which were implemented and tested in the Smoldyn simulator, are accurate, easy to implement, and computationally efficient. They allow longer time steps and better address reversible processes than an algorithm that Erban and Chapman recently presented (Physical Biology 4:16-28, 2007). This article also presents a method for simulating unbounded diffusion in a limited spatial domain using a partially absorbing boundary, as well as new solutions to the diffusion differential equation with reversible Robin boundary conditions.
Cell modeling studies typically represent molecular interactions with mass action methods. These methods are simple, well supported by software (e.g. Copasi  and Virtual Cell ), and fully adequate for a wide variety of studies. However, they cannot represent the natural stochasticity that arises within and between cells , nor the intricate intracellular spatial organization that most cell systems exhibit. To investigate these topics, many scientists use more detailed modeling methods (reviewed in ).
Particle-based simulation is one of these more detailed methods. It represents individual molecules with point-like particles that behave in accordance with elementary processes, including diffusion and chemical reactions. For typical cell biology systems, particle-based simulation can support spatial resolution to scales of several nanometers [1, 14], can simulate up to tens of minutes of real time , and can represent hundreds of thousands of simulated molecules . Freely available particle-based simulators include Cell++ , ChemCell , MCell , and Smoldyn . Using these tools, researchers have shown, for example, that neural signaling may depend upon neurotransmitter release from sites that are away from the synapse  and that stable protein concentration gradients may arise in bacterial cells .
Recently, Erban and Chapman developed particle-based algorithms for simulating molecular adsorption to surfaces  (Figure 1A). They addressed the following problem: given a typical particle-based simulator design [3, 22, 35] (Figure 2), in which molecules diffuse in discrete steps and either adsorb to (with adsorption probability Pa) or reflect from any surface they encounter, how should the simulator compute this adsorption probability from the surface's physical adsorption coefficient? Their result (Eq. 10 of ) is
where κ is the adsorption coefficient (defined below and in the glossary), Δt is the simulation time step, and D is the molecular diffusion coefficient. Bartol and coworkers  derived essentially the same equation 16 years earlier (Eq. 6a of ) for the closely related problem of simulating reactions between solution-phase and surface-bound molecules. Also, Singer and coworkers recently generalized Eq. 1 to multi-dimensional systems with molecular drift and anisotropic diffusion matrices .
All of these derivations of Eq. 1 assumed well-mixed systems, in which molecular concentrations are spatially uniform throughout the solution. An assumption about the average molecular positions is required for the adsorption probability to be independent of the precise position that each molecule diffused from. However, the well-mixed assumption is inconsistent with the problem definition because irreversible molecular adsorption to surfaces necessarily creates a local concentration gradient. Regardless of the starting state, this gradient develops rapidly in the immediate vicinity of the surface and approaches the time-invariant steady state (see section 3.3.1 of  and section 5, below). As a result, molecular concentrations near adsorbing surfaces, whether physical or simulated, are likely to spend the vast majority of their time much closer to the steady state than to the well-mixed state.
In this article, we re-investigate the problem that Erban and Chapman posed, but for steady-state systems rather than well-mixed systems. The relationship we find between the adsorption probability and the adsorption coefficient enables substantially more accurate simulations than Eq. 1 does. We also develop related algorithms for simulating irreversible desorption, reversible adsorption, and either irreversible or reversible transmission through permeable surfaces. Together, these algorithms enable quantitative simulation of most elementary molecule-surface interactions. Examples include peripheral membrane protein binding and release, membrane permeability, and heterogeneous catalysis. Where coarse-grained approximations are appropriate, our algorithms also support simulation of ligand binding to receptors or molecular transport through channels or pores. In addition, our adsorption algorithm enables a method for efficiently modeling unbounded diffusion in a limited spatial domain, which we discuss near the end of this article.
A specific example that illustrates the need for accurate adsorption algorithms arises in the yeast pheromone response system (reviewed in ), a model eukaryotic signaling system. Upon pheromone stimulation, pheromone-sensitive cells execute a series of actions, of which a central one is that the Ste5 protein translocates from the cytoplasm to the plasma membrane . Colman-Lerner and coworkers recently recorded microscopy image sequences that show this translocation . Because their data are two-dimensional images of three-dimensional cells, including in-focus and out-of-focus portions, they require further analysis to determine the adsorption coefficient between Ste5 and the membrane, which could lend insights into the interaction mechanisms. This analysis could be performed relatively easily by adjusting a simple computer model until it agreed with the microscope images; the simulator that ran this model would require accurate adsorption algorithms. Afterwards, the adsorption coefficient between Ste5 and the plasma membrane could be used in new spatial simulations of the pheromone response system. The resulting simulations, which would again require accurate adsorption algorithms, might lend further insights into the yeast signaling system dynamics or signal processing (see e.g. [24, 28]).
We explore one model system for adsorption and desorption and another for partial transmission (Figure 1). The adsorption/desorption system extends over all x > 0 while the partial transmission system extends over all space. Both active surfaces are at x = 0. While these models are three-dimensional for convenient visualization, they can be described equally well with fewer dimensions because their average molecular concentrations are translationally invariant parallel to the surfaces; the mathematics throughout this article uses one dimension.
We assume that solution-phase molecules move solely by diffusion. This means we treat water and other un-modeled species only through their implicit contributions to the Brownian dynamics of the modeled molecules. It also means we ignore intermolecular interactions between modeled molecules, which is typically valid when their concentrations are low. These approximations allow us to describe diffusion with Fick's laws [11, 16]:
where J is the average molecular flux towards positive x and C(x,t) is the average concentration at position x and time t. We call C(x,t), which is technically a distribution function , the concentration profile.
The net adsorption rate in the adsorption/desorption system is
where Ca(t) is the average concentration of molecules that are adsorbed to the surface at time t, and k is the desorption rate constant. The first equality expresses the conservation of molecules between the solution and surface phases. More precisely, it states that the rate of adsorption onto the surface equals the flux of solution-phase molecules towards the surface, measured adjacent to the surface (see Eq. 2). The second equality expresses the adsorption rate in terms of the solution and surface concentrations. The first term on the right side of the equation is the Robin or radiation boundary condition [12, 16, 19], in which the adsorption rate is proportional to the concentration of solution-phase molecules at the surface. The proportionality constant, κ, is the adsorption coefficient. In this model system, κ can range from zero, for a completely inert surface, to infinity, for a surface that adsorbs molecules immediately upon contact. Physically, κ is limited to about the thermal velocity of potential adsorbents , which is (kBT/m)½, where kB is Boltzmann's constant, T is the temperature, and m is the molecule's mass, because this velocity determines how frequently molecules collide with the surface (for a 50 kDa protein at 37°C, κ < 7×106 μm/s). Note that the Robin boundary condition does not depend on the surface-bound concentration, so it does not account for surface saturation (in contrast to the Langmuir equation ). Also, note that authors frequently use the term “adsorption coefficient” to mean other things, such as the adsorption equilibrium constant [9, 10, 40] (see below) or an exponent for a variant of the Langmuir equation . The second term on the right side of Eq. 4 is the desorption rate. Its proportionality constant, k, can also range from zero to infinity, now representing no desorption and instantaneous desorption, respectively. According to Chou and D'Orsogna , desorption rates are typically between 10−4 s−1 and 104 s−1.
The net flux of molecules through the permeable surface, from front to back, is
In analogy to prior definitions, κF and κB are the permeability coefficients from the front of the surface (which faces the positive x-axis) to the back, and from the back to the front, respectively. Eq. 5 specifies the surface sides with +0 for the front and –0 for the back because the concentration profile is typically discontinuous at x = 0. As before, the first equality expresses the conservation of molecules and the second equality represents the Robin boundary condition. Again, this equation ignores saturation; for example, it does not account for ion channels operating at capacity. In general, permeability coefficients are large for molecules that tend to partition into the membrane as well as for molecules with high diffusion coefficients within the membrane (see ). Paula et al.  found that permeability coefficients across a 27 Å phosphatidylcholine bilayer are about 3.5×10−8 μm/s for potassium ions, 0.014 μm/s for urea, 0.027 μm/s for glycerol, and 150 μm/s for water.
In Step 1, the simulator diffuses molecules for a fixed amount of time with Gaussian-distributed random displacements, while ignoring any surfaces or bimolecular interactions. The standard deviation of these Gaussian distributions, called the root mean square (rms) step length, is
This is a characteristic length for the simulation, against which we previously compared intermolecular interaction distances  and surface curvature radii . Molecules move about an rms step length at each time step, so the planar surfaces in our model systems accurately represent membranes that have curvature radii which are much larger than the rms step length, and are less accurate for more tightly curved membranes. For either model system, a diffusive step transforms a starting concentration profile, C(x,t), to
The star represents convolution and the term on its right is a Gaussian with unit area and standard deviation s.
In Step 2, the simulator addresses any solution-phase molecules that diffused across the surface in Step 1, as well as the possible desorption of any surface-bound molecules. The former molecules adsorb to the surface with probability Pa, transmit through it with probability Pt, and otherwise reflect off of it using ballistic reflection . The latter molecules desorb with probability Pd; the simulator displaces desorbed molecules away from the surface according to the probability density p(x). We solve for Pa, Pt, Pd, and p(x) below. The simulator does not perform multiple surface interactions, such as both adsorption and desorption, on the same molecule so as to improve computational efficiency and to make their order of execution irrelevant (however, the algorithms do account for the multiple surface interactions that physical molecules may undergo). Step 2 transforms the adsorption/desorption system concentration profile to
It also transforms the adsorbed concentration to
The analogs of these equations for the partial transmission system are not shown because they are more complicated and less useful.
This design for Step 2 incorporates two simplifications for good computational efficiency. First, Pa and Pt do not depend on how far molecules are from surfaces, but instead depend only on static properties, such as diffusion and adsorption coefficients. This allows a simulator to compute these probabilities just once for a simulation, rather than once for each molecule-surface interaction. Secondly, the simulator ignores all molecules that diffused across the surface during Step 1, but that returned back to the side where they started before the end of the time step (see section 3.2 of  and section 2.2 of ). In preliminary work, accounting for these additional surface contacts proved computationally intensive, complicated to implement, and minimally more accurate.
In Step 3, the simulator performs chemical reactions . We assume that any reactions are slow enough that their effects on concentrations near surfaces are negligible, where the rms step length provides the necessary length scale. In Step 4, the simulator increments the simulation time. This updates the concentration profile and surface concentration to
Eq. 10 applies to both model systems, while Eq. 11 only applies to the adsorption/desorption system. By definition, C(x,t+Δt) equals the starting concentration profile, C(x,t), if and only if they represent a steady state. Ordinarily, the adsorbed concentration, Ca(t), would also have to be constant at steady-state. In fact, we require this time independence for reversible adsorption but not for irreversible adsorption. We do not require it for irreversible adsorption because (i) we wish to investigate a constant adsorption flux, and (ii) the adsorbed concentration does not affect other system aspects. Finally, the simulator performs any observations of the system in Step 5. It then returns to Step 1 to start the next time step.
In this section, we relate the molecule-surface interaction probabilities to their respective coefficients for accurate steady-state simulation. We also derive the necessary probability densities for desorbed molecule displacements. These results are presented as five algorithms, along with several variants. Figure 5 presents the main results graphically, the Supplementary Information presents them in tabular form, and several equations below approximate them with interpolating functions. Also, the Supplementary Information includes the C language library files that Smoldyn uses to perform the necessary table look-up and interpolation. Except as noted, most results in this section were computed with Mathematica software .
Unitless reduced variables, shown with a prime symbol, simplify several results. We use the simulation time step as the standard time unit and the rms step length as the standard length unit. From these, the reduced adsorption coefficient, desorption rate constant, and diffusion coefficient are:
Reduced permeability coefficients, κF′ and κB′, are analogous to the reduced adsorption coefficient (Eq. 12).
Consider the adsorption/desorption system and assume irreversible adsorption (κ > 0 and k = 0). Solving Eqs. 3 and 4 shows that the model steady-state concentration profile, C(x), increases linearly with distance away from the surface and the surface concentration increases at a constant rate (Figure 3A),
where C0 is the solution-phase concentration at x = 0. (Eq. 3 and the outer portions of Eq. 4 are the equations of motion for C(x,t) and Ca(t), respectively, while the latter two portions of Eq. 4 give the necessary boundary condition.) Eq. 15 becomes physically unreasonable as the system increases in size because C(x) is unbounded; however, this model is likely to be physically correct near the surface, so it is adequate for our work. As before, note that the concentration profile and adsorption flux are time-invariant at steady state.
In an exact algorithm, meaning one in which the simulated system is theoretically indistinguishable from the underlying model system, the average simulated steady-state concentration profile would equal Eq. 15. We determine whether this is possible by seeing if our simulator design can preserve the model steady-state profile over one time step. Step 1 (Eq. 7) transforms Eq. 15 to
for x > 0. Eq. 18 is the concentration profile at the end of the time step. Comparing Eqs. 15 and 18 shows that they cannot equal each other for arbitrary κ′ and x values. Thus, Eq. 15 is a steady state for the model system but not for our simulator design, which means that our simulator design cannot treat irreversible adsorption exactly. (The simplifications in Step 2, which we included for computational efficiency, cause this inexactness.)
Knowing that simulated concentration profiles would necessarily be inexact, I computed the adsorption probabilities that yield the correct adsorption rates with a “simulation emulator” (see Figure 4A and the Supplementary Information). Each emulator run started with a fixed adsorption probability and a well-mixed numerical concentration profile (which extended from −6 to +10 rms step lengths and was tabulated at 401 equally spaced points, plus additional points immediately above and below 0 for higher accuracy). The emulator iterated Steps 1 and 2, using Eqs. 7-9, until the system was essentially at steady state (as determined by the adsorption rate changing by less than 0.01% between sequential iterations). It addressed the region with x < −6s by extrapolating the concentration profile with an error function, and it addressed the region with x > 10s by fixing these concentrations to 1. Once at steady state, it fit a straight line to the final concentration profile (x ranged from 3s to 7s to avoid edge effects), from which it calculated the effective reduced adsorption coefficient using
Eq. 19 derives from Eq. 15, with m as the concentration profile slope (C0κ/D) and C0 as the intercept. Because the emulator could not actually reach the steady state, it performed the above routine twice, starting from different states; it started with C(x,0) equal to either 0 or 1 so that it would approach the steady-state profile from both above and below. Finally, the emulator averaged the two results to yield a final κ′ value. These κ′ values are almost certainly accurate to within ±10–6 because (i) the two intermediate solutions never differed by more than 10–7, and (ii) more conservative choices for emulator parameters (number of tabulated points, size of x-domain, and number of iterations) affected results by even smaller amounts.
Figure 5A shows the emulator results, oriented so that Pa depends on κ′. Figure 5A also shows stochastic simulation results, generated by Smoldyn, for systems that were similar to those that the emulator quantified (Supplementary Information). The close agreement between the two sets of results lends confidence in the emulator software, the Smoldyn software, and the mathematical analysis. When reduced adsorption coefficients are small, the emulator results agree with Eq. 1. This makes sense because small κ′ values imply either: (i) small concentration gradients (because κ is small or D is large), so the steady state is similar to the well-mixed state, or (ii) a small simulation time step, so molecules diffuse over short distances and the concentration gradient affects the simulated adsorption rate minimally. On the other hand, the emulator results are up to a factor of two smaller than Eq. 1 for larger κ′ values. This implies that simulators basing their adsorption probabilities on Eq. 1 would simulate adsorption too quickly in many situations. Finally, Figure 5A shows that κ′ values larger than about 0.86 cannot be modeled with our simulator design because the adsorption probability cannot exceed 1. This means that large adsorption coefficients can only be simulated quantitatively by using small time steps. (This limitation arises from simplifications that we made for computational efficiency in the simulator's Step 2.)
Polynomial fits to the emulator's relationship between κ′ and Pa, which may be more convenient than the raw data, are:
The initial terms were constrained with Eq. 1 to make these fits exact as κ′ or Pa tend toward 0. These equations agree with the tabulated data to within 1% over their entire domains.
Continuing with the adsorption/desorption system, we now consider irreversible desorption (κ = 0 and k > 0). Desorption is a simple first order process, like a first order chemical reaction, so a molecule should desorb during a time step with probability [3, 42]
This is exact for all possible desorption rate constants and simulation time steps.
The more difficult question concerns where the simulator should place desorbed molecules. In one algorithm variant, the simulator places them adjacent to the surface. This ignores diffusion occurring between when a molecule desorbs and the end of the time step, which does not affect the accuracy of the total surface- and solution-phase concentrations, but means that this algorithm variant yields inexact concentration profiles.
A better variant displaces desorbed molecules away from the surface. To calculate the probability density for this displacement, consider a molecule that desorbs at time t, where t is between the start of the time step at time 0 and the end at time Δt. The probability density for this molecule's position at time Δt is 
We denote the probability density that the molecule desorbs at time t, given that it desorbs between 0 and Δt, as pt(t). Treating pt(t) as a weighting factor for Eq. 23, the spatial probability density for a molecule that desorbs at an unknown time between 0 and Δt is
To be consistent with Eq. 22, pt(t) should equal k/Pd exp(–kt). However, this makes p(x,Δt) depend on k, which is inconvenient. So, we make the approximation that the molecule is equally likely to desorb at any time during the time step, which simplifies pt(t) to 1/Δt. This approximation is valid when kΔt << 1 because the exact exponential function approaches 1/Δt in this limit. Because Eq. 22 already ensures accurate surface and solution concentrations and we are just refining the solution concentration profile here, this approximation is also likely to be acceptable for substantially larger values of kΔt as well. Using the approximation, Eq. 24 integrates to
Figure 5B shows this displacement probability density with a solid line.
A final task makes Eq. 25 useful. We need a relationship with which a simulator can transform random numbers that are uniformly distributed between 0 and 1 into random numbers that obey the probability density in Eq. 25. We do this by integrating and inverting Eq. 25 . The integral, from 0 to x, is
This cannot be solved for x analytically, so I inverted it numerically. A least-squares rational function fit to the result is
Eq. 27 converts a uniform deviate, P, into the desired deviate, x. Compared to the numerical solution, it is in error by less than 0.02s for nearly all P values.
If adsorption is reversible (κ > 0 and k > 0), the adsorption and desorption rates equal each other for the system to be at steady state. Adding this condition to Eqs. 3 and 4 yields the steady-state concentrations (Figure 3A),
where C∞ is the solution concentration. In contrast to the preceding investigations of irreversible processes, this steady state is a well-mixed state. This steady state is also an equilibrium state, because there are no net molecular fluxes. The ratio κ/k, which clearly equals Ca/C∞, is the adsorption equilibrium constant; it is analogous to the familiar equilibrium constant for reversible reactions .
As in section 4.1, we determine whether our simulator design can treat reversible adsorption exactly by seeing if it can preserve the model concentration profile over one time step (see Figure 4B). Transforming Eq. 28 with Eq. 7 yields
The adsorption and reflection portions of Step 2 (see Eq. 8) then transform the concentration profile to
Finally, the desorption portion of Step 2 transforms this result back to the starting condition, if and only if:
Thus, our simulation design can simulate reversible adsorption exactly under equilibrium conditions. Eq. 32, which relates the adsorption and desorption probabilities, arises from the constraint that equal numbers of molecules must adsorb and desorb during each time step, on average. Eq. 33, which is the probability density for desorbed molecule displacements and is shown in Figure 5B with a dashed line, causes desorbed molecules to exactly replace those that adsorbed. This density is the ratio of the function C(x,t + Δt) – C(2a)(x,t) and its integral (the “desorbed molecules” region of Figure 4B), so that the probability density has unit area.
We integrate and invert Eq. 33 to make it useful. Integration yields
A rational function fit to the numerically calculated inverse function is
As before, P is a uniform deviate between 0 and 1 and x is the initial displacement for a desorbed molecule. This fit is in error by less than 0.02s for nearly all P values.
The remaining question is how Pa or Pd should be chosen. Their solutions for irreversible adsorption (Eq. 21) and irreversible desorption (Eq. 22) do not simultaneously satisfy the equilibrium condition of Eq. 32 because the prior results did not consider molecules that both adsorb and desorb during a single time step. We address them with the following conceptual scheme: we label all molecules as being solution-phase or surface-bound at time 0, the system evolves over one time step, and we then use the labels to see how many molecules had a net adsorption transition by time Δt. We solve this problem in Appendix A. The answer, Eq. A12, represents the model adsorbed molecules at the end of a time step. Meanwhile, the expectation simulated adsorbed molecules for one time step (see Figure 4B and Eq. 9) is
Equating these amounts and solving for Pa, yields the adsorption probability,
Figure 5C shows this adsorption probability and the corresponding desorption probability (from Eq. 32) for different κ′ and k′ values. As before, (i) the adsorption probability approaches Eq. 1 as κ′ is reduced to zero, (ii) Eq. 1 simulates adsorption too quickly with larger κ′ values, and (iii) κ′ values that exceed some cut-off, which now depends on k′, cannot be simulated quantitatively with our simulator design because adsorption probabilities cannot exceed 1. In the limit that κ′ is reduced to zero, Pd approaches its value for irreversible desorption (Eq. 22), as one would expect. In contrast though, Pa does not approach the irreversible adsorption probability (Eq. 21) as k′ is reduced to zero. This is because the concentration profile remains well-mixed for the reversible adsorption case, even as k′ is reduced towards zero, which is a consequence of the fact that the reversible adsorption steady state is an equilibrium state (as k′ is reduced towards zero, the adsorption equilibrium constant and the adsorbed concentration increase towards infinity); in contrast, the concentration profile is sloped for the irreversible adsorption case.
The partially transmitting system (Figure 1B), with irreversible transmission from the front face to the back face (κF > 0 and κB = 0), is essentially identical to the irreversible adsorption system. Its steady-state concentration profile is
where C+0 represents the concentration at the front of the surface and is used if x > 0, and C–0 represents the concentration at the back of the surface and is used if x < 0. As we found before: our simulation design cannot simulate concentration profiles exactly for steady-state irreversible transmission, but it can achieve accurate concentrations on each side of the surface; Figure 5A shows the transmission probability as a function of the transmission coefficient; and Eqs. 20 and 21 are polynomial fits to the relationship between Pt and κF.
The question remains where the simulator should place transmitted molecules. It can, of course, place them at the back side of the surface, although this ignores diffusion occurring between transmission and the end of the time step. A better algorithm variant accounts for this diffusion using the probability density given in Eq. 25. Whereas this density was inexact before, it is exact now because a steady-state irreversibly transmitting surface acts as a constant molecular source to the region with x < 0. In yet a third algorithm variant, the simulator does not move transmitted molecules, but simply leaves them where they ended up after Step 1 of the simulation process. The resulting probability density (which can be approximated with the x < 0 portion of Eq. 17) is very similar to the exact result. This is the most computationally efficient variant.
Finally, consider the reversible partially transmitting system (κF > 0 and κB > 0). From Eqs. 3 and 5, the model steady state is well-mixed on each side of the surface but is likely discontinuous at the surface,
This steady state is an equilibrium state, for which each side of Eq. 42 represents the chemical partition coefficient . Physically, the concentration discontinuity can be maintained actively, such as by trans-membrane molecular pumps, or passively, such as by chemical partitioning between different solvents.
Our simulator design can simulate reversible transmission exactly under equilibrium conditions. The starting concentration profile (Eq. 41) diffuses to
where the first term represents molecules that diffuse from x < 0 and the second term represents molecules that diffuse from x > 0. In Step 2, the simulator transmits molecules that diffused across the surface with probabilities PtF and PtB for the front and back faces, respectively. The number of molecules that the simulator transmits from the front to the back is the integral of the second term of Eq. 43 from –∞ to 0 (see Eq. 36), times PtF. This result, and the analogous one for transmission the other way, is
These are equal at steady state. Combining them with Eq. 42 yields
To return the concentration profile to the steady state by the end of the time step, the simulator can displace transmitted molecules away from the surface according to Eq. 25 or just leave them where they ended up after Step 1. These produce the same steady-state concentration profiles and so are equally accurate at steady state, although the latter variant is more computationally efficient.
We find PtF much as we found Pa for the reversible adsorption problem. Now, we solve: given an initially well-mixed concentration profile for x > 0, no concentration for x < 0, and a reversible partially transmitting boundary, what is the total number of molecules in the region x < 0 at time Δt? Appendix B presents the result. The solution, Eq. B12, represents the number of molecules that the simulator should transmit over one time step. We equate it to Eq. 44, which represents the number that the simulator does transmit, and solve for the transmission probabilities. They are:
Figure 5D presents these probabilities.
The results are familiar: transmission probabilities approach Eq. 1 as κF′ and κB′ approach zero; Eq. 1 simulates transmission too quickly with larger κF′ and κB′ values; and, in most cases, a simulator can only quantitatively simulate transmission up to finite κF′ and κB′ values. The exception to the last result occurs if κF′ = κB′ (the bold dashed line in Figure 5B), which is typical for passive membranes.
All of the algorithms presented above rest on firm theoretical foundations that ensure high accuracy for steady-state systems. Here, we investigate their accuracy away from steady state. Figure 6 shows results from several algorithm tests, run with the Smoldyn simulator. Each test started with 20,000 molecules that were well-mixed within box-shaped regions and that had 5 μm2/s diffusion coefficients (see Figure 1; x extended from 0 to 2 μm and y and z were each 1 μm wide). These are reasonably typical parameters for proteins in bacteria [17, 43]. Initially, no molecules were adsorbed to a surface at x = 0 and there were no molecules in the region x < 0. Simulation time steps were 1 ms, which is larger than the 0.25 μs to 0.1 ms time steps that modelers typically use [14, 25, 29-31], both because this tested the algorithms where they are likely to be weakest and because a goal of this work is to enable accurate simulations with longer time steps. Depending on the specific test, molecules interacted with the surface with irreversible adsorption, reversible adsorption, or reversible transmission. These tests used either “fast” rates, with the maximum possible adsorption or transmission coefficients that the time step allowed, or “slow” rates, with coefficients that led to adsorption or transmission probabilities to be about 0.1 (the Figure 6 caption presents the coefficients). Figure 6A compares these stochastic simulation test results with exact calculations, where the latter are:
(Eq. 3.37 of ) for irreversible adsorption, Eq. A12 for reversible adsorption, and Eq. B12 for reversible transmission. In all cases, the simulation results agreed essentially perfectly with the model ones. This is true at all time points, despite the fact that the system was far from steady state for the initial ones. Additional tests (not shown) performed equally well for smaller adsorption and transmission coefficients.
Figure 6B shows several snapshots of molecular spatial distributions from the irreversible adsorption test that used the “fast” rate. Comparison with the model profiles,
(Eq. 3.35 of ), again showed excellent agreement at all time points. As before, additional tests that used smaller adsorption coefficients performed equally well. Note that even though the system shown in Figure 6B was well-mixed initially, its concentration profile evolved rapidly in the immediate vicinity of the surface (the rms step length for this simulation was 0.1 μm) toward a steady state profile; this is the reason why we developed our algorithms for steady-state systems.
These tests indicate that our irreversible adsorption, reversible adsorption, and reversible transmission algorithms are accurate not only for steady-state systems, but also for those away from steady state. Furthermore, they are accurate over the entire range of adsorption or transmission coefficients accessible for a given simulation time step. Considering the other two algorithms, the irreversible desorption probability presented above (Eq. 22) is accurate for all systems because it does not rely on the steady-state assumption, and the irreversible transmission algorithm is essentially identical to the irreversible adsorption algorithm, which we just verified.
In contrast, simulations based on Eq. 1 performed less well. First, they could not be tested with the “fast” set of adsorption or transmission coefficients because Eq. 1 does support such large coefficients (Eq. 1 only allows κ′ values up to 0.399, whereas they can reach 0.859 for our irreversible adsorption algorithm; see Figure 5A). For the “slow” set of coefficients, adsorption was about 6% too fast and partial transmission about 10% too fast, at all time points. With larger adsorption or transmission coefficients, use of Eq. 1 led simulations to irreversibly adsorb molecules up to 2.1 times too fast, or to reversibly adsorb or transmit molecules with even larger errors (the maximum errors depended on the rates of the reverse processes).
Our final algorithm, which applies to all spatial simulation methods, builds on the irreversible adsorption algorithm. It addresses the following situation: suppose a modeler wants to investigate biochemical interactions in unbounded space, such as between extracellular pheromones and cell-surface receptors, but wants to limit the simulated spatial domain for computational efficiency . The modeler needs to know how much space to simulate explicitly and how to treat its boundary. Clearly, the boundary should be neither completely reflective, which would keep molecules from diffusing away permanently, nor completely absorbing, which would inaccurately represent those molecules that diffuse away and then return. We propose the intermediate solution that the boundary should absorb molecules according to whatever absorption coefficient causes the simulated concentrations to equal those for the unbounded model system. We solve for this absorption coefficient here.
We derive the absorption coefficient for three-dimensional systems with stationary sources that emit molecules at constant rates, without significant obstacles to diffusion, and without significant molecular turnover through reactions. We number the sources with index j: source j is located at rj and emits molecules with rate qj. For our unbounded model, the steady-state molecular concentration at location r is 
The flux at r is computed from the concentration gradient using Fick's law, as in Eq. 2,
The Robin boundary (Eq. 4) relates the flux into the surface with the concentration at the surface, which we generalize to more than one dimension to yield
where ñ(r) is the unit outward normal for the surface at r; also, note that κ, the absorption coefficient, is now made a function of r. We equate the fluxes in Eqs. 52 and 53 and solve for the adsorption coefficient,
Thus, if a surface absorbs molecules with this absorption coefficient, the concentration profile will be the same as if the system were unbounded.
Figure 7 illustrates this algorithm. Using the Smoldyn program, I defined two molecular sources and bounded the system with a cube that was 5 μm on each side. Each cube face was tiled with 25 panels, each 1 μm square. Each of these panels absorbed molecules according to the coefficient given in Eq. 54, calculated with r set to the panel center. The upper portion of Figure 7 shows a slice through the middle of the system and the lower portion compares average steady-state molecular concentrations in a band through the system with model concentrations from Eq. 51. The simulated and model systems agree well.
Erban and Chapman's adsorption algorithm for particle-based simulators , based on Eq. 1, is simple to understand, easy to implement, and fast to simulate. These qualities make it ideal for “quick-and-dirty” simulations. However, it is only quantitatively accurate if (i) the modeled system is kept well-mixed, even in the immediate vicinity of adsorbing surfaces, or (ii) simulation time steps or adsorption coefficients are very small (more precisely, κ′ << (2π)−1/2). We built on their result here with algorithms that yield more accurate adsorbed concentrations and solution-phase concentration distributions, and also with new algorithms for desorption and partial transmission. We derived these algorithms for steady-state systems, but also showed that they are very accurate for systems that start well-mixed. They require essentially the same computational effort as Erban and Chapman's method at each simulation time step but may enable faster simulations overall because they maintain their accuracy as time steps are made longer.
This work shares two important themes with our prior work on simulating bimolecular reactions . First, we designed both sets of algorithms to yield correct rates when the systems are at steady state rather than when they are well-mixed. We did this because the adsorption and reaction processes we wished to model cause systems to evolve rapidly away from well-mixed states (or other initial states) and towards steady states. As a result, portions of these systems close to adsorbing surfaces or reactive molecules, which are the regions that determine adsorption and reaction rates, spend nearly all of their time close to steady state. As we saw above and in prior work , steady-state investigations do not lead to the simple analytical equations that other researchers found using well-mixed assumptions [18, 22, 41], but are nevertheless necessary for accurate simulations. Secondly, we found that reversible interactions need to be considered as special cases, rather than as independent forward and reverse processes. We saw that here for reversible adsorption and reversible transmission, and previously for reversible bimolecular chemical reactions [1, 3]. Treating the processes as independent forward and reverse processes leads to incorrect transition rates and, more importantly, incorrect equilibrium constants.
The algorithms presented here will, hopefully, make it easier for researchers to quantify biological adsorption, desorption, and transmission rates by using computational modeling. This will help address the current paucity of quantitative data for these interactions. Combining these algorithms with the better data should also enable particle-based biological simulations that are faster and more accurate than current ones. These investigations may help scientists better understand dynamic protein localization, intracellular molecular gradients, and related topics. To promote these research directions, all algorithms presented here are implemented in a code library that is supplied with the Supplementary Information, and in the general-purpose Smoldyn biochemical simulator (available at http://www.smoldyn.org).
This publication was made possible by Grant Number P50 HG002370 from the NHGRI, and also by MITRE Corp. contract number 77459; R. Brent was the principle investigator for both grants. Its contents are solely the responsibility of the author and do not necessarily represent the official views of the NIH. I thank Attila Szabo for help with the Appendix results; Radek Erban and Jon Chapman for helpful discussions; Nathan Addy, Esther Mecking, and Orna Resnekov for editing assistance; and anonymous reviewers for exceptionally helpful reviews.
Units are given in parentheses, assuming a 3-dimensional system; L = length, T = time, and 1 = unitless.
This section solves the following problem. In the adsorption/desorption model system, defined by Eqs. 3 and 4, suppose the solution is well-mixed with concentration C∞ and there are no adsorbed molecules at time 0. What is the adsorbed concentration at time t?
A “hat” symbol denotes a Laplace transformed function and z is the Laplace conjugate of the time (we use z instead of s, which is conventional, because we use s for the rms step length). In addition, three boundary conditions are
Combination with Eq. A6 simplifies the result to
We solve this for b to find
Using the method of partial fractions, Eq. A11 inverse Laplace transforms to yield the surface concentration as a function of time
where c1 and c2 are variables that the partial fraction expansion introduced, but which do not have a clear physical interpretation. This solution applies for all physically reasonable (i.e. non-negative) values of the parameters, κ, k, D, and t. However, it may not be particularly convenient for many parameter choices because the constants c1 and c2 are complex for many reasonable values of κ and k. Note that the solution-phase concentration profile can be found from results presented in section 12.4.III.iii of .
This section is very similar to Appendix A. In the partial transmission model system, defined by Eqs. 3 and 5, suppose the concentration is 0 for all x < 0 and is C∞ for all x > 0 at time 0. How much substance is in the x < 0 portion of the system at time t?
CB(x,t) and CF(x,t), along with their Laplace transformed versions, are the concentration profiles on the back and front sides of the surface, respectively. Some additional boundary conditions are
The first two of these boundary conditions simplify the differential equations, Eqs. B1 and B2. They are then solved (see Eq. A6) and simplified further with the latter two boundary conditions to yield
In these equations, a and b are constants of integration. These general solutions are combined with the Robin boundary condition in Eq. B3 and simplified to yield solutions for a and b,
Substitution of Eq. B10 into Eqs. B8 and B9 yields the Laplace transforms of the concentration profiles. The inverse Laplace transform, which would yield the time-dependent concentration profiles appears to be intractable. Instead, we work towards the desired solution by integrating the concentration profile on the back side of the surface to find the total number of transmitted molecules. Its Laplace transform is
This inverse Laplace transforms to yield the desired solution
This solution applies to all physically reasonable values of the parameters (κF > 0, κB > 0, t > 0, and D > 0).
The Supplementary Information for this article includes (i) details about the simulation emulator, (ii) tables that relate adsorption coefficients to adsorption probabilities for irreversible and reversible adsorption, (iii) details about stochastic simulations used in Figures 5A, ,66 and and7,7, and (iv) C source code that converts between molecule-surface interaction coefficients and their simulation probabilities, and also includes the simulation emulator, and (v) documentation for the source code.