Numerical integration in Grind
Rob J. de Boer
2024-03-22
This document was made in the RStudio environment, and expects that
you have installed the deSolve, rootSolve,
coda and FME R-packages beforehand. This can
be done by Install Packages in the Tools menu
of RStudio.
To start a session in Grind first source the grind.R
script:
source("grind.R")## grind.R (15-03-2024) was sourcedGrind needs to sourced to define its functions. This document only
uses the function run() to explain a few of its
features.
The Lotka Volterra model in Grind
Consider the Lotka Volterra model \(R'=rR(1-R/K)-aRN\) and \(N'=caRN-dN\) with parameter values,
\(r=1,K=1,a=1,c=1\), and \(d=0.5\), and the initial condition \(R=1\) and \(N=0.01\). Using Grind’s default names for
the model, parameters, and state, i.e., model(),
p, and s, this should be written as
model <- function(t, state, parms) {
with(as.list(c(state,parms)), {
dR <- r*R*(1 - R/K) - a*R*N
dN <- c*a*R*N - d*N
return(list(c(dR, dN)))
})
}
p <- c(r=1,K=1,a=1,c=1,d=0.5)
s <- c(R=1,N=0.01)See the basic tutorial for more explanation of the syntax of this function.
Numerical integration is performed with run()
One can simulate the model for 50 time steps with output generated at every 10 time steps by calling
run(50,10,table=TRUE)
## time R N
## 1 0 1.0000000 0.0100000
## 2 10 0.6332281 0.4720938
## 3 20 0.5002436 0.4911588
## 4 30 0.4992284 0.5006415
## 5 40 0.5000431 0.5000221
## 6 50 0.5000031 0.4999936which makes a timeplot and returns a dataframe containing the output.
Adding events and noise during numerical integration
Run() has a simple option after to handle events, change
parameters, or variables after each time step. For instance one can draw
a random value of the carrying capacity by the option after, which is a
text string that will evaluated after each time step.
run(after="parms[\"K\"]<-rnorm(1,mean=1,sd=0.1)")
## R N
## 0.4634025 0.4660910Likewise,
run(state=c(R=0.5,N=0.5),after="state<-state+rnorm(2,mean=0,sd=0.05)",ymax=1)
## R N
## 0.5667427 0.4735298runs the model starting at the steady state \((R,N)=(0.5,0.5)\) while adding randomly
distributed values around zero to the state after each time step. Note
that the state is called state, and the parameters are
called parms, when accessing them via the after string.
One can also use after with arrest to define events at
specific time points. Arrest halts the integration at predefined points,
e.g.
run(50,arrest=33.14,after="if(t==33.14)state[\"N\"]<-0")
## R N
## 1 0Arrest can also be a vector of parameters,
e.g. arrest=c("t1","t2"). Arrest should also be used to
halt the integrator at time points at which something discontinuous
happens within a model (e.g., due to an ifelse(t<t1,..)
statement).
The deSolve library that Grind uses also has its own option to handle events changing the state:
fun <- function(t, state, parms) {
state["N"] <- 0
return(state)
}
run(events=list(func=fun,time=33.14))## Warning in checkevents(events, times, Ynames, dllname): Not all event times
## 'events$time' are in output 'times' so they are automatically included.
## R N
## 1 0See the manual of deSolve. Note that with after can one call functions returning either a new state or new parameters.
Tweaking the output of run()
The option tweak can be used to change the dataframe
provided by the numerical solution. This can be useful for plotting or
fitting additional observables. For instance, one can add a column
containing the total \(T=R+N\):
run(tweak="nsol$T<-nsol$R+nsol$N")
## R N
## 0.5 0.5For large models it may be useful to use R’s function
sum() for this, e.g.
run(tweak="nsol$T<-apply(nsol[,2:3],1,sum))
Vectors of equations
Since R is a parallel language handling vectors like scalers, one can easily define models based upon vectors of equations. Consider a model with \(n\) prey populations, \(R_i\), that are competing with each other via a logistic term. Each prey has its own predator, \(N_i\). For example, consider \(R_i'=b_iR_i(1-\sum R_j) - d_1R_i - aR_iN_i\) and \(N_i'= aR_iN_i -d_2N_i\), and write that as
model <- function(t, state, parms){
with(as.list(c(state,parms)),{
R <- state[1:n]
N <- state[(n+1):(2*n)]
S <- sum(R)
dR <- b*R*(1-S) - d1*R - a*R*N
dN <- a*R*N - d2*N
return(list(c(dR,dN)))
})
}Next define a 3-dimensional system. Draw random prey birth rates, \(b_i\), from a normal distribution, set all other parameters, and define \(R_i=0.1/n\) \(N_i=0.01/n\) as the initial condition:
n <- 3
b <- rnorm(n,mean=1,sd=0.1) # b is a global parameter
p <- c(d1=0.1,d2=0.2,a=1) # define all other parameters
R <- rep(0.1/n,n) # initial condition of R
names(R) <- paste("R",seq(1,n),sep="") # Name them R1, R2, ... Rn
N <- rep(0.01/n,n) # initial condition of N
names(N) <- paste("N",seq(1,n),sep="") # Name them N1, N2, ... Nn
s <- c(R,N) # combine R and N into s
run()
## R1 R2 R3 N1 N2 N3
## 0.34400497 0.08738242 0.15794242 0.38934263 0.33864169 0.12016641Note that b is a global parameter and that the other
parameters are passed on via the parameter vector p. The
variables are named R1, R2, R3, N1, N2 and N3.
The Fisher equation as an example involving diffusion
To spatially model a population that is growing logistically on a one-dimensional spatial domain one can just add a diffusion term, and write a partial differential equation (PDE) like, \[\frac{\partial N}{\partial t}=rN\Bigl(1-\frac{N}{K}\Bigr) + D\frac{\partial^2N}{\partial x^2} \ ,\] where the parameter \(D\) is a diffusion constant. To study such a PDE numerically one needs to discretize it into something like \[N_i'=rN_i\Bigl(1-\frac{N_i}{K}\Bigr) + D(N_{i-1}+N_{i+1}-2N_i)\ ,\quad \hbox{for}\ i=1,2,\dots,n\ \] where \(i\) defines a location of a (small) compartment in space, and \(D\) describes the flow of individuals between neighboring compartments. Note that \(D\) should be sufficiently small, i.e., one needs to define many small compartments to accurately describe the PDE. Wrapping the space, this can be defined as follows
model <- function(t, state, parms) {
with(as.list(c(state,parms)), {
N <- state[1:n]
Nleft <- rep(0,n); Nright <- rep(0,n);
Nleft[2:n] <- N[1:(n-1)]; Nleft[1] <- N[n]
Nright[1:(n-1)] <- N[2:n]; Nright[n] <- N[1]
dtN <- r*N*(1-N/K) + D*(Nleft+Nright-2*N)
return(list(dtN))
})
} Define \(n=100\) compartments and start with \(N_{25}=1\):
n <- 100
N <- rep(0,n)
names(N) <- paste("N",seq(1,n),sep="")
N[25] <- 1
s <- N
p <- c(r=1,K=1,D=0.1)
f <- run(10,legend=FALSE)
plot(seq(n),f,xlab="Location",ylab="Density",main="Distribution at t=10")
Mutations after each time step
Combining vectors and events one can model a series of evolving
bacterial strains with increasing replication rates. Consider the
following model, \(N_i'=b_iN_i(1-\sum N_j)
- dN_i\) for \(i=1,2,\dots,n\).
Let strain \(N_{i+1}\) evolve from
strain \(N_i\) at a mutation rate \(\mu\). When the expected number of mutants,
\(\mu N_i\), is smaller than a single
bacterium we calculate the probability that a single mutant appears, and
add a single cell to \(N_{i+1}\) (and
subtract it from \(N_i\)). Otherwise
the expected number of mutants is added to strain \(N_{i+1}\) (and subtracted from \(N_i\)). This is realized by using the
option after to call a function mutate().
model <- function(t, state, parms){
with(as.list(c(state,parms)),{
N <- state[1:n]
S <- sum(N)
dN <- b*N*(1-S/K) - d*N
return(list(c(dN)))
})
}
mutate <- function(t, state, parms){
nmut <- rep(0,n+1)
emut <- parms["mu"]*state # Expected number of mutants
nmut[2:(n+1)] <- ifelse(emut>1,emut,ifelse(runif(n)<emut,1,0))
state <- state + nmut[1:n] - nmut[2:(n+1)]
return(state)
}Consider \(n=10\) variants and let the birth rates \(b_i\) of the strains increase linearly from \(b_i=1\) to \(b_n=2\):
n <- 10
b <- seq(1,2,length=n)
p <- c(d=0.1,K=5000,mu=0.001)
s <- rep(0,n); s[1] <- 1
names(s) <- paste("N",seq(1,n),sep="")
run(1000,ymax=5000,after="state<-mutate(t,state,parms)")
## N1 N2 N3 N4 N5 N6
## 4.355982e-09 1.208047e-06 1.025996e-04 9.391411e-03 2.622610e-01 5.126638e+00
## N7 N8 N9 N10
## 5.782755e+01 3.174352e+02 1.345584e+03 3.013610e+03Maps
One can study maps by simply switching to Euler integration. For instance, the Logistic map, \(N_{t+1}=rN_t(1-N_t)\), can be studied with:
model <- function(t, state, parms) {
with(as.list(c(state,parms)), {
dN <- r*N*(1 - N) - N
return(list(c(dN)))
})
}
p <- c(r=3.75)
s <- c(N=0.01)
run(100,method="euler")
## N
## 0.2729081Make the famous bifurcation diagram by running for various values of
r Allow for a transient “burnin” time of 500 steps and plot
the next 100 steps:
plot(range(2,4),range(0,1),type='n',xlab="r",ylab="N")
npoints <- 100; burnin <- 500
for (r in seq(2,4,0.01)) {
p["r"] <- r
data <- run(burnin+npoints,method="euler",table=TRUE,timeplot=FALSE)
points(rep(r,npoints),data$N[(burnin+1):(burnin+npoints)],pch=".")
}
Make a Takens reconstruction, i.e., make a time series by running the model and plot \(N_{t+1}\) as a function of \(N_t\):
p <- c(r=3.75)
s <- c(N=0.01)
data <- run(1000,method="euler",table=TRUE,timeplot=FALSE)
plot(data$N[1:999],data$N[2:1000],pch=".")
Delay differential equations (DDEs)
One can study DDEs because the deSolve package
implements a general solver dede() using the same syntax as
its general ODE solver ode(). For instance, a Lotka
Volterra model with delayed growth of the predator, \(R(t)'=rR(t)(1-R(t)/K)-aR(t)N(t)\) and
\(N(t)'=caR(t-\Delta)N(t-\Delta)-dN(t)\)
would look like:
model <- function(t, state, parms) {
with(as.list(c(state,parms)), {
tlag <- t - Delta
if (tlag < 0) lags <- c(0,0) # no initial predation
else lags <- lagvalue(tlag) # returns lags of R and N
dR <- r*R*(1 - R/K) - a*R*N
dN <- a*lags[1]*lags[2] - d*N
return(list(c(dR, dN)))
})
}
p <- c(r=1,K=1,a=1,c=1,d=0.5,Delta=10)
s <- c(R=1,N=0.1)
run(delay=TRUE)
## R N
## 0.4951771 0.5084769where the option delay tells Grind to use the
dede() solver. This would correspond to the situation where
a prey at carrying capacity starts to be eaten by a predator at time
zero, but the predator will only start to grow \(\Delta\) time steps later. The function
lagvalue() stores previous values of \(R\) and \(N\) in a vector (called lags[]
in this example). This vector can be indexed to obtain the time lagged
\(R_{t-\Delta}\) and \(N_{t-\Delta}\) values. Read the manual for
further documentation and/or type ?dede for help. Grind’s
steady state functions, newton(), continue(), and
plane() should not be called with DDE models containing
calls to lagvalue().