In this tutorial, Gui explains how to perform Linear Stability Analysis in the way that many classical theoretical papers in ecology deal with this type of stability in the context of food webs.
Aims
Our main goal is to understand the elements of the linear stability analysis of a multidimensional dynamical system. This requires an intuition of elements such as the Jacobian and the eigenvalues of a matrix. Our task is divided in 3 steps:
-
Understanding the reasoning of the procedure by looking at the 1D system.
-
Developing the 2D form analogously to the 1D form.
-
Calculating the eigenvalues, which will allow us to perform the stability analysis in 2D.
1D Dynamical system
Our variable of interest is $x(t)$, which is a curve parameterised by the temporal variable $t$. The time-derivative of $x$, $\dot{x}=\frac{dx}{dt}$, is another curve parameterised by $t$. In a dynamical system, we interpret $x$ and $\dot{x}$ as two different functions of $t$ inhabiting the same space. The usual dynamical law to determine $x$ is to define an operation that transforms $x$ into $\dot{x}$:
\[f:x\to\dot{x}\]
such that
\[\dot{x}=f(x).\]
Note that this equation doesn’t explicitly involve $t$, which can be taken as an external parameterisation. We can then interpret $x$ and $\dot{x}$ as 1D vectors with the operator $f$ as a geometrical shift from $x$ to $\dot{x}$, independently of $t$.
If we are interested in the behaviour around an equilibrium point $x_0$, we can linearise $\dot{x}$ as a function of $x$ around that point. This is the same as the curved surface of Earth being ‘linearised’ around us in the very small region delimited by our line of sight. For this, we consider the Taylor series expansion of the operator around the equilibrium $x_0$:
\[\dot{x}=f(x)=f(x_0)+\frac{df}{dx}(x_0)[x-x_0]+\frac{1}{2}\frac{d^f}{dx^2}(x_0)[x-x_0]^2+...=\sum_n\frac{d^nf}{dx^n}(x_0)\frac{[x-x_0]^n}{n!}.\]
All derivatives are evaluated at equilibrium, so they are numbers, constant coefficients. Therefore, this expansion is a mapping from $f$ to an infinite polynomial, which is a recipe allowing us to write a function as a polynomial. At equilibrium, we have $\dot{x}(x_0)=f(x_0)=0$. Then, since we are very close to the equilibrium, the distance $[x-x_0]$ is very close to zero. The process of linearisation is to approximate to zero all higher orders $[x-x_0]^n$ for $n>1$. Thus, the linearized system, which is valid only if we are close to equilibrium, becomes:
\[\dot{x}=f(x)=\frac{df}{dx}(x_0)[x-x_0]=\lambda[x-x_0].\]
We write \(\frac{df}{dx}(x_0)=\lambda\) because it is a constant coefficient anyway. Then, we can define $z=[x-x_0]$ (with $\dot{z}=\dot{x}$, because $x_0$ is a constant) as the displacement from the equilibrium. Then, we solve this system for $z(t)$, which is a simple exponential solution:
\[\frac{dz}{dt}=\lambda z\]
\[z(t)=z(0)e^{\lambda t}.\]
Then, if we perturb the equilibrium with a displacement $z(0)$, will the system return to equilibrium or will it diverge away from it? It all depends on the derivative $\lambda$. If $\lambda<0$, then the equilibrium will be restored, making this equilibrium stable (what happens if $\lambda=0$?). Therefore, we have the condition for stability of the equilibrium $x_0$:
\[\frac{df}{dx}(x_0)<0,\]
which means that the derivative of the curve defining $\dot{x}$ as a function of $x$, when evaluated at the equilibrium, determines its stability.
As an example, consider the logistic growth:
\[\dot{x}=rx(1-\frac{x}{K}).\]
It has an equilibrium $x_0=K$. The linearised system around this equilibrium is
\[\dot{x}=-r[x-K].\]
Then, if $-r<0$, this equilibrium is stable. That is $r>0$.
2D Dynamical system
Now, we do the same for a 2D system. Consider the vector
\[\begin{align}
v &= \begin{bmatrix}
x \\
y \\
\end{bmatrix}.
\end{align}\]
A dynamical law connects $v$ with its derivative $\dot{v}$:
\[\begin{align}
\dot{v} &= \begin{bmatrix}
f(x,y) \\
g(x,y) \\
\end{bmatrix}.
\end{align}\]
We can use the same expansion on both entries to linearise the system around the equilibrium $(x_0,y_0)$. But now the maps $f$ and $g$ are functions of two variables. Therefore, they can change by varying both $x$ and $y$.
Now we have to introduce partial derivatives. Suppose we want to produce a variation $df$. For that, we cause a variation $dx$ and a variation $dy$, giving for $f$ two independent sources of variation. The partial derivatives are the coefficient of each source, geometrically meaning the tangents of the shadows of $f$ projected on $f-x$ and $f-y$ planes. Thus,
\[df = \frac{\partial f}{\partial x}dx + \frac{\partial f}{\partial y}dy.\]
This means that the polynomial expansion for two variables will consider both partial derivatives, computing perturbations away from equilibrium for both dimensions:
\[\begin{align}
\dot{v} &= \begin{bmatrix}
f(x_0,y_0)+\frac{\partial f}{\partial x}(x_0,y_0)[x-x_0]+\frac{\partial f}{\partial y}(x_0,y_0)[y-y_0]+... \\
g(x_0,y_0)+\frac{\partial g}{\partial x}(x_0,y_0)[x-x_0]+\frac{\partial g}{\partial y}(x_0,y_0)[y-y_0]+... \\
\end{bmatrix}.
\end{align}\]
Now we impose that $\dot{v}$ is zero at equilibrium, so $f(x_0,y_0)=g(x_0,y_0)=0$, then we define the vector of displacement from equilibrium
\[\begin{align}
z &= \begin{bmatrix}
[x-x_0] \\
[y-y_0] \\
\end{bmatrix}.
\end{align}\]
Then, the linearised system becomes
\[\begin{align}
\dot{z} &= \begin{bmatrix}
\frac{\partial f}{\partial x}(x_0,y_0)[x-x_0]+\frac{\partial f}{\partial y}(x_0,y_0)[y-y_0] \\
\frac{\partial g}{\partial x}(x_0,y_0)[x-x_0]+\frac{\partial g}{\partial y}(x_0,y_0)[y-y_0] \\
\end{bmatrix}=Jz,
\end{align}\]
where we define the Jacobian matrix of the system at this point as
\[\begin{align}
J &= \begin{bmatrix}
\frac{\partial f}{\partial x}(x_0,y_0)\quad \frac{\partial f}{\partial y}(x_0,y_0) \\
\frac{\partial g}{\partial x}(x_0,y_0)\quad \frac{\partial g}{\partial y}(x_0,y_0) \\
\end{bmatrix}.
\end{align}\]
The Jacobian is the generalisation of the derivative for multidimensional systems. Instead of a single number, it is a matrix comprising all partial derivatives of the system. Here, we present the Jacobian evaluated at the equilibrium point, which is also called the community matrix, in the context of community ecology. Note that the Jacobian itself is a function of $(x,y)$, but when evaluated at the equilibrium it becomes a simple constant matrix.
Now, in the same way as the 1D system, we can expect the solution of the linearised system $\dot{z}=Jz$ to be an exponential:
\[z(t) = ce^{\lambda t},\]
where $c$ is a vector of initial conditions and $\lambda$ is a number determining the behaviour of the linearized system around the equilibrium. For this case, how can we know the identities of $c$ and $\lambda$? As we will see, they are not single objects, they are the eigenvectors and eigenvalues of $J$.
Eigenvalue calculation
If we derive the exponential solution and identify with the dynamical rule $Jz$, we obtain the following equation:
\[Jc = \lambda c.\]
This equation is telling us that the operator $J$ applied to the vector $c$ does not project it on another direction, but only changes its size by a factor $\lambda$. From this, we can build a projection to zero:
\[(J-\lambda I)c=0,\]
where $I$ is the identity matrix. Since the vector $c$ is not zero, it follows that the matrix operator $(J-\lambda I)$ is doing something special: it’s taking $c$ to zero, regardless of what $c$ is, which is like an algebraic black hole. Therefore, this matrix is a singular matrix, it causes an indetermination, promoting ‘loss of information’ as it brings the entire vector space to zero. The condition for singular matrices is to have a zero determinant, which is the equivalent of a scalar operator being zero in a 1D system.
Why the determinant? We can think of the determinant of a matrix as a type of scalar substitute for it. More precisely, it is the area scaling of the transformation carried out by the matrix, how much the transformation expands and contracts the vector space. We can even derive the determinant formula for a 2D system from the area scaling of a paralelogram defined by two vectors being transformed by the matrix. The determinant for 2D matrices is:
\[\begin{align}
det[\begin{bmatrix}
\quad a\quad\quad b\quad \\
\quad c\quad\quad d\quad \\
\end{bmatrix}]=ad-bc.
\end{align}\]
Then, the task of calculating $\lambda$ becomes finding the roots of the following equation:
\[\begin{align}
det[\begin{bmatrix}
\quad a-\lambda\quad\quad b\quad \\
\quad c\quad\quad d-\lambda\quad \\
\end{bmatrix}]=\lambda^2-\lambda(a+d)+ad-bc=0,
\end{align}\]
where $\frac{\partial f}{\partial x}(x_0,y_0)=a$, $\frac{\partial f}{\partial y}(x_0,y_0)=b$, $\frac{\partial g}{\partial x}(x_0,y_0)=c$, $\frac{\partial g}{\partial y}(x_0,y_0)=d$. We can rewrite this equation in a way that hints to its generalization to higher dimensions involving the trace and the determinant of $J$:
\[\lambda^2-Tr(J)\lambda+det(J)=0\]
Since this is a second-order equation, there are two solutions: $\lambda_1$ and $\lambda_2$. For each of these eigenvalues, there is an associated eigenvector that we calculate from $Jc_1=\lambda_1c_1$ (and the same for $\lambda_2$ and $c_2$), of which only the direction is important, not the magnitude.
Once we have the eigenvalues and eigenvectors, the complete solution for the system of differential equations $\dot{z}=Jz$ is a linear combination of the solutions given by each pair $(\lambda,c)$. You can verify that the linear combination is also a solution by deriving it (or by noting that it has to be, since it’s a linear system). Then:
\[z(t) = a_1c_1e^{\lambda_1 t} + a_2c_2e^{\lambda_2 t},\]
where $a_1$ and $a_2$ are constant numbers we determine from the initial state system $z(0)=a_1c_1+a_2c_2$.
Given the final solution, we can assess the stability of the equilibrium. From an initial displacement $z(0)$, if all eigenvalues $(\lambda_1,\lambda_2)$ are negative, the system returns to the equilibrium. If any of the eigenvalues is positive, even if the other components return to the equilibrium value, at least one diverges away. Therefore, the equilibrium is stable if all eigenvalues are negative.
Let’s show an example with a simplified version of the classic Lotka-Volterra predator-prey system:
\[\begin{align}
\dot{v} &= \begin{bmatrix}
x-\alpha xy \\
\alpha xy -y \\
\end{bmatrix}.
\end{align}\]
The Jacobian of the system is
\[\begin{align}
J(x,y)=\begin{bmatrix}
\quad 1-\alpha y\quad\quad -\alpha x\quad \\
\quad \alpha y\quad\quad 1-\alpha x\quad \\
\end{bmatrix}.
\end{align}\]
This system has an equilibrium point $x_0=y_0=1/\alpha$. The Jacobian at this point is then
\[\begin{align}
J=\begin{bmatrix}
\quad 0\quad\quad -1\quad \\
\quad 1\quad\quad 0\quad \\
\end{bmatrix}.
\end{align}\]
The characteristic equation for the eigenvalues is then
\[\lambda^2+1=0.\]
The solutions are imaginary: $\lambda_1=i$ and $\lambda_2=-i$. What happens if the eigenvalues have imaginary components? They represent oscillations, so it’s neither stable nor unstable. The associated eigenvectors are the directions
\[\begin{align}
c_1=\begin{bmatrix}
1 \\
i \\
\end{bmatrix}, \quad c_2=\begin{bmatrix}
i \\
1 \\
\end{bmatrix}.
\end{align}\]
Therefore, we can write the solution for linearized displacements from the equilibrium as
\[\begin{align}
z(t)=a_1\begin{bmatrix}
1 \\
i \\
\end{bmatrix}e^{it} + a_2\begin{bmatrix}
i \\
1 \\
\end{bmatrix}e^{-it}.
\end{align}\]
How can the solution be complex if $z(t)$ is real? The answer is that the solution is, in fact, real. It can be rewritten in terms of real sines and cosines (using Euler’s formula), and this also explains why imaginary parts of eigenvalues translate into oscillations.
Complex eigenvalues can be divided into real and imaginary parts. The real parts will determine the stability, so the criteria actually states that the largest real part of an eigenvalue has to be negative for the equilibrium to be stable. The imaginary parts determine oscillations. That’s because, in the end, the real solution features real parts of the eigenvalues leading the exponentials and the imaginary parts transformed into sines and cosines. In the case of the Lotka-Volterra above, the real parts are all zero, so the equilibrium is actually neutral. Thus, the system oscillates around the equilibrium with neutral amplitude (it depends on the size of the perturbation).
In this tutorial, Marco explains how to implement the main routines to replicate the effects of the paper on Noise-induced effects in collective dynamics and inferring local interactions from data
Part 1 Gillespie algorithm
The code contains a function for exactly simulating runs from the master equation through the Gillespie algorithm. The Gillespie algorithm is divided in four main parts:
1.- Find the velocity (or propensity) for every reaction. This is done by multiplying the reaction rate (which is conditioned on the reagents encountering each other) with the rate that the reagents do encounter each other. The rate that two reagents encounter each other can be found through the mass action assumption. Usually this is done simply by multiplying the concentration of the reagents, but this is an approximation that holds only for big population sizes because reagents do not interact with themselves, and the concentration of the species with and without one molecule can be quite different in small population sizes and omega. Hence it is possible to use the stochastic law of mass action, equal for one reaction to $ k \Omega \prod \frac{x_i!}{(x_i-1)! \Omega^{s_i}} $ (for details see David Schnoerr et al (2017) J. Phys. A: Math. Theor. 50, 093001). Since there are some factorials, it is hard to compute in practice, and has been rearranged in the code implementation.
2.- Find the waiting time before the next reaction occurs by sampling from an exponential distribution with mean rate equal to the sum of the reaction velocities.
3.- Find the reaction that happens by choosing from the possible reactions with a probability proportional to the reaction velocity.
4.- Update the states based on the reaction that happened.
The function has not been optimized but works for simulating every type of system starting from the microscopic reactions. It performs only one update and must be inserted inside a for loop for producing the entire simulation. As input, it requires:
1.- The current number of molecules/individuals in a specific state through a vector of size equal to the number of possible states and each entry the number of individuals in that state.
2.- A vector containing the reaction rates.
3.- A matrix indicating the stoichiometry of the reagents, with the number of rows equal to the number of reactions, and number of columns equal to the number of possible states, and entries equal to the number of species of that state that must encounter to produce the specific reaction.
4.- A matrix indicating the stoichiometry of the products, with the number of rows equal to the number of reactions, and number of columns equal to the number of possible states, and entries equal to the number of species of that state that are produced from that specific reaction.
5.- A parameter omega indicating the spatial scale of the system.
As output, it gives a vector containing the waiting time, the number of the reaction the occurred, and the new number of individuals/molecules in each state.
It could be a good test to modify the code to simulate clonal logistic growth by using these reactions: A -> 0 with rate 0.05, A -> 2A with rate 0.2, 2A -> A with rate 0.05. You can start multiple replicates by having only one individual and set omega to 50.
Part 2 Replication of the results of the paper “Noise-induced effects in collective dynamics and inferring local interactions from data”
Jhawar J, Guttal V. (2020) Phil. Trans. R. Soc. B 375: 20190381. http://dx.doi.org/10.1098/rstb.2019.0381
The replication of the results consists of two parts, one applied to the voter model, and one applied to the higher order interactions model (see the paper for details). The only difference between the two is that the higher order interactions also has a reaction of the type 2A + B -> 3A, which leads to a non-linear per capita rate of change in the ordinary differential equation and consensus formation with multi-stability in the limit of infinite population size (not present in the voter model). For each model it is shown examples of the dynamics through 20 replicates of the time series obtained iterating the Gillespie algorithm and the probability density of the population being in certain states. Interestingly, also the voter model shows consensus and multi-stability! This is because of the constructive effect of intrinsic noise. To better understand this result, it is possible to fit a Langevin equation to the time series we have simulated and find the functional form of the deterministic (drift) and stochastic (diffusion) terms. Despite the Langevin equation being an approximation, it explicitly defines a term for stochasticity, enabling us to investigate noise, quite cool.
To fit the data, we first have to find the autocorrelation function for figuring out the best time scale (see paper for details).
\[ACF(\tau) = \frac{\langle (M(t) - \langle M(t) \rangle ) (M(t+\tau)- \langle M(t) \rangle) \rangle}{\langle (M(t) - \langle M(t) \rangle)^{2} \rangle }\]
When I see an equation my first though is “how do I calculate that?”. In my opinion Tidyverse has some cool operations which can be used to translate the equation into a specific value. For example, the averages in the ACF are calculated over different things, which can be captured by grouping.
The last part is to find the deterministic and stochastic component, which produce respectively the first and second moments (see paper for details). The first moment is the average “jump” from a specific population state. The second moment is the variance of such “jumps”. In this case there are only two possible individual states, only one degree of freedom to describe the state of the entire population, but in more complex cases a “jump” could be a vector in a multidimensional state space.
Comparing the two models we see that the deterministic component “pushes” the population towards non-consensus in the voter model (equal number of individuals for both options), while there are two stable equilibria separated by and unstable equilibrium in the higher interaction model (equilibria are found when the “push” is equal to zero, and equilibria are stable when the “push” goes towards the equilibrium). For both models the noise is stronger when the population is undecided because there is more uncertainty in the reactions that can occur (the reactions have similar velocities), which is the definition of intrinsic noise. Hence, the voter model is pushed away from its non-consensus equilibrium condition, leading to noise induced consensus. Noise constructs order!
Marco has kindly share his code with us here.
