Delay oscillators¶
(c) 2018 Justin Bois. With the exception of pasted graphics, where the source is noted, this work is licensed under a Creative Commons Attribution License CC-BY 4.0. All code contained herein is licensed under an MIT license.
This document was prepared at Caltech with financial support from the Donna and Benjamin M. Rosen Bioengineering Center.
This tutorial was generated from an Jupyter notebook. You can download the notebook here.
As usual, we will first import the modules we need.
# Numerical workhorses
import numpy as np
import scipy.optimize
# Numba for just-in-time compilations (speed)!
import numba
# Plotting modules
import bokeh.io
import bokeh.plotting
bokeh.io.output_notebook()
Simple delay oscillator¶
We start with a simple delay oscillator. The dimensionless delay differential equation describing its dynamics is
\begin{align} \frac{\mathrm{d}x(t)}{\mathrm{d}t} = \frac{\beta}{1 - (x(t-\tau))^n} - x(t). \end{align}
We solve the delayed differential equation (DDE) using simple Euler integration. There are more sophisticated ways to do this; see, e.g., Bellen and Zennaro, Numerical Methods for Delay Differential Equations, Oxford University Press, 2003.
We will use Numba to perform just-in-time compilation, which will greatly speed the calculation. (When I tested it, I got about a 180 fold speed up.) This is not necessary, unless you really care about speed, and the @numba decorators can be removed and the code will still function as expected.
@numba.jit(nopython=True)
def f(x, n):
"""
Hill-like repression
"""
return 1 / (1 + x**n)
@numba.jit(nopython=True)
def dx_dt(x_past, x_current, beta, tau, n):
"""
Dynamics of x, dependent on a past and current value of x.
"""
return beta * f(x_past, n) - x_current
@numba.jit(nopython=True)
def solve_delay_oscillator(x_0, beta, tau, n, dt=0.001, t_stop=30.0):
"""
Use Euler integration to solve ODEs
"""
# Number of indices to go back in time
i_time = int(tau / dt)
# Time points
t = np.linspace(-tau, t_stop, int((t_stop + tau) / dt))
# Initialize output array
x = x_0 * np.ones_like(t)
# Do Euler stepping
for i in range(i_time, len(t) - 1):
x[i+1] = x[i] + dt * dx_dt(x[i - i_time], x[i], beta, tau, n)
return t, x
Now that we have the machinery in place, we can go ahead and solve.
# Specify parameters
x_0 = 0
beta = 4
n = 2
tau = 10
# Perform the solution
t, x = solve_delay_oscillator(x_0, beta, tau, n, dt=0.001, t_stop=200)
# Plot the result
p = bokeh.plotting.figure(plot_width=600,
plot_height=300,
x_axis_label='time',
y_axis_label='x')
p.line(t, x, line_width=2)
bokeh.io.show(p)
Linear stability diagram of delay oscillator¶
In lecture, we found that the unique fixed point, $x_0$ is given by
\begin{align} \beta = x_0(1+x_0^n). \end{align}
The bifurcation line between a stable fixed point and an oscillatory instability is given by
\begin{align} \tau = \frac{\pi - \cos^{-1}\left(a_\tau^{-1}\right)}{\sqrt{a_\tau^2 - 1}}, \end{align}
where
\begin{align} a_\tau = \frac{nx_0^n}{1+x_0^n}. \end{align}
We will compute and plot the bifurcation line in the $\tau$-$\beta$ plane for various values of ultrasensitivity $n$.
def fp_rhs(x, beta, n):
return beta / (1 + x**n)
def bifurcation_line(beta, n):
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(plot_width=500,
plot_height=300,
x_axis_label='β',
y_axis_label='Ï„',
x_axis_type='log',
y_axis_type='log',
x_range=bokeh.models.Range1d(0.8, 100),
y_range=bokeh.models.Range1d(0.1, 25))
# Compute bifurcation line and plot
colors = bokeh.palettes.Blues8[::-1]
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='seashell',
alpha=0.7,
level='underlay')
p.line(beta, tau, color=colors[i+2], legend=f'n={n}', line_width=2)
p.text(x=10, y=1, text=['unstable'])
p.text(x=1, y=0.2, text=['stable'])
bokeh.io.show(p)
Delay due to mRNA dynamics¶
We can make a linear stability diagram for our analysis of delay due to slow mRNA dynamics. As derived in lecture, the fixed point is given by
\begin{align} m_0 &= r_0 = \gamma p_0,\\[1em] \frac{\beta}{\gamma} &= p_0(1+p_0^n). \end{align}
The linear stability matrix is
\begin{align} \mathsf{A} = \begin{pmatrix} -1 & 0 & -a \\ 1 & -1 & 0 \\ 0 & 1 & -\gamma \end{pmatrix}, \end{align}
where
\begin{align} a = \frac{\beta n p_0^{n-1}}{(1+p_0^n)^2} = \frac{\gamma n p_0^n}{1+p_0^n}. \end{align}
We will code up functions to compute the fixed point and to determine whether or not there is a positive real part to the complex conjugate eigenvalues.
def delay_fixed_point(beta, gamma, n):
"""
Computes fixed point of delayed oscillator. Uses np.roots because
it has better convergence properties than scipy.optimize.fixed_point
for integer n.
"""
coeffs = np.zeros(n+1)
coeffs[0] = 1
coeffs[-2] = 1
coeffs[-1] = -beta / gamma
p_0 = np.roots(coeffs)
return p_0[(p_0.imag == 0) & (p_0 > 0)].real[0]
def lin_stab(beta, gamma, n):
"""
Returns True is delay system is linearly stable and False otherwise.
"""
p_0 = delay_fixed_point(beta, gamma, n)
# Linear stability matrix
a = gamma * n * p_0**n / (1 + p_0**n)
L = np.array([[-1, 0, -a],
[1, -1, 0],
[0, 1, -gamma]])
# Compute eigenvalues
lam = np.linalg.eigvals(L)
return (lam.real < 0).all()
Now, we'll compute the stability diagram.
# Ultrasensitivity
n = 10
# Compute the linear stability diagram
beta_vals = np.logspace(0, 4, 400)
gamma_vals = np.logspace(-1, 2, 400)
stab = np.empty((len(beta_vals), len(gamma_vals)))
for i, beta in enumerate(beta_vals):
for j, gamma in enumerate(gamma_vals):
stab[i,j] = lin_stab(beta, gamma, n)
# Plot the result (shaded region has oscillatory instability)
p = bokeh.plotting.figure(plot_height=300,
plot_width=350,
x_axis_label='γ',
y_axis_label='β',
x_axis_type='log',
y_axis_type='log',
x_range=bokeh.models.Range1d(0.1, 100),
y_range=bokeh.models.Range1d(1, 10000))
p.image([stab], x=0.1, y=1, dw=100, dh=10000, palette='Greys3', alpha=0.5)
p.text(x=0.45, y=100, text=['unstable'])
p.text(x=10, y=2, text=['stable'])
bokeh.io.show(p)
Now that we know the region in parameter space to get oscillations, we can solve the dynamical equations and look at the results.
def rhs(y, t, beta, gamma, n):
"""
Right-hand-side for dimensionless delay-due-to-mRNA processing dynamics.
"""
m, r, p = y
return np.array([beta / (1 + p**n) - m,
m - r,
r - gamma * p])
# Set parameters
t = np.linspace(0, 50, 2000)
y0 = np.array([0, 0, 0])
beta = 10
gamma = 1
n = 10
# Perform solution
y = scipy.integrate.odeint(rhs, y0, t, args=(beta, gamma, n))
m, r, p = y.transpose()
# Plot the result
plot = bokeh.plotting.figure(plot_height=250,
plot_width=600,
x_axis_label='time')
plot.line(t, m, line_width=2, legend='m')
plot.line(t, r, line_width=2, color='orange', legend='r')
plot.line(t, p, line_width=2, color='mediumorchid', legend='p')
bokeh.io.show(plot)
We can also look at the trajectories in phase space.
# Compute fixed point
p_0 = delay_fixed_point(beta, gamma, n)
r_0 = gamma * p_0
m_0 = r_0
# Make plots
plots = []
for x, y, x0, y0, x_label, y_label in zip([m, m, r],
[r, p, p],
[m_0, m_0, r_0],
[r_0, p_0, p_0],
['m', 'm', 'r'],
['r', 'p', 'p']):
plot = bokeh.plotting.figure(plot_height=200,
plot_width=250,
x_axis_label=x_label,
y_axis_label=y_label)
plot.circle(x0, y0, color='tomato', size=10)
plot.line(x, y, line_width=2)
plots.append(plot)
bokeh.io.show(bokeh.layouts.gridplot(plots, ncols=2))