Abstract
Background
In this paper we apply a novel agentbased simulation method in order to model intracellular reactions in detail. The simulations are performed within a virtual cytoskeleton enriched with further crowding elements, which allows the analysis of molecular crowding effects on intracellular diffusion and reaction rates. The cytoskeleton network leads to a reduction in the mobility of molecules. Molecules can also unspecifically bind to membranes or the cytoskeleton affecting (i) the fraction of unbound molecules in the cytosol and (ii) furthermore reducing the mobility. Binding of molecules to intracellular structures or scaffolds can in turn lead to a microcompartmentalization of the cell. Especially the formation of enzyme complexes promoting metabolic channeling, e.g. in glycolysis, depends on the colocalization of the proteins.
Results
While the colocalization of enzymes leads to faster reaction rates, the reduced mobility decreases the collision rate of reactants, hence reducing the reaction rate, as expected. This effect is most prominent in diffusion limited reactions. Furthermore, anomalous diffusion can occur due to molecular crowding in the cell. In the context of diffusion controlled reactions, anomalous diffusion leads to fractal reaction kinetics. The simulation framework is used to quantify and separate the effects originating from molecular crowding or the reduced mobility of the reactants. We were able to define three factors which describe the effective reaction rate, namely f ^{diff }for the diffusion effect, f ^{volume }for the crowding, and f ^{access }for the reduced accessibility of the molecules.
Conclusions
Molecule distributions, reaction rate constants and structural parameters can be adjusted separately in the simulation allowing a comprehensive study of individual effects in the context of a realistic cell environment. As such, the present simulation can help to bridge the gap between in vivo and in vitro kinetics.
Background
The complex structured and crowded intracellular conditions [1] have a tremendous impact on intracellular reactions. Accordingly, the in vivo rate constants or even the structure of the kinetic rate expression can significantly differ from those obtained in in vitro assays [2]. First of all, the crowded conditions squeeze all molecules closer together which favors the formation of more compact complexes [35]. Associations or more general bimolecular reactions are governed by the occurrence of collisions of the respective molecules. The frequency of the collisions, in turn, depends on the mobility of the molecules. Molecular crowding and especially the cytoskeleton structure lead to a reduction in the diffusion rate, which depends on the size of the molecules [6]. Via the collision based principle of (diffusionlimited) reactions this also translates into reduced reaction rates [7,8]. In this context, it is also worth noting that anomalous (timedependent) diffusion, which was observed in crowded systems [6,9], leads to timedependent  fractal  reaction rate constants [7,1012].
In order to investigate the effects of a given intracellular state on the reaction rate, we have developed an agentbased simulation which tracks individual molecules through a virtual cell containing a model cytoskeleton (see Figure 1) [6,13,14]. The irregular cellular architecture requires an offlattice continuous space Monte Carlo method in which all structures are modeled explicitly as static obstacles. As long as no active transport e.g. by motor proteins is introduced, the molecules of interest move solely by diffusion, which translates into a random walk in the present simulation. Obviously, steps into the obstacles are prohibited which enforces to a tortuous way of the mobile molecules around the obstacles. The resulting effective diffusion has been explored for example in [6,1517]. Since the molecules still move with their initial 'speed'  just on a detouric way, the measurement of the displacement will return D_{0 }if sampled on short distances/times and a reduced D_{eff }on longer distances. Therefore, D_{eff }is transiently converging to a fixed long time diffusion coefficient. The corresponding crossover time/distance depends inter alia on the level of crowding [9,18].
Figure 1. Intracellular structure. Comparison of the 3D intracellular structures: in vivo (left) and in silico (right) cytoskeleton. The 815 × 870 × 97 nm section shows actin filaments and ribosomes.The in vivo image is reprinted from Medalia et.al. (2002), Science 298:12091213 [14] with permission from AAAS. The structures in the in silico cell are randomly arranged based on a uniform distribution. Actin volume fraction: 8.5% (modeled by cylinders with diameter of 7 nm and length of 700 nm), 'ribosome' volume fraction: 5.8% (modeled by spheres with diameter of 20 nm).
The mobility of the reactants is not the only factor determining the effective in vivo reaction rate. Research on the interactome, describing the reactions between proteins, metabolites, and further biomolecules, has revealed a multitude of interactions for each molecule species [19]. Proteins, for example, can unspecifically bind to cytoskeleton structures, which in turn leads to a further reduction in their mobility [1,20]. The decreased effective diffusion coefficient reduces the collision probability between reactants and can thus reduce the reaction rate if the reaction is significantly diffusion limited.
While the adsorption (for instance to the cytoskeleton) and subsequent immobilization hampers the reactions [21] (in the worst case the active site of the enzyme faces the cytoskeleton and is therefore blocked), it is the prerequisite to assign the proteins to specific locations in the cell. This is of particular importance in the following cases:
• Metabolic channeling [2124] describes the concept of enzyme colocalization along the cytoskeleton [21,24,25], which can be highly regulated [2426]. It has been suggested that metabolites are processed in this channel like in an assembly line, for example in glycolysis (see Figure 2a). Binding to actin filaments can also lead to an allosteric activation of the enzymes [24].
Figure 2. Metabolic channeling and scaffolds in signal transduction. (a) Metabolic channeling in contrast to unconnected metabolic reactions with unbound enzymes [24]. (b) Signal transduction through a localized scaffold versus unbound signaling components that can also activate other pathways (crosstalk)  here Kss1 in addition to Fus3 in yeast Saccharomyces cerevisiae [27].
• Regulation of signal transduction: Cells are subject to many and sometimes contradictory signals. The information is carried from the receptors in the plasma membrane towards the nucleus by signaling molecules. Especially in multistage cascades, for example in MAPK (mitogen activated protein kinase) signaling, this transfer can be regulated by scaffolds. The scaffolds integrate several stages of the cascade in one place (see Figure 2b) [27,28]. Scaffolding proteins can regulate signal transduction [28] and furthermore the crosstalk with other signaling pathways [29]. It was also found that the subcellular localization of the molecules matters in signal transduction [30]. Signaling molecules can likewise bind to actin filaments, which was reported for NFκB and its inhibitor IκB [31]. Thus the molecules are sequestered from the cytosol. Only with the right fraction of unbound molecules the correct signal is transmitted [32].
• Pharmacokinetics and drug detoxification: If drug molecules bind to proteins or membranes, they are likewise sequestered from the cytoplasm. This reduces both their action and their degradation, for example by enzymes of the CYP family in the liver [33,34].
All these details are omitted in models based on differential equations in which only the number/concentration of the species is tracked. The general compartmentalization of the organism/cell can be included in the model if the respective compartments and transport rates are defined [35]. Spatial aspects like the transport from the plasma membrane towards the nucleus in case of signaling molecules were also investigated with partial differential equations [36,37]. However, all approaches based on differential equations are based on (local) mean values and neglect the detailed, microscopic aspects in the cell, for instance a microcompartmentalization along the cytoskeleton.
A particle or agent based simulation allows exploring the effects introduced by the spatial organization of the cell including (transient) binding to the cytoskeleton. The in silico simulation environment also enables to change just one parameter of the complex setup for a comprehensive study of the individual crowding, structuring, or binding effects within a realistic environment. This paper presents the results of a Brownian dynamics particle based simulation covering diffusion and reaction rates in a virtual cell. The effects of the cytoskeleton and transient binding on the mobility of tracer molecules are evaluated, and the respective effective reaction rates are analyzed. Eventually, different spatial distributions of enzymes in the cell are tested in order to compare the effect of the formation of enzyme channels with homogeneously distributed enzymes. The simulation framework is described in the Methods section.
Results and discussion
The simulations are performed in a small model cell with a diameter of 7 μm (see Additional file 1, Section 1). The cytoskeleton structure is created by 25,000 randomly arranged cylinders with a length of 2.5 μm and a diameter of 35 nm. In addition, 100,000 immobile crowding spheres with a diameter of 60 nm are placed in the cell. Cylinders and spheres together occupy 24% of the volume in agreement with experimental results [38]. For tracer molecules with a radius of 2.5 nm, the excluded volume fraction ∈ increases to 30.5% due to their own volume. The respective volume fractions were sampled using a MonteCarlo testing method.
Additional file 1. Contains further simulation results and details of the model setup. Simulation. The simulation is available upon request from Michael Klann, http://mklann@ee.ethz.ch webcite.
Format: PDF Size: 2.7MB Download file
This file can be viewed with: Adobe Acrobat Reader
Effective diffusion
The mean squared displacement (MSD) of diffusing molecules should increase linearly with time according to
where d is the given dimension [16]. The effective diffusion coefficient in a given intracellular condition can hence be calculated from the resulting displacement. The effective diffusion for tracer molecules with a radius of 2.5 nm through the present sample cell was accordingly calculated from the MSD in the simulation to D_{eff }/D_{0 }= 0.77 ± 0.01. This slowdown is in agreement with previous studies of the impact of the cytoskeleton structure on the diffusion of inert (i.e. molecules that do not interact with other molecules) tracer molecules for the given excluded volume fraction [6,15].
The slowdown of the diffusion can also be explained by the local confinement of the obstacles [39]. In our model cell the random cytoskeleton/crowding structures create a multitude of randomly arranged confinement boundaries. The uniformly distributed test molecules in the cell sample the average effect of these boundaries in their joint MSD measure. Initially, a nonlinear diffusion can be observed, due to the fact that most particles can move a few steps before they hit an obstacle [6]. Accordingly the MSD grows proportional to the original diffusion coefficient D_{0 }in the beginning, and later on proportional to the average D_{eff}.
If tracers bind transiently to the cytoskeleton and are therefore temporarily immobilized, the effective diffusion coefficient is further reduced (see Figure 3a and Figure 3b). The simulation shows that this reduction is proportional to the steady state fraction of unbound molecules (fu: = unbound molecules/total molecules of the respective species):
Figure 3. Transient binding. (a) Transient binding of molecules to the cytoskeleton. (b) Diffusion rises proportionally to the fraction of unbound molecules. (c) Nonlinearity in 〈r^{2}〉 induced by the dynamics of the reversible binding process. k_{diss }= 2.5 1/s; 10 1/s; 25 1/s; 100 1/s; from top to bottom, the corresponding. k_{binding }= 2.6 1/s; 10.4 1/s; 26 1/s; 104 1/s; leads to fu = 0.51.
Figure 3c shows, how the equilibrium fu develops if initially all molecules are unbound. The mean squared displacement of the molecules shows a nonlinear behavior during this transition phase, which depends on the rate constants of the binding/dissociation reaction.
Effective reaction rates
Diffusioncontrolled reactions
The theory of diffusion controlled reactions requires to take into account the following points [40]:
• Diffusion Limit: The maximal reaction rate constant for a bimolecular reaction of two spherical molecules i and j with radius r_{i }and r_{j }is: k_{D }= 4π(r_{i }+ r_{j})(D_{i }+ D_{j}) (in 3D) [41]. It equals the collision rate of the molecules. (The fact that initially some nearby pairs will react faster leads to an initially time dependent reaction rate constant .)
• Microscopic Reaction Rate Constant: If not every encounter between two reactants leads to a reaction, the microscopic reaction rate constant k_{micro }determines the fraction of collisions which lead to subsequent reactions.
• Effective Macroscopic/Bulk Reaction Rate Constant: The resulting reaction rate constant which is observed on the macroscopic level, corresponding to the rate constant of ODE models is determined as [40]:
Test setup
The test molecules in the simulation are enzyme E and substrate S molecules with the following properties: radius r_{E }= r_{S }= 2.5 nm and diffusion coefficient of D = 1 μm^{2}/s. This leads to a diffusion limit of the reaction rate of k_{D }= 7.57 × 10^{7 }l/(mol·s). The chosen macroscopic reaction rate for a test simulation can be given as a fraction of k_{D }allowing a dimensionless survey of effects on the reaction rate. In the following, rates of k_{macro }= 7.57 × 10^{5 }l/(mol·s) (1% of k_{D}), k_{macro }= 7.57 × 10^{6 }l/(mol·s) (10% of k_{D}), and k_{macro }= 2.27 × 10^{7 }l/(mol·s) (30% of k_{D}) are tested.
The resulting reaction rate can be accessed from the change in the number of molecules. The noise in a stochastic simulation, however, hampers the identification of the current reaction rate. If the considered species is, in turn, created and destroyed by two reactions, it will accumulate to a dynamic equilibrium steady state. Averaging the steady state number over time reduces the stochastic noise in the result. This situation can be found in vivo in the sequence of enzymatic reactions, for example in glycolysis as shown in Figure 4.
Figure 4. Model setup. Description of an enzymatic reaction in a metabolic pathway based on mass action kinetics.
In order to reduce the complexity, the considered substrate species S is created in a zero order reaction with rate constant k_{1}. It is consumed in the enzymatic reaction S + E → P + E, which is modeled here with mass action kinetics based on the rate constant k_{2}. The number of enzymes E (concentration c_{E}) is not affected by this reaction. c_{E }is set to 2 × 10^{7 }mol/l (20600 molecules). The macroscopic balance equation for the substrate concentration is in this model
which leads to the equilibrium steady state
Detailed simulation vs. ODEmodel
This section compares the outcome of the detailed stochastic simulation with the result of the macroscopic ODEmodel of Figure 4. In order to elucidate the differences between diluted in vitro and crowded in vivo conditions, one simulation is conducted in an empty 'in vitro' virtual cell and one in the crowded model cell described above.
The steady state of the 'in vitro' reaction rate for both the predicted (Equation (5)) and simulated molecule numbers is given in Table 1, showing that the simulation is able to reproduce the macroscopic reaction rates correctly under the 'in vitro' conditions.
Table 1. Results for the 'in vitro' setup
The situation is quite different in the 'in vivo' case. While r_{1 }= k_{1 }is held constant in the simulation, the bimolecular reaction is affected by the crowded intracellular conditions. The rate for the second reaction becomes
The steady state shifts accordingly to
In order to understand the corresponding change in the reaction rate, the overall effect (f ^{eff}) is broken down into the following three factors:
1. The first factor arises from the reduced free volume fraction , which leads to an increased effective concentration of the reactants (given that c_{0 }is calculated as number of molecules per cell/total cell volume). This factor has to be added only once (for c_{E}) in the mass action reaction framework dc_{S}/dt = k_{2}c_{E}c_{S }because c_{S }appears both on the right and the left side of the equation. Instead of using effective concentrations, the respective factor can also be applied to the reaction rate constant, which leads to an apparent reaction rate of . Accordingly the in vitro reaction rate has to be multiplied with a factor
2. The reduced effective diffusion has an influence on the reaction rate because it reduces the collision rate. For the present analysis it is assumed that the molecules react in the same way in vivo and in vitro, i.e. the microscopic reaction rate constant (which describes the reaction probability upon a collision) stays the same (cf. Equation (3). The effect of the reduced diffusion on the reaction rate can be calculated using β: = k_{2}/k_{D }(see Methods section), leading to
3. The hindered accessibility of the molecules due to steric effects of nearby obstacles contributes a further reduction f ^{access }of the reaction rate (see Figure 5 for an explanation). Using a MonteCarlo sampling method of the respective volume fraction this factor was estimated to f ^{access }= 0.966 ± 0.001 in the given virtual cell.
Figure 5. Inaccessible volume fraction. The restricted volume close to all structures in the cell reduces the interaction volume (green) for the reaction. In order to estimate the effect of the reduced interaction volume on the reaction rate, the fraction of the accessible reaction volume has to be averaged over all possible molecule positions in the given cell. In the complex environment of the given random intracellular architecture the calculation of the corresponding f ^{access}factor is only possible with a MonteCarlo sampling method, averaging the accessible (green) volume fraction of the interaction volume of all possible molecular positions in the cell.
In combination the effective macroscopic bimolecular reaction rate is accordingly:
Table 2 contains the resulting steady state molecule numbers of three simulations with β: = k_{2}/k_{D }= 0.01, β = 0.1 and β = 0.3 as well as the model prediction for the virtual cell based on Equation (7) and the steady state molecule numbers of the 'in vitro' simulation . While the model prediction and the simulation results are in a perfect agreement in the case without diffusion limitation (β = 0.01), the values for β ≥ 0.1 show a significant and increasing deviation. The number of molecules in the in vivo simulation is smaller than predicted by the ODE model  which means that the reaction rate r_{2 }is faster. For comparison also a simulation in a homogenized cell was conducted. This cell does not contain any hindering obstacles but the size is reduced by a factor of 0.695 so that the effective concentration of molecules matches the effective concentration in the detailed virtual cell. Also the diffusion is reset in order to match the respective effective diffusion  but only after the microscopic reaction rate was set based on the in vitro diffusion coefficient. The model prediction and the latter simulation show a good agreement (see Table 2). This leads to the conclusion that in the detailed and crowded virtual cell the local properties differ from the average properties, and that the reaction rate depends on the local effective diffusion. In turn, the reaction rate could also be used to probe the local effective diffusion.
Table 2. Results for the 'in vivo' setup
In order to understand this result, it is necessary to recall the transient anomalous diffusion in the crowded environment. At short distances, the molecules still move with their original (fast) diffusion coefficient. Only on longer distances the tortuous way around the obstacles leads to a reduced mobility. The results indicate that the diffusion limited bimolecular reaction senses an intermediary effective diffusion coefficient which is slower than D_{0 }but faster than the long term effective D_{eff }(cf. Additional file 1, Section 2). This argument is supported by the stronger deviation of the result of the more diffusion limited reaction k_{2}/k_{D }= 0.3.
Future work could investigate this effect with respect to the local confinement, and also include the influence of unspecific and transient binding on the reaction rate, i.e. the nullified mobility of one or both of the reactants due to binding to cellular structures. In addition, also the combined action of individual reaction rate constants for different substates of a molecule (free, bound, phosphorylated at site x, etc.) could be analyzed in a more complex model.
Enzyme colocalization and metabolic channeling
This section considers more than one reaction of the enzymatic conversion of the metabolites in the cell. The pathway is simplified to a sequence of 4 reactions as shown in Figure 6a. First the substrate (e.g. glucose) is transported across the plasma membrane, then it is converted enzymatically, and finally the (virtual) product is transported out of the cell (like e.g. lactate).
Figure 6. Metabolic channeling model. (a) Metabolic pathway and (b) possible enzyme distributions. Note, that the enzymes are immobilized (D_{E }= 0) in the simulation in order to maintain the distribution. (c) Comparison of the intracellular metabolite 5 level for different enzyme distributions and different diffusion coefficients of the metabolites (left: D = 1 μm^{2}/s; right: D = 10 μm^{2}/s)).
The aim of the present simulations is to elucidate the effect of the localization of the enzymes. The most optimal setup promises to be the enzyme channel (A) of Figure 6b in which the enzymes are colocalized (cf. Figure 2) and aligned in the order of the reactions. Since the substrate enters the cell through the plasma membrane, the enzymes should be localized there (along actin filaments which mainly polymerize next to the plasma membrane [42]). Since it is not clear, whether the enzymes arrange in a sequential complex or are just randomly attached to the actin filaments [21,24,25] close to the surface, also a dense but unstructured enzyme layer (B) is modeled for comparison. These structured setups are furthermore compared to a uniform random distribution (C) of the enzymes in the cell and a well mixed (D) model based on ODE (in order to keep the models comparable, the corresponding stochastic solution of the ODE model is evaluated based on the Gillespie method [43]). This means that the spatial aspects are completely neglected and only the molecule numbers are tracked in (D).
All enzymes are immobilized in the simulation, i.e. they stay at their fixed initial position with an artificial D_{E }= 0. Therefore it is not necessary to include the cytoskeleton filaments as structuring elements as well as the enzyme binding to the cytoskeleton, which would increase the complexity of the model. This nicely underscores the advantage of in silico simulations in which the different factors affecting the reaction rates in the cell can be separated.
The metabolites are moving with a diffusion coefficient of (i) D = 1 μm^{2}/s and (ii) D = 10 μm^{2}/s for a comparison of the mobility effects. The macroscopic reaction rate constant is set to k = 3.78 × 10^{6 }l/(mol·s), which is fairly fast but not extremely diffusion limited (k_{D }= 3.78 × 10^{7 }l/(mol·s) for D = 1 μm^{2}/s, and k_{D }= 3.78 × 10^{8 }l/(mol·s) for D = 10 μm^{2}/s)  see Additional file 1, Section 3.
Since all setups are conducted with the same number of enzymes and the same reaction rates, they produce similar results (see Figure 6c). The deviations are stronger for the more diffusion limited case  i.e. the lower diffusion D = 1 μm^{2}/s. The overall development of the individual metabolite pools is shown in the Additional file 1, Section 3. It is worth noting that the export rate of metabolite 5 is diffusion controlled in the setups (A)(C). Especially in the randomly distributed setup (C), metabolite 5 can be created deep inside the cell and has to diffuse all the way back to the plasma membrane. This leads to a stronger accumulation of the metabolite in the cell compared to the wellmixed ODEmodel (D).
Likewise the setups where the enzymes are close to the plasma membrane (i.e. channel (A) or layer (B)) lead to a faster formation of the product because the substrate enters the cell through the plasma membrane (note, that this setup also promotes a faster export of the metabolite because it is produced next to the plasma membrane). This is clearly visible in the excerpt of the initial phase shown in Figure 6c. Due to the optimal location of the enzymes, the product formation is even faster than in the wellmixed ODEmodel. Locally the enzyme concentration is much higher than the average concentration in the cell  right where it is needed. Interestingly, the product formation rate is much faster for the slower diffusion  in this case the metabolites stay longer in the vicinity of the plasma membrane where they can interact with the enzymes. As such, the cell has become compartmentalized although no explicit compartments are defined.
The setup where the respective enzymes are colocalized in an enzyme channel indeed shows the fastest (initial) product formation rate. It can be expected that the differences between the enzyme distributions tested in this paper will increase if the MichaelisMenten enzyme kinetics is used. The high local enzyme concentration close to the surface leads to a locally higher V_{max }value  right were the highest substrate concentration is found in the given setup. Accordingly, the reactions in the channeled and layered setups should become even faster.
The improved reaction rate in the colocalized channel structure is in agreement with the findings of Bauler et al. [44], which found that a sequential reaction is faster if the two active zones of the respective enzymes are aligned to face each other in a Brownian dynamics simulation. However, it should be noted that the faster reaction rate arises solely from the fact that the intermediate metabolites have a higher probability to hit the next enzyme within the next diffusion steps in the simulation but are still allowed to diffuse away. In reality, the metabolites might be directed through the entire enzyme complex, which leads to a completely different description of the overall reaction. The analysis of the latter effect, in turn, requires a molecular dynamics study based on the molecular multi enzyme structure of the proposed channeling complex. Experimental observations [21,22,45] indicate that indeed enzyme channels might be formed in vivo. However, the theoretical analysis did not yet identify channeling as the only possibility to describe the observed kinetics [21].
The present study shows that the location of the reactants (here enzymes and metabolites) can play an important role. This work focused on the influence of spatial aspects and the possible enzyme colocalization, which allows a more realistic study of enzyme channeling properties. Future work could include more advanced reaction kinetics in order to verify channeling, as well as a study of the resulting control properties on the metabolic flux [45].
Conclusions
The present simulation allows a detailed analysis of the effects of the intracellular properties on the reaction rates. The results have shown that the in vivo reaction rate differs from the in vitro reaction rate. The difference is related to the crowded conditions in the cell and three factors, namely (i) the increased effective concentration, (ii) the reduced accessibility of the reactants, and (iii) the reduced mobility of the reactants could be determined.
In addition, the influence of the subcellular localization of the reactants was tested. The results show that the colocalization of enzymes in a metabolicchanneling framework can improve the product formation. It is worth noting that the advantage of a specific location of the enzymes is accompanied by the disadvantage of the reduced enzyme mobility. Hence the reaction rate will be reduced in the diffusion limited case. This reduction could even outbalance the superior product formation rate of enzyme channeling. Since this depends on the actual diffusion and reaction rate constants, further simulations are required in order to quantify the advantage of the channel configuration  especially in the context of a more advanced kinetics within the multienzymecomplex.
Thus the present simulation framework is a promising tool to investigate intracellular reactions and signal transduction processes in the detailed spatial organization of the cell [46]. If all intracellular factors are put together correctly, the simulation will return a prediction for the in vivo reaction rate. As such, the present simulation method enables to reconstruct the in vivo properties in silico.
On the way to a model which is in agreement with living cells, several parameters like the correct cytoskeleton structure, molecular crowding, and additional unspecific interactions which can for example transiently immobilize the molecules have to be adjusted, giving a deeper insight into the cell [5]. The resulting detailed in silico cell can also be visualized by advanced and interactive visualization tools, thus providing a powerful virtual microscope [13,47,48].
Methods
Description of the agentbased 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 D_{0 }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 t_{n+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 FokkerPlanck 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 r_{0}. 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 eventbased stochastic framework can be derived from the macroscopic (ordinary differential equation) description
For simplicity, the units of the concentrations should be [molecules/μm^{3}] (note: the units of k_{ij }are then [μm^{3}/molecules/s]. If c_{i }is given in [mol/liter], the factor N_{A}/10^{15 }[liter/(μm^{3}mol)], 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 [reactionvolume/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 ∝ Δt^{1/3 }according to Equation (20), compensating the Δtdependency 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 nonoverlapping 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 volumebased framework developed above.
• Solution: reaction probability in the interaction volume: Both concepts, the flux/surfacebased description and the macroscopic, volumebased 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:
 the corresponding reaction volume should be in analogy to Equation(18)
but this will (most likely) not match the collision volume 4π = 3(r_{i }+ r_{j})^{3}
 The mismatch is adjusted by introducing the reaction probability
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.
• Resulting reaction algorithm: Two particles i and j will react if the distance between them is smaller than and a random number of the interval [0,1] is smaller than (which on average leads to a reaction with the probability ).
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 ABpairs 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 agentbased offlattice methods can be used. Therefore we leave out the spatial Gillespie method as well as derivatives thereof, and also the ECell plugin 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 Greensfunction 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 distancedependent reaction probability derived from the FokkerPlanck Equation as outlined above based on Equations (12) and (13) (cf. [51,52]) requires a lookup 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 lookup 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 FokkerPlanck 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 reactiondiffusion 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 ODEmodels (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
Authors' contributions
MK developed, designed, and performed the simulations and drafted the manuscript. AL calculated the reaction probability based on the FokkerPlanck equation and revised the core functionality of the simulation. MR is the group leader, has initiated the program and revised the manuscript. All authors read and approved the final manuscript.
Acknowledgements
The authors would like to thank Martin Falk and Thomas Ertl for the valuable discussions about the visualization method as well as on the possible parallelization of the simulation and acknowledge the funding by the State of BadenWürttemberg/Center Systems Biology in Stuttgart.
References

LubyPhelps K: Cytoarchitecture and Physical Properties of Cytoplasm: Volume, Viscosity, Diffusion, Intracellular Surface Area.
Int Rev Cytol 2000, 192:189221. PubMed Abstract

De Silva A, Fraenkel D: The 6Phosphogluconate Dehydrogenase Reaction in Escherichia coli.
Journal of Biological Chemistry 1979, 254(20):1023710242. PubMed Abstract  Publisher Full Text

Ellis R: Macromolecular crowding: an important but neglected aspect of the intracellular environment.
Curr Opin Struct Biol 2001, 11:114119. PubMed Abstract  Publisher Full Text

Minton AP: The Inuence of Macromolecular Crowding and Macromolecular Confinement on Biochemical Reactions in Physiological Media.
J Bio Chem 2001, 276(14):1057710580. Publisher Full Text

Minton A: How can biochemical reactions within cells differ from those in test tubes?
Journal of Cell Science 2006, 119(14):2863. PubMed Abstract  Publisher Full Text

Klann M, Lapin A, Reuss M: Stochastic Simulation of Signal Transduction: Impact of the Cellular Architecture on Diffusion.
Biophys J 2009, 96(12):51225129. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Torney D, McConnel H: DiffusionLimited Reaction Rate Theory for TwoDimensional Systems.
Proc R Soc Lond A 1983, 387(1792):147170. Publisher Full Text

Rice S: Diffusionlimited reactions. Elsevier Amsterdam; 1985.

Weiss M, Elsner M, Kartberg F, Nilsson T: Anomalous Subdiffusion Is a Measure for Cytoplasmic Crowding in Living Cells.
Biophys J 2004, 87:35183524. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Schnell S, Turner T: Reaction kinetics in intracellular environments with macromolecular crowding: simulations and rate laws.
Progress in biophysics and molecular biology 2004, 85(23):235260. PubMed Abstract  Publisher Full Text

Grima R, Schnell S: A systematic investigation of the rate laws valid in intracellular environments.
Biophys chem 2006, 124:110. PubMed Abstract  Publisher Full Text

Haugh J: Analysis of ReactionDiffusion Systems with Anomalous Subdiffusion.
Biophys J 2009, 97(2):435442. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Falk M, Klann M, Reuss M, Ertl T: Visualization of Signal Transduction Processes in the Crowded Environment of the Cell.
Proceedings of IEEE Pacific Visualization Symposium 2009 (PacificVis '09) 2009, 169176.

Medalia O, Weber I, Frangakis AS, Nicastro D, Gerisch G, Baumeister W: Macromolecular Architecture in Eukaryotic Cells Visualized by Cryoelectron Tomography.
Science 2002, 298:12091213. PubMed Abstract  Publisher Full Text

Novak I, P K, Slepchenko B: Diffusion in Cytoplasm: Effects of Excluded Volume Due to Internal Membranes and Cytoskeletal Structures.
Biophys J 2009, 97(3):758767. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Trinh S, Arce P, Locke BR: Effective Diffusivities of PointLike Molecules in Isotropic Porous Media by Monte Carlo Simulation.
Transport in Porous Media 2000, 38:241259. Publisher Full Text

Netz P, Dorfmüller T: Computer simulation studies of anomalous diffusion in gels: Structural properties and probesize dependence.
J Chem Phys 1995, 103(20):90749082. Publisher Full Text

Chang R, Jagannathan K, Yethiraj A: Diffusion of hard sphere fluids in disordered media: A molecular dynamics simulation study.

Welch G: The 'fuzzy' interactome.
Trends in Biochemical Sciences 2009, 34:12. PubMed Abstract  Publisher Full Text

Saxton M: A Biological Interpretation of Transient Anomalous Subdiffusion. I. Qualitative Model.
Biophys J 2007, 92:11781191. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Welch G, Easterby J: Metabolic channeling versus free diffusion: transitiontime analysis.
Trends in biochemical sciences 1994, 19(5):193197. PubMed Abstract  Publisher Full Text

Winkel B: Metabolite Channeling and Multienzyme Complexes.
Plantderived Natural Products: Synthesis, Function, and Application 2009, 195208.

Ovadi J, Sreret P: Macromolecular compartmentation and channeling.

AlHabori M: Microcompartmentation, metabolic channelling and carbohydrate metabolism.
International Journal of Biochemistry and Cell Biology 1995, 27(2):123132. PubMed Abstract  Publisher Full Text

Masters C: Interactions between glycolytic enzymes and components of the cytomatrix.
Journal of Cell Biology 1984, 99:222225. Publisher Full Text

Ovadi J: Old pathwaynew concept: control of glycolysis by metabolitemodulated dynamic enzyme associations.
Trends in biochemical sciences 1988, 13:486450. PubMed Abstract  Publisher Full Text

Klipp E, Liebermeister W: Mathematical modeling of intracellular signaling pathways.
BMC neuroscience 2006, 7(Suppl 1):S10. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

McKay M, Ritt D, Morrison D: Signaling dynamics of the KSR1 scaffold complex.
PNAS 2009, 106(27):1102211027. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Force T, Woulfe K, Koch W, Kerkela R: Molecular Scaffolds Regulate Bidirectional Crosstalk Between Wnt and Classical SevenTransmembrane Domain Receptor Signaling Pathways.

Harding A, Tian T, Westbury E, Frische E, Hancock J: Subcellular localization determines MAP kinase signal output.
Current Biology 2005, 15(9):869873. PubMed Abstract  Publisher Full Text

Are A, Galkin V, Pospelova T, Pinaev G: The p65/RelA Subunit of NFκB Interacts with ActinContaining Structures.
Experimental Cell Research 2000, 256:533544. PubMed Abstract  Publisher Full Text

Pogson M, Holcombe M, Smallwood R, Qwarnstrom E: Introducing Spatial Information into Predictive NFκB Modelling  An AgentBased Approach.
PLoS ONE 2008, 3(6):e2367. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Paine S, Parker A, Gardiner P, Webborn P, Riley R: Prediction of the pharmacokinetics of atorvastatin, cerivastatin, and indomethacin using kinetic models applied to isolated rat hepatocytes.
Drug Metabolism and Disposition 2008, 36(7):13651374. PubMed Abstract  Publisher Full Text

Gao H, Steyn S, Chang G, Lin J: Assessment of in silico models for fraction of unbound drug in human liver microsomes.
Expert Opinion on Drug Metabolism & Toxicology 2010, 2010(5):533542. PubMed Abstract  Publisher Full Text

Maier K, Hofmann U, Bauer A, Niebel A, Vacun G, Reuss M, Mauch K: Quantification of statin effects on hepatic cholesterol synthesis by transient 13Cflux analysis.

Kholodenko B: Fourdimensional organization of protein kinase signaling cascades: the roles of diffusion, endocytosis and molecular motors.
J Exp Biol 2003, 206(12):20732082. PubMed Abstract  Publisher Full Text

Markevich N, Tsyganov M, Hoek J, Kholodenko B: Longrange signaling by phosphoprotein waves arising from bistability in protein kinase cascades.

LubyPhelps K: Cytoarchitecture and Physical Properties of Cytoplasm: Volume, Viscosity, Diffusion, Intracellular Surface Area.
Int Rev Cytol 2000, 192:189221. PubMed Abstract

Hall D, Hoshino M: Effects of macromolecular crowding on intracellular diffusion from a single particle perspective.
Biophys Rev 2010, 2:3953. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Berg O, von Hippel P: Diffusioncontrolled macromolecular interactions.
Annual Review of Biophysics and Biophysical Chemistry 1985, 14:131158. PubMed Abstract  Publisher Full Text

Smoluchowski M: Versuch einer mathematischen Theorie der Koagulationskinetik kolloider Lösungen.

Alberts B, Johnson A, Lewis J, Raff M, Roberts K, Walter P: Molekularbiologie der Zelle. WILEYVCH Verlag; 2004.

Gillespie D: A General Method for Numerically Simulating the Stochastic Time Evolution of Coupled Chemical Reactions.
J Comp Phys 1976, 22:403434. Publisher Full Text

Bauler P, Huber G, Leyh T, McCammon J: Channeling by proximity: the catalytic advantages of active site colocalization using Brownian dynamics.
The Journal of Physical Chemistry Letters 2010, 1(9):13321335. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Rohwer J, Postma P, Kholodenko B, Westerhoff H: Implications of macromolecular crowding for signal transduction and metabolite channeling.
PNAS 1998, 95(18):1054710552. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Takahashi K, Arjunan S, Tomita M: Space in systems biology of signaling pathways  towards intracellular molecular crowding in silico.
FEBS Lett 2005, 579:17831788. PubMed Abstract  Publisher Full Text

Falk M, Klann M, Reuss M, Ertl T: 3D visualization of concentrations from stochastic agentbased signal transduction simulations.
Biomedical Imaging: From Nano to Macro, 2010 IEEE International Symposium on, IEEE 2010, 13011304.

Shillcock J: Insight or illusion? Seeing inside the cell with mesoscopic simulations.
HFSP Journal 2008, 2:16. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Tolle D, Le Novere N: ParticleBased Stochastic Simulation in Systems Biology.
Curr Bioinf 2006, 1(3):315320. Publisher Full Text

Fox RO: Computational Models for Turbulent Reacting Flows. Cambridge University Press; 2003.

Lapin A, Klann M, Reuss M: Stochastic Simulations of 4Dspatial Temporal Dynamics of Signal Transduction Processes.
Proceedings of the FOSBE 2007, Stuttgart, Germany 2007, 421425.

Lapin A, Klann M, Reuss M: MultiScale SpatioTemporal Modeling: Lifelines of Microorganisms in Bioreactors and Tracking Molecules in Cells.
In Biosystems Engineering II, Advances in Biochemical Engineering Biotechnology Edited by Wittmann C, Krull R. 2010., 21

Pogson M, Smallwood R, Qwarnstrom E, Holcombe M: Formal agentbased modelling of intracellular chemical interactions.
Biosystems 2006, 85:3745. PubMed Abstract  Publisher Full Text

Andrews S, Bray D: Stochastic simulation of chemical reactions with spatial resolution and single molecule detail.
Phys Biol 2004, 1:137151. PubMed Abstract  Publisher Full Text

Morelli M, Ten Wolde P: Reaction Brownian dynamics and the effect of spatial fluctuations on the gain of a pushpull network.
J chem phys 2008, 129:054112. PubMed Abstract  Publisher Full Text

Ridgway D, Broderick G, LopezCampistrous A, Ru'aini M, Winter P, Hamilton M, Boulanger P, Kovalenko A, Ellison M: Coarsegrained molecular simulation of diffusion and reaction kinetics in a crowded virtual cytoplasm.
Biophys J 2008, 94(10):37483759. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Arjunan S, Tomita M: A new multicompartmental reactiondiffusion modeling method links transient membrane attachment of E. coli MinE to Ering formation.
M Syst Synth Biol 2010, 4:3553. Publisher Full Text

van Zon J, Ten Wolde P: Simulating biochemical networks at the particle level and in time and space: Green's function reaction dynamics.

Byrne M, Waxham M, Kubota Y: Cellular Dynamic Simulator: An Event Driven Molecular Simulation Environment for Cellular Physiology.
Neuroinformatics 2010, 8(2):6382. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Collins F, Kimball G: Diffusioncontrolled reaction rates.
J Colloid Sci 1949, 4(4):425437. Publisher Full Text