Figure 4.37:

Deterministic simulation of reaction A + B <-> C compared to stochastic simulation.

Code for Figure 4.37

Text of the GNU GPL.

main.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
# Converted from stochnon.m
# A + B <-> C stochastic simulation
import numpy as np
from scipy.integrate import solve_ivp
from misc import stairs, octave_save

nmolec = 400
k1 = 1.0; k2 = 1.0
k = np.array([k1/nmolec, k2])
stoi = np.array([[-1, -1, 1], [1, 1, -1]])
stoiT = stoi.T
nspec = 3
x_arr = np.zeros((3, nmolec*4+1))
x_arr[0, 0] = nmolec
x_arr[1, 0] = 0.9*nmolec
x_arr[2, 0] = 0
nsim = nmolec*4
time = np.zeros(nsim+1)

rng = np.random.RandomState(0)
for n in range(nsim):
    r = np.array([k[0]*x_arr[0, n]*x_arr[1, n],
                  k[1]*x_arr[2, n]])
    rtot = np.sum(r)
    if rtot == 0:
        time[n+1] = time[n]; x_arr[:, n+1] = x_arr[:, n]; continue
    pp = rng.rand(2)
    tau = -np.log(pp[0]) / rtot
    time[n+1] = time[n] + tau
    m = int(np.sum(np.cumsum(r) <= pp[1]*rtot))
    if m >= 2: m = 1
    x_arr[:, n+1] = x_arr[:, n] + stoiT[:, m]

ts, xs = stairs(time, x_arr.T)
ts = ts[:, 0]
x0_0 = x_arr[0, 0]
xscale = xs / x0_0

ca0 = x_arr[0, 0] / x0_0
cb0 = x_arr[1, 0] / x0_0
cc0 = x_arr[2, 0] / x0_0

def f_det(t, x):
    ca = x[0]
    cb = ca - ca0 + cb0
    cc = cc0 + ca0 - ca
    return [-k1*ca*cb + k2*cc]

sqeps = np.sqrt(np.finfo(float).eps)
sol = solve_ivp(f_det, [time[0], time[-1]], [ca0], method='Radau',
                t_eval=time, rtol=sqeps, atol=sqeps)
ca = sol.y[0]
cb = ca - ca0 + cb0
cc = cc0 + ca0 - ca

table1 = np.column_stack([ts, xscale])
table2 = np.column_stack([time, ca, cb, cc])
octave_save('stochnon.dat', ('table1', table1), ('table2', table2))