## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width=7,
  fig.height=5
)

## ----setup--------------------------------------------------------------------
library(dynafluxr)
#pkgload::load_all("~/projs/dynafluxr/dev/dynafluxr")

## -----------------------------------------------------------------------------
fmeas=system.file("dataglyco/data.tsv", package="dynafluxr")
cat(head(readLines(fmeas)), sep="\n")

## -----------------------------------------------------------------------------
mf=read.delim(fmeas, comment.char="#")
print(head(mf))

## -----------------------------------------------------------------------------
fsto=system.file("dataglyco/network.txt", package="dynafluxr")
cat(readLines(fsto), sep="\n")

## -----------------------------------------------------------------------------
ddir=tempfile(pattern="glyco")
dir.create(ddir)
file.copy(c(fmeas, fsto), ddir)

## -----------------------------------------------------------------------------
fmeas=file.path(ddir, "data.tsv")
fsto=file.path(ddir, "network.txt")
mf=read.delim(fmeas, comment.char="#")
nm=colnames(mf)[-1L]
msp=fitsmbsp(mf$Time, mf[, nm], n=4, nki=5)
matplot(mf$Time, msp(mf$Time), t="l", ylab="Concentration [mM]", xlab="Time [min]", lwd=2)
coltr=apply(col2rgb(1:6)/255, 2, function(v) do.call(rgb, c(as.list(v), alpha=0.25)))
matpoints(mf$Time, mf[, nm], pch="o", cex=0.75, col=coltr)
legend("topright", legend=nm, lty=1:5, col=1:6, cex=0.75, lwd=2)

## -----------------------------------------------------------------------------
res=cli(c("-m", fmeas, "-s", fsto, "--decr", "GLC", "--incr", "FBP", "--skip", "24"))

## -----------------------------------------------------------------------------
list.files("glyco/data/")

## -----------------------------------------------------------------------------
res=cli(c("-m", fmeas, "-s", fsto, "--decr", "GLC", "--incr", "FBP", "--skip", "24", "-o", ""))

## -----------------------------------------------------------------------------
print(res$chi2tab)

## -----------------------------------------------------------------------------
nm=colnames(res$mf)[-1]
matplot(res$tpp, res$isp(res$tpp, nm), t="l", xlab="Time [min]", ylab="Concentration [mM]", lwd=2)
matpoints(res$tp, res$mf[,nm], pch="o", cex=0.5, col=coltr)
legend("topright", legend=nm, lty=1:5, col=1:6, cex=0.75, lwd=2)

## -----------------------------------------------------------------------------
res=cli(c("-m", fmeas, "-s", fsto, "--decr", "GLC", "--incr", "FBP", "--skip", "24", "--sf=FBP,F6P", "-o", "", "--sderr=GLC=0.05,G6P=0.05,FBP=0.1,F6P=0.05"))
print(res$chi2tab)

## -----------------------------------------------------------------------------
matplot(res$tpp, res$vsp(res$tpp), t="l", ylab="Rate [1/min]", xlab="Time [min]", lwd=2)
legend("topright", colnames(bsppar(res$vsp)$qw), lty=1:5, col=1:6, cex=0.75, lwd=2)

## -----------------------------------------------------------------------------
nm="G6P"
jnz=names(which(res$stofull[nm,] != 0)) # non-zero coeffs, i.e. reactions involved in G6P mass balance
fl=t(res$stofull[nm,jnz]*t(res$vsp(res$tpp, jnz))) # each rate v_j is multiplied by S_ij
fl=cbind(res$fsp(res$tpp, nm), fl) # add total flux
matplot(res$tpp, fl, t="l", ylab="Flux [mM/min]", xlab="Time [min]", main="G6P flux components", lwd=2)
abline(h=0)
legend("topright", legend=c("Total", jnz), lty=1:5, col=1:6, lwd=2)