Example Simulation
simAKAP79.RmdThis article provides code to simulate the AKAP79 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).
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] «/home/andrei/R/library/uqsa/extdata/AKAP79/AKAP79_Compound.tsv» belongs to Document «AKAP79»
#> I'll take this as the Model Name.
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/andrei/R/library/uqsa/extdata/AKAP79/AKAP79_gvf.c' `pkg-config --libs gsl`
comment(modelName)
#> [1] "./AKAP79.so"
# model related functions, in R, e.g. AKAP79_default() parameters
source(uqsa_example("AKAP79",pat='^AKAP79[.]R$',full.names=TRUE))
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)
This also includes instructions for the simulator.
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 prameters for simulation
nInput <- length(experiments[[1]]$input)
p <- head(AKAP79_default(),-nInput)
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
This will make a function s, which will always simulate the scenarios described in the experiments list, but for user supplied parameters.
s <- simulator.c(experiments,modelName,noise=TRUE)
#> Loading required package: rgsl
y <- s(p)Plot
E <- 10 # which experiment to plot
out <- experiments[[E]]$outputValues$AKAR4pOUT
err <- experiments[[E]]$errorValues$AKAR4pOUT
tm <- experiments[[E]]$outputTime
par(bty='n',xaxp=c(80,120,4))
plot(tm,
y[[E]]$func[1,,1],
type='l',
ylim=c(90,130), ylab="AKAR4p",
xlab="t",
main=sprintf("Experiment %i",E),
lwd=1.5,
col="purple"
)
points(tm,out)
arrows(x0=tm,x1=tm,y0=out,y1=out+err,angle=90,length=0.025)
arrows(x0=tm,x1=tm,y0=out,y1=out-err,angle=90,length=0.025)
gg-Plot
require(ggplot2)
#> Loading required package: ggplot2
D<-data.frame(time=experiments[[E]]$outputTime,
AKAR4p=experiments[[E]]$outputValues$AKAR4pOUT,
AKAR4pERR=experiments[[E]]$errorValues$AKAR4pOUT,
sim=y[[E]]$func[1,,1])
ggplot(D) +
geom_linerange(mapping=aes(x=time,y=AKAR4p,ymin=AKAR4p-AKAR4pERR,ymax=AKAR4p+AKAR4pERR),na.rm=TRUE) +
geom_point(mapping=aes(x=time,y=AKAR4p),na.rm=TRUE) +
geom_line(mapping=aes(x=time,y=sim),color="purple",lwd=1.2)