Description of the agent-based simulation
Only the molecules of interest are tracked in the simulation in order to reduce the complexity, which allows modeling the whole cell [49]. Since no solvent molecules are present, the stochastic force pushing the tracer molecules around (leading to diffusion) is implicitly included in a random walk model. The position of the tracer molecules is updated in every step Δt according to [50]
where D0 is the diffusion coefficient. ζ is a Gaussian random number with mean 0 and variance 1. If the motion step would end in an obstacle, it is rejected and the molecule stays at the previous position waiting for the next timestep bringing a new chance to move [16].
The simulation only requires defining the particle radius and diffusion coefficient for each species as well as the initial number and distribution of the molecules. The particles are initially placed in the virtual cell at positions which are not restricted by the cytoskeleton or crowding molecules. Likewise all reactions have to be defined (educts, products, rate constants).
Reactions between molecules
A reaction between two molecules can only occur, if the reactants are close enough together. The reaction probability between two molecules is therefore given by the probability of the collision and the probability that a reaction occurs given that a collision is occurring. The claim of the simulation is that it can reproduce the macroscopic (mass action kinetics) rate constant k
ij
for a bimolecular reaction between species i and j in homogeneous conditions. This means that the combination of the collision and reaction probability has to yield the given macroscopic reaction rate.
The discrete time simulation framework complicates the estimation of the reaction probability. The position of the molecules is only known at t
n
and tn+1= t
n
+ Δt. All the collisions that happen within the interval [t, t + Δt) are not directly accessible. Furthermore, the number of tractable collisions depends on Δt. A smaller Δt leads to a finer sampling of the original time interval and reveals more collisions. The reaction probability per collision therefore has to be adjusted with Δt.
The gap between t and t + Δt could be closed if the probability density distribution of the relative motion of a pair of molecules is tracked using the Fokker-Planck equation [51, 52]
with the initial probability density function is
(
is the initial separation) and the partially absorbing boundary condition at the collision distance.
The flux across the boundary within [t, t + Δt) is determined by the surface reaction rate constant k'. The flux accordingly returns the diffusion controlled reaction probability for two molecules which are initially separated by the distance r0. By this approach the requested true number of reactions can be estimated [51, 52]. In order to increase the performance of the simulation, a simpler method was developed based on the approach of Pogson et al. [53].
Derivation of the simulation method
As shown by Pogson et al. [53], the calculation of the bimolecular reaction properties in an event-based stochastic framework can be derived from the macroscopic (ordinary differential equation) description
For simplicity, the units of the concentrations should be [molecules/μ m3] (note: the units of k
ij
are then [μ m3/molecules/s]. If c
i
is given in [mol/liter], the factor N
A
/1015 [liter/(μ m3mol)], where N
A
is the Avogadro number, is required for the respective conversion). The balance equation can be discretized
and converted according to c
i
× V
cell
= N
i
in order to track molecule numbers:
ΔN describes the change in the number of molecules and thus the number of reactions within Δt. Within the timestep Δt the fraction
of the i and j molecules respectively will react. According to mass action kinetics, the number of reactions is proportional to the number of reaction partners and the rate constant k
ij
.
The corresponding 'reaction' volume
can be introduced [53], which indeed has the units of volume: the units of k
ij
are in the given framework [volume/molecules × 1/s], and the unit of Δt is [s]. So (k
ij
Δt) gives just [volume/molecules] which can be further specified to [reaction volume/molecule]. Replacing Δt × Δk
ij
in Equation (17) by ΔV leads to
So in the completely homogeneous framework the fraction of reacting molecules corresponds to a fraction of the volume in which all molecules react, while they do not react in the remaining volume.
From the perspective of the i molecules, the reaction volume is located at the j molecules with which they react. Accordingly, it can be wrapped around these molecules. For symmetry reasons the reaction volume should be spherical. The corresponding 'reaction' radius of the spherical reaction volume is then
This reaction radius is used by Pogson et. al. [53]. Since the reaction volume is wrapped around the j molecules from the perspective of the i molecules, it is multiplied with N
j
as required by Equation (19) leading to the right number of reactions in the homogeneous case. From the perspective of j molecules it is accordingly multiplied by N
i
. Since the reaction volume is now bound to the molecules and not arbitrarily distributed in the cell, the approach is not limited to uniform particle distributions.
A molecule of species i will react with one molecule of species j, if both are closer than the critical reaction distance
(and vice versa j will react with i). This distance grows ∝ Δt1/3 according to Equation (20), compensating the Δt-dependency in the sampled collision frequency in the given framework, which was discussed above. Note, that the particles have to overlap in order to react. This is allowed in the present simulation framework for the moving agents. Since the concentration of each signaling molecule species is remarkably low, so the differences between overlapping and non-overlapping molecules are negligible. In contrast, in the case of high concentrations, a particle based simulation tracking each molecule individually becomes unreasonable [46].
From macroscopic theory to one microscopic reaction: impact of diffusion
Initially, a uniform random distribution of molecules is assumed. On average, ΔN
i
molecules are closer to j molecules than the reaction distance and will therefore react. Subsequently, the reaction rate depends on the flux of the remaining i molecules towards the remaining j molecules. The number of molecules entering the reaction volume is accordingly determined by the combined diffusion coefficient D
i
+ D
j
and the size of the surface of the reaction volume. A comparison with the theory of diffusion limited reactions leads to the following conclusions [40]:
-
The collision rate constant between the reactants is given by k
D
= 4π(r
i
+ r
j
)(D
i
+ D
j
) [41].
-
If not all collisions lead to a reaction but only a fraction which is determined by a microscopic reaction rate k
micro
, the macroscopic reaction rate is given by [40]:
(cf. the results section on effective reaction rates). Obviously the microscopic reaction rate constant k
micro
has to be used in the present collision based simulation. Since only macroscopic reaction rate constants are given in the literature, k
ij
is used as input parameter. The corresponding micro rate constant is then
-
The collision radius (r
i
+ r
j
) will most likely not match the critical reaction radius determined above in Equation (20). Using an artificially smaller radius would reduce the collision rate constant k
D
- requiring a different k
micro
. In other words, the flux through the reaction surface would be reduced due to the smaller surface area. In turn, a higher fraction of this flux has to react in order to yield the original macroscopic reaction rate. A different reaction rate constant would however lead to a different critical reaction radius in the volume-based framework developed above.
-
Solution: reaction probability in the interaction volume: Both concepts, the flux/surface-based description and the macroscopic, volume-based framework can be brought into agreement in the following way:
-
1.
The true collision radius
is used in the simulation as critical reaction radius.
-
2.
The microscopic reaction rate constant k
micro
is calculated as shown in Equation (22) and subsequently used to determine the fraction of the collisions which lead to a reaction within Δt:
but this will (most likely) not match the collision volume 4π = 3(r
i
+ r
j
)3
which effectively reduces the reaction volume determined by the collision radius to the reaction volume given by Equation (23) while it retains the correct interaction surface.
This approach also reflects the nature of reactions in a probabilistic framework: the overall, macroscopic reaction probability is now determined by the probability to collide and the probability to react, given that a collision has occurred.
Adsorption to cellular structures
The association with the cytoskeleton can be described in the same way, and also the adsorption to surfaces like the plasma membrane. Since these objects are impenetrable, the reaction volume has to be outside of the cellular structures - leading to a reaction layer around them. The height of this layer is given by k
binding
× Δt, and the reaction probability is 1. The total reaction volume is determined by multiplying the height with the total surface of the structures (the curvature can be neglected because k
binding
× Δt ≪ r
structure
). For simple binding reactions not the molecule type is changed, but its mobility is set to zero at the given position. A reverse first order dissociation reaction restores the mobility. The corresponding dissociation probability within [t, t + Δt) is P
diss
= k
diss
× Δt for each bound molecule (cf. [54]).
Remarks on reversible reactions
The binding and unbinding process to the cytoskeleton is a diffusion controlled process by itself. The details of the effective rates in diffusion controlled reversible processes have been studied for instance in [55, 56]. In the present study we adjusted the binding and unbinding rate constants to approximately achieve the desired steady state fu. Future work will analyze the details of the diffusion controlled binding, unbinding, and rebinding process.
Parallelization
Parallelization of the simulation is possible and benefits from the fact that all agents are updated simultaneously with a global Δt. Since the mobile agents can overlap with each other their random walk steps could be computed in parallel in order to reduce the computation time (either on a multicore CPU or GPU [13, 47]). The pairwise testing of agents for the reactions however requires to make sure that the following situation is treated correctly: (i) Assume A and B are supposed to react together to form complex C. (ii) At some point three molecules A,B, and another A will be close enough to react. (iii) A parallel implementation could now find two AB-pairs simultaneously, while B can only react with one of the two A molecules.
Review and comparison with other simulation methods
For the given purpose of the simulation environment, i.e. modeling of the cell including a realistic intracellular environment such as a detailed cytoskeleton structure, only agent-based off-lattice methods can be used. Therefore we leave out the spatial Gillespie method as well as derivatives thereof, and also the E-Cell plug-in of Arjunan and Tomita [57].
We also require that the treatment of bimolecular reactions is implemented efficiently and as correctly as possible. As such, the Greens-function reaction dynamics allows to do large steps if reactants are far apart from each other which would suit this need [58]. However, in the crowded environment the next obstacle is always close, which cancels out the advantage of the method.
The distance-dependent reaction probability derived from the Fokker-Planck Equation as outlined above based on Equations (12) and (13) (cf. [51, 52]) requires a look-up table and is therefor slower than a method which uses a fixed critical reaction radius (and if necessary a fixed reaction probability). The reaction method from Ridgway et al. [56] is likewise derived from Equation (12) and (13) but aims at such a simplified description with a critical radius. However it remains unclear how their reaction probability corresponds to the macroscopic reaction rate because it depends on Δx. Likewise the relationship between the macroscopic reaction rate and the reaction parameters in Smoldyn [54] and the cellular dynamics simulator from Byrne et al. [59] is only established by an iterative search or a look-up table. In contrast, the approach of Pogson et al. [53] promises to give a simple answer how to get the reaction radius from the macroscopic reaction rate (cf. the derivation of the present method above). However we found, that it does not work correctly for diffusion controlled reactions. The present method extends Pogson's method by matching the collision and reaction radius which correspond to the microscopic reaction rate. The new method was tested for different time steps and reaction rate constants. The results are in agreement both with the results obtained with the more detailed reaction scheme based on the Fokker-Planck equation [51, 52], as well as with ODE models for well mixed conditions. The only limitation is that the step length Δt should be chosen such that the reaction probability
(cf. Equation 24). Thus we claim that we have found a simple, efficient, and sufficiently accurate description of the reaction-diffusion process in Brownian dynamics simulations for simulations on the cell level.
Quantifying the influence of the reduced diffusion on a bimolecular reaction
The macroscopic reaction rate k
macro
: = k
ij
for ODE-models (see above) is related to the microscopic reaction rate, which states the reaction probability upon a collision between molecules in the detailed simulation, by Equation (3). The corresponding collision rate constant k
D
(which is the diffusion limit of the given reaction) is determined in 3D to [7, 8, 41, 60])
k
D
is accordingly calculated based on the collision distance and the combined diffusion coefficient of the reactants.
For a given macroscopic reaction rate, the microscopic reaction rate constant is accordingly given by Equation (3):
(given that the user does not try to exceed the diffusion limit with the macroscopic reaction rate, i.e. k
macro
< k
D
). Now we are interested in the effective macroscopic reaction rate constant, given that the microscopic reaction rate constant is held constant but the diffusion is reduced in the crowded intracellular conditions. This leads to a reduced
and the effective macroscopic reaction rate can now be calculated based on Equation (3)
Inserting Equation (26) into Equation (28) leads to
If also the definition of k
D
and k
D, eff
are inserted, this becomes
The initial (unperturbed) macroscopic reaction rate can be set into relation with the diffusion limit, defining
which leads to a simplification of Equation (30)
From this equation it can be deduced that the effective macroscopic reaction rate constant is reduced by the factor