Generating and Plotting Data from ODE Models: Part 1
Overview
Our work has taken us into the realm of Ordinary Differential Equations. This tutorial is an exploration into how to plot output from these models - both for simulations (to increase understanding of how the models work), and for reporting of empirical results. Much is borrowed directly from http://desolve.r-forge.r-project.org/slides/tutorial.pdf.
Outline
In this session we cover …
A. Analytic Solution of Differential Equation
B. Numerical Solution of Differential Equation (4 Steps)
C. Expansion to Multiple Scenarios (Persons)
D. Oscillator Model Example (VanDerPol Equation)
E. Oscillator Model Example (Univariate Linear Oscillator)
F. Bivariate System - Predator Prey Model
G. Bivariate System - Linear Oscillator Model
Prelim - Loading libraries used in this script.
library(deSolve) #for the differential equation solvers
library(ggplot2) #for plottingA. Analytic Solution of Differential Equation
Let’s start with a very simple model - the logistic growth model. In differential form, the model is … (using the notation often used in ecology) …
\[\frac{dN}{dt} = rN\left( \frac{1 - N}{K} \right)\]
The analytical solution is derived by integrating and is … \[ N = \frac{KN_0exp(rt)} {(K + N_0(exp(rt) - 1))}\]
(see Ram & Grimm, 2015 for discussion from a growth curve modeling perspective).
We can make a function of the integrated “solution” and plot a simulated curve.
logistic <- function(t, r, K, N0){
K * N0 * exp(r*t) / (K + N0 * (exp(r*t) - 1))
}
plot(0:100, logistic(t = 0:100, r = 0.1, K = 10, N0 = 0.1))
B. Numerical Solution of Differential Equation (4 Steps)
We can also use the deSolve library to obtain a numerical solution. Four steps are needed.
1. Write the differential equation as a function.
difflogistic <- function(time, y, parms){
with(as.list(c(y,parms)),{
dN <- r * N * (1 - N / K)
list(dN)
})
}2. Specify vectors with the initial values, model parameters, and times over which to integrate.
#Initial Values
initvalues.1 <- c(N = 0.1)
#Model Parameters
parms.1 <- c(r = 0.1, K = 10)
#Times over which to integrate
times <- seq(0, 100, 1)3. Solve the differential equation using numerical methods.
This is done using the ode() function, which is itself a wrapper of the lsoda() function (and a collection of similar functions for different types of integrations (e.g., Runge-Katta), some of which work better with particular types of models (e.g., stiff and non-stiff systems)).
out <- ode(y = initvalues.1, times = times, func = difflogistic, parms = parms.1) This produces an object of class deSolve that can then be examined, summarized, used for diagnosis, and plotted.
head(out)## time N
## [1,] 0 0.1000000
## [2,] 1 0.1104022
## [3,] 2 0.1218708
## [4,] 3 0.1345160
## [5,] 4 0.1484538
## [6,] 5 0.1638111summary(out)## N
## Min. 0.100000
## 1st Qu. 1.095729
## Median 5.998616
## Mean 5.395638
## 3rd Qu. 9.480876
## Max. 9.955256
## N 101.000000
## sd 3.902511diagnostics(out)##
## --------------------
## lsoda return code
## --------------------
##
## return code (idid) = 2
## Integration was successful.
##
## --------------------
## INTEGER values
## --------------------
##
## 1 The return code : 2
## 2 The number of steps taken for the problem so far: 105
## 3 The number of function evaluations for the problem so far: 211
## 5 The method order last used (successfully): 5
## 6 The order of the method to be attempted on the next step: 5
## 7 If return flag =-4,-5: the largest component in error vector 0
## 8 The length of the real work array actually required: 36
## 9 The length of the integer work array actually required: 21
## 14 The number of Jacobian evaluations and LU decompositions so far: 0
## 15 The method indicator for the last succesful step,
## 1=adams (nonstiff), 2= bdf (stiff): 1
## 16 The current method indicator to be attempted on the next step,
## 1=adams (nonstiff), 2= bdf (stiff): 1
##
## --------------------
## RSTATE values
## --------------------
##
## 1 The step size in t last used (successfully): 1
## 2 The step size to be attempted on the next step: 1
## 3 The current value of the independent variable which the solver has reached: 100.8645
## 4 Tolerance scale factor > 1.0 computed when requesting too much accuracy: 0
## 5 The value of t at the time of the last method switch, if any: 0
## These diagnostics output may be useful for checking the kind of integration method used (ode() chooses on its own), and if there were problems. (I have not run enough of these yet to know how, specifically, to use the diagnostics)
4. Plot the output trajectory (solution).
plot(out, main="Logistic Growth", lwd=2, xlab="Time", ylab="N")
Given our proclivity for plotting using ggplot2 function, we can also make the plot there, after conversion of the deSolve object into a data frame.
outplot=as.data.frame(out)
ggplot(data=outplot, aes(x=time, y=N)) +
geom_point() +
geom_line() +
xlab("Time") +
ylab("N") +
ggtitle("Logistic Growth")
C. Expansion to Multiple Scenarios (Persons)
We can expand by plotting multiple scenarios (e.g., simualted trajectories for a set of selected or randomly chosen parameters (or starting values)).
#Use a list if many unknown number of outputs
outlist <- list()
plist <- cbind(r = runif(30,min=0.1,max=5),
K = runif(30,min=8,max=15))
#Loop through the matrix of parameter values
for(i in 1:nrow(plist)){
outlist[[i]] <- ode(y = initvalues.1, times = times, func = difflogistic, parms = plist[i,])
}
plot(out,outlist, lwd=2)
Ok - that is all good and bascially replicates most of what we do in individual-level growth modeling, with access to a large set of interesting models - basically the first 4 chapters of Banks (1994).
D. Oscillator Model Example
OK - lets move on to the oscillator models. The classic example used in the deSolve vignettes is the Van der Pol oscillator …
\[\frac{d^2y}{dt^2} = \mu(1-y^2)\frac{dy}{dt} - y\]
A key difference with the logistic model is that this model is a 2nd order ODE. It must be transformed into two 1st order equations …. \[\frac{dy_1}{dt} = y_2\] \[\frac{dy_2}{dt} = \mu(1-y_{1}^2)y_{2} - y_{1}\] Then we can implment in a similar way with our 4 steps.
1. Write the model as a function
#using code from vignette
model.vdpol <- function(time, y, parms){
with(as.list(c(y,parms)),{
list(c(
y[2],
mu * (1-y[1]^2) * y[2] - y[1]
))
})
}2. setting initial values, model parameters, times
initvalues.vdpol <- c(y1 = 2, y2 = 0)
parms.vdpol <- c(mu = 1000)
times <- seq(0, 3000, 1)3. Solving to get trajectory
stiff <- ode(y = initvalues.vdpol, times=times, func=model.vdpol, parms=parms.vdpol)
diagnostics(stiff)##
## --------------------
## lsoda return code
## --------------------
##
## return code (idid) = 2
## Integration was successful.
##
## --------------------
## INTEGER values
## --------------------
##
## 1 The return code : 2
## 2 The number of steps taken for the problem so far: 3970
## 3 The number of function evaluations for the problem so far: 5310
## 5 The method order last used (successfully): 4
## 6 The order of the method to be attempted on the next step: 4
## 7 If return flag =-4,-5: the largest component in error vector 0
## 8 The length of the real work array actually required: 52
## 9 The length of the integer work array actually required: 22
## 14 The number of Jacobian evaluations and LU decompositions so far: 232
## 15 The method indicator for the last succesful step,
## 1=adams (nonstiff), 2= bdf (stiff): 2
## 16 The current method indicator to be attempted on the next step,
## 1=adams (nonstiff), 2= bdf (stiff): 2
##
## --------------------
## RSTATE values
## --------------------
##
## 1 The step size in t last used (successfully): 1
## 2 The step size to be attempted on the next step: 1
## 3 The current value of the independent variable which the solver has reached: 3000.728
## 4 Tolerance scale factor > 1.0 computed when requesting too much accuracy: 0
## 5 The value of t at the time of the last method switch, if any: 2421.491
## 4. Plot the trajectory (note our interest is in y1)
plot(stiff, type="l", which = "y1", lwd=2,ylab="y", main = "van der Pol oscillator, mu = 1000")
The output has been called stiff as a note that there are some nuances in the background about “stiff” and “non-stiff” systems. If one changes the parmaters one can see that ode() chooses a different integration method when mu = 1000 (a “stiff” case) and when mu = 1 (a “non-stiff” case)). Models are “stiff” when there are flat areas and fast change areas (a harder problem to solve). Models are “non-stiff” when there is similar change across the whole series (an easier problem to solve).
The coding here for how the initial state vector and parameters are brought into the model function are not as elegantly matched to the equations as they might be. So, in the below, we do a bit of recoding to get a nicer naming structure in the function for the 2nd order system that matches the equations a bit better. Still, care must be taken with consistency in naming across the model function (model.abcd), the initial values vector (initvalues.abcd), and the parameters vector (parms.abcd).
1. Write the model as a function
In writing out the function we follow the equations above, except that we use x instead of y. This is because deSolve uses y = for the vector of initial values, and we did not want to get confused with the embedded object naming.
#with closer notation match (but using x)
model.vdpol = function(times, y, parms) {
## t = current time
## y = current state vector: in this case the 0th and 1st derivatives, c(x, dx)
## parms = vector of model parameters.
with(as.list(c(y,parms)),{
derivs = c(
dx = dx,
d2x = mu*(1-x^2)*dx - x
)
## Return
list(derivs)
})
}2. setting initial values, model parameters, times
initvalues.vdpol <- c(x = 2, dx = 0)
parms.vdpol <- c(mu = 5)
times <- seq(0, 20, .2)3. Solving to get trajectory
stiff <- ode(y = initvalues.vdpol, times=times, func=model.vdpol, parms=parms.vdpol)
summary(stiff)## x dx
## Min. -2.01902696 -6.2945456
## 1st Qu. -1.68715548 -0.2225785
## Median 0.75820689 -0.1405331
## Mean 0.09955845 -0.1748458
## 3rd Qu. 1.67364437 0.1618086
## Max. 2.01818619 6.1384929
## N 101.00000000 101.0000000
## sd 1.61408535 1.39119494. Plot the trajectory
#plotting x
plot(stiff, type="l", which = "x", lwd=2,ylab="x", main = "van der Pol oscillator, mu = 5")
#Plotting both x and dx together
matplot(stiff[,1], stiff[,-1], type="l",xlab="time", ylab="x,dx",lwd=2, main = "damped linear oscillator")
legend("topright", c("x", "dx"), col = 1:2, lty = 1:2)
E. Oscillator Model Example (Univariate Linear Oscillator)
Lets adapt the above example (and the cleaned-up naming convention) for the linear oscillator model we have been using in our work.
\[\frac{d^2x}{dt^2} = \eta x + \zeta \frac{dx}{dt}\]
We implement in a similar way with our 4 steps.
1. Write the model as a function
#writing model as function
model.damposc = function(times, y, parms) {
## t = current time
## y = current state vector in this case the 0th and 1st derivatives, c(x, dx)
## parms = vector of model parameters
with(as.list(c(y,parms)),{
derivs = c(
dx = dx,
d2x = eta*x + zeta*dx
)
## Return
list(derivs)
})
}2. setting initial values, model parameters, times
initvalues.damposc <- c(x = 2, dx = 0)
parms.damposc <- c(eta=-.14, zeta = -.1)
times <- seq(0, 50, .2)3. Solving to get trajectory
out <- ode(y = initvalues.damposc, times=times, func=model.damposc, parms=parms.damposc)
summary(out)## x dx
## Min. -1.30887690 -0.616292863
## 1st Qu. -0.29864644 -0.140599628
## Median -0.01144032 -0.005469619
## Mean 0.02792894 -0.036804782
## 3rd Qu. 0.33804279 0.092591884
## Max. 2.00000000 0.403551544
## N 251.00000000 251.000000000
## sd 0.65785793 0.2330056094. Plot the trajectory (note our interest is in y1)
plot(out, type="l", which = "x", lwd=2,ylab="x", main = "damped linear oscillator")
#Plotting both x and dx together
matplot(out[,1], out[,-1], type="l",xlab="time", ylab="x,dx",lwd=2, main = "damped linear oscillator")
legend("topright", c("x", "dx"), col = 1:2, lty = 1:2)
F. Bivariate System - Predator Prey Model
Moving on to some bivariate systems …
First lets do the classic Predator Prey model - a 1st order system
1. Specify the model as a function
#Writing out model
model.LV <- function(times, y , parms) {
with(as.list(c(y, parms)), {
Ingestion <- rIng * Prey * Predator
GrowthPrey <- rGrow * Prey * (1 - Prey/K)
MortPredator <- rMort * Predator
derivs = c(
dPrey <- GrowthPrey - Ingestion,
dPredator <- Ingestion * assEff - MortPredator
)
return(list(derivs))
})
}2. setting model parameters, initial values, times
parms.LV <- c(rIng = 0.2, # /day, rate of ingestion
rGrow = 1.0, # /day, growth rate of prey
rMort = 0.2, # /day, mortality rate of predator
assEff = 0.5, # -, assimilation efficiency
K = 10) # mmol/m3, carrying capacity
initvalues.LV <- c(Prey = 1, Predator = 2)
times <- seq(0, 200, by = 1)3. Solving to get trajectory
out <- ode(y=initvalues.LV, times = times, func = model.LV, parms = parms.LV)
summary(out)## Prey Predator
## Min. 1.0000000 1.8632829
## 1st Qu. 1.9999571 3.9999115
## Median 2.0000000 4.0000000
## Mean 2.0317905 3.9606228
## 3rd Qu. 2.0000751 4.0000418
## Max. 4.2001812 4.9787222
## N 201.0000000 201.0000000
## sd 0.3138139 0.34890794. Plotting the implied trajectory (solution)
matplot(out[ , 1], out[ , 2:3], type = "l", xlab = "time", ylab = "Conc",
main = "Lotka-Volterra", lwd = 2)
legend("topright", c("Prey", "Predator"), col = 1:2, lty = 1:2)
G. Bivariate System - Linear Oscillator Model
Now lets do our 2nd order Coupled linear oscillator model
Our coupled equations are …
\[\frac{d^2x}{dt^2} = \eta_x x + \zeta_x\frac{dx}{dt} + \gamma_1z\] \[\frac{d^2z}{dt^2} = \eta_z z + \zeta_z \frac{dz}{dt} + \gamma_2 x\] We do the four steps as usual
1. specifying the model as a function
#Writing out the model as a function
model.coupleosc <- function(times, y , parms) {
with(as.list(c(y, parms)), {
derivs = c(
dx = dx,
d2x = eta_x*x + zeta_x*dx + gamma1*z,
dz = dz,
d2z = eta_z*z + zeta_z*dz + gamma2*x
)
return(list(derivs))
})
}2. setting model parameters, initial values, times
parms.coupleosc <- c(eta_x = -0.2, # frequency of x
zeta_x = -.1, # damping of x
eta_z = -0.05, # frequency of z
zeta_z = -0.3, # damping of z
gamma1 = -.1, # coupling parameter 1 (z on d2x)
gamma2 = +.1) # coupling parameter 2 (x on d2z)
initvalues.coupleosc <- c(x = 3, dx = 0, z = 2, dz = 0)
times <- seq(0, 200, by = 1)3. Solving to get trajectory
out <- ode(y=initvalues.coupleosc, times = times, func = model.coupleosc, parms = parms.coupleosc)
summary(out)## x dx z dz
## Min. -14.79261987 -5.83560949 -9.26433637 -3.60227910
## 1st Qu. -5.16354504 -1.92409340 -2.93865183 -1.22593568
## Median 0.03717895 -0.03126275 0.16523351 -0.07442915
## Mean -0.13493234 0.04198308 0.02922775 -0.05662182
## 3rd Qu. 5.03297122 2.10343092 3.00129769 1.19886562
## Max. 14.19363853 6.07324780 8.71366181 3.43320050
## N 201.00000000 201.00000000 201.00000000 201.00000000
## sd 6.53258288 2.67200730 3.97986042 1.585087034. Plotting the implied trajectory (solution)
matplot(out[ , 1], out[ , c(2,4)], type = "l", xlab = "time", ylab = "x,z",
main = "Coupled Oscillator", lwd = 2)
legend("topright", c("x", "z"), col = 1:2, lty = 1:2)
##ggplot version
#convert to data frame
outplot=as.data.frame(out)
ggplot(data=outplot, aes(x=time)) +
#geom_point() +
geom_line(aes(y=x), color="red") +
geom_line(aes(y=z), color="black") +
xlab("Time") +
ylab("x = red, z = black") +
ggtitle("Coupled Oscillator")
Yay!
We have an approach here in place to plot the results from our model fitting procedures. A next step is to figure out how to get parameters from the lme model fitting of the model straight into the ode() and plotting routine. We will look to do that in Part 2.