# Delay oscillators and DDEs¶

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.

In [1]:
# Numerical workhorses
import numpy as np

# Numba for just-in-time compilations (speed)!
import numba

# Plotting modules
import bokeh.io
import bokeh.plotting

bokeh.io.output_notebook()


## The Swaminathan delay oscillator¶

The dynamical equations for the delay oscillator are

\begin{align} &\frac{\mathrm{d}x}{\mathrm{d}t} = \frac{\beta R}{1 + (x(t-\tau)/k)^n} - \gamma_{0} \,x - \gamma_{1}\, \frac{x}{K_m + x + y}, \\[1em] &\frac{\mathrm{d}y}{\mathrm{d}t} = \frac{\beta F}{1 + (x(t-\tau)/k)^n} - \gamma_{0} \,y - \gamma_{1}\, \frac{y}{K_m + x + y}, \end{align}

The parameter values suggested by Anandh Swaminathan are as follows.

Paramter Value Units
β 300 per minute
γ₀ 0.02 per minute
γ₁ 90 molecules per minute
R 12 copies
F 60 copies
k (treRL) 0.03 same as x
k (LacI) 0.5 same as x
Km 100 same as x
n 2
τ 5 minutes

We solve the delay 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.) The basic idea is that we take small steps forward in time along the derivative in each variable. Since this is a delay differential equation, the derivative at time $t$ depends on that at time $t-\tau$. We therefore look back into the past to compute the current derivative.

Throughout, we will use Numba to perform just-in-time compilation, which will greatly speed the calculation. (When I tested it, I got about a 20 fold speed up.) This is not necessary, unless you really care about speed (which you will if you care about interactive plotting, below), and the @numba decorators can be removed and the code will still function as expected.

First, it is convenient to code up a Hill function.

In [2]:
@numba.njit
def f(x, k, n):
"""
Hill-like repression
"""
return 1 / (1 + (x/k)**n)


Next, we write a function to compute the derivatives, dx/dt and dy/dt. For convenience, we pass in the values of x and y together as an array, which we will call u. We need to give both u at time t-τ and at time t. Finally, we pass in the parameters.

In [3]:
@numba.njit
def du_dt(u_past, u_current, beta, gamma_0, gamma_1, R, F, k, K_m, n):
"""
Dynamics of u, dependent on a past and current value of u.
"""
# Unpack inputs
x_past, y_past = u_past
x, y = u_current

# Compute derivatives
dx_dt = (
beta * R * f(x_past, k, n) - gamma_0 * x - gamma_1 * x / (K_m + x + y)
)

dy_dt = (
beta * F * f(x_past, k, n) - gamma_0 * y - gamma_1 * y / (K_m + x + y)
)

return np.array([dx_dt, dy_dt])


Finally, we write a function to perform Euler integration. We pass the parameters as keyword arguments, allowing us to set defaults. The default values hold for both treRL and LacI, except that for LacI, k=0.5

In [4]:
@numba.njit
def solve_delay_oscillator(
u_0,
beta=300.0,
gamma_0=0.02,
gamma_1=90.0,
R=12,
F=60,
k=0.03,
K_m=100.0,
n=2,
tau=5,
dt=0.001,
t_stop=300,
):
"""
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
u = np.empty((len(t), len(u_0)))
for i in range(i_time + 1):
u[i, :] = u_0

# Do Euler stepping
for i in range(i_time, len(t) - 1):
u_deriv = du_dt(
u[i - i_time], u[i], beta, gamma_0, gamma_1, R, F, k, K_m, n
)

u[i + 1] = u[i] + dt * u_deriv

return t, u


Now that we have the machinery in place, we can go ahead and solve.

In [5]:
# Set up initial conditions (doesn't really matter)
u_0 = np.array([1.0, 1.0])

# Perform the solution
t, u = solve_delay_oscillator(u_0)

# Plot the result (normalized by highest value)
p = bokeh.plotting.figure(
frame_width=550,
frame_height=250,
x_axis_label="time (min.)",
y_axis_label="fluoresence (a.u.)",
)

# Only plot every 300th point, since we don't need too many points for display
p.line(t[::300], u[::300, 1] / u[:, 1].max(), line_width=2, line_join="bevel")
bokeh.io.show(p)


Indeed, we see sustained oscillations.

## Computing environment¶

In [6]:
%load_ext watermark
%watermark -v -p numpy,numba,bokeh,jupyterlab

CPython 3.7.7
IPython 7.13.0

numpy 1.18.1
numba 0.49.0
bokeh 2.0.2
jupyterlab 1.2.6