In bmarkslash7/QPot: Quasi-Potential Analysis for Stochastic Differential Equations

Stochastic simulation

Example 1 from Moore et al. (2016)

Creating a model with stochasticity

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.

Running the simulation

Specifically, we created TSTraj (for Time Series Trajectory), 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)

Visualizing the results

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)

Discrepancies in the notion of stability in stochastic systems

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