Figure 5.12:

Solution to master equation for A + B <-> C starting with 200 A molecules, 1000 B molecules and 0 C molecules, k_1=1/200, k_{-1}=3.

Code for Figure 5.12

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
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import expm, null_space

# Stochastic reaction simulation
# A +B <-> C
# Determine equilibrium via chemical master equation

# Set state vector x: [A; B; C]
Vol = 200
x = np.array([1, 5, 0]) * Vol
nmol = x[0]
# Set the reaction rate vector k
k = np.array([1 / Vol, 3])
extent = np.linspace(1, 0, nmol+1)

# Construct A of the chemical master equation
# dP/dt = A*P
A = np.zeros((nmol+1, nmol+1))
for i in range(nmol+1):
    n = i
    nu = nmol - n
    a = n
    b = x[1] - nu
    c = nu
    A[i, i] = -(k[0] * a * b + k[1] * c)
    if i > 0:
        A[i, i - 1] = k[1] * (c + 1)
    if i < nmol:
        A[i, i + 1] = k[0] * (a + 1) * (b + 1)

v = null_space(A)
v = v / np.sum(v)
plt.figure()
plt.plot(v)

# Set initial condition for p
p = np.zeros((nmol+1, 1))
p[-1] = 1 * nmol

deltat = 0.01
tfinal = 1
time = np.arange(0, tfinal + deltat, deltat)
Ad = expm(A * deltat)
p_all = np.zeros((nmol+1, len(time)))
p_all[:, 0] = p.flatten()
for i in range(1, len(time)):
    p_all[:, i] = Ad @ p_all[:, i-1]
    p_all[:, i] = p_all[:, i] / np.sum(p_all[:, i]) * nmol

# Decide on density of extents for plotting
elines = 50
eskip = int(np.ceil(nmol / elines))
eindex = np.arange(0, nmol+1, eskip)
tlines = 30
tskip = int(np.ceil((len(time) - 1) / tlines))
tstart = 2
tindex = np.arange(tstart, len(time), tskip) - 2
plt.figure()

X, Y = np.meshgrid(time[tindex], extent[eindex])
plt.pcolormesh(X, Y, p_all[eindex, :][:, tindex])
plt.colorbar()
plt.show(block=False)

# Save output for plotting in readable .dat format
with open('prob_dis_large.dat', 'w') as f:
    for j in range(len(eindex)):
        for i in range(len(tindex)):
            output = [time[tindex[i]], extent[eindex[j]], p_all[eindex[j],tindex[i]]]
            np.savetxt(f, [output], fmt='%g', newline='\n')
        if j < len(eindex) - 1:
            f.write('\n')