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 (25-06-2024) was sourced

Grind 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.4999936

which 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.4798607 0.5224990

Likewise,

run(state=c(R=0.5,N=0.5),after="state<-state+rnorm(2,mean=0,sd=0.05)",ymax=1)

##         R         N 
## 0.3903761 0.6660798

runs 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 0

Arrest 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 0

See 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.5

For 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.1294354 0.2368022 0.2429884 0.3081778 0.2775318 0.2636773

Note 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  <- c(N[n], N[1:(n-1)])
    Nright <- c(N[2: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){
  emut <- parms["mu"]*state         # Expected number of mutants
  nmut <- ifelse(emut>1,emut,ifelse(runif(n)<emut,1,0))
  left <- c(0,nmut[1:(n-1)])
  state <- state - nmut + left
  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 
## 5.315495e-09 2.283427e-06 3.970162e-04 2.312944e-02 6.491735e-01 1.170369e+01 
##           N7           N8           N9          N10 
## 8.843790e+01 4.711609e+02 1.551049e+03 2.615319e+03

Maps

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.2729081

Make 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.5084769

where 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().

More information

These are just a few examples of solving models by numerical integration. The manual explains all options that can be given to run(). Check Grind’s homepage for other examples.