Next Article in Journal
Mitogen-Activated Protein (MAP) Kinases in Plant Metal Stress: Regulation and Responses in Comparison to Other Biotic and Abiotic Stresses
Next Article in Special Issue
Homology Modeling and Analysis of Structure Predictions of the Bovine Rhinitis B Virus RNA Dependent RNA Polymerase (RdRp)
Previous Article in Journal
Plasma Depolymerization of Chitosan in the Presence of Hydrogen Peroxide
Previous Article in Special Issue
Target Molecular Simulations of RecA Family Protein Filaments
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Review

Spatial Simulations in Systems Biology: From Molecules to Cells

Automatic Control Laboratory, ETH Zurich, Physikstrasse 3, 8092 Zurich, Switzerland
*
Authors to whom correspondence should be addressed.
Int. J. Mol. Sci. 2012, 13(6), 7798-7827; https://doi.org/10.3390/ijms13067798
Submission received: 2 May 2012 / Revised: 8 June 2012 / Accepted: 12 June 2012 / Published: 21 June 2012
(This article belongs to the Special Issue Advances in Biomolecular Simulation)

Abstract

:
Cells are highly organized objects containing millions of molecules. Each biomolecule has a specific shape in order to interact with others in the complex machinery. Spatial dynamics emerge in this system on length and time scales which can not yet be modeled with full atomic detail. This review gives an overview of methods which can be used to simulate the complete cell at least with molecular detail, especially Brownian dynamics simulations. Such simulations require correct implementation of the diffusion-controlled reaction scheme occurring on this level. Implementations and applications of spatial simulations are presented, and finally it is discussed how the atomic level can be included for instance in multi-scale simulation methods.

1. Introduction

At the beginning of this century Tomita [1] wrote the article “Whole-cell simulation: a grand challenge of the 21st century”. With the present contribution, we want to review how far we progressed after one decade. Clearly, the cell is not an unstructured bag of enzymes but highly organized in space [2]. Some structures or compartments can be readily seen under the microscope [3], others emanate more subtle, for instance by the clustering of molecules leading to a microcompartmentalization of the cell [46]. This organization gives rise to spatial-temporal dynamics in the cell [79].
Depending on the chosen description of the cell, the cellular state is encoded by the amount of molecules of each molecular species, but also by the distribution of molecules in the cell and the actual state (e.g., phosphorylation level) or conformation of each molecule [3]. Thus also on the molecular level the spatial organization of the atoms or the electron density respectively controls the state of the cell [13]. If we want to track how a conformational change of a single molecule can influence the state of the cell, we therefore have to bridge scales from 10−10 to 10−5 m and 10−12 to 103 s [14].
For each individual level many simulation methods exist, as reviewed in [1519]. For instance, the overall temporal dynamics in the cell can be simulated with ordinary differential equations (ODE), and the spatial component can be included if partial differential equations (PDE) are used [7]. In order to include stochasticity due to low particle numbers, stochastic differential equations (SDE) have to be used, and Gillespie [20] developed an efficient way to simulate the outcome of the chemical master equation, describing the evolution of the population of molecule species. On the molecular level, the dynamics can be investigated in detail using molecular dynamics simulations (MDS) [21,22].
Despite the increasing computational power of workstations and supercomputers, simulation of the whole cell at the atomic level remains prohibitively expensive. For example a small yeast cell contains “only” 50 million proteins [23], however there are more relevant molecules (DNA, RNA, lipids, metabolites), and especially all the ions and water molecules in the cell. In addition, each macromolecule itself consists of thousands of atoms, such that the number of states for each object in the simulation has to be reduced for a simulation [24]. The cell could for instance be simulated at a more coarsely grained level, while only relevant parts could still be tracked on the atomic level in a multi-scale simulation [25].
Particle-based methods can be a useful tool to analyze cellular dynamics. The particles in these methods can represent atoms as in MDS, molecular subunits [26,27], whole molecules [28], sub-volumes of the cell [29,30] or even a whole cell each [31,32]. The particle interactions and actions then have to be defined according to the chosen level of granularity. Especially if the actions of the objects in the simulation are determined by their own internal state or dynamics, e.g., if they represent individual cells, the term cellular-automata or agent-based becomes more appropriate [33].
In the present review we first describe a particle-based mesoscopic level which can serve as a basis for whole cell simulations (cf. Figure 1). Then we show possible applications of such simulations and finally we discuss how it can be extended in several directions in a multi-scale approach.

2. The Mesoscale Level

In order to track the dynamics for instance of signaling molecules in the cell, each molecule is represented by a particle in the simulation. This level is often called mesoscale, because it is between the microscopic atomic resolution of MDS and the macroscopic resolution of ODE/SDE/PDE methods. The following sections describe how the molecules move and interact in the cell at this level.

2.1. Diffusion in the Cell

Obviously the solvent cannot be included explicitly in a whole cell simulation [34]. Langevin dynamics describe the motion of particle j (position vector xj) with mass Mj in such a case based on [22]
M j d 2 x j d t 2 = - U ( x j ) - γ j M j d x j d t + 2 k B T γ j M j R ( t )
within the interaction potential U(x), at temperature T and subject to friction γj with the Boltzmann constant kB and (zero mean) white noise vector R(t), which represents the force induced by collisions with solvent molecules. According to the fluctuation-dissipation theorem the energy added by R is dissipated by γj such that the system reaches and fluctuates around T [22].
The overdamped case (i.e., the damping is so strong that inertia can be neglected) reduces Langevin dynamics to Brownian dynamics
d x j d t = - U ( x j ) / ζ j + 2 D j R ( t )
where the damping is expressed by ζj = γjMj and the diffusion coefficient Dj = kBT/ζj. Especially if the particles are connected, e.g., because they represent sub-segments of a cellular filament, great care has to be taken on the definition of the corresponding interaction potential U between them [26,27,31,35].
The authors understand the importance of all intramolecular forces, and especially electrostatic and hydrodynamic interactions [12,3638]. Still, if all forces are neglected, the Brownian dynamics motion can be easily integrated with a discrete time EulerMaruyama scheme as a random walk [39]:
x j ( t + Δ t ) = x j ( t ) + 2 D j Δ t ξ
with ξ a three-dimensional zero mean Gaussian random variable with unit variance (the difference of the Wiener process R(t + Δt) − R(t) ~ Ijms 13 07798f3(0,Δt); for convenience the Δt has been included in the square root in Equation (3) such that readily available standard normal random numbers can be used).
As long as intramolecular potentials/forces are included, they prevent the overlap of particles [40]. Otherwise the new position will have to be rejected if the particle jumped to an excluded position [4143]. Thus intracellular structures and all the other molecules hamper the diffusion of each molecule. While the expected mean squared displacement (MSD)
( x ( t ) - x ( 0 ) ) 2 = 2 d D t
should depend linearly on D, t and the dimension d, it is found that it only grows with a reduced D′ if diffusion is hindered by other objects. Several Monte Carlo or Brownian dynamics simulations evaluated the effect of fixed crowding structures on the diffusion of tracer molecules in the cytoplasm or cellular compartments [30,42,4448]. Especially the intracellular matrix built by the cytoskeleton and proteins bound to it (microtrabecular lattice) looks like porous media and can cover up to 20% of the intracellular volume [49,50]. Diffusion through such structures leads to similar effective diffusion coefficients [42,43,51], like the predicted D/D′ from Weissberg [52] for porous media (Note that the comparison can require volume averaging [53]). In general, bigger molecules are hindered more strongly in their mobility, and at the same excluded volume fraction many small obstacles lead to a stronger hindrance compared to a few big ones [37,42,49]. Large diffusing objects with respect to the mesh size of the cytoskeleton network can be caged/trapped in these meshes and are restricted to their subvolume, which means that their MSD is limited [42,49]. Several theories have been developed to describe diffusion through polymer solutions like the cytoplasm [5458]. The complexity of the cytoplasm requires to use approximations or to rely on empirical formulas to compute the expected effective diffusion depending on the molecular radius (and the shape), the free/excluded volume fraction, and the size (distribution, shape...) of the crowding objects [49,59,60].
Other studies also revealed subdiffusion or at least transient anomalous behavior [58]. This means, that the MSD did grow only α tα with α < 1 [40,42,61], a fact that can also be observed with Fluorescence Correlation Spectroscopy (FCS) [40] The degree of subdiffusion also is a measure of the crowding level in the cell [40,62]. Such studies with mobile crowding objects require much more computation time than simulations with static objects [63]. In addition, the obstacles hinder the test molecules, while the test molecules influence the mobility of the obstacles. This feedback requires to model molecular crowding as exactly as possible in order to analyze the effective diffusion in the cell, for instance by using an obstacle size and abundance distribution similar to the in vivo conditions [41,60]. Table 1 shows relations between molecular weight and hydrodynamic radius of the molecules as used in such detailed simulations. Thus mesoscale simulations can be used to calculate the mobility of biomolecules on the cellular level, although techniques like MDS [59], dissipative particle dynamics [64] or multiparticle collision dynamics [65] have been employed as well. Connected Brownian particles can also represent polymer chains like in the Brownmove package [34]. Such simulations allow the analysis of filament shapes and diffusion [26,66] and the rheology of biopolymer networks [27,29,35,67]. Simulations with motor proteins furthermore show how network patterns emerge [68].
In addition to undirected diffusion, cells also contain motor proteins that can pull molecules along the cytoskeleton [73]. This effect can either be included in the model as a general drift term in U(xj) or by switching the transport mode if molecules bind to a cytoskeleton filament and moving the molecules linearly in their direction [74,75]. Of course also the action and properties of the motors have been studied extensively in detail using MDS [76,77]. Eventually, also the effect of the dynamic cytoskeleton and cell shape has to be included in the simulations [78].
Most biomolecules are not inert, which means that they interact with the others. This interaction is for instance required to bind to the motor proteins for the directed transport. But molecules can also unspecifically bind to the cytoskeleton [79]. This sequesters molecules from the cytosol making it less crowded, forms the microtrabecular lattice in the cell and can co-localize molecules which belong together [4]. Transient binding of course reduces the mobility of the molecules, because bound molecules are immobile [79,80].

2.2. Reactions in the Cell

Reactions between molecules require that these molecules come into contact. Since the motion of the molecules is mostly determined by diffusion, all reactions between molecules are diffusion controlled. Smoluchowski [81] calculated the collision “rate constant” of spherical molecules (here between molecule i and j with the respective radii ri, rj and diffusion coefficients Di, Dj ):
γ i j = 4 π ( r i + r j ) ( D i + D j )
and obviously no reaction can occur with a faster rate than this collision process. For example two molecules with ri = rj = 2.5 nm (corresponding to a molecular weight of about 25 kDa, cf. Table 1) and diffusion coefficient Di = Dj = 80 μm2/s have a diffusion limit for their reaction of γij = 6.05 × 109 M−1s−1. Probably enzymes cannot react that efficiently, and thus this limit is in agreement with observed biomolecular reaction rate constants.
Let kij denote the macroscopic (observable) reaction rate constant between molecule i and j, i.e., the rate constant which is also used in ODE/PDE models, assuming (locally) well mixed conditions. Obviously, the microscopic situation of two molecules being in contact is fundamentally different. The corresponding microscopic rate constant κij describes the fraction of collisions that lead to reactions. These microscopic rates can be observed e.g., in MDS or BD where the reaction dynamics is solely described by interaction potentials [8284]. In three dimensions, the macroscopic rate constant related to the microscopic one is given by [85,86]
1 k i j = 1 γ i j + 1 κ i j
Therefore, if diffusion is fast enough (1ij → 0), the microscopic rate equals the macroscopic one. The possible difference between these rates should however be considered, when comparing MDS and ODE models: while ODE models use the macroscopic descriptions, in MDS the simulation volume is so small that the molecules are basically always in contact, i.e., the microscopic description is more applicable [82].
The physiological conditions in the cell can have a strong effect on such reactions, and render them different from in vitro kinetics obtained in test tubes [87,88]. Of course all kinds of co-factors and allosteric modifiers can alter the functionality of an enzyme and the “pressure” by the surrounding crowding molecules can change the molecular conformation [87,89]. But in addition molecular crowding alters the collision frequency of the molecules [88].
Let us assume we have Ni and Nj molecules in the cell with volume V and the reaction is described by mass action kinetics. Then the (ODE) reaction rate is given based on the concentrations cm = Nm/V, where m is species i or j: rij = kijcicj (the stochastic propensity is aij = kijNiNj/V respectively). This rate/propensity changes due to excluded volume effects, because only a smaller volume is accessible for the molecules. The molecules effectively have a higher concentration, but based on the original volume, the rate constant will appear to be increased (compared to the in vitro case in diluted conditions) [80,88].
In contrast, the reduced diffusion leads to a reduced collision frequency, as can be seen from Equation (5). As long as we assume that the microscopic reaction rate constant remains constant (no allosteric change due to crowding), the macroscopic reaction rate constant appears to be decreased (cf. Equation (6)) [80]. Moreover, if diffusion becomes nonlinear due to crowding, the reaction rate constant can appear to be time dependent. Several in silico studies have observed such fractal kinetics for instance due to compartmentalization or crowding, when the obstacles hamper the restoration of the well mixed molecule distribution that is critical for spatially non-resolved population level descriptions of the dynamics [17,90,91]. Several studies therefore targeted reaction diffusion processes under crowding [37,41,63,80,92]. Further studies also investigated reactions within tubular interconnected reaction compartments such as the ER [93,94] or reactions that are enhanced by transport with motor proteins [95,96].
Finally it should be noted that the dimensionality of the system plays an important role for reactions. The collision rate constant γij has different expressions in one, two, and three dimensions [97,98]. Thus it was for instance discussed whether signaling molecules are bound to the plasma membrane to make use of the resulting 2D kinetics [5,99]. Also the membranes can be crowded and structured, and in 2D caging effects occur much faster, leading to complex reaction diffusion behavior [100].

2.3. Reactions in the Simulation: Implementation Issues

Especially for the description of biological processes several nonlinear reaction schemes and kinetics have been developed. Examples are Michaelis–Menten and Briggs–Haldane enzyme kinetics [101] but also Hill kinetics to describe cooperativity and additional effects [102]. These schemes provide one reaction rate for the whole process (under the quasi-steady state assumptions [103]), while the individual steps of the process are described by mass action kinetics [101]. The effective rate constant for Michaelis–Menten kinetics describing the process
E + S k r k f E S k c a t E + P
or in terms of balance Equations
- d [ S ] d t = d [ P ] d t = k ( [ S ] ) [ E ] [ S ]
(with [E] = const) depends on the overall current substrate concentration [S] [101]
k ( [ S ] ) = k c a t K m + [ S ]
and the two parameters kcat and Km. This means that information is lost from the original description with three reactions (Km = (kr + kcat)/kf ). With respect to particle based simulations, where reactions are triggered by collisions of molecules, this means that each colliding pair would have to know the current overall concentration [S] in order to determine the rate applicable to its own reaction, which is unphysical. Therefore the authors suggest to only use plain mass action kinetics in particle based simulations on the molecular level.
The corresponding reaction rates and rate constants for mass action kinetics count/determine the (average) number of events per time unit. These rates can also be observed in MDS or BD where the dynamics is solely described by interaction potentials [8284]. However with respect to a whole cell simulation, an event based algorithm using the macroscopic rate constants as parameters is more desirable.
The highest order of reaction which can be reasonably modeled are second order/bimolecular reactions, because third order reactions mean that three particles collide simultaneously. This is extremely unlikely and mostly such reactions form intermediary dimers which collide with the third particle instead. Accordingly third order reactions should be modeled as a set of second order reactions.
First order reactions with the rate ri = kici or propensity ai = kiNi require to execute (on average) ai events per time unit on the population level. In a discrete time particle based simulation, each particle will react in the time interval (t, tt] with probability [28]
P i = 1 - exp  ( k i Δ t ) k i Δ t
In the simulation therefore in each time step and for each particle involved in this reaction a uniform random number ξ1 has to be compared with Pi. The particle will react if ξ1Pi with ξ1 ~ Ijms 13 07798f4[0, 1] (and of course Pi < 1).
A lot of random numbers and computation time can be saved, if instead the time to the next event is pre-computed when the molecule is created in analogy to Gillespie’s stochastic simulation algorithm [20]. The waiting times ti for all molecules are exponentially distributed (Exp(ki)), and in this context it is worth noting that the minimum time min (ti) ~ Exp(kiNi). Thus this description is compatible with the chemical master equation description on the population level. The trade-off is that all individual waiting times have to be stored and ordered, executed in their sequence, and especially updated if other processes interfere with the assigned reaction channel [104106].
Zero order reactions ∅ → I with rate ri = ki(0) or propensity ai = ki can be implemented similarly. New molecules of the product species I are created (out of nothing) after the waiting time min (ti) ~ Exp(ki), where the position of the new molecule is drawn from a given spatial distribution. Alternatively they can be generated by a first order reaction based on a dummy species D with fixed concentration cD and the scheme DD + I and rate constant k′ = ki(0)/cD [75]. Since the number of molecules is changing, a buffer of empty particles is required for the simulation [107].
Second order reactions require a more elaborate description [39]. As stated above, bimolecular reactions are triggered by (diffusion-controlled) collisions of molecules, i.e., the simulation has to identify the pairs currently in contact. In a discrete time simulation the number of pairs which will be found in a time interval thus depends on Δt, i.e., how frequently it is checked. The true number of collisions remains unknown, because the path of the molecules between t and tt is undetermined. Clifford and Green [108] suggested a Brownian bridge to interpolate the paths of the molecules in order to find all collisions. Lapin et al. [109] in turn used the Fokker–Planck equation in order to determine the probability of a reaction within Δt given the current distance ||xi(t) − xj(t)||, microscopic rate constant κij and diffusion coefficient Dij = Di + Dj (see Figure 2). Note that in this framework the new particle position is not updated according to Equation (3) but according to the probability density distribution (pdf) such that the particles will never overlap. Obviously this description has to be applied only if the particles are within the interaction range. For separated particles, the pdf converges to the normal distribution.
Again, most steps will not result in reactions. Similarly to the treatment of first order reactions using Gillespie’s approach also here the time to the event can be computed instead of checking and wasting random numbers each step. The reaction-diffusion equation for the probability density function (see Figure 2) can be solved analytically using Green’s functions from which the probability of a reaction between particle i at xi(t) and j at xj(t) and if so also the time t + τ* and position x* of it can be deduced [110]. The resulting Green’s function reaction dynamics (GFRD) is in agreement with the corresponding analytical solutions for diffusion limited reactions [111113].
The challenge in the underlying theory is that it can only be solved analytically for pairs, but in the cell we have much more molecules [115]. Therefore the step length or event horizon respectively in GFRD has to be set such that no additional molecule enters the vicinity of a pair [110]. For low particle concentrations, in turn, GFRD allows a great speedup of the simulation.
Several methods aim at a simpler model based on a critical radius rij* (kij, Δt) such that molecules will react if their distance is smaller than r*ij [28,80,116118]. Such an algorithm will compute reactions faster than GFRD or the Fokker–Planck method because it does not need complex look-up tables. However, in order to reach the same accuracy, a shorter time step than in GFRD might have to be chosen, which requires more steps in total to complete the simulation [111]. In addition, the computation of a suitable critical radius r*ij(kij,Δt, . . . ) is not straightforward, requiring numerical calculations before the simulation can start and can be implementation dependent [28,39,118].
A stable and rather simple derivation is given in [80] for particles that can overlap:
  • set the critical reaction radius to the physical collision radius
    r i j * = r i + r j
  • and execute reactions for particles with ||xi(t) − xj(t)|| ≤ r*ij with probability
    P i j * = κ i j Δ t 4 π ( r i + r j ) 3 / 3
Obviously P* ij < 1, which limits Δt. From tests we found that this approach works reliably up to P*ij < 0.2 even for significant degrees of diffusion control. Theoretical limits of this approach or respective correction factors for the reaction probability are calculated in [39].
It is especially important for this approach to work correctly, that the critical radius is the same as the radius used to determine the diffusion limit γij Equation (5) in order to obtain the expected macroscopic rate constant. Note, that bimolecular rate constants (for reactions in 3D volumes) have units volume/time (conversion from M−1s−1 to μm3/s using 1015/6.022 × 1023 Mμm3). This can be interpreted as reaction volume per molecule and time [116]. The present approach uses the ratio of the microscopic reaction volume to the interaction volume as reaction probability given that a collision has happened, which is a mechanistic analogy to diffusion controlled reaction scheme where the microscopic rate constant describes how efficiently collisions are turned into reactions. Based on this formalism also reactions with other geometries can be defined, for instance binding/adsorption to membranes or cytoskeleton filaments [75]. If the reacting objects are not allowed to overlap, obviously the reaction volume has to be wrapped around the collision radius rather than being distributed within their collision radius. If the thickness of this reaction volume layer compared to the collision radius is negligible, the constraint that the critical reaction radius should match the collision radius in Equation (5) will be satisfied without using a reaction probability [75].
Reversible reactions like A + BC require a special treatment in this context [28,111,119,120]. Not only the association but also dissociation can be diffusion limited [121]. For strong association and slow diffusion, a pair of just dissociated particles can hardly diffuse away from each other before they will recombine. In order to obey detailed balance therefore also a microscopic dissociation reaction has to be introduced such that the dissociation constant
K d = [ A ] [ B ] [ C ] = k c k a b = κ c κ a b
is preserved [111]. Similar to the microscopic binding reaction rate constant, which was introduced above, the microscopic rate constant κc counts all dissociation events, while the macroscopic rate constant kc only includes the successful events where the molecules could escape from each other and reach a macroscopically observable distance. From Equation (10) follows that the backward first order reaction rate constant has the same scaling like the forward second order rate constant, which is determined by Equation (6) in 3D: κc = kc × γab/(γab − kab).
In order to reach the correct geminate recombination rate also the initial condition of the new pair A and B upon dissociation has to be chosen carefully. In the Smoldyn framework of Andrews et al. [28,119,120] as well as in GFRD [110,111] therefore an unbinding radius outside of the critical binding/reaction radius has to be used. Since space and time are related, the onset of the geminate recombination process can likewise be delayed by blocking reactions for particles that just reacted. If the particles A and B are placed at identical positions xA = xB (= xC from which they were created from) this delay term is
τ i = 1 10 ( r a + r b ) 2 / ( D a + D b ) , i = a , b
which is the maximum likelihood of the (diffusion-controlled) escape process [122]. Placing the molecules overlapping can be necessary in simulations with crowding objects, because the dissociation CA + B might occur at a position where no valid separated positions for A and B exists. Note that such crowding structures which reduce diffusion keep the molecules together and therefore push equilibrium towards the associated C state.
Depending on the implementation for parallelized execution in multi-threaded CPU or GPU applications anyway a flag might be necessary to execute each reaction for each molecule not more than once [107]. This flag can be implemented using the same memory like the necessary functional delay term Equation (11) by preventing that molecules react in any reaction, which have a delay > t assigned, and setting the minimum delay for each particle that reacted to t + 0.5Δt. Note that the global delay for reversible reactions does not affect the overall reaction process because it is so small [122].
Such a delay term can also be used in order to mimic Michaelis–Menten enzyme kinetics[75,107]. In the reduced reaction scheme E +SE +P the enzyme E is blocked for a certain time τ after each reaction, leading to a saturation of the reaction rate. Let us use the standard notation kcat and kr for the first order reactions and kES instead of kf for the second order E +SES reaction [101]. If only kcat and Km are given, kES can be set to kES = kcat/Km [75,107].
The delay term introduced here for reversible reactions is also in agreement with the overdamped Langevin regime discussed here. Upon dissociation we would expect that the two molecules have to pick up and maintain impulses which separate them for a given time. The reversion of these impulses such that the molecules can re-associate certainly requires additional time, which is not covered by the plain random walk implementation. In the detailed dissociation process, the molecule will first have to accumulate the necessary energy in order to overcome the potential energy of the bond and to dissociate. Once it has this energy, the fate of the molecule is already dissociation, but only with a delay the bond will be actually broken. A similar discussion, however in terms of the adjoint dissociation distance, is found in Andrews et al. [120].
Future work has to find ways to include more detailed description on such a simplified simulation scheme or at least bridge the detailed and simplified approach. Especially nonspherical molecules and specific binding sites on the surface of the molecules are relevant for biomolecular reaction kinetics [123128].

2.4. Performance and Accuracy

The performance of a particle based simulation strongly depends on the efficiency of the used pair finding and collision detection algorithms [129]. In order to cover a sufficient time span with the simulation, also a large Δt is desirable to reduce the number of steps. Constraints on the accuracy limit the choice of Δt.
Within a crowded intracellular environment the random walk steps have to comply with the boundaries of the obstacles. Steps can either be rejected (completely or retried) if they end in an obstacle [42,43,75,80] or reflected [110,130]. However, if complex surfaces are considered, calculation of the reflected endpoint can become time consuming and suffer from numerical imprecisions. Especially for densely crowded regions, one step might interfere with more than one obstacle. Thus obstacle density and size give upper bounds for the step length, similar to the condition in GFRD that a particle must interact with at maximum one particle within the next time step [110]. As such, the step size could be adjusted individually for each particle depending on its distance to the closest object.
The random number distribution for diffusion simulations (Equation (3)) should be a normal distribution with μ = 0, σ2 = 1. These random numbers can be obtained by several algorithms based on uniformly distributed pseudo random numbers [131]. From our observations, a normally distributed random number costs at least 20% more than a uniformly distributed one. However, the repeated application of the uniform distribution will lead to a normal distribution due to the central limit theorem [132], and the convergence sufficient for simulation occurs within 5 steps. The authors also are very faithful that a repeated rejection sampling of steps from the uniform distribution will converge to a reasonable distribution within a crowded intracellular environment, even without using complex reflection calculations [130].
In addition, the uniform distribution does not have long tails. For a σ2 = 1 uniform random variable the distribution has to have a width of 12. Taking into account the scaling factor of Equation (3), the maximum step length Δ x m a x , u = 1 / 2 12 × 2 D Δ t u = 6 D Δ t u. The maximum step length in space becomes important for instance if only the end point of a step is checked with all obstacles/boundaries. In order to not just jump across boundaries, Δxmax has to be smaller than the smallest object (similar to microscopes, where the wavelength determines the spatial resolution). The normal distribution, in contrast, should not be truncated to less than 3σ, which means that Δ x m a x , n = 3 2 D Δ t n = 18 D Δ t n. If the step length decision is determined by the spatial aspects as indicated above, then a normal distribution requires Δtn ≤ 1/tu (depending on the truncation). Taking also the higher costs of normally distributed random numbers into account, a simulation with uniformly distributed random numbers could run about 4 times faster along the simulated time.
Depending on the implementation of the reaction scheme(s), further constraints on Δt can occur. For instance the reaction probability from Equation (9) has to be smaller than 1 (or better 0.2), which leads to
Δ t r min ( 4 π ( r i + r j ) 2 / 3 κ i j )
Except for strongly diffusion controlled reactions Δtr can be much bigger than the constraints from the random walk steps Δtn > Δtu. In order to save computation time on the expensive pair finding for bimolecular reactions, reactions can thus be executed with a lower frequency, while diffusion steps executed with uniformly distributed random walk steps at higher frequency can converge to the normal distribution and correctly sample the details of the spatial structures in the cell. Likewise this approach allows refilling of reaction volumes between two reaction steps to the local concentration, which is necessary to obtain the correct reaction rate [75].

3. Applications and Results of Spatial Simulations

3.1. Current Spatial Simulation Frameworks for the Cellular Level

Table 2 lists current spatial simulation packages. The methods employed are mostly particle based (Brownian/Smoluchowski or Green’s function reaction dynamics), or lattice based using ODE/PDE/SDE methods. Due to its relevance, we will describe Gillespie’s stochastic simulation algorithm for the chemical master equation (CME) or the spatial version thereof (reaction diffusion master equation, RDME) in the next section.
Further methods that have been used for cell simulations are for instance: dissipative particle dynamics [151], an Ising model to calculate the ability of gradient sensing in the presence of multiple competitive ligands [152], cellular Potts models [153], or lattice based methods where the lattice sites have the size of molecules [105]. Special multi-scale methods have been developed to accommodate the heterogeneous process of the cells like crowding with small and large molecules [154] or slow and fast reactions [155], which also involves for instance local mean field closures for interactions [156] and coupling of several simulation methods [14].
Parallelization for all these methods has been investigated, in order to increase the performance [107,133,157,158]. Note that the computational particle object does not necessarily have to coincide with an atom, molecule, or other cellular object. It can also represent a subvolume of the cell or the membrane, exchanging mass or other content with neighbours, interact and evolute as described by additional rules [150]. Thus particle or agent based simulations can be used for many simulation tasks in biology [30].

3.2. The Reaction Diffusion Master Equation and Gillespie’s Algorithm

The chemical master equation CME or reaction diffusion master equation RDME if space is considered describes the changes in a reaction system consisting of M molecular species and K reaction channels with rate constants k = (k1, . . ., kK)T [20]. The reactions take place in reaction compartment Ω. For a spatially resolved description Ω can be subdivided into U subvolumes of volume V1, . . ., VU [144,157,159161]. Here we denote the number of particles in subvolume ν with Nν(t) = (N1ν (t), . . ., NMν (t))T. Time evolution on this level is driven by Markovian population dynamics. More specifically, the probability distribution P(N1(t) = n1, . . ., NU(t) = nU) = p(n1, . . ., nU, t) satisfies the RDME
t p ( n 1 , , n U , t ) = ( R + D ) p ( n 1 , , n U , t )
where and Ijms 13 07798f5 are the reaction and diffusion operators, respectively. The definition of the reaction operator follows from the classical master equation and thus reads
R = ν = 1 U j = 1 K ( E ν - δ j - 1 ) a j ( n ν )
where δj denotes the stoichiometric change vector associated with the j-th reaction. This change is performed by the shift operator Eν that applies to everything to its right-side. It is defined as Eνδj f(n1, . . ., nU, t) = f(n1, . . ., nνδj, . . ., nU, t) for any function f of appropriate domain—in particular f(n1, . . ., nU, t) = aj(nν)p(n1, . . ., nU, t). Diffusion of species i with diffusion coefficient Di into another volume μ is translated into a first-order transport reactions with effective propensity kiνμniν. The corresponding rate constant can be expressed as
k i ν μ = D i l 2 = D i × l 2 l × l 3 = D i S ν , μ l V ν
where Sν,μ is the surface/interface area of the cubic subvolumes. Thus, the diffusion operator Ijms 13 07798f5 reads
D = ν = 1 U μ N ( ν ) i = 1 M ( E ν Δ i E μ - Δ i - 1 ) k i ν μ n i ν
where Δi is the vector of change caused by the transport of a single molecule of type i from one volume to another one. The set Ijms 13 07798f3(ν) denotes the volumes μ in the neighborhood of volume ν.
In terms of simulating the reactions the waiting time τjν for every reaction j within volume ν is distributed exponentially with parameter aj(nν), called the propensity of reaction j within ν. The waiting time τν for any reaction to occur in ν is exponentially distributed according to parameter a 0 ( n ν ) = j = 1 K a j ( n ν ). Starting from a given time t, the next event of reaction j in ν is according to Gillespie’s algorithm [20] at
t j ν = t + τ j ν ; τ j ν ~ Exp ( a j )
A transport reaction into μ at time t′ causes a state change Nμ(t) → Nμ(t′), which according to Gillespie’s algorithm would require updating the precomputed waiting times tμ in μ. Anderson [162] proved that the remaining fraction of the time to the next reaction can simply be stretched according to the changed propensity
t + μ = t + ( t - μ - t ) × a 0 ( N μ ( t ) ) a 0 ( N μ ( t ) )
Of course there exist other ways in implementing the simulation, because the next reaction time Equation (14) can either be calculated cumulatively for all reaction based on a0 or individually for ajν for each subvolume and reaction channel. Note that the minimum waiting time of all individual waiting times will be distributed as Exp(a0) such that both descriptions are equivalent, and actually any partitioning/grouping of reaction channels is possible in order to improve the execution performance [155]. Individual reactions have to be executed in their order in time, while the resulting changes in Nν(t) can require updates in other waiting times as in Equation (15). If the waiting time is calculated based on the cumulative a0 instead, the reaction channel and compartment of the next reaction have to be found based on their individual probabilities ajν/a0. Simulators based on Gillespie’s algorithm for the RDME are for example MesoRD [143] or STEPS [146] (cf. Table 2). The discretization into the subvolumes should not be too small such that the chance of two reactants being in the same subvolume goes to zero. In that case simulations of the RDME become diffusion limited [163165].
Note that for a conversion from particle based simulations to population level of the ith species in ν all particles of that species with x(t) ∈ Vν have to be counted. Conversely, the underlying assumption of the population dynamics model is that the Niν (t) molecules are uniformly distributed in Vν.

3.3. Rule-based Modeling

The different phosphorylation levels as well as the localization of signaling molecules creates a multitude of states in a species/population based description of the system. The number of reactions and rate constants that have to be determined increases even more dramatic, because two interacting proteins each with n possible states can have up to n2 different reaction rate constants. In order to be able to generate and analyze detailed models of signal transduction despite this combinatorial explosion, rule based modeling strategies have been developed [166,167]. The syntax/input scheme of such models minimizes/structures the information that has to be entered by using state dependent rules. Based on these rules the models can also be further abstracted, decomposed or reduced for a comprehensive analysis of the complex system [168,169].
The rule-based concept can be easily extended for spatial models, where the particle interactions are now determined by their internal state and the rules assigned [170,171]. Spatial rule-based models can be for instance used to track how simple molecular units self-assemble into large structures [172].

3.4. Applications and Results

The following list shows the importance of the mesoscale level in the cell. These applications cover time and length scales above molecular dynamics simulations, though not necessary on the cellular level. Thus they can be used to analyze important questions in molecular and cell biology, ranging from more detailed descriptions of reaction process to the spatial organization of the cell.
  • Binding kinetics and binding sites: depending on the description level, protein-protein association can become quite complex [36]. For instance if multiple binding sites and diffusion-controlled reactions are considered. Biomolecules can have several binding sites for the same ligand, for instance receptors forming multimers or antibodies [128]. Kang et al. [173] analysed this and Park et al. [174] developed a theory for reversible reactions under these circumstances. For instance, two binding sites on a molecule would mean that the microscopic reaction rate constant κij is doubled, while the reaction radius is the same as for a molecule with just one binding site. Equation (6) shows that the macroscopic rate constant will not necessarily double under these circumstances.
  • Scaffolds and Channeling: Both in signaling and metabolic pathways co-localization of related molecules has been observed. Obviously co-localization has advantages because the local high concentration boosts the reaction rate [4,80,106,175,176]. Specific and even nonspecific binding interactions which modify the localization properties of molecules can thus enhance reactions [177]. Note, that the localization requires that molecules do not diffuse around/away, such that there is a trade-off between advantages due to co-localization and disadvantages due to the reduced mobility [75,80].
  • Protein DNA interactions: Transcription factors have to find their target site on the DNA amongst millions of binding sites, and they do it surprisingly efficiently, e.g., by combining 1D sliding and 3D diffusion [178]. For instance nonspecific interactions could enhance association rates respectively [177]. Note that even DNA is well organized in space [6]. The spatial organization of DNA strands plays an important role, but long DNA strands can obviously not be modeled with full atomic detail in a MDS simulation such that multi-scale approaches have to be employed [25]. The observed bursting kinetics of transcription rates is likewise explained using open and closed chromatin states, which involve large-scale transitions of the DNA state [179].
  • Assembly and fusion processes: Large polymer structures such as the cytoskeleton filaments play an important role for the spatial organization of the cell. Guo et al. modeled the actin assembly using Brownian dynamics [180]. Langevin dynamics have been used to simulate the assembly of virus polymers [181]. A rule-based description was likewise used to analyze the emergence of complex structures [172]. Likewise the fusion of membrane enclosed structures like vesicles is important for the functionality of the cell [64,151,182]. Note that the interplay of cytoskeleton filaments, motor proteins and vesicles can enhance their fusion process [150], while the cytoskeleton structure is for instance organized by the aforementioned growth processes but also motors pulling them together and creating spatial patterns [68].
  • Non-uniform molecule distributions in space: In order to grow/move in specific directions cells have to polarize into front and back, which is associated with nonsymmetric particle distributions across the cell [134,183]. In addition receptors on the membrane can cluster together [9,184], and the output of spatial simulations shows the importance of the spatial organization in the cell [9]. Note that again reversible binding and/or unspecific binding interactions influence these reaction rates [177,184].

4. Towards Multi-Scale Simulations from Atoms to Cells

Simulations on the molecular and cellular level are available and have been used extensively to better understand biological processes [1519,39]. At least theoretically the different levels can be included in a multi-scale simulation, where all levels are executed in parallel. Starting from a coarse grained Brownian dynamics simulation, more detailed levels could be invoked as soon as two particles are closer than a threshold, e.g., switching from spheres to more realistic shapes [60,185] and eventually resolving each collision on the MDS level.
A sequential execution would probably be easier to implement and thus more relevant in order to model a signaling pathway “ab initio”, starting from the protein sequences of all involved molecules. The corresponding molecular structures can be predicted (at least for small signaling molecules) if not yet known from structure analysis by molecular simulation/analysis tools based on the atomic level [22,186189].
If possible from these structures the hydrodynamic radii [60,190], binding rates and reaction kinetics can be calculated/estimated [82,188,189,191198]. Especially important with respect to signal transduction are the conformational changes and their rates associated with “active” and “inactive” states, mostly regulated by the phosphorylation level of the proteins [199]. Based on the diffusion model and the crowding level of the cell, the effective diffusion coefficient can then be estimated (e.g., if crowding should only be modeled implicitly) [12]. The number/concentration of signaling molecules can nowadays be quantified [200] and likewise serves as an input for the simulation [12].
Thus all necessary parameters are available to simulate the spatial and temporal dynamics of a signaling pathway in the cell. The relation of transport and reaction rates beforehand gives an estimator how well mixed the system will be and if the spatial dynamics have to be resolved explicitly [7]. If not, even a population based ODE/SDE model with the effective reaction rates might be sufficient for simulation (of the respective well-mixed subset of the system). But of course the compartmentalization and thus transport rates through nuclear pores/membranes should be included [201].
Given the challenges involved with parameter estimation and model discrimination from experimental results [179,202205], such a bottom-up simulation approach might bring invaluable insights for the signaling pathway under consideration. Simulations tracking all particle positions inside the cell allow analysis of particle distributions, e.g., to find receptor clusters in the membrane [9,184,206]. Based on such simulations the cell can also be visualized interactively at the desired resolution [11,74,107,129,207], for instance showing local variations in molecule concentrations [207] or even showing the molecules with their atomic structure as in Figure 1. The described bottom-up simulation approach as such cannot include all relevant effects on all levels. Leaving out solvent molecules is the most obvious fact, but intermediary effects like crowding and hydrodynamic forces [12,37] seem to be important. Thus especially efficient yet correct multi-scale simulation algorithms for the mesoscale resolution have to be further developed, and a consensus on the necessary levels of details has to be found.

Acknowledgments

The authors thank Peter R. ten Wolde, Alexei Lapin, Martin Falk and Pablo de Heras Ciechomski for valuable discussions and ScienceVisuals for rendering our simulated particle distributions. We acknowledge funding through the SystemsX.ch initiative and the Swiss Confederation’s Commission for Technology and Innovation (CTI) project 12532.1 PFLS-LS. H.K. acknowledges the support from the Swiss National Science Foundation, grant no. PP00P2 128503.

References

  1. Tomita, M. Whole-cell simulation: A grand challenge of the 21st century. Trends Biotechnol 2001, 19, 205–210. [Google Scholar]
  2. Mathews, C.K. The cell-bag of enzymes or network of channels? J. Bacteriol 1993, 175, 6377–6381. [Google Scholar]
  3. Alberts, B.; Johnson, A.; Lewis, J.; Raff, M.; Roberts, K.; Walter, P. Molecular Biology of the Cell; Garland Science: New York, NY, USA, 2002. [Google Scholar]
  4. Ovádi, J.; Srere, P.A. Macromolecular compartmentation and channeling. Int. Rev. Cyt 1999, 192, 255–280. [Google Scholar]
  5. Bray, D. Signaling complexes: Biophysical constraints on intracellular communication. Ann. Rev. Biophys. Biomol. Struct 1998, 27, 59–75. [Google Scholar]
  6. Lieberman-Aiden, E.; van Berkum, N.L.; Williams, L.; Imakaev, M.; Ragoczy, T.; Telling, A.; Amit, I.; Lajoie, B.R.; Sabo, P.J.; Dorschner, M.O.; et al. Comprehensive mapping of long-range interactions reveals folding principles of the human genome. Science 2009, 326, 289–293. [Google Scholar]
  7. Kholodenko, B.N. Cell-signalling dynamics in time and space. Nat. Rev. Mol. Cell Biol 2006, 7, 165–176. [Google Scholar]
  8. Kholodenko, B.N.; Hancock, J.B.; Kolch, W. Signalling ballet in space and time. Nat. Rev. Mol. Cell Biol 2010, 11, 414–426. [Google Scholar]
  9. Costa, M.N.; Radhakrishnan, K.; Wilson, B.S.; Vlachos, D.G.; Edwards, J.S. Coupled stochastic spatial and non-spatial simulations of ErbB1 signaling pathways demonstrate the importance of spatial organization in signal transduction. PLoS One 2009, 4. [Google Scholar] [CrossRef]
  10. ScienceVisuals. Available online: http://www.sciencevisuals.com accessed on 08 June 2012.
  11. De Heras Ciechomski, P.; Mange, R.; Peternier, A. Two-Phased Real-Time Rendering of Large Neuron Databases. Proceedings of the 2008 International Conference on Innovations in Information Technology, Al Ain, United Arab Emirates, 16–18 December 2008; pp. 712–716.
  12. Ando, T.; Skolnick, J. Crowding and hydrodynamic interactions likely dominate in vivo macromolecular motion. Proc. Natl. Acad. Sci. USA 2010, 107, 18457–18462. [Google Scholar]
  13. Rafelski, S.M.; Marshall, W.F. Building the cell: Design principles of cellular architecture. Nat. Rev. Mol. Cell Biol 2008, 9, 593–602. [Google Scholar]
  14. Bittig, A.T.; Uhrmacher, A.M. Spatial Modeling in Cell Biology at Multiple Levels. Proceedings of the Winter Simulation Conference (WSC), Baltimore MD, USA, 5–8 December 2010; pp. 608–619.
  15. Takahashi, K.; Arjunan, S.N.V.; Tomita, M. Space in systems biology of signaling pathways— Towards intracellular molecular crowding in silico. FEBS Lett 2005, 579, 1783–1788. [Google Scholar]
  16. Tolle, D.P.; le Novere, N. Particle-based stochastic simulation in systems biology. Curr. Bioinf 2006, 1, 315–320. [Google Scholar]
  17. Turner, T.E.; Schnell, S.; Burrage, K. Stochastic approaches for modelling in vivo reactions. Comput. Biol. Chem 2004, 28, 165–178. [Google Scholar]
  18. Ridgway, D.; Broderick, G.; Ellison, M.J. Accommodating space, time and randomness in network simulation. Curr. Opin. Biotechnol 2006, 17, 493–498. [Google Scholar]
  19. Burrage, K.; Burrage, P.; Leier, A.; Marquez-Lago, T.; Nicolau, D., Jr. Stochastic simulation for spatial modelling of dynamic process in a living cell. Des. Anal. Biomol. Circuits: Eng. Approaches Syst. Synth. Biol 2011, 43–62. [Google Scholar]
  20. Gillespie, D.T. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comp. Phys 1976, 22, 403–434. [Google Scholar]
  21. Leach, A.R. Molecular Modelling: Principles and Applications; Pearson Education Ltd: Harlow, England, 2001. [Google Scholar]
  22. Schlick, T. Molecular Modeling and Simulation: An Interdisciplinary Guide; Springer Verlag: Berlin, Germany, 2010. [Google Scholar]
  23. Bionumbers. Available online: http://bionumbers.hms.harvard.edu/bionumber.aspx?&id=106198&ver=2 accessed on 8 June 2012.
  24. Tozzini, V. Coarse-grained models for proteins. Curr. Opin. Struct. Biol 2005, 15, 144–150. [Google Scholar]
  25. Villa, E.; Balaeff, A.; Mahadevan, L.; Schulten, K. Multiscale method for simulating protein-DNA complexes. Multiscale Model. Simul 2004, 2, 527–553. [Google Scholar]
  26. Chandran, P.L.; Mofrad, M.R.K. Averaged implicit hydrodynamic model of semiflexible filaments. Phys. Rev. E 2010, 81, 031920:1–031920:17. [Google Scholar]
  27. Cyron, C.J.; Wall, W.A. Numerical method for the simulation of the Brownian dynamics of rod-like microstructures with three-dimensional nonlinear beam elements. Int. J. Numer. Methods Eng 2012, 90, 955–987. [Google Scholar]
  28. Andrews, S.S.; Bray, D. Stochastic simulation of chemical reactions with spatial resolution and single molecule detail. Phys. Biol 2004, 1, 137–151. [Google Scholar]
  29. Sandersius, S.A.; Newman, T.J. Modeling cell rheology with the Subcellular Element Model. Phys. Biol 2008, 5. [Google Scholar] [CrossRef]
  30. Sbalzarini, I.F.; Walther, J.H.; Bergdorf, M.; Hieber, S.E.; Kotsalis, E.M.; Koumoutsakos, P. PPM–A highly efficient parallel particle–mesh library for the simulation of continuum systems. J. Comp. Phys 2006, 215, 566–588. [Google Scholar]
  31. Newman, T.J. Grid-free models of multicellular systems, with an application to large-scale vortices accompanying primitive streak formation. Curr. Topics Dev. Biol 2008, 81, 157–182. [Google Scholar]
  32. Holcombe, M.; Adra, S.; Bicak, M.; Chin, S.; Coakley, S.; Graham, A.; Green, J.; Greenough, C.; Jackson, D.; Kiran, M.; et al. Modelling complex biological systems using an agent-based approach. Integr. Biol 2012, 4, 53–64. [Google Scholar]
  33. Walker, D.C.; Southgate, J. The virtual cell a candidate co-ordinator for middle-outmodelling of biological systems. Brief. Bioinf 2009, 10, 450–461. [Google Scholar]
  34. Geyer, T. Many-particle Brownian and Langevin Dynamics Simulations with the Brownmove package. BMC Biophys 2011, 4. [Google Scholar] [CrossRef]
  35. Kim, T.; Hwang, W.; Lee, H.; Kamm, R.D. Computational analysis of viscoelastic properties of crosslinked actin networks. PLoS Comput. Biol 2009, 5. [Google Scholar] [CrossRef]
  36. Gabdoulline, R.R.; Wade, R.C. Protein-protein association: Investigation of factors influencing association rates by Brownian dynamics simulations. J. Mol. Biol 2001, 306, 1139–1155. [Google Scholar]
  37. Sun, J.; Weinstein, H. Toward realistic modeling of dynamic processes in cell signaling: Quantification of macromolecular crowding effects. J. Chem. Phys 2007, 127, 155105:1–155105:10. [Google Scholar]
  38. Schmidt, R.R.; Cifre, J.G.H.; de la Torre, J.G. Comparison of Brownian dynamics algorithms with hydrodynamic interaction. J. Chem. Phys 2011, 135, 084116:1–084116:10. [Google Scholar]
  39. Erban, R.; Chapman, S.J. Stochastic modelling of reaction–diffusion processes: Algorithms for bimolecular reactions. Phys. Biol 2009, 6, 046001. [Google Scholar]
  40. Weiss, M.; Elsner, M.; Kartberg, F.; Nilsson, T. Anomalous subdiffusion is a measure for cytoplasmic crowding in living cells. Biophys. J 2004, 87, 3518–3524. [Google Scholar]
  41. Ridgway, D.; Broderick, G.; Lopez-Campistrous, A.; Ru’aini, M.; Winter, P.; Hamilton, M.; Boulanger, P.; Kovalenko, A.; Ellison, M.J. Coarse-grained molecular simulation of diffusion and reaction kinetics in a crowded virtual cytoplasm. Biophys. J 2008, 94, 3748–3759. [Google Scholar]
  42. Klann, M.T.; Lapin, A.; Reuss, M. Stochastic simulation of signal transduction: Impact of the cellular architecture on diffusion. Biophys. J 2009, 96, 5122–5129. [Google Scholar]
  43. Trinh, S.; Arce, P.; Locke, B.R. Effective diffusivities of point-like molecules in isotropic porous media by monte carlo simulation. Trans. Porous Media 2000, 38, 241–259. [Google Scholar]
  44. Długosz, M.; Trylska, J. Diffusion in crowded biological environments: Applications of Brownian dynamics. BMC Biophys 2011, 4. [Google Scholar] [CrossRef]
  45. Chang, R.; Jagannathan, K.; Yethiraj, A. Diffusion of hard sphere fluids in disordered media: A molecular dynamics simulation study. Phys. Rev. E 2004, 69. [Google Scholar] [CrossRef]
  46. Ölveczky, B.P.; Verkman, A.S. Monte carlo analysis of obstructed diffusion in three dimensions: Application to molecular diffusion in organelles. Biophys. J 1998, 74, 2722–2730. [Google Scholar]
  47. Verkman, A.S. Solute and macromolecule diffusion in cellular aqueous compartments. Trends Biochem. Sci 2002, 27, 27–33. [Google Scholar]
  48. Lipkow, K.; Andrews, S.S.; Bray, D. Simulated diffusion of phosphorylated CheY through the cytoplasm of Escherichia coli. J. Bacteriol 2005, 187, 45–53. [Google Scholar]
  49. Luby-Phelps, K. Cytoarchitecture and physical properties of cytoplasm: Volume, viscosity, diffusion, intracellular surface area. Int. Rev. Cytol 2000, 192, 189–221. [Google Scholar]
  50. Jacobson, K.; Wojcieszyn, J. The translational mobility of substances within the cytoplasmic matrix. Proc. Natl. Acad. Sci. USA 1984, 81, 6747–6751. [Google Scholar]
  51. Blum, J.J.; Lawler, G.; Reed, M.; Shin, I. Effect of cytoskeletal geometry on intracellular diffusion. Biophys. J 1989, 56, 995–1005. [Google Scholar]
  52. Weissberg, H.L. Effective diffusion coefficient in porous media. J. Appl. Phys 1963, 34, 2636– 2639. [Google Scholar]
  53. Whitaker, S. The Method of Volume Averaging; Springer: Berlin, Germany, 1998. [Google Scholar]
  54. Fan, T.H.; Dhont, J.K.G.; Tuinier, R. Motion of a sphere through a polymer solution. Phys. Rev. E 2007, 75. [Google Scholar] [CrossRef]
  55. Ogston, A.G.; Preston, B.N.; Wells, J.D. On the transport of compact particles through solutions of chain-polymers. Proc. R. Soc. Lond. Ser. A Math. Phys. Sci 1973, 333, 297–316. [Google Scholar]
  56. Cukier, R.I. Diffusion of Brownian spheres in semidilute polymer solutions. Macromolecules 1984, 17, 252–255. [Google Scholar]
  57. Han, J.; Herzfeld, J. Macromolecular diffusion in crowded solutions. Biophys. J 1993, 65, 1155–1161. [Google Scholar]
  58. Bruna, M.; Chapman, S.J. Excluded-volume effects in the diffusion of hard spheres. Phys. Rev. E 2012, 85. [Google Scholar] [CrossRef]
  59. Sakha, F.; Fazli, H. Three-dimensional Brownian diffusion of rod-like macromolecules in the presence of randomly distributed spherical obstacles: Molecular dynamics simulation. J. Chem. Phys 2010, 133, 234904:1–234904:6. [Google Scholar]
  60. Ando, T.; Skolnick, J. Brownian Dynamics Simulation of Macromolecule Diffusion in a Protocell. Proceedings of the International Conference of the Quantum Bio-Informatics IV, Tokyo, Japan, 10–13 March 2010; Georgia Institute of Technology and World Scientific Publishing: Atlanta, GA, USA, 2011; 28, pp. 413–426. [Google Scholar]
  61. Saxton, M.J. A biological interpretation of transient anomalous subdiffusion. I. Qualitative model. Biophys. J 2007, 92, 1178–1191. [Google Scholar]
  62. Hiroi, N.; Lu, J.; Iba, K.; Tabira, A.; Yamashita, S.; Okada, Y.; Flamm, C.; Oka, K.; Köhler, G.; Funahashi, A. Physiological environment induces quick response–slow exhaustion reactions. Frontiers Physiol 2011, 2. [Google Scholar] [CrossRef]
  63. Echeveria, C.; Tucci, K.; Kapral, R. Diffusion and reaction in crowded environments. J. Phys. Condens. Matter 2007, 19. [Google Scholar] [CrossRef]
  64. Shillcock, J. Insight or illusion? Seeing inside the cell with mesoscopic simulations. HFSP J 2008, 2, 1–6. [Google Scholar]
  65. Kapral, R. Multiparticle Collision Dynamics: Simulation of Complex Systems on Mesoscales. In Advances in Chemical Physics; Rice, S.A., Ed.; John Wiley and Sons, Inc: Hoboken, NJ, USA, 2008; Volume 140, pp. 89–146. [Google Scholar]
  66. Cyron, C.J.; Wall, W.A. Consistent finite-element approach to Brownian polymer dynamics with anisotropic friction. Phys. Rev. E 2010, 82, 66705:1–66705:12. [Google Scholar]
  67. Lee, H.; Pelz, B.; Ferrer, J.M.; Kim, T.; Lang, M.J.; Kamm, R.D. Cytoskeletal deformation at high strains and the role of cross-link unfolding or unbinding. Cell. Mol. Bioeng 2009, 2, 28–38. [Google Scholar]
  68. Karsenti, E.; Nédélec, F.; Surrey, T. Modelling microtubule patterns. Nat. Cell Biol 2006, 8, 1204–1211. [Google Scholar]
  69. Renkin, E.M. Multiple pathways of capillary permeability. Circ. Res 1977, 41, 735–743. [Google Scholar]
  70. Taylor, A.E.; Granger, D.N. Exchange of Macromolecules across the Microcirculation. In Handbook of Physiology: The Cardiovascular System: Microcirculation; American Physiological Society: Bethesda, MD, USA, 1984; Volume 4, pp. 467–520. [Google Scholar]
  71. Zimmerman, S.B.; Trach, S.O. Estimation of macromolecule concentrations and excluded volume effects for the cytoplasm of Escherichia coli. J. Mol. Biol 1991, 222, 599–620. [Google Scholar]
  72. Niederalt, C. Bayer Technology Services. PK-Sim/MoBi from Bayer Technology Services. Personal communication, 2011. [Google Scholar]
  73. Vale, R.D. The molecular motor toolbox for intracellular transport. Cell 2003, 112, 467–480. [Google Scholar]
  74. 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), Beijing, China, 20–23 April 2009; pp. 169–176.
  75. Klann, M. Development of a Stochastic Multi-Scale Simulation Method for the Analysis of Spatiotemporal Dynamics in Cellular Transport and Signaling Processes. Ph.D. Dissertation, Universität Stuttgart, Germany, November 2011. [Google Scholar]
  76. Li, G.; Cui, Q. Mechanochemical coupling in myosin: A theoretical analysis with molecular dynamics and combined QM/MM reaction path calculations. J. Phys. Chem. B 2004, 108, 3342–3357. [Google Scholar]
  77. Kawakubo, T.; Okada, O.; Minami, T. Molecular dynamics simulations of evolved collective motions of atoms in the myosin motor domain upon perturbation of the ATPase pocket. Biophys. Chem 2005, 115, 77–85. [Google Scholar]
  78. Otten, M.; Nandi, A.; Arcizet, D.; Gorelashvili, M.; Lindner, B.; Heinrich, D. Local motion analysis reveals impact of the dynamic cytoskeleton on intracellular subdiffusion. Biophys. J 2012, 102, 758–767. [Google Scholar]
  79. Gershon, N.D.; Porter, K.R.; Trus, B.L. The cytoplasmic matrix: Its volume and surface area and the diffusion of molecules through it. Proc. Natl. Acad. Sci. USA 1985, 82, 5030–5034. [Google Scholar]
  80. Klann, M.T.; Lapin, A.; Reuss, M. Agent-based simulation of reactions in the crowded and structured intracellular environment: Influence of mobility and location of the reactants. BMC Syst. Biol 2011, 5. [Google Scholar] [CrossRef]
  81. Smoluchowski, M. Versuch einer mathematischen Theorie der Koagulationskinetik kolloider Lösungen. Z. Phys. Chem 1917, 92, 129–168. [Google Scholar]
  82. Zhang, Y.; McCammon, J.A. Studying the affinity and kinetics of molecular association with molecular-dynamics simulation. J. Chem. Phys 2003, 118, 1821:1–1821:7. [Google Scholar]
  83. Gabdoulline, R.R.; Wade, R.C. Biomolecular diffusional association. Curr. Opin. Struct. Biol 2002, 12, 204–213. [Google Scholar]
  84. Northrup, S.H.; Erickson, H.P. Kinetics of protein-protein association explained by Brownian dynamics computer simulation. Proc. Natl. Acad. Sci. USA 1992, 89, 3338–3342. [Google Scholar]
  85. Rice, S.A. Diffusion-Limited Reactions; Elsevier: Amsterdam, The Netherlands, 1985. [Google Scholar]
  86. Collins, F.C.; Kimball, G.E. Diffusion-controlled reaction rates. J. Colloid. Sci 1949, 4, 425– 437. [Google Scholar]
  87. Ellis, R.J. Macromolecular crowding: An important but neglected aspect of the intracellular environment. Curr. Opin. Struct. Biol 2001, 11, 114–119. [Google Scholar]
  88. Minton, A.P. The influence of macromolecular crowding and macromolecular confinement on biochemical reactions in physiological media. J. Biol. Chem 2001, 276, 10577–10580. [Google Scholar]
  89. Al-Habori, M. Microcompartmentation, metabolic channelling and carbohydrate metabolism. Int. J. Biochem. Cell Biol 1995, 27, 123–132. [Google Scholar]
  90. Schnell, S.; Turner, T.E. Reaction kinetics in intracellular environments with macromolecular crowding: Simulations and rate laws. Prog. Biophys. Mol. Biol 2004, 85, 235–260. [Google Scholar]
  91. Grima, R.; Schnell, S. A systematic investigation of the rate laws valid in intracellular environments. Biophys. Chem 2006, 124, 1–10. [Google Scholar]
  92. Nicolau, D.V., Jr; Burrage, K. Stochastic simulation of chemical reactions in spatially complex media. Comput. Mathe. Appl 2008, 55, 1007–1018. [Google Scholar]
  93. Means, S.; Smith, A.J.; Shepherd, J.; Shadid, J.; Fowler, J.; Wojcikiewicz, R.J.H.; Mazel, T.; Smith, G.D.; Wilson, B.S. Reaction diffusion modeling of calcium dynamics with realistic ER geometry. Biophys. J 2006, 91, 537–557. [Google Scholar]
  94. Bergdorf, M.; Sbalzarini, I.F.; Koumoutsakos, P. A Lagrangian particle method for reaction–diffusion systems on deforming surfaces. J. Math. Biol 2010, 61, 649–663. [Google Scholar]
  95. Loverdo, C.; Benichou, O.; Moreau, M.; Voituriez, R. Enhanced reaction kinetics in biological cells. Nat. Phys 2008, 4, 134–137. [Google Scholar]
  96. Chaudhuri, A.; Bhattacharya, B.; Gowrishankar, K.; Mayor, S.; Rao, M. Spatiotemporal regulation of chemical reactions by active cytoskeletal remodeling. Proc. Natl. Acad. Sci. USA 2011, 108, 14825–14830. [Google Scholar]
  97. Hardt, S.L. Rates of diffusion controlled reactions in one, two and three dimensions. Biophys. Chem 1979, 10, 239–243. [Google Scholar]
  98. Torney, D.C.; McConnel, H.M. Diffusion-Limited reaction rate theory for two-dimensional systems. Proc. R. Soc. Lond. A 1983, 387, 147–170. [Google Scholar]
  99. Kholodenko, B.N.; Hoek, J.B.; Westerhoff, H.V. Why cytoplasmic signalling proteins should be recruited to cell membranes. Trends Cell Biol 2000, 10, 173–178. [Google Scholar]
  100. Lampoudi, S.; Gillespie, D.T.; Petzold, L.R. Effect of excluded volume on 2D discrete stochastic chemical kinetics. J. Comp. Phys 2009, 228, 3656–3668. [Google Scholar]
  101. Bisswanger, H. Enzyme Kinetics; Wiley VCH: Chichester, England, 2002. [Google Scholar]
  102. Rohwer, J.; Hanekom, A.; Hofmeyr, J.H. A Universal Rate Equation for Systems Biology. Proceedings of the 2nd International Symposium on Experimental Standard Conditions of Enzyme Characterizations (ESEC 2006), Rüdesheim, Germany, 19–23 March 2006; pp. 175–187.
  103. Segel, L.A.; Slemrod, M. The quasi-steady-state assumption: A case study in perturbation. SIAM Rev 1989, 31, 446–477. [Google Scholar]
  104. Byrne, M.J.; Waxham, M.N.; Kubota, Y. Cellular dynamic simulator: An event driven molecular simulation environment for cellular physiology. Neuroinformatics 2010, 8, 63–82. [Google Scholar]
  105. Arjunan, S.N.V.; Tomita, M. A new multicompartmental reaction-diffusion modeling method links transient membrane attachment of E. coli MinE to E-ring formation. M. Syst. Synth. Biol 2010, 4, 35–53. [Google Scholar]
  106. Klann, M.; Ganguly, A.; Koeppl, H. Improved Reaction Scheme for Spatial Stochastic Simulations with Single Molecule Detail. Proceedings of the 8th International Workshop on Computational Systems Biology (WCSB 2011), Zurich, Switzerland, 6–8 June 2011.
  107. Falk, M.; Klann, M.; Ott, M.; Ertl, T.; Koeppl, H. Parallelized Agent-Based Simulation on CPU and Graphics Hardware for Spatial and Stochastic Models in Biology. Proceedings of the 9th International Conference on Computational Methods in Systems Biology, Paris, France, 21–23 September 2011; pp. 73–82.
  108. Clifford, P.; Green, N.J.B. On the simulation of the Smoluchowski boundary condition and the interpolation of brownian paths. Mol. Phys 1986, 57, 123–128. [Google Scholar]
  109. Lapin, A.; Klann, M.; Reuss, M. Stochastic Simulations of 4D Spatial Temporal Dynamics of Signal Transduction Processes. Proceedings of the FOSBE, Stuttgart, Germany, 9–12 September 2007; pp. 421–425.
  110. Van Zon, J.S.; ten Wolde, P.R. Simulating biochemical networks at the particle level and in time and space: Green’s function reaction dynamics. Phys. Rev Lett 2005, 94. [Google Scholar] [CrossRef]
  111. Morelli, M.J.; ten Wolde, P.R. Reaction Brownian dynamics and the effect of spatial fluctuations on the gain of a push-pull network. J. Chem. Phys 2008, 129. [Google Scholar] [CrossRef]
  112. Kim, H.; Shin, K.J. Exact solution of the reversible diffusion-influenced reaction for an isolated pair in three dimensions. Phys. Rev. Lett 1999, 82, 1578–1581. [Google Scholar]
  113. Gopich, I.V.; Szabo, A. Kinetics of reversible diffusion influenced reactions: The self-consistent relaxation time approximation. J. Chem. Phys 2002, 117, 507:1–507:11. [Google Scholar]
  114. Lapin, A.; Klann, M.; Reuss, M. Multi-Scale Spatio-Temporal Modeling: Lifelines of Microorganisms in Bioreactors and Tracking Molecules in Cells. In Biosystems Engineering II; Springer-Verlag: Berlin, Germany, 2010; Volume 121, pp. 23–43. [Google Scholar]
  115. Park, S.; Agmon, N. Theory and simulation of diffusion-controlled michaelis-menten kinetics for a static enzyme in solution. J. Phys. Chem. B 2008, 112, 5977–5987. [Google Scholar]
  116. Pogson, M.; Smallwood, R.; Qwarnstrom, E.; Holcombe, M. Formal agent-based modelling of intracellular chemical interactions. Biosystems 2006, 85, 37–45. [Google Scholar]
  117. Pogson, M.; Holcombe, M.; Smallwood, R.; Qwarnstrom, E. Introducing spatial information into predictive NF-κB modelling—An agent-based approach. PLoS One 2008, 3. [Google Scholar] [CrossRef]
  118. Lipková, J.; Zygalakis, K.C.; Chapman, S.J.; Erban, R. Analysis of Brownian dynamics simulations of reversible bimolecular reactions. SIAM J. Appl. Math 2011, 71, 714–730. [Google Scholar]
  119. Andrews, S.S.; Addy, N.J.; Brent, R.; Arkin, A.P. Detailed simulations of cell biology with Smoldyn 2.1. PLoS Comp. Biol 2010, 6. [Google Scholar] [CrossRef]
  120. Andrews, S.S. Serial rebinding of ligands to clustered receptors as exemplified by bacterial chemotaxis. Phys. Biol 2005, 2, 111–122. [Google Scholar]
  121. Berg, O.G. On diffusion-controlled dissociation. Chem. Phys 1978, 31, 47–57. [Google Scholar]
  122. Klann, M.; Koeppl, H. Escape times and geminate recombinations in spatial simulations of chemical reactions. Biophys. J 2012.
  123. Wade, R.C.; Luty, B.A.; Demchuk, E.; Madura, J.D.; Davis, M.E.; Briggs, J.M.; McCammon, J.A. Simulation of enzyme–substrate encounter with gated active sites. Nat. Struct. Mol. Biol 1994, 1, 65–69. [Google Scholar]
  124. Shoup, D.; Lipari, G.; Szabo, A. Diffusion-controlled bimolecular reaction rates. The effect of rotational diffusion and orientation constraints. Biophys. J 1981, 36, 697–714. [Google Scholar]
  125. Dudko, O.K.; Berezhkovskii, A.M.; Weiss, G.H. Rate constant for diffusion-influenced ligand binding to receptors of arbitrary shape on a cell surface. J. Chem. Phys 2004, 121, 1562–1565. [Google Scholar]
  126. Traytak, S.D. Diffusion-controlled reaction rate to an active site. Chem. Phys 1995, 192, 1–7. [Google Scholar]
  127. Wu, Y.T.; Nitsche, J.M. On diffusion-limited site-specific association processes for spherical and nonspherical molecules. Chem. Eng. Sci 1995, 50, 1467–1487. [Google Scholar]
  128. Bongini, L.; Fanelli, D.; Piazza, F.; de los Rios, P.; Sanner, M.; Skoglund, U. A dynamical study of antibody–antigen encounter reactions. Phys. Biol 2007, 4, 172–180. [Google Scholar]
  129. Pettré, J.; Ciechomski, P.H.; Maïm, J.; Yersin, B.; Laumond, J.P.; Thalmann, D. Real-time navigating crowds: Scalable simulation and rendering. Comput. Animat. Virtual Worlds 2006, 17, 445–455. [Google Scholar]
  130. Behringer, H.; Eichhorn, R. Hard-wall interactions in soft matter systems: Exact numerical treatment. Phys. Rev. E 2011, 83, 065701:1–065701:4. [Google Scholar]
  131. Press, W.H.; Flannery, B.P.; Teukolsky, S.A.; Vetterling, W.T. Numerical Recipes; Cambridge University Press: Cambridge, UK, 1986; Volume 547. [Google Scholar]
  132. Trotter, H.F. An elementary proof of the central limit theorem. Arch. Math 1959, 10, 226–234. [Google Scholar]
  133. Dematte, L. Smoldyn on graphics processing units: Massively parallel brownian dynamics simulation. IEEE/ACM Trans. Comput. Biol. Bioinf 2012, 9, 655–667. [Google Scholar]
  134. Jilkine, A.; Angenent, S.B.; Wu, L.F.; Altschuler, S.J. A density-dependent switch drives stochastic clustering and polarization of signaling molecules. PLoS Comp. Biol 2011, 7. [Google Scholar] [CrossRef]
  135. Plimpton, S.; Slepoy, A. ChemCell: A Particle-Based Model of Protein Chemistry and Diffusion in Microbial Cells; Sandia National Laboratories Technical Report 2003-4509; Sandia National Laboratories: Albuquerque, USA, 2003. [Google Scholar]
  136. Plimpton, S.J.; Slepoy, A. Microbial cell modeling via reacting diffusive particles. J. Phys. Conf. Ser 2005, 16. [Google Scholar] [CrossRef]
  137. Takahashi, K.; Ishikawa, N.; Sadamoto, Y.; Sasamoto, H.; Ohta, S.; Shiozawa, A.; Miyoshi, F.; Naito, Y.; Nakayama, Y.; Tomita, M. E-Cell 2: Multi-platform E-Cell simulation system. Bioinformatics 2003, 19, 1727–1729. [Google Scholar]
  138. Tomita, M.; Hashimoto, K.; Takahashi, K.; Shimizu, T.S.; Matsuzaki, Y.; Miyoshi, F.; Saito, K.; Tanida, S.; Yugi, K.; Venter, J.C.; et al. E-CELL: Software environment for whole-cell simulation. Bioinformatics 1999, 15, 72–84. [Google Scholar]
  139. Takahashi, K.; Kaizu, K.; Hu, B.; Tomita, M. A multi-algorithm, multi-timescale method for cell simulation. Bioinformatics 2004, 20, 538–546. [Google Scholar]
  140. Takahashi, K.; Tănase-Nicola, S.; ten Wolde, P.R. Spatio-temporal correlations can drastically change the response of a MAPK pathway. Proc. Natl. Acad. Sci. USA 2010, 107, 2473–2478. [Google Scholar]
  141. Van Zon, J.S.; Morelli, M.J.; Tanase-Nicola, S.; tenWolde, P.R. Diffusion of transcription factors can drastically enhance the noise in gene expression. Biophys. J 2006, 91, 4350–4367. [Google Scholar]
  142. Stiles, J.R.; Bart, T.M. Monte Carlo Methods for Simulating Realistic Synaptic Microphysiology Using MCell. In Computational Neuroscience—Realistic Modeling for Experimentalists; de Schutter, E., Ed.; CRC Press: Boca Raton, FL, USA, 2001. [Google Scholar]
  143. Hattne, J.; Fange, D.; Elf, J. Stochastic reaction-diffusion simulation with MesoRD. Bioinformatics 2005, 21, 2923–2924. [Google Scholar]
  144. Elf, J.; Doncic, A.; Ehrenberg, M. Mesoscopic reaction-diffusion in intracellular signaling. Proc. SPIE 2003, 5110, 114–124. [Google Scholar]
  145. Ander, M.; Beltrao, P.; di Ventura, B.; Ferkinghoff-Borg, J.; Foglierini, M.; Kaplan, A.; Lemerle, C.; Tomas-Oliveira, I.; Serrano, L. SmartCell, a framework to simulate cellular processes that combines stochastic approximation with diffusion and localisation: Analysis of simple networks. Syst. Biol 2004, 1, 129–138. [Google Scholar]
  146. Wils, S.; de Schutter, E. STEPS: Modeling and simulating complex reaction-diffusion systems with Python. Frontiers Neuroinf 2009, 3. [Google Scholar] [CrossRef]
  147. Stoma, S.; Fröhlich, M.; Gerber, S.; Klipp, E. STSE: Spatio-temporal simulation environment dedicated to biology. BMC Bioinf 2011, 12. [Google Scholar] [CrossRef]
  148. Moraru, I.I.; Schaff, J.C.; Slepchenko, B.M.; Loew, L.M. The virtual cell. Ann. N. Y. Acad. Sci 2002, 971, 595–596. [Google Scholar]
  149. Slepchenko, B.M.; Schaff, J.C.; Macara, I.; Loew, L.M. Quantitative cell biology with the virtual cell. Trends Cell Biol 2003, 13, 570–576. [Google Scholar]
  150. Klann, M.; Koeppl, H.; Reuss, M. Spatial modeling of vesicle transport and the cytoskeleton: The challenge of hitting the right road. PLoS One 2012, 7. [Google Scholar] [CrossRef]
  151. Shillcock, J.C.; Lipowsky, R. The computational route from bilayer membranes to vesicle fusion. J. Phys. Condens. Matt 2006, 18. [Google Scholar] [CrossRef]
  152. Liou, S.H.; Chen, C.C. Cellular ability to sense spatial gradients in the presence of multiple competitive ligands. Phys. Rev. E 2012, 85, 011904:1–011904:5. [Google Scholar]
  153. Angermann, B.R.; Klauschen, F.; Garcia, A.D.; Prustel, T.; Zhang, F.; Germain, R.N.; Meier-Schellersheim, M. Computational modeling of cellular signaling processes embedded into dynamic spatial contexts. Nat. Methods 2012, 9, 283–289. [Google Scholar]
  154. Jeschke, M.; Uhrmacher, A.M. Multi-Resolution Spatial Simulation for Molecular Crowding. Proceedings of the 2008 Winter Simulation Conference, Miami, FL, USA, 7–10 December 2008; pp. 1384–1392.
  155. Pahle, J. Biochemical simulations: Stochastic, approximate stochastic and hybrid approaches. Brief. Bioinf 2009, 10, 53–64. [Google Scholar]
  156. Chatterjee, A.; Vlachos, D.G. Multiscale spatial monte carlo simulations: Multigriding, computational singular perturbation, and hierarchical stochastic closures. J. Chem. Phys 2006, 124, 064110. [Google Scholar]
  157. Jeschke, M.; Ewald, R.; Park, A.; Fujimoto, R.; Uhrmacher, A.M. A parallel and distributed discrete event approach for spatial cell-biological simulations. ACM SIGMETRICS Perform. Eval. Rev 2008, 35, 22–31. [Google Scholar]
  158. Xing, F.; Yao, Y.P.; Jiang, Z.W.; Wang, B. Fine-grained parallel and distributed spatial stochastic simulation of biological reactions. Adv. Mater. Res 2012, 345, 104–112. [Google Scholar]
  159. Rodríguez, J.V.; Kaandorp, J.A.; Dobrzyński, M.; Blom, J.G. Spatial stochastic modelling of the phosphoenolpyruvate-dependent phosphotransferase (PTS) pathway in Escherichia coli. Bioinformatics 2006, 22, 1895–1901. [Google Scholar]
  160. Lampoudi, S.; Gillespie, D.T.; Petzold, L.R. The multinomial simulation algorithm for discrete stochastic simulation of reaction-diffusion systems. J. Chem. Phys 2009, 130, 094104. [Google Scholar]
  161. Stundzia, A.B.; Lumsden, C.J. Stochastic simulation of coupled reaction-diffusion processes. J. Comp. Phys 1996, 127, 196–207. [Google Scholar]
  162. Anderson, D.F. A modified next reaction method for simulating chemical systems with time dependent propensities and delays. J. Chem. Phys 2007, 127, 214107:1–214107:10. [Google Scholar]
  163. Gillespie, D.T. A diffusional bimolecular propensity function. J. Chem. Phys 2009, 131, 164109:1–164109:13. [Google Scholar]
  164. Fange, D.; Berg, O.G.; Sjöberg, P.; Elf, J. Stochastic reaction-diffusion kinetics in the microscopic limit. Proc. Natl. Acad. Sci. USA 2010, 107, 19820–19825. [Google Scholar]
  165. Gillespie, D.T. Deterministic limit of stochastic chemical kinetics. J. Phys. Chem. B 2009, 113, 1640–1644. [Google Scholar]
  166. Danos, V.; Feret, J.; Fontana, W.; Harmer, R.; Krivine, J. Rule-based Modelling of Cellular Signalling. Proceedings of the Eighteenth International Conference on Concurrency Theory, CONCUR’07, Lisbon, Portugal, 3–8 September 2007, Caires, L., Vasconcelos, V.T., Eds.; Springer: Berlin, Germany, 2007; Volume 4703, pp. 17–41. [Google Scholar]
  167. Faeder, J.R.; Blinov, M.L.; Hlavacek, W.S. Rule-based modeling of biochemical systems with BioNetGen. Methods Mol. Biol 2009, 500, 113–167. [Google Scholar]
  168. Camporesi, F.; Feret, J.; Koeppl, H.; Petrov, T. Automatic Reduction of Stochastic Rules-Based Models in a Nutshell. AIP Conference Proceedings. International Conference of Numerical Analysis and Applied Mathematics (ICNAAM 2010), Rhodes, Greece, 19–25 September 2010; 1281, pp. 1330–1334.
  169. Petrov, T.; Ganguly, A.; Koeppl, H. Model decomposition and stochastic fragments. Theor. Comput. Sci 2012, (in press). [Google Scholar]
  170. Tolle, D.P.; le Novère, N. Meredys, a multi-compartment reaction-diffusion simulator using multistate realistic molecular complexes. BMC Syst. Biol 2010, 4. [Google Scholar] [CrossRef]
  171. Yang, J.; Meng, X.; Hlavacek, W.S. Rule-based modelling and simulation of biochemical systems with molecular finite automata. Syst. Biol. IET 2010, 4, 453–466. [Google Scholar]
  172. Gruenert, G.; Ibrahim, B.; Lenser, T.; Lohel, M.; Hinze, T.; Dittrich, P. Rule-based spatial modeling with diffusing, geometrically constrained molecules. BMC Bioinf 2010, 11. [Google Scholar] [CrossRef]
  173. Kang, A.; Kim, J.H.; Lee, S.; Park, H. Diffusion-influenced reactions involving a reactant with two active sites. J. Chem. Phys 2009, 130, 094507:1–094507:8. [Google Scholar]
  174. Park, S.; Agmon, N. Multisite reversible geminate reaction. J. Chem. Phys 2009, 130, 074507:1–074507:11. [Google Scholar]
  175. Bauler, P.; Huber, G.; Leyh, T.; McCammon, J.A. Channeling by proximity: The catalytic advantages of active site colocalization using Brownian dynamics. J. Phys. Chem. Lett 2010, 1, 1332–1335. [Google Scholar]
  176. Locasale, J.W.; Shaw, A.S.; Chakraborty, A.K. Scaffold proteins confer diverse regulatory properties to protein kinase cascades. Proc. Natl. Acad. Sci. USA 2007, 104. [Google Scholar] [CrossRef]
  177. Zhou, H.X.; Szabo, A. Enhancement of association rates by nonspecific binding to DNA and cell membranes. Phys. Rev. Lett. 2004, 93, 178101:1–178101:4. [Google Scholar]
  178. Halford, S.E. An end to 40 years of mistakes in DNA-protein association kinetics? Biochem. Soc. Trans 2009, 37, 343–348. [Google Scholar]
  179. Zechner, C.; Ruess, J.; Krenn, P.; Pelet, S.; Peter, M.; Lygeros, J.; Koeppl, H. Moment-based inference predicts bimodality in transient gene expression. Proc. Natl. Acad. Sci. USA 2012, (in press). [Google Scholar]
  180. Guo, K.; Shillcock, J.; Lipowsky, R. Treadmilling of actin filaments via Brownian dynamics simulations. J. Chem. Phys 2010, 133. [Google Scholar] [CrossRef]
  181. Mahalik, J.P.; Muthukumar, M. Langevin dynamics simulation of polymer-assisted virus-like assembly. J. Chem. Phys 2012, 136, 135101:1–135101:13. [Google Scholar]
  182. Noguchi, H.; Takasu, M. Fusion pathways of vesicles: A Brownian dynamics simulation. J. Chem. Phys 2001, 115, 9547–9551. [Google Scholar]
  183. Mogilner, A.; Allard, J.; Wollman, R. Cell polarity: Quantitative modeling as a tool in cell biology. Science 2012, 336, 175–179. [Google Scholar]
  184. Mugler, A.; Bailey, A.G.; Takahashi, K.; Wolde, P.R. Membrane clustering and the role of rebinding in biochemical signaling. Biophys. J 2012, 102, 1069–1078. [Google Scholar]
  185. Cichocki, B.; Ekiel-Jezewska, M.L.; Wajnryb, E. Communication: Translational Brownian motion for particles of arbitrary shape. J. Chem. Phys 2012, 136, 071102:1–071102:4. [Google Scholar]
  186. Lee, D.; Redfern, O.; Orengo, C. Predicting protein function from sequence and structure. Nat. Rev. Mol. Cell Biol 2007, 8, 995–1005. [Google Scholar]
  187. Zhou, H.; Skolnick, J. Protein structure prediction by pro-Sp3-TASSER. Biophys. J 2009, 96, 2119–2127. [Google Scholar]
  188. Accelrys Discovery Studio. Available online: http://accelrys.com/products/discovery-studio/index.html accessed on 8 June 2012.
  189. Molsoft. Available online: http://www.molsoft.com accessed on 8 June 2012.
  190. García de la Torre, J.; Huertas, M.L.; Carrasco, B. Calculation of hydrodynamic properties of globular proteins from their atomic-level structure. Biophys. J 2000, 78, 719–730. [Google Scholar]
  191. Moal, I.H.; Bates, P.A. Kinetic rate constant prediction supports the conformational selection mechanism of protein binding. PLoS Comp. Biol 2012, 8. [Google Scholar] [CrossRef]
  192. De Jong, D.H.; Schäfer, L.V.; de Vries, A.H.; Marrink, S.J.; Berendsen, H.J.C.; Grubmüller, H. Determining equilibrium constants for dimerization reactions from molecular dynamics simulations. J. Comp. Chem 2011, 32. [Google Scholar] [CrossRef]
  193. Lee, J.; Yang, S.; Kim, J.; Lee, S. An efficient molecular dynamics simulation method for calculating the diffusion-influenced reaction rates. J. Chem. Phys 2004, 120, 7564–7575. [Google Scholar]
  194. Thomas, A.S.; Elcock, A.H. Direct measurement of the kinetics and thermodynamics of association of hydrophobic molecules from molecular dynamics simulations. J. Phys. Chem. Lett 2011, 2, 19–24. [Google Scholar]
  195. Thomas, A.S.; Elcock, A.H. Direct observation of salt effects on molecular interactions through explicit-solvent molecular dynamics simulations: Differential effects on electrostatic and hydrophobic interactions and comparisons to Poisson-Boltzmann theory. J. Am. Chem. Soc 2006, 128, 7796–7806. [Google Scholar]
  196. Garcia-Viloca, M.; Gao, J.; Karplus, M.; Truhlar, D.G. How enzymes work: Analysis by modern rate theory and computer simulations. Science 2004, 303, 186–195. [Google Scholar]
  197. Gabdoulline, R.R.; Wade, R.C. Simulation of the diffusional association of barnase and barstar. Biophys. J 1997, 72, 1917–1929. [Google Scholar]
  198. Schlosshauer, M.; Baker, D. Realistic protein–protein association rates from a simple diffusional model neglecting long-range interactions, free energy barriers, and landscape ruggedness. Protein Sci 2004, 13, 1660–1669. [Google Scholar]
  199. Peter, E.; Dick, B.; Baeurle, S.A. A novel computer simulation method for simulating the multiscale transduction dynamics of signal proteins. J. Chem. Phys 2012, 136, 124112:1–124112:14. [Google Scholar]
  200. Ghaemmaghami, S.; Huh, W.K.; Bower, K.; Howson, R.W.; Belle, A.; Dephoure, N.; O’Shea, E.K.; Weissman, J.S. Global analysis of protein expression in yeast. Nature 2003, 425, 737–741. [Google Scholar]
  201. Fujioka, A.; Terai, K.; Itoh, R.E.; Aoki, K.; Nakamura, T.; Kuroda, S.; Nishida, E.; Matsuda, M. Dynamics of the Ras/ERK MAPK cascade as monitored by fluorescent probes. J. Biol. Chem 2006, 281, 8917–8926. [Google Scholar]
  202. Gutenkunst, R.N.; Waterfall, J.J.; Casey, F.P.; Brown, K.S.; Myers, C.R.; Sethna, J.P. Universally sloppy parameter sensitivities in systems biology models. PLoS Comp. Biol 2007, 3. [Google Scholar] [CrossRef]
  203. Chen, W.W.; Schoeberl, B.; Jasper, P.J.; Niepel, M.; Nielsen, U.B.; Lauffenburger, D.A.; Sorger, P.K. Input–output behavior of ErbB signaling pathways as revealed by a mass action model trained against dynamic data. Mol. Syst. Biol 2009, 5. [Google Scholar] [CrossRef]
  204. Hengl, S.; Kreutz, C.; Timmer, J.; Maiwald, T. Data-based identifiability analysis of non-linear dynamical models. Bioinformatics 2007, 23, 2612–2618. [Google Scholar]
  205. Koh, G.; Hsu, D.; Thiagarajan, P.S. Component-based construction of bio-pathway models: The parameter estimation problem. Theor. Comput. Sci 2011. [Google Scholar] [CrossRef]
  206. Gopalakrishnan, M.; Forsten-Williams, K.; Nugent, M.A.; Täuber, U.C. Effects of receptor clustering on ligand dissociation kinetics: Theory and simulations. Biophys. J 2005, 89, 3686–3700. [Google Scholar]
  207. Falk, M.; Klann, M.; Reuss, M.; Ertl, T. 3D Visualization of Concentrations from Stochastic Agent-based Signal Transduction Simulations. Proceedings of the IEEE International Symposium on Biomedical Imaging: From Nano to Macro, 2010, Rotterdam, The Netherlands, 14–17 April 2010; pp. 1301–1304.
Figure 1. (a) Visualization of the cytoplasm from a Brownian dynamics simulation including cytoskeleton filaments and just the signaling molecules of one pathway. For visualization, all molecules and cytoskeleton filaments have been replaced by their atomic structure and rendered by raytracing (ScienceVisuals [10,11]); (b) Physiological level of crowding, i.e., a representative molecular size distribution and abundance, modeled either with spheres or real molecule shapes by Ando and Skolnick [12] in a cubic subvolume of the cytoplasm. Reproduced with permission of PNAS. Due to the high density of molecules it is impossible to see through the cytoplasm. These crowding conditions affect diffusion and reactions in the cell.
Figure 1. (a) Visualization of the cytoplasm from a Brownian dynamics simulation including cytoskeleton filaments and just the signaling molecules of one pathway. For visualization, all molecules and cytoskeleton filaments have been replaced by their atomic structure and rendered by raytracing (ScienceVisuals [10,11]); (b) Physiological level of crowding, i.e., a representative molecular size distribution and abundance, modeled either with spheres or real molecule shapes by Ando and Skolnick [12] in a cubic subvolume of the cytoplasm. Reproduced with permission of PNAS. Due to the high density of molecules it is impossible to see through the cytoplasm. These crowding conditions affect diffusion and reactions in the cell.
Ijms 13 07798f1
Figure 2. (a) Two diffusing molecules can collide and will be reflected if they do not react; (b) Corresponding probability density function (pdf) for the distance of j relative to i as described by the Fokker–Planck equation; (c) Reaction probability depending on the initial distance; (d)–(f) Fokker–Planck equation and boundary conditions. The pdf for the distance r between two diffusing molecules as described by (d) starting from W(0) = δ(xixj) is shown in (b). In this description, particles react if they diffuse “into” the reaction partner, which is accounted for by the flux across the collision surface and depends on the microscopic rate constant κ. Therefore the boundary condition (e) is partially reflecting. For κ = 0 (e) becomes completely reflective, describing two inert particles. Note the “blister” in (b) which deforms the normal distribution at the boundary caused by the reflection. Due to incomplete reflection, the total probability ∫ WdV < 1. The loss corresponds to the reaction probability. Hence the reaction probability (c) for a given initial distance is found as P F P r e a c t i o n = 1 - r * W d V [75,109,114].
Figure 2. (a) Two diffusing molecules can collide and will be reflected if they do not react; (b) Corresponding probability density function (pdf) for the distance of j relative to i as described by the Fokker–Planck equation; (c) Reaction probability depending on the initial distance; (d)–(f) Fokker–Planck equation and boundary conditions. The pdf for the distance r between two diffusing molecules as described by (d) starting from W(0) = δ(xixj) is shown in (b). In this description, particles react if they diffuse “into” the reaction partner, which is accounted for by the flux across the collision surface and depends on the microscopic rate constant κ. Therefore the boundary condition (e) is partially reflecting. For κ = 0 (e) becomes completely reflective, describing two inert particles. Note the “blister” in (b) which deforms the normal distribution at the boundary caused by the reflection. Due to incomplete reflection, the total probability ∫ WdV < 1. The loss corresponds to the reaction probability. Hence the reaction probability (c) for a given initial distance is found as P F P r e a c t i o n = 1 - r * W d V [75,109,114].
Ijms 13 07798f2
Table 1. Empiric approximations for the hydrodynamic radius rh based on the molecular weight MW in kDa. (i) is a fit to experimental data, e.g., from [69,70]. The other equations assume that the mass is (re-)distributed in a sphere, for instance with a specific volume of 1 cm3/g in (ii) [41,71]. Due to the in general nonspherical shape and the “holes” of the molecule, an exponent larger than 1/3 as in (i) is reasonable. The hydrodynamic radii reported by [60] fall between (i) and (ii).
Table 1. Empiric approximations for the hydrodynamic radius rh based on the molecular weight MW in kDa. (i) is a fit to experimental data, e.g., from [69,70]. The other equations assume that the mass is (re-)distributed in a sphere, for instance with a specific volume of 1 cm3/g in (ii) [41,71]. Due to the in general nonspherical shape and the “holes” of the molecule, an exponent larger than 1/3 as in (i) is reasonable. The hydrodynamic radii reported by [60] fall between (i) and (ii).
Hydrodynamic Radius [nm]Reference
(i)rh = 0.6169 × MW0.4226[72]
(ii)rh = 0.7468 × MW1/3[41]
(iii)rh = 0.5429 × MW1/3[40]
Table 2. Spatial simulations on the cellular level.
Table 2. Spatial simulations on the cellular level.
Name/AuthorsFeaturesWebsite/References
Smoldyn
S. Andrews et al.
Particle based simulator for reaction diffusion processes in arbitrarily shaped compartments. (point particles, no crowding). www.smoldyn.org
[28,48,119,120,133,134]

ChemCellParticle based simulator within realistic cell shapes.chemcell.sandia.gov
[93,135,136]

E-CellComplete software environment for simulations on several levels.
Contains further analysis tools.
www.e-cell.org
[137139]
(GFRD,eGFRD) ten Wolde et al.Green’s function reaction dynamics will be included in the E-Cell project[110,111,140,141]

FLAMEAgent-based multi-scale simulation (also beyond the cellular level). www.flame.ac.uk
[32,116,117]

MCellMonte Carlo simulator of reaction diffusion processes. Reactions can only happen at membranes www.mcell.cnl.salk.edu
[142]

MesoRDSpatial derivative of Gillespie’s algorithm to solve the Reaction-Diffusion Master Equation (RDME) with the “next subvolume method”mesord.sourceforge.net
[143,144]

SmartCell
Serrano et al.
Spatial derivative of Gillespie’s algorithm in arbitrarily shaped compartments.software.crg.es/smartcell
[145]

STEPSTetrahedral mesh based spatial derivative of Gillespie’s algorithmsteps.sourceforge.net/STEPS
[146]

STSE
S.Stoma
PDE based simulator with compartments and direct linking to microscope images. www.stse-software.org
[147]

V CellODE/PDE or SDE based simulator within realistic cell shapes. www.nrcam.uchc.edu
[148,149]

M. Klann et al.Agent-based Brownian dynamics including cytoskeleton, crowding and vesicle transport. www.bison.ethz.ch/research/spatial simulations
[75,80,106,107,122,150]

Share and Cite

MDPI and ACS Style

Klann, M.; Koeppl, H. Spatial Simulations in Systems Biology: From Molecules to Cells. Int. J. Mol. Sci. 2012, 13, 7798-7827. https://doi.org/10.3390/ijms13067798

AMA Style

Klann M, Koeppl H. Spatial Simulations in Systems Biology: From Molecules to Cells. International Journal of Molecular Sciences. 2012; 13(6):7798-7827. https://doi.org/10.3390/ijms13067798

Chicago/Turabian Style

Klann, Michael, and Heinz Koeppl. 2012. "Spatial Simulations in Systems Biology: From Molecules to Cells" International Journal of Molecular Sciences 13, no. 6: 7798-7827. https://doi.org/10.3390/ijms13067798

APA Style

Klann, M., & Koeppl, H. (2012). Spatial Simulations in Systems Biology: From Molecules to Cells. International Journal of Molecular Sciences, 13(6), 7798-7827. https://doi.org/10.3390/ijms13067798

Article Metrics

Back to TopTop