10d. Linear stability analysis of delay differential equations
[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 scipy.optimize
import bokeh.io
import bokeh.plotting
bokeh.io.output_notebook()
In this technical appendix, we describe a general procedure for linear stability analysis of delay differential equations, and then apply the technique to a simple example of a single gene the features delayed autorepression.
General approach for linear stability analysis of a DDEs
Consider a system of delay differential equations,
\begin{align} \frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t} &= \mathbf{f}(\mathbf{u}(t), \mathbf{u}(t-\tau)), \end{align}
where for simplicity we consider a single delay \(\tau\). The procedure generalizes to multiple delays. We have written that the time derivative explicitly depends on the value of \(\mathbf{u}\) at the present time and also at some time \(\tau\) time units ago in the past. A steady state of the delay differential equations, \(\mathbf{u}_0\) satisfies \(\mathbf{u}_0(t) = \mathbf{u}_0(t-\tau) \equiv \mathbf{u}_0\) by the definition of a steady state. So, the steady state \(\mathbf{u}_0\) satisfies \(\mathbf{f}(\mathbf{u_0}, \mathbf{u}_0) = 0\). We will perform a linear stability analysis about this steady state.
To linearize, we have to remember that \(\mathbf{f}\) is a function of two sets of variables, present and past \(\mathbf{u}\), \(\mathbf{u}(t)\) and \(\mathbf{u}(t-\tau)\), respectively. We write a Taylor series expansion with respect to both of these variables to linearize. (For multiple delays, we expand in terms of \(\mathbf{u}(t)\), \(\mathbf{u}(t-\tau_1)\), \(\mathbf{u}(t-\tau_2)\), etc.)
\begin{align} \mathbf{f}(\mathbf{u}(t), \mathbf{u}(t-\tau)) \approx \mathbf{f}(\mathbf{u}_0, \mathbf{u}_0) + \mathsf{A}\cdot(\mathbf{u}(t) - \mathbf{u}_0) + \mathsf{A}_{\tau}\cdot(\mathbf{u}(t-\tau) - \mathbf{u}_0), \end{align}
where
\begin{align} \mathsf{A} = \begin{pmatrix} \frac{\partial f_1}{\partial u_1(t)} & \frac{\partial f_1}{\partial u_2(t)} & \cdots \\[0.5em] \frac{\partial f_2}{\partial u_1(t)} & \frac{\partial f_2}{\partial u_2(t)} & \cdots \\ \vdots & \vdots & \ddots \end{pmatrix} \end{align}
and
\begin{align} \mathsf{A}_{\tau} = \begin{pmatrix} \frac{\partial f_1}{\partial u_1(t-\tau)} & \frac{\partial f_1}{\partial u_2(t-\tau)} & \cdots \\[0.5em] \frac{\partial f_2}{\partial u_1(t-\tau)} & \frac{\partial f_2}{\partial u_2(t-\tau)} & \cdots \\ \vdots & \vdots & \ddots \end{pmatrix}, \end{align}
where all derivatives in both matrices are evaluated at \(\mathbf{u}(t) = \mathbf{u}(t-\tau) = \mathbf{u}_0\). Defining \(\delta\mathbf{u}(t) = \mathbf{u}(t) - \mathbf{u}_0\) and \(\delta\mathbf{u}(t-\tau) = \mathbf{u}(t-\tau) - \mathbf{u}_0\), our linearized differential equations are
\begin{align} \frac{\mathrm{d}\delta\mathbf{u}}{\mathrm{d}t} = \mathsf{A}\cdot \delta\mathbf{u}(t) + \mathsf{A}_{\tau}\cdot \delta\mathbf{u}(t-\tau). \end{align}
To solve this linear system, we insert the ansatz \(\delta\mathbf{u} = \mathbf{w}\,\mathrm{e}^{\lambda t}\), giving
\begin{align} \mathbf{w}\,\lambda \mathrm{e}^{\lambda t} = \mathsf{A}\cdot\mathbf{w}\,\mathrm{e}^{\lambda t} + \mathsf{A}_{\tau}\cdot \mathbf{w}\,\mathrm{e}^{\lambda (t-\tau)}, \end{align}
or, upon dividing both sides by \(\mathrm{e}^{\lambda t}\),
\begin{align} \lambda \mathbf{w} = (\mathsf{A} + \mathrm{e}^{-\lambda \tau} \mathsf{A}_{\tau})\cdot \mathbf{w}. \end{align}
So, \(\lambda\) is an eigenvalue and \(\mathbf{w}\) is the corresponding eigenvector for the matrix \((\mathsf{A} + \mathrm{e}^{-\lambda \tau} \mathsf{A}_{\tau})\). In general, there are infinitely many values of \(\lambda\) that solve the above equation. Nonetheless, if we can show that the real parts of all possible values of \(\lambda\) are less than zero, then the fixed point is stable. Otherwise, if the real parts of any values of \(\lambda\) are positive, the fixed point is locally unstable. If \(\lambda\) has an imaginary part for the eigenvalues with a positive real part, the instability is oscillatory.
Linear stability analysis of delayed autorepression
As a specific example, we consider now the simple case of an autorepressor with delay.
\begin{align} \frac{\mathrm{d}x(t)}{\mathrm{d}t} = \frac{\beta}{1 + (x(t-\tau))^n} - x(t). \end{align}
The steady state \(x_0\) is given by setting the time derivative to zero,
\begin{align} \beta = x_0(1+x_0^n). \end{align}
The \(x_0\) that satisfies this equality is unique because the right hand is monotonically increasing from zero, so it only crosses a value of \(\beta\) once. So, we will consider the stability of this unique steady state.
We linearize the right hand side of the ODE around the fixed point. The matrices \(\mathsf{A}\) and \(\mathsf{A}_{\tau}\) are scalars in this case because we have a single species. Remember that the eigenvalue of a scalar is just the scalar itself. So, we have
\begin{align} &\mathsf{A} = -1 \\[1em] &\mathsf{A}_{\tau} = -\frac{\beta n x_0^{n-1}}{(1+x_0^n)^2} = -\frac{nx_0^n}{1+x_0^n} \equiv -a_{\tau}, \end{align}
where we have defined \(a_{\tau}\) (which must be positive) for convenience. Thus, our characteristic equation is
\begin{align} \lambda = -1 - a_{\tau}\mathrm{e}^{-\lambda \tau}. \end{align}
There are infinitely many \(\lambda\) that satisfy this equation in general. Nonetheless, we can still make progress. We write \(\lambda = a + ib\), where \(a\) is the real part at \(b\) is the imaginary part. Then, the characteristic equation is
\begin{align} a + ib = -a_{\tau} \mathrm{e}^{-a\tau} \mathrm{e}^{-ib\tau} - 1 = -a_{\tau}\mathrm{e}^{-a\tau}(\cos b\tau - i\sin b\tau) - 1. \end{align}
This can be written as two equations by equating the real and imaginary parts of both sides of the equation.
\begin{align} a &= -(1+a_{\tau} \mathrm{e}^{-a\tau} \cos b\tau ) \\[1em] b &= a_{\tau} \mathrm{e}^{-a\tau} \sin b\tau. \end{align}
Right away, we can see that if \(a\) is positive, we must have \(a_{\tau} > 1\), since \(\left|\mathrm{e}^{-a\tau}\cos b\tau\right| < 1\). Recall our expression for \(a_{\tau}\),
\begin{align} a_{\tau} = \frac{nx_0^n}{1+x_0^n}. \end{align}
Because \(x_0^n/(1+x_0^n) < 1\), we must have \(n>1\) in order to have the eigenvalue have a positive real part and therefore have an instability. So, ultrasensitivity is a requirement for a simple delay oscillator.
To investigate when we get an oscillatory instability, we will compute the values of \(a_{\tau}\) and \(\tau\) that lie on the bifurcation line between a stable steady state and an oscillatory instability, a so-called Hopf bifurcation. To do this, we solve for the line with \(a = 0\) and \(b \ne 0\). In this case, the characteristic equation gives
\begin{align} -a_{\tau} \cos bt &= 1 \\[1em] a_{\tau} \sin bt &= b. \end{align}
Squaring both equations and adding gives
\begin{align} a_{\tau}^2 = 1 + b^2, \end{align}
Thus,
\begin{align} b = \sqrt{a_{\tau}^2 - 1}, \end{align}
which can only hold for \(a_{\tau} > 1\), which we already found was a requirement for an oscillatory instability. To find \(\tau\), we have
\begin{align} &-a_{\tau} \cos b\tau = -a_{\tau} \cos \sqrt{a_{\tau}^2 - 1}\,t = 1 \\[1em] &\Rightarrow\; \tau = \frac{\pi - \cos^{-1}\left(a_{\tau}^{-1}\right)}{\sqrt{a_{\tau}^2 - 1}}. \end{align}
The region of stability is below the bifurcation line in the \(\tau\)-\(\beta\) plane, since a smaller \(\tau\) serves to decrease the real part of \(\lambda\).
Now that we have worked it out analytically, we can compute and plot the linear stability diagram in the \(\tau\)-\(\beta\) plane.
[2]:
def fp_rhs(x, beta, n):
""" Right-hand side of the fixed point equation."""
return beta / (1 + x ** n)
def bifurcation_line(beta, n):
"""Solve for bifurcation line by numerically finding fixed points
and plugging in to expression for bifurcation."""
x_0 = scipy.optimize.fixed_point(fp_rhs, 1, args=(beta, n))
a_tau = n * x_0 ** n / (1 + x_0 ** n)
if a_tau <= 1:
return np.nan
return np.arccos(-1 / a_tau) / np.sqrt(a_tau ** 2 - 1)
# Set up values for plotting
n_vals = (1.2, 2, 5, 10)
beta = np.logspace(-0.2, 2, 2000)
# Set up plot
p = bokeh.plotting.figure(
frame_width=300,
frame_height=300,
x_axis_label="β",
y_axis_label="τ",
x_axis_type="log",
y_axis_type="log",
x_range=[0.8, 100],
y_range=[0.1, 25],
toolbar_location='above',
)
# Compute bifurcation lines and plot
colors = bokeh.palettes.Blues8[::-1]
lines = []
legend_items = []
for i, n in enumerate(n_vals):
tau = np.array([bifurcation_line(b, n) for b in beta])
p.patch(
np.append(beta, beta[-1]),
np.append(tau, np.nanmax(tau)),
color="gray",
alpha=0.2,
level="underlay",
)
lines.append(p.line(beta, tau, color=colors[i + 2], line_width=2))
legend_items.append((f"n = {n}", [lines[-1]]))
p.text(x=10, y=1, text=["unstable"])
p.text(x=1, y=0.2, text=["stable"])
# Instantiate legend
legend = bokeh.models.Legend(items=legend_items, click_policy='hide')
# Add the legend to the right of the plot
p.add_layout(legend, 'right')
bokeh.io.show(p)
The unstable region is above the respective curves for the various values of the Hill coefficient \(n\).
Computing environment
[3]:
%load_ext watermark
%watermark -v -p numpy,scipy,bokeh,jupyterlab
Python implementation: CPython
Python version : 3.11.3
IPython version : 8.12.0
numpy : 1.24.3
scipy : 1.10.1
bokeh : 3.2.0
jupyterlab: 3.6.3