Skip to content

Commit 44db10d

Browse files
dkaschekdkaschek
committed
EXAMPLE: Add bile acid PD model to vignette folder
1 parent 6b87993 commit 44db10d

17 files changed

+703
-440
lines changed

‎data/BAdata.tab‎

Lines changed: 229 additions & 229 deletions
Large diffs are not rendered by default.

‎inst/examples/errmodel.R‎

Lines changed: 14 additions & 82 deletions
Original file line numberDiff line numberDiff line change
@@ -3,12 +3,14 @@
33

44
library(dMod)
55
library(ggplot2)
6+
library(dplyr)
7+
68
setwd(tempdir())
79

810
# Set up reactions
9-
f <- NULL
10-
f <- addReaction(f, "A", "B", "k1*A", "Production of B")
11-
f <- addReaction(f, "B", "C", "k2*B", "Production of C")
11+
f <- eqnvec() %>%
12+
addReaction("A", "B", "k1*A", "Production of B") %>%
13+
addReaction("B", "C", "k2*B", "Production of C")
1214

1315
# Define observables and error model
1416
observables <- eqnvec(B_obs = "B + off_B")
@@ -24,13 +26,16 @@ e <- Y(errors, g, attach.input = FALSE,
2426

2527
# Generate parameter transformation
2628
innerpars <- getParameters(model, g, e)
27-
trafo <- repar("x~x", x = innerpars)
28-
trafo <- repar("x~1", x = "A", trafo)
29-
trafo <- repar("x~0", x = c("B", "C"), trafo)
30-
trafo <- repar("x~exp(x)", x = innerpars, trafo)
29+
covariates <- data.frame(Aini = 1:2, row.names = c("C1", "C2"))
3130

32-
p <- P(trafo, condition = "C1", modelname = "parfn", compile = FALSE) +
33-
P(trafo, condition = "C2", modelname = "parfn", compile = FALSE)
31+
p <-
32+
eqnvec() %>%
33+
define("x~x", x = innerpars) %>%
34+
define("x~0", x = c("B", "C")) %>%
35+
branch(table = covariates) %>%
36+
insert("A~Aini", Aini = Aini) %>%
37+
insert("x~exp(x)", x = innerpars) %>%
38+
P(modelname = "parfn", compile = FALSE)
3439

3540
compile(g, x, e, p, output = "errtest_total")
3641
#compile(g, x, e, p, cores = 4)
@@ -53,79 +58,6 @@ profiles <- profile(obj + constraintL2(myfit$argument, 10),
5358
plotProfile(profiles)
5459

5560

56-
## Test if fixing parameters works -----------------------------------
57-
58-
# Set up reactions
59-
f <- NULL
60-
f <- addReaction(f, "A", "B", "k1*A", "Production of B")
61-
f <- addReaction(f, "B", "C", "k2*B", "Production of C")
62-
63-
# Define observables and error model
64-
observables <- eqnvec(B_obs = "B + off_B")
65-
errors <- eqnvec(B_obs = "sqrt((sigma_rel*B_obs)^2 + sigma_abs^2)")
66-
67-
# Generate dMod objects
68-
model <- odemodel(f, modelname = "errtest", compile = FALSE, solver = "deSolve")
69-
x <- Xs(model, optionsSens = list(method = "lsoda"), optionsOde = list(method = "lsodes"))
70-
g <- Y(observables, x, parameters = c("k1", "k2"),
71-
compile = FALSE, modelname = "obsfn")
72-
e <- Y(errors, g, attach.input = FALSE,
73-
compile = FALSE, modelname = "errfn")
74-
75-
# Generate parameter transformation
76-
innerpars <- getParameters(model, g, e)
77-
trafo <- repar("x~x", x = innerpars)
78-
trafo <- repar("x~1", x = "A", trafo)
79-
trafo <- repar("x~0", x = c("B", "C"), trafo)
80-
trafo <- repar("x~exp(x)", x = innerpars, trafo)
81-
82-
p <- P(trafo, condition = "C1", modelname = "parfn", compile = FALSE) +
83-
P(trafo, condition = "C2", modelname = "parfn", compile = FALSE)
84-
85-
compile(g, x, e, p, output = "errtest_total")
86-
#compile(g, x, e, p, cores = 4)
87-
88-
89-
## Simulate data
90-
ptrue <- c(k1 = -2, k2 = -3, sigma_rel = log(.1), sigma_abs = log(.1))
91-
fixed <- c(off_B = -3)
92-
times <- seq(0, 50, 1)
93-
prediction <- (g*x*p)(times, ptrue, fixed = fixed, deriv = TRUE)
94-
datasheet <- subset(as.data.frame(prediction, errfn = e), name == "B_obs")
95-
datasheet$value <- datasheet$value + rnorm(length(datasheet$value), sd = datasheet$sigma)
96-
data <- as.datalist(datasheet)
97-
98-
## Fit data with error model
99-
obj <- normL2(data, g*x*p, e)
100-
myfit <- trust(obj, ptrue, rinit = 1, rmax = 10)
101-
fits <- mstrust(obj, center = ptrue, sd = 3, fits = 10, cores = 10)
102-
profiles <- profile(obj + constraintL2(myfit$argument, 10),
103-
myfit$argument, names(myfit$argument), limits = c(-5, 5), cores = length(myfit$argument))
104-
plotProfile(profiles)
105-
106-
107-
108-
## Fit externally
109-
out <- runbg({
110-
trust(obj, ptrue, rinit = 1, rmax = 10)
111-
}, machine = "localhost", filename = "test", input = c("obj", "ptrue"), compile = TRUE)
112-
113-
## Fit on grid
114-
out <- runbg_bwfor({
115-
trust(obj, ptrue, rinit = 1, rmax = 10)
116-
}, machine = "bwfor", filename = "test", input = c("obj", "ptrue"), compile = TRUE, nodes = 2, cores = 1, walltime = "00:01:00")
117-
118-
119-
120-
## Plotting
121-
out <- as.data.frame((g*x*p)(times = seq(0, 50, len = 100), pars = myfit$argument), errfn = e)
122-
ggplot(out, aes(x = time, y = value, ymin = value-sigma, ymax = value+sigma, color = condition, fill = condition)) +
123-
facet_wrap(~name, scales = "free") +
124-
geom_line() + geom_ribbon(alpha = .2, lty = 0) +
125-
geom_point(data = as.data.frame(data)) +
126-
theme_dMod() + scale_color_dMod() + scale_fill_dMod()
127-
128-
12961
## Compute prediction profile
13062
datapoint <- datapointL2(name = "A", time = 10, value = "d1", sigma = .05, condition = "C1")
13163
par <- trust(normL2(data, g*x*p, e) + datapoint, c(ptrue, d1 = 0), rinit = 1, rmax = 10)$argument

‎vignettes/BA_PDmodel.so‎

-12.2 KB
Binary file not shown.

‎vignettes/BA_PDmodel_s.so‎

-37.2 KB
Binary file not shown.

‎vignettes/BA_basemodel.so‎

-12.1 KB
Binary file not shown.

‎vignettes/BA_basemodel_s.so‎

-16.4 KB
Binary file not shown.

‎vignettes/BA_publication_part1.R‎

Lines changed: 89 additions & 55 deletions
Original file line numberDiff line numberDiff line change
@@ -1,122 +1,156 @@
11

2+
3+
# Parameter Estimation in an ODE Model of Bile Acid Transport
4+
5+
## Preliminaries
6+
27
library(dMod)
38
library(dplyr)
49
library(ggplot2)
10+
library(pander)
511
theme_set(theme_dMod())
6-
rm(list = ls())
712

8-
## Read the data ---------------------------------------------
13+
setwd("/tmp")
14+
15+
16+
## The model scheme
17+
18+
## Read the data
19+
920

1021
data(BAdata)
1122
data <- BAdata %>%
1223
filter(experiment == "exp1") %>%
1324
as.datalist()
1425
covtable <- covariates(data)
15-
print(covtable)
16-
1726
plot(data)
1827

28+
## Setting up the base model
29+
30+
### Setting up the prediction function
1931

20-
## Set up a base model ----------------------------------------
2132

22-
# Set up base reactions
2333
reactions <- eqnlist() %>%
34+
#addReaction(from, to, rate, description) %>%
2435
addReaction("Tca_buffer", "Tca_cyto", "import_Tca*Tca_buffer", "Basolateral uptake") %>%
2536
addReaction("Tca_cyto", "Tca_buffer", "export_Tca_baso*Tca_cyto", "Basolateral efflux") %>%
2637
addReaction("Tca_cyto", "Tca_canalicular", "export_Tca_cana*Tca_cyto", "Canalicular efflux") %>%
2738
addReaction("Tca_canalicular", "Tca_buffer", "transport_Tca*Tca_canalicular", "Transport bile") %>%
28-
addReaction("0", "cations", "0", "Cation concentration in buffer")
39+
#addReaction("0", "cations", "0", "Cation concentration in buffer") %>%
40+
{.}
2941

30-
# Translate to ODEs and let cation concentration change transport rate
3142
odes <- reactions %>%
3243
as.eqnvec() %>%
33-
reparameterize("x ~ (x*crit/(crit + cations))", x = "transport_Tca")
34-
35-
# Set up events to control the experiment
36-
events <- data.frame(
37-
var = c("Tca_buffer", "cations"),
38-
time = c("change_buffer", "change_cations"),
39-
value = c("value_buffer", "value_cations"),
40-
method = c("replace", "replace")
41-
)
42-
43-
# Make ODE model available as prediction function
44-
noSens <- c("change_buffer", "change_cations", "value_cations", "cations")
44+
insert("transport_Tca ~ (transport_Tca*crit/(crit + cations))")
45+
46+
events <- eventlist() %>%
47+
#addEvent(var, time, value, method = "replace") %>%
48+
addEvent("Tca_buffer", "t_addTca", "amount_Tca") %>%
49+
addEvent("Tca_buffer", "t_removeTca", "0") %>%
50+
addEvent("cations", "t_removeCa", "0")
51+
52+
noSens <- c("t_removeCa", "t_addTca", "t_removeTca", "cations")
4553
x <- odemodel(odes, events = events, fixed = noSens, modelname = "BA_basemodel") %>%
4654
Xs()
4755

48-
# Make observation function available
49-
g <- eqnvec(
50-
buffer = "s*Tca_buffer",
51-
cellular = "s*(Tca_cyto + Tca_canalicular)") %>%
56+
## Demonstration of the usage of x()
57+
n <- getParameters(x)
58+
pars <- runif(length(n), min = 0, max = 100)
59+
names(pars) <- n
60+
61+
times <- 0:200
62+
pred <- x(times, pars)
63+
plot(pred) + theme(legend.position = "none")
64+
65+
### Setting up the observation function
66+
67+
g <- eqnvec(buffer = "s*Tca_buffer",
68+
cellular = "s*(Tca_cyto + Tca_canalicular)") %>%
5269
Y(f = x, modelname = "obsfn", compile = TRUE)
5370

54-
# Parameterize the model
71+
# Demonstration of the usage of g
72+
pars["s"] <- 1
73+
pred <- (g*x)(times, pars)
74+
plot(pred) + theme(legend.position = "none")
75+
76+
77+
### Parameterization
78+
79+
80+
5581
parameters <- getParameters(g, x)
56-
transformation <- eqnvec() %>%
57-
reparameterize("x~x", x = parameters) %>%
58-
reparameterize("x~0", x = c("Tca_cyto", "Tca_canalicular", "value_buffer", "value_cations")) %>%
59-
reparameterize("x~1", x = "cations")
60-
82+
p <- eqnvec() %>%
83+
define("x ~ x", x = parameters) %>%
84+
define("x ~ 0", x = c("Tca_cyto", "Tca_canalicular", "Tca_buffer")) %>%
85+
define("cations ~ 1") %>%
86+
define("t_addTca ~ 0") %>%
87+
define("t_removeTca ~ 30") %>%
88+
branch(covtable) %>%
89+
define("t_removeCa ~ time", time = ifelse(changeCa == "yes", 30, 1e3)) %>%
90+
insert("x ~ exp(X)", x = parameters, X = toupper(parameters)) %>%
91+
P()
92+
6193

62-
p <- P()
63-
for (c in rownames(covtable)) {
64-
p <-p + transformation %>%
65-
reparameterize("x~cations", x = "change_cations", cations = covtable[c, "cations"]) %>%
66-
reparameterize("x~buffer", x = "change_buffer", buffer = covtable[c, "drug_time"]) %>%
67-
reparameterize("x~exp(x)", x = parameters) %>%
68-
P(condition = c)
69-
}
94+
## Parameter estimation
7095

96+
set.seed(12837)
7197
estimate <- getParameters(p)
72-
pars <- structure(rep(-1, length(estimate)), names = estimate)
98+
pars <- structure(runif(length(estimate), min = -2, max = 0), names = estimate)
7399
times <- seq(0, 200, 1)
74100
(g*x*p)(times, pars) %>% plot(data = data, time <= 180)
75101

76-
obj <- normL2(data, g*x*p) + constraintL2(pars, sigma = 10)
77-
myfit <- trust(obj, pars, rinit = 1, rmax = 10, iterlim = 500)
78102

79-
(g*x*p)(times, myfit$argument) %>% plot(data = data)
103+
obj <- normL2(data, g*x*p) + constraintL2(0*pars, sigma = 10)
104+
myfit <- trust(obj, pars, rinit = 1, rmax = 10, iterlim = 500)
80105

106+
(g*x*p)(times, myfit$argument) %>% plot(data = data, time <= 180)
107+
### Identifiability
81108
profiles <- profile(obj, myfit$argument, names(myfit$argument), cores = 4)
82109
profiles %>% plotProfile(mode == "data")
83110

84111

85112
# Add additional data
86-
data$exp1_30_no_0_0_30 <- rbind(data$exp1_30_no_0_0_30, data.frame(
87-
name = "Tca_buffer", time = 0, value = 21.5, sigma = 2
88-
))
113+
data$exp1_yes_no_0 <- data$exp1_yes_no_0 %>%
114+
rbind(data.frame(name = "Tca_buffer", time = 1, value = 21.5, sigma = 2))
89115

90-
obj <- normL2(data, g*x*p) + constraintL2(pars, sigma = 10)
116+
obj <- normL2(data, g*x*p) + constraintL2(0*pars, sigma = 10)
91117
myfit <- trust(obj, pars, rinit = 1, rmax = 10, iterlim = 500)
92118

93-
(g*x*p)(times, myfit$argument) %>% plot(data = data)
119+
pred <- (g*x*p)(times, myfit$argument)
120+
pred %>% plot(data = data, name %in% c("buffer", "cellular", "Tca_buffer"))
94121

95-
profiles <- profile(obj, myfit$argument, c("s", "Tca_buffer"), cores = 2)
122+
123+
profiles <- profile(obj, myfit$argument, c("S", "AMOUNT_TCA"), cores = 2)
96124
profiles %>% plotProfile(mode == "data")
97125

126+
98127
# Explore parameter space
99128
fits <- mstrust(obj, pars, fits = 30, rinit = 1, rmax = 10, samplefun = "runif", min = -5, max = 5, cores = 4)
100129
parframe <- as.parframe(fits)
101-
plotValues(parframe[1:29,])
102-
plotPars(parframe[1:29,])
103-
subframe <- unique(parframe[1:29,])
130+
plotValues(parframe[1:28,])
131+
plotPars(parframe[1:28,])
132+
subframe <- unique(parframe[1:28,])
104133

105134
prediction <- predict(g*x*p, times = times, pars = subframe, data = data)
106135

107-
ggplot(prediction, aes(x = time, y = value, color = as.factor(.value))) + facet_wrap(~name*cations, scales = "free") +
136+
ggplot(prediction, aes(x = time, y = value, color = as.factor(round(.value)))) + facet_wrap( ~ name*changeCa, scales = "free") +
108137
geom_line() +
109-
geom_point(data = as.data.frame(data), color = "black")
138+
geom_point(data = as.data.frame(data), color = "black") +
139+
scale_color_discrete(name = "obj. value")
140+
110141

111142
profiles <- lapply(1:nrow(subframe), function(i) {
112-
profile(obj, as.parvec(subframe, i), names(myfit$argument), cores = 4, algoControl = list(gamma = 10))
143+
profile(obj, as.parvec(subframe, i), names(myfit$argument), limits = c(-4, 4), cores = 4)
113144
})
114145
names(profiles) <- c("obj = 23", "obj = 24", "obj = 27")
115146
profiles %>% plotProfile(mode == "data")
116147

148+
117149
partable <- confint(profiles[[1]], level = 0.9, val.column = "value")
118150

151+
119152
# Save objects from this session
120153
save.image(file = "part1.RData")
121154

122-
155+
156+

‎vignettes/BA_publication_part2.R‎

Lines changed: 10 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -24,34 +24,24 @@ plot(data)
2424

2525
## Add new condition to parameterization --------------------------
2626

27-
28-
# Parameterize the model
2927
parameters <- getParameters(g, x)
30-
transformation <- eqnvec() %>%
31-
reparameterize("x~x", x = parameters) %>%
32-
reparameterize("x~0", x = c("Tca_cyto", "Tca_canalicular", "Tca_buffer", "value_cations")) %>%
33-
reparameterize("x~1", x = "cations")
34-
35-
36-
for (c in rownames(covtable)[3]) {
37-
p <-p + transformation %>%
38-
reparameterize("x~cations", x = "change_cations", cations = covtable[c, "cations"]) %>%
39-
reparameterize("x~buffer", x = "change_buffer", buffer = covtable[c, "tca_time"]) %>%
40-
reparameterize("x~exp(x)", x = parameters) %>%
41-
P(condition = c)
42-
}
43-
44-
estimate <- getParameters(p)
28+
p_add <- getEquations(p, conditions = "exp1_no_no_0") %>%
29+
define("t_addTca ~ 60") %>%
30+
define("x ~ 1e3", x = c("t_removeTca", "t_removeCa")) %>%
31+
P(condition = "exp3_no_no_0")
32+
33+
estimate <- getParameters(p, p_add)
4534
pars <- structure(rep(-1, length(estimate)), names = estimate)
4635
pars[partable$name] <- partable$value
4736

4837
times <- seq(0, 200, 1)
49-
(g*x*p)(times, pars) %>% plot(data = data, time <= 180)
38+
(g*x*(p + p_add))(times, pars) %>% plot(data = data, time <= 180)
5039

51-
obj <- normL2(data, g*x*p) + constraintL2(pars, sigma = 10)
40+
y <- g * x * (p + p_add)
41+
obj <- normL2(data, y) + constraintL2(pars, sigma = 10)
5242
myfit <- trust(obj, pars, rinit = 1, rmax = 10, iterlim = 500)
5343

54-
(g*x*p)(times, myfit$argument) %>% plot(data = data)
44+
y(times, myfit$argument) %>% plot(data = data)
5545

5646
profiles_part2 <- profile(obj, myfit$argument, names(myfit$argument), cores = 4)
5747
profiles_part2 %>% plotProfile(mode == "data")

0 commit comments

Comments
 (0)