Figure 4.34

Figure 4.34:

Envelope versus time for hepatitis B virus model; deterministic and average stochastic models.

Code for Figure 4.34

hbv_common.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
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118 # Converted from hbv_common.m
import numpy as np
from scipy.integrate import solve_ivp

def hbv_deriv(t, x, p):
    k = p['k']
    return [
        k[1]*x[1] - k[3]*x[0],
        k[0]*x[0] - k[1]*x[1] - k[5]*x[1]*x[2],
        k[2]*x[0] - k[4]*x[2] - k[5]*x[1]*x[2]
    ]

def _load_stochsim_lib():
    """Try to load the compiled C stochsim shared library via ctypes."""
    import ctypes, os
    so = os.path.join(os.path.dirname(__file__), 'stochsim_py.so')
    if not os.path.exists(so):
        return None
    try:
        lib = ctypes.CDLL(so)
        lib.stochsim.restype = None
        lib.stochsim.argtypes = [
            ctypes.POINTER(ctypes.c_double), ctypes.c_int,  # x0, nx
            ctypes.POINTER(ctypes.c_double),                 # k
            ctypes.POINTER(ctypes.c_double), ctypes.c_int,  # stoiT, nrxns
            ctypes.POINTER(ctypes.c_double), ctypes.c_int,  # tout, nts
            ctypes.POINTER(ctypes.c_double),                 # xout
            ctypes.c_ulong,                                  # seed
        ]
        return lib
    except Exception:
        return None

_stochsim_lib = _load_stochsim_lib()

def hbv_stochsim(x0, k, stoiT, tout, rng):
    """Gillespie algorithm for HBV model. Uses compiled C extension if available."""
    import ctypes
    x0 = np.array(x0, dtype=float)
    nts = len(tout)
    nx = len(x0)

    if _stochsim_lib is not None:
        # Fast path: call compiled C library
        xout = np.zeros((nx, nts), dtype=np.float64, order='C')
        stoiT_cm = np.ascontiguousarray(stoiT, dtype=np.float64)  # nx x nrxns
        tout_c   = np.ascontiguousarray(tout,  dtype=np.float64)
        k_c      = np.ascontiguousarray(k,     dtype=np.float64)
        x0_c     = np.ascontiguousarray(x0,    dtype=np.float64)
        seed = rng.randint(0, 2**31)
        nrxns = stoiT.shape[1]
        _stochsim_lib.stochsim(
            x0_c.ctypes.data_as(ctypes.POINTER(ctypes.c_double)), nx,
            k_c.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
            stoiT_cm.ctypes.data_as(ctypes.POINTER(ctypes.c_double)), nrxns,
            tout_c.ctypes.data_as(ctypes.POINTER(ctypes.c_double)), nts,
            xout.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
            ctypes.c_ulong(seed),
        )
        return xout

    # Fallback: pure Python Gillespie
    x = x0.copy()
    t = tout[0]
    n_rxns = stoiT.shape[1]
    tout_idx = 0
    results = np.zeros((nx, nts))
    results[:, 0] = x0

    while tout_idx < nts - 1:
        r = np.zeros(n_rxns)
        r[0] = k[0]*x[0]; r[1] = k[1]*x[1]; r[2] = k[2]*x[0]
        r[3] = k[3]*x[0]; r[4] = k[4]*x[2]; r[5] = k[5]*x[1]*x[2]
        rtot = np.sum(r)
        if rtot == 0:
            for idx in range(tout_idx+1, nts):
                results[:, idx] = x
            break
        p2 = rng.rand(2)
        timenew = t - np.log(p2[0]) / rtot
        while tout_idx < nts - 1 and tout[tout_idx+1] <= timenew:
            tout_idx += 1
            results[:, tout_idx] = x
        if tout_idx >= nts - 1:
            break
        cumr = np.cumsum(r)
        m = np.searchsorted(cumr, p2[1]*rtot)
        if m < n_rxns:
            x = np.maximum(x + stoiT[:, m], 0)
        t = timenew
    return results

def hbv_common():
    k = np.array([1., 0.025, 1000., 0.25, 2., 7.5e-6])
    p = {'k': k}
    stoi = np.array([
        [0, 1, 0], [1, -1, 0],
        [0, 0, 1], [-1, 0, 0],
        [0, 0, -1], [0, -1, -1]
    ], dtype=float)
    x0 = np.array([1., 0., 0.])
    stoiT = stoi.T
    tfin = 200.0
    nts  = 200
    tout = np.linspace(0, tfin, nts)

    p['x0'] = x0; p['stoiT'] = stoiT; p['tout'] = tout; p['tfin'] = tfin; p['nts'] = nts

    sqeps = np.sqrt(np.finfo(float).eps)
    sol = solve_ivp(lambda t, x: hbv_deriv(t, x, p), [0, tfin], x0,
                    method='Radau', t_eval=tout, rtol=sqeps, atol=sqeps)
    tsolver = sol.t
    x = sol.y.T
    return p, tsolver, x

if __name__ == '__main__':
    p, tsolver, x = hbv_common()
    print(f'hbv_common ran ok, x shape = {x.shape}')

hbv_binary.npz

hbv_binary.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 # Converted from hbv_binary.m
# [depends] hbv_common.py stochsim_py.c
# [makes] hbv_binary.npz
# [rule] copymat
# Runs 500 stochastic simulations and saves results as hbv_binary.npz.
import numpy as np
from hbv_common import hbv_common, hbv_stochsim

p, tsolver, xdet = hbv_common()

rng = np.random.RandomState(0)

nsim = 10  # change to 500 for figure in text
nts   = p['nts']
x0    = p['x0']
nspec = len(x0)
xavg  = np.zeros((nspec, nts))

for i in range(nsim):
    xsim = hbv_stochsim(x0, p['k'], p['stoiT'], p['tout'], rng)
    xavg += xsim / nsim

xnull = hbv_stochsim(x0, p['k'], p['stoiT'], p['tout'], np.random.RandomState(121))

np.savez('hbv_binary.npz',
         xavg=xavg,
         xout=xsim,
         tout=p['tout'],
         x_det=xdet.T,
         xnull=xnull)
print('hbv_binary.npz saved')

stochsim_py.c (Python C extension source)


  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
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103 /*
 * stochsim_py.c -- Gillespie stochastic simulator for HBV model.
 * Called from Python via ctypes.
 *
 * Interface:
 *   void stochsim(const double *x0, int nx,
 *                 const double *k,
 *                 const double *stoiT,  // shape (nx, nrxns), column-major
 *                 int nrxns,
 *                 const double *tout, int nts,
 *                 double *xout,         // shape (nx, nts), column-major (output)
 *                 unsigned long seed)
 */

#include <math.h>
#include <stdlib.h>
#include <string.h>

/* Simple LCG random number generator (thread-safe, seedable per call) */
static unsigned long rng_state;

static void rng_seed(unsigned long seed) {
    rng_state = seed;
}

static double rng_rand(void) {
    rng_state = rng_state * 6364136223846793005ULL + 1442695040888963407ULL;
    return (double)((rng_state >> 33) & 0x7FFFFFFF) / (double)0x7FFFFFFF;
}

void stochsim(const double *x0, int nx,
              const double *k,
              const double *stoiT,   /* nx * nrxns, column-major */
              int nrxns,
              const double *tout, int nts,
              double *xout,          /* nx * nts, column-major (output) */
              unsigned long seed)
{
    rng_seed(seed);

    double *x = (double *)malloc(nx * sizeof(double));
    memcpy(x, x0, nx * sizeof(double));

    double time = tout[0];
    int iout = 0;

    /* HBV-specific rates (6 reactions, 3 species) */
    double r[6];

    while (iout < nts) {
        double x0v = x[0], x1v = x[1], x2v = x[2];

        r[0] = k[0] * x0v;
        r[1] = k[1] * x1v;
        r[2] = k[2] * x0v;
        r[3] = k[3] * x0v;
        r[4] = k[4] * x2v;
        r[5] = k[5] * x1v * x2v;

        double rtot = r[0]+r[1]+r[2]+r[3]+r[4]+r[5];

        if (rtot == 0.0) {
            while (iout < nts) {
                for (int s = 0; s < nx; s++)
                    xout[s * nts + iout] = x[s];
                iout++;
            }
            break;
        }

        double p1 = rng_rand();
        double p2 = rng_rand();
        double timenew = time - log(p1) / rtot;

        /* record output times that fall before timenew */
        while (iout < nts && tout[iout] <= timenew) {
            for (int s = 0; s < nx; s++)
                xout[s * nts + iout] = x[s];
            iout++;
        }
        if (iout >= nts)
            break;

        /* choose reaction */
        double rcum = 0.0;
        double target = p2 * rtot;
        int m = 0;
        for (m = 0; m < nrxns - 1; m++) {
            rcum += r[m];
            if (rcum > target)
                break;
        }

        /* apply stoichiometry */
        for (int s = 0; s < nx; s++) {
            double xnew = x[s] + stoiT[s * nrxns + m];
            x[s] = xnew < 0.0 ? 0.0 : xnew;
        }
        time = timenew;
    }

    free(x);
}

main.py


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15 # Converted from hbv_stoch_env.m
# [depends] hbv_common.py hbv_binary.npz hbv_binary.py stochsim_py.c
# This py-file loads data file 'hbv_binary.npz'.
import numpy as np
from misc import save_ascii

data  = np.load('hbv_binary.npz')
tout  = data['tout']
x_det = data['x_det'].T
xavg  = data['xavg']
xout  = data['xout']
xnull = data['xnull']

table = np.column_stack([tout, x_det, xavg.T, xout.T, xnull.T])
save_ascii('hbv_stoch_env.dat', table)