added BA example update

simonbeyer1 · simonbeyer1 · commit 8e5416cd75b7 · 2025-01-31T15:26:56.000+01:00
diff --git a/inst/examples/BA_transport.R b/inst/examples/BA_transport.R
@@ -0,0 +1,345 @@
+## Load Libraries -----
+rm(list=ls())
+library(dMod)
+library(dplyr)
+library(ggplot2)
+theme_set(theme_dMod())
+setwd(tempdir()) # change if you want to keep the files that dMod creates internally
+set.seed(5555)
+
+## Set up the ODE model -----------------------------------------------------------------------
+
+# Define via reactions
+
+reactions <- eqnlist() %>%
+  addReaction("TCA_buffer", "TCA_cell",  rate = "import*TCA_buffer", description = "Uptake") %>%
+  addReaction("TCA_cell", "TCA_buffer",  rate = "export_sinus*TCA_cell", description = "Sinusoidal export") %>%
+  addReaction("TCA_cell", "TCA_cana",    rate = "export_cana*TCA_cell", description = "Canalicular export") %>%
+  addReaction("TCA_cana", "TCA_buffer",  rate = "reflux*TCA_cana", description = "Reflux into the buffer")
+
+# Define via ODE system (equivalent)
+
+# f <- eqnvec(TCA_buffer = "-import * TCA_buffer + export_sinus * TCA_cell + reflux * TCA_cana",
+#             TCA_cana = "export_cana * TCA_cell - reflux * TCA_cana",
+#             TCA_cell = "import * TCA_buffer - export_sinus * TCA_cell - export_cana * TCA_cell")
+
+# Translate reactions into ODE model object
+mymodel <- odemodel(reactions, modelname = "bamodel", compile = T)
+
+
+# Generate trajectories for the default condition
+x <- Xs(mymodel, condition = NULL)
+
+
+# For demonstration define parameters (initials and dynamic parameters)
+pars <- c(reflux = 0.1, TCA_buffer = 1, TCA_cell = 0, TCA_cana = 0, import = 0.2, export_sinus = 0.2, export_cana = 0.04)
+
+
+# Plot trjectories
+times <- seq(0, 50, 0.01) 
+out <- x(times, pars)
+plot(out) # Note this is a ggplot object
+
+# Define observables buffer and cellular
+observables <- eqnvec(buffer = "s*TCA_buffer", cellular = "s*(TCA_cana + TCA_cell)")
+g <- Y(observables, f = x, condition = NULL, compile = TRUE, modelname = "obsfn")
+
+
+
+## Simulate data -------------------------------------------------------------------------------------------------------------
+
+# Fix parameters to steady state values (calculated on paper / by hand), replace buffer -> TCA_buffer = 0
+pars["TCA_cell"] <- 0.3846154
+pars["TCA_cana"] <- 0.1538462
+pars["TCA_buffer"] <- 0
+pars["s"] <- 1e3
+
+# Simulate trajectories 
+# Evaluate the expression separately
+out <- (g * x)(times, pars, conditions = "closed")
+plot(out)
+
+
+
+
+
+# simulate data for time points
+timesD <- c(0.1, 1, 3, 7, 11, 15, 20, 41)
+datasheet <- subset(as.data.frame(out), time %in% timesD & name %in% names(observables))
+datasheet <- within(datasheet, {
+  sigma <- 0.1*value
+  value <- rnorm(length(value), value, sigma)
+})
+data <- as.datalist(datasheet)
+plot(out, data)
+
+
+# Define parameter transformations using define(), insert() and branch(). Old function repar also avaiable!
+innerpars <- getParameters(x,g)
+trafo <- NULL %>% 
+  define("x~x", x = innerpars) %>% # identity
+  define("TCA_buffer~0") %>% 
+  insert("x~exp(log(10)*x)", x = .currentSymbols)
+
+
+# # Explicit trafo
+# trafo <- eqnvec(TCA_buffer = "0", 
+#                 TCA_cell = "exp(log(10)*TCA_cell)",
+#                 TCA_cana = "exp(log(10)*TCA_cana)", 
+#                 import = "exp(log(10)*import)",
+#                 export_sinus = "exp(log(10)*export_sinus)", 
+#                 export_cana = "exp(log(10)*export_cana)",
+#                 reflux = "exp(log(10)*reflux)", 
+#                 s = "exp(log(10)*s)")
+
+
+p <- P(trafo, condition = "closed")
+# Set every parameter to -1 in the log-space
+outerpars <- getParameters(p)
+pouter <- structure(rep(-1, length(outerpars)), names = outerpars)
+plot((g*x*p)(times, pouter),data)
+
+
+## Use simulate data to calibrate outer model parameters -------------------------------------------------------------
+# # One Fit
+obj <- normL2(data, g * x * p) + constraintL2(pouter, sigma = 20)
+
+# Fit on time (starting from pouter)
+myfit_1 <- trust(obj, pouter, rinit = 0.1, rmax = 5)
+mypred_1 <- (g * x * p)(times, myfit_1$argument)
+plot(mypred_1, data)
+
+## Handling different experimental conditions
+
+# Define reflux and open condition according to "Dynamic Modelling in R p. 19-20"
+pars["reflux"] <- 1e3
+out <- (g * x)(times, pars, conditions = "open")
+datasheet <- subset(as.data.frame(out), time %in% timesD & name %in% names(observables))
+datasheet <- within(datasheet, {
+  sigma <- 0.1*value
+  value <- rnorm(length(value), value, sigma)
+})
+data <- data + as.datalist(datasheet)
+plotData(data)
+
+# Parameter Trafo, usage of "+" operator for trafo functions (output of P())
+trafo <- getEquations(p, conditions = "closed")
+trafo["reflux"] <- "exp(log(10)*reflux_open)"
+p <- p + P(trafo, condition = "open")
+
+outerpars <- getParameters(p)
+pouter <- structure(rep(-1, length(outerpars)), names = outerpars)
+plot((g*x*p)(times, pouter),data)
+
+
+# Objective function
+obj <- normL2(data, g * x * p) + constraintL2(pouter, sigma = 20)
+
+# Evaluation of obj at pouter
+obj(pouter)
+
+
+# Fit 50 times, sample with sd=4 around pouter
+out_frame <- mstrust(obj, pouter, sd = 4, studyname = "bamodel", cores=8, fits=50, iterlim = 200)
+out_frame <- as.parframe(out_frame)
+plotValues(out_frame) # Show "Waterfall" plot
+plotPars(out_frame) # Show parameter plot
+bestfit <- as.parvec(out_frame)
+
+
+# Plot predictions along data
+plot((g * x * p)(times, bestfit), data)
+
+# Plot sensis
+plot(getDerivs((g * x * p)(times, bestfit)))
+
+# Calculate Parameter Profiles and plot different contributions (for identifiablility only "data" is of interest)
+profiles <- profile(obj, bestfit, names(bestfit), method = "optimize", cores = 12)
+plotProfile(profiles)
+plotProfile(profiles,mode %in% c("data", "prior"))
+plotPaths(profiles, whichPar = "TCA_cana")
+plotPaths(profiles, whichPar = "export_cana")
+plotPaths(profiles, whichPar = "s")
+
+# Tighten model assumptions with steady state constraint
+pSS <- NULL
+equations <- getEquations(p) # getEquations retrieves the "trafo" from p 
+conditions <- names(equations)
+
+# Set up steady state constraint
+for (n in conditions) {
+  equations[[n]]["TCA_cana"] <- "exp(log(10)*export_cana)*exp(log(10)*TCA_cell)/exp(log(10)*reflux)"
+  pSS <- pSS + P(equations[[n]], condition = n)
+}
+outerpars <- getParameters(pSS)
+pouter <- structure(rep(-1, length(outerpars)), names = outerpars)
+plot((g*x*pSS)(times, pouter),data)
+
+# Objective function
+obj <- normL2(data, g * x * pSS, times = times, attr.name = "data") + constraintL2(pouter, sigma = 5, attr.name = "prior")
+
+# Fit 50 times again
+out_frame <- mstrust(obj,pouter,studyname = "bamodel", cores=10,fits=50)
+out_frame <- as.parframe(out_frame)
+plotValues(out_frame) # Show "Waterfall" plot
+plotPars(out_frame) # Show parameter plot
+bestfit <- as.parvec(out_frame)
+plot((g * x * pSS)(times, bestfit), data)
+
+
+# Calculate Parameter Profiles
+profiles <- profile(obj, bestfit, names(bestfit), cores = 10, limits = c(lower = -10, upper = 10),
+                    stepControl = list(stepsize = 1e-4, min = 1e-4, max = Inf, atol = 1e-2, rtol = 1e-2, limit = 100, stop = "data"))
+plotProfile(profiles,mode %in% c("data", "prior"))
+plotPaths(profiles, whichPar = "s")
+
+
+# s and initial value TCA_cell render model structurally non identifiable
+#  --> fix s to 1 (on log scale to 0)
+
+trafo <- getEquations(pSS)
+trafo <- repar("x~0", x = "s", trafo)
+p <- P(trafo)
+
+
+# Objective function
+outerpars <- getParameters(p)
+pouter <- structure(rep(-1, length(outerpars)), names = outerpars)
+obj <- normL2(data, g * x * p, times = times, attr.name = "data") + constraintL2(pouter, sigma = 20, attr.name = "prior")
+
+# Fit 50 times again
+out_frame <- mstrust(obj,pouter,studyname = "bamodel", cores=10,fits=50)
+out_frame <- as.parframe(out_frame)
+plotValues(out_frame) # Show "Waterfall" plot
+plotPars(out_frame) # Show parameter plot
+bestfit <- as.parvec(out_frame)
+plot((g * x * p)(times, bestfit), data)
+
+# remove the prior
+obj <- normL2(data, g * x * p, times = times, attr.name = "data")
+
+# refit
+refit <- trust(obj, bestfit, rinit = 1, rmax = 10)
+plot((g * x * p)(times, refit$argument), data)
+
+profiles <- profile(obj, refit$argument, names(refit$argument), cores = 10, limits = c(lower = -5, upper = 5), method = "optimize",
+                    stepControl = list(stop = "data"))
+plotProfile(profiles, mode == "data")
+plotPaths(profiles, whichPar = "reflux_open")
+
+# Note: The parameter "reflux_open" is still practical non-identifiable 
+# The profile is open to the left -> A possible reduction of the model would be the assumption of an immediate reflux, i.e. the limit case reflux_open -> infinity
+# An pragmatic but ugly way to circumvent the reformulation of the ODE system is to fix the reflux_open parameter to a high value. E.g. 1e3
+
+# One could also check the models ability to produce reliable predictions, by the calculation of prediction uncertainty with profile likelihood
+# The calculation of prediction confidence intervals is done at next
+
+## Prediction uncertainty taken from validation profile --------------------------------------------------------------------------
+
+# choose sigma below 1 percent of the prediction in order to pull the prediction strongly towards d1
+obj.validation <- normL2(data, g * x * p, times = times, attr.name = "data") +
+  datapointL2(name = "TCA_cell", time = 10, value = "v", sigma = 1, attr.name = "validation", condition = "closed")
+
+# If sigma is not known, and you therefore decide to calculate prediction confidence intervals, just choose a very small sigma, in order to "pull strongly" on the trajectory
+
+# refit
+myfit <- trust(obj.validation, parinit = c(v = 180, bestfit), rinit = 1, rmax = 10, iterlim = 1000)
+
+# Calculate profile
+validation_profile <- profile(obj.validation, myfit$argument, "v", cores = 4, method = "optimize",
+                              optControl = list(rinit = .1, rmax = 10, iterlim = 100))
+
+
+# plotProfile(validation_profile) # This also plos the prediction colums, which is a bug in the code.
+plotProfile(validation_profile, mode %in% c("validation", "data")) # Plots only the two contributions validation and data, along with the sum (total)
+# Is the contribution of validation small?? # If yes: The total aligns with a prediction profile
+
+# Confidence Interval of the prediction
+confint(validation_profile, val.column = "value")
+
+
+
+## Prediction band (prediction uncertainty for several time points) --------------------------------------------------------------
+# Here we calculate a prediction CI for different timepoints. In the end we interpolate to a "prediction band"
+prediction_band <- do.call(rbind, lapply(seq(10, 50, 10), function(t) {
+  
+  cat("Computing prediction profile for t =", t, "\n")
+  
+  obj.validation <- normL2(data, g * x * p, times = times, attr.name = "data") + 
+    datapointL2(name = "TCA_cell", time = t, value = "v", sigma = 1, attr.name = "validation", condition = "closed")
+  
+  refit <- trust(obj.validation, parinit = c(v = 180, bestfit), rinit = 1, rmax = 10, iterlim = 1000)
+  
+  profile_prediction <- profile(obj.validation, myfit$argument, "v", cores = 4, method = "optimize")
+  
+  d1 <- confint(profile_prediction, val.column = "value")
+  
+  # Output
+  data.frame(time = t, condition = "closed", name = "TCA_cell",  d1[-1])
+  
+}))
+
+prediction <- (g * x * p)(times, refit$argument) %>%
+  as.data.frame() 
+
+
+prediction_band_spline <- data.frame(
+  time = prediction$time[prediction$time>=10],
+  value = prediction$value[prediction$time>=10],
+  condition = "closed",
+  name = "TCA_cell",
+  lower = spline(prediction_band$time, prediction_band$lower, xout = prediction$time[prediction$time>=10])$y,
+  upper = spline(prediction_band$time, prediction_band$upper, xout = prediction$time[prediction$time>=10])$y
+)
+
+# Create the ggplot
+ggplot(prediction, aes(x = time, y = value, color = condition)) +
+  geom_line() +  # Line connecting the points for each condition
+  geom_ribbon(data = prediction_band_spline, aes(x = time, ymin = lower, ymax = upper, fill = condition), 
+              lty = 0, alpha = .3, show.legend = F) +  # Show ribbon in the legend
+  geom_point(data = prediction_band, aes(x = time, y = lower, color = condition), shape = 4, show.legend = F) + 
+  geom_point(data = prediction_band, aes(x = time, y = upper, color = condition), shape = 4, show.legend = F) + 
+  facet_wrap(~ name, scales = "free_y") +  # Facet by 'name' column
+  labs(
+    x = "Time",
+    y = "Value",
+    color = "Condition"
+  ) +
+  dMod::theme_dMod() +  # Apply dMod theme
+  dMod::scale_color_dMod() +  # Apply dMod color scale to lines
+  dMod::scale_fill_dMod()   # Apply the same color scale to the fill
+
+
+# ## Alternative implementation of the Steady state via implicit parameter transformation of the steady state  -----------------------------------------------
+# 
+# # Redefine reactions in order to control the standard and open condition by events
+# reactions <- eqnlist() %>% 
+#   addReaction("TCA_buffer", "TCA_cell", rate = "import*TCA_buffer", description = "Uptake")%>% 
+#   addReaction("TCA_cell", "TCA_buffer", rate = "export_sinus*TCA_cell", description = "Sinusoidal export")%>% 
+#   addReaction("TCA_cell", "TCA_cana", rate = "export_cana*TCA_cell", description = "Canalicular export")%>% 
+#   addReaction("TCA_cana", "TCA_buffer", rate = "(reflux*(1-switch) + reflux_open*switch)*TCA_cana", description = "Reflux into the buffer")%>% 
+#   addReaction("0", "switch", rate = "0", description = "Create a switch")
+# 
+# events <- NULL
+# events <- addEvent(events, var = "TCA_buffer", time = 0, value = 0)
+# events <- addEvent(events, var = "switch" , time = 0, value = "OnOff")
+# mymodel <- odemodel(reactions, modelname = "bamodel2", events = events)
+# x <- Xs(mymodel)
+# 
+# 
+# 
+# # Replace one reaction with a analytical expression for the conserved quantity: TCR_tot
+# f <- as.eqnvec(reactions)[c("TCA_buffer", "TCA_cana", "TCA_cell")]
+# f["TCA_cell"] <- "TCA_buffer + TCA_cana + TCA_cell - TCA_tot"
+# pSS <- P(f, method = "implicit", compile = TRUE, modelname = "pfn") 
+# 
+# observables <- eqnvec(buffer = "s*TCA_buffer", cellular = "s*(TCA_cana + TCA_cell)")
+# 
+# innerpars <- unique(c(getParameters(mymodel), getSymbols(observables), getSymbols(f)))
+# trafo <- repar("x~x" , x = innerpars)
+# trafo <- repar("x~0" , x = reactions$states, trafo)
+# 
+# trafo <- repar("x~exp(log(10)*x)", x = setdiff(innerpars, "OnOff"), trafo)
+# p <- P(repar("OnOff~0", trafo), condition = "closed") + P(repar("OnOff~1", trafo), condition = "open")
+# g <- Y(observables, f = x, compile = TRUE, modelname = "obsfn2")
+# 
