Smoldyn

A spatial stochastic simulator for chemical reaction networks

Written and maintained by Steve Andrews

Introduction

Smoldyn is a computer program for simulating reactions and diffusion at a mesoscopic size scale. This scale lies between the more traditional microscopic size scale, in which all or most atoms are considered, and the macroscopic size scale, in which molecules are treated with concentrations rather than as individuals. Smoldyn is based on and named for Smoluchowski dynamics, which is a very simple description of chemical reaction dynamics. In Smoldyn, every molecule of interest is represented as an individual point, while those that are not of interest (water, non-reactive molecules, etc.) are not represented. Simulated space is continuous, meaning that space is neither a lattice, nor binned into voxels. The use of continuous space provides good spatial resolution and avoids artifacts that can arise from lattice geometries. Space can be 1, 2, or 3-dimensional.

Some of the algorithms that Smoldyn executes are described below.

Diffusion

Smoldyn does not consider inter-molecular forces, so each molecule diffuses independently of each other. If diffusion is isotropic, this is simulated by displacing each molecule, on each axis, by an amount that is chosen from a Gaussian probability density that has a standard deviation of s = (2D∆t)^1/2; s is called the rms step length, D is the diffusion coefficient, and ∆t is the simulation time step. Anisotropic diffusion is simulated by choosing a Gaussian distributed random displacement vector, exactly as in the isotropic case, and then multiplying it by the appropriate anisotropy matrix. Anistropic diffusion occurs in liquid crystals, aligned polymer media, and other ordered systems.

The pictures below show isotropic and anistropic diffusion, respectively, from Smoldyn simulations. In the former, red, green, and blue molecules were all started at the system center and allowed to diffuse outwards. Red molecules were assigned a high diffusion coefficient, green were assigned an intermediate value, and blue a very small value. For the latter picture, red and green molecules were again started at the system center and allowed to diffuse outwards, but this time both were assigned anisotropic diffusion matrices.

Surfaces

Surfaces in Smoldyn are built of panels, where each panel is a rectangle, triangle, sphere, cylinder, hemisphere, or a 1- or 2-dimensional analog of these shapes. Triangular panels are particularly useful for constructing complex mesh surfaces, while the others can be more useful for designing surfaces manually. Included in the Smoldyn download package is a small utility program called wrl2smol which converts wrl type triangle mesh data to Smoldyn format input.

When a molecule collides with a surface, it can be reflected off of the surface, absorbed by the surface, or transmitted through the surface. Reflection is performed with ballistic-type reflection, despite the fact that molecules are assumed to move exclusively through Brownian motion, because it yields the correct molecule position probability distribution. For absorption, Smoldyn used to absorb a molecule into an absorbing surface if the molecule either crossed the surface during a time step, or it was determined to have both crossed and recrossed the surface during the time step. It was found that the latter determination was difficult to work with, computationally intensive, and made little difference to results, so it is no longer implemented. Following are some images of molecules that are diffusing among surfaces in 1, 2, and 3 dimensions, and in a complex triangulated mesh (which was generously provided by Dylan Morris and made Smoldyn-readable with wrl2smol).

Reactions

Smoldyn can simulate three types of fundamental reactions: zeroth-order, unimolecular, and bimolecular. In zeroth-order reactions, product molecules are produced spontaneously, which is unphysical, but can be useful for simulating chemical inputs to the system. They are simulated by adding a Poisson-distributed random number of product molecules to the system at each time step. For unimolecular reactions, the probability that a molecule reacts is constant over time, which means that the probability that a reactant is still in existance after some time decreases exponentially. To simulate this, a reaction is considered to have happened during a time step if a uniformly distributed random number is greater than exp(-k₁∆t), where k₁ is the first-order reaction rate constant.

A particular strong point of Smoldyn is its ability to accurately and efficiently simulate bimolecular reactions. Following is a typical sequence of events.

Reaction diagram

A) Red and green molecules are diffusing.
B) Through diffusion, the reactants approach each other to within a binding radius, σ_b, which causes a reaction to occur. The size of the binding radius depends on the diffusion coefficients, the reaction rate constant, and the simulation time step. This radius is almost always smaller than the sum of the physical radii of the two reactants.
C) The reaction product is a blue molecule, which is placed at the reaction location.
D) The blue molecule diffuses.
E) The blue molecule dissociates into a red and a green molecule which are initially separated by the unbinding radius, σ_u. This radius is typically slightly larger than the sum of the physical radii of the two reactants.
F) The red and green product molecules of the dissociation reaction diffuse away from each other. Alternatively, they could diffuse back together and react again, which is called a geminate recombination.

The binding and unbinding radii are chosen by Smoldyn so that the simulated steady-state reaction rates and geminate recombination probabilities will equal values that the user requests.

The following images show simulations of the Lotka-Volterra reaction system in a thin 3-D system and in a crowded 2-D system, respectively. Both systems have periodic boundaries. Using R for red and G for green, the reactions are simply: R → 2 R, R + G → 2 G, and G → nothing.

lotka-volterra crowded lotka-volterra

Papers that used or discuss Smoldyn:

Key: •• = research paper that used Smoldyn, • = Smoldyn referenced, but not used.

2004
•• Andrews, Steven S. and Dennis Bray. "Stochastic simulation of chemical reactions with spatial resolution and single molecule detail." Phys. Biol. 1:137-151, 2004. (Description and analysis of algorithms).
• Ander, M., P. Beltrao, B. Di Ventura, J. Ferkinghoff-Borg, M. Foglierini, A. Kaplan, C. Lemerle, I. Tomás-Oliveira, and L. Serrano. "SmartCell, a framework to simulate cellular processes that combines stochastic approximation with diffusion and localisation: analysis of simple networks" Syst. Biol., 1:129-138, 2004. (Alternate software for spatial stochastic simulation.)
• Rivas, Germán, Frank Ferrone, and Judith Herzfeld. "Life in a crowded world." EMBO Reports 5:23-27, 2004. (Meeting report.)

2005
•• Andrews, Steven S. “Serial rebinding of ligands to clustered receptors as exemplified by bacterial chemotaxis.” Phys. Biol. 2:111-122, 2005. (Different code, but identical algorithm.)
•• Lipkow, Karen, Steven S. Andrews, and Dennis Bray. “Simulated diffusion of phosphorylated CheY through the cytoplasm of Escherichia coli.” J. Bact., 187:45-53, 2005. (Used Smoldyn.)
• Hazelbauer, Gerald L. “Myriad molecules in motion: Simulated diffusion as a new tool to study molecular movement and interaction in a living cell.” J. Bact., 187:23-25, 2005. (Commentary.)
• Takahashi, Kouichi, Satya Nanda Vel Arjunan, Masaru Tomita. “Space in systems biology of signaling pathways – towards intracellular molecular crowding in silico” FEBS Lett., 579:1783-1788, 2005. (Review.)
• Plimpton, Steven J. and Alex Slepoy. “Microbial cell modeling via reacting diffusive particles” J. Phys.: Conf. Ser., 16:305-309, 2005. (About ChemCell.)
•• DePristo, Mark A., Lynne Chang, Karen Lipkow, Ronald D. Vale, and Shahid Khan. “FRAP kinetics of the Escherichia coli chemoreceptor cluster analyzed by simulated diffusion” Submitted, 2005. (Used Smoldyn.)
• Le Novère, Nicolas and Dominic Tolle. “Particle-based stochastic simulations” Proceedings of the 4^th Workshop on Computation of Biochemical Pathways and Genetic Networks. Logos, Berlin, pp. 41-45. (Review.)
• Lemerle, Caroline, Barbara Di Ventura, and Luis Serrano. "Space as the final frontier in stochastic simulations of biological systems." FEBS Lett. 579:1789-1794, 2005. (Review.)

2006
•• Lipkow, Karen. "Changing cellular location of CheZ predicted by molecular simulations." PLoS Comp. Biol. 2:0301-0310, 2006. (Used Smoldyn).
•• Andrews, Steven S. and Adam P. Arkin. "Simulating cell biology" Curr. Biol. 16:R523-R527, 2006 (Used Smoldyn).
• Tolle, Dominic P. and Nicolas Le Novère. "Particle-based stochastic simulation in systems biology" Curr. Bioinf. 1:315-320, 2006. (Review.)
• Ridgway, Douglas, Gordon Broderick, and Michael J. Ellison. "Accomodating space, time and randomness in network simulation." Curr. Opinion Biotechnol. 17:493-498, 2006. (Review.)
• Tournier, Alexander L., Paul W. Fitzjohn, and Paul A. Bates. "Probability-based model of protein-protein interactions on biological timescales." Algorithms Molec. Biol. 1:25, 2006. (Investigated an alternative reaction algorithm to the one that Smoldyn uses.)
• Fange, David and Johan Elf. "Noise-induced Min phenotypes in E. coli." PLoS Comput. Biol. 2:637-648, 2006. (Used Meso-RD simulation program.)

2007
• Grima, Ramon and Santiago Schnell. "A mesoscopic simulation approach for modeling intracellular reactions." J. Stat. Phys. 128:139-164, 2007. (Developed an alternative algorithm.)
• Bray, Dennis, Matthew D. Levin, and Karen Lipkow. "The chemotactic behavior of computer-based surrogate bacteria." Curr. Biol. 17:12-19, 2007. (Used a non-spatial simulation method.)

Smoldyn download

Platform: Compiles fairly easily for Macintosh or Linux (instructions are in the user manual and Makefiles are provided). An executable application is included in the download package for Windows. Because Smoldyn uses the OpenGL library for graphics and the libtiff library for saving tiff files, these need to be supported by your computer. Mac users should already have OpenGL installed but will need libtiff, Linux users may need both, and Windows users don't need either but may need several .dll files instead. These upgrade instructions are in the manual. Here's the libtiff link: libtiff.

Current version: 1.85.

Documentation: The manual is in two parts. Part I (pdf) is for Smoldyn users and Part II (in the download package) is for Smoldyn programmers. Some algorithms are presented in a research paper that is available on my publications page.

Download: A downloadable Smoldyn package with the source code, documentation, makefiles, test files, binaries, and the wrl2smol utility is available as a .tgz format for Mac and Linux users, Smoldyn.tgz, or as a .zip format for Windows users, Smoldyn.zip. There is a sourceforge website for Smoldyn here but it's not active.

Old versions: 1.72 1.73 1.74 1.75 1.76 1.77 1.78 1.79 1.80 1.81 1.82 1.83 1.84

Collaborators: The Computational Research Laboratories, located in Pune, India, are adding parallel processing capability to Smoldyn. Also, Upi Bhalla is adding the Smoldyn simulation core into the MOOSE program.

History:
1.5 (7/2003) - First public release.
1.51 (9/5/03) - Unimolecular reactions can now have much lower reaction rates, a memory error upon termination was fixed, parameter checking was improved, several minor bugs were fixed, and two new interpreter commands were added.
1.52 (10/24/03) - A significant bug was fixed in the equation for absorbing walls (thanks to Dan Gillespie for catching it), it is now possible to pause simulations while using graphics, and a few minor bugs were fixed.
1.53 (2/9/04) - Cleaned up code for the runtime command interpreter, it is now possible to save file stacks, minor display and user interface improvements, and several new commands for text output and excluded volumes.
1.55 (8/20/04) - Minor graphics changes, some improvements in the design of the command queue for the runtime command interpreter, and it is now possible for commands to return error messages.
1.56 (1/14/05) - Improved graphics and a few other minor changes.
1.57 (2/17/05) - A new runtime command.
1.58 (7/22/05) - A new runtime command, a randomized Brownian dynamics lookup table to remove artifacts that can arise from random number generator problems, and a more versitile mol statement.
1.59 (?) - Displays the random number seed.
1.70 (5/17/06) - Reasonably good surface support, where surfaces can be reflective, absorbing or transparent, and they are comprised of panels that can be rectangular, triangular, or spherical.
1.71 (12/8/06) - Better surface support with cylinder and hemisphere panels, and improved reflective surfaces.
1.72 (2/26/07) - Surfaces are finally totally functional for both diffusion and reactions (or so I thought for a while).
1.73 (9/25/07) - Minor surface improvements, added exponential time-stepping option for run-time commands, and improved multi-platform compiling.
1.74 (10/22/07) - Fixed yet more bugs in surfaces, added disk-shaped surface panels, and made some changes in surface input statement formats. Only the code got updated for this release, so the documentation is slightly out of date.
1.75 (11/6/07) - Added support of surface-bound molecules, which still need some debugging. Also, changed several surface input statement formats. The documentation is current, but incomplete.
1.76 (11/7/07) - Main source file was split to two files and assorted function rearranging. More command-line options are now available.
1.77 (11/18/07) - Fixed a few minor bugs and modified several functions to make Smoldyn easier to parallelize.
1.78 (11/30/07) - Fixed lots of bugs with surface-bound molecules and commands, and improved input formats of several statements.
1.79 (12/6/07) - Fixed more bugs with commands. Also added basic compartments.
1.80 (12/22/07) - 'Q' now quits the program, time steps can be changed mid-simulation, and basic porting between Smoldyn and MOOSE was added.
1.81 (1/23/08) - Surface-bound molecules now work, lots of minor bugs were fixed, multiple molecule lists were implemented, porting works well, and documentation was improved.
1.82 (2/28/08) - Fixed a bug in which reaction products ended up across surfaces. Also, reactions were completely overhauled and can now be entered with a vastly improved format. Allosteric reactions were added.
1.83 (3/14/08) - Improved unimolecular reaction rates for multiple reaction channels, improved allosteric reactions, and added molecular drift.
1.84 (4/11/08) - First Windows release. Also, added reactions that are only active in specified compartments, added several commands, and fixed the change time step command.
1.85 (6/3/08) - The Mersenne Twister is now used for random numbers, added define statements to parser, renamed allosteric reactions to conformational spread reactions, and fixed lots of minor bugs and one significant one (reaction products ended up across surfaces).

Known bugs: No significant ones; a few very minor ones are listed at end of part II of the documentation. If you have a configuration file that causes crashes or incorrect behavior, please send it to me so I can find the bugs. My e-mail address is steven.s.andrews@gmail.com.