9b. Linear stability analysis

[1]:

# Colab setup ------------------
import os, sys, subprocess
if "google.colab" in sys.modules:
    cmd = "pip install --upgrade watermark"
    process = subprocess.Popen(cmd.split(), stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    stdout, stderr = process.communicate()
# ------------------------------

import numpy as np

import bokeh.plotting
import bokeh.io

bokeh.io.output_notebook()

Loading BokehJS ...

Thus far, we have assessed the stability of fixed points of a dynamical system graphically. We have drawn plots of production rates and degradation rates of gene products, with the crossing of these respective plots being fixed points. We then investigated when degradation was faster than production, and vice versa, to determine whether a system is attracted toward or repelled away from a fixed point. We also did similar nullcline-based analysis.

In this technical appendix, we present a more general, and very useful, approach called linear stability analysis. We first give a minimal introduction to the technique, and then apply it to the repressilator. You can read about the technique in greater depth in Strogatz’s book.

A minimal introduction to linear stability analysis

The main idea behind linear stability analysis is to locally approximate a nonlinear dynamical system by its Taylor series to first order near the fixed point, and then examine the behavior of the resulting simpler linear system. The Hartman-Grobman theorem (which we will not derive here) ensures that the linearized system faithfully represents the phase portrait of the full nonlinear system near the fixed point.

Say we have a dynamical system with variables \(\mathbf{u}\) with

\begin{align} \frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t} = \mathbf{f}(\mathbf{u}), \end{align}

where \(\mathbf{f}(\mathbf{u})\) is a vector-valued function, i.e.,

\begin{align} \mathbf{f}(\mathbf{u}) = (f_1(u_1, u_2, \ldots), f_2(u_1, u_2, \ldots), \ldots). \end{align}

Assuming there is a fixed point \(\mathbf{u}_0\), linear stability analysis proceeds with the following steps:

1. Linearize about \(\mathbf{u}_0\), defining \(\delta\mathbf{u} = \mathbf{u} - \mathbf{u}_0\). To do this, expand \(f(\mathbf{u})\) in a Taylor series about \(\mathbf{u}_0\) to first order.

\begin{align} \mathbf{f}(\mathbf{u}) = \mathbf{f}(\mathbf{u}_0) + \nabla \mathbf{f}(\mathbf{u}_0)\cdot \delta\mathbf{u} + \cdots, \end{align}

where \(\nabla \mathbf{f}(\mathbf{u}_0) \equiv \mathsf{A}\) is the Jacobi matrix evaluated at \(\mathbf{u}_0\),

\begin{align} \nabla \mathbf{f}(\mathbf{u}_0) \equiv \mathsf{A} = \left.\begin{pmatrix} \frac{\partial f_1}{\partial u_1} & \frac{\partial f_1}{\partial u_2} & \cdots \\[0.5em] \frac{\partial f_2}{\partial u_1} & \frac{\partial f_2}{\partial u_2} & \cdots \\ \vdots & \vdots & \ddots \end{pmatrix}\right|_{\mathbf{\, u}_0} \end{align}

Thus, we have

\begin{align} \frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t} = \frac{\mathrm{d}\mathbf{u}_0}{\mathrm{d}t} + \frac{\mathrm{d}\delta\mathbf{u}}{\mathrm{d}t} = \mathbf{f}(\mathbf{u}_0) + \mathsf{A} \cdot \delta\mathbf{u} + \text{higher order terms}. \end{align}

Since

\begin{align} \frac{\mathrm{d}\mathbf{u}_0}{\mathrm{d}t} = \mathbf{f}(\mathbf{u}_0) = 0, \end{align}

we have, to linear order,

\begin{align} \frac{\mathrm{d}\delta\mathbf{u}}{\mathrm{d}t} = \mathsf{A} \cdot \delta\mathbf{u}. \end{align}

This is now a system of linear ordinary differential equations. If this were a one-dimensional system, the solution would be an exponential \(\delta u \propto \mathrm{e}^{\lambda u}\). In the multidimensional case, the growth rate \(\lambda\) is replaced by a set of eigenvalues of \(\mathsf{A}\) that represent the growth rates in different directions, given by the eigenvectors of \(\mathsf{A}\)

2. Compute the eigenvalues, \(\lambda\) of \(\mathsf{A}\).

3. Determine the stability of the fixed point using the following rules.

If \(\mathrm{Re}(\lambda) < 0\) for all \(\lambda\), then the fixed point \(\mathbf{u}_0\) is linearly stable.
If \(\mathrm{Re}(\lambda) > 0\) for any \(\lambda\), then the fixed point \(\mathbf{u}_0\) is linearly unstable.
If \(\mathrm{Im}(\lambda) \ne 0\) for a linearly unstable fixed point, the trajectories spiral out, potentially leading to oscillatory dynamics.
If \(\mathrm{Re}(\lambda) = 0\) for one or more \(\lambda\), with the rest having \(\mathrm{Re}(\lambda) < 0\), then the fixed point \(\mathbf{u}_0\) lies at a bifurcation.

So, if we can assess the dynamics of the linearized system near the fixed point, we can get an idea what is happening with the full system.

Handy Taylor series of Hill functions

Because we model dynamics using Hill function so often, we will write the linearization of them here for you to memorize.

\begin{align} \frac{x^n}{1+x^n} &= \frac{x_0^n}{1+x_0^n} + \frac{n x_0^{n-1}}{(1+x_0^n)^2}\,\delta x + \text{higher order terms}, \\[1em] \frac{1}{1+x^n} &= \frac{1}{1+x_0^n} - \frac{n x_0^{n-1}}{(1+x_0^n)^2}\,\delta x + \text{higher order terms}. \end{align}

In the following analysis of the repressilator, we only need the second, repressing case.

Linear stability analysis for the repressilator

In Technical Appendix 9a, we wrote the dimensionless dynamical equations of the repressilator as

\begin{align} &\frac{\mathrm{d}x_1}{\mathrm{d}t} = \frac{\beta}{1 + x_3^n} - x_1, \\[1em] &\frac{\mathrm{d}x_2}{\mathrm{d}t} = \frac{\beta}{1 + x_1^n} - x_2, \\[1em] &\frac{\mathrm{d}x_3}{\mathrm{d}t} = \frac{\beta}{1 + x_2^n} - x_3, \end{align}

and showed that there exists a unique fixed point that satisfies

\begin{align} x_1 = x_2 = x_3 \equiv x_0 = \frac{\beta}{1 + x_0^n}, \end{align}

or

\begin{align} \beta = x_0(1+x_0^n). \end{align}

To perform linear stability analysis for the repressilator, we begin by writing the linearized system, performing a Taylor series expansion about the fixed point \(x_1 = x_2 = x_3 \equiv x_0\).

\begin{align} \frac{\mathrm{d}\delta x_1}{\mathrm{d}t} &\approx -\frac{\beta n x_0^{n-1}}{(1+x_0^n)^2}\,\delta x_3 - \delta x_1, \\[1em] \frac{\mathrm{d}\delta x_2}{\mathrm{d}t} &\approx -\frac{\beta n x_0^{n-1}}{(1+x_0^n)^2}\,\delta x_1 - \delta x_2, \\[1em] \frac{\mathrm{d}\delta x_3}{\mathrm{d}t} &\approx -\frac{\beta n x_0^{n-1}}{(1+x_0^n)^2}\,\delta x_2 - \delta x_3. \end{align}

Defining

\begin{align} a = \frac{\beta n x_0^{n-1}}{(1+x_0^n)^2}, \end{align}

we can write this in matrix form as

\begin{align} \frac{\mathrm{d}}{\mathrm{d}t}\begin{pmatrix} \delta x_1 \\ \delta x_2 \\ \delta x_3 \end{pmatrix} = \mathsf{A}\cdot\begin{pmatrix} \delta x_1 \\ \delta x_2 \\ \delta x_3 \end{pmatrix}, \end{align}

with

\begin{align} \mathsf{A} = -\begin{pmatrix} 1 & 0 & a \\ a & 1 & 0 \\ 0 & a & 1 \end{pmatrix}. \end{align}

To compute the eigenvalues of \(\mathsf{A}\), we compute the characteristic polynomial using cofactors,

\begin{align} (1+\lambda)(1+\lambda)^2 + a(a^2) = (1+\lambda)^3 + a^3 = 0. \end{align}

This is solved to give

\begin{align} \lambda = -1 + a \sqrt[3]{-1}. \end{align}

Recalling that there are three cube roots of \(-1\), we get our three eigenvalues.

\begin{align} &\lambda = -1 - a, \\[1em] &\lambda = -1 + \frac{a}{2}(1 + i\sqrt{3}),\\[1em] &\lambda = -1 + \frac{a}{2}(1 - i\sqrt{3}). \end{align}

The first eigenvalue is always real and negative. The second two have a positive real part if \(a > 2\);

\begin{align} a = \frac{\beta n x_0^{n-1}}{(1 + x_0^n)^2} > 2. \end{align}

We previously derived that the fixed point \(x_0\) satisfies

\begin{align} \beta = x_0(1+x_0^n), \end{align}

so

\begin{align} a = \frac{\beta n x_0^{n-1}}{(1 + x_0^n)^2} = \frac{n x_0^n}{1 + x_0^n}. \end{align}

So, \(a>2\) only if \(n > 2\), meaning that we must have ultrasensitivity for the fixed point to be unstable.

At the bifurcation,

\begin{align} a = \frac{n x_0^n}{1+x_0^n} = 2, \end{align}

so

\begin{align} x_0^n = \frac{2}{n-2}. \end{align}

Again using \(\beta = x_0(1+x_0^n)\), we can write

\begin{align} \beta = \frac{n}{2}\left(\frac{n}{2} - 1\right)^{-\frac{n+1}{n}} \end{align}

at the bifurcation. So, for \(n > 2\) and

\begin{align} \beta > \frac{n}{2}\left(\frac{n}{2} - 1\right)^{-\frac{n+1}{n}}, \end{align}

we have imaginary eigenvalues with positive real parts. This is therefore an oscillatory instability.

Linear stability diagram

When analyzing stability of dynamical systems, it is useful to make a linear stability diagram, which is a map of parameter space highlighting stable and unstable regions. We know the bifurcation line is

\begin{align} \beta = \frac{n}{2}\left(\frac{n}{2} - 1\right)^{-\frac{n+1}{n}}. \end{align}

We can plot this line and delineate the regions of stability and instability.

[2]:

# Get bifurcation line
n = np.linspace(2.01, 5, 200)
beta = n / 2 * (n / 2 - 1) ** (-(1 + 1 / n))

# Build the plot
p = bokeh.plotting.figure(
    height=300,
    width=400,
    x_axis_label="n",
    y_axis_label="β",
    y_axis_type="log",
    x_range=[2, 5],
    y_range=[1, 2000],
)
p.patch(
    np.append(n, n[-1]), np.append(beta, beta[0]), color="lightgray", alpha=0.7
)
p.line(n, beta, line_width=4, color="black")
p.text(x=2.1, y=2, text=["stable"])
p.text(x=2.5, y=100, text=["unstable (limit cycle oscillations)"])

bokeh.io.show(p)

References

Strogatz, S. H., Nonlinear Dynamics and Chaos With Applications to Physics, Biology, Chemistry, and Engineering, 2nd Ed., CRC Press, 2015. (link)

Computing environment

[3]:

%load_ext watermark
%watermark -v -p numpy,bokeh,jupyterlab

Python implementation: CPython
Python version       : 3.10.10
IPython version      : 8.12.0

numpy     : 1.23.5
bokeh     : 3.1.0
jupyterlab: 3.5.3