In PabRod/rolldown: Computation of approximate potentials for weakly non-gradient fields
knitr::opts_chunk$set(echo = TRUE)
knitr::include_graphics('img/logo.png')
Installation
To install the package, please type
devtools::install_github("PabRod/rolldown", ref = "master")
in your R console.
Libraries
In this vignette we will use the following libraries:
library(rolldown)
# To calculate some trajectories
library(deSolve)
# To plot our results
library(ggplot2)
# To arrange our plots in panels
library(latticeExtra)
library(gridExtra)
# For nicer plots
library(colorRamps)
One dimensional examples
Allee effect
A single-species population dynamics model with Allee effect is governed by the following differential equation:
$$
\frac{dN}{dt} = r N \left( \frac{N}{A} - 1 \right) \left( 1 - \frac{N}{K} \right)
$$
It is easy to see that this differential equation has three equilibrium points, $N = 0$, $N = K$ and $N = A$, being all of them stable but the latter one, which is unstable. We'll use the parameters $r = 1$, $A = 0.5$ and $K = 1$.
r <- 1
A <- 0.5
K <- 1
f <- function(x) { r * x * (x/A - 1) * (1 - x/K) }
We can use our method approxPot1D to approximate the potential function at a set of points. First, we have to create the points.
xs <- seq(0, 1.25, by = 0.01)
and then pass them to our algorithm:
Vs <- approxPot1D(f, xs)
By plotting the result, we clearly see that the two stable equilibria appear at $N = 0$ and $N = K = 1$, and the unstable one at $N = A = 0.5$, as we expected.
plot(xs, Vs,
type = 'l', xlab = 'N', ylab = 'V')
Two dimensional examples
In this section we'll apply our method to a collection two dimensional systems.
Synthetic examples
We generated some abstract, synthetic examples in order to test our method. Here we present some of them.
A gradient system: the four well potential
In this section we'll deal with the two-dimensional differential equation given by:
$$
\begin{cases}
\frac{dx}{dt} = f(x,y) = -x(x^2-1) \
\frac{dy}{dt} = g(x,y) = -y(y^2-1)
\end{cases}
$$
This is a gradient system (because $\frac{\partial f}{\partial x} = \frac{\partial g}{\partial y}$ everywhere). This means that the gradient - curl decomposition will have zero curl term everywhere and, thus, there exists a well defined potential. Particularly, the potential can be analytically proven to be:
$$
V(x,y) = \frac{x^2}{4}(x^2 - 2) + \frac{y^2}{4}(y^2 - 2) + V_0
$$
Let's try to compute it using our algorithm. First, we'll code our function as a vector:
f <- function(x) {c(-x[1]*(x[1]^2 - 1),
-x[2]*(x[2]^2 - 1))}
Our region of interest is now two-dimensional. We need, thus, two vectors to create our grid of points:
xs <- seq(-1.5, 1.5, by = 0.025)
ys <- seq(-1.5, 1.5, by = 0.025)
Now we are ready to apply approxPot2D:
result <- approxPot2D(f, xs, ys)
# Transform result into dataframe
data <- expand.grid(X = xs, Y = ys)
data$V <- as.vector(result$V)
data$err <- as.vector(result$err)
# Input equilibrium points (calculated externally)
eqPoints <- data.frame(x_eq = c(-1, -1, 0, 1, 1),
y_eq = c(-1, 1, 0, -1, 1),
equilibrium = factor(c('stable', 'stable', 'unstable', 'stable', 'stable')))
result is a list that contains two fields:
result$V contains the estimated values of the potentials at each grid pointresult$err contains the estimated error at each grid point
By plotting them we see:
nbins <- 15
plotV <- ggplot() +
geom_tile(data = data, aes(x = X, y = Y, fill = V)) +
geom_contour(data = data, aes(x = X, y = Y, z = V), colour = 'white', alpha = 0.5, bins = nbins) +
geom_point(data = eqPoints, aes(x = x_eq, y = y_eq, color = equilibrium)) +
coord_fixed() +
scale_fill_gradientn(colours = colorRamps::matlab.like(nbins)) +
xlab("x") + ylab("y") + ggtitle("Approximate potential") +
theme_bw()
plotErr <- ggplot() +
geom_tile(data = data, aes(x = X, y = Y, fill = err)) +
# geom_contour(data = data, aes(x = X, y = Y, z = err), colour = 'white', alpha = 0.5, bins = nbins) +
geom_point(data = eqPoints, aes(x = x_eq, y = y_eq, color = equilibrium)) +
coord_fixed() +
scale_fill_gradientn(colours = colorRamps::green2red(nbins), limits = c(0,1)) +
xlab("x") + ylab("y") + ggtitle("Error map") +
theme_bw()
grid.arrange(plotV, plotErr, ncol = 2)
Provided our example is a purely gradient system, it should not surprise us that the error is zero everywhere.
max(result$err) == 0
A non-gradient system
In this example we will apply our algorithm to the system given below:
$$
\begin{cases}
\frac{dx}{dt} = f(x,y) = -y \
\frac{dy}{dt} = g(x,y) = x
\end{cases}
$$
This is an extreme case. The gradient - curl decomposition will give us zero gradient part everywhere. Let's feed our algorithm with these example to see what happens:
First, we code the dynamics in vector form:
f <- function(x) {c(-x[2],
x[1])}
Secondly, we define our region of interest:
xs <- seq(-2, 2, by = 0.05)
ys <- seq(-2, 2, by = 0.05)
And then we are ready to apply our algorithm:
result <- approxPot2D(f, xs, ys)
# Transform result into dataframe
data <- expand.grid(X = xs, Y = ys)
data$V <- as.vector(result$V)
data$err <- as.vector(result$err)
# Input equilibrium points (calculated externally)
eqPoints <- data.frame(x_eq = c(0),
y_eq = c(0),
equilibrium = factor(c('unstable')))
The resulting approximate potential is plotted below. Being this a purely non-gradient system we expect our pseudopotential to not be trustworthy. By calculating the error we can see that actually that's the case.
nbins <- 15
plotV <- ggplot() +
geom_tile(data = data, aes(x = X, y = Y, fill = V)) +
geom_contour(data = data, aes(x = X, y = Y, z = V), colour = 'white', alpha = 0.5, bins = nbins) +
geom_point(data = eqPoints, aes(x = x_eq, y = y_eq, color = equilibrium)) +
coord_fixed() +
scale_fill_gradientn(colours = colorRamps::matlab.like(nbins)) +
xlab("x") + ylab("y") + ggtitle("Approximate potential") +
theme_bw()
plotErr <- ggplot() +
geom_tile(data = data, aes(x = X, y = Y, fill = err)) +
# geom_contour(data = data, aes(x = X, y = Y, z = err), colour = 'white', alpha = 0.5, bins = nbins) +
geom_point(data = eqPoints, aes(x = x_eq, y = y_eq, color = equilibrium)) +
coord_fixed() +
scale_fill_gradientn(colours = colorRamps::green2red(nbins), limits = c(0,1)) +
xlab("x") + ylab("y") + ggtitle("Error map") +
theme_bw()
grid.arrange(plotV, plotErr, ncol = 2)
The fact that the underlying equations are non-gradient has been captured by the algorithm.
Biological examples
Here we apply our methods to some dynamical equations well known in biology. While the abstract equations in the previous sections can be manipulated to increase or decrease their curl to gradient ratio, equations describing natural dynamical systems don't allow such a manipulation. Once again, the error map will let us know if our system allows a pseudopotential or not.
Simple regulatory gene network
A bistable network cell fate model can be described by the set of equations:
$$
\begin{cases}
\frac{dx}{dt} = f(x,y) = b_x - r_x x + \frac{a_x}{k_x + y^n} \
\frac{dy}{dt} = g(x,y) = b_y - r_y y + \frac{a_y}{k_y + x^n}
\end{cases}
$$
Such a system represents two genes ($x$ and $y$) that inhibit each other. This circuit works as a toggle switch with two stable steady states, one with dominant $x$ , the other with dominant $y$.
We can code it in vector form:
# Parameters
bx <- 0.2
ax <- 0.125
kx <- 0.0625
rx <- 1
by <- 0.05
ay <- 0.1094
ky <- 0.0625
ry <- 1
n <- 4
# Dynamics
f <- function(x) {c(bx - rx*x[1] + ax/(kx + x[2]^n),
by - ry*x[2] + ay/(ky + x[1]^n))}
This set of equations is, in general, not gradient (because $\frac{\partial f}{\partial y} \neq \frac{\partial g}{\partial x}$). Anyways, we can use the method approxPot2D to compute the approximate potential.
First, we need to define our region of interest:
xs <- seq(0, 4, by = 0.05)
ys <- seq(0, 4, by = 0.05)
And then we are ready to apply our algorithm:
result <- approxPot2D(f, xs, ys)
# Transform result into dataframe
data <- expand.grid(X = xs, Y = ys)
data$V <- as.vector(result$V)
data$err <- as.vector(result$err)
# Input equilibrium points (calculated externally)
#
# Estimated with Wolfram Alpha
# Prompt: 0 = 0.2 - x + 0.125/(0.0625 + y^4); 0 = 0.05 - y + 0.1094/(0.0625 + x^4)
eqPoints <- data.frame(x_eq = c(0.213416, 0.559865, 2.19971),
y_eq = c(1.74417, 0.730558, 0.0546602),
equilibrium = factor(c('stable', 'unstable', 'stable')))
The resulting approximate potential is plotted below. Being this not a gradient system it is advisable to plot also the estimated error. The areas in green represent small approximation error, so the potential can be safely used in those regions.
nbins <- 25
plotV <- ggplot() +
geom_tile(data = data, aes(x = X, y = Y, fill = V)) +
geom_contour(data = data, aes(x = X, y = Y, z = V), colour = 'white', alpha = 0.5, bins = nbins) +
geom_point(data = eqPoints, aes(x = x_eq, y = y_eq, color = equilibrium)) +
coord_fixed() +
scale_fill_gradientn(colours = colorRamps::matlab.like(nbins)) +
xlab("x") + ylab("y") + ggtitle("Approximate potential") +
theme_bw()
plotErr <- ggplot() +
geom_tile(data = data, aes(x = X, y = Y, fill = err)) +
geom_contour(data = data, aes(x = X, y = Y, z = err), colour = 'white', alpha = 0.5, bins = nbins) +
geom_point(data = eqPoints, aes(x = x_eq, y = y_eq, color = equilibrium)) +
coord_fixed() +
scale_fill_gradientn(colours = colorRamps::green2red(nbins), limits = c(0,1)) +
xlab("x") + ylab("y") + ggtitle("Error map") +
theme_bw()
grid.arrange(plotV, plotErr, ncol = 2)
Selkov equation
The Selkov model for glycolysis reads like:
$$
\begin{cases}
\frac{dx}{dt} & = -x + ay + x^2 y \
\frac{dy}{dt} & = b - a y - x^2 y
\end{cases}
$$
where $x$ and $y$ represent the concentrations of two chemicals. If we fix $a = 0.1$, such a system shows a limit cycle for $b \in [0.42, 0.79]$. At each side of this interval, a Hopf bifurcation happens.
We can code it in vector form:
# Parameters
a <- 0.1
b <- 0.5
# Dynamics
f <- function(x) {c(-x[1] + a*x[2] + x[1]^2*x[2],
b - a*x[2] - x[1]^2*x[2])}
This system is interesting because the Jacobian is independent of $b$, so our error map will remain constant along the bifurcation. Let's see what happens:
# Package desolve requires a slightly different syntax
f_dyn <- function(t, state, parameters) {
with(as.list(c(state, parameters)),{
# rate of change
df <- f(state)
dX <- df[1]
dY <- df[2]
# return the rate of change
list(c(dX, dY))
}) # end with(as.list ...
}
roi <- c(0, 2.5, 0, 2.5)
init_state <- c(1, .05)
ts <- seq(0, 1000, by = 0.01)
bs <- c(0.1, 0.6, 1.3)
for (b in bs) {
out <- ode(y = init_state, times = ts, func = f_dyn, parms = c(a = a, b = b))
colnames(out) <- c("time", "x", "y")
out <- as.data.frame(out)
xs <- seq(roi[1], roi[2], by = 0.05)
ys <- seq(roi[3], roi[4], by = 0.05)
result <- approxPot2D(f, xs, ys)
# Get the limit cycle attractor
attr <- dplyr::filter(as.data.frame(out), time > 0)
# Transform result into dataframe
data <- expand.grid(X = xs, Y = ys)
data$V <- as.vector(result$V)
data$err <- as.vector(result$err)
nbins <- 15
plotV <- ggplot() +
geom_tile(data = data, aes(x = X, y = Y, fill = V)) +
geom_contour(data = data, aes(x = X, y = Y, z = V), colour = 'white', alpha = 0.5, bins = nbins) +
geom_path(data = attr, aes(x = x, y = y)) +
coord_fixed() +
scale_fill_gradientn(colours = colorRamps::matlab.like(nbins)) +
xlab("x") + ylab("y") + ggtitle("Approximate potential") +
theme_bw()
plotErr <- ggplot() +
geom_tile(data = data, aes(x = X, y = Y, fill = err)) +
geom_contour(data = data, aes(x = X, y = Y, z = err), colour = 'white', alpha = 0.5, bins = nbins) +
geom_path(data = attr, aes(x = x, y = y)) +
coord_fixed() +
scale_fill_gradientn(colours = colorRamps::green2red(nbins), limits = c(0,1)) +
xlab("x") + ylab("y") + ggtitle(sprintf("Error map. b = %.3f ", b)) +
theme_bw()
grid.arrange(plotV, plotErr, ncol = 2)
}
Lotka-Volterra predator prey dynamics
Here we will use a variation of the classical Lotka-Volterra predator prey model. Particularly, the Rosenzweig-MacArthur model, that adds a function $g$ accounting for predator saturation and a carrying capacity $k$ to the prey growth.
The dynamics, being $x$ the prey biomass and $y$ the predator biomass, look like:
$$
\begin{cases}
\frac{dx}{dt} = f(x,y) = r x (1 - \frac{x}{k}) - g(x) x y \
\frac{dy}{dt} = g(x,y) = e g(x) x y - m y
\end{cases}
$$
with $g(x)$ being a saturation function:
$$
g(x) = \frac{1}{h + x}
$$
We can code it in vector form:
# Parameters
r <- 1
k <- 10
h <- 2
e <- 0.2
m <- 0.1
# Auxiliary function
g <- function(x) {1/(h + x)}
# Dynamics
f <- function(x) {c(r*x[1]*(1 - x[1]/k) -g(x[1])*x[1]*x[2],
e*g(x[1])*x[1]*x[2] - m*x[2])}
Such a system has a limit cycle attractor, as we can see simulating one of its trajectories (particularly, the one beginning at $x = 1$ and $y = 2$):
# Package desolve requires a slightly different syntax
f_dyn <- function(t, state, parameters) {
with(as.list(c(state, parameters)),{
# rate of change
df <- f(state)
dX <- df[1]
dY <- df[2]
# return the rate of change
list(c(dX, dY))
}) # end with(as.list ...
}
parms <- c(r =r,
k = k,
h = h,
e = e,
m = m)
init_state <- c(1,2)
ts <- seq(0, 300, by = 0.01)
out <- ode(y = init_state, times = ts, func = f_dyn, parms = parms)
colnames(out) <- c("time", "x", "y")
out <- as.data.frame(out)
plot(out$x, out$y, type = 'l', asp = 1,
main = 'Trajectory', xlab = 'x (prey biomass)', ylab = 'y (predator biomass)')
For such a system, we expect our pseudopotential to have a high error over large portions of the phase plane. Let's check it. First, we need to define our region of interest:
xs <- seq(0, 10, by = 0.05)
ys <- seq(0, 5, by = 0.05)
And then we are ready to apply our algorithm:
result <- approxPot2D(f, xs, ys)
# Get the limit cycle attractor
attr <- dplyr::filter(as.data.frame(out), time > 200)
# Transform result into dataframe
data <- expand.grid(X = xs, Y = ys)
data$V <- as.vector(result$V)
data$err <- as.vector(result$err)
# Input equilibrium points (calculated externally)
eqPoints <- data.frame(x_eq = c(0),
y_eq = c(0),
equilibrium = factor(c('unstable')))
Even for highly non-gradient systems, our method will compute some pseudopotential. By plotting the estimated error, we can see that it is high almost everywhere. This means that our algorithm noticed that the system is highly non-gradient, and thus, the previously computed pseudopotential is of very limited use.
nbins <- 15
plotV <- ggplot() +
geom_tile(data = data, aes(x = X, y = Y, fill = V)) +
geom_contour(data = data, aes(x = X, y = Y, z = V), colour = 'white', alpha = 0.5, bins = nbins) +
geom_point(data = eqPoints, aes(x = x_eq, y = y_eq, color = equilibrium)) +
geom_path(data = attr, aes(x = x, y = y)) +
coord_fixed() +
scale_fill_gradientn(colours = colorRamps::matlab.like(nbins)) +
xlab("x") + ylab("y") + ggtitle("Approximate potential") +
theme_bw()
plotErr <- ggplot() +
geom_tile(data = data, aes(x = X, y = Y, fill = err)) +
geom_contour(data = data, aes(x = X, y = Y, z = err), colour = 'white', alpha = 0.5, bins = nbins) +
geom_point(data = eqPoints, aes(x = x_eq, y = y_eq, color = equilibrium)) +
geom_path(data = attr, aes(x = x, y = y)) +
coord_fixed() +
scale_fill_gradientn(colours = colorRamps::green2red(nbins), limits = c(0,1)) +
xlab("x") + ylab("y") + ggtitle("Error map") +
theme_bw()
grid.arrange(plotV, plotErr, ncol = 2)
PabRod/rolldown documentation built on Dec. 22, 2020, 6:52 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo = TRUE)
knitr::include_graphics('img/logo.png')
Installation
To install the package, please type
devtools::install_github("PabRod/rolldown", ref = "master")
in your R console.
Libraries
In this vignette we will use the following libraries:
library(rolldown) # To calculate some trajectories library(deSolve) # To plot our results library(ggplot2) # To arrange our plots in panels library(latticeExtra) library(gridExtra) # For nicer plots library(colorRamps)
One dimensional examples
Allee effect
A single-species population dynamics model with Allee effect is governed by the following differential equation:
$$ \frac{dN}{dt} = r N \left( \frac{N}{A} - 1 \right) \left( 1 - \frac{N}{K} \right) $$
It is easy to see that this differential equation has three equilibrium points, $N = 0$, $N = K$ and $N = A$, being all of them stable but the latter one, which is unstable. We'll use the parameters $r = 1$, $A = 0.5$ and $K = 1$.
r <- 1 A <- 0.5 K <- 1 f <- function(x) { r * x * (x/A - 1) * (1 - x/K) }
We can use our method approxPot1D to approximate the potential function at a set of points. First, we have to create the points.
xs <- seq(0, 1.25, by = 0.01)
and then pass them to our algorithm:
Vs <- approxPot1D(f, xs)
By plotting the result, we clearly see that the two stable equilibria appear at $N = 0$ and $N = K = 1$, and the unstable one at $N = A = 0.5$, as we expected.
plot(xs, Vs, type = 'l', xlab = 'N', ylab = 'V')
Two dimensional examples
In this section we'll apply our method to a collection two dimensional systems.
Synthetic examples
We generated some abstract, synthetic examples in order to test our method. Here we present some of them.
A gradient system: the four well potential
In this section we'll deal with the two-dimensional differential equation given by:
$$ \begin{cases} \frac{dx}{dt} = f(x,y) = -x(x^2-1) \ \frac{dy}{dt} = g(x,y) = -y(y^2-1) \end{cases} $$
This is a gradient system (because $\frac{\partial f}{\partial x} = \frac{\partial g}{\partial y}$ everywhere). This means that the gradient - curl decomposition will have zero curl term everywhere and, thus, there exists a well defined potential. Particularly, the potential can be analytically proven to be:
$$ V(x,y) = \frac{x^2}{4}(x^2 - 2) + \frac{y^2}{4}(y^2 - 2) + V_0 $$
Let's try to compute it using our algorithm. First, we'll code our function as a vector:
f <- function(x) {c(-x[1]*(x[1]^2 - 1), -x[2]*(x[2]^2 - 1))}
Our region of interest is now two-dimensional. We need, thus, two vectors to create our grid of points:
xs <- seq(-1.5, 1.5, by = 0.025) ys <- seq(-1.5, 1.5, by = 0.025)
Now we are ready to apply approxPot2D:
result <- approxPot2D(f, xs, ys)
# Transform result into dataframe data <- expand.grid(X = xs, Y = ys) data$V <- as.vector(result$V) data$err <- as.vector(result$err) # Input equilibrium points (calculated externally) eqPoints <- data.frame(x_eq = c(-1, -1, 0, 1, 1), y_eq = c(-1, 1, 0, -1, 1), equilibrium = factor(c('stable', 'stable', 'unstable', 'stable', 'stable')))
result is a list that contains two fields:
result$Vcontains the estimated values of the potentials at each grid pointresult$errcontains the estimated error at each grid point
By plotting them we see:
nbins <- 15 plotV <- ggplot() + geom_tile(data = data, aes(x = X, y = Y, fill = V)) + geom_contour(data = data, aes(x = X, y = Y, z = V), colour = 'white', alpha = 0.5, bins = nbins) + geom_point(data = eqPoints, aes(x = x_eq, y = y_eq, color = equilibrium)) + coord_fixed() + scale_fill_gradientn(colours = colorRamps::matlab.like(nbins)) + xlab("x") + ylab("y") + ggtitle("Approximate potential") + theme_bw() plotErr <- ggplot() + geom_tile(data = data, aes(x = X, y = Y, fill = err)) + # geom_contour(data = data, aes(x = X, y = Y, z = err), colour = 'white', alpha = 0.5, bins = nbins) + geom_point(data = eqPoints, aes(x = x_eq, y = y_eq, color = equilibrium)) + coord_fixed() + scale_fill_gradientn(colours = colorRamps::green2red(nbins), limits = c(0,1)) + xlab("x") + ylab("y") + ggtitle("Error map") + theme_bw() grid.arrange(plotV, plotErr, ncol = 2)
Provided our example is a purely gradient system, it should not surprise us that the error is zero everywhere.
max(result$err) == 0
A non-gradient system
In this example we will apply our algorithm to the system given below:
$$ \begin{cases} \frac{dx}{dt} = f(x,y) = -y \ \frac{dy}{dt} = g(x,y) = x \end{cases} $$ This is an extreme case. The gradient - curl decomposition will give us zero gradient part everywhere. Let's feed our algorithm with these example to see what happens:
First, we code the dynamics in vector form:
f <- function(x) {c(-x[2], x[1])}
Secondly, we define our region of interest:
xs <- seq(-2, 2, by = 0.05) ys <- seq(-2, 2, by = 0.05)
And then we are ready to apply our algorithm:
result <- approxPot2D(f, xs, ys)
# Transform result into dataframe data <- expand.grid(X = xs, Y = ys) data$V <- as.vector(result$V) data$err <- as.vector(result$err) # Input equilibrium points (calculated externally) eqPoints <- data.frame(x_eq = c(0), y_eq = c(0), equilibrium = factor(c('unstable')))
The resulting approximate potential is plotted below. Being this a purely non-gradient system we expect our pseudopotential to not be trustworthy. By calculating the error we can see that actually that's the case.
nbins <- 15 plotV <- ggplot() + geom_tile(data = data, aes(x = X, y = Y, fill = V)) + geom_contour(data = data, aes(x = X, y = Y, z = V), colour = 'white', alpha = 0.5, bins = nbins) + geom_point(data = eqPoints, aes(x = x_eq, y = y_eq, color = equilibrium)) + coord_fixed() + scale_fill_gradientn(colours = colorRamps::matlab.like(nbins)) + xlab("x") + ylab("y") + ggtitle("Approximate potential") + theme_bw() plotErr <- ggplot() + geom_tile(data = data, aes(x = X, y = Y, fill = err)) + # geom_contour(data = data, aes(x = X, y = Y, z = err), colour = 'white', alpha = 0.5, bins = nbins) + geom_point(data = eqPoints, aes(x = x_eq, y = y_eq, color = equilibrium)) + coord_fixed() + scale_fill_gradientn(colours = colorRamps::green2red(nbins), limits = c(0,1)) + xlab("x") + ylab("y") + ggtitle("Error map") + theme_bw() grid.arrange(plotV, plotErr, ncol = 2)
The fact that the underlying equations are non-gradient has been captured by the algorithm.
Biological examples
Here we apply our methods to some dynamical equations well known in biology. While the abstract equations in the previous sections can be manipulated to increase or decrease their curl to gradient ratio, equations describing natural dynamical systems don't allow such a manipulation. Once again, the error map will let us know if our system allows a pseudopotential or not.
Simple regulatory gene network
A bistable network cell fate model can be described by the set of equations:
$$ \begin{cases} \frac{dx}{dt} = f(x,y) = b_x - r_x x + \frac{a_x}{k_x + y^n} \ \frac{dy}{dt} = g(x,y) = b_y - r_y y + \frac{a_y}{k_y + x^n} \end{cases} $$
Such a system represents two genes ($x$ and $y$) that inhibit each other. This circuit works as a toggle switch with two stable steady states, one with dominant $x$ , the other with dominant $y$.
We can code it in vector form:
# Parameters bx <- 0.2 ax <- 0.125 kx <- 0.0625 rx <- 1 by <- 0.05 ay <- 0.1094 ky <- 0.0625 ry <- 1 n <- 4 # Dynamics f <- function(x) {c(bx - rx*x[1] + ax/(kx + x[2]^n), by - ry*x[2] + ay/(ky + x[1]^n))}
This set of equations is, in general, not gradient (because $\frac{\partial f}{\partial y} \neq \frac{\partial g}{\partial x}$). Anyways, we can use the method approxPot2D to compute the approximate potential.
First, we need to define our region of interest:
xs <- seq(0, 4, by = 0.05) ys <- seq(0, 4, by = 0.05)
And then we are ready to apply our algorithm:
result <- approxPot2D(f, xs, ys)
# Transform result into dataframe data <- expand.grid(X = xs, Y = ys) data$V <- as.vector(result$V) data$err <- as.vector(result$err) # Input equilibrium points (calculated externally) # # Estimated with Wolfram Alpha # Prompt: 0 = 0.2 - x + 0.125/(0.0625 + y^4); 0 = 0.05 - y + 0.1094/(0.0625 + x^4) eqPoints <- data.frame(x_eq = c(0.213416, 0.559865, 2.19971), y_eq = c(1.74417, 0.730558, 0.0546602), equilibrium = factor(c('stable', 'unstable', 'stable')))
The resulting approximate potential is plotted below. Being this not a gradient system it is advisable to plot also the estimated error. The areas in green represent small approximation error, so the potential can be safely used in those regions.
nbins <- 25 plotV <- ggplot() + geom_tile(data = data, aes(x = X, y = Y, fill = V)) + geom_contour(data = data, aes(x = X, y = Y, z = V), colour = 'white', alpha = 0.5, bins = nbins) + geom_point(data = eqPoints, aes(x = x_eq, y = y_eq, color = equilibrium)) + coord_fixed() + scale_fill_gradientn(colours = colorRamps::matlab.like(nbins)) + xlab("x") + ylab("y") + ggtitle("Approximate potential") + theme_bw() plotErr <- ggplot() + geom_tile(data = data, aes(x = X, y = Y, fill = err)) + geom_contour(data = data, aes(x = X, y = Y, z = err), colour = 'white', alpha = 0.5, bins = nbins) + geom_point(data = eqPoints, aes(x = x_eq, y = y_eq, color = equilibrium)) + coord_fixed() + scale_fill_gradientn(colours = colorRamps::green2red(nbins), limits = c(0,1)) + xlab("x") + ylab("y") + ggtitle("Error map") + theme_bw() grid.arrange(plotV, plotErr, ncol = 2)
Selkov equation
The Selkov model for glycolysis reads like:
$$ \begin{cases} \frac{dx}{dt} & = -x + ay + x^2 y \ \frac{dy}{dt} & = b - a y - x^2 y \end{cases} $$
where $x$ and $y$ represent the concentrations of two chemicals. If we fix $a = 0.1$, such a system shows a limit cycle for $b \in [0.42, 0.79]$. At each side of this interval, a Hopf bifurcation happens.
We can code it in vector form:
# Parameters a <- 0.1 b <- 0.5 # Dynamics f <- function(x) {c(-x[1] + a*x[2] + x[1]^2*x[2], b - a*x[2] - x[1]^2*x[2])}
This system is interesting because the Jacobian is independent of $b$, so our error map will remain constant along the bifurcation. Let's see what happens:
# Package desolve requires a slightly different syntax f_dyn <- function(t, state, parameters) { with(as.list(c(state, parameters)),{ # rate of change df <- f(state) dX <- df[1] dY <- df[2] # return the rate of change list(c(dX, dY)) }) # end with(as.list ... } roi <- c(0, 2.5, 0, 2.5) init_state <- c(1, .05) ts <- seq(0, 1000, by = 0.01) bs <- c(0.1, 0.6, 1.3) for (b in bs) { out <- ode(y = init_state, times = ts, func = f_dyn, parms = c(a = a, b = b)) colnames(out) <- c("time", "x", "y") out <- as.data.frame(out) xs <- seq(roi[1], roi[2], by = 0.05) ys <- seq(roi[3], roi[4], by = 0.05) result <- approxPot2D(f, xs, ys) # Get the limit cycle attractor attr <- dplyr::filter(as.data.frame(out), time > 0) # Transform result into dataframe data <- expand.grid(X = xs, Y = ys) data$V <- as.vector(result$V) data$err <- as.vector(result$err) nbins <- 15 plotV <- ggplot() + geom_tile(data = data, aes(x = X, y = Y, fill = V)) + geom_contour(data = data, aes(x = X, y = Y, z = V), colour = 'white', alpha = 0.5, bins = nbins) + geom_path(data = attr, aes(x = x, y = y)) + coord_fixed() + scale_fill_gradientn(colours = colorRamps::matlab.like(nbins)) + xlab("x") + ylab("y") + ggtitle("Approximate potential") + theme_bw() plotErr <- ggplot() + geom_tile(data = data, aes(x = X, y = Y, fill = err)) + geom_contour(data = data, aes(x = X, y = Y, z = err), colour = 'white', alpha = 0.5, bins = nbins) + geom_path(data = attr, aes(x = x, y = y)) + coord_fixed() + scale_fill_gradientn(colours = colorRamps::green2red(nbins), limits = c(0,1)) + xlab("x") + ylab("y") + ggtitle(sprintf("Error map. b = %.3f ", b)) + theme_bw() grid.arrange(plotV, plotErr, ncol = 2) }
Lotka-Volterra predator prey dynamics
Here we will use a variation of the classical Lotka-Volterra predator prey model. Particularly, the Rosenzweig-MacArthur model, that adds a function $g$ accounting for predator saturation and a carrying capacity $k$ to the prey growth.
The dynamics, being $x$ the prey biomass and $y$ the predator biomass, look like:
$$ \begin{cases} \frac{dx}{dt} = f(x,y) = r x (1 - \frac{x}{k}) - g(x) x y \ \frac{dy}{dt} = g(x,y) = e g(x) x y - m y \end{cases} $$
with $g(x)$ being a saturation function:
$$ g(x) = \frac{1}{h + x} $$
We can code it in vector form:
# Parameters r <- 1 k <- 10 h <- 2 e <- 0.2 m <- 0.1 # Auxiliary function g <- function(x) {1/(h + x)} # Dynamics f <- function(x) {c(r*x[1]*(1 - x[1]/k) -g(x[1])*x[1]*x[2], e*g(x[1])*x[1]*x[2] - m*x[2])}
Such a system has a limit cycle attractor, as we can see simulating one of its trajectories (particularly, the one beginning at $x = 1$ and $y = 2$):
# Package desolve requires a slightly different syntax f_dyn <- function(t, state, parameters) { with(as.list(c(state, parameters)),{ # rate of change df <- f(state) dX <- df[1] dY <- df[2] # return the rate of change list(c(dX, dY)) }) # end with(as.list ... } parms <- c(r =r, k = k, h = h, e = e, m = m) init_state <- c(1,2) ts <- seq(0, 300, by = 0.01) out <- ode(y = init_state, times = ts, func = f_dyn, parms = parms) colnames(out) <- c("time", "x", "y") out <- as.data.frame(out) plot(out$x, out$y, type = 'l', asp = 1, main = 'Trajectory', xlab = 'x (prey biomass)', ylab = 'y (predator biomass)')
For such a system, we expect our pseudopotential to have a high error over large portions of the phase plane. Let's check it. First, we need to define our region of interest:
xs <- seq(0, 10, by = 0.05) ys <- seq(0, 5, by = 0.05)
And then we are ready to apply our algorithm:
result <- approxPot2D(f, xs, ys)
# Get the limit cycle attractor attr <- dplyr::filter(as.data.frame(out), time > 200) # Transform result into dataframe data <- expand.grid(X = xs, Y = ys) data$V <- as.vector(result$V) data$err <- as.vector(result$err) # Input equilibrium points (calculated externally) eqPoints <- data.frame(x_eq = c(0), y_eq = c(0), equilibrium = factor(c('unstable')))
Even for highly non-gradient systems, our method will compute some pseudopotential. By plotting the estimated error, we can see that it is high almost everywhere. This means that our algorithm noticed that the system is highly non-gradient, and thus, the previously computed pseudopotential is of very limited use.
nbins <- 15 plotV <- ggplot() + geom_tile(data = data, aes(x = X, y = Y, fill = V)) + geom_contour(data = data, aes(x = X, y = Y, z = V), colour = 'white', alpha = 0.5, bins = nbins) + geom_point(data = eqPoints, aes(x = x_eq, y = y_eq, color = equilibrium)) + geom_path(data = attr, aes(x = x, y = y)) + coord_fixed() + scale_fill_gradientn(colours = colorRamps::matlab.like(nbins)) + xlab("x") + ylab("y") + ggtitle("Approximate potential") + theme_bw() plotErr <- ggplot() + geom_tile(data = data, aes(x = X, y = Y, fill = err)) + geom_contour(data = data, aes(x = X, y = Y, z = err), colour = 'white', alpha = 0.5, bins = nbins) + geom_point(data = eqPoints, aes(x = x_eq, y = y_eq, color = equilibrium)) + geom_path(data = attr, aes(x = x, y = y)) + coord_fixed() + scale_fill_gradientn(colours = colorRamps::green2red(nbins), limits = c(0,1)) + xlab("x") + ylab("y") + ggtitle("Error map") + theme_bw() grid.arrange(plotV, plotErr, ncol = 2)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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.