Skip to content

Commit 690c297

Browse files
dkaschekdkaschek
committed
VIGNETTE: Preparations of QSP workshop.
1 parent f091dad commit 690c297

16 files changed

+663
-0
lines changed

‎data/BAdata.tab‎

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

‎man/unique.parframe.Rd‎

Lines changed: 23 additions & 0 deletions
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

‎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.R‎

Lines changed: 115 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,115 @@
1+
2+
library(dMod)
3+
library(dplyr)
4+
library(ggplot2)
5+
theme_set(theme_dMod())
6+
7+
## Read the data ---------------------------------------------
8+
9+
data(BAdata)
10+
data <- BAdata %>%
11+
filter(experiment == "exp1") %>%
12+
as.datalist(split.by = c("experiment", "cations", "compound", "dose", "tca_time", "drug_time"))
13+
covtable <- covariates(data)
14+
print(covtable)
15+
16+
plot(data)
17+
18+
19+
## Set up a base model ----------------------------------------
20+
21+
# Set up base reactions
22+
reactions <- eqnlist() %>%
23+
addReaction("Tca_buffer", "Tca_cyto", "import_Tca*Tca_buffer", "Basolateral uptake") %>%
24+
addReaction("Tca_cyto", "Tca_buffer", "export_Tca_baso*Tca_cyto", "Basolateral efflux") %>%
25+
addReaction("Tca_cyto", "Tca_canalicular", "export_Tca_cana*Tca_cyto", "Canalicular efflux") %>%
26+
addReaction("Tca_canalicular", "Tca_buffer", "transport_Tca*Tca_canalicular", "Transport bile") %>%
27+
addReaction("0", "cations", "0", "Cation concentration in buffer")
28+
29+
# Translate to ODEs and let cation concentration change transport rate
30+
odes <- reactions %>%
31+
as.eqnvec() %>%
32+
reparameterize("x ~ (x*crit/(crit + cations))", x = "transport_Tca")
33+
34+
# Set up events to control the experiment
35+
events <- data.frame(
36+
var = c("Tca_buffer", "cations"),
37+
time = c(30, "time_remove_cations"),
38+
value = c(0, 0),
39+
method = c("replace", "replace")
40+
)
41+
42+
# Make ODE model available as prediction function
43+
x <- odemodel(odes, events = events, fixed = c("time_remove_cations", "cations"), modelname = "BA_basemodel") %>%
44+
Xs()
45+
46+
# Make observation function available
47+
g <- eqnvec(
48+
buffer = "s*Tca_buffer",
49+
cellular = "s*(Tca_cyto + Tca_canalicular)") %>%
50+
Y(f = x, modelname = "obsfn", compile = TRUE)
51+
52+
# Parameterize the model
53+
parameters <- getParameters(g, x)
54+
transformation <- eqnvec() %>%
55+
reparameterize("x~x", x = parameters) %>%
56+
reparameterize("x~0", x = c("Tca_cyto", "Tca_canalicular")) %>%
57+
reparameterize("x~1", x = "cations")
58+
59+
60+
p <- P()
61+
for (c in rownames(covtable)) {
62+
p <-p + transformation %>%
63+
reparameterize("x~cations", x = "time_remove_cations", cations = covtable[c, "cations"]) %>%
64+
reparameterize("x~exp(x)", x = parameters) %>%
65+
P(condition = c)
66+
}
67+
68+
estimate <- getParameters(p)
69+
pars <- structure(rep(-1, length(estimate)), names = estimate)
70+
times <- seq(0, 200, 1)
71+
(g*x*p)(times, pars) %>% plot(data = data, time <= 180)
72+
73+
obj <- normL2(data, g*x*p) + constraintL2(pars, sigma = 10)
74+
myfit <- trust(obj, pars, rinit = 1, rmax = 10, iterlim = 500)
75+
76+
(g*x*p)(times, myfit$argument) %>% plot(data = data)
77+
78+
profiles <- profile(obj, myfit$argument, names(myfit$argument), cores = 4)
79+
profiles %>% plotProfile(mode == "data")
80+
81+
82+
# Add additional data
83+
data$exp1_30_no_0_0_30 <- rbind(data$exp1_30_no_0_0_30, data.frame(
84+
name = "Tca_buffer", time = 0, value = 21.5, sigma = 2
85+
))
86+
87+
obj <- normL2(data, g*x*p) + constraintL2(pars, sigma = 10)
88+
myfit <- trust(obj, pars, rinit = 1, rmax = 10, iterlim = 500)
89+
90+
(g*x*p)(times, myfit$argument) %>% plot(data = data)
91+
92+
profiles <- profile(obj, myfit$argument, names(myfit$argument), cores = 4)
93+
profiles %>% plotProfile(mode == "data")
94+
95+
# Explore parameter space
96+
fits <- mstrust(obj, pars, rinit = 1, rmax = 10, samplefun = "runif", min = -5, max = 5, cores = 4)
97+
parframe <- as.parframe(fits)
98+
subframe <- unique(parframe[1:19,])
99+
100+
prediction <- predict(g*x*p, times = times, pars = subframe, data = data)
101+
102+
ggplot(prediction, aes(x = time, y = value, color = as.factor(.value))) + facet_wrap(~name*cations, scales = "free") +
103+
geom_line() +
104+
geom_point(data = as.data.frame(data), color = "black")
105+
106+
profiles <- lapply(1:nrow(subframe), function(i) {
107+
profile(obj, as.parvec(subframe, i), names(myfit$argument), cores = 4, algoControl = list(gamma = 10))
108+
})
109+
names(profiles) <- c("obj = 23", "obj = 24", "obj = 27")
110+
profiles %>% plotProfile(mode == "data")
111+
112+
# Save objects from this session
113+
save(g, x, p, transformation, data, file = "part1.RData")
114+
115+

‎vignettes/BA_publication_part1.R‎

Lines changed: 122 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,122 @@
1+
2+
library(dMod)
3+
library(dplyr)
4+
library(ggplot2)
5+
theme_set(theme_dMod())
6+
rm(list = ls())
7+
8+
## Read the data ---------------------------------------------
9+
10+
data(BAdata)
11+
data <- BAdata %>%
12+
filter(experiment == "exp1") %>%
13+
as.datalist()
14+
covtable <- covariates(data)
15+
print(covtable)
16+
17+
plot(data)
18+
19+
20+
## Set up a base model ----------------------------------------
21+
22+
# Set up base reactions
23+
reactions <- eqnlist() %>%
24+
addReaction("Tca_buffer", "Tca_cyto", "import_Tca*Tca_buffer", "Basolateral uptake") %>%
25+
addReaction("Tca_cyto", "Tca_buffer", "export_Tca_baso*Tca_cyto", "Basolateral efflux") %>%
26+
addReaction("Tca_cyto", "Tca_canalicular", "export_Tca_cana*Tca_cyto", "Canalicular efflux") %>%
27+
addReaction("Tca_canalicular", "Tca_buffer", "transport_Tca*Tca_canalicular", "Transport bile") %>%
28+
addReaction("0", "cations", "0", "Cation concentration in buffer")
29+
30+
# Translate to ODEs and let cation concentration change transport rate
31+
odes <- reactions %>%
32+
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")
45+
x <- odemodel(odes, events = events, fixed = noSens, modelname = "BA_basemodel") %>%
46+
Xs()
47+
48+
# Make observation function available
49+
g <- eqnvec(
50+
buffer = "s*Tca_buffer",
51+
cellular = "s*(Tca_cyto + Tca_canalicular)") %>%
52+
Y(f = x, modelname = "obsfn", compile = TRUE)
53+
54+
# Parameterize the model
55+
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+
61+
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+
}
70+
71+
estimate <- getParameters(p)
72+
pars <- structure(rep(-1, length(estimate)), names = estimate)
73+
times <- seq(0, 200, 1)
74+
(g*x*p)(times, pars) %>% plot(data = data, time <= 180)
75+
76+
obj <- normL2(data, g*x*p) + constraintL2(pars, sigma = 10)
77+
myfit <- trust(obj, pars, rinit = 1, rmax = 10, iterlim = 500)
78+
79+
(g*x*p)(times, myfit$argument) %>% plot(data = data)
80+
81+
profiles <- profile(obj, myfit$argument, names(myfit$argument), cores = 4)
82+
profiles %>% plotProfile(mode == "data")
83+
84+
85+
# 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+
))
89+
90+
obj <- normL2(data, g*x*p) + constraintL2(pars, sigma = 10)
91+
myfit <- trust(obj, pars, rinit = 1, rmax = 10, iterlim = 500)
92+
93+
(g*x*p)(times, myfit$argument) %>% plot(data = data)
94+
95+
profiles <- profile(obj, myfit$argument, c("s", "Tca_buffer"), cores = 2)
96+
profiles %>% plotProfile(mode == "data")
97+
98+
# Explore parameter space
99+
fits <- mstrust(obj, pars, fits = 30, rinit = 1, rmax = 10, samplefun = "runif", min = -5, max = 5, cores = 4)
100+
parframe <- as.parframe(fits)
101+
plotValues(parframe[1:29,])
102+
plotPars(parframe[1:29,])
103+
subframe <- unique(parframe[1:29,])
104+
105+
prediction <- predict(g*x*p, times = times, pars = subframe, data = data)
106+
107+
ggplot(prediction, aes(x = time, y = value, color = as.factor(.value))) + facet_wrap(~name*cations, scales = "free") +
108+
geom_line() +
109+
geom_point(data = as.data.frame(data), color = "black")
110+
111+
profiles <- lapply(1:nrow(subframe), function(i) {
112+
profile(obj, as.parvec(subframe, i), names(myfit$argument), cores = 4, algoControl = list(gamma = 10))
113+
})
114+
names(profiles) <- c("obj = 23", "obj = 24", "obj = 27")
115+
profiles %>% plotProfile(mode == "data")
116+
117+
partable <- confint(profiles[[1]], level = 0.9, val.column = "value")
118+
119+
# Save objects from this session
120+
save.image(file = "part1.RData")
121+
122+

‎vignettes/BA_publication_part2.R‎

Lines changed: 66 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,66 @@
1+
2+
library(dMod)
3+
library(dplyr)
4+
library(ggplot2)
5+
theme_set(theme_dMod())
6+
rm(list = ls())
7+
8+
## Load contents of session 1 --------------------------------
9+
10+
load("part1.RData")
11+
loadDLL(g, x)
12+
13+
## Read the data ---------------------------------------------
14+
15+
data(BAdata)
16+
data <- data + BAdata %>%
17+
filter(experiment == "exp3" & compound == "no") %>%
18+
as.datalist()
19+
covtable <- covariates(data)
20+
print(covtable)
21+
22+
plot(data)
23+
24+
25+
## Add new condition to parameterization --------------------------
26+
27+
28+
# Parameterize the model
29+
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)
45+
pars <- structure(rep(-1, length(estimate)), names = estimate)
46+
pars[partable$name] <- partable$value
47+
48+
times <- seq(0, 200, 1)
49+
(g*x*p)(times, pars) %>% plot(data = data, time <= 180)
50+
51+
obj <- normL2(data, g*x*p) + constraintL2(pars, sigma = 10)
52+
myfit <- trust(obj, pars, rinit = 1, rmax = 10, iterlim = 500)
53+
54+
(g*x*p)(times, myfit$argument) %>% plot(data = data)
55+
56+
profiles_part2 <- profile(obj, myfit$argument, names(myfit$argument), cores = 4)
57+
profiles_part2 %>% plotProfile(mode == "data")
58+
59+
plotProfile(list(old = profiles[[1]], new = profiles_part2), mode == "data")
60+
61+
partable <- confint(profiles_part2, level = 0.68)
62+
63+
# Save objects from this session
64+
save.image(file = "part2.RData")
65+
66+

0 commit comments

Comments
 (0)