QPot: An R Package for Stochastic Differential Equation Quasi-Potential Analysis
In QPot: Quasi-Potential Analysis for Stochastic Differential Equations

Analyzing the deterministic skeleton

Example 1 from Moore et al. (2016)

A system of equations

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:

Visualizing dynamics

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"))

Identifying and classifying equilibria

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.

Finding equilibria

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)

Classifying equilibria

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.