# Michaelis-Menten - well mixed

In [2]:
%matplotlib inline
import os
import numpy as np

In [3]:
readdy.__version__

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

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.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.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)
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