Abstract
Background
Realistic biochemical simulators aim to improve our understanding of many biological processes that would be otherwise very difficult to monitor in experimental studies. Increasingly accurate simulators may provide insights into the regulation of biological processes due to stochastic or spatial effects.
Results
We have developed GridCell as a threedimensional simulation environment for investigating the behaviour of biochemical networks under a variety of spatial influences including crowding, recruitment and localization. GridCell enables the tracking and characterization of individual particles, leading to insights on the behaviour of low copy number molecules participating in signaling networks. The simulation space is divided into a discrete 3D grid that provides ideal support for particle collisions without distance calculation and particle search. SBML support enables existing networks to be simulated and visualized. The user interface provides intuitive navigation that facilitates insights into species behaviour across spatial and temporal dimensions. We demonstrate the effect of crowing on a MichaelisMenten system.
Conclusion
GridCell is an effective stochastic particle simulator designed to track the progress of individual particles in a threedimensional space in which spatial influences such as crowding, colocalization and recruitment may be investigated.
Background
One of the main goals of computational cell biology aims to accurately simulate large biological systems at molecular resolution. Stochastic effects and spatial constraints are increasingly being recognized as important factors in the normal functioning of molecular networks [1]. The efficiency of biochemical networks is enhanced by component colocalization [2], and certain signaling networks are thought to be facilitated by transport and colocalization [3]. In addition, molecular crowding has been shown to affect biochemical systems [46]. Modeling and simulation of these kinds of networks requires new kinds of stochastic simulators.
We developed GridCell to simulate biological models with specific consideration for stochasticity, locality, and collision. GridCell is based on a simplified model for molecular movement and interaction. It uses a discrete threedimensional cubic grid based on the D3Q27 model often used in the application of the LatticeBoltzmann Method (LBM) [7]. Each voxel has access to itself and its 26 neighbours and is independent of voxels outside this immediate surrounding. Figure 1 shows the 27 possible locations accessible to a voxel from a D3Q27 grid. The integeraddressed 3D grid avoids floatingpoint computation and distance calculations, resulting in an efficient implementation. Molecules are represented as particles that move and react stochastically within discrete volumes in discrete timesteps. Collisions and molecular crowding are enforced since only one particle can occupy a given location at any time. GridCell stores the coordinates of all the particles on the 3D grid at every turn, thereby enabling particle tracking in both space and time.
Figure 1. D3Q27 cubic grid structure. The 27 possible locations of a D3Q27 cubic grid that a given voxel can access.
The simulation space is visualized via a 3D interface and 2D graphs, and surface plots summarize molecule concentrations over space and time. GridCell supports models specified by the Systems Biology Markup Language (SBML) (please see Availability & requirements for further information). SBML models can be obtained from public repositories such as EBI's Biomodels database (please see Availability & requirements for further information) or designed using software such as SBMLeditor or JDesigner (please see Availability & requirements for further information).
Implementation
Algorithm
The simulation employs a twophase process in which particles (1) attempt to move and then (2) attempt to react every turn.
Movement Phase
A particle can move at most once per timestep. Since a particle only has access to its immediate surrounding, a particle can only move in one of 27 nearest locations, including the current location. The selection of the movement direction is made randomly; therefore the particles follow a Brownian random walk. Figure 2 shows an example of 4 different Brownian random walks of particles starting from the same location. In GridCell, any particle attempting to move to an occupied location will generate a collision. A collision prevents the particle attempting to move from moving during that turn and does not affect the other particle.
Figure 2. Random Brownian walk. Random walks of 4 different particles in GridCell after 1000 timesteps.
Diffusion
Particles following a Brownian random walk should also follow the wellknown EinsteinSmoluchowski equation
where <r^{2}> is the meansquare displacement, d is the dimensionality, D the diffusion coefficient and t the elapsed time. Figure 3 shows the meansquare displacement in units of voxel <v^{2}> (averaged over 1000 iterations) versus the number of timesteps n_{ts }when the probability of movement of the particles at every timestep is equal to 1. As expected for an uncrowded case, the meansquare displacement increases linearly with the number of timesteps. This leads to the following relation
Figure 3. Diffusion in GridCell. Meansquare distance (over 1000 iterations) of particles versus the number of timesteps in GridCell.
where A is the slope of the graph. By substituting and where s_{vox }is the length of the sides of the voxels in meters and l_{ts }is the length of the timestep in seconds, we get
Since the probability of movement at each timestep of the particle is equal to 1, D can be substituted for the maximum diffusion speed D_{max }that GridCell can support for a given timestep and voxel size. This upper limit on diffusion speed is caused by the design decision of restraining particle movement to its immediate neighbourhood (the D3Q27 grid). By calculating the slope of the graph and setting the dimensionality d equal to 3, D_{max }can be calculated as
Smaller diffusion speeds are simulated by applying a different probability of movement such that
where p_{m }is the probability of movement of a particle at every timestep. As long as the diffusion speeds of the particles are smaller than D_{max}, diffusion will be modeled correctly. If a larger diffusion speed is needed, one can reduce the timestep or increase the size of the voxels.
Reaction Phase
A particle may react only once per turn and only with its immediate surrounding. The reaction phase is completely independent from the movement phase, therefore it does not matter if a particle previously moved or collided with another particle. Common interactions include aggregation events such as molecular complex formation/dissolution or conversion events such as chemical reactions. Only the simplest reactions involving 3 or less particles are directly supported. Complex reactions involving more than 3 particles are decomposed into several elementary reactions. The probability of reaction per timestep is derived from the overall rate of reaction and is very similar to the approach taken by ChemCell [8]. Only 3 different reactions involve 3 or less participants: 1 reactant and 1 product, 1 reactant and 2 products and 2 reactants and 1 product. Let's consider the two reactions that involve a single reactant
Both reactions have a forward rate of reaction k in units of time^{1 }and timestep is t second. Assuming a wellmixed approximation and N particles of type A in the system, then in both cases the expected number of reaction per turn is given by N(1  e^{kt}). Considering each particle individually, each particle has a probability equal to 1  e^{kt }to react during each timestep. In our stochastic model, a uniform random number R_{n }between 0 and 1 is generated for each particle, and the reaction occurs if R_{n }< 1  e^{kt}. In a reaction with only 1 reactant and 1 product, the reactant is simply replaced by the product. In a reaction with 1 reactant and 2 products, a search is first conducted in the surrounding area. If there is at least 1 free voxel in the surrounding area of the particle, the reaction takes place, and the second product is positioned in that free location while the first product is placed at the position of the initial reactant. The reaction is blocked if no free position is found. This limitation only modifies the overall reaction rate of the reaction in a situation where the whole cell is completely filled which would prevent any movement and reaction to take place.
Consider the following reaction with 2 reactants:
with a rate constant k in units of (molarity*time)^{1 }and a timestep between each iteration of t second. Assuming N_{a }particle of type A, N_{b }particle of type B, a Volume V and the Avogadro's number being A_{v}, then the total number of reactions N_{r }in a wellmixed system is given by
On average, the desired number of reactions in our system should be equivalent to the result of the above equation. In our system, particles can only react with their immediate surrounding locations. In a wellmixed system, the number of A,B pairs that are close enough to each other to generate a reaction is given by N = N_{a}N_{b}V_{c}/V where V_{c }is the volume of the cube containing the 26 "neighbouring" voxels and V is still the total volume of the simulation. If each of those pairs react with probability P, then N_{r }= NP. Setting the 2 equations N_{r }= NP = kN_{a}N_{b}t/(A_{v}V) gives the equation
The formula is independent of V, N_{a }and N_{b }as expected. Also, for a given rate constant k, it is possible to have a set of parameter t/V_{c }such that P is greater than 1. If that is the case, a smaller timestep or larger voxel size (proportional to V_{c}) has to be selected. A random number R_{n }is generated. If R_{n }<P, then the first reactant will search its surrounding area for the second reactant. If it is found, the reaction takes place and the product is placed at the location of the first reactant. If no reactant is found, the reaction is aborted. Note that the second reactant does not try to react with the first one, doing so still enables us to get the right rate of reaction and reduces the number of operations needed to update the simulation. In the case where a particle can participate in many different reactions, a random number is generated to select which reaction is going to be tested first. If the first reaction does not take place, then the next reaction on the list is tested. The procedure will go on until either a reaction takes place or all possible reactions have been tried. This ensures that, on average, all reactions are tested equally.
When 2 reactants of the same species form a product such as
the individual rate of reaction of particle A needs to be modified to ensure that the overall rate of reaction is respected since each reactant will attempt to react with the other one. Assuming y is the overall probability of reaction and x is the individual probability of reaction of species A, then
More complex reactions are implemented by creating a cascade of several elementary equations. This process, done automatically by the software, will break the complex reactions into a series of simpler reactions by introducing "temporary" species. For example, consider the following reaction with 1 reactant and 5 products.
where k is the rate of reaction in time^{1}. For each product exceeding 2, a temporary species is created. In this case, 3 temporary species are created. It follows that the reaction is broken down into:
where T_{1}, T_{2 }and T_{3 }are respectively the first, second and third temporary species. By setting the rate of reaction of equation 17 equal to k and the probability of reaction of equations containing any temporary species on the reactant side equal to 1, we reduce the artifacts due to the creation of the "temporary" species to a minimum. Indeed, the temporary species disappear from the system as quickly as possible and the overall rate of reaction is identical.
Shown below is the case where more than 2 reactants merge into a single product:
The procedure is similar to the previous case, 1 temporary species is created for each reactant above 2.
where T_{1}, T_{2 }and T_{3 }are respectively the first, second and third temporary species. In order to obtain the same overall probability of reaction and to reduce the impact of the temporary species on the system to a minimum, the probability of reaction of any reaction containing temporary species on the reactant side (equation 21 and 22) is set to 1. Assuming that P is the probability of reaction of the reaction presented in equation 18 and P_{1 }and P_{2 }are the probability of the first and second simple reactions A + B → T_{1 }and C + D → T_{2 }then, we set P = P_{1}P_{2}. We also set P_{1 }= P_{2}. Equating the 2 equations gives P_{1 }= P_{2 }= .
In general, the probability of the simple reactions P_{n }containing no temporary species is equal to
where P is the probability of reaction and N_{reactants }is the number of reactants of the initial reaction. Each temporary particle has a parameter lifetime which indicates the number of turns the particle has to live in the system before reverting back to its previous state. The short lifetime of temporary particles is important for 2 reasons. First, it makes sure that temporary particles are effectively temporary and never stay in the system for a long period of time. It also makes sure that all the reactants are to be close to each other in order for the reaction to complete. Usually, a lifetime of 2–3 turns is reasonable since it gives enough time to react with the neighbouring particles while making sure temporary particles do not constitute the bulk of the system.
Reversible reactions are handled by creating 2 different separate reactions, 1 for the forward reaction with the forward reaction rate and 1 for the backward reaction with the corresponding backward reaction rate. Assuming the following reaction
with forward reaction rate k_{f }and backward reaction rate k_{b}. This reversible reaction is then split into
with a reaction rate k_{f }and
with reaction rate k_{b}.
Temporary particles involved in a reversible reaction have a flag mentioning if they are participating in a forward or backward reaction such that they can revert back to the proper reactants when their lifetime reaches zero.
Performance analysis
Preliminary tests have been conducted to determine how the software reacts to different system sizes. The tests have been executed on a standalone microprocessor: a 3.2 GHz P4 with 2GB of RAM. The current algorithm is computed serially. As it can be shown in Table 1, the time required to compute a timestep increases linearly with the number of particles and voxels present in the system. Tables 2 and 3 demonstrate how the performance is affected by independently modifying the number of voxels or the number of particles. The maximum number of particles that can be currently simulated is equal to the maximum number of voxels that can be supported, which is 10^{7}. Table 4 shows that the number of reactions occurring at each timestep has a negligible effect on the performance. The reason is that all reactions have to be tested, regardless of whether or not they actually react. There are no practical limitations to the number of chemical species or the number of different reactions present in the system beyond the absolute limit on the number of voxels.
Table 1. GridCell performance versus system size
Table 2. GridCell performance versus number of voxels
Table 3. GridCell performance versus number of particles
Table 4. GridCell performance versus the average number of reactions
User Interface Features
The rendering is implemented in OpenGL, and most userinterface functions are written using the PLIB library, which is available online http://plib.sourceforge.net/ webcite. GridCell's user interface (Figure 4) consists of a) a menu system, b) an interactive 3D simulation space, c) a species panel, d) a 2D plot of concentration versus time, and e) a 2D plot of concentration versus space.
Figure 4. GridCell user interface. GridCell user interface with (a) menu, (b) 3D space, (c) species panel, (d) 2D plot of concentration versus time, and (e) 2D surface plot of concentration versus space. Simulated model involves the translocation of particles through a membrane with embedded enzymes.
The menu system (Figure 4a) provides the ability to load SBML models, set parameters and control the simulation. Userdesignated simulation parameters include the number of times to run the simulation, the timestep, the total simulation time, the sampling rate which is the frequency that the 2D graphs are updated and the results saved to file, and the frame rate which designates the frequency of updating the 3D visualization. GridCell computes the means and the standard deviations of the concentration over time if the user chooses to run multiple iterations of the simulation. These preferences may be saved and used later in any simulation. GridCell saves the particle concentrations and the 2D surface plot data in userspecified tabdelimited files. Spatial information such as specific compartment geometries or colocalization of particles is specified in an optional configuration file.
A key feature of the GridCell user interface is the ability to interact with the threedimensional simulation volume (Figure 4b). Users can navigate into the 3D scene with mouse and keyboard controls to rotate, pan and zoom. Buttons are present to i) start/pause simulations, ii) change the particle representation from cubes to points for faster rendering, iii) turn off the visualization for optimal simulation performance, and iv) hide or show all particle types.
The species panel (Figure 4c) contains the current amount of each species, and allows species selection for the visualization plots. A second column specifies which species to render in the 2D surface plot of concentration versus space (Figure 4e). Particle colours are automatically selected from a predefined colour palette.
Finally, two plots to summarize particle concentrations with respect to time (Figure 4d) and space (Figure 4e) are provided in realtime. The 2D spatial plot displays increasing concentration with increasing brightness along a selected Cartesian axis.
Results and discussion
MichaelisMenten reaction
The MichaelisMenten equations are used to describe most enzymatic reactions. Its kinetics is given by the following equation:
The enzyme E reacts with the substrate S to form the complex ES with a rate of reaction k_{1}. ES decomposes into the enzyme E and a new product P with a rate k_{2}, or reverts back to its original form E + S with rate k_{r}.
Crowding
One of the main differences between GridCell and other simulators is its ability to simulate particle crowding. Molecular crowding occurs when particle density affects movement and reactivity. Crowding is typically ignored in most models since kinetics are often based on controlled, in vitro conditions that are not crowded. In addition, simulators do not typically support this feature since it is computationally expensive to keep track of all particle positions and their excluded volume, and to implement collisiondetection algorithms. Some simulators (e.g. Smoldyn [9]) have shown crowding effects by explicitly introducing cubic obstacles [10] in the model. In contrast, GridCell implicitly exhibits molecular crowding effects by allowing interparticle collisions. We demonstrate the effect of crowding by adding inert particles to a MichaelisMenten system. Inert particles do not react with other molecules but reduce their movement and affect the overall number of reactions. The simulation parameters are described in Table 5. Figure 5 shows the number of products over time for a wide range of concentrations of inert particles averaged over 20 iterations. The individual simulations provided almost identical results to one another with a relative standard deviation smaller than 3.5% at the end of the simulation. The number indicated in the legend signifies the percentage of the voxels occupied by inert particles. In this specific example, with a voxel size of 3.2^{20 }litres, this amounts to approximately 30000 inert particles per step of 10%.
Table 5. Simulation parameters
Figure 5. Effect of crowding on MichaelisMenten product formation using GridCell. Effect of increasing the number of inert particles on product formation of a MichaelisMenten system using GridCell. The mean has been calculated over 20 iterations. Percentage of voxels occupied by inert particles.
Interestingly, the maximum rate of reaction is obtained when the inert particles occupy 20% of the volume, which agrees with the fact that macromolecular crowding may enhance reaction rates, as the particles have to search a smaller volume to find each other [11]. However, above 30%, the reaction rates decrease linearly as more and more inert particles are added. Under wellmixed and uncrowded systems, GridCell provides similar results to other ODE simulators and stochastic algorithms such as the stochastic simulation algorithm (SSA) from Gillespie [12] with the exception of a small stochastic noise [13].
Related Work
GridCell is related to a family of Monte Carlo (MC) simulators (Table 6). SmartCell [14] and MesoRD [15] subdivide the simulation space into smaller subvolumes (voxels) that can contain many particles. Each subvolume is composed of a wellmixed solution, and particles can diffuse to adjacent subvolumes. This approach permits quicker simulations, but it is impossible to track individual particles, and molecular crowding has no effect on movement and reaction rates. Cell++ [16] combines a cellular automata engine with Brownian dynamics in order to simulate large quantities of small molecules on a discretized grid, while large molecules exhibit stochastic behaviour and move in a continuous space. Both spaces are then superimposed onto each other, and reactions can take place between the two different spaces. Currently, only collisions between particles (both small and large) and a fixed membrane separating two compartments are supported. Therefore, molecular crowding effects are not simulated. Unlike Cell++, MCell [17] tracks all individual particles in a continuous 3D space, and the diffusing particles follow Brownian dynamics. Particles may collide and interact with effector sites and 2D membrane surfaces, but not with other particles. ChemCell [8] calculates the probability of reaction at every timestep, and particles follow Brownian motion in a continuous space that requires a computationally expensive search algorithm to find nearby particles, and the use of dimensionless particles removes the ability of particles to collide with one another. In contrast, GridCell enables stochasticitybased investigations of SBML networks while considering spatial effects of recruitment, localization and crowding.
Table 6. Spatial simulators
Future Directions
GridCell performance is tightly linked to the number of voxels in the simulation space. The simulator can currently support a maximum of 10^{7 }voxels/particles which is not enough to simulate at a molecular resolution structures as complex as a complete cell, the long term goal of GridCell. However, the simple and regular algorithm of GridCell, which does not require any searches or complex operations, is a prime candidate for acceleration by parallelization to achieve performance speedup and simulate largescale systems.
Conclusion
GridCell is a stochastic simulator that uses a 3D grid and accounts for locality, very low concentration stochastic effects and particle collisions. Its userinterface makes it easy to use while providing several tools to analyze the system. GridCell reproduces the results obtained with ODEs and the Stochastic Simulation Algorithm (SSA) for simple systems when crowding and locality do not affect the system. We also show that particle collisions can have a significant impact on the speed of reaction and that the wellmixed assumption and dimensionless particles can induce a significantly different response in a biological system. The discrete 3D grid and the nearestneighbour interactions remove the need to do any distance calculation, particle search and floatingpoint arithmetic. The regularity and simplicity of the algorithm makes it a good candidate for acceleration with a parallel architecture which will open the door to the simulation of even more complex systems.
Availability and Requirements
The software is available at http://iml.ece.mcgill.ca/GridCell webcite and runs under the Windows XP operating system. This package includes sample SBML and GridCell configuration files. GridCell requires the Systems Biology Markup Language Library (libSBML 2.3.4Xerces; http://sbml.org/software/libsbml/ webcite) and the OpenGL utility toolkit (GLUT 3.7.6; http://www.xmission.com/~nate/glut.html webcite). EBI's Biomodels database: http://www.ebi.ac.uk/biomodels/ webcite. SBMLeditor: http://www.ebi.ac.uk/compneursrv/SBMLeditor.html webcite. JDesigner: http://sbw.kgi.edu/software/jdesigner.htm webcite. Systems Biology Markup Language (SBML): http://sbml.org webcite.
Authors' contributions
LB designed and programmed the GridCell simulator. SAA contributed to the programming of the GridCell simulator and its graphical user interface. MD and WJG contributed to the conception and design of the GridCell simulator. LB and SAA drafted the manuscript and MD and WJG revised the manuscript. All authors read and approved the final manuscript.
Acknowledgements
The authors gratefully acknowledge funding provided by the Natural Sciences and Engineering Research Council of Canada.
References

Lemerle C, Ventura BD, Serrano L: Space as the final frontier in stochastic simulations of biological systems.
FEBS Letters 2005, 578:17891794. Publisher Full Text

Jorgensen K, et al.: Metabolon formation and metabolic channeling in the biosynthesis of plant natural products.
Curr Opin Plant Biol 2005, 8:280291. PubMed Abstract  Publisher Full Text

Kholodenko BN: Fourdimensional organization of protein kinase signaling cascades: the roles of diffusion, endocytosis and molecular motors.
Exp Biol 2003, 206:20732082. Publisher Full Text

Saxton M: Anomalous Subdiffusion in Fluorescence Photobleaching Recovery: A Monte Carlo Study.
Biophysical Journal 2001, 81:22262240. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Turner T, Schnell S, Burrage K: Stochastic approaches for modelling in vivo reactions.
Computational Biology and Chemistry 2004, 28:165178. Publisher Full Text

Schnell S, Turner T: Reaction kinetics in intracellular environments with macromolecular crowding: simulations and rate laws.
Progress in Biophysics and Molecular Biology 2004, 85:235260. Publisher Full Text

Williams S, Carter J, Oliker L, Shalf J, Yelick K: Lattice Boltzmann Simulation Optimization on Leading Multicore Platforms.
International Parallel and Distributed Processing Symposium (IPDPS) 2008.

Plimpton S, Slepoy A: ChemCell: A ParticleBased Model of Protein Chemistry and Diffusion in Microbial Cells.

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

Lipkow K, Andrews SS, Bray D: Simulated diffusion of phosphorylated CheY through the cytoplasm of Escherichia coli.
J Bact 2005, 187:4553. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Zimmerman S, Minton A: Macromolecular crowding: biochemical, biophysical, and physiological consequences.
Annu Rev Biophys Biomol Struct 1993, 22:2765. PubMed Abstract  Publisher Full Text

Gillespie DT: A general method for numerically simulating the stochastic time evolution of couple chemical reactions.
Journal of Computational Physics 1976, 22:403434. Publisher Full Text

Boulianne L, Dumontier M, Gross WJ: A Stochastic ParticleBased Biological System Simulator.
Proceedings of the Summer Computer Simulation Conference, San Diego, California (USA) 2007, 794801.

Ander M, et al.: SmartCell, a framework to simulate cellular processes that combines stochastic approximation with diffusion and localisation: analysis of simple networks.
Syst Biol 2004, 1:129138. Publisher Full Text

Hattne J, Fange D, Elf J: Stochastic reactiondiffusion simulation with MesoRD.
Bioinfomatics 2005, 21:29232924. Publisher Full Text

Sanford C, Yip MLK, White C, Parkinson J: Cell++simulating biochemical pathways.
Bioinfomatics 2006, 22:29182925. Publisher Full Text

Stiles JR, Bartol TM: Monte Carlo methods for simulating realistic synaptic microphysiology using MCell. In Computational Neuroscience: Realistic Modeling for Experimentalists. Edited by DeSchutter E. Boca Raton: CRC Press; 2001:87127.