(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()
```

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 |

K_{m} |
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.

In [6]:

```
%load_ext watermark
%watermark -v -p numpy,numba,bokeh,jupyterlab
```