In the previous vignette we used an example of a coupled consumer-resource equation developed for plankton and their consumers. We now write the equation as a stochastic differential equation:

$$\mathrm{d}X(t) = \left(\alpha X(t)\left(1 - \frac{X(t)}{\beta}\right) - \frac{\delta X^2(t)Y(t)}{\kappa + X^2(t)}\right)dt + \sigma \mathrm{d}W_1(t)$$ $$\mathrm{d}Y(t) = \left(\frac{\gamma X^2(t)Y(t)}{\kappa + X^2(t)} - \mu Y^2(t)\right) + \sigma \mathrm{d}W_2(t)$$

There are several existing packages in R that will run stochastic simulations (see references in CRAN Task View: Differential Equations), but we include a less-computationally-efficient function for less-intensive simulations that are sufficient for our purposes and takes arguments in the form that we later use for the quasi-potential analysis.

Specifically, we created `TSTraj`

(for **T**ime **S**eries **Traj**ectory), to take a equations as strings, with the option to specify the parameter values as a separate list. `TSTraj`

allows users to add stochasticity to the model (`sigma`

), control the time step ($\Delta t$), set upper and lower bounds (e.g., a lower bound of 0 may be useful for biologists studying populations since a population of < 0 is yet to be discovered), and returns a matrix with the numbers of rows equal to the total time steps ($\Delta t \times T$) and two columns for each state variable (e.g., $X$ and $Y$).

Specifying the above equations can be separate from the parameters:

var.eqn.x <- "(alpha * x) * (1 - (x / beta)) - ((delta * (x^2) * y) / (kappa + (x^2)))" var.eqn.y <- "((gamma * (x^2) * y) / (kappa + (x^2))) - mu * (y^2)" model.parms <- c(alpha = 1.54, beta = 10.14, delta = 1, gamma = 0.476, kappa = 1, mu = 0.112509)

or together, using `QPot::Model2String`

parms.eqn.x <- Model2String(model = var.eqn.x, parms = model.parms, supress.print = T) parms.eqn.y <- Model2String(model = var.eqn.y, parms = model.parms, supress.print = T)

For this simulation, we start our simulation at (1, 2), add Gaussian noise with a mean of 0 and a standard deviation of 0.05 every $\Delta t = 0.025$, for a total time of $T = 1000$.

model.state <- c(x = 1, y = 2) model.sigma <- 0.05 model.time <- 1000 model.deltat <- 0.025 set.seed(6174) ts.ex1 <- TSTraj(y0 = model.state, time = model.time, deltat = model.deltat, x.rhs = var.eqn.x, y.rhs = var.eqn.y, sigma = model.sigma, parms = model.parms)

Or alternatively, one could also use TSTraj to combine equation strings and parameter values.

ts.ex1 <- TSTraj(y0 = model.state, time = model.time, deltat = model.deltat, x.rhs = parms.eqn.x, y.rhs = parms.eqn.y, parms = model.parms, sigma = model.sigma)

We have two functions that allow users to see the results from `TSTraj`

: `TSPlot`

and `TSDensity`

. First, `TSPlot`

plots the time series for each state variable if `dim = 1`

(default) and shows the trajectory in state space if `dim = 2`

. For `dim = 1`

, we provide an option (default) to plot the density of each state variable adjacent to the time series plot.

TSPlot(mat = ts.ex1, deltat = model.deltat) TSPlot(mat = ts.ex1, deltat = model.deltat, dim = 2)

Second, `TSDensity`

takes the simulation results and plots it as either a single (`dim = 1`

) dimension or in two dimensional state space (`dim = 2`

).

TSDensity(mat = ts.ex1, dim = 1) TSDensity(mat = ts.ex1, dim = 2)

We can see from these simulations that the system spends a great deal of time around the stable focus at x = 1.405 and y = 2.808 than the stable node at x = 4.904 and y = 4.062. This realization is typical of this system, regardless of $\Delta t$, $T$, or $\sigma$. We should therefore describe the behavior of system in a way that captures this behavior.

The typical way of describing stability is through linear stability analysis, which respectively yields the following results for stable focus and stable node:

stability(deriv = model.ex1, y.star = eqs[1,], parameters = model.parms, summary = F)$eigenvalues stability(deriv = model.ex1, y.star = eqs[3,], parameters = model.parms, summary = F)$eigenvalues

The largest real eigenvalue of the Jacobian matrix for the stable focus is -0.0473848 and for the stable node is -0.37737660. Because the stable node is larger (i.e., more negative), **we would conclude that the stable node is more stable than the stable focus---in direct contrast to what our simulation shows**. This type of stability---known as asymptotic stability---is not necessarily sufficient for describing stability in stochastic systems, with continual perturbations.

This is why we need another tool that will better describe the behavior of how a system will behave with continual perturbations. For this, Nolting and Abbott (2016) have argued that the quasi-potential should be used.

For a more in-depth description and mathematical details, we again encourage readers to see

B. C. Nolting and K. C. Abbott. Balls, cups, and quasi-potentials: Quantifying stability in stochastic systems. Ecology, 97(4):850–864, 2016.

and specifically **$\S$ A Path Through the Quagmire of Stability Concepts**.

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.