ReaDDy - A particle-based
reaction-diffusion simulator

Michaelis-Menten - well mixed

In [2]:
%matplotlib inline
import os
import numpy as np
import readdy
In [3]:
readdy.__version__
Out[3]:
'v1.0.1-0'

Setup ReaDDy system

The intrinsic forward rate $\lambda_\mathrm{f}$ is obtained from the macroscopic forward rate $k_\mathrm{f} = 5.88 \,\mathrm{\mu M}^{-1} \mathrm{s}^{-1} = 0.98 \times 10^{-2} \mathrm{\mu m}^{3} \mathrm{s}^{-1}$ by application of

$$ k_\mathrm{f} = 4\pi (D_E + D_S) R \left[1 - \frac{\tanh(\kappa R)}{\kappa R}\right] $$

where

$$ \kappa = \sqrt{\frac{\lambda_\mathrm{f}}{D_E + D_S}} $$

as described in [1]. Note that this will only work in the well mixed regime. The equation cannot be solved explicitly for $\kappa$. Hence, the intrinsic rate was numerically approximated by solving a minimization problem.

[1] R. Erban and J. Chapman, “Stochastic modelling of reaction-diffusion processes: algorithms for bimolecular reactions.,” Phys. Biol., vol. 6, no. 4, p. 46001, Jan. 2009.

In [4]:
system = readdy.ReactionDiffusionSystem(
    box_size=[0.3, 0.3, 0.3],
    unit_system={"length_unit": "micrometer", "time_unit": "second"}
)

system.add_species("E", diffusion_constant=10.)
system.add_species("S", diffusion_constant=10.)
system.add_species("ES", diffusion_constant=10.)
system.add_species("P", diffusion_constant=10.)

system.reactions.add("fwd: E +(0.03) S -> ES", rate=86.78638438)
system.reactions.add("back: ES -> E +(0.03) S", rate=1.)
system.reactions.add("prod: ES -> E +(0.03) P", rate=1.)

Simulate the system

In [5]:
simulation = system.simulation(kernel="CPU")
simulation.output_file = "out.h5"
simulation.reaction_handler = "UncontrolledApproximation"

n_particles_e = 909
n_particles_s = 9091
edge_length = system.box_size[0]
initial_positions_e = np.random.random(size=(n_particles_e, 3)) * edge_length - .5*edge_length
initial_positions_s = np.random.random(size=(n_particles_s, 3)) * edge_length - .5*edge_length
simulation.add_particles("E", initial_positions_e)
simulation.add_particles("S", initial_positions_s)

simulation.observe.number_of_particles(stride=1, types=["E", "S", "ES", "P"])
In [6]:
if os.path.exists(simulation.output_file):
    os.remove(simulation.output_file)

dt = 8e-4
n_steps = int(10. / dt)

simulation.run(n_steps=n_steps, timestep=dt)
Configured kernel context with:
--------------------------------
 - kBT = 2.43614
 - periodic b.c. = (true, true, true)
 - box size = (0.3, 0.3, 0.3)
 - particle types:
     *  particle type "P" with D=10
     *  particle type "ES" with D=10
     *  particle type "E" with D=10
     *  particle type "S" with D=10
 - unimolecular reactions:
     * Fission ES -> E + S with a rate of 1, a product distance of 0.03, and weights 0.5 and 0.5
     * Fission ES -> E + P with a rate of 1, a product distance of 0.03, and weights 0.5 and 0.5
 - bimolecular reactions:
     * Fusion E + S -> ES with a rate of 86.7864, an educt distance of 0.03, and weights 0.5 and 0.5


Obtain simulation output

In [7]:
traj = readdy.Trajectory(simulation.output_file)
time, counts = traj.read_observable_number_of_particles()
counts = counts / system.box_volume.magnitude
time = time * dt

Analytical solution (integrated ODE)

In [8]:
from scipy.integrate import odeint

def f(x, t0, kf, kr, kcat):
    """
    x: state vector with concentrations of E, S, ES, P
    """
    return np.array([
        -kf * x[0] * x[1] + (kr+kcat)*x[2],
        -kf * x[0] * x[1] + kr * x[2],
        kf * x[0] * x[1]- (kr+kcat)*x[2],
        kcat*x[2]
    ])

init_state = np.array([n_particles_e, n_particles_s, 0., 0.]) / system.box_volume.magnitude
ode_time = np.linspace(0.,10,100000)
ode_result = odeint(f, y0=init_state, t=ode_time, args=(0.98e-2, 1., 1.))
In [9]:
import matplotlib.pyplot as plt

stride = 100
plt.plot(time[::stride], counts[:,0][::stride], "-", label='E')
plt.plot(time[::stride], counts[:,1][::stride], "-", label='S')
plt.plot(time[::stride], counts[:,2][::stride], "-", label='ES')
plt.plot(time[::stride], counts[:,3][::stride], "-", label='P')
plt.plot(ode_time, ode_result[:,0], "--", color="grey", label="LMA solution")
plt.plot(ode_time, ode_result[:,1], "--", color="grey")
plt.plot(ode_time, ode_result[:,2], "--", color="grey")
plt.plot(ode_time, ode_result[:,3], "--", color="grey")
plt.legend()
plt.xlabel(r"time in $\mathrm{s}$")
plt.ylabel(r"concentration in $\mathrm{\mu m}^{-3}$")
plt.show()

Performance data

In [10]:
perf = simulation._simulation.performance_root()
In [11]:
print(perf)
simulation: time 104.921 s, count 1
	initNeighborList: time 0.003026 s, count 1
		set_up: time 0.002545 s, count 1
			setUpCellNeighbors: time 0.001843 s, count 1
			setUpBins: time 0.000693 s, count 1
				allocate: time 1.3e-05 s, count 1
				fillBins parallel: time 0.000673 s, count 1
		update: time 0.000469 s, count 1
			setUpBins: time 0.000465 s, count 1
				allocate: time 1.1e-05 s, count 1
				fillBins parallel: time 0.000446 s, count 1
	forces: time 0.000563 s, count 12501
	integrator: time 5.49103 s, count 12500
	neighborList: time 3.85435 s, count 25000
		update: time 3.80827 s, count 25000
			setUpBins: time 3.76443 s, count 25000
				allocate: time 0.120096 s, count 25000
				fillBins parallel: time 3.56029 s, count 25000
	reactionScheduler: time 93.8283 s, count 12500
	evaluateTopologyReactions: time 7.8e-05 s, count 12500
	clearNeighborList: time 0 s, count 1

© Copyright 2018 CMB group