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

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.


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 Gprotein coupled receptors (GPCR) [3] and cyclic adenosine monophosphate (cAMP) [4]. These are related to various diseases [5][6][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][11][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][24][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][28][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][34][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][50][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).

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: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 2D i ðTÞ where x i ðtÞ 2 R 3 contains the particle position at time t, D i (T) is its diffusion coefficient, k B is the Boltzmann constant, and T the system temperature. The particle moves according to the deterministic force f i and the stochastic velocity ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 2D i ðTÞ p ξ i 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 6 ¼ j.
In ReaDDy 2 the default assumption is that the diffusion coefficients D i (T) are given for the simulation temperature T. Additionally, we offer the option to define diffusion coefficients for a reference temperature T 0 = 293K and then generate the diffusion coefficients at the simulation temperature T by employing the Einstein-Smoluchowski model for particle diffusion in 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. 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 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi This means that the mobility is preserved if the temperature changes and (1) is recovered for T = T 0 .
The simplest integration scheme for (1), (2) is Euler-Maruyama, according to which the particle positions evolve as: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Where τ > 0 is a finite time step size and η t � N ð0; 1Þ is a normally-distributed random variable. The diffusion constant D i 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 acting on every component j of the single particle position x i , where k is the force constant, x ðjÞ origin þ x ðjÞ extent � the cuboid in which there is no repulsion contribution of the potential, and d(�, W j ) ≔ inf{d(�, w): w 2 W j } the shortest distance to the set W j . 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: 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 R u around the educt's position x 0 . This is achieved by randomly selecting an orientation n 2 R 3 , distance d � R u , and weights w 1 � 0,w 2 � 0, s.t. w 1 + w 2 = 1. The products are placed at x 1 = x 0 + dw 1 n and x 2 = x 0 − dw 2 n. Per default, w 1 = w 2 = 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 w 1 = 0, w 2 = 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 R b or less, where R b 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 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. 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].
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. 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 tobe 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 D A = 143.1 μm 2 s −1 , D B = 71.6 μm 2 s −1 , and D C = 68.82 μm 2 s −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 R 1 = 4.5 nm, and C * A + B with microscopic dissociation rate constant λ off = 5 × 10 −5 ns −1 and dissociation radius R 2 = R 1 . Particles are subject to an harmonic repulsion interaction potential which reads 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 r A = 1.5 nm, r B = 3 nm, and r C = 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 Javabased 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. 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%.
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 sixcore 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.

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 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 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.
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 k on = 3.82 × 10 −1 nm 3 s −1 and k off = 5 × 10 −5 s −1 resembling a cellular system. The  Number of time steps per day for benchmark system of Sec. Performance using the machine described in Fig 3. In case of the parallelized implementation the peak performance with respect to the number of threads is shown.
https://doi.org/10.1371/journal.pcbi.1006830.t001 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, where V eff ¼ R R 0 exp ðÀ bUÞ4pr 2 dr is the accessible reaction volume, R the reaction radius, β the inverse thermal energy, and U the pair potential. The harmonic repulsion potential reduces V eff with respect to the volume of the reactive sphere. The expression (8) originates from an approximation for k on 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 k on 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 k on and k off , a solution can be calculated for the massaction 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 massaction 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 D 0 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 /D 0 , and energy is given in units of k B T. 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 where x is the position of a particle and t is time (Fig 6). For long timescales t � 10 1 transport is obstructed by walls, resulting in finite size saturation. 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] UðrÞ ¼ where U 0 = k B T is a repulsion energy, σ is the length scale, λ = 8 is the screening parameter, and r c = 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 D 0 . At intermediate timescales t 2 [10 −2 , 10 −1 ), particle transport is subdiffusive due to crowding. At long timescales, t 2 [10 −1 , 10 1 ), 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 � 10 1 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 r c = 5σ, where U(r c ) < 10 −14 k B T. The particle suspension is equilibrated for at least t eq � 3 with a time-step size of τ = 10 −5 . We observe the mean squared displacement until t obs = 4.5 and measure the diffusion coefficient as the slope of a linear function for t 2 [4, 4.5). We obtain D/D 0 = 0.54 ± 0.01, which agrees with the reference value [66] D ref /D 0 = 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  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 r C = 4σ and shifted to avoid a discontinuity. The rescaled temperature is T � = k B T� −1 = 3. We perform simulations of the equilibrated Lennard-Jones system for 10 6 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 OðtÞ. On the other hand HALMD uses a Velocity-Verlet method [70], which has a discretization error of second order in the time step size Oðt 2 Þ.

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] where l p = 4lk(k B T) −1 is the persistence length, R max = (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.
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.
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 hn(t)i 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 t begin 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 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 where s 1 , s 2 2 {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, U r is the repulsion energy, and s s 1 s 2 ¼ r s 1 þ r s 2 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 = (N A + N B )/2 + N C . In case of a reversibly associating fluid described by the law of mass action, the equilibrium constant is given by K = k off /k on , where k on is the macroscopic association rate constant of (12) and k off the respective dissociation rate constant. In a well-mixed (i.e., reaction-limited) and sufficiently diluted system, k on can be approximated as in Sec. Reaction kinetics and detailed balance. However, for a diffusion-influenced process which we consider here, k on is typically understood as a harmonic mean of encounter and formation rates [78][79][80][81], i.e., k À 1 on ¼ k À 1 enc þ k À 1 form . At low densities, only two-body interactions between A and B determine the on-rate constant, in this limit, k on 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, manybody interactions have an influence on k on , 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 [2 × 10 1 , 1.5 × 10 4 ] in a constant volume which then are allowed to relax into an equilibrium state subject to detailedbalance 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 K dilute ¼ k off =k dilute on , where k off = λ off and k dilute on 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].
We show that the reference value K dilute 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 n ≳ 10 2 which then is followed by a roughly stable regime up to n � 5 × 10 3 . 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. (PDF)