This document was made in the RStudio environment, and expects that you have installed the deSolve, rootSolve, coda and FME R-packages beforehand (which 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-05-2024) was sourced

Grind needs to sourced to define its functions. This document uses run(), plane(), newton() and continue().

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=K=a=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, 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)

where the arguments (t, state, parms) provide the time, the state \((R,N)\), and the parameters (as defined by the vector p) to the equations. The with(as.list(c(state,parms)),..) function allows one to identify parameters and variables by their name (instead of having to write parms["K"] and state["R"], and so on). Note that the derivatives are returned as a list, and that the order of the variables in this list has to be the same as their order in the state!

Numerical integration and phase plane analysis

One can simulate the model for 50 time steps by calling run(50)

run(50)

##         R         N 
## 0.5000029 0.4999935

One can sketch nullclines by calling plane(),

plane()

which uses the first variable R for the horizontak axis, and the second N for the vertical axis. One ask for a vector field in plane() by setting the option vector=TRUE, and one can add a trajectory by using the traject option in run()

plane(vector=TRUE)
run(traject=TRUE)

##   R   N 
## 0.5 0.5

or make a phase portrait by adding the portrait option to plane()

plane(portrait=TRUE,tstep=0.1)

Finding and continuing steady states

The newton() function searches a steady state in the neighborhood of the state provided as its first argument, and provides the eigenvalues of the Jacobian of that state. Such a state can be continued along one of the parameters with the function continue()

s1 <- newton(c(R=0.5,N=0.5))
##   R   N 
## 0.5 0.5 
## Stable spiral point, eigenvalues:
## [1] -0.25+0.4330127i -0.25-0.4330127i
s2 <- newton(c(R=1,N=0))
## R N 
## 1 0 
## Saddle point (unstable), eigenvalues:
## [1] -1.0  0.5
continue(s1,x="K",xmin=0.1,xmax=2,y="N",ymin=-0.01,step=0.001)
## Starting at K = 1 with:
##   R   N 
## 0.5 0.5
## Bifurcation at K = 0.5
## NULL
continue(s2,step=0.001,add=TRUE)

## Starting at K = 1 with:
## R N 
## 1 0 
## Bifurcation at K = 0.498
## NULL

To see what happens at this bifurcation, we set the \(K\) parameter to its critical value \(K=0.5\) and sketch the corresponding phase plane

p["K"] <- 0.5
plane()

More information

These examples provide the basics of running models, drawing nullclines, and finding steady states. More information can be found on the webpages on running and fitting. Check Grind’s homepage for even more info.