Simulation of a Stochastic Model for the AKAR4 reaction network

require(SBtabVFGEN)
#> Loading required package: SBtabVFGEN
library(uqsa)
library(GillespieSSA2)

This article provides code to simulate the AKAR4 stochastic model (one time, no sampling). We are plotting the model with default parameters which are not expected to fit the data (this is the starting point).

The Stochastic Model

When the copy number of molecular species in a reaction network system (AKAP79 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. The likelihood of this model is hard to compute; however, we can easily sample trajectories from this stochastic model using the Gillespie 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).

Load the Model

This model is included with the package. To load your own model, see the user model article.

modelFiles <- uqsa_example("AKAR4",full.names=TRUE)
SBtab <- SBtabVFGEN::sbtab_from_tsv(modelFiles)
#> [tsv] file[1] «AKAR4_100nM.tsv» belongs to Document «AKAR4»
#>  I'll take this as the Model Name.
#> AKAR4_100nM.tsv  AKAR4_25nM.tsv  AKAR4_400nM.tsv  AKAR4_Compound.tsv  AKAR4_Experiments.tsv  AKAR4_Output.tsv  AKAR4_Parameter.tsv  AKAR4_Reaction.tsv
modelName <- checkModel("AKAR4",uqsa_example("AKAR4",pat="_gvf[.]c$"))
#> building a shared library from c source, and using GSL odeiv2 as backend (pkg-config is used here).
#> cc -shared -fPIC `pkg-config --cflags gsl` -o './AKAR4.so' '/home/andreikr/.local/R/library/uqsa/extdata/AKAR4/AKAR4_gvf.c' `pkg-config --libs gsl`
comment(modelName)
#> [1] "./AKAR4.so"

# model related functions, in R, e.g. AKAP79_default() parameters
source(uqsa_example("AKAR4",pat='^AKAR4[.]R$'))
#> Loading required package: deSolve
print(AKAR4_default())
#> kf_C_AKAR4 kb_C_AKAR4 kcat_AKARp 
#>      0.018      0.106     10.200

Load Experiments (data)

experiments <- sbtab.data(SBtab)

# for example, these is the initial state of experiment 1:
print(experiments[[1]]$initialState)
#>   AKAR4 AKAR4_C  AKAR4p       C 
#>     0.2     0.0     0.0     0.4

# pick parameters for simulation
p <- SBtab$Parameter[["!DefaultValue"]]
par_names <- SBtab$Parameter[["!Name"]]
names(p) <- par_names
print(p)
#> kf_C_AKAR4 kb_C_AKAR4 kcat_AKARp 
#>      0.018      0.106     10.200

Simulate

Function simulator.stoch will output a function s, which will always simulate the scenarios described in experiment e (i.e., same initial conditions, same inputs), but for user supplied parameters.

exp_ind <- 1 #index of the experiment to consider
e <- experiments[[exp_ind]]

# generate a function (s) that simulates a trajectory given a parameter in input
s <- simulator.stoch(experiment = e, model.tab = SBtab, vol = 4e-16, nStochSim = 3)
# vol indicates the volume in m^3 where the reactions take place
# nStochSim indicates how many times we simulate the stochastic model; the output trajectory is the average of all (in this case nStochSim = 3) the simulated trajectories

# simulate a trajectory (y) given parameter p
y <- s(p)

Plot

par(bty='n',xaxp=c(80,120,4))
plot(e$outputTimes, e$outputValues$AKAR4pOUT,ylim=c(90,200), ylab="AKAP79",
     xlab="t",
     main=sprintf("Experiment %i",exp_ind),
     lwd=1.5)
lines(e$outputTimes, y$output, col="blue")

gg-Plot

require(ggplot2)
#> Loading required package: ggplot2

D<-data.frame(time=experiments[[exp_ind]]$outputTime,
              AKAP79=experiments[[exp_ind]]$outputValues$AKAR4pOUT,
              AKAP79ERR=experiments[[exp_ind]]$errorValues$AKAR4pOUT,
              sim=y$output)
ggplot(D) +
  geom_linerange(mapping=aes(x=time,y=AKAP79,ymin=AKAP79-AKAP79ERR,ymax=AKAP79+AKAP79ERR),na.rm=TRUE) +
  geom_point(mapping=aes(x=time,y=AKAP79),na.rm=TRUE) +
  geom_line(mapping=aes(x=time,y=sim),color="purple",lwd=1.2)