|
| 1 | +# The main functions of cppDE on small models: sensitivities of first and second |
| 2 | +# order, reverse-mode gradients, events, switches, forcings, steady states, the |
| 3 | +# cvode() backend and cppFUN(). |
| 4 | + |
| 5 | +library(cppDE) |
| 6 | +library(ggplot2) |
| 7 | + |
| 8 | +theme_set(theme_bw(base_size = 11)) |
| 9 | + |
| 10 | +# A trajectory and its sensitivities as one long data frame, one panel per |
| 11 | +# state and per derivative of a state. |
| 12 | +longFormat <- function(out, states = colnames(out$variable), sens = TRUE) { |
| 13 | + nt <- length(out$time) |
| 14 | + value <- data.frame(time = out$time, panel = rep(states, each = nt), |
| 15 | + value = as.vector(out$variable[, states])) |
| 16 | + if (!sens) return(value) |
| 17 | + tangent <- out$tangent[, states, , drop = FALSE] |
| 18 | + grid <- expand.grid(time = out$time, state = states, |
| 19 | + par = dimnames(tangent)[[3]], stringsAsFactors = FALSE) |
| 20 | + rbind(value, data.frame(time = grid$time, |
| 21 | + panel = sprintf("d%s/d%s", grid$state, grid$par), |
| 22 | + value = as.vector(tangent))) |
| 23 | +} |
| 24 | + |
| 25 | + |
| 26 | +# Values and sensitivities ------------------------------------------------------ |
| 27 | +# A decay chain A -> B. cppODE() generates and compiles the solver; solveODE() |
| 28 | +# returns the states and, by default, their derivatives in every initial value |
| 29 | +# and parameter. |
| 30 | +eqns <- c(A = "-k1*A", B = "k1*A - k2*B") |
| 31 | +parms <- c(A = 1, B = 0, k1 = 0.5, k2 = 0.2) |
| 32 | +times <- seq(0, 10, by = 0.1) |
| 33 | + |
| 34 | +chain <- cppODE(eqns, modelname = "chain") |
| 35 | +out <- solveODE(chain, times, parms) |
| 36 | + |
| 37 | +dim(out$tangent) # [time, state, parameter] |
| 38 | +out$tangent[101, "B", ] # dB(10) / d(A, B, k1, k2) |
| 39 | + |
| 40 | +ggplot(longFormat(out), aes(time, value)) + |
| 41 | + geom_line() + |
| 42 | + facet_wrap(~panel, scales = "free_y", ncol = 5) |
| 43 | + |
| 44 | + |
| 45 | +# Second order ----------------------------------------------------------------- |
| 46 | +# deriv2 = TRUE adds the Hessian of every state in the initial values and |
| 47 | +# parameters. |
| 48 | +chain2 <- cppODE(eqns, deriv2 = TRUE, modelname = "chain2") |
| 49 | +out2 <- solveODE(chain2, times, parms) |
| 50 | +out2$hessian[101, "B", , ] # Hessian of B(10) |
| 51 | + |
| 52 | + |
| 53 | +# Reverse mode ----------------------------------------------------------------- |
| 54 | +# derivMode = "reverse" returns the gradient of a weighted sum of the states, |
| 55 | +# sum(cotangent * x), in one backward sweep. The cotangent below weights B at |
| 56 | +# t = 10 only, so the gradient is the row dB(10) / d(A, B, k1, k2) from above. |
| 57 | +chainR <- cppODE(eqns, derivMode = "reverse", modelname = "chainR") |
| 58 | +w <- matrix(0, length(times), 2, dimnames = list(NULL, c("A", "B"))) |
| 59 | +w[101, "B"] <- 1 |
| 60 | +solveODE(chainR, times, parms, cotangent = w)$cotangent |
| 61 | + |
| 62 | +# The gradient of a least-squares objective: the cotangent holds the residuals. |
| 63 | +data <- data.frame(time = times, B = 0.4 * exp(-0.3 * times) * (1 - exp(-times))) |
| 64 | +res <- out$variable[, "B"] - data$B |
| 65 | +w[] <- 0 |
| 66 | +w[, "B"] <- 2 * res |
| 67 | +solveODE(chainR, times, parms, cotangent = w)$cotangent |
| 68 | +2 * colSums(res * out$tangent[, "B", ]) # the same from forward mode |
| 69 | + |
| 70 | + |
| 71 | +# Events ----------------------------------------------------------------------- |
| 72 | +# Events are a data frame: which state, when or on which root, which value, and |
| 73 | +# how it is applied (replace, add, multiply). The sensitivities take each jump, |
| 74 | +# also where a parameter moves the event time. |
| 75 | +events <- data.frame(var = c("A", "A"), time = c("t_e", NA), value = c("dose", "dose"), |
| 76 | + method = c("add", "add"), root = c(NA, "A - Amin")) |
| 77 | +dosed <- cppODE(eqns, events = events, modelname = "dosed") |
| 78 | +outE <- solveODE(dosed, times, c(parms, t_e = 2, dose = 1, Amin = 0.3), maxroot = 3) |
| 79 | + |
| 80 | +ggplot(longFormat(outE, "A"), aes(time, value)) + |
| 81 | + geom_line() + |
| 82 | + facet_wrap(~panel, scales = "free_y", ncol = 4) |
| 83 | + |
| 84 | + |
| 85 | +# Switches --------------------------------------------------------------------- |
| 86 | +# Comparisons and piecewise() switch the right-hand side. The solver locates the |
| 87 | +# switch as a root, so a short pulse is not stepped over. Below, a pulse of |
| 88 | +# inflow into A between t1 and t2. |
| 89 | +pulse <- cppODE(c(A = "u*(time >= t1 && time < t2) - k1*A", B = "k1*A - k2*B"), |
| 90 | + modelname = "pulse") |
| 91 | +outP <- solveODE(pulse, times, c(A = 0, B = 0, k1 = 0.5, k2 = 0.2, u = 1, t1 = 1, |
| 92 | + t2 = 1.5)) |
| 93 | + |
| 94 | +ggplot(longFormat(outP, "B"), aes(time, value)) + |
| 95 | + geom_line() + |
| 96 | + facet_wrap(~panel, scales = "free_y", ncol = 4) |
| 97 | + |
| 98 | +# A switch on a state. Above xc, x decays; below, an inflow p > k*xc pushes it |
| 99 | +# back up. Both branches move onto x = xc, so x slides along it, and from then |
| 100 | +# on its only sensitivity is the one in xc. |
| 101 | +slide <- cppODE(c(x = "-k*x + p*(x < xc)"), modelname = "slide") |
| 102 | +outS <- solveODE(slide, times, c(x = 2, k = 1, p = 3, xc = 1)) |
| 103 | + |
| 104 | +ggplot(longFormat(outS), aes(time, value)) + |
| 105 | + geom_line() + |
| 106 | + facet_wrap(~panel, scales = "free_y", ncol = 5) |
| 107 | + |
| 108 | + |
| 109 | +# Forcings --------------------------------------------------------------------- |
| 110 | +# A forcing is an input given as data when solving. Between its points it is |
| 111 | +# interpolated monotonically. |
| 112 | +forced <- cppODE(c(A = "u - k1*A", B = "k1*A - k2*B"), forcings = "u", |
| 113 | + modelname = "forced") |
| 114 | +u <- list(u = data.frame(time = seq(0, 10, by = 0.5), |
| 115 | + value = 1 + sin(seq(0, 10, by = 0.5)))) |
| 116 | +outF <- solveODE(forced, times, parms, forcings = u) |
| 117 | + |
| 118 | +ggplot(longFormat(outF, sens = FALSE), aes(time, value)) + |
| 119 | + geom_line() + |
| 120 | + facet_wrap(~panel, scales = "free_y") |
| 121 | + |
| 122 | + |
| 123 | +# Steady states ---------------------------------------------------------------- |
| 124 | +# rootfunc = "equilibrate" integrates until every derivative is negligible and |
| 125 | +# returns the steady state with its sensitivities. |
| 126 | +steady <- cppODE(c(A = "k0 - k1*A", B = "k1*A - k2*B"), rootfunc = "equilibrate", |
| 127 | + modelname = "steady") |
| 128 | +outQ <- solveODE(steady, c(0, 1e4), c(parms, k0 = 1), abstol = 1e-10, |
| 129 | + reltol = 1e-10, roottol = 1e-10) |
| 130 | +tail(outQ$time, 1) |
| 131 | +tail(outQ$variable, 1) # k0/k1 and k0/k2 |
| 132 | +outQ$tangent[nrow(outQ$variable), , ] |
| 133 | + |
| 134 | +# ptc() solves for the steady state directly, by pseudo-transient continuation, |
| 135 | +# with the conservation laws of the system as constraints. The right-hand side |
| 136 | +# is a cppFUN() object. |
| 137 | +enzyme <- c(E = "-kon*E*S + (koff + kcat)*ES", |
| 138 | + S = "-kon*E*S + koff*ES + kback*P", |
| 139 | + ES = " kon*E*S - (koff + kcat)*ES", |
| 140 | + P = " kcat*ES - kback*P") |
| 141 | +enz <- cppFUN(enzyme, variables = names(enzyme), |
| 142 | + parameters = c("kon", "koff", "kcat", "kback"), modelname = "enzyme", |
| 143 | + compile = TRUE) |
| 144 | +C <- rbind(Etot = c(E = 1, S = 0, ES = 1, P = 0), |
| 145 | + Stot = c(E = 0, S = 1, ES = 1, P = 1)) |
| 146 | +ss <- ptc(enz, x = c(E = 1, S = 1, ES = 1, P = 1), |
| 147 | + parms = c(kon = 1, koff = 0.5, kcat = 0.3, kback = 0.1), |
| 148 | + C = C, total = c(1, 10)) |
| 149 | +ss$x |
| 150 | + |
| 151 | + |
| 152 | +# The cvode() backend ----------------------------------------------------------- |
| 153 | +# The same equations on SUNDIALS CVODES, if cppDE was installed with it. |
| 154 | +if (codegenAvailable("cvode")) { |
| 155 | + chainC <- cvode(eqns, deriv = TRUE, modelname = "chainC") |
| 156 | + outC <- solveODE(chainC, times, parms, abstol = 1e-10, reltol = 1e-10) |
| 157 | + outN <- solveODE(chain, times, parms, abstol = 1e-10, reltol = 1e-10) |
| 158 | + max(abs(outC$tangent - outN$tangent)) |
| 159 | +} |
| 160 | + |
| 161 | + |
| 162 | +# Algebraic functions ---------------------------------------------------------- |
| 163 | +# cppFUN() compiles functions of variables and parameters with their Jacobian, |
| 164 | +# here an observable of the decay chain. The Jacobian is in the variables and |
| 165 | +# parameters together. |
| 166 | +g <- cppFUN(c(y = "s*(A + B)", z = "log(B)"), variables = c("A", "B"), |
| 167 | + parameters = "s", modelname = "observable", compile = TRUE) |
| 168 | +g$func(A = 0.2, B = 0.5, s = 2) |
| 169 | +g$jac(A = 0.2, B = 0.5, s = 2)[1, , ] |
0 commit comments