# GridCell: a stochastic particle-based biological system simulator

- Laurier Boulianne
^{1}, - Sevin Al Assaad
^{1}, - Michel Dumontier
^{2}and - Warren J Gross
^{1}Email author

**2**:66

**DOI: **10.1186/1752-0509-2-66

© Boulianne et al; licensee BioMed Central Ltd. 2008

**Received: **15 February 2008

**Accepted: **23 July 2008

**Published: **23 July 2008

## 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 three-dimensional 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 Michaelis-Menten system.

### Conclusion

GridCell is an effective stochastic particle simulator designed to track the progress of individual particles in a three-dimensional space in which spatial influences such as crowding, co-localization 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 co-localization [2], and certain signaling networks are thought to be facilitated by transport and co-localization [3]. In addition, molecular crowding has been shown to affect biochemical systems [4–6]. Modeling and simulation of these kinds of networks requires new kinds of stochastic simulators.

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 two-phase process in which particles (1) attempt to move and then (2) attempt to react every turn.

#### Movement Phase

#### Diffusion

Particles following a Brownian random walk should also follow the well-known Einstein-Smoluchowski equation

<*r*^{2}> = 2*dDt*,

*r*

^{2}> is the mean-square displacement,

*d*is the dimensionality,

*D*the diffusion coefficient and

*t*the elapsed time. Figure 3 shows the mean-square 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 mean-square displacement increases linearly with the number of timesteps. This leads to the following relation

<*v*^{2}> = *An*_{
ts
}

*A*is the slope of the graph. By substituting $<{v}^{2}>=\frac{<{r}^{2}>}{{s}_{vox}^{2}}$ and ${n}_{ts}=\frac{t}{{l}_{ts}}$ 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

*D*_{
max
}= 0.335*s*_{
vox
}^{2}/*t*_{
s
}.

Smaller diffusion speeds are simulated by applying a different probability of movement such that

*D* = *p*_{
m
}*D*_{
max
}

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

*A* → *B*

*A* → *B* + *C*.

Both reactions have a forward rate of reaction *k* in units of time^{-1} and timestep is *t* second. Assuming a well-mixed 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:

*A* + *B* → *C*

*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 well-mixed system is given by

*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

*A* + *A* → *B*

*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.

*A* → *B* + *C* + *D* + *E* + *F*

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:

*T*_{1} → *B* + *C*

*T*_{2} → *D* + *E*

*T*_{3} → *F* + *T*_{1}

*A* → *T*_{2} + *T*_{3}

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:

*A* + *B* + *C* + *D* + *E* → *F*.

The procedure is similar to the previous case, 1 temporary species is created for each reactant above 2.

*A* + *B* → *T*_{1}

*C* + *D* → *T*_{2}

*E* + *T*_{1} → *T*_{3}

*T*_{2} + *T*_{3} → *F*

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} = $\sqrt{P}$.

*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

*A* + *B* + *C* + *D* + *E* ↔ *F*.

with forward reaction rate *k*_{
f
}and backward reaction rate *k*_{
b
}. This reversible reaction is then split into

*A* + *B* + *C* + *D* + *E* → *F*

with a reaction rate *k*_{
f
}and

*F* → *A* + *B* + *C* + *D* + *E*

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

^{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.

GridCell performance versus system size

Number of Voxels | 1e3 | 1e4 | 1e5 | 1e6 |
---|---|---|---|---|

Number of Particles | 3e2 | 3e3 | 3e4 | 3e5 |

Time (s) | 1.62e-4 | 1.58e-3 | 1.6e-2 | 1.7e-1 |

GridCell performance versus number of voxels

Number of Voxels | 1e3 | 1e4 | 1e5 | 1e6 |
---|---|---|---|---|

Time (s) | 1.6e-5 | 1.6e-4 | 2.14e-3 | 2.06e-2 |

GridCell performance versus number of particles

Number of Particles | 1e3 | 1e4 | 1e5 | 5e5 |
---|---|---|---|---|

Time (s) | 21.3e-2 | 26.4e-2 | 68.1e-2 | 22.8e-1 |

GridCell performance versus the average number of reactions

Average Number of Reactions | 0 | 16.5e3 | 29.5e3 | 39.5e3 | 47.5e3 |
---|---|---|---|---|---|

Time (s) | 7.4e-2 | 7.3e-2 | 7.25e-2 | 7.2e-2 | 7.1e-2 |

### User Interface Features

The menu system (Figure 4a) provides the ability to load SBML models, set parameters and control the simulation. User-designated 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 user-specified tab-delimited files. Spatial information such as specific compartment geometries or co-localization of particles is specified in an optional configuration file.

A key feature of the GridCell user interface is the ability to interact with the three-dimensional 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 real-time. The 2D spatial plot displays increasing concentration with increasing brightness along a selected Cartesian axis.

## Results and discussion

### Michaelis-Menten reaction

The Michaelis-Menten equations are used to describe most enzymatic reactions. Its kinetics is given by the following equation:

*E* + *S* ↔ *ES* → *E* + *P*.

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

^{-20}litres, this amounts to approximately 30000 inert particles per step of 10%.

Simulation parameters

Volume (litres) | 10 |
---|---|

Number of | 3000 |

Number of | 1000 |

| 10 |

| 1 |

| 1 |

Simulation time (s) | 10 |

Timestep (s) | 10 |

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 well-mixed and un-crowded 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

Spatial simulators

GridCell | SmartCell | MesoRD | Cell++ | MCell | Smoldyn | ChemCell | |
---|---|---|---|---|---|---|---|

Molecule representation | Particle | Population | Population | Large particles and populations of small particles | Particle | Particle | Particle |

Stochastic | Yes | Yes | Yes | Large particles only | Yes | Yes | Yes |

Space | Discretized | Discretized | Discretized | Continuous (large particles) and discretized (small particles) | Continuous | Continuous | Continuous |

Particle-collision support | Yes | No | No | No | No | No | No |

Diffusion support | Yes | Yes | Yes | Yes | Yes | Yes | Yes |

SBML support | Yes | Yes | Yes | No | No | No | No |

Web availability | Yes | Yes | Yes | Yes | Yes | Yes | No |

### 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 large-scale 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 user-interface 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 well-mixed assumption and dimensionless particles can induce a significantly different response in a biological system. The discrete 3D grid and the nearest-neighbour interactions remove the need to do any distance calculation, particle search and floating-point 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 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.4-Xerces; http://sbml.org/software/libsbml/) and the OpenGL utility toolkit (GLUT 3.7.6; http://www.xmission.com/~nate/glut.html). EBI's Biomodels database: http://www.ebi.ac.uk/biomodels/. SBMLeditor: http://www.ebi.ac.uk/compneur-srv/SBMLeditor.html. JDesigner: http://sbw.kgi.edu/software/jdesigner.htm. Systems Biology Markup Language (SBML): http://sbml.org.

## Declarations

### Acknowledgements

The authors gratefully acknowledge funding provided by the Natural Sciences and Engineering Research Council of Canada.

## Authors’ Affiliations

## References

- Lemerle C, Ventura BD, Serrano L: Space as the final frontier in stochastic simulations of biological systems. FEBS Letters. 2005, 578: 1789-1794. 10.1016/j.febslet.2005.02.009.View ArticleGoogle Scholar
- Jorgensen K: Metabolon formation and metabolic channeling in the biosynthesis of plant natural products. Curr Opin Plant Biol. 2005, 8: 280-291. 10.1016/j.pbi.2005.03.014View ArticlePubMedGoogle Scholar
- Kholodenko BN: Four-dimensional organization of protein kinase signaling cascades: the roles of diffusion, endocytosis and molecular motors. Exp Biol. 2003, 206: 2073-2082. 10.1242/jeb.00298.View ArticleGoogle Scholar
- Saxton M: Anomalous Subdiffusion in Fluorescence Photobleaching Recovery: A Monte Carlo Study. Biophysical Journal. 2001, 81: 2226-2240.PubMed CentralView ArticlePubMedGoogle Scholar
- Turner T, Schnell S, Burrage K: Stochastic approaches for modelling
*in vivo*reactions. Computational Biology and Chemistry. 2004, 28: 165-178. 10.1016/j.compbiolchem.2004.05.001.View ArticlePubMedGoogle Scholar - Schnell S, Turner T: Reaction kinetics in intracellular environments with macromolecular crowding: simulations and rate laws. Progress in Biophysics and Molecular Biology. 2004, 85: 235-260. 10.1016/j.pbiomolbio.2004.01.012.View ArticlePubMedGoogle Scholar
- 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). 2008Google Scholar
- Plimpton S, Slepoy A: ChemCell: A Particle-Based Model of Protein Chemistry and Diffusion in Microbial Cells. Sandia Technical Report SAND2003-4509. 2003Google Scholar
- Andrews SS, Bray D: Stochastic simulation of chemical reactionswith spatial resolution and single molecule detail. Phys Biol. 2004, 1: 137-151. 10.1088/1478-3967/1/3/001View ArticlePubMedGoogle Scholar
- Lipkow K, Andrews SS, Bray D: Simulated diffusion of phosphorylated CheY through the cytoplasm of
*Escherichia coli*. J Bact. 2005, 187: 45-53. 10.1128/JB.187.1.45-53.2005PubMed CentralView ArticlePubMedGoogle Scholar - Zimmerman S, Minton A: Macromolecular crowding: biochemical, biophysical, and physiological consequences. Annu Rev Biophys Biomol Struct. 1993, 22: 27-65. 10.1146/annurev.bb.22.060193.000331View ArticlePubMedGoogle Scholar
- Gillespie DT: A general method for numerically simulating the stochastic time evolution of couple chemical reactions. Journal of Computational Physics. 1976, 22: 403-434. 10.1016/0021-9991(76)90041-3.View ArticleGoogle Scholar
- Boulianne L, Dumontier M, Gross WJ: A Stochastic Particle-Based Biological System Simulator. Proceedings of the Summer Computer Simulation Conference, San Diego, California (USA). 2007, 794-801.Google Scholar
- Ander M: 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. 10.1049/sb:20045017.View ArticleGoogle Scholar
- Hattne J, Fange D, Elf J: Stochastic reaction-diffusion simulation with MesoRD. Bioinfomatics. 2005, 21: 2923-2924. 10.1093/bioinformatics/bti431.View ArticleGoogle Scholar
- Sanford C, Yip MLK, White C, Parkinson J: Cell++simulating biochemical pathways. Bioinfomatics. 2006, 22: 2918-2925. 10.1093/bioinformatics/btl497.View ArticleGoogle Scholar
- Stiles JR, Bartol TM: Monte Carlo methods for simulating realistic synaptic microphysiology using MCell. Computational Neuroscience: Realistic Modeling for Experimentalists. Edited by: DeSchutter E. 2001, 87-127. Boca Raton: CRC PressGoogle Scholar

## Copyright

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.