Simulation of a Stochastic Model for the AKAR4 reaction network

library(uqsa)
library(errors)

This article provides code to simulate the AKAR4 stochastic model (one time, no sampling).

The Stochastic Model

When the copy number of molecular species in a reaction network system (AKAR4 in our case) is low, we cannot model the amount of molecules deterministically (e.g., with an ODE model), because the stochasticity in the reactions that take place cannot be ignored. In particular, the time at which reactions take place is random, as well as the specific reactions that take place (i.e., what pair of molecules react). To model such system we can use the master equation: we model the (integer) number of each molecule species in the system (e.g., proteins) and how the number of each molecule type (randomly) evolves in time. We can easily sample trajectories from this stochastic model using the Gillespie’s Stochastic Simulation algorithm. Given the current amount of each molecular species in the system at a given time point, we can sample the time at which the next reaction takes place and we can sample the type of reaction (i.e., what pair of molecules react).

To obtain the stochastic model for the reaction network, we need to determine the reaction propensities. These can be derived from the reaction rate coefficients (which are the parameters that we usually use in our models and on which we perform uncertainty quantification). To derive the reaction propensities we referred to this article.

Load the Model

The AKAR4 model is included with the package. To load your own model, see the article “Build and simulate your own Model” article.

f <- uqsa_example("AKAR4")
m <- model_from_tsv(f)
o <- as_ode(m,cla=FALSE)
sm <- makeGillespieModel(m)
C <- generateGillespieCode(sm,LV=1e8)

where LV is Avogadro’s number multiplied by the system’s volume. Higher numbers mean more particles (more reactions per time-span).

Load Experiments (data)

ex <- experiments(m,o) # same as with ODE models

Simulate

Function simstoch will output a function s, which will always simulate the scenarios described in the experiments ex (i.e., same initial conditions, same inputs), but for user supplied parameters, as with simulator.c for ODEs

We auto-generate C code, which uses the kinetic laws, but rescales all parameters internally to have the unit 1/s and all compound species are rescaled to be in particle counts. So, the internal simulator can be used with the list of experiments that we use for ODE models, with no changes. This rescaling is performed as described in this article.

We can also write code and compile that model to a shared library manually:

print(length(C))
#> [1] 142
cat(C,sep="\n",file="./AKAR4_st.c")

In the shell:

echo $0
cc -shared -fPIC -O2 -o akar4-st.so ./AKAR4_st.c -lm
#> sh

And next we simulate this compiled model:

p <- values(m$Parameter)
st <- simstoch(ex,model.so="./akar4-st.so")
y <- st(p)

Plot the results

tm <- ex[[1]]$outputTimes

par(bty="n")
plot(
    as.errors(tm),
    ex[[1]]$data,
    xlab="time",
    ylab="AKAR4p",
    main=names(ex)[1],
    ylim=c(90,210)
)
lines(
    tm,
    y[[1]]$func,
    lwd=2,
    col="red",
    type='s'
)

Benchmark

Stochastic simulations are slower than ODE simulations, this is a benchmark of the builtin stochastic solver:

if (require(rbenchmark)){
    BM <- benchmark(
        uqsa = {y <- st(p)},
        replications = 5000
    )
    print(BM)
}
#> Loading required package: rbenchmark
#>   test replications elapsed relative user.self sys.self user.child sys.child
#> 1 uqsa         5000   0.424        1     0.418    0.006          0         0