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] «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"
# 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.00000
#> AKAR4_ConservedConst CaN_ConservedConst
#> 0.20000 1.50000
#> Rii_C_ConservedConst cAMP_ConservedConst
#> 0.63000 0.00000
#> Rii_ConservedConst
#> 6.30000Conservation laws (determined by sbtab_to_vfgen and
saved as RDS file):
clf <- uqsa_example("AKAP79",pat='Laws[.]RDS$')
cl <- readRDS(clf)With conservation laws, some species are calculated algebraically. Their initial values are turned into input parameters (using the found law):
With this hypothetical relationship (\(c\) is a constant):
\[ A+B = c\,, \]
we can determine that
\[ c = A_0 + B_0\,. \]
And thus we can replac either of the two species:
\[ A(t) = A_0 + B_0 - B(t) \]
And \(A_0+B_0\) are turned into an
input called A_ConservedConst (the \(c\) above) with the value determined from
the stated initial condition.
Load Experiments (data)
This also includes instructions for the simulator.
experiments <- sbtab.data(SBtab,cl) # with conservation laws
# for example, these are the input
# and initial state of experiment 1:
print(experiments[[1]]$input)
#> b_AKAP AKAR4_ConservedConst CaN_ConservedConst
#> 1.00 0.20 1.50
#> Rii_C_ConservedConst cAMP_ConservedConst Rii_ConservedConst
#> 0.63 0.00 6.30
print(experiments[[1]]$initialState)
#> RiiP RiiP_cAMP RiiP_C RiiP_C_cAMP C
#> 0 0 0 0 0
#> Rii_cAMP Rii_C_cAMP RiiP_CaN RiiP_cAMP_CaN AKAR4_C
#> 0 0 0 0 0
#> AKAR4p
#> 0
# 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)