Model prediction function for ODE models.

Share:

Description

Interface to combine an ODE and its sensitivity equations into one model function x(times, pars, deriv = TRUE) returning ODE output and sensitivities.

Usage

1
2
3
Xs(odemodel, forcings = NULL, events = NULL, names = NULL,
  condition = NULL, optionsOde = list(method = "lsoda"),
  optionsSens = list(method = "lsodes"))

Arguments

odemodel

object of class odemodel

forcings

data.frame with columns name (factor), time (numeric) and value (numeric). The ODE forcings.

events

data.frame of events with columns "var" (character, the name of the state to be affected), "time" (numeric, time point), "value" (numeric, value), "method" (character, either "replace", "add" or "multiply"). See events. Within Xs() a data.frame of additional events is generated to reset the sensitivities appropriately, depending on the event method. ATTENTION: The addional events are not dynamically recalculated. If you call the prediction function with alternative events, the prediction is fine but the sensitivities can be wrong.

names

character vector with the states to be returned. If NULL, all states are returned.

condition

either NULL (generic prediction for any condition) or a character, denoting the condition for which the function makes a prediction.

optionsOde

list with arguments to be passed to odeC() for the ODE integration.

optionsSens

list with arguments to be passed to odeC() for integration of the extended system

Value

Object of class prdfn. If the function is called with parameters that result from a parameter transformation (see P), the Jacobian of the parameter transformation and the sensitivities of the ODE are multiplied according to the chain rule for differentiation. The result is saved in the attributed "deriv", i.e. in this case the attibutes "deriv" and "sensitivities" do not coincide.

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
## Not run: 

################################################################  
## Walk through the most frequently used functions in
## connection with ODE models
################################################################  
  
library(deSolve)

## Model definition (text-based, scripting part)
r <- NULL
r <- addReaction(r, "A", "B", "k1*A", "translation")
r <- addReaction(r, "B",  "", "k2*B", "degradation")

f <- as.eqnvec(r)
observables <- eqnvec(Bobs = "s1*B")

innerpars <- getSymbols(c(names(f), f, observables))
trafo0 <- structure(innerpars, names = innerpars)
trafo0 <- replaceSymbols(innerpars, 
                         paste0("exp(", innerpars, ")"), 
                         trafo0)

conditions <- c("a", "b")

trafo <- list()
trafo$a <- replaceSymbols("s1", "s1_a", trafo0)
trafo$b <- replaceSymbols("s1", "s1_b", trafo0)

events <- list()
events$a <- data.frame(var = "A", time = 5, value = 1, method = "add")
events$b <- data.frame(var = "A", time = 5, value = 2, method = "add")


## Model definition (generate compiled objects 
## and functions from above information) 

# ODE model
model <- odemodel(f)

# Observation function
g <- Y(observables, f, compile = TRUE, modelname = "obsfn")

# Parameter transformation for steady states
pSS <- P(f, method = "implicit", condition = NULL)

# Condition-specific transformation and prediction functions
p0 <- x <- NULL
for (C in conditions) {
  p0 <- p0 + P(trafo[[C]], condition = C)
  x <- x + Xs(model, events = events[[C]], condition = C)
}

## Process data

# Data
data <- datalist(
  a = data.frame(time = c(0, 2, 7),
                 name = "Bobs",
                 value = c(.3, .3, .3),
                 sigma = c(.03, .03, .03)),
  b = data.frame(time = c(0, 2, 7),
                 name = "Bobs",
                 value = c(.1, .1, .2),
                 sigma = c(.01, .01, .02))
)

## Evaluate model/data

# Initialize times and parameters
times <- seq(0, 10, .1)
outerpars <- attr(p0, "parameters")
pouter <- structure(rnorm(length(outerpars)), 
                    names = outerpars)

# Plot prediction
plot((g*x*p0)(times, pouter))
plot((g*x*pSS*p0)(times, pouter))
plotFluxes(pouter, g*x*p0, times, getFluxes(r)$B)

# Fit data with L2 constraint
constr <- constraintL2(mu = 0*pouter, sigma = 5)
myfit <- trust(normL2(data, g*x*p0) + constr, 
               pouter, rinit = 1, rmax = 10)
plot((g*x*p0)(times, myfit$argument), data)

# List of fits
obj <- normL2(data, g*x*p0) + constr
mylist <- mstrust(normL2(data, g*x*p0) + constr, 
                  pouter, fits = 10, cores = 4, sd = 4)
plot(mylist)
pars <- as.parframe(mylist)
plotValues(pars)

bestfit <- as.parvec(pars)

# Compute profiles of the fit
profile <- profile(normL2(data, g*x*p0) + constr, 
                   bestfit, "k1", limits = c(-5, 5))
plotProfile(profile)

# Add a validation objective function
vali <- datapointL2("A", 2, "mypoint", .1, condition = "a")
# Get validation point parameters
mypoint <- trust(normL2(data, g*x*p0) + constr + vali, 
                 c(mypoint = 0), rinit = 1, rmax = 10, 
                 fixed = bestfit)$argument
# Compoute profile and plot
profile <- profile(normL2(data, g*x*p0) + constr + vali, 
                   c(bestfit, mypoint), "mypoint")
plotProfile(profile)


## Using the controls() function to manipulate objects ------------------------

# List the available controls of an object
controls(x)

# List a specific control
controls(x, condition = "a", name = "optionsSens")

# Change specific control
controls(x, condition = "a", name = "optionsSens") <- 
  list(method = "lsoda", rtol = 1e-8, atol = 1e-8)

# Almost every function (objfn, parfn, prdfn, obsfn) saves the arguments 
# by which it was invoked in a list named "controls", which can be manipulated
# by the controls() function.


## New example: condition-specific observation parameters ---------------------

# Condition-specific observation parameters
f <- eqnvec(A = "-k1*A", B = "k1*A - k2*B")
observables <- eqnvec(Bobs = "s1*B")
conditions <- c("scale1", "scale2")

dynpars <- getSymbols(c(names(f), f))
obspars <- setdiff(getSymbols(observables), dynpars)
trafo0 <- structure(c(dynpars, obspars), names = c(dynpars, obspars))

obspars_all <- outer(obspars, conditions, paste, sep = "_")
trafo_all <- structure(c(dynpars, obspars_all), 
                       names = c(dynpars, obspars_all))
trafo_all <- replaceSymbols(dynpars, 
                            paste0("exp(", dynpars, ")"), 
                            trafo_all)
trafo_all <- replaceSymbols(obspars_all, 
                            paste0("exp(", obspars_all, ")"), 
                            trafo_all)

pObs <- NULL
for (C in conditions) {
  pObs <- pObs + P(replaceSymbols(obspars, 
                                  paste(obspars, C, sep = "_"), 
                                  trafo0), 
                   condition = C)
}


x <- Xs(model)
p <- P(trafo_all)
y <- g*pObs*x*p

times <- seq(0, 10, .1)
outerpars <- attr(p, "parameters")
pouter <- structure(rnorm(length(outerpars)), names = outerpars)

# Plot prediction
plot((g*pObs*x*p)(times, pouter))



## End(Not run)

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker.