ReaDDy

What it simulates

ReaDDy simulates particles that diffuse and react in 3D space. Each particle is a point. It moves by overdamped Brownian dynamics. A reaction happens between two particles when they come closer than a reaction radius. Reaction rates are microscopic rates with units of inverse time. The per-step probability is p = 1 - exp(-rate * dt).

You describe species, diffusion constants and reactions in Python, place particles, and run. ReaDDy can add potentials, polymer bonds and topologies. The example on this page uses none of those. It runs one reversible reaction, A + B <-> C, in a cubic box of edge 10. It returns the count of A, B and C over time.

The key assumption is the spatial resolution. Particles have no volume unless you add an interaction potential. The fusion reaction is evaluated on particle pairs that are within the radius.

When to use it

Use ReaDDy when space matters in a reaction network:

• Does crowding or slow diffusion lower the binding rate below the well-mixed value?
• How many complexes form in a small compartment with a few hundred molecules?
• What do reactions do on polymers, membranes or other bonded structures?

GillesPy2 on biosimhub fits better for well-mixed stochastic models, where space does not matter and runs are much cheaper. Cytosim fits filament mechanics.

Install it yourself

The README gives conda as the preferred route, on Mac and Linux only. The dry run micromamba create -n t --dry-run -c conda-forge -c bioconda readdy=2.0.14 resolves readdy 2.0.14 on linux-64.

conda create -n readdy -c conda-forge readdy=2.0.14
conda activate readdy

The recipe's image uses the same package and runs the single-threaded kernel. The minimal Python script follows the ReaDDy docs for the system and the simulation:

import numpy as np
import readdy

system = readdy.ReactionDiffusionSystem(box_size=[10., 10., 10.], unit_system=None)
for name in ("A", "B", "C"):
    system.add_species(name, diffusion_constant=1.)
system.reactions.add("fusion: A +(1) B -> C", rate=1.)
system.reactions.add("fission: C -> A +(1) B", rate=0.1)

sim = system.simulation(kernel="SingleCPU")
sim.output_file = "traj.h5"
sim.add_particles("A", np.random.uniform(-5, 5, size=(100, 3)))
sim.add_particles("B", np.random.uniform(-5, 5, size=(100, 3)))
sim.observe.number_of_particles(stride=25, types=["A", "B", "C"], callback=print)
sim.run(5000, 0.01)

With unit_system=None the system is unitless and kbt is 1. The default unit system is nanometers, nanoseconds and kJ/mol. In the recipe's image the box has periodic_boundary_conditions of [True, True, True].

Reading the output

The recipe's helper script run_readdy.py wraps the script above. It builds the system from the command-line values, scatters the particles uniformly in the box, and records counts every steps // 200 steps. With 5000 steps that is every 25 steps, which gives 201 rows in particle_counts.csv (4052 bytes).

The file has the columns step,time,A,B,C. It starts at 0,0.0,100,100,0. At step 25 (time 0.25) it reads 91, 91 and 9. By time 2.5 the line is 58, 58, 42. A and B always match because each reaction removes one of each and starts with 100 of each.

readdy_log.txt begins with the helper's line final counts A,B,C: [40, 40, 60]. ReaDDy's own progress lines follow, one per 100 steps, in the form Simulation progress: 100 / 5000 steps. At time 50 the system holds 60 complexes.

The helper also writes a trajectory file to /tmp inside the image.

Your own runs

Copy the docker run line from the recipe and change values. All runs used 5000 steps of 0.01.

1. Faster fusion. Set --fusion_rate 10. The final counts changed from 40, 40, 60 to 17, 17, 83. More A and B ended up bound.
2. Faster fission. Set --fission_rate 1. The final counts were 78, 78, 22. C splits apart almost as fast as it forms.
3. Faster diffusion. Set --diffusion 10. At time 1.5 the example run has 30 complexes. The default seed with --diffusion 10 gave 43. The final counts were 45, 45, 55 against 40, 40, 60.

A different --seed 7 gave the same final counts of 40, 40, 60 as seed 42, but a different early path. At time 0.25 it had 7 complexes, against 9 for seed 42.

Pitfalls

• The seed is partial. The --seed value only sets the initial positions. ReaDDy's own random generator is not seedable, so repeat runs of the same seed can differ.
• Rates are per unit time of your system. With unit_system=None the numbers have no units. Keep diffusion, radius, rate and time step in one consistent set.
• Time step against radius. A particle should move a small fraction of the reaction radius per step. At D = 1 and dt = 0.01 the root-mean-square step is about 0.14 per axis (sqrt(2 D dt)), against a radius of 1.
• Box and counts set the concentration. 100 particles in a box of edge 10 is 0.1 per unit volume. Doubling the edge divides the concentration by 8.
• Linux amd64 only. The recipe's image runs as linux/amd64. On Apple silicon it runs under emulation and is slower.

Further reading

• Documentation: https://readdy.github.io/
• System configuration: https://readdy.github.io/system.html
• Simulation and observables: https://readdy.github.io/simulation.html
• Hoffmann M, Frohner C, Noe F (2019). ReaDDy 2: Fast and flexible software framework for interacting-particle reaction dynamics. PLoS Comput Biol 15:e1006830. https://doi.org/10.1371/journal.pcbi.1006830