Skip to contents

Uncertainty Quantification and global Sensitivity Analysis on AKAR4 (deterministic)

uqsaAKAR4.Rmd

Load Model and Data

First we load the model and data, for detailed explanations see previous examples (Simulate AKAR4 deterministically).

f  <- uqsa_example("AKAR4")
m <- model_from_tsv(f)
o <- as_ode(m)
#> Loading required namespace: pracma
C <- generateCode(o)
#> Warning in generateCode(o): This function will start a background yacas process via Ryacas.
#> There is currently no working way to reset/restart that process.
#>  It is therefore not advisable to generate the code for two different models in the same R session.
#> The definitions for the two models will be mixed up.
modelName <- write_and_compile(C)
#> './AKAR4_gvf.c' was created.
#> 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` -O2 -o './AKAR4.so' './AKAR4_gvf.c' `pkg-config --libs gsl`
ex <- experiments(m,o)

Create simulator

We construct a simulator function and and test it on default parameter values.

s <- simfi(ex,modelName,parMap=log10ParMap,omit=2) # or simulator.c
p <- log10(values(m$Parameter))
y <- s(p)

Uncertainty quantification: Sample the posterior distribution

For details on how to sample from the posterior distribution, see, e.g., the deterministic AKAR4 example. First we construct the prior. For this article, we override the information in the file and pick an iid normal prior with \(\sigma=2\) for all parameters:

sd <- rep(2,length(p))        # sigma = 2
dprior <- dNormalPrior(p,sd)
rprior <- rNormalPrior(p,sd)
mh <- mcmc(                   # metropolis hastings
    metropolis_update(        # random walk
        s,
        dprior=dprior,
        Sigma=diag(1/sd^2)
    )
)

We obtain a high quality sample by following the steps described in details in the deterministic AKAR4 example). This will take a minute or so on a laptop.

set.seed(137)
x <- mcmc_init(beta=1.0,p,s,ll,dprior)
h <- tune_step_size(mh,x)
#> acceptance rate: 0.66, step-size: 0.01;
#> acceptance rate: 0.46, step-size: 0.0174905;
#> acceptance rate: 0.28, step-size: 0.0270046;
#> acceptance rate: 0.33, step-size: 0.030052;
#> acceptance rate: 0.26, step-size: 0.0381874;
N <- 3e4

nCores <- parallel::detectCores()

bigSample <- parallel::mclapply(
    seq(nCores),
    function(i) {
        return(mh(x,N,h))
    },
    mc.cores=nCores
)

S_ <- Reduce(\(a,b) {rbind(a,b)},bigSample,init=NULL)
L_ <- Reduce(\(a,b) {c(a,attr(b,"logLikelihood"))},bigSample,init=NULL)
colnames(S_) <- names(p)

We should estimate the auto-correlation of this sample, to reduce the size effectively, for plotting and post-processing:

if (require(hadron)){
    res <- uwerr(data=L_)
    tau <- ceiling(res$tauint + res$dtauint)
} else {
    A <- acf(L_,lag.max=1000)
    i <- which(A$acf < exp(-3))[1]
    tau <- ceiling(sum(A$acf[seq(i)]))
}
#> Loading required package: hadron
#> 
#> Attaching package: 'hadron'
#> The following object is masked from 'package:base':
#> 
#>     kappa

print(tau)
#> [1] 7

if (1 < tau && tau < length(L_)){
    N <- length(L_)
    S_ <- S_[seq(1,N,by=tau),]
    L_ <- L_[seq(1,N,by=tau)]
    N <- length(L_)
}

Pairs plot (either with hexagonal binning or normal)

if (require(hexbin)){
    hexbin::hexplom(S_)    # binning
} else {
    pairs(
        S_,                # the sample
        pch=16,            # filled circle
        ces=4,             # so that points overlap
        col=rgb(0,0,1,0.1) # transparent blue
    )
}
#> Loading required package: hexbin

Sensitivity analysis on the posterior distribution

First we need to (re-)create the simulation result corresponding to the posterior parameter sample. This will take a minute or so on a laptop.

y <- s(t(S_))

Next we calculate the first order sensitivity indexes. Here for compute them for a specific experimental setting and time point:

E <- 2                      # Experiment idx
T <- 60                     # Time idx
fM <- y[[E]]$state[,T,]
SIappr <- gsa_binning(S_, t(fM))

and plot the result:

cols=rainbow(3)
par(mfrow = c(1, 1))
fM <- y[[E]]$state[,T,]
barplot(
    t(SIappr),
    col=cols,
    border="white",
    space=0.04,
    cex.axis=0.7,
    legend.text=names(p)
)

Here for compute the first order sensitivity indexes for all three experiments, looking at the different compounds in the system (on the x axis):

T <- 200                    #Time idx
par(mfrow = c(1,3))
for (E in seq(3,1,-1)){
fM <- y[[E]]$state[,T,]
SIappr <- gsa_binning(S_, t(fM))
lgd <- list(names(p), NULL,NULL)
barplot(
    t(SIappr),
    col = cols,
    cex.names = 0.7,
    las = 2,
    border = "white",
    space = 0.04,
    cex.axis = 1,
    main = names(ex)[E],
    legend.text = lgd[[E]])
}

Now we investigate the first order sensitivity indexes at different time points (time is on the x axis) for each of the compounds and one experiment.

E <- 3
timePts <- seq(2,40,by=4)  #Avoid first time point
nStates=dim(y[[E]]$state)[1]
cols=rainbow(3)
par(mfrow = c(1, nStates))
lgd <- rep(list(NULL),nStates)
lgd[[nStates-1]] <- names(p)
for (Cm in seq(nStates)) {
    fM <- y[[E]]$state[Cm,timePts,]
    SIappr <-gsa_binning(S_, t(fM))
    names(SIappr) <- timePts
    barplot(
        t(SIappr),
        col=cols,
        names.arg = timePts,
        border="white",
        cex.axis=1,
        main=dimnames(y[[E]]$state)[[1]][Cm],
        xlab="Time",
        ylab="SI",
        legend.text=lgd[[Cm]],
        args.legend = list(x = "topright", inset=c(-0.1, 0), cex=0.7)
    )
}
mtext(names(ex[E]), side=3, line = - 1.5, outer=TRUE)