2b. Numerical solutions to ODEs with SciPy

---

# 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 colorcet

import bokeh.io
import bokeh.plotting

bokeh.io.output_notebook()

Loading BokehJS ...

In this chapter, we numerically solve the ODE for a negative autoregulatory circuit,

\begin{align} \frac{\mathrm{d}x}{\mathrm{d}t} = \beta_0\,\frac{(s/k_s)^{n_s}}{1 + (s/k_s)^{n_s}}\, \frac{1}{(1 + (x/k)^n)} - \gamma x. \end{align}

In this technical appendix, we solve this ODE numerically to generate the plots shown in the chapter.

The scipy.integrate module

The SciPy Library is a Python library for scientific computing. It contains many modules, including scipy.stats, scipy.special, and scipy.optimize, which respectively include functions to perform statistical calculations, “special functions,” and optimization routines, among many others. We will use the scipy.integrate module to integrate systems of ODEs.

There are three main APIs for solving real-valued initial value problems in the module. They are solve_ivp(), ode(), and odeint(). According to the SciPy developers, solve_ivp() is the preferred method, with the others labeled as having an “old” API. The solve_ivp() function has the flexibility of allowing choice of multiple numerical algorithms for solving ODEs. However, for the kinds of problems we encounter in this class, we favor the generic LSODA algorithm developed by Linda Petzold and Alan Hindmarsh that handles both stiff and non-stiff problems with variable time stepping. This is also the only solver offered in the odeint() function. If we compare the two solvers, solve_ivp() and odeint(), the former has a large overhead, which can lead to performance issues for small problems (for large problems, the overhead is negligible). Since most of our problems are small and we need to solve them rapidly for interactive graphics, we will use odeint(), which performs better for these problems and through a distinct, but still intuitive, API.

The basic call signature for odeint() is

scipy.integrate.odeint(func, y0, t, args=())

There are many other keyword arguments to set algorithmic parameters, but we will generally not need them (and you can read about them in the documentation). Importantly, func is a vector-valued function with call signature func(y, t, *args) that specifies the right hand side of the system of ODEs to be solved. t is a scalar time point and y is a one-dimensional array (though multidimensional arrays are possible). y0 is an array with the initial conditions.

As is often the case, use of this function is best seen by example, and we will now apply it to the negative autoregulation circuit.

Solving for a constant input

We first seek to solve the ODE where we have a sudden step in \(s\) at time \(t = 0\). Since we are only simulating our system starting at \(t=0\), we can treat \(s\) as a single, constant parameter, rather than a time-varying variable. So, we need seven parameters for the right hand side of our ODEs: \(\beta_0\), \(\gamma\), \(k\), \(n\), \(k_s\), \(n_s\), and \(s\).

We now define the function for the right hand side of the ODEs.

def neg_auto_rhs(x, t, beta0, gamma, k, n, ks, ns, s):
    """
    Right hand side for negative autoregulation motif with s dependence.
    Return dx/dt.
    """
    # Compute dx/dt
    return (
        beta0 * (s / ks) ** ns / (1 + (s / ks) ** ns) / (1 + (x / k) ** n) - gamma * x
    )

We can now define the initial conditions, our parameters (taking \(n_s = 10\)), and the time points we want and use scipy.integrate.odeint() to solve.

# Time points we want for the solution
t = np.linspace(0, 10, 200)

# Initial condition
x0 = 0.0

# Parameters
beta0 = 100
gamma = 1.0
k = 1.0
n = 1.0
s = 100.0
ns = 10.0
ks = 0.1

# Package parameters into a tuple
args = (beta0, gamma, k, n, ks, ns, s)

# Integrate ODES
x = scipy.integrate.odeint(neg_auto_rhs, x0, t, args=args)

That’s it! The integration is done. Let’s take a quick look at the output.

x.shape

(200, 1)

We see that the output of odeint() has each column corresponding to a given species and each row to a given time point. It is often convenient to transpose the output so that the each species is more easily index-able.

# Extract time course for first (in this case only) species
x = x.transpose()[0]

When we plot the result, we would like to compare it both to the limiting result of Rosenfeld et al.,

\begin{align} x(t) \approx x_\mathrm{ss} \sqrt{1-\mathrm{e}^{-2 \gamma t}}, \end{align}

and also to the unregulated case,

\begin{align} x(t) = \frac{\beta_0}{\gamma}\left(1 - \mathrm{e}^{-\gamma t}\right), \end{align}

so we need to code up those solutions as well for plotting.

# Unregulated solution
x_unreg = beta0 / gamma * (1 - np.exp(-gamma * t))

# Limiting analytical solution
x_limiting = x[-1] * np.sqrt(1 - np.exp(-2 * gamma * t))

Now let’s plot the results. As a reminder, you can read about making plots with Bokeh in Appendix B.

# Set up color palette for this notebook
colors = colorcet.b_glasbey_category10

# Set up figure
p = bokeh.plotting.figure(
    frame_width=325,
    frame_height=250,
    x_axis_label="time",
    y_axis_label="x",
    x_range=[np.min(t), np.max(t)],
)

cds = bokeh.models.ColumnDataSource(
    dict(t=t, x=x, x_limiting=x_limiting, x_unreg=x_unreg)
)

# Populate glyphs
p.line(source=cds, x="t", y="x", line_width=2, color=colors[0], legend_label="numerical solution")
p.line(source=cds, x="t", y="x_limiting", line_width=2, color=colors[1], legend_label="Rosenfeld limiting solution")
p.line(source=cds, x="t", y="x_unreg", line_width=2, color=colors[2], legend_label="unregulated")

# Aesthetic tweaks
p.legend.location = "center_right"
p.legend.click_policy = "hide"
p.title.text = "Constant-input dynamics"

bokeh.io.show(p)

Solving for a time-varying input

We consider now a pulsatile input signal,

\begin{align} s(t) = \exp\left[-\frac{4(t-t_0)^2}{\tau^2}\right]. \end{align}

In order to incorporate these dynamics into our model, we can write a function for the pulse.

def s_pulse(t, t_0, tau):
    """
    Returns s value for a pulse centered at t_0 with duration tau.
    """
    # Return 0 is tau is zero, otherwise Gaussian
    return 0 if tau == 0 else np.exp(-4 * (t - t_0) ** 2 / tau ** 2)

If we want to solve the ODEs for a varying input, we need to have a way to pass a function defining the variation as a parameter. Fortunately, we can pass functions as arguments in Python. So, we write a new function for the right-hand-side of our ODE that takes s_fun, the function describing \(s(t)\) as an argument, as well as s_args, the set of parameters passed into s_fun.

def neg_auto_rhs_s_fun(x, t, beta0, gamma, k, n, ks, ns, s_fun, s_args):
    """
    Right hand side for negative autoregulation function, with s variable.
    Returns dx/dt.

    s_fun is a function of the form s_fun(t, *s_args), so s_args is a tuple
    containing the arguments to pass to s_fun.
    """
    # Compute s
    s = s_fun(t, *s_args)

    # Correct for x possibly being numerically negative as odeint() adjusts step size
    x = np.maximum(0, x)

    # Plug in this value of s to the RHS of the negative autoregulation model
    return neg_auto_rhs(x, t, beta0, gamma, k, n, ks, ns, s)

Now that we have this new function in hand, we can numerically integrate our ODEs as we did before. We’ll start with a pulse that is on from roughly \(t=2\) to \(t=6\), as above.

# Set up parameters for the pulse
s_args = (4.0, 2.0)

# Package parameters into a tuple
args = (beta0, gamma, k, n, ks, ns, s_pulse, s_args)

# Integrate ODEs
x = scipy.integrate.odeint(neg_auto_rhs_s_fun, x0, t, args=args).transpose()[0]

# Plot the normalized values
x /= x.max()

# Also calculate the pulse for plotting purposes
s = s_pulse(t, *s_args)

# Plot the results
p = bokeh.plotting.figure(
    frame_width=450,
    frame_height=250,
    x_axis_label="time",
    y_axis_label="normalized concentration",
    x_range=[0, 10],
)

# Populate glyphs
p.line(t, s, line_width=2, color=colors[0], legend_label="s")
p.line(t, x, line_width=2, color=colors[1], legend_label="x")

# Plot aesthetics
p.legend.location = "top_right"
p.legend.click_policy = "hide"

# Show plot
bokeh.io.show(p)

Computing environment

%load_ext watermark
%watermark -v -p numpy,scipy,bokeh,colorcet,jupyterlab

Python implementation: CPython
Python version       : 3.10.10
IPython version      : 8.10.0

numpy     : 1.23.5
scipy     : 1.10.0
bokeh     : 3.1.0
colorcet  : 3.0.1
jupyterlab: 3.5.3