Meredys, a multi-compartment reaction-diffusion simulator using multistate realistic molecular complexes
© Tolle and Le Novère; licensee BioMed Central Ltd. 2010
Received: 15 May 2009
Accepted: 16 March 2010
Published: 16 March 2010
Most cellular signal transduction mechanisms depend on a few molecular partners whose roles depend on their position and movement in relation to the input signal. This movement can follow various rules and take place in different compartments. Additionally, the molecules can form transient complexes. Complexation and signal transduction depend on the specific states partners and complexes adopt. Several spatial simulator have been developed to date, but none are able to model reaction-diffusion of realistic multi-state transient complexes.
Meredys allows for the simulation of multi-component, multi-feature state molecular species in two and three dimensions. Several compartments can be defined with different diffusion and boundary properties. The software employs a Brownian dynamics engine to simulate reaction-diffusion systems at the reactive particle level, based on compartment properties, complex structure, and hydro-dynamic radii. Zeroth-, first-, and second order reactions are supported. The molecular complexes have realistic geometries. Reactive species can contain user-defined feature states which can modify reaction rates and outcome. Models are defined in a versatile NeuroML input file. The simulation volume can be split in subvolumes to speed up run-time.
Meredys provides a powerful and versatile way to run accurate simulations of molecular and sub-cellular systems, that complement existing multi-agent simulation systems. Meredys is a Free Software and the source code is available at http://meredys.sourceforge.net/.
The influence of geometry and space on the functioning of cellular processes, the vast quantity of potential interactions due to molecular complex formation, and the stochasticity caused by low copy numbers of molecular species are all recognised features of many biological systems [1–3]. Within the field of computational biology these systems are best modelled using a particle-based stochastic approach . Here we present Meredys (MEsoscopic REaction DYnamics Simulator), a stochastic, particle-based simulation software designed to model and simulate reaction-diffusion systems at the mesoscopic level. The software is derived from an idea initially developed by Dan Mossop and Fred Howell in the Abstracted Protein Simulator (APS) . It is implemented in the Java programming language and uses Java3D as visualization framework for rendering to the screen. The input to the software is a model of a reaction diffusion system encoded in a Meredys specific implementation of the NeuroML model description language . The specification includes entries for molecule geometry and position, feature states of molecular entities, position of reaction sites, as well as types of reactions occurring and the biophysical properties of the diffusion landscapes. During a simulation, the software implements a Brownian Dynamics algorithm  to simulate the evolution of the system through time. Among the features Meredys was designed to tackle are the accurate simulation of reaction-diffusion systems operating in three dimensions, the potential for multi-state, multi-component molecular species, and the effect of multiple molecular states on the rate and outcome of the reactions these molecular species undergo.
Random Number Generation
Representation of molecules entities
Particles are the basic building blocks for the construction of compound objects. They contain the sites of all bi-molecular and some uni-molecular reactions. A particle's centre of mass is described by a position vector relative to the centre of mass of its parent entity (see P1 and P2 relative to E1 in Figure 5A).
Entities are permanent objects that never dissociate into their component particles during run-time. An entity's particle make up is defined by the modeller in the simulation input file. An entity maintains its identity throughout a simulation, even when it is part of a cluster comprising two or more entities. Entities have a centre of mass encoded as a position vector relative to the centre of mass of its parent cluster (see E1 and E2 relative to C1 in Figure 5B).
The final member of the component hierarchy of the objects used in Meredys is the cluster. Clusters are the run-time instantiations of one or more entities. Entities which undergo binding reactions are considered part of the same cluster (Figure 1B). This association can be transient. When two bonded entities separate, they each form an independent cluster. Additionally, a cluster's hydro-dynamic radius, used in the calculation of the cluster diffusion constant, is determined by the hydro-dynamic radii of all its member particles. The member particles are assumed to be spheres of a volume calculated from their user-defined hydro-dynamic radius. The parent entity is assumed to be a sphere of volume equal to the sum of the volumes of its child particles. The parent cluster is assumed to be a sphere of volume equal to the sum of the volumes of its child entities. The sphere's radius is taken as the cluster's hydro-dynamic radius (see Figure 5B). Clusters have a centre of mass which is a position vector relative to the centre of the simulation volume.
Clusters are the software objects which diffuse through the simulation volume at each time step. The diffusion properties of the cluster are determined by the particles composing the cluster's entities. Each particle belongs to a specific diffusion landscape. An entity then contains a set of diffusion landscapes composed of the diffusion landscapes of all the entity's child particles. A cluster, in turn, contains a set of diffusion landscapes constructed from the union of the sets of its child entities. The diffusion landscape which determines the cluster's diffusion properties is the most limiting landscape from the set of landscapes. Currently there are five landscapes to chose from: unrestricted, membrane, above membrane, below membrane, and static. The 'unrestricted' diffusion landscape allows for unrestricted diffusion over the whole simulation volume. Clusters in the 'membrane' diffusion landscape have their movement restricted to two dimensions, with no movement along the y-axis. This is used for movement of membrane-bound molecules. The membrane position point is the y-coordinate at which the membrane is located. It is user-defined. The membrane is therefore represented by a plane in x-z, which crosses the simulation volume at a user defined y-coordinate. The 'above membrane' diffusion landscape allows unrestricted three-dimensional diffusion in the sub-volume above the membrane position point. Conversely, the 'below membrane' diffusion landscape allows unrestricted three-dimensional diffusion in the sub-volume below the membrane position point. The 'static' diffusion landscape disallows any kind of movement. The 'static' diffusion landscape is set as the most limiting diffusion landscape, followed by the membrane landscape as the second most limiting followed by the remaining landscapes, which are considered of equal precedence.
A cluster's diffusion landscape can change during a simulation, as the cluster incorporates more entities containing a different set of diffusion landscapes, or when a cluster separates into two clusters and the inheritance of possible diffusion landscapes is unequal due to unequal entity composition. For example, a cluster can change from an unrestricted diffusion landscape to a static diffusion landscape by binding a different, static cluster. In addition to determining the limits of the cluster diffusion within the simulation volume, diffusion landscapes also influence the actual displacement a cluster experiences at each time step, by affecting the cluster's Diffusion coefficient, D.
where γ is Euler's constant, and h is the thickness of the plasma membrane (5 nm in Meredys), μ is the viscosity of the membrane, and r is the radius of a cylindrical particle in the membrane. The membrane landscape can be further sub-divided by defining membrane domains. These are circular sub-domains within the membrane. Membrane domains may have different viscosities from the membrane landscape. The user can assign specific boundary conditions to the boundaries between membrane domains and the membrane landscape.
In addition to translational motion, clusters also undergo rotational motion during each time step. As a cluster's rotational motion is much faster than its translational motion, clusters assume a random orientation after each time-step. Rotation is restricted for clusters diffusing in the membrane diffusion landscape.
Simulations take place in a simulation volume of user defined size delimited by the simulation volume boundaries. Additionally, specific membrane domains can be described which are separated from the canonical membrane environment by user-defined boundaries. As a consequence, types of behaviour need to be specified to resolve interaction of diffusing clusters with the available boundaries. There are four types of boundary interactions that can be simulated: Open, absorbing, periodic, and reflective. Boundary interactions are invoked whenever a cluster crosses the boundary. A cluster is said to have crossed a boundary if it is found on a different side of the boundary at the end of the movement step of the iteration cycle, compared to the start of the movement step. Open boundary interactions do not obstruct cluster diffusion at all. A cluster is freely allowed to cross an open boundary. Clusters crossing an absorbing boundary are removed from the simulation. Periodic boundary interaction allows the translation of the cluster across the domain volume to emerge at the opposite side. A reflective boundary interaction reflects the molecule according to the law of reflection. Any described boundary can have a number of boundary conditions associated with it. A boundary condition is defined as a boundary interaction type and an associated probability. The sum of all the probabilities of a domains boundary conditions must equal to one. The probabilities determine what type of interaction occurs when a cluster comes in contact with a boundary. Additionally, a boundary can posses a different set of boundary conditions, depending on the directionality of the crossing. For example, a boundary between two domains A and B may be open to molecules crossing from A into B, but reflective for molecules attempting to cross from B into A.
Meredys is capable of simulating zeroth-order reactions, uni-molecular reactions and bi-molecular reactions. Reactions involving three reacting partners simultaneously, tertiary reactions, cannot be simulated.
Where k is the reaction rate in units of Molar per second, Ms-1, δt is the time-step in seconds, V is the volume of the landscape the entities are created in and N Avogadro is Avogadro's number. At program initialisation, a Poisson distribution with mean λ is used to determine the time elapsed until creation of one entity. This is repeated until the total elapsed time is equal to the time-step of one iteration, δt. All the resulting new molecules are stored in a table and indexed by the iteration at which they are created. This process is repeated until the iteration step reached equals the total simulation run time. Since this process occurs at program initialisation, during the simulation run time only the relevant table entry needs to be queried at specific iteration steps, thus avoiding the need for computationally expensive random number generation during run time. The comparison of Meredys' zeroth order reactions and their analytic equivalent is provided in additional file 1.
k is the first order reaction rate constant, in units of s-1, and t is the elapsed time in seconds. This time is added to the elapsed simulation time to calculate the iteration step at which the reaction will occur. The reaction with associated time of occurrence is termed a reaction event. The reaction event is stored in a table indexed by the iteration step at which the event occurs. Additionally, reactions occurring during an iteration step are executed according to the order within that iteration step. Whenever an entity, reaction site or bond is created, during program initialisation or as result of a reaction, for example, Meredys determines the uni-molecular reactions the reactant can undergo and creates a reaction event for each reaction and adds it to the scheduler. The event is executed at its determined iteration step. If the reactant undergoes a state change or other reaction which affects a previously determined reaction event, then the affected reaction event is removed from the scheduler and a new reaction event determined if need be. The comparison of Meredys' first order reactions and their analytic equivalent is provided in additional file 1.
The simulated bi-molecular reactions take place on reaction sites. A reaction site software object is contained within the particle software object. A particle may contain more than one reaction site. The reaction sites are roughly analogous to biological binding sites or enzyme active sites. As active sites they determine the site of reactions for a particle and its parent entity. As binding sites, they determine the site of binding and geometry of the binding between two entities. The reaction site is described by a set of three points (see Figure 5). The first point, the centre point, gives the centre of the reaction site. The first and second points together describe a vector, the normal vector, through the centre of the reaction site. The first and the third points together describe a vector perpendicular to the normal vector, called the plane vector. The centre point is used as the centre of the sphere describing the reaction radius of the reaction surface for the purpose of bi-molecular reactions. The two vectors are used to determine the geometry of binding during binding reactions.
Meredys implements the bi-molecular reaction algorithm outlined in Andrews and Bray . This algorithm is based on the Smoluchowski model for reaction-diffusion systems . Within a physical system, a collision occurs when the reactant centres are separated by a distance equal to the sum of the molecular radii. Not every collision in a physical system leads to a reaction, as not every collision overcomes the reactions activation barrier. In order to take account of this, the algorithm replaces the sum of the molecular radii by an effective binding radius, σ. Andrews and Bray  determine the binding radius by deriving the simulated reaction rate constant in terms of the binding radius, equating this to the experimentally observed rate constant and then inverting the result to get σ. Bi-molecular reactions occur when two reaction site centre points come within a distance determined by the reacting pairs binding radius following the molecular displacement step of the iteration cycle.
The accuracy of Brownian Dynamics based simulators of chemical reactions is depended on the chosen time step length of the iteration step . In the case of the bi-molecular reaction algorithm used in Meredys, the authors of Smoldyn present a practical heuristic to allow the determination of an acceptable time step length . In brief, simulations are run with a trial time step at first instance. The simulations are then repeated using a time step half as long as the initial time step. If no significant differences are found between the results obtained using the two different time steps, the longer time step is sufficient. Otherwise the procedure is repeated.
Many molecular biological species interact to form transient complexes, such as the protein-protein interactions that dominate cellular signalling networks. Two procedures exist to simulate bi-molecular reactions resulting in bond formation between the two participating entities. Although the modeller has the option of encoding binding reactions by using the aforementioned bi-molecular reaction scheme and treating the bonded product as a new entity altogether, Meredys does allow for binding reactions where the identities of the participating entities are retained. This is particularly useful for the modelling of transient, reversible interactions, such as ligand binding to a receptor, as it eases the tracking of individual molecular species. The reaction scheme for binding reactions is that of the general bi-molecular reaction scheme given above. However, the reaction outcome differs, as a new cluster needs to be formed from the existing reacting partners. The structural rearrangements required for binding are encoded in the set of three points describing the reaction surface (see Figure 5). First, the centre points of the partner reaction sites are superimposed. Then the reaction partners are rotated to make their normal vectors anti-parallel. Finally, the reactants are rotated perpendicular to the normal vector to superimpose their plane vectors. The hydro-dynamic radii and diffusion landscapes of the reaction partners determine the relative contribution of each partner to the rotational movements required to bring them into the right orientation. Reactants with a overall hydro-dynamic radius rotate less, and both the membrane diffusion landscape as well as the static diffusion landscape restrict the amount of rotation a reactant can undergo. The comparison of Meredys' bi-molecular reactions and the equivalent simulation in an ODE solver is provided in additional file 1.
Many biological molecules can assume different states. Common examples include post-translational mod-ifications, ligand occupation or conformation states. These states often influence the molecules overall biophysical properties, including the reactions the molecule partakes in. At times, a modeller would like to keep track of a molecules different states but still maintain the identity of the original molecule; that is, avoid creating new entities every time a state change occurs. Meredys supports this concept, by allowing user defined feature states for simulation entities. An entity's state can have a direct effect on the reaction probability of any reaction the entity is capable of undergoing. Equally, any reaction the molecule undergoes can effect a state change. An entity's states are defined by describing a particular entity feature, such as channel gating, or phosphorylation site, and an enumeration of the possible states the feature can assume.
Sometimes the spontaneous creation through the zero-order reaction mechanism is not sufficient for the addition of new entities. It is possible to load pre-defined entities directly into the simulation at a given time point. These non-movement and non-reaction occurrences are termed events and can be specified within the input XML file.
Output is described in the NeuroML output class. This class defines for which entity we want output, at what time points this output is created and what kind of output to produce. The options are position, orientation, feature state and entity count. Output is printed whenever the iteration step number is evenly divisible by the number given in the timepoints attribute. Output is only printed if an iteration step is executed. The format is text with each line showing one cluster - entities user ids as well as internal entity unique identifiers of entities belonging to the cluster are given. In addition the position, orientation or feature states are given, if specified as output by the user. Simple text output allows the data to be further analysed by user created scripts.
Results and Discussion
Meredys is a stochastic, particle-based software designed to simulate reaction-diffusion systems at the meso-scopic level. Molecular composition and states are maintained for each simulated molecular species within the simulation. The random walk algorithm effectively models the molecular diffusion of molecular species in two- and three-dimensions. The random walk is influenced by the diffusion landscape the moving cluster operates in. The membrane diffusion landscape can further be subdivided into domains such as the extra-synaptic or synaptic membrane. Meredys offers algorithms for simulating uni-molecular reactions, such as receptor-channel state changes, and bi-molecular reactions, such as ligand binding reactions, as well as zeroth-order reactions to model creation processes without the need to fully specify all of the entities involved in the creation process. In the case of binding reactions between two molecular species, Meredys allows for the control of the geometry of the interaction. Additionally, interacting partners maintain their identities throughout the interaction. Networks with over 4000 components have been simulated using Meredys (data not shown). Run time is highly dependent on the compute power, number of molecules involved in the simulation and number of interactions the molecules can undergo.
There are a number of simulation engines available to the biochemical modeller [11, 15–19]. Often, the simulators are classified into groups based on the amount of spatial detail which they are able to simulate . Depending on the system of interest, particle-based simulators may prove to be the better choice over other simulation engines. Particle-based stochastic simulators like Meredys are capable of capturing a host of features which population based simulations engines are not capable of capturing due to the lack of spatial information presented in the model . Examples include analyses of glutamate release location relative to post-synaptic receptor location during neurotransmitter release in the glutamate synapse , or investigating the effect of signal molecule locations on signalling during bacterial chemotaxis . A feature which sets Meredys apart from these other particle-based simulation software such as MCell  and Smoldyn  is Meredys' multi-component, multi-state feature clusters. Many biological entities form large, interacting multi-component clusters . Meredys' multi-component cluster formation allows for the identification and tracking of specific members of the clusters during the whole of the simulation.
Availability and requirements
We thank Dan Mossop and Fred Howell for providing the Abstract Protein Simulator, the software from which Meredys was derived. We also thank Simone Zorzan for helpful discussions of the Meredys source code, Sarah Birch for careful reading of the manuscript and Anton Enright for aid with the creation of the images. Both DPT and NL were supported by the European Molecular Biology Laboratory.
- Lemerle C, Ventura BD, Serrano L: Space as the final frontier in stochastic simulations of biological systems. FEBS Letters. 2005, 579: 1789-1794. 10.1016/j.febslet.2005.02.009View ArticlePubMedGoogle Scholar
- Lok L, Brent R: Automatic generation of cellular reaction networks with Moleculizer 1.0. Nat Biotechnol. 2005, 23: 131-136. 10.1038/nbt1054View ArticlePubMedGoogle Scholar
- Halling PJ: Do the laws of chemistry apply to living cells?. Trends Biochem Sci. 1989, 14 (8): 317-318. 10.1016/0968-0004(89)90158-8View ArticlePubMedGoogle Scholar
- Tolle DP, Le Novère N: Particle-Based Stochastic Simulation in Systems Biology. Curr Bioinformatics. 2006, 1: 315-320. 10.2174/157489306777827964.View ArticleGoogle Scholar
- Mossop D, Howell F: Abstracted Protein Simulator. 2001, http://www.neurogems.org/protsim1/Google Scholar
- Goddard NH, Hucka M, Howell F, Cornelis H, Shankar K, Beeman D: Towards NeuroML: model description methods for collaborative modelling in neuroscience. Philos Trans R Soc Lond B Biol Sci. 2001, 356 (1412): 1209-1228. 10.1098/rstb.2001.0910PubMed CentralView ArticlePubMedGoogle Scholar
- Ermak DL, McCammon JA: Brownian dynamics with hydrodynamic interactions. J Chem Phys. 1978, 69: 1352-1360. 10.1063/1.436761.View ArticleGoogle Scholar
- Knuth D: The Art of Computer Programming. 1999, Reading, Massachusetts, Addison Wesley,Google Scholar
- Devroye L: Non-Uniform Random Variate Generation. 1986, New York, Springer-Verlag,View ArticleGoogle Scholar
- Berg HC: Random walks in biology, expanded edition. 1993, Princeton: Princeton University Press,Google Scholar
- Andrews SS, Bray D: Stochastic simulation of chemical reactions with spatial resolution and single molecule detail. Phys Biol. 2004, 1 (3-4): 137-151. 10.1088/1478-3967/1/3/001View ArticlePubMedGoogle Scholar
- Stiles JR, Helden DV, Bartol TM, Salpeter EE, Salpeter MM: Miniature endplate current rise times less than 100 microseconds from improved dual recordings can be modeled with passive acetylcholine diffusion from a synaptic vesicle. Proc Natl Acad Sci USA. 1996, 93 (12): 5747-5752. 10.1073/pnas.93.12.5747PubMed CentralView ArticlePubMedGoogle Scholar
- Saffman PG, Delbrück M: Brownian motion in biological membranes. Proc Natl Acad Sci USA. 1975, 72 (8): 3111-3113. 10.1073/pnas.72.8.3111PubMed CentralView ArticlePubMedGoogle Scholar
- v Smoluchowski M: Versuch einer mathematischen Theorie der Koagulationskinetik kolloider Lösungen. Zeitschrift für physikalische Chemie. 1916, XCII: 129-168.Google Scholar
- Stiles J, Bartol T: Computational Neuroscience: Realistic Modeling for Experimentalists. 2001, CRC Press, Boca Raton,Google Scholar
- Schaff J, Fink CC, Slepchenko B, Carson JH, Loew LM: A general computational framework for modeling cellular structure and function. Biophys J. 1997, 73 (3): 1135-1146. 10.1016/S0006-3495(97)78146-3PubMed CentralView ArticlePubMedGoogle Scholar
- Hattne J, Fange D, Elf J: Stochastic reaction-diffusion simulation with MesoRD. Bioinformatics. 2005, 21 (12): 2923-2924. 10.1093/bioinformatics/bti431View ArticlePubMedGoogle Scholar
- Tomita M, Hashimoto K, Takahashi K, Shimizu TS, Matsuzaki Y, Miyoshi F, Saito K, Tanida S, Yugi K, Venter JC, Hutchison CA: E-CELL: software environment for whole-cell simulation. Bioinformatics. 1999, 15: 72-84. 10.1093/bioinformatics/15.1.72View ArticlePubMedGoogle Scholar
- Le Novère N, Shimizu TS: STOCHSIM: modelling of stochastic biomolecular processes. Bioinformatics. 2001, 17 (6): 575-576. 10.1093/bioinformatics/17.6.575View ArticlePubMedGoogle Scholar
- Franks KM, Stevens CF, Sejnowski TJ: Independent sources of quantal variability at single glutamatergic synapses. J Neurosci. 2003, 23 (8): 3186-3195.PubMed CentralPubMedGoogle Scholar
- Lipkow K, Andrews SS, Bray D: Simulated diffusion of phosphorylated CheY through the cytoplasm of Escherichia coli. J Bacteriol. 2005, 187: 45-53. 10.1128/JB.187.1.45-53.2005PubMed CentralView ArticlePubMedGoogle Scholar
- Bartol TM, Land BR, Salpeter EE, Salpeter MM: Monte Carlo simulation of miniature endplate current generation in the vertebrate neuromuscular junction. Biophys J. 1991, 59 (6): 1290-1307. 10.1016/S0006-3495(91)82344-XPubMed CentralView ArticlePubMedGoogle Scholar
- Bray D: Genomics. Molecular prodigality. Science. 2003, 299 (5610): 1189-1190. 10.1126/science.1080010View ArticlePubMedGoogle Scholar