10c. Numerical solution 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.integrate
import scipy.interpolate
import scipy.optimize
import bokeh.io
import bokeh.plotting
bokeh.io.output_notebook()
Numerical solution of DDEs
In general, we may write a system of DDEs as
\begin{align} \frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t} = &= \mathbf{f}(t, \mathbf{u}(t), \mathbf{u}(t-\tau_1), \mathbf{u}(t-\tau_2), \ldots, \tau_n), \end{align}
where \(\tau_1, \tau_2, \ldots, \tau_n\) are the delays with \(\tau_1 < \tau_2 < \ldots < \tau_n\). We define the longest delay to be \(\tau_n = \tau\) for notational convenience. If we wish to have a solution for time \(t \ge t_0\), we need to specify initial conditions that go sufficiently far back in time,
\begin{align} y(t) = y_0(t) \text{ for }t_0-\tau < t \le t_0. \end{align}
The approach we will take is the method of steps. Given \(y_0(t)\), we integrate the differential equations from time \(t_0\) to time \(t_0 + \tau\). We store that solution and then use it as the past values when we integrate again from time \(t_0+\tau\) to time \(t_0 + 2\tau\). We continue this process until we integrate over the desired time span.
Because we need to access \(y(t-\tau_i)\) at time \(t\), and \(t\) can be arbitrary, we often need an interpolation of the solution from the previous time steps. We will use SciPy’s functional interface for B-spline interpolation, using the scipy.interpolate.splrep() and scipy.interpolate.splev() functions.
We write the function below to solve DDEs with similar call signature as scipy.integrate.odeint().
ddeint(func, y0, t, tau, args=(), y0_args=(), n_time_points_per_step=200)
The input function has a different call signature than odeint() because it needs to allow for a function that computes the values of y in the past. So, the call signature of the function func giving the right hand side of the system of DDEs is
func(y, t, y_past, *args)
Furthermore, instead of passing in a single initial condition, we need to pass in a function that computes \(y_0(t)\), as defined for initial conditions that go back in time. The function has call signature
y0(t, *y0_args)
We also need to specify tau, the longest time delay in the system. Finally, in addition to the additional arguments to func and to y0, we specify how many time points we wish to use for each step to build our interpolants. Note that even though we specify tau in the arguments of ddeint(), we still need to explicitly pass all time delays as extra args to func. This is because we could in general have more than one delay, and only the longest is pertinent for the
integration.
[2]:
def ddeint(func, y0, t, tau, args=(), y0_args=(), n_time_points_per_step=None):
"""Integrate a system of delay differential equations.
"""
if np.isscalar(tau):
tau = np.array([tau])
else:
tau = np.array(tau)
if (tau <= 0).any():
raise RuntimeError("All tau's must be greater than zero.")
tau_short = np.min(tau)
tau_long = np.max(tau)
if n_time_points_per_step is None:
n_time_points_per_step = max(
int(1 + len(t) / (t.max() - t.min()) * tau_long), 20
)
t0 = t[0]
# Past function for the first step
y_past = lambda time_point: y0(time_point, *y0_args)
# Integrate first step
t_step = np.linspace(t0, t0 + tau_short, n_time_points_per_step)
y = scipy.integrate.odeint(func, y_past(t0), t_step, args=(y_past,) + args)
# Store result from integration
y_dense = y.copy()
t_dense = t_step.copy()
# Get dimension of problem for convenience
n = y.shape[1]
# Integrate subsequent steps
j = 1
while t_step[-1] < t[-1]:
t_start = max(t0, t_step[-1] - tau_long)
i = np.searchsorted(t_dense, t_start, side="left")
t_interp = t_dense[i:]
y_interp = y_dense[i:, :]
# Make B-spline
tck = [scipy.interpolate.splrep(t_interp, y_interp[:, i]) for i in range(n)]
# Interpolant of y from previous step
y_past = (
lambda time_point: np.array(
[scipy.interpolate.splev(time_point, tck[i]) for i in range(n)]
)
if time_point > t0
else y0(time_point, *y0_args)
)
# Integrate this step
t_step = np.linspace(
t0 + j * tau_short, t0 + (j + 1) * tau_short, n_time_points_per_step
)
y = scipy.integrate.odeint(func, y[-1, :], t_step, args=(y_past,) + args)
# Store the result
y_dense = np.append(y_dense, y[1:, :], axis=0)
t_dense = np.append(t_dense, t_step[1:])
j += 1
# Interpolate solution for returning
y_return = np.empty((len(t), n))
for i in range(n):
tck = scipy.interpolate.splrep(t_dense, y_dense[:, i])
y_return[:, i] = scipy.interpolate.splev(t, tck)
return y_return
Now that we have a function to solve DDEs, let’s solve the simple delay oscillator.
[3]:
# Right hand side of simple delay oscillator
def delay_rhs(x, t, x_past, beta, n, tau):
return beta / (1 + x_past(t - tau) ** n) - x
# Initially, we just have no gene product
def x0(t):
return 0
# Specify parameters (All dimensionless)
beta = 4
n = 2
tau = 10
# Time points we want
t = np.linspace(0, 200, 400)
# Perform the integration
x = ddeint(delay_rhs, x0, t, tau, args=(beta, n, tau))
# Plot the result
p = bokeh.plotting.figure(
frame_width=450,
frame_height=200,
x_axis_label="time",
y_axis_label="x",
x_range=[t.min(), t.max()],
)
p.line(t, x[:, 0], line_width=2)
bokeh.io.show(p)
The ddeint() function is available as biocircuits.ddeint() in the biocircuits package for future use.
Computing environment
[4]:
%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