#80 Fundamentals--Linear Stability Analysis 常識集--線性穩定分析

(This is the second edition of this post because the previous one disappeared from the beginning to the end after I pressed ctrl+Z......Unimaginable.)

After a struggling semester, I could finally work back on 3min biophysics and share something interesting with you guys. Recently I am preparing for a series on the energy trade-off in sensory adaptation. However, I have to make sure that every reader understands some basic concepts, such as linear stability analysis, before we move on. So today we will introduce linear stability analysis using a classic system -- the Lotka-Volterra prey and predator model, with the aids of python.

Set up a system of ODE

As a biophysicist, it is tempting to describe a system using system of ODEs. How could we formulate the interactions between preys and predators using ODEs? Let's first consider what this system supposed to be. The preys should reproduce themselves otherwise they will all be eaten eventually, and there should have some terms describing the predations. The growth of predators relies on how many preys they could catch, and they will die if there is no preys at all. So the ODEs should look like this:

$\\ \dot{\textup{prey}}=\textup{reproduction}-\textup{predation}\\ \dot{\textup{predator}}=\textup{predation}-\textup{die}$

So there are many possible ways to describe the precise dynamics. One famous possibility is the one proposed by Lotka in 1910:

$\\ \dot{N}=rN-aPN\\ \dot{P}=afPN-mP$

Determine the stability of fixed points

In such systems of ODEs, there are usually some fixed points, where the derivatives of all variables vanish. It is easy to solve this point by set the above 2 equations equal to 0, and the result should be:

$(N_0,P_0)=(\frac{m}{af},\frac{r}{a})$

However, is that fixed point stable? Before answering this question, we don't even know what do we mean by "being stable." In fact, there are different ways to define stability, each with some vigorous mathematical definition. However, we will just use our intuitive thoughts on "local stability" to think through it instead. Let's say now our system is slightly deviated from our fixed point. That is:

$(N',P')=(N_0,P_0)+(\Delta N, \Delta P)$

How will this point evolve? Will it move away from the fixed point or will it move toward? This question could be easily answered by plugging it back to our system of ODEs.

$\\ \dot{N}=\dot{\Delta N}=r(N+\Delta N)-a(P+\Delta P)(N+\Delta N)\\ \dot{P}=\dot{\Delta P}=af(P+\Delta P)(N+\Delta N)-m(P+\Delta P)$

After some simplification, you should arrive at something just like this:

$\\ \dot{\Delta N}=-\frac{m}{f}\Delta P-a\Delta N\Delta P\\ \dot{\Delta P}=rf\Delta N +af\Delta N\Delta P$

~~(I want something just like this~~Oh~~Ohhh~~)~~

Since we have stated that the system is only slightly deviated from the fixed point, so we could remove all higher order terms and leave the linear terms alone. We could then restate our system in terms of matrix:

$\begin{bmatrix}\dot{\Delta N}\\\dot{\Delta P}\end{bmatrix}=\begin{bmatrix} 0 & -m/f\\rf & 0\end{bmatrix}\begin{bmatrix}\Delta N\\\Delta P\end{bmatrix}$

And the matrix in the middle determines the evolution of deviations from the fixed point. For a 2-by-2 matrix, there are usually some intrinsic vectors associated with this matrix, or the "eigenvectors." These vectors have a special properties that:

$\lambda\bold{\upright{x}}=\bold{\upright{J}}\bold{\upright{x}}$

That is, the operation of matrix J (the Jacobian matrix) on vector x will not change its directions. If we could find 2 non-parallel eigenvectors, any 2 dimensional vectors could then be described in terms of the linear combinations of these 2 eigenvectors. For example, if our deviation vectors could be described as:

$\begin{bmatrix}\Delta N\\\Delta P\end{bmatrix}=\alpha\bold{\upright{x_1}}+\beta\bold{\upright{x_2}}$

Then the evolution of the deviations from fixed point could be easily solved:

$\begin{bmatrix}\Delta N(t)\\\Delta P(t)\end{bmatrix}=\alpha\bold{\upright{x_1}}e^{\lambda_1t}+\beta\bold{\upright{x_2}}e^{\lambda_2t}$

So that λ, or the eigenvalues, determine the linear stability of fixed points.

Now back to our Lotka-Volterra problem. What is the eigenvalues of our Jacobian matrix? From the definition of eigenvectors, the following must hold for an eigenvector:

$(\bold{\upright{J}}-\lambda\bold{\upright{I}})\bold{\upright{x}}=\bold{\upright{0}}$

One obvious solution is x = 0, but that would be meaningless. For a meaningful x to exist:

$\textup{det}(\bold{\upright{J}}-\lambda\bold{\upright{I}})=0$

must be true. In our case, the eigenvalues could be solve easily:

$\lambda=\pm \sqrt{mr}i$

2 pure imaginary numbers. That means the deviation will not vanish in terms of linear stability and the system will circle around the fixed point at least temporarily. The bold parts of the previous sentence could not be omitted because we arrive at such conclusion by excluding higher order terms. Let's see if that is true. By writing some simple python codes:

The evolution of the system on phase plane would be

The system orbits around the fixed point on the phase plane.

Some modifications guided by linear stability analysis

It seems a little bit boring because no matter how we adjust the parameters, the dynamics of the system and the properties of the fixed point will be left unchanged. Let's do some modifications by introducing capacity and "logistic growth." Now we define our system as:

$\\ \dot{N}=rN(1-\frac{N}{K})-aPN\\ \dot{P}=afPN-mP$

The introduction of capacity will limit the growth rate of preys. Would it change the dynamics of the system and the properties of the fixed point? To answer that question we have to turn to our new Jacobian matrix...(you could derive it either through the same way as presented above or using differentiation and plugging in the fixed point, they will arrive at something exactly the same)

$\bold{\upright{J}}=\begin{bmatrix} -rm/afK & -m/f \\ rf-rm/aK & 0 \end{bmatrix}$

The eigenvalues could be solved:

$\lambda=\frac{-rm}{2fK}\pm \sqrt{-rm+\frac{rm^2}{afK}+\frac{r^2m^2}{4a^2f^2K^2}}$

Since K is usually a large number, the eigenvalues could be approximated as:

$\lambda\sim\frac{-rm}{2fK}\pm i\sqrt{rm}$

There is an imaginary part in the eigenvalues so the system would circle around the fixed point. However, there is a negative real part so the system will in fact spiral toward the fixed point. Let's see if it is the case. By slightly modify our code

We will find

However, this is not always the case since we have some space for adjustment in our eigenvalues. We should ignore the term with 1/K^2 because it would be way too small. Yet the 1/K term could balance out -rm if the following inequality is true:

$m>afK$

So we modify the constants to satisfy the above inequality and let's see how we could destroy this system:

The result would be

All predators extincted because they can't assimilate enough energy from preys to reproduce themselves and die too quickly. That really make sense!

Take home message

So today we use Lotka-Volterra model and its modification to demonstrate linear stability analysis and how this process shed light on the change of system dynamics with respect to some adjustment in key parameters. A more formal descriptions could be easily found in any engineering mathematics textbooks. I hope you enjoy this episode and please stay tuned!

Three Minute Biophysics

Search This Blog

Featured Post

#77 Brownian Carnot engine-I 微觀卡諾引擎-I