In Moore et al. (2016) we use an example of a coupled consumer-resource equation developed for plankton, $x$, and their consumers, $y$:
$$\frac{dx(t)}{dt} = \alpha x(t)\left(1 - \frac{x(t)}{\beta}\right) - \frac{\delta x^2(t)y(t)}{\kappa + x^2(t)}$$ $$\frac{dy(t)}{dt} = \frac{\gamma x^2(t)y(t)}{\kappa + x^2(t)} - \mu y^2(t)$$
For this system, we use a specific parametrization to generate two non-trivial stable equilibria. Specifically, we use the parameters:
Parameter | Value | Biological description ---------- | :-------- | :---------------------------------------- $\alpha$ | 1.54 | Maximal growth rate of plankton $\beta$ | 10.14 | Carrying capacity of plankton $\delta$ | 1.0 | Maximal feeding rate of the consumers $\gamma$ | 0.476 | Conversion rate of plankton to consumer $\mu$ | 0.112509 | Death rate of the consumer
We can create a vector field of the deterministic skeleton by using package phaseR. First, we will load phaseR and its dependency, deSolve (phaseR makes use of deSolve::ode
):
library(package = "deSolve") library(package = "phaseR")
Second, we write our equations above in the pseudo-code format:
model <- function(time, initial conditions, parameters){ assign state variables to initial conditions assign parameters create an object to store output equations a list as the output }
in R as
model.ex1 <- function(t, y, parameters) { x <- y[1] y <- y[2] alpha <- parameters["alpha"] beta <- parameters["beta"] delta <- parameters["delta"] kappa <- parameters["kappa"] gamma <- parameters["gamma"] mu <- parameters["mu"] dy <- numeric(2) dy[1] <- (alpha*x)*(1-(x/beta)) - ((delta*(x^2)*y)/(kappa + (x^2))) dy[2] <- ((gamma*(x^2)*y)/(kappa + (x^2))) - mu*(y^2) list(dy) }
Then, we plot the direction field and the zero-growth isoclines (i.e., nullclines):
xlims <- c(0, 6) ylims <- c(0, 6) flowField(deriv = model.ex1, x.lim = xlims, y.lim = ylims, parameters = model.parms, points = 30, add = F) nullclines(model.ex1, x.lim = xlims, y.lim = ylims, parameters = model.parms, points = 250, col = c("blue","red"))
We can see from the field that trajectories may take many paths to different areas or points in phase space (i.e., there are multiple basins of attraction). But more reliably, we can see that the nullclines cross several times, which means that at those points we have equilibria. For the purposes of our example, we are interested in interior points (i.e., $x$ and $y$ have populations > 0).
There are several ways to find solutions and classify equilibrium points. Users are encouraged to familiarize themselves with CRAN Task View: Differential Equations for a summary of available packages. Here, we use package rootSolve to find equilibria. At each equilibrium point, we also classify the behavior about the point by finding the eigenvalues of the Jacobian matrix, known as linear or local stability analysis.
For simple models, equilibria can be found analytically. But for most non-linear models, a solver must be used. Because we have an idea where the equilibria are, we create an area to find steady-state equilibria using rootSolve::stode
. Because it's a small space, we'll first use a for
loop over the x- and y-area we wish to sample, then we find unique values of (x, y) to give us a matrix of unique equilibria.
library(package = "rootSolve") xspace <- seq(from = 1, to = 5, length.out = 10) yspace <- seq(from = 2.5, to = 4, length.out = 10) l.xspace <- length(x = xspace) l.yspace <- length(x = yspace) space.mat <- matrix(data = NA, nrow = l.xspace*l.yspace, ncol = 2) for (i in 1:l.xspace){ for (j in 1:l.yspace){ y <- c(x = xspace[i], y = yspace[j]) STO <- stode(y = y, func = model.ex1, parms = model.parms, positive = T) space.mat[(((i-1)*l.xspace)+j),] <- STO$y } } eqs <- unique(x = round(x = space.mat, digits = 3))
And, for further confirmation, if the plot window is still open, we can add them to ensure they look correct:
points(x = eqs[,1], y = eqs[,2], cex = 1.5)
Once we have equilibria, we can classify them using phaseR::stability
. There's a good deal of information beyond the $classification
that we call and with summary = T
.
for (i in 1:nrow(eqs)){ print(x = paste0("x = ", eqs[i, 1], ", y = ", eqs[i, 2], " is a ", stability(deriv = model.ex1, y.star = eqs[i,], parameters = model.parms, summary = F)$classification)) }
This will reveal our three interior equilibria---two stable and one unstable equilibrium. Our current focus is on determining the relative stability around the two interior equilibria when stochasticity is added to the system. In the next vignette, we heuristically run some stochastic simulations to better understand how the system behaves with stochasticity.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.