Delay oscillators and DDEs

(c) 2020 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.


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()
Loading BokehJS ...

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