Simulation of a Stochastic Model for the AKAP79 reaction network
simAKAP79stochastic.Rmd
require(SBtabVFGEN)
#> Loading required package: SBtabVFGEN
library(uqsa)
library(GillespieSSA2)This article provides code to simulate the AKAP79 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("AKAP79",full.names=TRUE)
SBtab <- SBtabVFGEN::sbtab_from_tsv(modelFiles)
#> [tsv] file[1] «AKAP79_Compound.tsv» belongs to Document «AKAP79»
#> I'll take this as the Model Name.
#> AKAP79_Compound.tsv AKAP79_Experiments.tsv AKAP79_Expression.tsv AKAP79_Input.tsv AKAP79_Output.tsv AKAP79_Parameter.tsv AKAP79_Reaction.tsv X0uM_cAMPCaN_AKAP79_0_nM_cAMP.tsv X0uM_cAMPCaN_only_0_nM_cAMP.tsv X0uM_cAMPno_CaN_0_nM_cAMP.tsv X1000nM_cAMPCaN_AKAP79_1_uM_cAMP.tsv X1000nM_cAMPCaN_only_1_uM_cAMP.tsv X1000nM_cAMPno_CaN_1_uM_cAMP.tsv X100nM_cAMPCaN_AKAP79_100_nM_cAMP.tsv X100nM_cAMPCaN_only_100_nM_cAMP.tsv X100nM_cAMPno_CaN_100_nM_cAMP.tsv X2000nM_cAMPCaN_AKAP79_2_uM_cAMP.tsv X2000nM_cAMPCaN_only_2_uM_cAMP.tsv X2000nM_cAMPno_CaN_2_uM_cAMP.tsv X200nM_cAMPCaN_AKAP79_200_nM_cAMP.tsv X200nM_cAMPCaN_only_200_nM_cAMP.tsv X200nM_cAMPno_CaN_200_nM_cAMP.tsv X500nM_cAMPCaN_AKAP79_500_nM_cAMP.tsv X500nM_cAMPCaN_only_500_nM_cAMP.tsv X500nM_cAMPno_CaN_500_nM_cAMP.tsv
modelName <- checkModel("AKAP79",uqsa_example("AKAP79",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 './AKAP79.so' '/home/andreikr/.local/R/library/uqsa/extdata/AKAP79/AKAP79_gvf.c' `pkg-config --libs gsl`
comment(modelName)
#> [1] "./AKAP79.so"
# In the terminal run the following:
# ``ode.sh -R --maxima AKAP79.tar.gz > AKAP79.R```
# (Read article "..." in the documentation)
system2(command = "ode_uqsa", args = c("-R", "--maxima", "AKAP79.tar.gz"), stdout = "AKAP79.R")
#> Warning in system2(command = "ode_uqsa", args = c("-R", "--maxima",
#> "AKAP79.tar.gz"), : error in running command
# model related functions, in R, e.g. AKAP79_default() parameters
source(uqsa_example("AKAP79",pat='^AKAP79[.]R$'))
print(AKAP79_default())
#> kf_Rii_C__RiiP_C kf_RiiP_CxcAMP__RiiP_C_cAMP
#> 33.00000 0.49600
#> kf_RiiP_cAMPxC__RiiP_C_cAMP kb_RiiP_cAMPxC__RiiP_C_cAMP
#> 0.00545 0.01560
#> kb_RiiPXcAMP__RiiP_cAMP kf_RiiPXcAMP__RiiP_cAMP
#> 0.00160 0.01500
#> kf_RiiPxC__RiiP_C kb_RiiPxC__RiiP_C
#> 0.03800 0.00260
#> kf_cAMPxRii__Rii_cAMP kb_cAMPxRii__Rii_cAMP
#> 0.01500 0.00160
#> kf_Rii_CxcAMP__Rii_C_cAMP kb_Rii_CxcAMP__Rii_C_cAMP
#> 0.49600 1.41300
#> kf_RiixC__Rii_C kf_Rii_cAMPxC__Rii_C_cAMP
#> 2.10000 0.29840
#> kb_Rii_cAMPxC__Rii_C_cAMP kf_Rii_C_cAMP__RiiP_C_cAMP
#> 0.01800 33.00000
#> kb_RiixC__Rii_C AKAPoff_1
#> 0.00030 2.60000
#> AKAPoff_3 AKAPon_1
#> 20.00000 0.45000
#> AKAPon_3 kf_C_AKAR4
#> 2.00000 0.01800
#> kb_C_AKAR4 kcat_AKARp
#> 0.10600 10.20000
#> kmOFF kmON
#> 100.00000 1.00000
#> KD_T b_AKAP
#> 0.70000 0.00000Load Experiments (data)
experiments <- sbtab.data(SBtab)
# for example, these are the input and initial state of experiment 1:
print(experiments[[1]]$input)
#> [1] 1
print(experiments[[1]]$initialState)
#> Rii cAMP RiiP Rii_C RiiP_cAMP
#> 6.30 0.00 0.00 0.63 0.00
#> RiiP_C RiiP_C_cAMP C Rii_cAMP Rii_C_cAMP
#> 0.00 0.00 0.00 0.00 0.00
#> CaN RiiP_CaN RiiP_cAMP_CaN AKAR4 AKAR4_C
#> 1.50 0.00 0.00 0.20 0.00
#> AKAR4p
#> 0.00
# pick parameters for simulation
nInput <- length(experiments[[1]]$input)
p <- SBtab$Parameter[["!DefaultValue"]]
par_names <- SBtab$Parameter[["!Name"]]
names(p) <- par_names
print(p)
#> kf_Rii_C__RiiP_C kf_RiiP_CxcAMP__RiiP_C_cAMP
#> 33.00000 0.49600
#> kf_RiiP_cAMPxC__RiiP_C_cAMP kb_RiiP_cAMPxC__RiiP_C_cAMP
#> 0.00545 0.01560
#> kb_RiiPXcAMP__RiiP_cAMP kf_RiiPXcAMP__RiiP_cAMP
#> 0.00160 0.01500
#> kf_RiiPxC__RiiP_C kb_RiiPxC__RiiP_C
#> 0.03800 0.00260
#> kf_cAMPxRii__Rii_cAMP kb_cAMPxRii__Rii_cAMP
#> 0.01500 0.00160
#> kf_Rii_CxcAMP__Rii_C_cAMP kb_Rii_CxcAMP__Rii_C_cAMP
#> 0.49600 1.41300
#> kf_RiixC__Rii_C kf_Rii_cAMPxC__Rii_C_cAMP
#> 2.10000 0.29840
#> kb_Rii_cAMPxC__Rii_C_cAMP kf_Rii_C_cAMP__RiiP_C_cAMP
#> 0.01800 33.00000
#> kb_RiixC__Rii_C AKAPoff_1
#> 0.00030 2.60000
#> AKAPoff_3 AKAPon_1
#> 20.00000 0.45000
#> AKAPon_3 kf_C_AKAR4
#> 2.00000 0.01800
#> kb_C_AKAR4 kcat_AKARp
#> 0.10600 10.20000
#> kmOFF kmON
#> 100.00000 1.00000
#> KD_T
#> 0.70000Simulate
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 <- 9 #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)
# vol indicates the volume in m^3 where the reactions take place
# 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,180), 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)