 Methodology Article
 Open Access
 Published:
Potential based, spatial simulation of dynamically nested particles
BMC Bioinformatics volume 20, Article number: 607 (2019)
Abstract
Background
To study cell biological phenomena which depend on diffusion, active transport processes, or the locations of species, modeling and simulation studies need to take space into account. To describe the system as a collection of discrete objects moving and interacting in continuous space, various particlebased reaction diffusion simulators for cellbiological system have been developed. So far the focus has been on particles as solid spheres or points. However, spatial dynamics might happen at different organizational levels, such as proteins, vesicles or cells with interrelated dynamics which requires spatial approaches that take this multilevelness of cell biological systems into account.
Results
Based on the perception of particles forming hollow spheres, MLForce contributes to the family of particlebased simulation approaches: in addition to excluded volumes and forces, it also supports compartmental dynamics and relating dynamics between different organizational levels explicitly. Thereby, compartmental dynamics, e.g., particles entering and leaving other particles, and bimolecular reactions are modeled using pairwise potentials (forces) and the Langevin equation. In addition, forces that act independently of other particles can be applied to direct the movement of particles. Attributes and the possibility to define arbitrary functions on particles, their attributes and content, to determine the results and kinetics of reactions add to the expressiveness of MLForce. Its implementation comprises a rudimentary rulebased embedded domainspecific modeling language for specifying models and a simulator for executing models continuously. Applications inspired by cell biological models from literature, such as vesicle transport or yeast growth, show the value of the realized features. They facilitate capturing more complex spatial dynamics, such as the fission of compartments or the directed movement of particles, and enable the integration of nonspatial intracompartmental dynamics as stochastic events.
Conclusions
By handling all dynamics based on potentials (forces) and the Langevin equation, compartmental dynamics, such as dynamic nesting, fusion and fission of compartmental structures are handled continuously and are seamlessly integrated with traditional particlebased reactiondiffusion dynamics within the cell. Thereby, attributes and arbitrary functions allow to flexibly describe diverse spatial phenomena, and relate dynamics across organizational levels. Also they prove crucial in modeling intracellular or intracompartmental dynamics in a nonspatial manner, and, thus, to abstract from spatial dynamics, on demand which increases the range of multicompartmental processes that can be captured.
Background
Space plays an important role in cell biological dynamics, such as cell signalling [1]. To study cell biological phenomena which depend on diffusion, active transport processes, or the locations of species, modeling and simulation studies need to take space into account. Over the last decade, a variety of spatial modeling and simulation methods and tools has been developed to support these simulation studies. They offer different approaches how to describe a spatial model, e.g., reactionbased, rulebased, or graphically, and they differ referring to what kind of spatial dynamics can be described, e.g., whether concentrations, populations, or individual particles are considered, whether a deterministic or stochastic approach is pursued, and whether movement takes place in discretized or continuous space [2, 3]. The kind of spatial dynamics that are supported by tools largely determines the results and questions that can be answered by a simulation study [4–6] and has been subject to different categorizations to structure the portfolio of spatial simulation methods from which the user can select.
TAKAHASHI et al. [7] distinguished spatial simulation approaches according to the representation of space. We adopt this structuring and slightly adapt it for the purpose of this paper as presented in Fig. 1. We distinguish between:
 1.
Stochastic nonspatial dynamics: To take the stochasticity of the modelled system into account, the propensity of reactions are sampled from exponential distributions and determine which reaction will occur at which time [9].
 2.
Compartmental dynamics: Many tools allow to constrain the dynamics of species to specific compartments. However also, compartments themselves can be subject to dynamics. For example, compartments can fuse and divide [10].
 3.
Individual particles moving in continuous space: Particles can be identified by their unique position in space, bimolecular reactions are triggered by collisions of particles, and typically particles diffuse by Brownian motion [11, 12].
 4.
Partial differential equations: Spatial gradients of concentrations are calculated deterministically [13].
 5.
Spatial stochastic dynamics: Multiple particles can occupy a position within a lattice space, and diffuse between neighboring positions [14].
Particlebased approaches (in the above categorization c)) are the subject of a further categorization suggested by SCHöNEBERG et al. [8] (Fig. 1, right lower corner). Spatial particlebased simulation tools are categorized into those that support
level 1 – free diffusion: basic 3D diffusion of point particles and reactions between them are considered,
level 2 – confined diffusion: diffusion can be constrained to compartments,
level 3 – excluded volume: point particles are replaced by volumetric entities, and
level 4 – potentials for particleparticle interaction: instead of a simple rejection of movements assuming rigid cells, potentials are used to determine the interaction of particles in terms of exclusion of movement or reactions.
According to this categorization [8], Smoldyn is characterised as level 2 particlebased approach: point particles which do not hamper each others movement are equipped with binding and unbinding radii, they can be constrained to compartments. A more recent version includes collision strategies to reach level 3 [15]. SpringSaLaD [16], ReaDDy [17], and SRSim [18] are based on forces to model the particleparticle interactions (level 4). Still, the levels introduced do not imply that a tool that supports level 3 provides all interesting features that a tool working at level 2 offers to study the system of interest. E.g., DONOVAN [19], whose spatial simulation works at level 2, supports arbitrarily complex 3D mesh geometries, whereas in MLSpace [20] (level 3) only rigid spheres are considered.
Merging both characterizations, approaches that combine different spatial representations can be characterized. The TwoRegime method [21] and KLANN et al. [22] allow to combine particle and RDME dynamics within one model: c(3)e. Similarly in [23], particles and partial differential equations are combined to focus on the spatial region of interest: c(3)d. [24] couples a partial differential equations solver and Smoldyn c(2)d. MLSpace combines particles at level 3 (as excluded volumes are considered), compartmental dynamics, and spatial, stochastic simulation (RDME): bc(3)e [20].
In addition, simulation environments offer the possibility to select different spatial semantics for one model specification. For example, in VCELL [13], rulebased models can be interpreted by a particlebased simulator (i.e., Smoldyn) or a partial differential equation solver, both confined to realistically geometrical compartmental structures (c(2), d). This list is far from being complete, and the characterization of the simulation tools may only depict a specific state in their development.
Against this background, we propose a new spatial particlebased modeling and simulation approach for cell biological systems, i.e., MLForce (Fig 1). It combines:
Potentialbased particle dynamics: The excluded volumes that are represented by the large amount of macromolecules within the cell affect the physicochemical kinetics of various intracellular processes [25]. To capture the effects of molecular crowding, space exclusions need to be considered. In addition, interactions based on potentials are an effective means to study the formation of clusters [18]. Additionally, potentials allow us to capture directed movements, such as the transport of vesicles [26] (see “Results” section). Similar to ReaDDy [17], SRSim [18], or SpringSaLaD [16], MLForce will use potentials for particleparticle interaction: c(4).
Compartmental dynamics: Intracellular space is further structured by compartments and vesicles, most of which are subject to frequent changes in terms of numbers, content, and interconnectivity. A prominent example is the endosomal system. Vesicles form at the membrane. Here they acquire their cargo and engulf protein receptor complexes. Those are transported towards the inner cell, where part is degraded and part is recycled. The vesicles themselves move through stages of early sorting, recycling to late endosomes, closely interacting with each other – processes which include frequent fission and fusion [26] (Fig. 2). MLForce uses the same forcebased approach that it uses for the interaction of particles and global force functions, to model the compartmental dynamics. Similar to MLSpace, MLForce supports both: compartmental and particle dynamics. However, unlike MLSpace it uses potentials bc(4).
Rulebased approach with arbitrary attributes and functions: As other spatial simulators, MLForce supports rulebased modeling [29]. However in MLForce, attributes are not constrained to a finite set of values. Similarly, as in MLRules [30], attributes are of arbitrary type and arbitrary functions are supported to determine the result and kinetics of reactions. The possibility to define attributed species and arbitrary functions increases significantly the expressiveness of a modeling formalism [31]. It allows to express spatial stochastic models (RDME) on a lattice (e) within a modeling approach that relies on a standard stochastic simulation algorithm for execution (a). By attributing species with indices that locate them in a specific point in the lattice (which represent a subvolume), e.g., A(x,y),B(x,y), reactions can be constrained to species within one subvolume, e.g., A(x_{a},y_{a})+B(x_{b},y_{b})→2 B(x_{b},y_{b})@ if x_{a}=x_{b}∧y_{a}=y_{b} then k else 0. Furthermore diffusion between neighboring subvolumes can be defined, e.g., B(x_{b},y_{b})→B(oneOf((x+1,y),(h−1,y),(x,y+1),(x,y−1)))@k_{d}. Thus, attributed species and arbitrary functions can be used to add space to otherwise nonspatial approaches. This is an observation also made in colored PetriNets [32] and in expressive rulebased modeling languages such as MLRules [30].
In MLForce, equipping particles with attributes and arbitrary functions allows e.g. to describe some part of the cellular dynamics in a nonspatial manner. For example, the reaction B+C→D @ k which takes place in particle A can be transformed into a nonspatial stochastic first order reaction (a). Therefore, the particle A is equipped with the corresponding attributes and a reaction is defined A(b,c)→A(d) @ b·c·k. The function b·c·k calculates the propensity of the reaction to take place. The possibility to constrain reactions based on arbitrary functions defined on the attributes, allows us to include nonspatial stochastic dynamics (a) which brings MLForce up to abc(4) (see Fig. 1).
Methods
MLForce adopts a rulebased approach to specify the spatial biochemical dynamics. Firstly the mathematical concepts and abstractions are presented. This is followed by a short description of its implementation and embedded domainspecific modeling language.
Mathematical concepts and abstractions
As the name suggests, forces are central in MLForce: they govern the interactions between and most of the activity of particles in a continuous manner. This is in contrast to the discrete event view chosen by simulators such as MLSpace or other agentbased approaches which interpret compartmental dynamics such as the proliferation of cells as discrete events [20, 33].
Representation of particles
MLForce is a particlebased approach. Particles represent hollow spheres. As particles can contain other particles, they can be used to model cellular compartments. Particles are spherically symmetric and deformable, their size, properties and nested content may vary. Particles p_{m} are categorized into classes called species: p_{m}∈S_{i}. Species are characterized by a set of arbitrary attributes, S(a_{1},…a_{n}), e.g., radius, phosphorylation sites, or phases. In our implementation of MLForce we require particles to have at least the attributes radius r (nonzero) and density ρ. The information about its radius together with information about the system’s temperature and the viscosity of its environment allows to derive the diffusion coefficient of a particle. In combination with the density, this information is used to determine the actual velocity (see “Spatial propagation” section). Particles may contain other particles, p_{m}(v_{1},…v_{n})[p_{j}+…+p_{k}]. All particles belonging to one class share the same attributes and attributetypes. Attribute values and content of particles may vary among particles of the same species. The interactions between and the activity of particles is governed by rules of behavior that operate on the particles, their attributes and their content. They rely on arbitrary functions for accessing and updating attributes and content of particles. Particles belonging to the same species share the same behavioral patterns.
Spatial propagation
For particles to interact they need to be in close proximity. In MLForce, we base our model of propagation on the LANGEVIN equation [34]
with the dot denoting time derivative and boldness indicating vectors. A Langevin equation is an ordinary differential equation with a random term added. This equation can be intuitively understood as an extension to NEWTON’s \(\boldsymbol {F}(\boldsymbol {x},t)=m\boldsymbol {\ddot {x}}(t)\) by adding both friction (friction coefficient γ) which is linear to the velocity and a random force f which describes the accumulated effect of the system interaction with the individual particle. This interaction between the particles and the system results in the diffusion of the particles. In particular f is commonly sampled from a Gaussian distribution to model the aggregate system behavior.
Values for the constants involved can generally be found using known equations like STOKESEINSTEIN relation (\( D_{i} =\frac {k_{B}T}{\gamma _{i}} =\frac {k_{B}T}{6\pi \eta r_{i}}\)), as well as system properties like viscosity η.
The F_{ext}(x,t) expression is the sum of all external forces on a particle. Two types of external forces are distinguished:
global forces are those, that act on one particle, independent of other particles. Global forces functions may be specified as arbitrary functions of the individual particle state (most importantly its position) and the global state of the simulator. This may be used to model, for example, some global current in a blood stream or some other directed cellular motion (as is done in “A model of vesicular transport” section below).
pairwise forces are those, that originate in the interaction of two particles like the nonreaction force (Eq. 2), reaction force (Eq. 7) or a division force (Eq. 6).
Nonreactive force
In order to model excluded volumes, we introduce a soft sphere potential between two nonreacting particles. This is a common and natural way to model excluded volumes in forcebased approaches. In addition, a hardsphere potential would imply more efforts numerically, as it would lead to a discontinuous force. An established choice for the modeling of elastic (softsphere) collision is the HERTZ collision theory [35]:
where d is the overlap of the particles, \(E^{*} = \left (\frac {1\nu _{1}^{2}}{E_{1}}+\frac {1\nu _{2}^{2}}{E_{2}}\right)^{1}\) results from the YOUNG’S modulus E and POISSON’S ratio ν and \(r =\left (\frac {1}{r_{1}}+\frac {1}{r_{2}}\right)^{1}\) from their radii. However, for many sparse applications, any kind of polynomial would suffice, as the exact dynamics of the elastic collision do not matter, when the mean time between collisions is sufficiently large.
Reactions in mLForce
Reactions between particles in MLForce include zero, first and second order reactions. MLForce supports a rulebased description of models. Whereas lower order reactions are sampled from an exponential distribution, second order reactions are calculated, similarly as the movement of particles, based on forces.
Lower order reactions
Simulation of zeroth and first order reaction, e.g., ⊥→p_{j}, or p_{i}→pi′, or p_{i}→p_{j}+p_{k} is a well understood problem (e.g., in SpringSaLaT [16]). The most basic approach is to sample the rate every (sufficiently small) timestep and check it against a random number.
In MLForce, rates can be defined by arbitrary functions that access the attributes and contents of the particles involved in a reaction, e.g., for a first order reaction:
here λ denotes a function defined by the user, k a parameter, which might refer to the system globally, e.g., in terms of temperature and viscosity, or to a kinetic constant. The result of the function λ serves as the parameter for the exponential random distribution, that describes the probability density of the time until the reaction occurs. Thus in MLForce, zero and first order reactions are treated semantically the same as in stochastic simulation algorithms [9]. In stochastic, nonspatial and spatial, simulation approaches (see Fig. 1a and e) exponentially distributed stochastic events form the basis of all, including second order, reactions.
In MLForce the rate is sampled every time step. Instead of propagating the simulation until the exact time point a reaction occurs, we check (stochastically), if the reaction occurs in the current time step. Sampling the rate every time step is useful if the information (v_{1},…,v_{n},[p_{1},…,p_{m}]), on which the rates depend, frequently changes. If this information remains largely constant over time, a scheduling approach may be more efficient computationally, as events do not have to be rescheduled frequently [36]. In executing MLForce models so far, zeroth and first order reactions have only been responsible for a small fraction of the required computingtime.
In addition to computational considerations, the error introduced by this procedure needs to be considered. Due to the discretization of time (explicit time steps) nomore than one reaction per particle may take place per time step. This is important to keep in mind for future models. In the application domain the mean time between reactions tend to be large compared to the very small time steps that result from the fast diffusion of particles. As long as this is the case, the error is negligible. With resonable technical effort, this problem could also be avoided by allowing multiple reactions per time step.
One kind of lower order reaction is the changing of particle size, in terms of a particle’s radius r. In MLForce arbitrary functions can be applied. The new value of the radius r, i.e., \(v^{\prime }_{1}\), of a particle p_{i}∈S_{j}(r,…) is calculated in dependence of some global parameters k, and the current state of the particle in terms of its attribute values v_{1},…v_{n}, and its content p_{1}…p_{m}, such that
By default, whenever a particle is created, it is first put into the system as a very small particle and then quickly grows continuously to its desired size. This continuous change of particle radii is also used, whenever the radius of a particle is changed in the model by any kind of reaction. Carrying the changes out suddenly (so not slowly over multiple integration steps) would potentially lead to integration errors, especially in systems with nesting. A sudden change in size could for example lead to a high particle overlap, which in turn would create a high unrealistic potential energy leading to a strong force and thus a disturbance in the system, e.g., numerical heating. Slow changes in size, allow for each artificial change in energy to slowly dissipate e.g., in friction due to the numerically very stable Langevin model. Naturally, it is also possible to create a particle within or at the surface of another particle. Those new positions are chosen randomly. At the surface this means that the new particle’s shell barely touches the old particle’s shell.
Compartment division
An interesting kind of first order reaction is that of a cellular compartment (or specifically a cell) dividing.
with y=k,v_{1},…v_{n},[p_{1},…,p_{m}] or in shorthand, omitting attributes and contained particles
In MLForce this is realized in a continuous and nesting compatible manner. As a first order reaction, the event that the division starts is sampled from an exponential distribution. The simulator temporarily introduces an “ignoring” relation between particles. This can be best explained using an example: Imagine a particle p_{i} containing many particles \(\hat {p}_{1} \ldots \hat {p}_{m}\). If now the reaction for p_{i} to split into p_{i1} and p_{i2} triggers, a few things happen. First of all, each particle \(\hat {p}_{1} \ldots \hat {p}_{m}\) gets assigned to remain in either p_{i1} or p_{i2}. By default, this assignment happens randomly. If required by the model, a specific splitting function may be specified by the user. Therefore, two arbitrary functions, i.e., λ_{c1},λ_{c2} are applied. The results of these applications form the new content of the newly generated particles, i.e., p_{i1} and p_{i2} respectively. The functions λ_{i1} to λ_{n1}, respectively λ_{i2} to λ_{n2} determine the new attribute values of the particles. If attribute values remain the same, simply the old values v_{i} can be used.
Each \(\hat {p}_{i}\) that will remain in p_{i1} is set to ignore the particle p_{i2} and vice versa. Also p_{i1} and p_{i2} are set to ignore one another. Now an external force is applied to push p_{i1} and p_{i2} apart
The force is calculated based on the acceleration a that is part of this particular reaction type (Eq.: 4, 5) in MLForce and the mass of the involved particles (which again depend on the respective radii and their density ρ). This is chosen, instead of a uniform force, to create a uniform motion of all particles. In principle however, any force expression suitable to the model could be used here, like a half sine wave for example. Once p_{i1} and p_{i2} are fully outside one another, all the newly introduced ignore relations are removed and regular propagation proceeds.
The principle of continuous change in particle size also applies here. At the beginning both particles have the same size. If p_{i1} and p_{i2} are to change their size throughout the division (see example in “The yeast model” section), this change is carried out throughout the division process. Initially, both p_{i1} and p_{i2} share position and radius. The further they move out, the closer their radii become to the desired ones. The fusion of compartments could be added in an identical (if reversed) fashion. However so far, it has not been implemented yet.
Bimolecular reactions
Handling bimolecular reactions in particlebased approaches is rather challenging. Whilst solutions exist when not considering excluded volumes (e.g. Smoldyn [12]), there has not yet been agreement among the community on the best way to approach this problem. This is evident from the many different approaches of recent tools [5, 12, 16, 17, 37]. Two processes contribute to the macroscopic rate,
the diffusion process describes how likely it is that two particle actually move into the vicinity of each other. This upper bound to the reaction rate is well described by the SmoluchowskiTheory [38] (elaborated in the Additional file 1)
the microscopic activation process describes how likely it is for two particle to react once they are close to each other [39, 40].
The rate is considered to be diffusion limited if every interaction results in a reaction, all other cases are considered to be activation limited. We propose to use forces for both processes. This means also the microscopic activation process is based on forces. In case a reaction can occur between two particles, the nonreactive force (Eq.: 2) has to be replaced by a reactive force. We used a rather simple approach which is calculated based on a degree of overlapping required d^{∗} and the work to cross the energy barrier between the two particles. Both, the degree of overlapping and the energy barrier, are specific for the reaction to occur and have to be defined by the modeler.
where d^{∗} denotes the required overlap of the particles and W the work required to cross the energy barrier. It should be noted due to the modular design (see “Simulator and implementation” section) within the simulator these functions can easily be changed. Suppose we were to look at the reaction p_{i}+p_{j}→p_{k}{d^{∗},W}. Whenever p_{i} and p_{j} overlap, they start to repel each other with the force F_{react}. Once they overlap more than d^{∗}, we consider the reaction to be triggered. Thus, the energy barrier and the required overlap of the microscopic interaction potential determines the likelihood of the reaction to take place. The higher the energy barrier, the less likely the particles will overlap sufficiently. If the energy barrier W is negligible, the diffusion is the only limit to the macroscopic reaction rate in the reaction process, it is diffusion limited.
This type of bimolecular reaction handling with forces, allows also for selfconsistent handling of dynamic nesting (and “unnesting”) reactions, i.e., one particle shuttles into (or out of) another particle. In principle the same concepts apply as for the regular bimolecular reaction. Looking at the reaction p_{i}+p_{j}→p_{i}[p_{j}]{d^{∗},W} (and the same for p_{i}[p_{j}]→p_{i}+p_{j}{d^{∗},W}), first the particles need to diffuse into proximity of one another. Once particles overlap, they repel each other. If p_{i}’s diameter has fully passed over p_{j}’s outline, the nesting reaction is triggered. However at this point p_{j} is already fully contained in p_{i} so there is no jerkiness in the simulation. For example, if there is no unnesting reaction defined for these particles, the only thing, that changes, is that, if p_{j} were to move towards again touching p_{i}, it would perceive the force F_{non−react} (pulling inwards). This ensures a softspherestyle excluded volume. Similar as for the regular p_{i}+p_{j}→p_{k} type bimolecular reaction, a high energy barrier makes for a less, and a low for a more likely reaction.
The advantage of this forcebased approach, is the smooth integration of the nesting process. On the other hand, it does require a very small timestep, compared to other methods. Furthermore, in the current situation, the parameterization of bimolecular reactions is not optimal. Ideally, the modeler would specify the macroscopic rate, instead of the microscopic barrier. In the future we plan on investigating this further and, using principles from thermodynamics, provide an automatic conversion method. This would mean, by modeling the system as a thermodynamic ensemble we could determine the statistical likelihood of a particle surpassing a certain energy threshold during interaction and put this into explicit mathematical relation to the reaction rate.
Simulator and implementation
The simulator has been implemented using C++14 and has been tested on both Windows and Linux (using CMakebuild system) and is available at https://git.informatik.unirostock.de/mosi/mlforcepublication.
The userinterface is built around a rudimentary embedded domain specific language (DSL). An example and explanation can be found in the provided code snippet below. In the long run, it would be advantageous to create an external DSL using techniques like transpiling, to still have the performance benefits of a compiling language. A more indepth comparison and discussion of DSLs can be found in the literature [41]. The developed embedded DSL bears some similarities with the ℓlanguage [42] which also supports an imperative rather than a declarative approach towards modeling.
As can be seen in Fig. 3, the software pursues a componentbased design with a clear separation of concern [43]. Its components are arranged around the main simulation engine. The simulation engine itself is very basic and designed to make as little assumptions about the underlying semantics as possible. This simulation engine interacts with the components via a well defined, small interface. Each of the key ingredients of the simulation is separated into exchangeable components, namely the integrator, collision detection engine, and the visualizer, as well as the overarching modeling layer which provides the interface between the simulation engine and DSL. They are invoked by the simulator and implement the functionality introduced in the previous section. The simulator is responsible for executing the lower and higher order reactions by calling the specified functions of the model. First it is checked whether lower order reactions trigger in this step. If so, those are executed. Based on the forces a numerical integrator calculates the new positions of the particles. To identify the particles that interact, a collision detection engines is invoked. Afterward the simulator calculates the bimolecular forces and executes the biomolecular reactions accordingly. The calculated state can be subject to visualizations or other reporting mechanisms.
Integrator Currently integration is carried out using a twostep kickdriftkick integrator [44] of the stochastic LANGEVIN equation. This integration includes velocity, and is thus fairly costly numerically. However as an alternative, one could select an integrator with immediate dissipation of momentum, that would allow for a considerably larger timestep. The choice for the more complex integrator was made as MLForce makes heavy use of various forces and carrying out the more complex integration, helps in terms of correctness. In the current implementation the time step is chosen adaptively based on the smallest particle and the fastest velocity currently in the system, as \(\Delta t = \frac {r_{\text {min}}}{v_{\text {max}}\,g}\) with g as the parameter of granularity, that describes the precision of the integration.
Collision detection engine Finding possible interactions of particles using the collision detection engine is one of the key drivers of runtime in the current implementation. Several collision detection engines have been implemented, using different algorithms and parallelization. For example, for highly nested particles, the collisiondetection problem lends itself to be solved by a specifically tailored algorithm and a massively parallel execution on the GPU [28].
Visualizer Currently 3 basic visualizers have been created. One to just create an output text file, and two based on CImg [45], one for live viewing and one for generating movie files. The visualizers are encapsulated in a separate thread to minimize overhead. One more advanced visualization has been implemented that maps the simulation to 3D space using OpenGL [46]. It offers some useful features for the case studies (see “Results” section), such as tracking individual particles in space.
Modeling layer The modeling layer makes heavy use of λfunctions in order to encapsulate model specific behavior into single function calls for the simulator. λfunctions here are a C++11 feature that allows to define anonymous functions, that can be used as parameters in a functionallike programming style. This is facilitated by various templates that have been defined in the DSL. For example, after each reaction, the simulator invokes the post reaction function (see e.g. afterFuncA in listing 2). In the modeling layer this function is specified by the user and denotes the results of applying a rule, such as particles changing state, being degraded or created.
This simplifies the implementation of the simulator and gives a unified and lean interface for implementing further features. The simulator also provides a set of units (via the units:: namespace) to the modeler, allowing for consistent parameterization. For simplicity the reactants in the DSL are always denoted as A and B. The particular type of A and B is then specified via the as_A etc. functions. In general the DSL is very expressive, as all rates and reactions are described as λfunctions. These can be of arbitrary complexity. However, it should be noted that the current DSL presents only a proof of concept rather than a language that allows a succinct description of MLForce models. Therefore, further efforts will be dedicated to improving the design of the language.
Results
In order to test the MLForce simulator and its implementation, we conduct a few case studies. These include simple test cases to investigate correctness and more complex models. The latter are inspired by realistic cell biological models, make use of the features in MLForce, and, thus, demonstrate their usefulness. Both, the basic test models and the more realistic ones, can be found in the Additional file 1.
Testing correctness
As other particlebased simulation approaches [16, 17, 37, 47] we test the correctness of our simulator, based on a set of simple reactions which we simulate to compare the achieved results to theory. To those tests belong the creation of a particle with a constant rate as a 0 order reaction (\(\emptyset \xrightarrow {\text {k}} A\)), the decay with a rate \(A \xrightarrow {\text {k}} \emptyset \) as first order reaction, and irreversible and reversible second order reactions A+B⇔C in the diffusion limited (intrinsic rate k_{int}=∞) case. A comparison of simulation results and theory shows overall a good agreement with the theory (see Additional file 1). However, those tests also reveal the decisive role of selected parameters to determine forces (“Bimolecular reactions” section) and time steps. For their selection, suitable computational support should be provided. As a form of preprocessing step, energy barriers respectively accelerations for individual reactions could be generated by automatically fitting simulation results to theory. A similar strategy can help selecting suitable time steps, so that the diffusion of particles is correctly simulated (see Additional file 1).
A model of vesicular transport
HEINRICH and RAPOPORT proposed a model of vesicular transport [48]. Although the model describes a spatial process, it was formalized by means of differential equations.
The vesicular transport model refers to two cellular compartments which exchange membranebound soluble Nethylmaleimide–sensitive factor attachment protein receptors (SNAREs) and cargo proteins with the help of differently coated vesicles. These are budding from the compartments and move, driven by motorproteins, to the other compartment where they are fused. Thereby the coat of vesicles defines, which type of SNAREs (X or Y), cargo proteins and motorproteins are bound to them. Since different types of motorproteins move in different direction, SNARE X and Y accumulate in different compartments.
By modeling this process in a spatial regime, the problem starts, as in each particlebased approach, with the description of the species. Each species needs a radius and a diffusion coefficient. The MLForce simulator determines the diffusion based on the temperature, solvent viscosity and particle radius according to the STOKESEINSTEIN equation (see Additional file 1). The cell, compartments and the vesicles are represented as particles while the SNAREs are attributes of particles, which change during the budding and fusion processes. As KLANN et al. [2] has already shown through the translation of the original ODE model into a spatial agentbased model, this transport relies on a directed movement of vesicles. Therefore Klann et al. added a cytoskeleton structure to the cell on which the motorproteins moved along. In MLForce the direction is determined by forces which take the type of particles and their attributes into account. An external force field is introduced which is shaped like a dipole field with the compartments as poles. Based on the coat of the vesicles a force along this field is added to their movement, which lead them to one of the compartments. To test the sorting mechanism of this simple vesicular transport model we study a system with 2 types of SNAREs (X and Y) with the same amount. As in [48], we initialize the model: the first compartment is nine times bigger (V_{1}≈ 0.118 μm^{3}) than the second one (V_{2}≈ 0.013 μm^{3}) and contains 90% of the SNAREs (X and Y) (\(N_{1}^{X/Y}\) = 90000 and \(N_{2}^{X/Y}\) = 10000). Since the budding process of the vesicles depends on the volume of the compartments, they should reach an equal size and the SNAREs should be sorted.
As shown in Fig. 4 both compartments reach an equal volume and SNARE X accumulates in compartment 1 while SNARE Y accumulates in compartment 2 (not shown).
To show that the directed motion is essential for the sorting, as a control we run the simulation without the external force field and with compartments with equal initial size (V_{1/2}≈ 0.065 μm^{3}) and equally initially distributed SNAREs (\(N_{1/2}^{X/Y}\) = 50000). In a purely Brownian model, i.e., a model where all particles perform Brownian motion, no sorting takes place and both compartments are just shrinking as shown in Fig. 5. It can also be seen that the vesicle’s slowly diffusion around the compartments and there behavior is independent of their coating. In contrast the vesicles in the model with a directed movement either refuse with the compartment of origin or move to the other compartment based on their coating. As a result vesicles in this model fuse with a specific compartment after a short time and by this they are sorting the SNAREs. This simplified model illustrates how external forces can simulate the directed motion in a biological system.
The yeast model
A more complex model that showcases more of MLForce’s abilities is based on an illustrative yeast model [30], which in turn uses the model of the cell cycle proposed by TYSON [49]. In the original model each yeast cell can grow and contains five different proteins whose amounts oscillate periodically and trigger the fission of a yeast cell. Each cell has a mating type (P or M) which can change during the cell fission. Depending on the type the yeast releases pheromones (Mfactor and Pfactor, respectively) to inhibit the cell cycle of cells with opposite mating type. In addition, the type M cells secrete a protease called Sxa2 which inactivates the Pfactor pheromone that stems from the type P cells.
With the MLForce model, we want to check whether cells of identical mating types accumulate in some areas and inhibit the cell cycle of cells with opposite mating type in this area. In MLForce, each species of the system can be modeled as a particle, as shown in the upper part of Fig. 6. The problem of this approach are the different time scales of the dynamics. The cell cycle lasts about 120 min, while the diffusion of the proteins needs a time step in the sub nanosecond regime. In addition, for our question, a spatial simulation of the intracellular dynamics is not required. Therefore, we represent the pheromones and the cells as particles and the proteins involved in the cell cycle (as well as the cell cycle) as attributes of the cells, as shown in the lower part of Fig. 6. Whereas cells only move during the fission process in the xy plane, pheromones and Sxa2 diffuse in the extracellular medium and disappear when they reach the borders of the test volume.
As shown in Fig. 7 after 900 min there are twice as many Ptype cells as Mtype cells. Furthermore some of the Mtype cells are old (large) and have a highly inhibited cell cycle, which makes it likely that these cell will enter apoptosis before they can divide.
Attributed particles and arbitrary functions allow to model the intracellular dynamics nonspatially as first order reaction.
Nevertheless these nonspatial dynamic can influence the particle, here by triggering the division of a cell.
It is also possible for a particle to influence the attributes (nonspatial) of particle.
For example this is the case in the yeast model when a pheromone reacts with a cell (both particles) and inhibits the cell cycle of the cell. This combination of spatial and nonspatial behaviors enables the simulation of systems which include processes on different temporal scales.
Lipid rafts model
The role of lipid rafts in inducing and promoting receptor accumulation within the cell membrane and the recruitment and binding of proteins from the cytosol has been subject to several computational studies [50, 51]. A small model of receptor lipid raft interaction which is inspired by a submodel of Wnt signaling [52] shall illustrate the ability of MLForce to support dynamic nesting. The model consists of two proteins (LRP 5/6 and CK1 γ) which diffuse in the membrane and may enter or leave lipid rafts. In the model, the mobility of proteins is reduced in lipid rafts by a 10 times larger viscosity in the rafts compared to the rest of the membrane. Both proteins have different affinities to enter the lipid raft which is modeled by different energy barriers in the nesting reactions. LRP 5/6 has to pass a barrier of 12·10^{−21} J which was estimated from the distribution of the kinetic energy for demonstration purposes. It shows how the energy barrier at a bimolecular reaction can be used to slow down the reaction. CK1 γ can enter the lipid raft freely (E_{barrier} = 0 J). By choosing a larger energy barrier for LRP 5/6 while they are entering the lipid raft we model their lower affinity to the lipid raft and can show how the energy barrier at a bimolecular reaction can be used to slow down the reaction. To show how these functions influence the accumulation of the proteins in the lipid raft we run a 2D simulation. In this simple experiment, our model contains a single lipid raft which covers 25% of the surface. The 10 times larger viscosity of the lipid raft causes a 10 times lower diffusion coefficient in the raft compared to the rest of the membrane. The simulation starts with 200 LRP 5/6 and 200 CK1 γ particles outside the lipid raft. To enter the lipid raft, the particles need to overcome the energy barriers described above. After a particle entered a lipid raft they can leave it again without any energy barrier.
As expected, Fig. 8 shows that due to the slower diffusion in the lipid raft CK1 γ accumulates in the lipid raft. The amount of LRP 5/6 which accumulates in the lipid raft is lower due to the energy barrier which hampers receptors to enter the lipid raft.
The model shows how the repulsive force during a bimolecular reaction can lower the reaction rate. Also the dynamic nesting of particles can change their behavior (here the diffusion of the nested particle) during simulation. It is also possible to restrict reactions to only occur if particles are nested inside a specific type of particle, such as the phosphorylation of LRP 5/6 which is constrained to lipid rafts.
Discussion
The benchmark models which range from the multicellular yeast model to the subcellular lipid raft model have shown the usefulness of the realized features. MLForce expands upon the state of the art. The combination of compartmental dynamics and particle simulation based on forces is a unique feature of MLForce. This feature has been used in the lipid raft model (“Lipid rafts model” section). MLForce provides a consistent forcebased semantics for spatially resolved dynamic nesting with excluded volumes. Arbitrary functions and attributes have proven a necessity for a concise and computationally feasible realization of the yeast model in MLForce (“The yeast model” section) where they are used for the nonspatial stochastic simulation of the cell cycle. In MLForce, intraparticle dynamics can be simulated either spatially resolved by nested particles or in a nonspatial stochastic manner. Applying lower and higher spatial resolutions on demand facilitates modeling and simulating complex spatial models. The vesicle transport model (“A model of vesicular transport” section) relied on the possibility to let particles grow and to define global force functions that apply to all particles independently of other particles to simulate the directed movement of vesicles. The global force functions provide an additional means for the modeler for abstraction (in this case from the cytoskeleton structure on which the motorproteins move). Along all other dynamics this is seamlessly integrated into the forcebased semantics of MLForce.
Conclusions
The particlebased modeling and simulation approach MLForce combines excluded volumes and forces with support for dynamic nesting (compartmental dynamics) and the ability of constraining cellular dynamics by arbitrary attributes and functions. Whereas other simulation approaches have treated compartmental dynamics, such as particles entering or leaving a compartment, or compartmental fission or fusion, discretely, in MLForce those are simulated in a smooth, continuous manner. Thereby, structural dynamics are seamlessly integrated into the particlebased simulation governed by the LANGEVIN equation. However, this adds to the requirements the integrator has to face and consequently the induced calculation efforts.
MLForce utilizes a rulebased modeling approach. In combination with attributed species of arbitrary types and arbitrary functions that work on particles, their attributes and content, this contributes to the expressiveness and flexibility of MLForce. In particular, part of a spatial model can be executed in a nonspatial stochastic manner. Abstracting from spatial details on demand reduces the complexity of the model and the induced calculation effort and allows applying MLForce to a wider range of cell biological systems. The possibility to define force functions that act on all particles independently of others provides another valuable means for abstraction in MLForce.
Currently, MLForce is being used in a combined invitro and insilico study to analyze the impact of external electrical fields on cellular membranes. Future work will be aimed at enhancing the readability of MLForce models by a more declarative expression of rules. In addition, user support for selecting time steps and suitable parameters for force calculation needs to be provided. Also advanced numerical integration schemes that are exploited by other particlebased simulators and their impact on compartmental dynamics shall be analyzed.
Availability of data and materials
The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.
Abbreviations
 CK1 γ :

casein kinase 1 γ
 DSL:

domain specific language
 LRP5/6:

lowdensity lipoprotein receptorrelated protein 5/6
 ML:

multilevel
 RDME:

reaction diffusion master equation
 SNARE:

soluble Nethylmaleimide–sensitive factor attachment protein receptors
References
 1
Scott J. D., Pawson T.Cell Signaling in Space and Time: Where Proteins Come Together and When They’re Apart. Science. 2009; 326(5957):1220–4. https://doi.org/10.1126/science.1175668. Accessed 23 July 2019.
 2
Klann M, Koeppl H. Spatial Simulations in Systems Biology: From Molecules to Cells. Int J Mol Sci. 2012; 13(6):7798–827. https://doi.org/10.3390/ijms13067798.
 3
Bittig AT, Uhrmacher AM. Spatial modeling in cell biology at multiple levels. In: Proceedings of the 2010 Winter Simulation Conference. Baltimore: 2010. p. 608–19. https://doi.org/10.1109/WSC.2010.5679125.
 4
Ridgway D, Broderick G, Ellison MJ. Accommodating space, time and randomness in network simulation. Curr Opin Biotechnol. 2006; 17(5):493–8. https://doi.org/10.1016/j.copbio.2006.08.004.
 5
Takahashi K, TănaseNicola S, Wolde PRt. Spatiotemporal correlations can drastically change the response of a MAPK pathway. Proc Natl Acad Sci. 2010; 107(6):2473–8. https://doi.org/10.1073/pnas.0906885107.
 6
Cowan AE, Moraru II, Schaff JC, Slepchenko BM, Loew LM. Chapter 8  Spatial Modeling of Cell Signaling Networks In: Asthagiri AR, Arkin AP, editors. Methods in Cell Biology, Computational Methods in Cell Biology, vol. 110, pp. 192221. Cambridge: Academic Press: 2012. https://doi.org/10.1016/B9780123884039.000084. http://www.sciencedirect.com/science/article/pii/B9780123884039000084. Accessed 23 July 2019.
 7
Takahashi K, Arjunan SNV, Tomita M. Space in systems biology of signaling pathways  towards intracellular molecular crowding in silico. FEBS Lett. 2005; 579(8):1783–8. https://doi.org/10.1016/j.febslet.2005.01.072.
 8
Schöneberg J, Ullrich A, Noé F. Simulation tools for particlebased reactiondiffusion dynamics in continuous space. BMC Biophysics. 2014; 7(1):11. https://doi.org/10.1186/s1362801400115.
 9
Gillespie DT. Stochastic Simulation of Chemical Kinetics. Annu Rev Phys Chem. 2007; 58(1):35–55. https://doi.org/10.1146/annurev.physchem.58.032806.104637. Accessed 23 July 2019.
 10
Regev A, Panina EM, Silverman W, Cardelli L, Shapiro E. BioAmbients: an abstraction for biological compartments. Theor Comput Sci. 2004; 325(1):141–67. https://doi.org/10.1016/j.tcs.2004.03.061. Accessed 23 July 2019.
 11
Kerr RA, Bartol TM, Kaminsky B, Dittrich M, Chang JCJ, Baden SB, Sejnowski TJ, Stiles JR. Fast Monte Carlo Simulation Methods For BIOLOGICAL ReactionDiffusion Systems In Solution And On Surfaces. SIAM J Sci Comput Publ Soc Ind Appl Math. 2008; 30(6):3126. https://doi.org/10.1137/070692017. Accessed 23 July 2019.
 12
Andrews SS, Addy NJ, Brent R, Arkin AP. Detailed Simulations of Cell Biology with Smoldyn 2.1. PLoS Comput Biol. 2010; 6(3). https://doi.org/10.1371/journal.pcbi.1000705.
 13
Blinov ML, Schaff JC, Vasilescu D, Moraru II, Bloom JE, Loew LM. Compartmental and Spatial RuleBased Modeling with Virtual Cell. Biophys J. 2017; 113(7):1365–72. https://doi.org/10.1016/j.bpj.2017.08.022.
 14
Hattne J, Fange D, Elf J. Stochastic reactiondiffusion simulation with MesoRD. Bioinformatics. 2005; 21(12):2923–4. https://doi.org/10.1093/bioinformatics/bti431. Accessed 23 July 2019.
 15
Andrews SS. Smoldyn: particlebased simulation with rulebased modeling, improved molecular interaction and a library interface. Bioinformatics. 2017; 33(5):710–7. https://doi.org/10.1093/bioinformatics/btw700.
 16
Michalski P, Loew L. SpringSaLaD: A Spatial, ParticleBased Biochemical Simulation Platform with Excluded Volume. Biophys J. 2016; 110(3):523–29. https://doi.org/10.1016/j.bpj.2015.12.026.
 17
Schöneberg J, Noé F. ReaDDy  A Software for ParticleBased ReactionDiffusion Dynamics in Crowded Cellular Environments. PLoS ONE. 2013; 8(9):74261. https://doi.org/10.1371/journal.pone.0074261.
 18
Gruenert G, Ibrahim B, Lenser T, Lohel M, Hinze T, Dittrich P.Rulebased spatial modeling with diffusing, geometrically constrained molecules. BMC Bioinformatics. 2010; 11(1):307. https://doi.org/10.1186/1471210511307.
 19
Donovan RM, Tapia JJ, Sullivan DP, Faeder JR, Murphy RF, Dittrich M, Zuckerman DM. Unbiased Rare Event Sampling in Spatial Stochastic Systems Biology Models Using a Weighted Ensemble of Trajectories. PLOS Comput Biol. 2016; 12(2):1004611. https://doi.org/10.1371/journal.pcbi.1004611.
 20
Bittig AT, Uhrmacher AM. MLSpace: Hybrid Spatial Gillespie and Particle Simulation of Multilevel Rulebased Models in Cell Biology. IEEE/ACM Transactions on Computational Biology and Bioinformatics. 2016; PP(99):1–16. https://doi.org/10.1109/TCBB.2016.2598162.
 21
Flegg MB, Chapman SJ, Erban R. The tworegime method for optimizing stochastic reactiondiffusion simulations. J R Soc Interf. 2012; 9(70):859–68. https://doi.org/10.1098/rsif.2011.0574.
 22
Klann M, Ganguly A, Koeppl H. Hybrid spatial Gillespie and particle tracking simulation. Bioinformatics. 2012; 28(18):549–55. https://doi.org/10.1093/bioinformatics/bts384.
 23
Kim Y, Stolarska MA, Othmer HG. A hybrid model for tumor spheroid growth in vitro i: theoretical development and early results. Mathematical Models and Methods in Applied Sciences. 2007; 17(supp01):1773–98. https://doi.org/10.1142/S0218202507002479. Accessed 23 July 2019.
 24
Schaff JC, Gao F, Li Y, Novak IL, Slepchenko BM. Numerical Approach to Spatial DeterministicStochastic Models Arising in Cell Biology. PLOS Comput Biol. 2016; 12(12):1005236. https://doi.org/10.1371/journal.pcbi.1005236. Accessed 3 May 2018.
 25
Hall D, Minton AP. Macromolecular crowding: qualitative and semiquantitative successes, quantitative challenges. Biochim Biophys Acta (BBA)  Protein Proteomics. 2003; 1649(2):127–39. https://doi.org/10.1016/S15709639(03)001675.
 26
Neefjes J, Jongsma MML, Berlin I. Stop or Go? Endosome Positioning in the Establishment of Compartment Architecture, Dynamics, and Function. Trends Cell Biol. 2017; 27(8):580–94. https://doi.org/10.1016/j.tcb.2017.03.002.
 27
Peng D, Warnke T, Haack F, Uhrmacher AM. Reusing simulation experiment specifications to support developing models by successive extension. Simul Model Pract Theory. 2016; 68:33–53. https://doi.org/10.1016/j.simpat.2016.07.006.
 28
Köster T, Perumalla K, Uhrmacher A. Efficient Simulation of Nested Hollow Sphere Intersections: For Dynamically Nested Compartmental Models in Cell Biology. In: Proceedings of the 2017 ACM SIGSIM Conference on Principles of Advanced Discrete Simulation SIGSIMPADS ’17, pp. 173–183. New York: ACM: 2017. https://doi.org/10.1145/3064911.3064920.
 29
Faeder JR, Blinov ML, Hlavacek WS. RuleBased Modeling of Biochemical Systems with BioNetGen. In: Maly, IV, (ed.) Systems Biology. Methods in Molecular Biology pp. 113167. Totowa: Humana Press: 2009. https://doi.org/10.1007/9781597455251_5. Accessed 23 July 2019.
 30
Maus C, Rybacki S, Uhrmacher AM. Rulebased multilevel modeling of cell biological systems. BMC Syst Biol. 2011; 5(1):166. https://doi.org/10.1186/175205095166.
 31
John M, Lhoussaine C, Niehren J, Uhrmacher AM. The Attributed PiCalculus with Priorities. In: Transactions on Computational Systems Biology XII. Lecture Notes in Computer Science, pp. 13?76. Berlin: Springer: 2010. https://link.springer.com/chapter/10.1007/9783642117121_2.
 32
Pârvu O., Gilbert D., Heiner M., Liu F., Saunders N., Shaw S.SpatialTemporal Modelling and Analysis of Bacterial Colonies with Phase Variable Genes. ACM Trans Model Comput Simul. 2015; 25(2):13–11325. https://doi.org/10.1145/2742546.
 33
Wang Z, Butner JD, Kerketta R, Cristini V, Deisboeck TS. Simulating Cancer Growth with Multiscale AgentBased Modeling. Semin Cancer Biol. 2015; 30:70–8. https://doi.org/10.1016/j.semcancer.2014.04.001.
 34
Coffey WT, Kalmykov YP, Waldron JT. The Langevin Equation: With Applications to Stochastic Problems in Physics, Chemistry and Electrical Engineering, Revised. Singapore; River Edge, NJ: World Scientific Pub Co Inc; 2004. https://www.worldscientific.com/worldscibooks/10.1142/8195.
 35
Hertz H. Über die Berührung fester elastischer Körper. Journal für die reine und angewandte Mathematik (Crelle’s Journal). 1882; 92:156–71. https://doi.org/10.1515/crll.1882.92.156.
 36
Jeschke M, Ewald R. LargeScale Design Space Exploration of SSA. In: Computational Methods in Systems Biology, Lecture Notes in Computer Science, pp. 211230. Rostock: Springer: 2008. https://link.springer.com/chapter/10.1007/9783540885627_17.
 37
Johnson ME, Hummer G. FreePropagator Reweighting Integrator for SingleParticle Dynamics in ReactionDiffusion Models of Heterogeneous ProteinProtein Interaction Systems. Phys Rev X. 2014; 4(3):031037. https://doi.org/10.1103/PhysRevX.4.031037.
 38
von Smoluchowski M. Versuch einer mathematischen Theorie der Koagulationskinetik kolloider Lösungen. Z für Phys Chem. 1917; 92:129–68.
 39
Collins FC, Kimball GE. Diffusioncontrolled reaction rates. J Colloid Sci. 1949; 4(4):425–37. https://doi.org/10.1016/00958522(49)900239.
 40
Noyes RM. Effects of diffusion rates on chemical kinetics. Prog React Kinet. 1961; 1:129–60. https://ci.nii.ac.jp/naid/10016743949/. Accessed 12 June 2018.
 41
Van Deursen A, Klint P, Visser J. Domainspecific languages: An annotated bibliography. SIGPLAN Not. 2000; 35(6):26–36. https://doi.org/10.1145/352029.352035.
 42
Zunino R, Nikolic D, Priami C, Kahramanogulları O, Schiavinotto T. l: An Imperative DSL to Stochastically Simulate Biological Systems. In: Bodei, C, Ferrari, G, Priami, C, (eds.) Programming Languages with Applications to Biology and Security: Essays Dedicated to Pierpaolo Degano on the Occasion Of His 65th Birthday, pp. 354357. Cham: Springer: 2015. https://doi.org/10.1007/9783319255279_23.
 43
Himmelspach J, Uhrmacher AM. Plug’n Simulate. In: Simulation Symposium, 2007. ANSS ’07. 40th Annual: 2007. p. 137–43. https://doi.org/10.1109/ANSS.2007.34.
 44
Mannella R. Numerical Stochastic Integration for QuasiSymplectic Flows. SIAM J Sci Comput. 2006. https://doi.org/10.1137/040620965.
 45
The CImg Library  C++ Template Image Processing Toolkit. http://cimg.eu/. Accessed 11 Feb 2019.
 46
OpenGL  The Industry Standard for High Performance Graphics. https://www.opengl.org/. Accessed 11 Feb 2019.
 47
Yogurtcu O. N., Johnson M. E.Theory of bimolecular association dynamics in 2d for accurate model and experimental parameterization of binding rates. J Chem Phys. 2015; 143(8):084117. https://doi.org/10.1063/1.4929390. Accessed 23 July 2019.
 48
Heinrich R, Rapoport TA. Generation of nonidentical compartments in vesicular transport systems. J Cell Biol. 2005; 168(2):271–80. https://doi.org/10.1083/jcb.200409087.
 49
Tyson JJ. Modeling the cell division cycle: cdc2 and cyclin interactions. Proc Natl Acad Sci. 1991; 88(16):7328–32. https://doi.org/10.1073/pnas.88.16.7328.
 50
Nicolau DV, Burrage K, Parton RG, Hancock JF. Identifying Optimal Lipid Raft Characteristics Required To Promote Nanoscale ProteinProtein Interactions on the Plasma Membrane. Mol Cell Biol. 2006; 26(1):313–23. https://doi.org/10.1128/MCB.26.1.313323.2006.
 51
Haack F, Burrage K, Redmer R, Uhrmacher AM. Studying the Role of Lipid Rafts on Protein Receptor Bindings with Cellular Automata. IEEE/ACM Trans Comput Biol Bioinformatics. 2013; 10(3):760–70. https://doi.org/10.1109/TCBB.2013.40.
 52
Haack F, Lemcke H, Ewald R, Rharass T, Uhrmacher AM. Spatiotemporal Model of Endogenous ROS and RaftDependent WNT/BetaCatenin Signaling Driving Cell Fate Commitment in Human Neural Progenitor Cells. PLOS Comput Biol. 2015; 11(3):1004106. https://doi.org/10.1371/journal.pcbi.1004106.
Acknowledgements
Felix Hauptmann implemented the OpenGL visualizer.
Funding
Financial support (design of the study, analysis, and interpretation of simulations and in writing the manuscript) was provided by the Deutsche Forschungsgemeinschaft (DFG) research grant ’ESCeMMo’ (UH66/13) and the DFG collaborative research centre ’ELAINE’ (CRC 1270). We acknowledge financial support by Deutsche Forschungsgemeinschaft (DFG) and Universität Rostock within the funding programme Open Access Publishing.
Author information
Affiliations
Contributions
TK conceived the method. TK designed and implemented the Software. PH carried out the Case study. TK, PH and AMU wrote the manuscript. AMU supervised the project. All authors have read and approved the manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This 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.
About this article
Cite this article
Köster, T., Henning, P. & Uhrmacher, A.M. Potential based, spatial simulation of dynamically nested particles. BMC Bioinformatics 20, 607 (2019). https://doi.org/10.1186/s128590193092y
Received:
Accepted:
Published:
Keywords
 Space
 Simulation
 Nesting
 Force
 Multilevel
 Modeling
 Attributed