# Bimolecular reaction - well mixed

In [2]:
import os
import numpy as np
import matplotlib.pyplot as plt

In [3]:
readdy.__version__

Out[3]:
'v1.0.1-0'

In [4]:
system = readdy.ReactionDiffusionSystem([20.,20.,20.], temperature=300.*readdy.units.kelvin)

lambda_on = 1.


Simulate the system

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

n_particles = 2000
initial_positions_a = np.random.random(size=(n_particles, 3)) * 20. - 10.
initial_positions_b = np.random.random(size=(n_particles, 3)) * 20. - 10.

simulation.observe.number_of_particles(stride=1, types=["A"])

In [6]:
if os.path.exists(simulation.output_file):
os.remove(simulation.output_file)


Configured kernel context with:
--------------------------------
- kBT = 2.49434
- periodic b.c. = (true, true, true)
- box size = (20, 20, 20)
- particle types:
*  particle type "C" with D=1
*  particle type "A" with D=1
*  particle type "B" with D=1
- bimolecular reactions:
* Fusion A + B -> C with a rate of 1, an educt distance of 1, and weights 0.5 and 0.5



In [7]:
traj = readdy.Trajectory(simulation.output_file)


Analytical solution

In ReaDDy, one defines the intrinsic rate constant. In a well-mixed setting, one can use the following equation (see [1]) to obtain the corresponding macroscopic rate (ODE rate equation/ law of mass action)

$$k_\mathrm{on} = 4\pi (D_A + D_B) R \left[1 - \frac{\tanh(\kappa R)}{\kappa R}\right]$$

where

$$\kappa = \sqrt{\frac{\lambda_\mathrm{on}}{D_A + D_B}}$$

Parameters:

• intrinsic rate constant $\lambda_\mathrm{on} = 1\,\mathrm{ns}^{-1}$
• diffusion coefficient $D_A=D_B=1\, \mathrm{nm}^2 \mathrm{ns}^{-1}$
• reaction radius $R = 1\,\mathrm{nm}$

Law of mass action ODE for the concentration of A particles $[A](t) = a(t)$

$$\frac{\mathrm{d}a}{\mathrm{d}t} = -k_\mathrm{on}\,a^2 \quad\text{, with } a(0) = a_0$$

which yields

$$a(t) = \frac{1}{a_0^{-1} + k_\mathrm{on}t}$$

[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 [8]:
kappa = np.sqrt(lambda_on / 2.)
k_on = 4. * np.pi * 2. * 1. * (1. - np.tanh(kappa * 1.) / (kappa * 1.) )

def a(t):
return 1. / ((system.box_volume.magnitude / n_particles) + k_on * t)

t_range = np.linspace(0., 5000 * 1e-3, 10000)

In [9]:
plt.plot(time[::200]*1e-3, counts[::200] / system.box_volume.magnitude, "x", label="ReaDDy")
plt.plot(t_range, a(t_range), label=r"analytical $a(t)$")
plt.legend(loc="best")
plt.xlabel("time in nanoseconds")
plt.ylabel(r"concentration in nm$^{-3}$")
plt.show()


Note that in the diffusion-limited regime, the law of mass action solution generally does not reflect what can be observed from the Doi model.