This document illustrates the usage of the Grind’s extension cube(), which allows one to study 3-dimensional phase spaces. The script cube.R can be downloaded from the Grind directory or via Grind’s home page. A short manual of cube() is provided on the webpage cube.html. This document was made in the RStudio environment, and expects that you have installed the plot3D, deSolve, rootSolve, coda and FME R-packages beforehand. This can be done by Install Packages in the Tools menu of RStudio.

For this example source both the grind.R`` and thecube.R` script:

source("grind.R")
## grind.R (03-05-2024) was sourced
source("cube.R")
## cube.R (03-05-2024) was sourced

Be reminded that cube() uses the plot3D R-package.

Chaos in a three-species food chain

Consider the following model proposed by Hastings and Powell (1991), with a resource, \(R\), consumer \(N\), and top-predator \(M\), i.e., \[R'= R(1 - R) - \frac{c_1NR}{1 + b_1R} \ ,\] \[N'= -a_N N + \frac{c_1NR}{1 + b_1R} - \frac{c_2MN}{1 + b_2N} \ ,\] \[M'= -a_M M + \frac{ c_2MN}{1 + b_2N} \ ,\] having logistic growth of the resource and saturated consumption terms of the consumer and the predator. This is written as

model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dR <- R*(1 - R/K) - c1*N*R/(1+b1*R)
    dN <- -aN*N + c1*N*R/(1+b1*R) - c2*M*N/(1+b2*N)
    dM <- -aT*M + c2*M*N/(1+b2*N)
    return(list(c(dR, dN, dM)))  
  }) 
}  
p <- c(b1=6,b2=2,c1=5,c2=0.1,aN=0.4,aT=0.01,K=1)

For these parameters the model behavior approaches a limit cycle in the absence of the predator:

par(mar=c(2.6,2.6,1.6,0.2),mgp=c(1.5,0.5,0)) 
s <- c(R=1,N=0.1,M=0)
plane()
s <- run(100,0.1,traject=TRUE)

run(state=s) 

##         R         N         M 
## 0.7843797 0.3495084 0.0000000

In the presence of the predator the 3D model model displays chaotic behavior for these parameters. Run for 2000 time steps to approach that attractor

s <- c(R=1,N=0.1,M=0.2)
s <- run(2000)

and save the end result as the new state.

Make a 3-dimensional phase space using cube(), and shade the \(N'=0\) nullcline along the x (\(R\)) axis and z (\(M\)) axis. Define narrow margins of the plot by par(), and plot the location of the steady state using points3D() from the plot3D library:

par(mar=c(0,0,0,0))
cube(ymax=0.5,zmax=6,theta=30,phi=30)
cube(shade="xz",grid=9,show="N",add=TRUE)
f <- newton(s)
##         R         N         M 
## 0.9025793 0.1250000 3.7929887 
## Unstable point, eigenvalues:
## [1] -0.80809098+0.00000000i  0.02421744+0.04333401i  0.02421744-0.04333401i
points3D(x=f[1],y=f[2],z=f[3],add=TRUE,pch=1)

where the red, blue, and green nullclines correspond to \(R'=0\), \(N'=0\), \(M'=0\), respectively, and the circle denotes the location of the steady state (an unstable spiral point).

One can open a new window displaying the same graphics, and zoom and rotate with the mouse after installing the rgl and the plot3Drgl R-packages and loading the corresponding libraries. If you don’t need this, just comment out the next 3 lines.

library(rgl)
library(plot3Drgl)
plotrgl()

Redraw the 3-D space from a different viewpoint and shade both the \(R'=0\) and \(N'=0\) nullcline

par(mar=c(0,0,0,0))
cube(ymax=0.5,zmax=6,theta=30,phi=20)
cube(shade="xz",show=c("R","N"),add=TRUE)
points3D(x=f[1],y=f[2],z=f[3],add=TRUE,pch=1)

and again pop-up a window allowing you to zoom and rotate (after installing rgl and plot3Drgl).

plotrgl()

Redraw the 3-D space and add a trajectory using run3d()

par(mar=c(0,0,0,0))
cube(zmax=6,theta=20,phi=30)
run3d(tmax=5000,tstep=0.2,add=TRUE)

##         R         N         M 
## 0.5187239 0.4945453 3.3352990

where the red, blue, and green nullclines correspond to \(R'=0\), \(N'=0\), \(M'=0\), respectively, and the black line is the trajectory.

Again pop-up a window allowing you to zoom and rotate (after installing rgl and plot3Drgl).

plotrgl()

Bifurcation diagram

Now that we have a chaotic trajectory in 3 dimensions, let’s make a bifurcation diagram using continue() and by following the attractor by plotting the minimum and maximum values of trajectories for increasing values of \(b_1\).
Start with a low value of \(b_1\)

p["b1"] <- 1
s <- c(R=1,N=0.1,M=0.2)
f <- newton(run(1e4,timeplot=FALSE))
##          R          N          M 
##  0.6123724  0.1250000 18.7372436 
## Stable point, eigenvalues:
## [1] -0.07413741+0.5666056i -0.07413741-0.5666056i -0.01708235+0.0000000i

we obtain a stable steady state (stable spiral) that we continue by changing \(b_1\). The bifurcation around \(b_1=2.12\) is a Hopf bifurcation. Starting slightly above that value of \(b_1\) we first run the model for a burn-in time of a 1000 time steps. Subsequent we run it for 500 time steps, save the data, and find all minima and maxima in that data. These are all plotted as dots at the corresponding value of \(r_1\):

par(mar=c(2.6,2.6,1.6,0.2),mgp=c(1.5,0.5,0)) 
continue(f,x="b1",xmin=0,xmax=4,y="R")
## Starting at b1 = 1 with:
##          R          N          M 
##  0.6123724  0.1250000 18.7372436
## Bifurcation at b1 = 2.12
## NULL
for (b1 in seq(2.12,4,0.01)) {
  p["b1"] <- b1
  f <- run(1000,state=f,timeplot=FALSE)
  data <- run(500,tstep=0.2,state=f,table=TRUE,timeplot=FALSE)
  nr <- nrow(data)
  Rt <- data$R
  Rmin <- Rt[1+which(Rt[2:(nr-1)]<Rt[1:(nr-2)] & Rt[2:(nr-1)]<Rt[3:nr])]
  Rmax <- Rt[1+which(Rt[2:(nr-1)]>Rt[1:(nr-2)] & Rt[2:(nr-1)]>Rt[3:nr])]
  points(rep(b1,length(Rmin)),Rmin,pch=".")
  points(rep(b1,length(Rmax)),Rmax,pch=".")
}

References

Hastings & Powell (1991) Ecology 72. 896-903.

More information

Check Grind’s homepage for more tutorials