Page contents

Citations to Smoldyn


Smoldyn citations

See Google Scholar for all papers that mention Smoldyn. PubMed is a poor way to find publications that reference Smoldyn because PubMed only searches article titles and abstracts, and authors typically don't mention Smoldyn in those places.

Citations to all spatial cell biology simulators over time

The figure shows citation counts to most spatial cell biology simulators using data from Google Scholar. The red line is for Smoldyn, blue lines for other particle-based simulators, and grey lines for lattice-based simulators. The green line is for Virtual Cell which performs both PDE simulations and particle-based simulations, of which the latter use Smoldyn.

Published research that used Smoldyn simulations

Hybrid method

Discrete-continuous reaction-diffusion model with mobile point-like sources and sinks

Kondrat, S., O. Zimmermann, W. Wiechert, and E. von Lieres, Eur. Phys. J. E 39:11, 2016

In many soft and biological physics applications, there are multiple distinct time and length scales. For example, enzymes are large, move slowly, and are few, while metabolites are small, fast, and abundant. The authors addressed this with a new hybrid simulation method, which is stochastic-deterministic and discrete-continuous. They demonstrate it by modelling enzyme-catalysed reactions with discrete enzymes and continuous metabolites (see figure). They validated their method by comparing simulation results from their new method against those from Smoldyn.

Excluded volume effects in on- and off-lattice reaction-diffusion models

Meinecke, L. and M. Eriksson, ArXiv 2016

Simulations with off-lattice Brownian dynamics (e.g. Smoldyn) are computationally expensive in crowded environments. This paper investigates the extent to which on-lattice simulations can simulate reactions and diffusion in the presence of crowders. The authors show that diffusion is slowed down in the off-lattice model since randomly distributed obstacles effectively exclude more volume than those ordered on an artificial grid. Crowded reaction rates can be both increased and decreased by the grid structure. Grid artifacts increase with increasing crowder density. The authors conclude that the computationally more efficient on-lattice simulations are accurate only for low crowder densities.

Multiscale modeling of diffusion in a crowded environment

Meinecke, L., ArXiv 2016

This paper presents a multiscale approach to modeling diffusion and reactions rates in a crowded environment. The method combines jumping according to local first exit times and jumping on a coarser Cartesian grid. Excluded volume is modeled by a diffusion equation with space-dependent diffusion coefficients. Simulations showed that crowding molecule shape and size play a crucial role in the effects on diffusive motion. They also showed that molecular crowding can enhance or inhibit chemical reactions depending on local obstacle density fluctuations.
Crowding reactions

Simulating macromolecular crowding with particle and lattice-based methods

Andrews, S.S., S.N.V. Arjunan, G. Balbo, A.T. Bittig, J. Feret, K. Kaizu, and F. Liu, In D. Gilbert, M. Heiner, K. Takahashi, and A.M. Uhrmacher, eds. Multiscale Spatial Computational Systems Biology 170-187, 2015

The authors investigated the effect of macromolecular crowding on reaction rates using 5 different simulators, both to develop a better understanding of the topic and to compare the simulators. The eGFRD algorithm was presumably very accurate but ran very slowly, Smoldyn and NL-space worked well although Smoldyn appeared to produce more accurate results, Spatiocyte was very fast but results had lattice artifacts, and KaSim worked, which was remarkable because it is typically a non-spatial method. The Smoldyn results (figure) agreed with qualitative expectations but showed that the theory that the authors investigated was incorrect.
Diffusion MRI

Sensitivity Analysis of Diffusion Tensor MRI in Simulated Rat Myocardium

Bates, J., I. Teh, P. Kohl, J.E. Schneider, and V. Grau, Lecture Notes In Computer Science 9126:120-128, 2015

Diffusion MRI (magnetic resonance imaging) is a non-invasive experimental method for visualizing diffusion in tissues, in this case in rat heart tissue. Here, the authors modeled this method using Smoldyn and found close correspondance between model and experimental results. Next, they used the simulations to test the method sensitivity. They found that the diffusivity has the greatest effect and the cross-sectional area and aspect ratio of cells are important, but the cell length and volume fraction of cells had no marked effect.
Bax model

Bax monomers form dimer units in the membrane that further self-assemble into multiple oligomeric species

Subburaj, Y., K. Cosentino, M. Axmann, E. Pedrueza-Villalmanzo, E. Hermann, S. Bleicken, J. Spatz, and A.J. García-Sáez, Nature Comm. 6:8042, 2015

Bax is a key regulator of apoptosis, mediating cytochrome c release to the cytosol via oligomerization in the outer mitochondrial membrane. These authors investigated the molecular mechanism of Bax assembly and regulation by other Bcl-2 members. They found that Bax binds to the membrane in a monomeric state and then rapidly self-assembles. They also show that other proteins, cBid and Bcl-xL, help drive Bax activity. Based on these experimental results, they developed a theoretical model, using Smoldyn, which presents a new mechanism for the molecular pathway of Bax assembly to form the apoptotic pore.
T cell motility

Linking morphodynamics and directional persistence of T lymphocyte migration

Xiaji Liu, Erik S. Welf, and Jason M. Haugh, J. R. Soc. Interface 12:20141412, 2015

These researchers experimentally investigated T lymphocyte motility, which is similar to ameboid motility. They found that cells typically turn by first creating a bifurcation of the lamellipodium at the cell's leading edge and then following one of the new projections. They proposed a model of the major interactions in this process. Simulation with Smoldyn (shown in the figure) agreed with their experimental results. The figure shows the cell rear in beige and bifurcated actin regions in green.

Bistability: Requirements on Cell-Volume, Protein Diffusion, and Thermodynamics

Robert G. Endres, PLoS ONE 10:e0121681, 2015

Many cell systems, including cell cycle ones and some signaling ones, are bistable. This has been studied most often using deterministic and/or non-spatial models. This paper used Smoldyn models to investigate the more physiologically relevent case, which is for stochastic chemical reactions in the confined 3-dimensional volume of a cell, considering both bacteria and eukaryotic cells. The author finds that bistability is fragile in these conditions, often requiring finely tuned parameters, small volumes, and fast diffusion coefficients. Switching can occur upon cell growth.
DNA target finding

An Integrated Model of Transcription Factor Diffusion Shows the Importance of Intersegmental Transfer and Quaternary Protein Structure for Target Site Finding

Hugo G. Schmidt, Sven Sewitz, Steven S. Andrews, and Karen Lipkow, PLoS ONE 9:e108575, 2014

Transcription factors and other DNA binding proteins find their DNA binding sites faster than simple 3D diffusion allows. These authors explored this acceleration with Smoldyn models. As in prior work, they found that non-specific DNA binding followed by 1D sliding reduces finding times. They also found that intersegmental transfer, in which a transcription factor that is non-specifically bound to DNA simultaneously binds a separate DNA loop and then releases the first binding site, is also effective. Finally, they found that DNA binding proteins are enriched in dimers and tetramers, perhaps because this favors intersegmental transfer.
chromosome domains

Nucleolar tethering mediates pairing between the IgH and Myc loci

Daniel E. Strongin, Mark Groudine, Joan C. Ritland Politz, Nucleus 5:474-481, 2014

The genome is spatially organized within the eukaryotic nucleus. One aspect of this is that gene loci on different chromosomes can preferentially colocalize. Using mouse strains that have different gene arragements on their chromosomes, these authors investigated the driving force behind the colocalization of the IgH and Myc loci. They found that it arose when both loci were on chromosomes that had nucleolar organizer regions (NORs). These are sites of ribosomal DNA repeat sequences, which nucleate nucleoli. Together with simulations, these results implied that chromosome tethering to nucleoli can help colocalize genes. The figure shows one chromosome in yellow, another in blue, and nucleoli in red, in a strain where colocalization occurs.
receptor stochasticity

Mosaic physiology from developmental noise: within-organism physiological diversity as an alternative to phenotypic plasticity and phenotypic flexibility

H. Arthur Woods, J. Experiment. Biol. 217:35-45, 2014

In this paper, the author proposes that stochastic variation during organismal development causes physiological diversity within a single individual. He calls this mosaic physiology. This article reviews known mechanisms by which stochastic effects arise in, and are controlled by, biological systems. It then goes on to show how then can give rise to mosaic physiology. This provides a set of diversified phenotypes within single organisms, which may help the organism to cope with novel environmental challenges. The figure shows stochastic effects in receptor-ligand binding, a simple example of biological stochasticity.
L. lactis diffusion

Impact of osmotic stress on protein diffusion in Lactococcus lactis

Jacek T. Mika, Paul E. Schavemaker, Victor Krasnikov, Bert Poolman, Mol. Microbiol. 94:857-870, 2014

Proteins diffuse many times slower inside cells than they do in aqueous solutions due to macromolecular crowding, non-specific binding, and hydrodynamic effects. This has been investigated most thoroughly in E. coli, a Gram-negative bacterium. The authors investigated diffusion in L. lactis, a Gram-positive bacterium here. They found similar diffusion coefficients and similar cell-to-cell variation in unstressed cells, but a smaller increase upon osmotic challenge (as expected due to a higher turgor pressure). They measured intracellular diffusion with FRAP methods, which they validated with Smoldyn simulations (figure).
Lambda switch

Stochastic cellular fate decision making by multiple infecting lambda phage

Matthew L. Robb and Vahid Shahrezaei, PLOS One 9:e103636, 2014

Bacteriophage lambda is a classic system for studying cellular decision making. Upon infection, lambda can adopt the lytic state in which it reproduces rapidly in the bacterium causing cell death and the release of virus, or it can adopt the lysogenic state in which it incorporates itself into the cell genome; later, it may stochastically return to the lytic state. This decision is made through a genetic switch, shown in the figure. These authors extended prior work by investigating the role of the cell volume and bacterial growth rate on the decision. They also investigated spatial effects, which arise primarily from the slow diffusion of mRNA across a bacterium (about 10 minutes), finding that spatial effects were minimal.
Diffusion barrier

The long and viscous road: uncovering nuclear diffusion barriers in closed mitosis

Eder Zavala and Tatiana T. Marquez-Lago, PLoS Comp. Biol. 10:e1003725, 2014

Yeast cells retain their nuclear membranes during cell division, in a processs called closed mitosis. Membrane-bound proteins segregate aysmmetrically in the process, with some getting localized in the mother cell and others in the bud (dots in the figure). These authors explored mechanisms by which yeast cells might prevent protein diffusion across the division plane, and hence maintain the localization. They found that a combination of protein rings and sphingolipid domains is necessary during early anaphase, but that sphingolipid domains alone are adequate during late anaphase (figure), due to the elongated nuclear neck.
Min system

Oscillations of Min-proteins in micropatterned environments: a three-dimensional particle-based stochastic simulation approach

Max Hoffmann and Ulrich S. Schwarz, Soft Matter 10:2388, 2014

One way in which E. coli bacteria find their mid-planes, so that cell division yields two equal size daughter cells, is with spatiotemporal oscillations of the Min proteins. These proteins oscillate from pole to pole, with minimal occupancy of the mid-plane. The simplicity and remarkable dynamics of this system has made it popular for spatial stochastic simulations. These authors investigated Min system operation in artificial micropatterned environments and in mutant filamentous cells, such as the one shown in the figure. This work highlights the robustness and variability of Min system oscillations, puts limits on the effect of putative division sites, and provides a computational framework for future studies.

Input-output relations in biological systems: measurement, information and the Hill equation

Steven A. Frank, Biology Direct 8:31, 2013

This paper discusses relations among Michaelis-Menten kinetics, the Hill equation, and biological information processing. A central question that the paper probes regards why Michaelis-Menten kinetics leads to linear input-output relations for low signal levels, but the Hill equation exhibits logarithmic sensitivity at these levels. This paper has numerous misconceptions and non-standard terminology but nevertheless presents some intriguing points. A particularly interesting result is that dose-response curve for the Michaelis-Menten reaction, which normally exhibits Hill equation behavior with Hill coefficient of 1, becomes ultrasensitive when reactions are substantially diffusion influenced, as shown in the figure.
TEM image

Dynamic transition states of ErbB1 phosphorylation predicted by spatial stochastic modeling

Meghan McCabe Pryor, Shalini T. Low-Nam, Ádám M. Halász, Diane S. Lidke, Bridget S. Wilson, and Jeremy S. Edwards, Biophys. J. 105:1533-1543, 2013

ErbB1 (epidermal growth factor receptor) is an important receptor for growth and development, and its overexpression can cause cancer. These authors built on their prior erbB1 single molecule tracking experiments to build a spatial stochastic model for erbB1 diffusion, dimerization, kinase activation, and phosphorylation. The figure shows erbB1 dimerization. The model yields new insight into the activation states of individual erbB1 monomers. For the most part, the authors used their own reimplementation of Smoldyn's algorithms, although they also used the Smoldyn for part of the research as well.
pore simulation

Membrane and actin reorganization in electropulse-induced cell fusion

Günther Gerisch, Mary Ecke, Ralph Neujahr, Jana Prassler, Andreas Stengl, Max Hoffmann, Ulrich S. Schwarz and Eberhard Neumann, J. Cell Science 127:4507-4517, 2013

Electric pulses induce Dictyostelium discoideum cells to fuse. The authors combined electron microscopy, fluorescence microscopy, and simulation to study the fusion pores and actin localization that arise in the membrane during cell fusion. They found that the plasma membranes of the contiguous cells become tangles of highly bent and interdigitated membranes. By imaging GFP diffusion from one cell to its neighbor, and then modeling this diffusion with Smoldyn simulations as shown in the figure, they found that membranes persist in a fusogenic state for up to 24 seconds before pores of about 3 nm are formed.

Anomalous diffusion and multifractional Brownian motion: simulating molecular crowding and physical obstacles in systems biology

T.T. Marquez-Lago, A. Leier, and K. Burrage, IET Syst. Biol. 6:134-142, 2012

The study compares both fractional Brownian motion and continuous time random walks and highlights how well they can represent different types of spatial crowding and physical obstacles. Although diffusion around immovable obstacles could be reasonably characterised by a single Hurst exponent, the authors found that diffusion in a crowded environment seemed to exhibit multifractional properties in the form of a different short- and long-time behaviors. The figure shows a 2-dimensional crowded environment, developed with the Smoldyn utility program SmolCrowd.
yeast cells

Nuclear envelope morphology constrains diffusion and promotes asymmetric protein segregation in closed mitosis

Barbara Boettcher, Tatiana T. Marquez-Lago, Mathias Bayer, Eric L. Weiss, and Yves Barral, J. Cell Biol. 197:921-937, 2012

The authors investigated the establishment and maintenance of asymmetric cell division in budding yeast. Here, unlike in mammalian cells, the nuclear envelope is maintained during cell division, and divides at about the same time as the cell membrane. Using photobleaching and Smoldyn simulations, they found that diffusion barriers compartmentalize the nuclear membranes, whereas protein diffusion was well explained by the dumbbell shape of the anaphase nucleus. The figure shows a dividing cell with the nucleoplasmic GFP in green and mCherry-Tub1 spindles in red.
mushroom spine

Spatiotemporal maps of CaMKII in dendritic spines

Khan, S., T.S. Reese, N. Rajpoot, and A. Shabbir, J. Comput. Neurosci. 33:123-139, 2012

The authors investigated the sequestration of calcium calmodulin dependent kinase (CaMKII) in neural dendritic spines, which is a key cellular mechanism for the formation and storage of memories. The figure shows confocal microscopy data of GFP-CamKII localization shortly the GFP was photo-activated at the asterisk; bright green regions, such as at the spine tip, indicate CamKII sequestration, while red regions indicate low CamKII concentrations. Smoldyn simulations showed that a major cause of sequestion is the high number of cytoskeletal binding sites at spine tips, rather than high binding affinities.
rate comparison

Correction factors for boundary diffusion in reaction-diffusion master equations

Leier, A. and T.T. Marquez-Lago, J. Chem. Phys. 135:134109, 2011

The two primary methods for simulating chemical reactions and diffusion are continuous-space single particle methods, like Smoldyn uses, and lattice-based reaction-diffusion master equation (RDME) methods. The latter are typically more efficient but can suffer from several artifacts. In this work, the authors presented an accurate RDME method for simulating reactive boundary conditions, which they validated using Smoldyn. The figure shows simulation results from their methods in black, analytical results in red, and Smoldyn results in gray.
low concentration high concentration

A density-dependent switch drives stochastic clustering and polarization of signaling molecules

Jilkine, A., S.B. Angenent, L.F. Wu, and S.J. Altschuler, PLoS Comput. Biol. 7:e1002271, 2011

Many cell types, including budding yeast, mammalian neutrophils, and amoeba, can spontaeously polarize in the absence of spatial cues. The authors propose that this polarization arises from positive feedback; it robustly maintains an off state at low concentrations of signaling molecules, and then switches to highly localized signaling clusters at higher concentrations. The figures show that in their model, increasing signaling molecule concentrations (right panel vs. left panel) leads to clusters of molecules that are in their active states (red).
neural synapse

Computational investigation of the changing patterns of subtype specific NMDA receptor activation during physiological glutametergic neurotransmission

Singh, P., A.J. Hockenberry, V. Tiruvadi, and D.F. Meaney, PLoS Comput. Biol. 7:e1002106, 2011

Neural NMDA receptors mediate many physiological functions, including the molecular basis for learning and memory. These receptors exist in various subtypes, which enables them to discriminate between different types of signals. Using Smoldyn models, the authors found that different receptors have different dynamic ranges, that specific subtypes dominate in long-term depression and long-term potentiation situations, and that the content of a specific subtype enhances response magnitude and fidelity during long-term potentiation. The yellow and red portions of the figure show pre- and post-synaptic regions, while dots represent receptors.
neural spine

Sequestration of CaMKII in dendritic spines in silico

Khan, S., Y. Zou, A. Amjad, A. Gardezi, C.L. Smith, C. Winters, T.S. Reese, J. Comput. Neurosci. 31:581-594, 2011

This work investigated the sequestration of calcium dependent kinase II (CamKII) in neural dendritic spines following synaptic stimulation. It looked at the effects of spine geometry on CamKII diffusion, binding of CamKII to the post-synaptic density, to the cytoskeleton, and to other CamKII proteins. Simulation results compared favorably with microscopy experiments. This work showed that self-aggregation of CamKII could provide a switch that amplifies CamKII sequestration and regulates its activity.

Ultrasensitivity in Multisite Phosphorylation of Membrane-Anchored Proteins

Dushek, O., P.A. van der Merwe, and V. Shahrezaei, Biophys. J. 100:1189-1197, 2011

The authors showed that multiple phosphorylation sites on a single protein can produce switch-like responses when reactions between proteins are limited by diffusion. This is the case for many reactions that occur on the plasma membrane. The figure shows that ultrasensitivity increases with as the number of protein phosphorylation sites increase. Smoldyn file.
NeuroRD comparison

The role of type 4 phosphodiesterases in generating microdomains of cAMP: Large scale stochastic simulations

Oliveira, R.F., A. Terrin, G. Di Benedetto, R.C. Cannon, W. Koh, M. Kim, M. Zaccolo, and K.T. Blackwell, PLoS ONE 5:e11725, 2010

This work investigates mechanisms that increase the local concentration of cyclic AMP (cAMP), which is essential for normal neural functioning, including for synaptic plasticity. To do so, the authors developed a new spatial simulator called NeuroRD, which is a compartment-based simulator. It combines the spatial Gillespie method with tau-leaping and a new diffusion algorithm. The authors validated NeuroRD by comparing results against mass action theory and against Smoldyn simulations. These comparisons are shown in the figure.
Bar1 model

Detailed Simulations of Cell Biology with Smoldyn 2.1

Andrews, S.S., N.J. Addy, R. Brent, and A.P. Arkin, PLoS Comp. Biol. 6:e1000705, 2010

The authors investigated mating pheromone signaling between yeast cells. A central "receiver" cell is covered with receptors that bind mating pheromone, which is secreted by surrounding "sender" cells (in the figure, unbound receptors are blue, bound receptors are red, and pheromone is green). This work showed that the central cell can better locate strong pheromone emitting cells when the central cell secretes a pheromone-degrading protease because the protease cloud sharpens the local pheromone gradient.
Simulated FRAP

Introducing simulated cellular architecture to the quantitative analysis of fluorescent microscopy

DePristo, M.A., L. Chang, R.D. Vale, S.M. Khan, and K. Lipkow, Prog. Biophys. Mol. Biol. 100:25-32, 2009

The authors used simulations to analyze and quantify experimental Fluorescence Recovery After Photobleaching (FRAP) data for proteins in the E. coli chemotaxis system. They quantified protein diffusion coefficients of 2 μm2/s and assessed turnover rates between cytoplasmic proteins and membrane-associated protein clusters. The figure shows simulated CheY-YFP fluorescence in a protein cluster as it is photobleached and then recovers from protein exchange. In general, simulations can be useful tools for the quantitative analysis of fluorescent microscopy experiments.
Protein gradients

Model for Protein Concentration Gradients in the Cytoplasm

Lipkow, K. and D.J. Odde, Cell. Mol. Bioeng. 1:84-92, 2008

Using simulations, the authors showed that several of the E. coli chemotaxis proteins (CheY and CheZ) likely have stable protein concentration gradients across the length of cells. This arises from protein activation that is localized to one cell pole and activation-dependent protein complexation that affects protein diffusion coefficients. These mechanisms, and similar ones, occur for many other proteins as well, so stable protein concentration gradients within cells are likely to be widespread.
Hair cell PMCA2

Rapid Turnover of Stereocilia Membrane Proteins: Evidence from the Trafficking and Mobility of Plasma Membrane Ca2+-ATPase 2

Grati, M., M.E. Schneider, K. Lipkow, E.E. Strehler, R.J. Wenthold, and B. Kachar, J. Neuroscience 26:6386-6395, 2006

This work combined experimental and modeling approaches to investigate membrane protein turnover in stereocilia (fine hairs in the inner ear that convert mechanical motion to electrical signals) to better understand stereocilia recovery from trauma, such as loud noises. The authors focused on the spatial distribution, mobility, and trafficking of the PMCA2 protein, which is both abundant and essential to stereocilia function. They found that PMCA2 exhibits rapid turnover, which supports other evidence that stereocilia undergo rapid continuous renewal. The figure shows PMCA2 (green dots) distribution on an unrolled stereocilia membrane over time.
Lotka Volterra scale Lotka Volterra dynamics

Simulating cell biology

Andrews, S.S. and A.P. Arkin, Current Biology 16:R523-R527, 2006

This is a primer on how to simulate cell biology systems. It presents the basics of ODE (ordinary differential equation) modeling, stochastic modeling, spatial modeling, and spatial-stochastic modeling. Of particular interest, it presents simulations of the same Lotka-Volterra predator-prey system that used each of the different methods. The results differ dramatically. Deterministic simulations show nearly sinusoidal oscillations, non-spatial stochastic simulations show oscillations with increasing amplitudes, and particle-tracking (Smoldyn; figure show here) simulations show a characteristic boom-and-bust pattern.
CheZ localization

Changing Cellular Location of CheZ Predicted by Molecular Simulations

Lipkow, K., PLoS Comp. Biol. 2:e39, 2006

In the E. coli chemotaxis system, the CheY messenger protein is phosphorylated at a receptor cluster near one cell pole and is dephosphorylated by the CheZ protein. The author used simulations to explore the signaling behaviors that arise with different CheZ localizations. She found that the model that best agreed with experiments included CheZ proteins that oligomerize and localize to the receptor cluster upon cellular stimulation. The black dots in the figure are CheZ oligomers (also, pink dots are CheY, red are CheYp, and green are non-oligomeric CheZ). This creates a negative feedback loop which improves system precision, robustness, and adaptation.
FliM occupancy

Simulated Diffusion of Phosphorylated CheY through the Cytoplasm of Escherichia coli

Lipkow, K., S.S. Andrews, and D. Bray, J. Bact. 187:45-53, 2005

The authors developed a spatial stochastic simulation of the central E. coli chemotaxis system. It includes the core proteins and reactions, a polar-localized receptor cluster, and several flagellar motors in the cell membrane. Increasing cytoplasmic crowding, due to ribosomes and other macromolecules, resulted in slower protein diffusion, steepened concentration gradients, and substantial differences in signal propagation delays to motors at different positions (figure).

Publications that describe Smoldyn software and algorithms

Andrews image

Rule-based modeling using wildcards

Andrews, S.S., BioRxiv

This article presents a new approach to rule-based modeling, which is a method for automatically generating new chemical species for modified and complexed molecules. This approach is based on wildcards that match to species names, much as wildcards can match to file names in computer operating systems. It is simple to use and remarkably powerful. This article demonstrates the approach through examples for: signaling systems, protein complexation, polymerization, nucleic acid sequence copying and mutation, the "SMILES" chemical notation, and others. It is implemented in Smoldyn for both the generate-first and on-the-fly expansion.
Andrews image Andrews image

Smoldyn: particle-based simulation with rule-based modeling, improved molecular interaction and a library interface

Andrews, S.S., Bioinformatics 33:710, 2017

This paper describes several additions to the Smoldyn software. Smoldyn now supports rule-based modeling with a convenient wildcard method, and also with the BioNetGen package with extensions for spatial simulation. New algorithms for simulating the diffusion of surface-bound molecules and molecules with excluded volume provide excellent accuracy (figure). In addition, Smoldyn supports single-molecule tracking simulations. Finally, the Smoldyn source code can be accessed through a C/C++ language library interface.
Seeliger image

Enabling surface dependent diffusion in spatial simulations using Smoldyn

Seeliger, C. and N. Le Novère, BMC Res. Notes 8:752, 2015

The authors added support for diffusion coefficients of surface-bound molecules to depend on the surface to which the molecule is bound. For example, in the figure, certain surface-bound molecules diffuse more slowly within the blue region of the surface than in other portions of the surface. This is useful for simulating lipid rafts or similar structures. Unfortunately, this addition is not included in the current Smoldyn release, but only in the authors' modification of Smoldyn 2.22, which is otherwise largely obsolete.
Hybrid Smoldyn

Multiscale reaction-diffusion simulations with Smoldyn

Robinson, M., S.S. Andrews, and R. Erban, Bioinformatics 31:2406-2408, 2015

The authors integrated a lattice-based simulation solver with Smoldyn's existing particle-based simulations. This is an adjacent-space type of hybrid method, meaning that different regions of space get simulated with different levels of detail. They used recently devloped methods (by Flegg and Erban) to address the boundary between two different simulation regions. This accelerates simulations several-fold. It is very easy to use.
Smoldyn example

Spatial and stochastic cellular modeling with the Smoldyn simulator

Andrews, S.S., In Bacterial Molecular Networks: Methods and Protocols, van Helden at al. (eds.) Methods in Molecular Biology 804:519-542, 2012

This paper describes the use of the Smoldyn simulation program. It focuses on the basic concepts, while leaving syntax details to the Smoldyn User's Manual. This figure was taken from an example simulation that presents the use of reactions, molecule-surface interactions, compartments, and other Smoldyn features. This paper presents an especially useful Notes section that gives suggestions on choosing simulation parameters (e.g. diffusion coefficient and reaction rates) and methods for improving computational performance.
GPU Smoldyn

Smoldyn on Graphics Processing Units: Massively Parallel Brownian Dynamics Simulation

Dematté, L., IEEE Trans. Comput. Biol. Bioinform. 9:655-667, 2012

This is a more thorough presentation of the work described below by the same author. In brief, he parallelized the core Smoldyn code to run on Graphics Processing Units (GPUs), thus accelerating it 130-fold.
GPU Smoldyn

Accelerating the Smoldyn spatial stochastic biochemical reaction network simulator using GPUs

Gladkov, D.V., S. Alberts, R.M. D'Souza, and S. Andrews, Proceedings of the 19th High Performance Computing Symposia, 2011

The authors accelerated Smoldyn simulations by about 200 fold by parallelizing the core code to run on Graphics Processing Units (GPUs). This work was performed independently of Dematté's work, listed below. The figure shows code acceleration for different Smoldyn algorithms. Diffusion has the best performace, while bimolecular reactions have the worst (but still excellent) performance.
GPU Smoldyn

Parallel particle-based reaction diffusion: A GPU implementation

Dematté, L., 2010 Ninth International Workshop on Parallel and Distributed Methods in Verification/2010 Second International Workshop on High Performance Computing 67-77, 2010

The author parallelized the core Smoldyn code to run on Graphics Processing Units (GPUs). This can accelerate simulations by factors of up to 130, thus enabling the efficient simulation of models that would be impractical to run on conventional computer hardware. The figure explains efficient memory organization for GPUs.
Smoldyn accuracy

Detailed Simulations of Cell Biology with Smoldyn 2.1

Andrews, S.S., N.J. Addy, R. Brent, and A.P. Arkin, PLoS Comp. Biol. 6:e1000705, 2010

This paper (which is also listed above) describes the Smoldyn software. It gives an overview of the software algorithms, features, and performance. It shows that Smoldyn is very accurate, is more than twice as fast as either MCell or ChemCell, and has a wide range of features. This figure shows the number of molecules adsorbed to a surface over time for several adsorption coefficients; close accuracy between simulation data (dots) and theoretical calculations (lines) show Smoldyn's accuracy. The inset shows that simulation noise matches theoretical calculations.

Accurate particle-based simulation of adsorption, desorption and partial transmission

Andrews, S.S., Phys. Biol. 6:046015, 2009

This paper describes methods for quantitatively simulating the adsorption of molecules to surfaces, the desorption of molecules from surfaces, and the partial transmission of molecules through surfaces. Central results of the paper include look-up tables that give a molecule's adsorption or transmission probability based on the surface's adsorption or transmission coefficient, the molecule's diffusion coefficient, the simulation time step, and the rate of the reverse desorption or transmission process. Smoldyn uses these algorithms for molecule-surface interactions.
Simulator comparison

Computational methods for diffusion-influenced biochemical reactions

Dobrzynski, M., J.V. Rodríguez, J.A. Kaandorp, and J.G. Blom, Bioinformatics 23:1969-1977, 2007

The authors compared the three major spatial stochastic simulation methods, which are particle-based methods (with Smoldyn), spatial Gillespie methods (with GMP), and Green's Function Reaction Dynamics (GFRD; with custom code). They found that GFRD is most accurate and that the relative accuracy between particle-based and spatial Gillepsie methods depends on the latter's spatial resolution; with fine resolution, spatial Gillespie methods can be better than particle-based methods. For systems with many molecules, GFRD is extremely slow, while Smoldyn takes about twice as long as GMP to yield results with comparable accuracy. The authors noted that particle-based methods are particularly flexible for creating realistic geometries.
Bimolecular reaction simulation

Stochastic simulation of chemical reactions with spatial resolution and single molecule detail

Andrews, S.S. and D. Bray, Phys. Biol. 1:137, 2004

The authors derived several of the central algorithms for particle-based simulation, all designed for high accuracy and good computational efficiency. Methods are presented for simulating diffusion and chemical reactions, where the reactions can be 0th, 1st, or 2nd order and they can be irreversible or reversible. This figure shows one of the paper's central results, which is a look-up table that enables simulated bimolecular chemical reaction rates (y-axis) to be determined from simulation parameters (x-axis) for both irreversible (bottom line) and reversible reactions (other lines).