Figure 4.31:

Stochastic simulation of the first-order reactions A-> B-> C starting with 4000 A~molecules.

Code for Figure 4.31

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
# Converted from stoch_big.m
import numpy as np
from misc import stairs, save_ascii

k = np.array([2.0, 1.0])
stoi = np.array([[-1, 1, 0], [0, -1, 1]])
stoiT = stoi.T
nsim = 8000
time = np.zeros(nsim+1)
x = np.zeros((3, nsim+1))
x[0, 0] = 4000

rng = np.random.RandomState(2)

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

ts, xs = stairs(time, x.T)
table = np.column_stack([ts[:, 0], xs])
save_ascii('stoch_big.dat', table)