Enabling surface dependent diffusion in spatial simulations using Smoldyn
 Christine Seeliger^{1}Email author and
 Nicolas Le Novère^{1, 2}Email author
Received: 9 February 2015
Accepted: 19 November 2015
Published: 8 December 2015
Abstract
Background
Spatial computer simulations are becoming more feasible and relevant for studies of signaling pathways due to technical advances in experimental techniques yielding better high resolution data. However, many common single particle simulation environments used in computational systems biology lack the functionality to easily implement spatially heterogeneous membrane environments.
Results
We introduce an extension to the single particle simulator Smoldyn that allows modeling of surfacedependent diffusion, without unnecessarily increasing molecular states or numbers, hence avoiding explosion of molecule and reaction definitions.
Conclusions
We demonstrate the usefulness of this approach studying AMPA receptor diffusion at the postsynaptic density and its spatial trapping without introducing hypothetical scaffold elements or membrane barriers.
Keywords
Findings
Computational modeling is increasingly used to study the dynamic behavior of complex signaling systems in time and space. Different modeling approaches, and with them simulation environments, exist to allow the development, implementation and analysis of signaling models. With improvements in experimental techniques, detailed spatial models become increasingly feasible (see [1] for review).
The composition of biological membranes is heterogeneous. This heterogeneity is a regulatory mechanism that controls and regulates signaling pathways via location and movement of their components and is reflected in the frequent observation of their subdiffusive behavior [2, 3].
An example of such a specialized environment is the postsynaptic density (PSD). The PSD is located at the postsynaptic side of neuronal synapses, receiving the biochemical signals originating from the presynaptic terminal during neurotransmission. The PSD is a highly organized structure visible in electron microscopy micrographs as an electron dense area in spine heads, the small synaptic protusions found on dendrites [4]. Techniques such as singleparticle tracking showed that molecules, including AMPA receptors, exhibit different macroscopic diffusive behavior when they are moving within the PSD or the extrasynaptic membrane [5–7]. The environment of the PSD acts as a trap for receptor molecules depending on synaptic activity, hence regulating synaptic strength [8, 9]. To study these subdiffusive effects in the context of signaling pathways, the membrane environments have to be explicitly represented in the geometry of spatial computer models. This integration should help to better understand the importance of spatially organized and regulated signaling pathways.
Currently, most available simulators [10, 11] require the introduction of an artificial boundary surface surrounding the desired surface area. In recent versions of Smoldyn, a molecule changes its type upon crossing this boundary. This allows a new type of behavior to be specified. However, this comes at the costs of defining the interactions with the boundary surface. This requires the definition of transition rates for molecules crossing the surface. In addition, reactions that take place on all types of surfaces have to be defined for all variants of molecules for each surface type. This results in an explosion of molecule and reaction definitions in the model. This aspect becomes more important as numbers of molecule and surface types increase.
We introduce an extension to the simulation environment Smoldyn that allows to model surfacedependent diffusion. This modified version allows the macroscopic modeling of different diffusive areas of biological membranes. Artificially introduced boundary surfaces intersecting the diffusive surface as described above are not necessary. The modified version Smoldyn is described in the following sections and validated against the original. Examples based on phosphoinositide phosphorylation and AMPAR diffusion in the dendritic spine illustrate its use.
Implementation
The main modifications to the original Smoldyn source code affect the structures that hold the surface and reaction data. Each surface now holds the information for its specific diffusion coefficients accessible by molecule id and state. Based on this, the functions that calculate reaction parameters such as binding and unbinding radii and diffusion steps can be calculated in a surface dependent manner using the original methods provided by Smoldyn. These modifications do not cause any changes in the control flow of the simulator.
In the following, the modified version of Smoldyn will be denoted SmoldynM to distinguish it from the original.
The modified files (smolmolec.c, smolreact.c, smolsurface.c) as well as the complete source code are available online and can be downloaded from GitHub.
Diffusion
First, diffusion was tested to show that the new version is working as expected by comparing both Smoldyn versions in scenarios that both simulators can handle. The testing geometry comprises a rectangular plane inside a 3dimensional box as shown in Fig. 1a. The plane is separated in two triangular panels, each defined as a surface in its own right. The random seeds were fixed to directly compare the diffusion between both Smoldyn versions. Single molecules were set to diffuse for a set amount of time (100 s). Their paths along the surfaces were compared and are identical between simulations using the same seed.
SmoldynM is expected to run slightly slower than the original version. For n molecular types and m surfaces, lookup tables increase by \(\text{ n } \cdot \text{ m }\). Additional operations are simple indexed lookups and should not increase runtimes too much. Running a diffusion simulation with 10,000 molecules over 300 s simulated time and the two triangular surfaces mentioned above, the realtime runtime increases from 617 s to 638 s when using SmoldynM (run on a standard laptop with an Intel i7 CPU and 8GB memory).
The mean square displacement (MSD) [12] of particles was calculated based on 20 simulations tracking one molecule. The MSD can be interpreted as the space explored by particles over time. The MSD of random diffusion is expected to be linear over time. If the molecules are allowed to cross between triangles defining a different diffusion coefficent, the MSD shows the anomalous behavior expected from such a system as shown in Fig. 1c for SmoldynM, indicating subdiffusion, a behavior that is observed for example in systems that exhibit molecular crowding [13].
Reactions
The same geometry as described before and illustrated in Fig. 1a was used to run a set of simulations to compare reactions in Smoldyn and SmoldynM. The modeled reaction is based on the phosphorylation of the phosphoinositide PIP2 by phosphatidylinositol3kinase (PI3K). The kinetic parameters where taken from [14]. Diffusion on the “fast” surface is again set up 100x faster than the slow one. Simulations are initiated with 3 kinases and 500 molecules of PIP2 for each triangle. PIP2 is not allowed to move between the fast and slow triangles in the first set of simulation runs and the boundaries of the triangles are reflective.
In this case, direct comparison between Smoldyn and SmoldynM is only possible if it is ensured that the same sequence of random numbers is used during simulation. This required SmoldynM simulations to be run separately for the “fast” and “slow” triangle, mirroring Smoldyn, although simulations using the two different diffusion coefficients for each triangle can run simultaneously using SmoldynM. However, the different number of events by having both triangles populated with molecules and their reactions changes the amount of numbers used for each timestep and would therefore alter the results.
Figure 1d shows the results of these simulations in purple. The development of PIP3 on the “fast” (dark purple) and “slow” (light purple) is identical for Smoldyn and SmoldynM for identical random number seeds. The number of PIP3 is larger on the “fast” surface due to increase chances of PIP2 encountering PI3K.
In addition, simulations were run in SmoldynM connecting the two surfaces, setting the two different diffusion coefficients and allowing PIP2 to cross between surfaces. The same random seed was used as before. The development of PIP3 is shown in blue for both surfaces on Fig. 1d. The trapping of PIP2 molecules on the “slow” diffusion surface causes an increase in the local concentration compared to the “fast” surface. Therefore, more PIP3 is produced on the slow than on the fast surface in this case compared to the other cases where crossing of molecules between surfaces was not possible.
Example: AMPAR trapping at the PSD
The increase in synaptic AMPAR numbers is one of the major processes of long term potentiation (LTP). Sources for new AMPAR at the synapse are exocytosis and lateral diffusion [15]. This section shows an application for surface dependent diffusion coefficients to study the latter (lateral diffusion) based on a model proposed by [9].
Different macroscopic diffusion environments can mimic scaffold binding
Binding partners and mechanisms of detainment of molecules in specific areas of the membrane are often unknown, only vaguely defined or very broad. This can render explicit modeling difficult. Molecular behavior could also be solely influenced by crowded environments. Different macroscopic diffusion environments can be used to abstract from this level of detail as shown for example in Fig. 1a.
Different diffusive environments are able to simulate receptor trapping at the synapse
Tolle et al. used explicit modeling of the PSD boundaries to show the effects of AMPAR receptor confinement to the PSD on receptor trapping [9]. They vary the probability that AMPAR are reflected back into the PSD upon encounter of this boundary instead of diffusing into the extra synaptic space. Increased confinement of AMPAR to the PSD decreases the time needed to bind most receptors to the scaffold due to increased probabilities of AMPAR encountering binding scaffold elements. As shown before, a different macroscopic diffusion environment resembles this behavior and enriches receptors at the PSD without the need of explicitly modeling a reflective boundary around it.
The effects of the diffusion environment on explicitly modeled scaffold binding reactions were examined by running SmoldynM simulations with scaffolds while modeling two different diffusive environments for AMPAR. Results are shown in Fig. 2c. The parameters for the model are taken from [9] (D \(_{\text{ ESM }}\) = 0.45 µm\(^2\) s\(^{1}\), n(AMPAR): 66, n(Scaffold) = 132). The extra synaptic diffusion coefficient is the same for all simulations (D\(_{\text{ ESM }}\) = 0.45 µm\(^2\) s\(^{1}\)). The diffusion coefficient for the PSD is changed tenfold from D\(_{\text{ PSD }}\) = 0.45 µm\(^2\) s\(^{1}\) to D\(_{\text{ PSD }}\) = 0.45 10\(^{5}\) µm\(^2\) s\(^{1}\). Scaffolding elements and AMPAR binding to them is modeled explicitly. The results indicate the positive effect of slower AMPAR diffusion at the PSD on the overall time that is needed to trap AMPAR by scaffold binding at the PSD. This improvement however has its limitations once diffusion at the PSD gets too slow to enable AMPAR to encounter the scaffold elements.
Conclusion
We described, validated and illustrated surfacedependent diffusion coefficients as an extension to the simulator Smoldyn.
Examples demonstrate that SmoldynM behaves as expected compared to Smoldyn when such a comparison is directly possible. The new extension enables simulations showing interesting aspects such as anomalous diffusion due to different diffusive areas or the enhancement of reaction speed by enriching substrate molecules.
Simulations of AMPAR trafficking at the PSD show that a change in diffusion coefficient can substitute for explicit modeling of binding reactions, especially in cases where binding partners are unknown, not well defined or an abundance of unspecific associations slows the movement of molecules down. Especially in case of the latter, modeling the behavior based on macroscopic diffusive behavior might be the easier and also biologically more realistic abstraction. It has to be kept in mind that the results are not the same with regards to the spatial distribution of the trapped molecules. Modeling molecular trapping by unspecified scaffolds defines the position of receptors via the positioning of the scaffold. Trapping them via a change in diffusivity initially concentrates them on the boundaries of the subarea. Achieving a random or uniform distribution of molecules via macroscopic diffusion in that subarea might take longer since the molecules have to distribute within this subarea via the slower diffusion coefficient.
In general, being able to easily model membranes as surfaces composed of subareas that exhibit different properties is an important feature to study the influence of special subdiffusive membrane environments like the PSD.
Availability and requirements

Project name: SmoldynM

Project source code The extended Smoldyn can be found at http://lenoverelab.org/documents/smoldynM2.22s.tar.gz

Operating System: Platform independent

Programming Language: C

Other requirements: None

License: GPL

Any restrictions to use by nonacademics: No
Declarations
Authors’ contributions
CS and NLN designed the work and wrote the manuscript. CS wrote the software and ran simulations. Both authors read and approved the final manuscript.
Acknowledgements
This work was funded by the European Molecular Biology Laboratory graduate programme. Open access fees are provided by RCUK.
Competing interests
The authors declare that they have no competing interests.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 Tolle DP, Le Novère N. Particlebased stochastic simulation in systems biology. Curr Bioinform. 2006;1(3):315–20.View ArticleGoogle Scholar
 Nicolau DV, Hancock JF, Burrage K. Sources of anomalous diffusion on cell membranes: a Monte Carlo study. Biophys J. 2007;92(6):1975–87.PubMed CentralView ArticlePubMedGoogle Scholar
 Espinoza FA, Wester MJ, Oliver JM, Wilson BS, Andrews NL, Lidke DS, Steinberg SL. Insights into cell membrane microdomain organization from live cell single particle tracking of the IgE high affinity receptor FcIuRI of mast cells. Bull Math Biol. 2012;74(8):1857–911.PubMed CentralView ArticlePubMedGoogle Scholar
 Kennedy MB. Signalprocessing machines at the postsynaptic density. Science. 2000;290(5492):750–4.View ArticlePubMedGoogle Scholar
 Triller A, Choquet D. Synaptic structure and diffusion dynamics of synaptic receptors. Biol Cell. 2003;95(7):465–76.View ArticlePubMedGoogle Scholar
 Jaskolski F, Henley J. Synaptic receptor trafficking: the lateral point of view. Neuroscience. 2009;158(1):19–24.PubMed CentralView ArticlePubMedGoogle Scholar
 Choquet D, Triller A. The dynamic synapse. Neuron 2013;80(3).Google Scholar
 Borgdorff AJ, Choquet D. Regulation of AMPA receptor lateral movements. Nature. 2002;417(6889):649–53.View ArticlePubMedGoogle Scholar
 Tolle DP, Le Novère N. Brownian diffusion of AMPA receptors is sufficient to explain fast onset of LTP. BMC Syst Biol. 2010;4:25.PubMed CentralView ArticlePubMedGoogle Scholar
 Tolle DP, Le Novère N. Meredys, a multicompartment reactiondiffusion simulator using multistate realistic molecular complexes. BMC Syst Biol. 2010;4(1):24.PubMed CentralView ArticlePubMedGoogle Scholar
 Andrews SS, Addy NJ, Brent R, Arkin AP. Detailed simulations of cell biology with Smoldyn 2.1. PLoS Comput Biol. 2010;6(3):1000705.View ArticleGoogle Scholar
 Berg HC. Diffusion: Microscopic theory. In: Random Walks in Biology. Princeton: Princeton University Press. 1993. p. 5–14.Google Scholar
 Banks DS, Fradin C. Anomalous diffusion of proteins due to molecular crowding. Biophys J. 2005;89(5):2960–71.PubMed CentralView ArticlePubMedGoogle Scholar
 Jain P, Bhalla US. Signaling logic of activitytriggered dendritic protein synthesis: an mTOR gate but not a feedback switch. PLoS Comp Biol 2009;5(2).Google Scholar
 Newpher TM, Ehlers MD. Glutamate receptor dynamics in dendritic microdomains. Neuron. 2008;58(4):472–97.PubMed CentralView ArticlePubMedGoogle Scholar