QPot (pronounced
Differential equations are an important modeling tool in virtually every scientific discipline. Most differential equation models are deterministic, meaning that they provide a set of rules for how variables change over time, and no randomness comes into play. Reality, of course, is filled with random events (i.e., noise or stochasticity). Unfortunately, many of the analytic techniques developed for deterministic ordinary differential equations are insufficient to study stochastic systems, where phenomena like noise-induced transitions between alternative stable states and metastability can occur. For systems subject to stochasticity, the quasi-potential is a tool that yields information about properties such as the expected time to escape a basin of attraction, the expected frequency of transitions between basins, and the stationary probability distribution. QPot [abbreviation of Quasi-Potential; Moore et al. (2016)] is an R package that allows users to calculate quasi-potentials, and this paper is a tutorial of its application. This package is intended for use by any researchers who are interested in understanding how stochasticity impacts differential equation models. QPot makes quasi-potential analysis accessible to a broad range of modelers, including those who have not previously encountered the topic. The key functions in package QPot are listed in Table 1.
Function | Main arguments | Description |
---|---|---|
TSTraj() |
Deterministic skeleton, |
Creates a realization (time series) of the stochastic differential equations. |
TSPlot() |
TSTraj() output |
Plots a realization of the stochastic differential equations, with an optional histogram side-plot. Plots can additionally be two-dimensional, which show realizations in |
TSDensity() |
TSTraj() output |
Creates a density plot of a trajectory in |
QPotential() |
Deterministic skeleton, stable equilibria, bounds, mesh (number of divisions along each axis) | Creates a matrix corresponding to a discretized version of the local quasi-potential function for each equilibrium. |
QPGlobal() |
Local quasi-potential matrices, unstable equilibria | Creates a global quasi-potential surface. |
QPInterp() |
Global quasi-potential, |
Evaluates the global quasi-potential at |
QPContour() |
Global quasi-potential | Creates a contour plot of the quasi-potential. |
VecDecomAll() |
Global quasi-potential, deterministic skeleton, bounds | Creates three vector fields: the deterministic skeleton, the negative gradient of the quasi-potential, and the remainder vector field. To find each field individually, the functions VecDecomVec() , VecDecomGrad() , or VecDecomRem() can be used. |
VecDecomPlot() |
Deterministic skeleton, gradient, or remainder field | Creates a vector field plot for the vector, gradient, or remainder field. |
Consider a differential equation model of the form
Consider System (2), with deterministic
skeleton (1). If there exists a function
The potential function is useful for understanding the stochastic
System (2). As in the deterministic case, the dynamics
of the stochastic system can be represented by a ball rolling on the
surface
Unfortunately, gradient systems are very special, and a generic system
of the form (1) will almost certainly not be a
gradient system. That is, there will be no function
Fortunately, quasi-potential functions generalize the concept of a
potential function for use in non-gradient systems. A non-gradient
system’s quasi-potential,
The mathematical definition of the quasi-potential is rather involved. We refer readers to Cameron (2012), Nolting and Abbott (2016), and the references therein for the technical construction.
In both this paper and in package QPot, the function that we refer to
as the quasi-potential is
QPot is an R package that contains tools for calculating and analyzing quasi-potentials (which, for the special case of gradient systems, are simply potentials). The following three examples show how to use the tools in this package. The first example is a simple consumer-resource model from ecology. This example is explained in detail, starting with the analysis of the deterministic skeleton, proceeding with simulation of the stochastic system, and finally demonstrating the calculation, analysis, and interpretation of the quasi-potential. The second and third examples are covered in less detail, but illustrate some special system behaviors. Systems with limit cycles, like Example 2, require a slightly different procedure than systems that only have point attractors. Extra care must be taken constructing global quasi-potentials for exotic systems, like Example 3.
Consider the stochastic version (sensu (2)) of a
standard consumer-resource model of plankton (
There are preexisting tools in R for analyzing the deterministic
skeleton of System (3), which will be described briefly in
this subsection. Many of the these tools can be found in the CRAN Task
View for Differential Equations
(https://CRAN.R-project.org/view=DifferentialEquations), but we use a
select few in our analysis. The first step is to find the equilibria for
the system and determine their stability with linear stability analysis.
Equilibria can be found using the package
rootSolve
(Soetaert and Herman 2008). rootSolve provides routines that allow users to
find roots of nonlinear functions, and perform equilibria and
steady-state analysis of ordinary differential equations (ODEs). In
Example 1, the equilibria are eigen()
in package base over the Jacobian matrix
(jacobian.full()
in package rootSolve), which determines the
asymptotic stability of the system. Model2String()
, which takes a function containing
equations and a list of parameters and their values, and returns the
equations in a string that is usable by QPot (see the help page for an
example).
The package phaseR
(Grayling 2014b,a) is an R package for the qualitative
analysis of one- and two-dimensional autonomous ODE systems using phase
plane methods (including the linear stability analysis described in the
preceding paragraph). We use phaseR to generate a stream plot of the
deterministic skeleton of System (3)
(Figure 1). Note that a stream plot is a phase plane
plot that displays solutions of a system of differential equations;
these solutions are also called streamlines. The deSolve package can
be used to find solutions corresponding to particular initial conditions
of the deterministic skeleton of System (3). During the
analysis of the deterministic skeleton of a system, it is important to
note several things. The first is the range of
QPot contains several tools for generating and visualizing
realizations of systems of the form (2). Examining
these realizations can help users understand qualitative features of the
system before computing and analyzing the quasi-potential.
Sim.DiffProc
(Guidoum and Boukhetala 2016) and
yuima (Brouste et al. 2014)
are two packages that offer a full suite of stochastic differential
equation simulation options, and many of the tools that they contain are
more efficient than those in QPot. Users interested in very
large-scale simulation are encouraged to seek out those packages. For
general exploration of a model’s behavior prior to quasipotential
analysis, however, TSTraj()
is extremely helpful and does not require
the use of a separate package.
Here, we show how to use QPot to obtain realizations of
System (3) for a specified level of noise intensity, TSTraj()
in QPot implements the Euler-Maruyama method.
All other code/function references hereafter are found in QPot, unless
specified otherwise. To generate a realization, the following arguments
are required: the right-hand side of the deterministic skeleton for both
equations, the initial conditions TSTraj()
accepts strings of equations with the
parameter values already included (supplied by the user or made with
Model2String()
) or can combine the equations with the parameter values
supplied as parms
. We supply the function Model2String()
to replace
parameters with their values in an equation, but the user can also input
the values themselves and may need to do so with complicated equations
(see the Model2String()
help page). Using Model2String()
will allow
a user to catch problems before they cause complications in the C code
within function QPotential()
.
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)
parms.eqn.x <- Model2String(var.eqn.x, parms = model.parms)
## Do not print to screen.
parms.eqn.y <- Model2String(var.eqn.y, parms = model.parms, supress.print = TRUE)
model.state <- c(x = 1, y = 2)
model.sigma <- 0.05
model.time <- 1000 # we used 12500 in the figures
model.deltat <- 0.025
ts.ex1 <- TSTraj(y0 = model.state, time = model.time, deltat = model.deltat,
x.rhs = parms.eqn.x, y.rhs = parms.eqn.y, sigma = model.sigma)
## 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 = var.eqn.x, y.rhs = var.eqn.y, parms = model.parms, sigma = model.sigma)
TSPlot()
, with dens = TRUE
, shows a histogram of the Figure 2 shows a realization for dim = 1
produces a time
series plot with optional histogram side-plot. The dim = 2
produces a
plot of a realization in TSDensity()
function
(e.g., Figure 3b).
[commandchars=\\\{\}]
TSPlot(ts.ex1, deltat = model.deltat) # Figure \ref{fig:ex1ts}
TSPlot(ts.ex1, deltat = model.deltat, dim = 2) # Figure \ref{fig:ex12d}a
TSDensity(ts.ex1, dim = 1) # like Figure \ref{fig:ex1ts} histogram
TSDensity(ts.ex1, dim = 2) # Figure \ref{fig:ex12d}b
Bounds can be placed on the state variables in all of the functions
described in this subsection. For example, it might be desirable to set
0 as the minimum size of a biological population, because negative
population densities are not physically meaningful. A lower bound can be
imposed on the functions described in this subsection with the argument
lower.bound
in the function TSTraj()
. Similarly, it might be
desirable to set an upper bound for realizations, and hence prevent
runaway trajectories (unbounded population densities are also not
physically meaningful). An upper bound can be imposed on the functions
described in this subsection with the argument upper.bound
.
TSPlot()
plotted in dim = 2
. (B) A density plot
obtained from a realization of System (3) using the function
TSDensity()
with dim = 2
. Red corresponds to high density, and blue
to low density.The next step is to compute a local quasi-potential for each attractor.
Because QPot deals with two-dimensional systems, attractor will be
used synonymously with stable equilibrium or stable limit cycle. A limit
cycle will be considered in example 2. For now, suppose that the only
attractors are stable equilibrium points,
In order to understand the local quasi-potential, it is useful consider
the analogy of a particle traveling according to
System (2). In the context of Example 1, the
coordinates of the particle correspond to population densities, and the
particle’s path corresponds to how those population densities change
over time. The deterministic skeleton of (2) can be
visualized as a force field influencing the particle’s trajectory.
Suppose that the particle moves along a path from a stable equilibrium
In the basin of attraction for
QPot calculates quasi-potentials using an adjustment developed by
Cameron (2012) to the ordered upwind algorithm
(Sethian and Vladimirsky 2001; Sethian and Vladimirsky 2003). The idea behind the algorithm is to
calculate
Calculating QPotential()
requires a
text string of the equations and parameter values, the stable
equilibrium points, the computation domain, and the mesh size. If the
equations do not contain the parameter values, the function
Model2String()
can be used to insert the values into the equations, as
presented above. For (3), this first means inputting the
equations:
In R:
## If not done in a previous step.
parms.eqn.x <- Model2String(var.eqn.x, parms = model.parms)
## Do not print to screen.
parms.eqn.y <- Model2String(var.eqn.y, parms = model.parms, supress.print = TRUE)
## Could also input the values by hand and use this version.
## parms.eqn.x <- "1.54 * x * (1.0 - (x / 10.14)) - (y * (x^2)) / (1.0 + (x^2))"
## parms.eqn.y <- "((0.476 * (x^2) * y) / (1 + (x^2))) - 0.112509 * (y^2)"
The coordinates of the points
eq1.x <- 1.40491
eq1.y <- 2.80808
eq2.x <- 4.9040
eq2.y <- 4.06187
Next, the boundaries of the computational domain need to be entered.
This domain will be denoted by
bounds.x <- c(-0.5, 20.0)
bounds.y <- c(-0.5, 20.0)
In some cases, it may be desirable to treat boundaries differently in the upwind algorithm. This is addressed below in Section 7.
Finally, the mesh size for the discretization of the domain needs to be
specified. Let
step.number.x <- 1000
step.number.y <- 1000 # we used 4100 in the figures
The update radii factors, k.x
and k.y
in
QPotential()
. These two parameters determine the neighborhood of
points that can be used to update a given point.
The R interface implements the QPotential()
algorithm using C code. By
default QPotential()
outputs a matrix that contains the
quasi-potentials to the R session. The time required to compute the
quasi-potential will depend on the size of the region and the fineness
of the mesh. This example with save.to.R
to TRUE
(default) outputs the matrix into the R session,
and setting the argument save.to.HD
to TRUE
saves the matrix to the
hard drive as a tab-delimited text file filename
in the current
working directory. For
eq1.local <- QPotential(x.rhs = parms.eqn.x, x.start = eq1.x, x.bound = bounds.x,
x.num.steps = step.number.x, y.rhs = parms.eqn.y, y.start = eq1.y,
y.bound = bounds.y, y.num.steps = step.number.y)
Step 3 should be repeated until local quasi-potentials
eq2.local <- QPotential(x.rhs = parms.eqn.x, x.start = eq2.x, x.bound = bounds.x,
x.num.steps = step.number.x, y.rhs = parms.eqn.y, y.start = eq2.y,
y.bound = bounds.y, y.num.steps = step.number.y)
Each local quasi-potential QPInterp()
. Inputs to
QPInterp()
include the (QPotential()
output (i.e., the
matrix with rows corresponding to QPInterp()
can be used for any of the local
quasi-potential or the global quasi-potential surfaces (see the next
subsection).
Recall that
ex1.global <- QPGlobal(local.surfaces = list(eq1.local, eq2.local),
unstable.eq.x = c(0, 4.2008), unstable.eq.y = c(0, 4.0039),
x.bound = bounds.x, y.bound = bounds.y)
This function QPGlobal()
calculates the global quasi-potential by
automatically pasting together the local quasi-potentials. This function
requires the input of all the discretized local quasi-potentials, and
the coordinates of all unstable equilibria. The output is a discretized
version of the global quasi-potential. The length of time required for
this computation will depend on the total number of mesh points; for the
parameters used in Example 1, it takes a couple of minutes. As with the
local quasi-potentials, the values of QPInterp()
.
c.parm
parameter in
QPContour()
can be used to generate non-equal contour spacing (e.g.,
for finer resolution near equilibria). The default creates evenly spaced
contour lines ((A); c.parm = 1
). In (B), contour lines are
concentrated at the bottom of the basin (c.parm = 5
). Default plot
colors are generated from package
viridis
(Garnier 2016), by setting
col.contour = viridis(n = 25, option = "D")
in
QPContour()
.To visualize the global quasi-potential, one can simply take the global
quasi-potential matrix from QPGlobal()
and use it to create a contour
plot using QPContour()
(Figure 4).
[commandchars=\\\{\}]
QPContour(surface = ex1.global, dens = c(1000, 1000), x.bound = bounds.x,
y.bound = bounds.y, c.parm = 5) # right side of Figure \ref{fig:ex1qp}
QPContour()
is based on the .filled.contour()
function from the base
package graphics. In most cases, the mesh sizes used for the
quasi-potential calculation will be much finer than what is required for
useful visualization. The argument dens
within QPContour()
reduces
the points used in the graphics generation. Although it might seem
wasteful to perform the original calculations at a mesh size that is
finer than the final visualization, this is not so. Choosing the mesh
size in the original calculations to be very fine reduces the
propagation of errors in the ordered upwind algorithm, and hence leads
to a more accurate numerical solution.
An additional option allows users to specify contour levels. R’s default
for the contour()
function creates contour lines that are equally
spaced over the range of values specified by the user. In some cases,
however, it is desirable to use a non-equidistant spacing for the
contours. For example, equally-spaced contours will not capture the
topography at the bottom of a basin if the changes in height are much
smaller than in other regions in the plot. Simply increasing the number
of equally-spaced contour lines does not solve this problem, because
steep areas of the plot become completely saturated with lines.
QPContour()
has a function for non-equidistant contour spacing that
condenses contour lines at the bottoms of basins. Specifically, for c.parm = 5
).
Finally, creating a 3D plot can be very useful for visualizing the
features of more complex surfaces. This is especially helpful when
considering the physical metaphor of a ball rolling on a surface
specified by a quasi-potential (Nolting and Abbott 2016). R has several packages
for 3D plotting, including static plotting with the base function
persp()
and with the package
plot3D
(Soetaert 2014). Interactive plotting is provided by
rgl (Adler et al. 2015). To
create an interactive 3D plot for Example 1 using rgl, use the code:
persp3d(x = ex1.global)
. Figure 5 shows a 3D plot of
example 1 using persp3D(z = ex1.global)
in plot3D that clearly
illustrates the differences between the two local basins. Users can also
export the matrix of quasi-potential values and create 3D plots in other
programs.
persp3D()
in package plot3D. 3D plotting
can further help users visualize the quasi-potential surfaces. Plot
colors are generated from package viridis by setting
col = viridis(n = 100, option = "A")
and
contour = TRUE
.Recall that the deterministic skeleton (1) can be
visualized as a vector field, as shown in Figure 1. In
gradient systems, this vector field is completely determined by the
potential function,
The remainder vector field can be interpreted as a force that causes
trajectories to circulate around level sets of the quasi-potential.
QPot enables users to perform this decomposition. The function
VecDecomAll()
calculates the vector field decomposition, and outputs
three vector fields: the original deterministic skeleton,
VecDecomVec()
,
VecDecomGrad()
, or VecDecomRem()
. These vector fields can be
visualized using the function VecDecomPlot()
. Code to create the
vector fields from VecDecomAll()
is displayed below; code for
generating individual vector fields can be found in the man pages
accessible by help()
for VecDecomVec()
, VecDecomGrad()
, or
VecDecomRem()
. The gradient and remainder vector fields are shown in
the left and right columns of Figure 6, respectively, with
proportional vectors (top row) and equal-length vectors (bottom row).
Three arguments within VecDecomPlot()
are important to creating
comprehensible plots: dens
, tail.length
, and head.length
. The
dens
parameter specifies the number of arrows in the plot window along
the tail.length
scales the length of
arrow tails. The argument head.length
scales the length of arrow
heads. The function arrows()
makes up the base of VecDecomPlot()
,
and arguments can be passed to it, as well as to plot()
. The code
below produces all three vector fields from the multi-dimensional array
returned by VecDecomAll()
:
## Calculate all three vector fields.
VDAll <- VecDecomAll(surface = ex1.global, x.rhs = parms.eqn.x, y.rhs = parms.eqn.y,
x.bound = bounds.x, y.bound = bounds.y)
## Plot the deterministic skeleton vector field.
VecDecomPlot(x.field = VDAll[, , 1], y.field = VDAll[, , 2], dens = c(25, 25),
x.bound = bounds.x, y.bound = bounds.y, xlim = c(0, 11), ylim = c(0, 6),
arrow.type = "equal", tail.length = 0.25, head.length = 0.025)
## Plot the gradient vector field.
VecDecomPlot(x.field = VDAll[, , 3], y.field = VDAll[, , 4], dens = c(25, 25),
x.bound = bounds.x, y.bound = bounds.y, arrow.type = "proportional",
tail.length = 0.25, head.length = 0.025)
## Plot the remainder vector field.
VecDecomPlot(x.field = VDAll[, , 5], y.field = VDAll[, , 6], dens = c(25, 25),
x.bound = bounds.x, y.bound = bounds.y, arrow.type = "proportional",
tail.length = 0.35, head.length = 0.025)
arrow.type = "proportional"
(top row) and
arrow.type = "equal"
(bottom row) arrow lengths using VecDecomPlot()
for System (3).Consider the following model:
The deterministic skeleton of this system has one equilibrium,
Figures 8 and 9a show a time series for a
realization of (4) with
TSPlot()
, with dens = TRUE
, shows a histogram of the [commandchars=\\\{\}]
var.eqn.x <- "- (y - beta) + mu * (x - alpha) * (1 - (x - alpha)^2 - (y - beta)^2)"
var.eqn.y <- "(x - alpha) + mu * (y - beta) * (1 - (x - alpha)^2 - (y - beta)^2)"
model.state <- c(x = 3, y = 3)
model.parms <- c(alpha = 4, beta = 5, mu = 0.2)
model.sigma <- 0.1
model.time <- 1000 # we used 2500 in the figures
model.deltat <- 0.005
ts.ex2 <- TSTraj(y0 = model.state, time = model.time, deltat = model.deltat,
x.rhs = var.eqn.x, y.rhs = var.eqn.y, parms = model.parms, sigma = model.sigma)
TSPlot(ts.ex2, deltat = model.deltat) # Figure \ref{fig:ex2ts}
TSPlot(ts.ex2, deltat = model.deltat, dim = 2, line.alpha = 25) # Figure \ref{fig:ex22d}a
TSDensity(ts.ex2, dim = 1) # Histogram
TSDensity(ts.ex2, dim = 2) # Figure \ref{fig:ex22d}b
dim = 2
in the function TSPlot()
) (B) A density plot
obtained from a realization of System (4) using TSDensity()
with dim = 2
. Red corresponds to high density, and blue to low
density.In this example, there are no stable equilibrium points. There is one
stable limit cycle, and this can be used to obtain a local
quasi-potential. Using
eqn.x <- Model2String(var.eqn.x, parms = model.parms)
eqn.y <- Model2String(var.eqn.y, parms = model.parms)
eq1.qp <- QPotential(x.rhs = eqn.x, x.start = 4.15611, x.bound = c(-0.5, 7.5),
x.num.steps = 4000, y.rhs = eqn.y, y.start = 5.98774, y.bound = c(-0.5, 7.5),
y.num.steps = 4000)
QPContour()
. Yellow corresponds to low
values of the quasi-potential, and purple to high
values.There is only one local quasi-potential in this example, so it is the
global quasi-potential,
Figure 10 shows a contour plot of the global quasi-potential.
QPContour(eq1.qp, dens = c(1000, 1000), x.bound = c(-0.5, 7.5),
y.bound = c(-0.5, 7.5), c.parm = 10)
In Example 1, the procedure for pasting local quasi-potentials together
into a global quasi-potential was a simple, two-step process. First, one
of the local quasi-potentials was translated so that the two surfaces
agreed at the saddle point separating the two basins of attraction.
Second, the global quasi-potential was obtained by taking the minimum of
the two surfaces at each point. A general algorithm for pasting local
quasi-potentials, as explained in Graham and Tél (1986) and Roy and Nauman (1995), is
slightly more complicated. This process is automated in QPGlobal()
,
but it is worth understanding the process in order to correctly
interpret the outputs.
To understand the full algorithm, consider the following model:
The deterministic skeleton of this system has five equilibria. These are
Figures 12 and 13a show a time series for a
realization of System (5) with
[commandchars=\\\{\}]
var.eqn.x <- "x * ((1 + alpha1) - (x^2) - x * y - (y^2))"
var.eqn.y <- "y * ((1 + alpha2) - (x^2) - x * y - (y^2))"
model.state <- c(x = 0.5, y = 0.5)
model.parms <- c(alpha1 = 1.25, alpha2 = 2)
model.sigma <- 0.8
model.time <- 5000
model.deltat <- 0.01
ts.ex3 <- TSTraj(y0 = model.state, time = model.time, deltat = model.deltat,
x.rhs = var.eqn.x, y.rhs = var.eqn.y, parms = model.parms, sigma = model.sigma)
TSPlot(ts.ex3, deltat = model.deltat) # Figure \ref{fig:ex3ts}
TSPlot(ts.ex3, deltat = model.deltat, dim = 2 , line.alpha = 25) # Figure \ref{fig:ex32d}a
TSDensity(ts.ex3, dim = 1) # Histogram of time series
TSDensity(ts.ex3, dim = 2 , contour.levels = 20 , contour.lwd = 0.1) # Figure \ref{fig:ex32d}b
TSPlot()
, with dens = TRUE
in the function TSPlot()
, shows a histogram of the TSPlot()
with dim = 2
. (B) A density plot
obtained from the realization of System (5) by using the
function TSDensity()
with dim = 2
, contour.levels = 20
, and
contour.lwd = 0.1
. Red corresponds to high density, and blue to low
density.Two local quasi-potentials need to be calculated,
equation.x <- Model2String(var.eqn.x, parms = model.parms)
equation.y <- Model2String(var.eqn.y, parms = model.parms)
bounds.x <- c(-3, 3); bounds.y <- c(-3, 3)
step.number.x <- 6000; step.number.y <- 6000
eq1.x <- 0; eq1.y <- -1.73205
eq2.x <- 0; eq2.y <- 1.73205
eq1.local <- QPotential(x.rhs = equation.x, x.start = eq1.x, x.bound = bounds.x,
x.num.steps = step.number.x, y.rhs = equation.y, y.start = eq1.y,
y.bound = bounds.y, y.num.steps = step.number.y)
eq2.local <- QPotential(x.rhs = equation.x, x.start = eq2.x, x.bound = bounds.x,
x.num.steps = step.number.x, y.rhs = equation.y, y.start = eq2.y,
y.bound = bounds.y, y.num.steps = step.number.y)
If one were to naively try to match the local quasi-potentials at
ex3.global <- QPGlobal(local.surfaces = list(eq1.local, eq2.local),
unstable.eq.x = c(0, -1.5, 1.5), unstable.eq.y = c(0, 0, 0), x.bound = bounds.x,
y.bound = bounds.y)
Figure 14 shows a contour plot of the global quasi-potential. Note that the surface is continuous, but not smooth. The lack of smoothness is a generic feature of global quasi-potentials created from pasting local quasi-potentials. Cusps usually form when switching from the part of solution obtained from one local quasi-potential to the other.
QPContour(ex3.global, dens = c(1000, 1000), x.bound = bounds.x, y.bound = bounds.y,
c.parm = 5)
QPContour()
. Yellow
corresponds to low values of the quasi-potential, and purple to high
values.It is important to consider the type of behavior that should be enforced
at the boundaries and on coordinate axes (QPotential()
, if bounce = "d"
corresponds to
this (d)efault behavior. A user can prevent the ordered upwind method
from passing trajectories through negative phase space by using the
option bounce = "p"
for (p)ositive values only. This option can be
interpreted as a reflecting boundary condition. It forces the front of
solutions obtained by the ordered upwind method to stay in the defined
boundaries, which is positive phase space in this case. A more generic
option is bounce = "b"
for (b)ounce, which allows users to supply
reflecting boundaries other than the coordinate axes. These are set with
x.bound
and y.bound
. Small numerical errors at a reflecting boundary
can cause the algorithm to terminate prematurely. To avoid this, the
option bounce.edge
adds a small amount of padding between the
reflecting boundary and the edge of the computational domain.
In the cases considered so far, the noise terms for the
Many models contain multiplicative noise terms. These are of the form:
QPot is an R package that provides several important tools for analyzing two-dimensional systems of stochastic differential equations. Future efforts will work toward extending QPot to higher-dimensional systems, but this is a computationally challenging task. QPot includes functions for generating realizations of the stochastic differential equations, and for analyzing and visualizing the results. A central component of QPot is the calculation of quasi-potential functions, which are highly useful for studying stochastic dynamics. For example, quasi-potential functions can be used to compare the stability of different attractors in stochastic systems, a task that traditional linear stability analysis is poorly suited for (Nolting and Abbott 2016). By offering an intuitive way to quantify attractor stability, quasi-potentials are poised to become an important means of understanding phenomena like metastability and alternative stable states. QPot makes quasi-potentials accessible to R users interested in applying this new framework.
K.C.A, C.M.M., B.C.N., and C.R.S. designed the project. M.K.C. wrote the C code for finding the quasi-potential; C.M.M. and C.R.S. wrote the R code and adapted the C code.
This work was supported by a Complex Systems Scholar grant to K.C.A. from the James S. McDonnell Foundation. M.K.C. was partially supported by NSF grant 1217118. We also thank the anonymous reviewers for their constructive comments and suggestions which helped us to improve the quality of our paper.
QPot, rootSolve, deSolve, phaseR, Sim.DiffProc, yuima, viridis, plot3D, rgl
DifferentialEquations, Epidemiology, Finance, HighPerformanceComputing, NumericalMathematics, Psychometrics, Spatial, SpatioTemporal, TimeSeries
This article is converted from a Legacy LaTeX article using the texor package. The pdf version is the official version. To report a problem with the html, refer to CONTRIBUTE on the R Journal homepage.
Text and figures are licensed under Creative Commons Attribution CC BY 4.0. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".
For attribution, please cite this work as
Moore, et al., "QPot: An R Package for Stochastic Differential Equation Quasi-Potential Analysis", The R Journal, 2016
BibTeX citation
@article{RJ-2016-031, author = {Moore, Christopher M. and Stieha, Christopher R. and Nolting, Ben C. and Cameron, Maria K. and Abbott, Karen C.}, title = {QPot: An R Package for Stochastic Differential Equation Quasi-Potential Analysis}, journal = {The R Journal}, year = {2016}, note = {https://rjournal.github.io/}, volume = {8}, issue = {2}, issn = {2073-4859}, pages = {19-38} }