Advertisement

  • Loading metrics

Open Access

Peer-reviewed

Research Article

ReaDDy 2: Fast and flexible software framework for interacting-particle reaction dynamics

Abstract

Interacting-particle reaction dynamics (iPRD) combines the simulation of dynamical trajectories of interacting particles as in molecular dynamics (MD) simulations with reaction kinetics, in which particles appear, disappear, or change their type and interactions based on a set of reaction rules. This combination facilitates the simulation of reaction kinetics in crowded environments, involving complex molecular geometries such as polymers, and employing complex reaction mechanisms such as breaking and fusion of polymers. iPRD simulations are ideal to simulate the detailed spatiotemporal reaction mechanism in complex and dense environments, such as in signalling processes at cellular membranes, or in nano- to microscale chemical reactors. Here we introduce the iPRD software ReaDDy 2, which provides a Python interface in which the simulation environment, particle interactions and reaction rules can be conveniently defined and the simulation can be run, stored and analyzed. A C++ interface is available to enable deeper and more flexible interactions with the framework. The main computational work of ReaDDy 2 is done in hardware-specific simulation kernels. While the version introduced here provides single- and multi-threading CPU kernels, the architecture is ready to implement GPU and multi-node kernels. We demonstrate the efficiency and validity of ReaDDy 2 using several benchmark examples. ReaDDy 2 is available at the https://readdy.github.io/ website.

Author summary

Biological cells are not well-mixed reaction containers. Cellular signaling strongly depends on crowding, space exclusion, association and dissociation of proteins and other macromolecules. These are often confined to complex geometries and cell compartments. Understanding the mechanisms is challenging, as experiments can only achieve either high spatial resolution or high temporal resolution. Computer simulations on the other hand can capture these levels of detail. In particular, reaction-diffusion simulations can describe processes in a cell on a mesoscopic scale.

Reaction-diffusion simulations possess the highest level of detail if they are particle based. The aforementioned spatial effects can be captured by including forces that act between particles. We introduce the simulation tool ReaDDy 2 which implements a reaction-diffusion model in which diffusing particles can react and interact via forces. Further, macromolecules can be modeled by complex multi-particle structures. Compared to its predecessor ReaDDy 1, the present version 2 is significantly faster, thus accessing longer simulation timescales, has new functionalities such as reactive multi-particle groups, and provides an easy-to-use Python interface.

Citation: Hoffmann M, Fröhner C, Noé F (2019) ReaDDy 2: Fast and flexible software framework for interacting-particle reaction dynamics. PLoS Comput Biol 15(2): e1006830. https://doi.org/10.1371/journal.pcbi.1006830

Editor: Qiang Cui, Boston University, UNITED STATES

Received: August 9, 2018; Accepted: January 16, 2019; Published: February 28, 2019

Copyright: © 2019 Hoffmann et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: All files required to generate the figures can be found under https://github.com/readdy/readdy2_publication_data. The software and its implementation are publicly available under https://github.com/readdy/readdy.

Funding: We gratefully acknowledge funding from Deutsche Forschungsgemeinschaft (SFB 958 / Project A04, TRR 186 / Project A12, SFB 1114 / Project C03, SFB 740 / Project D7), Einstein Foundation Berlin (ECMath Project CH17) and European Research Council (ERC CoG 772230 “ScaleCell”). The authors are grateful to the Center for Theoretical Biological Physics (CTBP, supported by NSF PHY-1427654) at Rice University for hosting their sabbatical visit, during which part of this work was performed. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

This is a PLOS Computational Biology Software paper.

Introduction

The physiological response of biological cells to stimuli can be a many-stage process. A widely studied example is the MAPK pathway [1, 2]. Many of such signaling pathways incorporate G-protein coupled receptors (GPCR) [3] and cyclic adenosine monophosphate (cAMP) [4]. These are related to various diseases [5–7]. An extracellular stimulus can activate the membrane bound GPCRs and lead to localized synthesis of cAMP as second messengers. Their transport through the cell is diffusive, however due to the geometry of cellular compartments cAMP molecules are non-uniformly distributed [8, 9]. Their presence needs to be resolved in space and time to understand their function.

Particle-based reaction dynamics (PBRD) simulations [10–12] are amongst the most detailed approaches to model reaction kinetics computationally as they simulate each reactive molecule as a particle and therefore can be used as a tool to investigate systems of the aforementioned kind. Reactions can occur when reactive particles are in proximity, resembling the physical process. PBRD is suitable when the spatial distribution of molecules does not equilibrate rapidly and must therefore be resolved, and some reactants are locally scarce, such that their discrete number must be kept track of [13, 14]. There is a wide range of simulation tools for PBRD [15], including Smoldyn [16], MCell [17], Cell++ [18], eGFRD [2], mesoRD [19], spatiocyte [20], SpringSaLaD [21], and SRSim [22]. A simulation tool that takes the molecular structure into account is SDA [23–25]. A recent review of particle-based stochastic simulators can be found in [26]. Alternatively when the spatial resolution is of less importance, one can apply more efficient tools like Lattice Microbes [27–29] which generates realizations of the reaction-diffusion master equation (RDME) [30, 31]. In case of large copy numbers of particles it can make sense to think of them in terms of concentrations and build hybrid models [32].

PBRD simulations usually contain purely reactive particles that are not subject to interaction forces, e.g., to model space exclusion with repulsive interactions or clustering with attractive interactions. On the other hand, molecular dynamics (MD) simulations are designed to model particle dynamics including complex interactions between the particles or particles and an external field. The particles in MD simulations are often atoms or groups of atoms and higher-order structures such as molecules are represented by topology graphs that define the bonding structure between particles and thus, together with a MD force field, imply which pair, triplets and quadruplets or particles interact by means of bond, angle and torsion potentials. While reactive force fields [33–35] include reactivity on the chemistry scale, and soft matter MD simulation tools include breakable bonds [36, 37], current MD models and simulation packages do not incorporate generic particle reactions.

Interacting-particle reaction dynamics (iPRD) was introduced in [38] to combine the benefits of PBRD and MD simulations by modeling particle-based reaction dynamics while enabling full-blown interactions between particles as well as particles and the environment. Available simulation tools that are capable of special cases of iPRD simulations are, e.g., the MD packages LAMMPS [39] which is capable of forming and breaking bonds dynamically and ESPResSo [36, 37] which additionally has an implementation of catalytic reactions. In comparison to the iPRD simulator ReaDDy [38], these do not support full iPRD and are built and optimized for particle numbers that stay roughly constant. Comparing iPRD and PBRD, the interaction potentials in iPRD can be used to induce structure on mesoscopic length scales, e.g., volume-exclusion in crowded systems [38, 40], clustering of weakly interacting macromolecules [41], restriction of diffusing particles to arbitrarily-shaped membranes [15, 38, 42]. Furthermore it allows to study the large-scale structure of oligomers [43], polymers and membranes [44]. When not only considering interactions but also reactions, a wide range of reactive biochemical systems are in the scope of the model. For example, the reaction dynamics of photoreceptor proteins in crowded membranes [45] including cooperative effects of transmembrane protein oligomers [42] have been investigated. Another example is endocytosis, in which different proteins interact in very specific geometries [46, 47]. The simulation tool Cytosim [48] is another software package that can be used to investigate mesoscopic biochemical systems, specifically geared towards the simulation of the cytoskeleton.

The price of resolving these details is that the computation is dominated by computing particle-particle interaction forces. Although non-interacting particles can be propagated quickly by exploiting solutions of the diffusion equation [11, 49–51], interacting particles are propagated with small time-steps [52, 53], restricting the accessible simulation timescales whenever parts of the system are dense. As this computational expense is not entirely avoidable when the particle interactions present in iPRD are needed to model the process of interest realistically, it is important to have a simulation package that can fully exploit the computational resources.

Here we introduce the iPRD simulation framework ReaDDy 2, which is significantly faster, more flexible, and more conveniently usable than its predecessor ReaDDy 1 [38, 54]. Specifically, ReaDDy 2 includes the following new features:

  • Computational efficiency and flexibility: ReaDDy 2 defines computing kernels which perform the computationally most costly operations and are optimized for a given computing environment. The current version provides a single-CPU kernel that is four to ten times (depending on system size) faster than ReaDDy 1, and a multi-CPU kernel that scales with 80% efficiency to number of physical CPU cores for large particle systems (Sec. Performance). Kernels for GPUs or parallel multi-node kernels can be readily implemented with relatively little additional programming work (Sec. Design and implementation).
  • Python user interface: ReaDDy 2 can be installed via the conda package manager and used as a regular python package. The python interface provides the user with functionality to compose the simulation system, define particle interactions, reactions and parameters, as well as run, store and analyze simulations.
  • C++ user interface: ReaDDy 2 is mainly implemented in C++. Developers interested in extending the functionality of ReaDDy 2 in a way that interferes with the compute kernels, e.g., by adding new particle dynamics or reaction schemes, can do that via the C++ user interface.
  • Reversible reaction dynamics: ReaDDy 2 can treat reversible iPRD reactions by using steps that obey detailed balance, as described in [55] (iPRD-DB), and thus ensure correct thermodynamic behavior for such reactions (Sec. Reaction kinetics and detailed balance).
  • Topologies: We enable building complex multi-particle structures, such as polymers, by defining topology graphs (briefly: topologies, see Sec. Topologies). As in MD simulations, topologies are an efficient way to encode which bonded interactions (bond, angle and torsion terms) should act between groups of particles in the same topology. Note that particles in topologies can still be reactive. For example, it is possible to define reactions that involve breaking or fusing polymers (Sec. Topology reactions).
  • Potentials and boundaries: Furthermore, the range of by default supported interaction potentials has been broadened, now including harmonic repulsion, a harmonic interaction potential with a potential well, Lennard–Jones interaction, and screened electrostatics. The simulation volume can also be equipped with partially or fully periodic boundary conditions.

This paper summarizes the features of ReaDDy 2 and the demonstrates its efficiency and validity of ReaDDy 2 using several benchmarks and reactive particle systems. With few exceptions, we limit our description to the general features that are not likely to become outdated in future versions. Please see https://readdy.github.io/ for more details, tutorials and sample code.

interacting-particle reaction dynamics (iPRD)

The ReaDDy 2 simulation system consists of particles interacting by potentials and reactions (Fig 1) at a temperature T. Such a simulation system is confined to a box with repulsive or periodic boundaries. A boundary always has to be either periodic or be equipped with repulsive walls so that particles cannot diffuse away arbitrarily. To simulate iPRD dynamics in complex architectures, such as cellular membrane environments with specific shapes, additional potentials can be defined that confine the particle to a sub-volume of the simulation box (Sec. Potentials).

thumbnail
Download:
Fig 1. The simulation model.

(a) Potentials: Particles are subject to position-dependent external potentials, such as boundary potentials or external fields and interaction potentials involving two, three or four particles. As in MD force fields, bonded potentials are defined within particle groups called “topologies” whose bonding structure is defined by a connectivity graph. (b) Reactions: Most reactions are unimolecular or bimolecular particle reactions. Topology reactions act on the connectivity graphs and particle types and therefore change the particle bonding structure. (c) Simulation box: The simulation box with edge lengths ℓx, ℓy, and ℓz. It can optionally be periodic in a combination of x, y, and z directions, applying the minimum image convention.

https://doi.org/10.1371/journal.pcbi.1006830.g001

Interacting particle dynamics

ReaDDy 2 provides a developer interface to flexibly design models of how particle dynamics are propagated in time. The default model, however, is overdamped Langevin dynamics with isotropic diffusion as this is the most commonly used PBRD and iPRD model. In these dynamics a particle i moves according to the stochastic differential equation: (1) where contains the particle position at time t, Di(T) is its diffusion coefficient, kB is the Boltzmann constant, and T the system temperature. The particle moves according to the deterministic force fi and the stochastic velocity in which ξi are independent, Gaussian distributed random variables with moments where I is the identity matrix. The stochastic terms ξi and ξj are uncorrelated for particles i ≠ j. In ReaDDy 2 the default assumption is that the diffusion coefficients Di(T) are given for the simulation temperature T. Additionally, we offer the option to define diffusion coefficients for a reference temperature T0 = 293K and then generate the diffusion coefficients at the simulation temperature T by employing the Einstein-Smoluchowski model for particle diffusion in liquids [56, 57]:

This way, simulations at different temperatures are convenient while only having to specify one diffusion constant. Using this model, the dynamics are (2)

This means that the mobility is preserved if the temperature changes and (1) is recovered for T = T0.

The simplest integration scheme for (1), (2) is Euler-Maruyama, according to which the particle positions evolve as: (3) Where τ > 0 is a finite time step size and is a normally-distributed random variable. The diffusion constant Di effects the magnitude of the random displacement. The particles’ positions are loosely bound to a cuboid simulation box with edge lengths ℓx, ℓy, ℓz (Fig 1). If a boundary is non-periodic it is equipped with a repulsive wall given by the potential (4) acting on every component j of the single particle position xi, where k is the force constant, the cuboid in which there is no repulsion contribution of the potential, and d(⋅, Wj) ≔ inf{d(⋅, w): w ∈ Wj} the shortest distance to the set Wj. The cuboid can be larger than the simulation box in the periodic directions. In non-periodic directions there must be at least one repulsive wall for which this is not the case.

Due to the soft nature of the walls particles still can leave the simulation box in non-periodic directions. In that case they are no longer subject to pairwise interactions and bimolecular reactions however still are subject to the force of the wall pulling them back into the box.

Other types of dynamical models and other integration schemes can be implemented in ReaDDy 2 via its C++ interface. For example, non-overdamped dynamics, anisotropic diffusion [53, 58], hydrodynamic interactions [59] or employing the MD-GFRD scheme to make large steps for noninteracting particles will all affect the dynamical model and can be realized by writing suitable plugins.

Potentials

The deterministic forces are given by the gradient of a many-body potential energy U (Fig 1a):

The potentials are defined by the user. ReaDDy 2 provides a selection of standard potential terms, additional custom potentials can be defined via the C++ interface and then included into a Python simulation script.

External potentials only depend on the absolute position of each particle. They can be used, e.g., to form softly repulsive walls (4) and spheres, or to attach particles to a surface, for example to model membrane proteins. Furthermore the standard potential terms enable the user to simulate particles inside spheres and exclude particles from a spherical volume. The mentioned potential terms can also be combined to achieve more complex geometrical structures. Pair potentials generally depend on the particle distance and can be used, e.g., to model space exclusion at short distances.

A fundamental restriction of ReaDDy 2 interaction potentials is that they have a finite range and can therefore be cut off. This means that, e.g., full electrostatics is not supported but screened electrostatic interactions are implemented (Sec. Nontrivial bimolecular association kinetics at high concentrations). Additionally a harmonic repulsion potential, a weak interaction potential made out of three harmonic terms, and Lennard–Jones interaction are incorporated.

ReaDDy 2 has a special way of treating interaction potentials between bonded particles. Topologies define graphs of particles that are bonded and imply which particle pairs interact via bond constraints, which triples interact via angle constraints, and which quadruplets interact via a torsions potential. See Sec. Topologies for details.

Reactions

Reactions are discrete events, that can change particle types, add, and remove particles (Fig 1b). Each reaction is associated with a microscopic rate constant λ > 0 which has units of inverse time and represents the probability per unit time of the reaction occurring. The integration time-steps used in ReaDDy 2 should be significantly smaller than the inverse of the largest reaction rate, and we therefore compute discrete reaction probabilities by: (5)

In the software it is checked whether the time step τ is smaller than the inverse reaction rate up to a threshold factor of 10, otherwise a warning is displayed as discretization errors might become too large. In general, ReaDDy 2 reactions involve either one or two reactants. At any time step, a particle that is subject to an unary reaction will react with probability p(λ; τ). If there are two products, they are placed within a sphere of specified radius Ru around the educt’s position x0. This is achieved by randomly selecting an orientation , distance d ≤ Ru, and weights w1 ≥ 0,w2 ≥ 0, s.t. w1 + w2 = 1. The products are placed at x1 = x0 + dw1n and x2 = x0 − dw2n. Per default, w1 = w2 = 0.5 and the distances d are drawn such that the distribution is uniform with respect to the volume of the sphere. When it is necessary to produce new particles, we suggest to define a producing particle A and use the unary reaction A ⇀ A + B with corresponding placement weights w1 = 0, w2 = 1 so that the A particle stays at its position.

The basic binary reaction scheme is the Doi scheme [60, 61] in which a reactive complex is defined by two reactive particles being in a distance of Rb or less, where Rb is a parameter, e.g., see Fig 1b Fusion or Enzymatic reaction. The reactive complex then forms with probability p(λ; τ) while the particles are within distance.

Optionally ReaDDy 2 can simulate reversible reactions using the reversible iPRD-DB scheme developed in [55]. This scheme employs a Metropolis-Hastings algorithm that ensures the reversible reaction steps to be made according to thermodynamic equilibrium by accounting for the system’s energy in the educt and product states.

Topologies

Topologies are a way to group particles into superstructures. For example, large-scale molecules can be represented by a set of particles corresponding to molecular domains assembled into a topology. A topology also has a set of potential energy terms such as bond, angle, and torsion terms associated. The specific potential terms are implied by finding all paths of length two, three, and four in the topology connectivity graph. The sequence of particle types associated to these paths then is used to gather the potential term specifics, e.g., force constant, equilibrium length or angle, from a lookup table (Fig 1a).

Reactions are not only possible between particles, but also between a topology and a particle (Fig 1b) or two topologies. In order to define such reactions, one can register topology types and then specify the consequences of the reaction on the topology’s connectivity graph. We distinguish between global and local topology reactions.

Global topology reactions are triggered analogously to unary reactions, i.e., they can occur at any time with a fixed rate and probability as given in (5). Any edge in the graph can be removed and added. Moreover, any particle type as well as the topology type can be changed, which may result in significant changes in the potential energy. If the reaction causes the graph to split into two or more components, these components are subsequently treated as separate topologies that inherit the educt’s topology type and therefore also the topology reactions associated with it. Such a reaction is the topology analogue of a particle fission reaction.

A local topology reaction is triggered analogously to binary reactions with probability p(λ;τ) if the distance between two particles is smaller than the reaction radius. At least one of the two particles needs to be part of a topology with a specific type. The product of the reaction is then either yielded by the formation of an edge and/or a change of particle and topology types. In contrast to global reactions only certain changes to particle types and graphs can occur:

  • Two topologies can fuse, i.e., an additional edge is introduced between the vertices corresponding to the two particles that triggered the reaction.
  • A topology and a free particle can fuse by formation of an edge between the vertex of the topology’s particle and a newly introduced vertex for the free particle.
  • Two topologies can react in an enzymatic fashion, i.e., particle types of the triggering particles and topology types can be changed.
  • Two topologies and a free particle can react in an enzymatic fashion analogously.

In all of these cases the involved triggering particles’ types and topology types can be changed.

Simulation setup and boundary conditions

Once the potentials, the reactions (Fig 1a and 1b), and a temperature T have been defined, a corresponding simulation can be set up. A simulation box can be periodic or partially periodic, see Fig 1c. Periodicity in a certain direction means that with respect to that direction particle wrapping and the minimum image convention are applied. Non-periodic directions require a harmonically repelling wall as given in (4).

In order to define the initial condition, particles and particle complexes are added explicitly by specifying their 3D position and type. A simulation can now be started by providing a time step size τ and a number of integration steps.

Design and implementation

ReaDDy 2 is mainly written in C++ and has Python bindings making usage, configuration, and extension easy while still being able to provide high performance. To encourage usage and extension of the software, it is Open Source and licensed under the BSD-3 license. It therefore can not only be used in other Open Source projects without them requiring to have a similar license, but also in a commercial context.

Design

The software consists of three parts. The user-visible toplevel part is the python user interface, see Fig 2a. It is a language binding of the C++ user interface (Fig 2b) and has additional convenience functionality. The workflow consists out of three steps:

  1. The user is creating a readdy.ReactionDiffusionSystem, including information about temperature, simulation box size, periodicity, particle species, reactions, topologies, and physical units. Per default the configurational parameters are interpreted in a unit set well suited for cytosolic environments (lengths in nm, time in ns, and energy in kJ/mol), e.g., particles representing proteins in solution. The initial condition, i.e., the positions of particles, is not yet specified.
  2. The system can generate one or many instances of readdy.Simulation, in which particles and particle complexes can be added at certain positions. When instantiating the simulation object, a compute kernel needs to be selected, in order to specify how the simulation will be run (e.g., single-core or multi-core implementation). Additionally, observables to be monitored during the simulation are registered, e.g., particle positions, forces, or the total energy. A simulation is started by entering a time step size τ > 0 in units of time and a number of integration steps that the system should be propagated.
  3. When a simulation has been performed, the observables’ outputs have been recorded into a file. The file’s contents can be loaded again into a readdy.Trajectory object that can be used to produce trajectories compatible with the VMD molecular viewer [62].
thumbnail
Download:
Fig 2. The software structure.

(a) Python user interface: Provides a Python binding to the “C++ user interface” with some additional convenience functionality. The user creates a “readdy.ReactionDiffusionSystem” and defines particle species, reactions, and potentials. From a configured system, a “readdy.Simulation” object is generated, which can be used to run a simulation of the system given an initial placement of particles. (b) C++ core library: The core library serves as an adapter between the actual implementation of the algorithms in a compute kernel and the user interface. (c) Compute kernel implementation: Implements the compute kernel interface and contains the core simulation algorithms. Different compute kernel implementations support different hard- or software environments, such as serial and parallel CPU implementations. The compute kernel is chosen when the “readdy.Simulation” object is generated and then linked dynamically in order to provide optimal implementations for different computing environments under the same user interface.

https://doi.org/10.1371/journal.pcbi.1006830.g002

Running a simulation based on the readdy.Simulation object invokes a simulation loop. The default simulation loop is given in Alg. 1. Individual steps of the loop can be omitted. This enables the user to, e.g., perform pure PBRD simulations by skipping the calculation of forces. Performing a step in the algorithm leads to a call to the compute kernel interface, see Fig 2b. Depending on the selected compute kernel the call is then dispatched to the actual implementation. Compute kernel implementations (Fig 2c) are dynamically loaded at runtime from a plugin directory. This modularity allows ReaDDy 2 to run across many platforms although not every computing kernel may run on a given platform, such as a CUDA-enabled computing kernel. ReaDDy version 2 includes two iPRD computing kernels: a single threaded default computing kernel, and a dynamically-loaded shared-memory parallel kernel.

Algorithm 1: ReaDDy 2 default simulation loop. Each of the calls are dispatched to the compute kernel, see Fig 2. Furthermore, the user can decide to switch off certain calls in the simulation loop while configuring the simulation.

Initialize compute kernel;

if has output file then

 Write simulation setup;

end

Set up neighbor list;

Compute forces;

Evaluate observables;

while continue simulating do

 Call integrator;

 Update neighbor list;

 Perform reactions;

 Perform topology reactions;

 Update neighbor list;

 Calculate forces;

 Evaluate observables;

end

Tear down compute kernel;

The computing kernels contain implementations for the single steps of the simulation loop. Currently, integrator and reaction handler are exchangeable by user-written C++ extensions. Hence, there is flexibility considering what is actually performed during one step of the algorithm or even what kind of underlying model is applied.

In comparison to the predecessor ReaDDy 1, the software is a complete rewrite and extension. The functionality of the Brownian dynamics integrator has been preserved, however the reaction handlers can behave slightly differently. In particular, if during an integration step a reaction conflict occurs, i.e., there are at least two reaction events which involve the same educt particles, only one of these events can be processed. One possibility of choosing the to-be processed event is the so-called “UncontrolledApproximation”, which draws the next reaction event uniformly from all events and prunes conflicting events. Another possibility is drawing the next reaction event from all events weighted by their respective reaction probability. Since this approach is loosely based on the reaction order in the Gillespie SSA, this reaction handler is named “Gillespie” in ReaDDy 2.

With respect to the microscopic evaluation of a reaction event, the ReaDDy 1 implementation places product particles of fission type reactions at a fixed distance, which is handled more flexibly in the current implementation, see Sec. Reactions.

Performance

To benchmark ReaDDy 2, we use a reactive system with three particle species A, B, and C introduced in [63] with periodic boundaries instead of softly repelling ones. The simulation temperature is set to T = 293K and the diffusion coefficients are given by DA = 143.1 μm2s−1, DB = 71.6 μm2s−1, and DC = 68.82 μm2s−1. Particles of these types are subject to the two reactions A + B ⇀ C with microscopic association rate constant λon = 10−3 ns−1 and reaction radius R1 = 4.5 nm, and C ⇀ A + B with microscopic dissociation rate constant λoff = 5 × 10−5 ns−1 and dissociation radius R2 = R1. Particles are subject to an harmonic repulsion interaction potential which reads (6) where σ is the distance at which particles start to interact and κ = 10 kJ mol−1 nm−2 is the force constant. The interaction distance σ is defined as sum of radii associated to the particles’ types, in this case rA = 1.5 nm, rB = 3 nm, and rC = 3.12 nm. All particles are contained in a cubic box with periodic boundaries. The edge length is chosen such that the initial number density of all particles is ρtot = 3141 nm−3. This total density is distributed over the species, such that the initial density of A is ρA = ρtot/4, the initial density of B is ρB = ρtot/4, and the initial density of C is ρC = ρtot/2. For the chosen microscopic rates these densities roughly resemble the steady-state of the system. The performance is measured over a simulation timespan of 300 ns which is much shorter than the equilibration time of this system. Thus the overall number of particles does not vary significantly during measurement and we obtain the computation time at constant density.

In the following the benchmark results are presented. A comparison between the sequential reference compute kernel, the parallel implementation, and the previous Java-based ReaDDy 1 [38] is made with respect to their performance when varying the number of particles in the system keeping the density constant. Since the particle numbers fluctuate the comparison is based on the average computation time per particle and per integration step (Fig 3). The sequential kernel scales linearly with the number of particles, whereas the parallelized implementation comes with an overhead that depends on the number of threads. The previous Java-based implementation does not scale linearly for large particle numbers, probably owing to Java’s garbage collection. The parallel implementation starts to be more efficient than the sequential kernel given sufficiently many particles.

thumbnail
Download:
Fig 3. Performance comparison.

Average computation time per particle and integration step for the benchmark system of Sec. Performance using a machine with an Intel Core i7 6850K processor, i.e., six physical cores at 3.8GHz, and 32GB DDR4 RAM at 2.4GHz (dual channel). The number of particles is varied, but the particle density is kept constant. The sequential kernel (orange) has a constant per-particle CPU cost independent of the particle number. For large particle numbers, the parallel kernels are a certain factor faster (see scaling plot Fig 4). For small particle numbers of a few hundred the sequential kernel is more efficient. ReaDDy 2 is significantly faster and scales much better than the previous Java-based ReaDDy 1 [38].

https://doi.org/10.1371/journal.pcbi.1006830.g003

Fig 4 shows the strong scaling behavior of the parallel kernel, i.e. the speedup and efficiency for a fixed number of particles as a function of the used number of threads. For sufficiently large particle numbers, the kernel scales linear with the number of physical cores and an efficiency of around 80%. In hyperthreading mode, it then continues to scale linear with the number of virtual cores with an efficiency of about 55–60%.

thumbnail
Download:
Fig 4. Speedup and efficiency.

Parallel speedup and efficiency of the benchmark system of Sec. Performance as a function of the number of cores using the machine described in Fig 3. (a) Speedup with different numbers of cores compared to one core. Optimally one would like to have a speedup that behaves like the identity (black dashed line). (b) Efficiency is the speedup divided by the number of threads, i.e., how efficiently the available cores were used.

https://doi.org/10.1371/journal.pcbi.1006830.g004

The number of steps per day for a selection of particle numbers and kernel implementation is displayed in Table 1. For a system with 13, 000 particles and a time step size of τ = 1 ns (e.g., membrane proteins [63]), a total of 17ms simulation time per day can be collected on a six-core machine (Fig 3 for details). The current ReaDDy kernels are thus suited for the detailed simulation of processes in the millisecond- to second timescale, which include many processes in sensory signalling and signal transduction at cellular membranes.

thumbnail
Download:
Table 1. Number of steps per day for the benchmark system.

https://doi.org/10.1371/journal.pcbi.1006830.t001

Results

In the following, several aspects of the model applied in ReaDDy 2 are validated and demonstrated by considering different application scenarios and comparing the results to analytically obtained results, simulations from other packages, or literature data.

Reaction kinetics and detailed balance

We simulate the time evolution of particle concentrations of the benchmark system described in Sec. Performance. In contrast to the benchmarks, the considered system initially only contains A and B particles at equal numbers. It then relaxes to its equilibrium mixture of A, B, and C particles (Fig 5). Since the number of A and B molecules remain equal by construction, only the concentrations of A and C are shown.

thumbnail
Download:
Fig 5. Reaction kinetics and detailed balance.

Concentration time series of a the reaction-diffusion system introduced in Sec. Performance with the reversible reaction A + B ⇋ C. Compared are cases with and without harmonic repulsion (6). Additionally we compare two different reaction mechanisms, the Doi reaction scheme and the detailed balance (iPRD-DB) method for reversible reactions. (a) 30% volume occupation and no interaction potentials. (b) 30% volume occupation with harmonic repulsion between all particles. (c) 60% volume occupation and no interaction potentials. (d) 60% volume occupation with harmonic repulsion between all particles.

https://doi.org/10.1371/journal.pcbi.1006830.g005

In addition we compare the solutions with and without harmonic repulsion potentials (6) between all particles, as well as two different methods for executing the reactions: The Doi reaction scheme as described in Sec. Reactions and the detailed-balance reaction scheme iPRD–DB described in [55].

In contrast to Sec. Performance, we construct a macroscopic reference system with rate constants kon = 3.82 × 10−1 nm3s−1 and koff = 5 × 10−5s−1 resembling a cellular system. The microscopic reaction rate constants λon and λoff are then chosen with respect to the reference system taking interaction potentials between A and B into account. In particular, (7) (8) where is the accessible reaction volume, R the reaction radius, β the inverse thermal energy, and U the pair potential. The harmonic repulsion potential reduces Veff with respect to the volume of the reactive sphere. The expression (8) originates from an approximation for kon in a sufficiently well-mixed (i.e., reaction–limited) and sufficiently diluted system. The derivation can be found in [55] based on calculating the total association rate constant kon for an isolated pair of A and B particles. In this case one obtains λoff = 5 × 10−5 ns−1 for the microsopic dissociation rate constant. The microscopic association rate constant reads λon = 10−3 ns−1 for the noninteracting system and λon = 2.89 × 10−3 ns−1 for the interacting system. Note that for non-reversible binary reactions without interaction potentials the formula provided by [10, 64] describes the relation between λ and k for slow diffusion encounter. In the case of non-reversible binary reactions with interaction potentials and slow diffusion encounter such a relation can still be numerically computed [65].

Using the macroscopic rate constants kon and koff, a solution can be calculated for the mass-action reaction rate equations (RRE). This solution serves as a reference for the noninteracting system (no potentials), because the system parameters put the reaction kinetics in the mass-action limit.

In the noninteracting system, the ReaDDy solution and the RRE solution indeed agree (Fig 5a and 5c). In the case of interacting particles, see Fig 5b and 5d, an exact reference is unknown. We observe deviations from the RRE solution that become more pronounced with increasing particle densities. A difference between the two reaction schemes can also be seen. The Doi reaction scheme shows faster equilibration compared to RRE for increasing density, whereas the iPRD-DB scheme shows slower equilibration, as it has a chance to reject individual reaction events based on the change in potential energy. Thus an increased density leads to more rejected events, consistent with the physical intuition that equilibration in a dense system should be slowed down. Furthermore the equilibrated states differ depending on the reaction scheme, showing a dependence on the particle density. For denser systems the iPRD-DB scheme favors fewer A and B particles than the Doi scheme, consistent with the density-dependent equilibria described in [55].

Diffusion

Next we simulate and validate the diffusive behavior of noninteracting particle systems and the subdiffusive behavior of dense interacting particle systems. The simulation box contains particles with diffusion coefficient D0 and is equipped with softly repelling walls, in order to introduce finite size effects. The observations are carried out with and without interaction potential. In the case without interaction potential we compare with an analytical solution and the case of an interaction potential is compared to the literature.

Length x is given in units of σ, time t is given in units of σ2/D0, and energy is given in units of kBT. The cubic box has an edge length of ℓ ≈ 28σ.

The noninteracting particle simulation has a mean-squared displacement of particles in agreement with the analytic solution given by Fick’s law for diffusion in three dimensions (9) where x is the position of a particle and t is time (Fig 6). For long timescales t ≥ 101 transport is obstructed by walls, resulting in finite size saturation.

thumbnail
Download:
Fig 6. Diffusion in crowded environments.

Mean squared displacement as a function of time. Multiple particles are diffusing with intrinsic diffusion coefficient D0 in a cubic box with harmonically repelling walls. Triangles were obtained by using the Yukawa repulsion potential (10) between all particles. The dashed line represents an effective diffusion coefficient from the literature [66] for the same Yukawa repulsion potential.

https://doi.org/10.1371/journal.pcbi.1006830.g006

Fig 6 also shows that more complex transport can be modeled, as, e.g., found in crowded systems. Particles interact via the Yukawa potential [67] (10) where U0 = kBT is a repulsion energy, σ is the length scale, λ = 8 is the screening parameter, and rc = 2.5σ the cutoff radius.

The particle density is nσ3 = 0.6 with n being the number density. In such a particle system, the mean-squared displacement differs significantly from the analytical result for free diffusion after an initial time t ≥ 10−2 in which particles travel their mean free path length with diffusion constant D0. At intermediate timescales t ∈ [10−2, 10−1), particle transport is subdiffusive due to crowding. At long timescales, t ∈ [10−1, 101), the particles are again diffusive with an effective diffusion coefficient D that is reduced to reflect the effective mobility in the crowded systems. We compare this to an effective diffusion coefficient obtained by Brownian dynamics simulations from Löwen and Szamel [66] and find that they qualitatively agree. For large timescales t ≥ 101 finite size saturation can explicitly be observed as almost every particle has been repelled at least once by the boundaries.

To quantitatively compare the long-time effective diffusion coefficient D, we set up 1100 particles in a periodic box without repelling walls with the edge length chosen to give the densired density nσ3 = 0.6. The cutoff of the potential (10) is set to rc = 5σ, where U(rc) < 10−14 kBT. The particle suspension is equilibrated for at least teq ≥ 3 with a time-step size of τ = 10−5. We observe the mean squared displacement until tobs = 4.5 and measure the diffusion coefficient as the slope of a linear function for t ∈ [4, 4.5). We obtain D/D0 = 0.54 ± 0.01, which agrees with the reference value [66] Dref/D0 = 0.55 ± 0.01.

Thermodynamic equilibrium properties

We validate that ReaDDy 2’s integration of equations of motion yields the correct thermodynamics of a Lennard-Jones colloidal fluid in an (N, V, T) ensemble. To this end, we simulate a system of N particles confined to a periodic box with volume V at temperature T. The results and comparisons with other simulation frameworks and analytical results are shown in Table 2. The particles interact via the Lennard-Jones potential with ε being the depth of the potential well and σ the diameter of particles. The potential is cut off at rC = 4σ and shifted to avoid a discontinuity. The rescaled temperature is T* = kBTϵ−1 = 3. We perform simulations of the equilibrated Lennard-Jones system for 106 integration steps with rescaled time step size τ* = 10−4. Time units are σ2/D and are determined by the self-diffusion coefficient D of the particles. We measure the rescaled pressure P* = Pσ3ε−1 by estimating the virial term from forces acting in the system as described in [68]. Additionally we measure the rescaled potential energy per particle u* = UN−1ε−1. Both pressure and potential energy are calculated every 100th time step. This sampling gives rise to the mean and its error of the mean given for the ReaDDy 2 results in Table 2. Comparing HALMD [69] and ReaDDy 2, the latter shows larger energy and pressure in the third decimal place for the lower density ρ* = 0.3. For the higher density ρ* = 0.6 pressure differs in the first decimal place and energy in the second. This can be explained by ReaDDy 2 using an Euler scheme (3) to integrate motion of particles, which has a discretization error of first order in the time step size . On the other hand HALMD uses a Velocity-Verlet method [70], which has a discretization error of second order in the time step size .

thumbnail
Download:
Table 2. Thermodynamic equilibrium properties of a Lennard–Jones colloidal fluid in a (N, V, T) ensemble.

https://doi.org/10.1371/journal.pcbi.1006830.t002

Topology reactions

We illustrate ReaDDy 2’s ability to model complex reactions between multi-particle complexes, called “topology reactions”. We model polymers as linear chains of beads, held together by harmonic bonds and stiffened by harmonic angle potentials.

When considering just one worm-like chain with a certain amount of beads n, its equilibrium mean-squared end-to-end distance should behave like [73] (11) where lp = 4lk(kBT)−1 is the persistence length, Rmax = (n − 1)l the chain contour length, l the bond length, and k the force constant of the harmonic angles. In order to verify that the considered chain model obeys the mechanics of a worm-like chain, the theoretical mean-squared end-to-end distance (11) can be compared to observations from simulations, see Fig 7. For each fixed number of beads, an isolated chain was relaxed into an equilibrium state without performing topology reactions, yielding a squared end-to-end distance at the end of the simulation. This experiment was repeated 51 times. From the figure it can be observed that there is good agreement between the theoretical and measured mean-squared end-to-end distances.

thumbnail
Download:
Fig 7. Mean-squared end-to-end distance of worm-like chains.

The theoretical mean-squared end-to-end distance of worm-like chains as a function of number of beads (11) is compared to simulation data obtained from linear chains of beads as described in Sec. Topologies. Error bars depict errors over the mean from multiple measurements.

https://doi.org/10.1371/journal.pcbi.1006830.g007

In a system with many of these chains, we introduce two different particle types for the beads. Either they are head particles and located at the ends of a polymer chain or they are core particles and located between the head particles, as shown in Fig 8a and 8c in blue and orange, respectively.

thumbnail
Download:
Fig 8. Topology reactions example.

Illustrative simulation of polymer assembly/disassembly using topology reactions. (a) Sketch of the involved topology reactions. Association: When two ends of different topologies come closer than R, there is a rate λ1 that an edge is formed. Dissociation: The inverse of association with a rate λ2 and a randomly drawn edge that is removed. (b) The number of beads in a polymer 〈n(t)〉 over time averaged over 15 realizations. (c) Two representative particle configurations showing the initial state and the end state at time tbegin and tend, respectively.

https://doi.org/10.1371/journal.pcbi.1006830.g008

We impose two different topology reactions in the system with many chains (Fig 8a):

  1. Association: Two nearby head particles (distance ≤R) can connect with rate λ1. The topology is changed by adding an edge between the connected particles, resulting in the addition of one bond and two angle potentials. Additionally, the particle types of the two connected particles change from “head” to “core”.
  2. Dissociation: A chain with n particles can dissociate with microscopic rate nλ2, such that longer chains have a higher probability to dissociate than shorter chains. When a dissociation occurs, a random edge between two core particles is removed. The particle types of the respective core particles are changed to “head”. As a result, the graph decays into two connected components which subsequently are treated as autonomous topology instances.

The temporal evolution of the average length of polymer chains is depicted in Fig 8b. The simulation was performed 15 times with an initial configuration of 500 polymers containing four beads each. After sufficient time 〈n(t)〉 reaches an equilibrium value. Over the course of the simulation the polymers diffuse and form longer polymers. This can also be observed from the two snapshots shown in Fig 8c, depicting a representative initial configuration at tbegin and a representative configuration at the end of the simulation at time. In that case, there are polymers of many different lengths.

Nontrivial bimolecular association kinetics at high concentrations

This section studies a biologically inspired system with three macromolecules A, B, and C, that resemble, e.g., proteins in cytosol. The macromolecules A and B can form complexes C that also can dissociate back into their original components, i.e., we introduce reactions (12)

This form of interaction has been studied for proteins bovine serum albumin and hen egg white lysozyme in coarse-grained atomistic detail in [74] and for barnase and barstar in [75]. Here, we consider the case where the association reaction of (12) does not preserve volume, i.e., the complex C is more compact.

The presence of ions in aqueous solutions has effects on protein interactions [76], therefore we assume the reversibly associating macromolecules to be weakly charged and thus subject to the Debye-Hückel interaction potential [77] including an additional repulsion term (13) where s1, s2 ∈ {A, B, C}, q are partial charges associated with the macromolecules, e is the elementary charge, ε0 is the vacuum permittivity, εr is the relative permittivity of an aqueous solution, κ is the screening parameter that describes shielding due to ions in the solution, Ur is the repulsion energy, and is the sum of two particle radii. Here, we do not take hydration effects into account.

We investigate the equilibrium constant K = [A][B]/[C] for different number densities n = (NA + NB)/2 + NC. In case of a reversibly associating fluid described by the law of mass action, the equilibrium constant is given by K = koff/kon, where kon is the macroscopic association rate constant of (12) and koff the respective dissociation rate constant. In a well-mixed (i.e., reaction-limited) and sufficiently diluted system, kon can be approximated as in Sec. Reaction kinetics and detailed balance. However, for a diffusion-influenced process which we consider here, kon is typically understood as a harmonic mean of encounter and formation rates [78–81], i.e., . At low densities, only two-body interactions between A and B determine the on-rate constant, in this limit, kon can be evaluated numerically as a function of the microscopic association rate constant λon in the presence of the interaction potential, based on solving the Smoluchowski diffusion equation with a sink term that accounts for the volume reaction model, see [65]. Furthermore, in dense reversibly associating fluids, many-body interactions have an influence on kon, in particular due to competition for reactants, clustering, volume exclusion, and caging [81].

Thus, it is challenging to find a consistent analytical description over multiple orders of magnitudes in density. In contrast, we perform an empirical evaluation by simulations as shown in Fig 9. To this end, we set up 6 simulations for different n ∈ [2 × 101, 1.5 × 104] in a constant volume which then are allowed to relax into an equilibrium state subject to detailed-balance and yield a measurement K(n). The exact simulation parameters can be found in S1 Table. The reference value for the dilute case is given by , where koff = λoff and is a function of the microscopic association rate constant λon as well as the interaction potential (13) and is numerically computed as described in [65].

thumbnail
Download:
Fig 9. Equilibrium constant transition from dilute to dense systems.

The equilibrium constant K is obtained by simulation for different choices of the number of particles n = (NA + NB)/2 + NC which corresponds to a density due to constant volume of the simulation box and compared to an analytically obtained equilibrium constant of a dilute system (dashed line). The number of particles n remains constant during the course of a simulation. The shaded areas are standard deviations from the recorded data.

https://doi.org/10.1371/journal.pcbi.1006830.g009

We show that the reference value Kdilute is recovered by the simulation for low densities. For increasing densities more complex behavior can be observed. In particular, there is a drop in the value of K for which then is followed by a roughly stable regime up to n ≈ 5 × 103. For even higher densities, the equilibrium state is dominated by the complexes C likely due to finite size of the simulation volume. This drop in the equilibrium constant is in accord with Le Châtelier’s principle [82], i.e., the system prefers the state of lower free energy.

Availability and future directions

We have described the iPRD simulation framework ReaDDy 2 for combined particle interaction dynamics and reaction kinetics, which permits to conduct highly realistic simulations of signal transduction in crowded cellular environments or chemical nanoreactors with complex geometries. ReaDDy 2 follows up upon and significantly extends the simulation package ReaDDy 1. ReaDDy 2 is significantly faster than its predecessor, it can be easily installed as a Python conda package, and it can be flexibly used and reconfigured via its Python interface.

In comparison to molecular dynamics software packages, ReaDDy 2 does not include long range interactions. The software comes with a set of default interaction potentials. These include, e.g., harmonic repulsion which can model steric repulsion, Lennard–Jones interaction, and screened electrostatics which provide a way to model charged interaction at short ranges. Furthermore, ReaDDy allows for implementation of any short-ranged potential via a C++ interface. It is possible to implement and subsequently use hydration models which are short-ranged [83, 84] in the ReaDDy 2 framework. Hydrodynamic interactions are currently not included. They can be added by, e.g., providing an appropriate integrator which represents these interactions by a particle pairwise friction tensor [59].

Currently all pair potentials implemented in ReaDDy 2 are isotropic, however anisotropic interactions can be emulated by using particle complexes, in particular allowing for patchy particles. If the particles and interactions should be anisotropic themselves, a new computation kernel or appropriate integrator can be implemented into the framework via the C++ interface.

We have conducted a set of numerical studies, showing that ReaDDy 2 produces quantitatively accurate results where references from analytical solutions or other simulation packages were available, and physically meaningful results where reference solutions were not available.

For a quick and easy start into simulating and developing with ReaDDy 2 step by step tutorials, sample code, and further details are available online (https://readdy.github.io/). The software itself is Open Source and available under a permissive licence in order to enable a broad group of people to run simulations without forcing them to make their own work public.

ReaDDy 2 has been designed to be easily extensible. Planned extensions include simulation kernels for specialized hardware platforms, such as graphics processors and highly parallel HPC environments. Also planned is a MD-GFRD integrator [49] to speed up computations in dilute systems, and a particle-based membrane model as described in [44] that reproduces mechanical properties of cellular membranes.

In its current state, membranes can be modeled in terms of external forces, i.e., constraining particles onto two-dimensional surfaces. As these constraints only apply to selected particle types, it is possible to, e.g., grow polymers against a static membrane, where one end is anchored.

Supporting information

S1 Table. Parameters of density-dependent reaction kinetics.

This table contains parameters for the Debye-Hückel system in the results section.

https://doi.org/10.1371/journal.pcbi.1006830.s001

(PDF)

Acknowledgments

We are grateful for inspiring discussions with Manuel Dibak, Luigi Sbailò, Mohsen Sadeghi, Felix Höfling and Christof Schütte.

References

  1. 1. Xing J, Ginty DD, Greenberg ME. Coupling of the RAS-MAPK Pathway to Gene Activation by RSK2, a Growth Factor-Regulated CREB Kinase. Science. 1996-08;273(5277):959–963. pmid:8688081
  2. 2. Takahashi K, Tanase-Nicola S, ten Wolde PR. Spatio-temporal correlations can drastically change the response of a MAPK pathway. Proceedings of the National Academy of Sciences. 2010-02;107(6):2473–2478.
  3. 3. Trzaskowski B, Latek D, Yuan S, Ghoshdastider U, Debinski A, Filipek S. Action of Molecular Switches in GPCRs—Theoretical and Experimental Studies. Current Medicinal Chemistry. 2012-03;19(8):1090–1109. pmid:22300046
  4. 4. Beavo JA, Brunton LL. Cyclic nucleotide research—still expanding after half a century. Nature Reviews Molecular Cell Biology. 2002-09;3(9):710–717. pmid:12209131
  5. 5. Thathiah A, De Strooper B. The role of G protein-coupled receptors in the pathology of Alzheimer’s disease. Nature Reviews Neuroscience. 2011-02;12(2):73–87. pmid:21248787
  6. 6. Dragun D, Philippe A, Catar R, Hegner B. Autoimmune mediated G-protein receptor activation in cardiovascular and renal pathologies. Thrombosis and Haemostasis. 2009-03;101(4):643–648. pmid:19350106
  7. 7. Abramovitch R, Tavor E, Jacob-Hirsch J, Zeira E, Amariglio N, Pappo O, et al. A Pivotal Role of Cyclic AMP-Responsive Element Binding Protein in Tumor Progression. Cancer Research. 2004-02;64(4):1338–1346. pmid:14973073
  8. 8. Agarwal SR, Clancy CE, Harvey RD. Mechanisms Restricting Diffusion of Intracellular cAMP. Scientific Reports. 2016-04;6(1):19577. pmid:26795432
  9. 9. Houslay MD. Underpinning compartmentalised cAMP signalling through targeted cAMP breakdown. Trends in Biochemical Sciences. 2010-02;35(2):91–100. pmid:19864144
  10. 10. Erban R, Chapman SJ. Stochastic modelling of reaction–diffusion processes: algorithms for bimolecular reactions. Physical biology. 2009;6(4):046001. pmid:19700812
  11. 11. van Zon JS, Ten Wolde PR. Simulating biochemical networks at the particle level and in time and space: Green’s function reaction dynamics. Physical review letters. 2005;94(12):128103. pmid:15903966
  12. 12. Donev A, Yang Cy, Kim C. Efficient reactive Brownian dynamics. The Journal of Chemical Physics. 2018-01;148(3):034103. pmid:29352800
  13. 13. Hoffmann M, Schwarz US. Oscillations of Min-proteins in micropatterned environments: a three-dimensional particle-based stochastic simulation approach. Soft Matter. 2014;10(14):2388–2396. pmid:24622920
  14. 14. Bhalla US. Signaling in Small Subcellular Volumes. I. Stochastic and Diffusion Effects on Individual Pathways. Biophysical Journal. 2004-08;87(2):733–744. pmid:15298882
  15. 15. Schöneberg J, Ullrich A, Noé F. Simulation tools for particle-based reaction-diffusion dynamics in continuous space. BMC Biophysics. 2014;7(1):11. pmid:25737778
  16. 16. Andrews SS. Smoldyn: Particle-based simulation with rule-based modeling, improved molecular interaction and a library interface. Bioinformatics. 2017;33(5):710–717. pmid:28365760
  17. 17. Kerr RA, Bartol TM, Kaminsky B, Dittrich M, Chang JCJ, Baden SB, et al. Fast Monte Carlo simulation methods for biological reaction-diffusion systems in solution and on surfaces. SIAM journal on scientific computing. 2008;30(6):3126–3149. pmid:20151023
  18. 18. Sanford C, Yip ML, White C, Parkinson J. Cell++—simulating biochemical pathways. Bioinformatics. 2006;22(23):2918–2925. pmid:17038347
  19. 19. Hattne J, Fange D, Elf J. Stochastic reaction-diffusion simulation with MesoRD. Bioinformatics. 2005;21(12):2923–2924. pmid:15817692
  20. 20. Arjunan SNV, Tomita M. A new multicompartmental reaction-diffusion modeling method links transient membrane attachment of E. coli MinE to E-ring formation. Systems and synthetic biology. 2010;4(1):35–53. pmid:20012222
  21. 21. Michalski PJ, Loew LM. SpringSaLaD: a spatial, particle-based biochemical simulation platform with excluded volume. Biophysical journal. 2016;110(3):523–529. pmid:26840718
  22. 22. Gruenert G, Ibrahim B, Lenser T, Lohel M, Hinze T, Dittrich P. Rule-based spatial modeling with diffusing, geometrically constrained molecules. BMC bioinformatics. 2010;11(1):307. pmid:20529264
  23. 23. Martinez M, Bruce NJ, Romanowska J, Kokh DB, Ozboyaci M, Yu X, et al. SDA 7: A modular and parallel implementation of the simulation of diffusional association software. Journal of computational chemistry. 2015;36(21):1631–1645. pmid:26123630
  24. 24. Gabdoulline RR, Wade RC. Brownian dynamics simulation of protein–protein diffusional encounter. Methods. 1998;14(3):329–341. pmid:9571088
  25. 25. Gabdoulline RR, Wade RC. Simulation of the diffusional association of barnase and barstar. Biophysical journal. 1997;72(5):1917–1929. pmid:9129797
  26. 26. Andrews SS. Particle-Based Stochastic Simulators. Encyclopedia of Computational Neuroscience. 2018;.
  27. 27. Roberts E, Stone JE, Luthey-Schulten Z. Lattice microbes: High-performance stochastic simulation method for the reaction-diffusion master equation. Journal of computational chemistry. 2013;34(3):245–255. pmid:23007888
  28. 28. Roberts E, Magis A, Ortiz JO, Baumeister W, Luthey-Schulten Z. Noise contributions in an inducible genetic switch: a whole-cell simulation study. PLoS computational biology. 2011;7(3):e1002010. pmid:21423716
  29. 29. Roberts E, Stone JE, Sepúlveda L, Hwu WMW, Luthey-Schulten Z. Long time-scale simulations of in vivo diffusion using GPU hardware. In: Parallel & Distributed Processing, 2009. IPDPS 2009. IEEE International Symposium on. IEEE; 2009. p. 1–8.
    • 30. Isaacson SA. A convergent reaction-diffusion master equation. The Journal of chemical physics. 2013;139(5):054101. pmid:23927237
    • 31. Isaacson SA. The reaction-diffusion master equation as an asymptotic approximation of diffusion to a small target. SIAM Journal on Applied Mathematics. 2009;70(1):77–111.
    • 32. Franz B, Flegg MB, Chapman SJ, Erban R. Multiscale reaction-diffusion algorithms: PDE-assisted Brownian dynamics. SIAM Journal on Applied Mathematics. 2013;73(3):1224–1247.
    • 33. Hofmann DW, Kuleshova L, D’Aguanno B. A new reactive potential for the molecular dynamics simulation of liquid water. Chemical Physics Letters. 2007;448(1-3):138–143.
    • 34. Kamerlin SC, Warshel A. The empirical valence bond model: theory and applications. Wiley Interdisciplinary Reviews: Computational Molecular Science. 2011;1(1):30–45.
    • 35. Van Duin AC, Dasgupta S, Lorant F, Goddard WA. ReaxFF: a reactive force field for hydrocarbons. The Journal of Physical Chemistry A. 2001;105(41):9396–9409.
    • 36. Arnold A, Lenz O, Kesselheim S, Weeber R, Fahrenberger F, Roehm D, et al. ESPResSo 3.1—Molecular Dynamics Software for Coarse-Grained Models. In: Griebel M, Schweitzer MA, editors. Meshfree Methods for Partial Differential Equations VI. vol. 89 of Lecture Notes in Computational Science and Engineering. Springer; 2013. p. 1–23.
      • 37. Limbach HJ, Arnold A, Mann BA, Holm C. ESPResSo—An Extensible Simulation Package for Research on Soft Matter Systems. Comp Phys Comm. 2006;174(9):704–727.
      • 38. Schöneberg J, Noé F. ReaDDy-a software for particle-based reaction-diffusion dynamics in crowded cellular environments. PLoS One. 2013;8(9):e74261. pmid:24040218
      • 39. Plimpton S, Crozier P, Thompson A. LAMMPS-large-scale atomic/molecular massively parallel simulator. Sandia National Laboratories. 2007;18:43.
      • 40. Höfling F, Franosch T. Anomalous transport in the crowded world of biological cells. arXiv. 2013; p. 1–55.
        • 41. Ullrich A, Böhme MA, Schöneberg J, Depner H, Sigrist SJ, Noé F. Dynamical Organization of Syntaxin-1A at the Presynaptic Active Zone. PLOS Computational Biology. 2015-09;11(9):e1004407. pmid:26367029
        • 42. Gunkel M, Schöneberg J, Alkhaldi W, Irsen S, Noé F, Kaupp UB, et al. Higher-Order Architecture of Rhodopsin in Intact Photoreceptors and Its Implication for Phototransduction Kinetics. Structure. 2015;23(4):628–638. pmid:25728926
        • 43. Klein HC, Schwarz US. Studying protein assembly with reversible Brownian dynamics of patchy particles. The Journal of chemical physics. 2014;140(18):05B612_1.
        • 44. Sadeghi M, Weikl TR, Noé F. Particle-based membrane model for mesoscopic simulation of cellular dynamics. The Journal of Chemical Physics. 2018-01;148(4):044901. pmid:29390800
        • 45. Schöneberg J, Heck M, Hofmann KP, Noé F. Explicit Spatiotemporal Simulation of Receptor-G Protein Coupling in Rod Cell Disk Membranes. Biophysical Journal. 2014;107(5):1042–1053. pmid:25185540
        • 46. Posor Y, Eichhorn-Gruenig M, Puchkov D, Schöneberg J, Ullrich A, Lampe A, et al. Spatiotemporal control of endocytosis by phosphatidylinositol-3,4-bisphosphate. Nature. 2013-07;499(7457):233–237. pmid:23823722
        • 47. Schöneberg J, Lehmann M, Ullrich A, Posor Y, Lo Wt, Lichtner G, et al. Lipid-mediated PX-BAR domain recruitment couples local membrane constriction to endocytic vesicle fission. Nature Communications. 2017;8(May):15873. pmid:28627515
        • 48. Nedelec F, Foethke D. Collective Langevin dynamics of flexible cytoskeletal fibers. New Journal of Physics. 2007;9(11):427.
        • 49. Sbailò L, Noé F. An efficient multi-scale Green’s function reaction dynamics scheme. The Journal of Chemical Physics. 2017;147(18):184106. pmid:29141429
        • 50. van Zon JS, Ten Wolde PR. Green’s-function reaction dynamics: a particle-based approach for simulating biochemical networks in time and space. The Journal of chemical physics. 2005;123(23):234910. pmid:16392952
        • 51. Donev A, Bulatov VV, Oppelstrup T, Gilmer GH, Sadigh B, Kalos MH. A First-Passage Kinetic Monte Carlo algorithm for complex diffusion-reaction systems. Journal of Computational Physics. 2010;229(9):3214–3236.
        • 52. Vijaykumar A, Bolhuis PG, ten Wolde PR. Combining molecular dynamics with mesoscopic Green’s function reaction dynamics simulations. The Journal of Chemical Physics. 2015-12;143(21):214102. pmid:26646864
        • 53. Vijaykumar A, Ouldridge TE, ten Wolde PR, Bolhuis PG. Multiscale simulations of anisotropic particles combining molecular dynamics and Green’s function reaction dynamics. The Journal of Chemical Physics. 2017;146(11):114106. pmid:28330367
        • 54. Biedermann J, Ullrich A, Schöneberg J, Noé F. ReaDDyMM: Fast interacting particle reaction-diffusion simulations using graphical processing units. Biophysical journal. 2015;108(3):457–61. pmid:25650912
        • 55. Fröhner C, Noé F. Reversible Interacting-Particle Reaction Dynamics. The Journal of Physical Chemistry B. 2018;.
        • 56. Von Smoluchowski M. Zur kinetischen theorie der brownschen molekularbewegung und der suspensionen. Annalen der physik. 1906;326(14):756–780.
        • 57. Einstein A. Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen. Annalen der physik. 1905;322(8):549–560.
        • 58. Munk T, Höfling F, Frey E, Franosch T. Effective Perrin theory for the anisotropic diffusion of a strongly hindered rod. EPL (Europhysics Letters). 2009;85(3):30003.
        • 59. Ermak DL, McCammon J. Brownian dynamics with hydrodynamic interactions. The Journal of chemical physics. 1978;69(4):1352–1360.
        • 60. Doi M. Stochastic theory of diffusion-controlled reaction. Journal of Physics A: Mathematical and General. 1976;9(9):1479.
        • 61. Teramoto E, Shigesada N. Theory of bimolecular reaction processes in liquids. Progress of Theoretical Physics. 1967;37(1):29–51.
        • 62. Humphrey W, Dalke A, Schulten K. VMD: visual molecular dynamics. Journal of molecular graphics. 1996;14(1):33–38. pmid:8744570
        • 63. Schöneberg J, Noé F. ReaDDy–a software for particle-based reaction-diffusion dynamics in crowded cellular environments. PloS one. 2013;8(9):e74261. pmid:24040218
        • 64. Doi M. Theory of diffusion-controlled reactions between non-simple molecules. I. Chemical Physics. 1975;ll:107–113.
        • 65. Dibak M, Fröhner C, Höfling F, Noé F. Diffusion influenced reaction rates in the presence of a pair potential. in preparation;.
          • 66. Löwen H, Szamel G. Long-time self-diffusion coefficient in colloidal suspensions: theory versus simulation. Journal of Physics: Condensed Matter. 1993;5(15):2295.
          • 67. Yukawa H. On the interaction of elementary particles. I. Proceedings of the Physico-Mathematical Society of Japan 3rd Series. 1935;17:48–57.
          • 68. Allen MP, Tildesley DJ. Computer Simulation of Liquids. New York: Oxford University Press; 1987.
            • 69. Colberg PH, Höfling F. Highly accelerated simulations of glassy dynamics using GPUs: Caveats on limited floating-point precision. Computer Physics Communications. 2011;182(5):1120–1129.
            • 70. Swope WC, Andersen HC, Berens PH, Wilson KR. A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters. The Journal of Chemical Physics. 1982;76(1):637–649.
            • 71. Johnson JK, Zollweg JA, Gubbins KE. The Lennard-Jones equation of state revisited. Molecular Physics. 1993;78(3):591–618.
            • 72. Ayadim A, Oettel M, Amokrane S. Optimum free energy in the reference functional approach for the integral equations theory. Journal of Physics Condensed Matter. 2009;21(11). pmid:21693909
            • 73. Rubinstein M, Colby RH. Polymer physics. vol. 23. Oxford University Press New York; 2003.
              • 74. Mereghetti P, Martinez M, Wade RC. Long range Debye-Hückel correction for computation of grid-based electrostatic forces between biomacromolecules. BMC biophysics. 2014;7(1):4. pmid:25045516
              • 75. Plattner N, Doerr S, De Fabritiis G, Noé F. Complete protein–protein association kinetics in atomic detail revealed by molecular dynamics simulations and Markov modelling. Nature chemistry. 2017;9(10):1005. pmid:28937668
              • 76. Record MT, Anderson CF, Lohman TM. Thermodynamic analysis of ion effects on the binding and conformational equilibria of proteins and nucleic acids: the roles of ion association or release, screening, and ion effects on water activity. Quarterly reviews of biophysics. 1978;11(2):103–178. pmid:353875
              • 77. Debye P, Hückel E. Zur Theorie der Elektrolyte. I. Gefrierpunktserniedrigung und verwandte Erscheinungen. Physikalische Zeitschrift. 1923;24(9):185–206.
              • 78. Shoup D, Szabo A. Role of diffusion in ligand binding to macromolecules and cell-bound receptors. Biophysical Journal. 1982;40(1):33. pmid:7139033
              • 79. Szabo A, Schulten K, Schulten Z. First passage time approach to diffusion controlled reactions. The Journal of chemical physics. 1980;72(8):4350–4357.
              • 80. Grote RF, Hynes JT. The stable states picture of chemical reactions. II. Rate constants for condensed and gas phase reaction models. The Journal of Chemical Physics. 1980;73(6):2715–2732.
              • 81. Northrup SH, Hynes JT. Short range caging effects for reactions in solution. I. Reaction rate constants and short range caging picture. The Journal of Chemical Physics. 1979;71(2):871–883.
              • 82. Atkins P, De Paula J, Keeler J. Atkins’ physical chemistry. Oxford university press; 2018.
                • 83. Schneck E, Sedlmeier F, Netz RR. Hydration repulsion between biomembranes results from an interplay of dehydration and depolarization. Proceedings of the National Academy of Sciences. 2012;109(36):14405–14409.
                • 84. Halle B, Davidovic M. Biomolecular hydration: from water dynamics to hydrodynamics. Proceedings of the National Academy of Sciences. 2003;100(21):12135–12140.