2. Introduction to Python for biological circuits¶

Design principle

A cascade can filter out high frequency fluctuations.

Techniques

Nondimensionalization.
Numerical solution of ODEs using scipy.integrate.odeint().
Interactive plotting.

[1]:

# Colab setup ------------------
import os, sys, subprocess
if "google.colab" in sys.modules:
    cmd = "pip install --upgrade colorcet watermark"
    process = subprocess.Popen(cmd.split(), stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    stdout, stderr = process.communicate()
# ------------------------------

import numpy as np
import scipy.integrate

import bokeh.io
import bokeh.layouts
import bokeh.models
import bokeh.plotting

import colorcet

# Set to True to have fully interactive plots
fully_interactive_plots = False
notebook_url = "localhost:8888"

bokeh.io.output_notebook()

Loading BokehJS ...

You have already installed Anaconda and have launched a Jupyter notebook. You have the tools you need to use your computer to analyze biological circuits. In this chapter, we will put these tools to use analyzing a genetic circuit. In doing so, we will demonstrate key techniques you will use over and over again when working with biological circuits.

Our model system¶

As we learn how to use some Python-based tools for analysis of genetic circuits, we will learn an important design principle about cascades. A three-component cascade is illustrated below.

\begin{align} \mathrm{X} \longrightarrow \mathrm{Y} \longrightarrow \mathrm{Z} \end{align}

Here, X is our input, which we will specify. X could be something like an externally imposed stimulus. We are interested in the response of Y and Z as a function of input X. We will expose the design principle that a cascade can filter out high frequency fluctuations. This makes sense intuitively. Imagine we have a brief pulse in X. Y will start to grow in concentration during the pulse, but will then fall back toward its basal level when the pulse stops. If there is not enough time for Y to accumulate enough to enhance production of Z, then the level of Z does not change appreciably. Thus, Z shows no response to a brief pulse.

The dynamical equations¶

We will model the circuit by writing down a system of ordinary differential equations and solving them. We will assume Hill-like behavior for the activation of Y by X and of Z by Y. (We will describe Hill functions, fortuitously named after their inventor and their shape, in the next chapter. For now, we take this regulatory function as given.) We also neglect leakage. We define the concentrations of X, Y, and Z, respectively as \(x\), \(y\), and \(z\). The system of ODEs describing this system is then

\begin{align} \frac{\mathrm{d}y}{\mathrm{d}t} &= \beta_y\,\frac{(x/k_x)^{n_x}}{1+(x/k_x)^{n_x}} - \gamma_y y, \\[1em] \frac{\mathrm{d}z}{\mathrm{d}t} &= \beta_z\,\frac{(y/k_y)^{n_y}}{1+(y/k_y)^{n_y}} - \gamma_z z. \end{align}

Note that the imposed \(x(t)\) is in general a function of time.

Nondimensionalization¶

As is generally a good idea for analysis of these systems, we will nondimensionalize the dynamical equations. The variables and parameters in the problem have dimension. For example, \(y\) has dimension of concentration, or number of particles per cubic length, written as

\begin{align} y = \left[\frac{N}{L^3}\right]. \end{align}

The parameter \(\gamma_y\) has dimension of inverse time, or \(\gamma_y = [T^{-1}]\). It is clear that each term in the dynamical equations has dimension of \(N/L^3 T\). In general, every term in an equation must have the same dimension; trying to add, for example, a meter to a kilogram is nonsensical.

The nondimensionalization procedure involves rewriting the equations such that every term is dimensionless, which means that every term has dimension of unity. This is beneficial for several reasons.

It reduces the number of parameters you need to consider. Only dimensionless ratios/products of parameters are considered.
It allows comparison of the relative magnitudes of terms in an equation.
It provides intuition by identifying which ratios and products of parameters determine the dynamics.

There is no one way to nondimensionalize a set of dynamical equations. We have found the following procedure to be generally effective. Let \(\theta\) be one of the variables, space, or time. In the above dynamical equations, \(\theta\) could be \(x\), \(y\), \(z\), or \(t\). (Here \(x\), \(y\), and \(z\) denote concentrations; there are no spatial variables in this problem.) By contrast, \(\beta_y\), \(\beta_z\), \(k_x\), \(k_y\), \(n_x\), \(n_y\), \(\gamma_y\), and \(\gamma_z\) are system parameters, not variables.

For each variable and time, define a dimensionless version (usually marked by a tilde) such that \(\theta = \theta_d \tilde{\theta}\), where \(\theta\) is some variable or time. The constant \(\theta_d\) imparts the dimension on \(\theta\) (hence the subscript d).
Substitute these expressions into the dynamical equations.
Rearrange the equations such that every term is dimensionless. As a result, you find ratios and products of the \(\theta_d\)’s and the system parameters. These ratios are dimensionless and are called dimensionless parameters.
Choose expressions for the \(\theta_d\)’s that minimize the number of dimensionless parameters.

As usual, this is best seen by example. We start be defining

\begin{align} &t = t_d\,\tilde{t},\\[1em] &x = x_d\,\tilde{x},\\[1em] &y = y_d\,\tilde{y},\\[1em] &z = z_d\,\tilde{z}. \end{align}

Inserting these into the dynamical equations gives

\begin{align} &\frac{y_d}{t_d}\,\frac{\mathrm{d}\tilde{y}}{\mathrm{d}\tilde{t}} = \beta_y\,\frac{(x_d\tilde{x}/k_x)^{n_x}}{1+(x_d\tilde{x}/k_x)^{n_x}} - \gamma_y y_d \tilde{y}, \\[1em] &\frac{z_d}{t_d}\frac{\mathrm{d}\tilde{z}}{\mathrm{d}\tilde{t}} = \beta_z\,\frac{(y_d\tilde{y}/k_y)^{n_y}}{1+(y_d\tilde{y}/k_y)^{n_y}} - \gamma_z z_d\tilde{z}. \end{align}

To make every term dimensionless, we can divide the top equation by \(y_d/t_d\) and the bottom equation by \(z_d/t_d\). The result is

\begin{align} &\frac{\mathrm{d}\tilde{y}}{\mathrm{d}\tilde{t}} = \frac{\beta_y\,t_d}{y_d}\,\frac{(x_d\tilde{x}/k_x)^{n_x}}{1+(x_d\tilde{x}/k_x)^{n_x}} - \gamma_y t_d \tilde{y}, \\[1em] &\frac{\mathrm{d}\tilde{z}}{\mathrm{d}\tilde{t}} = \frac{\beta_z\,t_d}{z_d}\,\frac{(y_d\tilde{y}/k_y)^{n_y}}{1+(y_d\tilde{y}/k_y)^{n_y}} - \gamma_z t_d\tilde{z}. \end{align}

We identify our dimensionless parameters as \begin{align} &\beta_y\,t_d/y_d, \\[1em] &\beta_z\,t_d/z_d, \\[1em] &x_d/k_x,\\[1em] &y_d/k_y,\\[1em] &\gamma_y t_d,\\[1em] &\gamma_z t_d, \end{align}

with \(n_x\) and \(n_y\) already being dimensionless. We see immediately that we can eliminate two of the dimensionless parameters if we choose \(x_d = k_x\) and \(y_d = k_y\). We can eliminate another by choosing \(t_d = 1/\gamma_y\). (We could have also chosen \(t_d = 1 / \gamma_z\).) Thus, we have \(\tilde{\gamma} \equiv \gamma_z / \gamma_y\) and \(\tilde{\beta} \equiv \beta_y/\gamma_y k_y\) as a dimensionless parameters. Finally, we can choose \(z_d = \beta_z/\gamma_z\) such that \(\beta_z\,t_d/z_d = \tilde{\gamma}\). This last choice does not add a dimensionless parameter, since we have already defined \(\tilde{\gamma}\).

Thus, our dimensionless equations are

\begin{align} &\frac{\mathrm{d}\tilde{y}}{\mathrm{d}\tilde{t}} = \tilde{\beta}\,\frac{\tilde{x}^{n_x}}{1+\tilde{x}^{n_x}} - \tilde{y}, \\[1em] &\frac{\mathrm{d}\tilde{z}}{\mathrm{d}\tilde{t}} = \tilde{\gamma}\,\frac{\tilde{y}^{n_y}}{1+\tilde{y}^{n_y}} - \tilde{\gamma} \tilde{z}. \end{align}

For notational convenience, and since we will always be working in dimensionless units, we will drop all of the tildes. Doing this and dividing the bottom equation by \(\tilde{\gamma}\) gives

\begin{align} \frac{\mathrm{d}y}{\mathrm{d}t} &= \beta \,\frac{x^{n_x}}{1+x^{n_x}} - y, \\[1em] \gamma^{-1}\,\frac{\mathrm{d}z}{\mathrm{d}t} &= \frac{y^{n_y}}{1+y^{n_y}} - z. \end{align}

Thus, in addition to the specifics of our input \(x(t)\), we have four parameters, \(\beta\), \(\gamma\), \(n_x\), and \(n_y\). We already know that \(n_x\) and \(n_y\) are Hill coefficients that parametrize the ultrasensitivity of the regulation by X and Y respectively. The parameter \(\beta = \beta_y/\gamma_y k_y\) is the ratio of the unregulated level of Y, \(\beta_y/\gamma\) to the level of Y that activates Z. The parameter \(\gamma\) is the ratio of the decay rate of Z to that of Y. If dilution is the primary mode of concentration depletion for both Y and Z, then \(\gamma \approx 1\). We now have a clear picture of which dimensionless parameters dictate the dynamics and what they mean physically.

The scipy.intergrate module¶

Now that we have our dimensionless dynamical equations, we will move to solving them.

TheSciPy Library is a Python library for scientific computing. It contains many modules, including scipy.stats, scipy.special, and scipy.optimize, which respectively have functions to perform statistical calculations, special functions, and optimization routines. There are many more. 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, I find that the generic LSODA algorithm developed by Linda Petzold and Alan Hindmarsh that handles both stiff and non-stiff problems with variable time stepping is the best option. This is 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, this is not a big deal). Since most of our problems are small, we will use odeint(). It has much better performance, and though its API is different, it is still intuitive.

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 cascade circuit.

Solving for a constant input X¶

We will first consider the case where we initially have no X, Y, or Z present. At time \(t = 0\), the concentration of X is suddenly increased to \(x_0\). So, we need five parameters for the right hand side of our ODEs, \(\beta\), \(\gamma\), \(n_x\), \(n_y\), and \(x_0\).

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

[2]:

def cascade_rhs(yz, t, beta, gamma, n_x, n_y, x):
    """
    Right hand side for cascade X -> Y -> Z.  Return dy/dt and dz/dt.
    """
    # Unpack y and z
    y, z = yz

    # Compute dy/dt
    dy_dt = beta * x ** n_x / (1 + x ** n_x) - y

    # Compute dz/dt
    dz_dt = gamma * (y ** n_y / (1 + y ** n_y) - z)

    # Return the result as a NumPy array
    return np.array([dy_dt, dz_dt])

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

[3]:

# Number of time points we want for the solutions
n = 1000

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

# Initial condition
yz_0 = np.array([0.0, 0.0])

# Parameters
beta = 1.0
gamma = 1.0
n_x = 2
n_y = 2
x0 = 2.0

# Package parameters into a tuple
args = (beta, gamma, n_x, n_y, x0)

# Integrate ODES
yz = scipy.integrate.odeint(cascade_rhs, yz_0, t, args=args)

That’s it! The integration is done. We can now look at what scipy.integrate.odeint()’s output looks like.

[4]:

yz.shape

[4]:

(1000, 2)

The first column of the output yz gives \(y(t)\) at the specified time points and the second column gives \(z(t)\). We would now like to plot the results.

Plotting results¶

We will use Bokeh to plot the results. The syntax is pretty self-explanatory from the example below. Note that you can save a plot as a PNG by clicking the disk icon next to the plot, which might be helpful for incorporating your plots into your presentations. (Note that for publications, you should usually save your figures in a vector graphics format, which Bokeh supports, but is not necessary, and in fact discouraged, for display of plots in the browser.)

Before building the plot, we will load in the color scheme we will use. The colorcet package is good for this. I like to use the Category10 color palette for categorical colors.

[5]:

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

Colors in place, we will write a function to set up the plot. Within this function, we will write the solver; not just the plotting function. It will become clear that this is useful for interactive graphics. Because transcription rates, decay rates, and input concentrations of X can vary over orders of magnitude, we will input those as their logarithms. Finally, we will allow for normalization of the responses. That is, we allow for rescaling of the Y and Z concentrations such that their maximum concentration is one.

[6]:

def cascade_response_plot(
    log_beta=0,
    log_gamma=0,
    n_x=2,
    n_y=2,
    log_x0=np.log10(2),
    t_max=10,
    normalize=False,
):
    """Generate a plot showing response of the cascade circuit
    to a step in x0.
    """
    # Package parameters into a tuple
    args = (
        10 ** log_beta,
        10 ** log_gamma,
        n_x,
        n_y,
        10 ** log_x0,
    )

    # Integrate ODES
    t = np.linspace(0, t_max, 1000)
    yz = scipy.integrate.odeint(cascade_rhs, yz_0, t, args=args)
    y, z = yz.transpose()

    # Normalize Y and Z responses
    if normalize:
        y /= y.max()
        z /= z.max()

    # Set up plot
    p = bokeh.plotting.figure(
        frame_width=325,
        frame_height=250,
        x_axis_label="dimensionless time",
        y_axis_label=f"{'normalized ' if normalize else ''}dimensionless y, z",
        x_range=[0, t_max],
    )

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

    # Place the legend
    p.legend.location = "bottom_right"

    return p

Now that we have our function, we call it to generate the graphic, and then we use bokeh.io.show() to display the graphic in the notebook. Note that at the top of the notebook, we called bokeh.io.output_notebook() which tells bokeh.io.show() to display the plot in the notebook instead of writing it out to a file.

[7]:

p = cascade_response_plot()

# Show plot
bokeh.io.show(p)

We see that the cascade acts as a delay for changes in \(z\) as a result of input \(x\).

Duration of input¶

Now imagine that the input is a pulse of duration \(\tau\). That is, the concentration of X instantaneously rises to \(x_0\), stays at \(x_0\) for time \(\tau\), and then instantaneously decreases back to zero. We can write a function for this and plot it.

[8]:

def x_pulse(t, t_0, tau, x_0):
    """
    Returns x value for a pulse beginning at t = 0
    and ending at t = t_0 + tau.
    """
    return np.logical_and(t >= t_0, t <= (t_0 + tau)) * x_0


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

# Populate glyphs
p.line(t, x_pulse(t, 1.0, 2.0, 2), line_width=2, color=colors[2])

# Show plot
bokeh.io.show(p)

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 that takes x_fun, the function describing \(x(t)\) as an argument, as well as x_args, the set of parameters passed into x_fun.

[9]:

def cascade_rhs_x_fun(yz, t, beta, gamma, n_x, n_y, x_fun, x_args):
    """
    Right hand side for cascade X -> Y -> Z.  Return dy/dt and dz/dt.

    x_fun is a function of the form x_fun(t, *x_args), so x_args is a tuple
    containing the arguments to pass to x_fun.
    """
    # Compute x
    x = x_fun(t, *x_args)

    # Return cascade RHS with this value of x
    return cascade_rhs(yz, t, beta, gamma, n_x, n_y, x)

With this in hand, we can now solve for a pulse. We will have a pulse during \(1 \le t \le 5\).

[10]:

# Set up parameters for the pulse (on at t = 1, off at t = 5, x_0 = 2)
x_args = (1.0, 4.0, 2.0)

# Package parameters into a tuple
args = (beta, gamma, n_x, n_y, x_pulse, x_args)

# Integrate ODEs
yz = scipy.integrate.odeint(cascade_rhs_x_fun, yz_0, t, args=args)

# Pluck out y and z
y, z = yz.transpose()

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

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

# Place the legend
p.legend.location = "top_right"

# Show plot
bokeh.io.show(p)

Let’s see what happens when we do a shorter pulse, this time with \(1 \le t \le 1.7\).

[11]:

# Set up parameters for the pulse (on at t = 1, off at t = 1.7, x_0 = 2)
x_args = (1.0, 0.7, 2.0)

# Package parameters into a tuple
args = (beta, gamma, n_x, n_y, x_pulse, x_args)

# Integrate ODEs
yz = scipy.integrate.odeint(cascade_rhs_x_fun, yz_0, t, args=args)

# Pluck out y and z
y, z = yz.transpose()

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

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

# Place the legend
p.legend.location = "top_right"

# Show plot
bokeh.io.show(p)

We see that Z basically does not respond strongly to a short pulse. The delay of the circuit allows short pulses to be ignored, but large pulses to be detected and responded to.

Really short pulses and a lesson about scipy.integrate.odeint()¶

Now, we will take a brief interlude to learn an important lesson about the algorithm of scipy.integrate.odeint() and its use in these applications. We will consider a very brief pulse, \(1 \le t \le 1.05\).

[12]:

# Set up parameters for the pulse (on at t = 1, off at t = 1.05, x_0 = 2)
x_args = (1.0, 0.05, 2.0)

# Package parameters into a tuple
args = (beta, gamma, n_x, n_y, x_pulse, x_args)

# Integrate ODEs
yz = scipy.integrate.odeint(cascade_rhs_x_fun, yz_0, t, args=args)

# Pluck out y and z
y, z = yz.transpose()

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

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

# Place the legend
p.legend.location = "top_right"

# Show plot
bokeh.io.show(p)

Uh oh! Something went wrong; the Y signal never went up at all! This exposes an important issue with the algorithm used by scipy.integrate.odeint(). The Hindmarsh-Petzold algorithm uses variable step sizes so that it takes long steps when the system is not changing much and short steps when it is. Therefore, if we have a long period of no changes (leading up to \(t = 1\)), the step sizes taken by the solver will increase, and we’ll step right over the pulse.

So, it is in general good practice to explicitly take into account discontinuities in the parameters over time. In this case, we would use scipy.integrate.odeint() to integrate to the pulse and use the end point of that as the initial condition of a new solution during the pulse. Then, at the end of the pulse, we start again. Let’s try again using this method.

[13]:

# Integrate prior to pulse
t_before_pulse = np.linspace(0, 1.0, 20)
args = (beta, gamma, n_x, n_y, 0.0)
yz_0 = np.array([0.0, 0.0])
yz_before_pulse = scipy.integrate.odeint(
    cascade_rhs, yz_0, t_before_pulse, args=args
)

# Integrate during pulse
t_during_pulse = np.linspace(1.0, 1.05, 50)
args = (beta, gamma, n_x, n_y, 2.0)
yz_0 = yz_before_pulse[-1]
yz_during_pulse = scipy.integrate.odeint(
    cascade_rhs, yz_0, t_during_pulse, args=args
)

# Integrate after pulse
t_after_pulse = np.linspace(1.05, 5, 50)
args = (beta, gamma, n_x, n_y, 0.0)
yz_0 = yz_during_pulse[-1]
yz_after_pulse = scipy.integrate.odeint(
    cascade_rhs, yz_0, t_after_pulse, args=args
)

# Piece together solution
t = np.concatenate((t_before_pulse, t_during_pulse[1:], t_after_pulse[1:]))
yz = np.vstack(
    (yz_before_pulse, yz_during_pulse[1:, :], yz_after_pulse[1:, :])
)
y, z = yz.transpose()

# Plot the results
p = bokeh.plotting.figure(
    frame_width=450,
    frame_height=250,
    x_axis_label="dimensionless time",
    y_axis_label="dimensionless y, z",
    x_range=[0, 5],
)

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

# Place the legend
p.legend.location = "top_right"

# Show plot
bokeh.io.show(p)

Much better. Dealing with discontinuities like this while solving systems of ODEs are facilitated by event handling. Unfortunately, scipy.integrate.odeint() does not allow for event handling, and we need to take the less elegant approach we just showed.

Note that scipy.integrate.solve_ivp() does have event handling capabilities. However, we will seldom need event handing for our purposes, and it is still advantageous to use scipy.integrate.odeint() due to its low overhead.

Periodic input¶

As our final variety of input, we will consider the case where we have periodic forcing of the circuit. We will do this for highly cooperative activation by Y, taking \(n_y = 10\). Recall that this gives a longer time delay.

We first write a function for the forcing, x_fun, which is periodic with frequency \(f\).

[14]:

def x_periodic(t, f, x_0):
    """
    Returns x value for periodic forcing of amplitude x_0 and frequency f.
    """
    if type(f) in [float, int]:
        return x_0 * (1 + np.sin(f * t))
    else:
        sin_sum = np.zeros_like(t)
        for freq, amp in zip(f, x_0):
            sin_sum += amp * (1 + np.sin(freq * t))
        return sin_sum


# Plot the forcing
t = np.linspace(0, 10, 500)

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

# Populate glyphs
p.line(t, x_periodic(t, 5, 2.0), line_width=2, color=colors[2])

# Show plot
bokeh.io.show(p)

Let’s see how the circuit responds to a low-frequency input.

[15]:

# Set up parameters for periodic forcing with f = 0.5 and x_0 = 2.
x_args = (0.5, 2.0)

# Package parameters into a tuple, now with high ultrasensitivity
n_y = 10
args = (beta, gamma, n_x, n_y, x_periodic, x_args)

# Time points
t = np.linspace(0, 40, 300)

# Initial condition
yz_0 = np.array([0.0, 0.0])

# Integrate ODES
yz = scipy.integrate.odeint(cascade_rhs_x_fun, yz_0, t, args=args)

# Pluck out y and z
y, z = yz.transpose()

# x
x = x_periodic(t, *x_args)
x /= x.max()

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

# Populate glyphs
p.line(t, y, line_width=2, color=colors[0], legend_label="y")
p.line(t, z, line_width=2, color=colors[1], legend_label="z")
p.line(
    t,
    x,
    line_width=2,
    color=colors[2],
    alpha=0.2,
    legend_label="x (normalized)",
    line_join="bevel",
)

# Place the legend
p.legend.location = "top_right"

# Show plot
bokeh.io.show(p)

We roughly follow the forcing with some lag. Now, for high-frequency forcing, we have a different response.

[16]:

# Set up parameters for periodic forcing with f = 5 and x_0 = 2.
x_args = (5.0, 2.0)

# Package parameters into a tuple
args = (beta, gamma, n_x, 10, x_periodic, x_args)

# Time points
t = np.linspace(0, 25, 600)

# Initial condition
yz_0 = np.array([0.0, 0.0])

# Integrate ODES
yz = scipy.integrate.odeint(cascade_rhs_x_fun, yz_0, t, args=args)

# Pluck out y and z
y, z = yz.transpose()

# x
x = x_periodic(t, *x_args)
x /= x.max()

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

# Populate glyphs
p.line(
    t, y, line_width=2, color=colors[0], legend_label="y", line_join="bevel"
)
p.line(
    t, z, line_width=2, color=colors[1], legend_label="z", line_join="bevel"
)
p.line(
    t,
    x,
    line_width=2,
    color=colors[2],
    alpha=0.2,
    legend_label="x (normalized)",
    line_join="bevel",
)

# Place the legend
p.legend.location = "top_right"

# Show plot
bokeh.io.show(p)

We see that Z does not really respond to high frequency forcing, even though the forcing is with the same amplitude. This gives us a design principle, that a cascade can filter out high frequency fluctuations. We can see this by adding another frequency to the signal.

[17]:

# Set up parameters for periodic forcing with f = 5 and x_0 = 2.
x_args = ((0.5, 10.0), (0.5, 2.0))

# Package parameters into a tuple
args = (beta, gamma, n_x, 10, x_periodic, x_args)

# Time points
t = np.linspace(0, 25, 600)

# Initial condition
yz_0 = np.array([0.0, 0.0])

# Integrate ODES
yz = scipy.integrate.odeint(cascade_rhs_x_fun, yz_0, t, args=args)

# Pluck out y and z
y, z = yz.transpose()

# x
x = x_periodic(t, *x_args)
x /= x.max()

# Plot the results
p = bokeh.plotting.figure(
    frame_width=450,
    frame_height=250,
    x_axis_label="dimensionless time",
    y_axis_label="dimensionless y, z",
)

# Populate glyphs
p.line(
    t, y, line_width=2, color=colors[0], legend_label="y", line_join="bevel"
)
p.line(
    t, z, line_width=2, color=colors[1], legend_label="z", line_join="bevel"
)
p.line(
    t,
    x,
    line_width=2,
    color=colors[2],
    alpha=0.2,
    legend_label="x (normalized)",
    line_join="bevel",
)

# Place the legend
p.legend.location = "top_right"

# Show plot
bokeh.io.show(p)

Even though the high frequency part of the forcing has a bigger amplitude, the signal in z responds predominantly to the low frequency part of the signal.

Interactive plotting and varying parameters¶

Plotting with Bokeh in Jupyter notebooks allows for interactivity with plots that can help to rapidly gain insights about how parameter values might affect the dynamics. We have found that this is a useful tool to rapidly explore parameter dependence on circuit behavior. We will construct an interactive plot for the cascade circuit to demonstrate how it is done.

To make an interactive plot in Bokeh, there are three major components.

The plot or plots themselves.
The widgets. Widgets for parameter values are primarily sliders, which enable you to vary parameter values by clicking and dragging. We will also make use of a toggle, which is a pushbutton whose True/False state dictates a feature of the plot.
The callback function. This is a function that is executed whenever a widget changes value. Most of the time, we use it to update the ColumnDataSource of the plot. A ColumnDataSource is a Bokeh object that dictates where the glyphs are placed on the plot. You may have more than one callback functions for different sliders, toggles, and also for changes in the range of the axis of the plot due to zooming.
The layout. This is the spatial arrangement of the plots and widgets.
The app. Bokeh will create an application that can be embedded in a notebook. To create, it you need to make a simple function that adds the layout you built to the document that Bokeh will make into an app. (This sounds a lot more complicated than it is; see the example below.)

We refer to a plot or set of plots with widgets for interactivity as a dashboard.

Note

The excellent package Panel allows for more declarative (and hence fewer lines of code) means of dashboarding. We intend to use Panel exclusively for this purpose in the future. At the present point in Panel and JupterLab’s development cycle, though, there are problems. As of Panel version 0.10.3, Panel graphics do not natively render in JupyterLab version 3.0.0 and above. Furthermore, when serving up dashboards with Panel (that is, serving them so they render by themselves outside of JupyterLab in their own browser tab), there are problems communicating with Bokeh’s JavaScript libraries.

We expect this to be fixed in short order, and in the future will likely replace the pure-Bokeh means of dashboarding presented here. Learning how to build the dashboards using only Bokeh is still quite valuable because it offers greater flexibility and performance. It is also not that much more effort to code up.

A simple contrived example¶

To demonstrate dashboard construction, we will first build a simple example, plotting a sine wave with a slider to adjust the frequency.

Step 1: Generate the plot (but don’t show it)¶

First, we will generate the plot. When we generate the plot, we specify the data in a ColumnDataSource, a Bokeh class that allows us to change the data on a plot. When we populate the plot with glyphs, we use the source keyword argument to specify the ColumnDataSource, and then specify x and y as strings indicating which columns of the ColumnDataSource are used to specify the \(x\) and \(y\) values, respectively.

We cannot show the plot here, since we will show it in the app. A Bokeh plot can only be in a single document, in our case the app (which as far as Bokeh is concerned is separate from showing it in the JupyterLab cell below).

[18]:

# x-y data for plotting
x = np.linspace(0, 1, 200)
y = np.sin(2 * np.pi * x)

# Place the data in a ColumnDataSource
cds = bokeh.models.ColumnDataSource(dict(x=x, y=y))

# Build the plot
p = bokeh.plotting.figure(
    frame_height=200,
    frame_width=400,
    x_axis_label="x",
    y_axis_label="sin(2π f x)",
    x_range=[0, 1],
    y_range=[-1.1, 1.1]
)
p.line(source=cds, x="x", y="y", line_width=2);

Step 2: Make the widgets¶

Now we will make our widget, in this case a slider.

[19]:

freq_slider = bokeh.models.Slider(
    title="freq", start=0, end=5, step=0.01, value=1, width=150
)

The slider is instantiated using bokeh.models.Slider(), with keyword arguments whose meaning should be obvious from their names. The value attribute of the slider is the present value of the slider; we started the example off at a value of 1.

Step 3: Make the callbacks¶

Next, we specify the callback function that will be used to update the plot as the slider changes. The callback function for a slider must take three arguments, the attribute that changes, its old value, and its new one. We will not directly use these arguments (though we could), but will rather directly read the value from the slider using its value attribute. Our callback function simply updates the y data values of the ColumnDataSource.

[20]:

def callback(attr, old, new):
    cds.data["y"] = np.sin(2 * np.pi * freq_slider.value * cds.data["x"])

    # Could also do: cds.data["y"] = np.sin(2 * np.pi * new * cds.data["x"])

Now, we need to alert Bokeh that it should trigger the callback function whenever the slider value changes.

Note

In the static HTML rendering of this notebook, the app will not be responsive. This is because the app needs a running Python engine for the interactivity, and no such engine is available in the browser. Also, Bokeh apps currently are not supported by Google Colab.

Therefore, in all chapters with Bokeh apps, we have a variable fully_interactive_plots, which is set in the first code cell along with the imports. This is set to False for static HTML rendering. In this case we do not link the sliders to callbacks (since there is no Python engine), and disable the slider.

In a later chapter, we will show how you can get interactivity in an HTML rendering of a notebook using JavaScript.

[21]:

if fully_interactive_plots:
    freq_slider.on_change("value", callback)
else:
    freq_slider.disabled = True

Step 4: Build the layout¶

Now that we have a slider and a plot and have linked the data source of the plot to the slider, we can lay out the dashboard. We will put the slider next to the plot, putting a spacer in between.

[22]:

layout = bokeh.layouts.row(p, bokeh.models.Spacer(width=30), freq_slider)

Step 5: Make the app¶

Now we are ready to make a function to produce the app. The function needs to have call signature app(doc), where doc represents the document that Bokeh will build into an app. The purpose of this function is to add the layout we have just build to the document.

[23]:

def app(doc):
    doc.add_root(layout)

Step 6: Enjoy your interactive plot!¶

And we are finally ready to see the app! When we call bokeh.io.show(), we pass the app() function as the first argument. We also need to use the notebook_url keyword argument to specify where the notebook is being hosted. This is usually "localhost:8888", but the number may change (e.g., 8889 or 8890). You can look in the navigation bar of your browser to make sure you get the right number. In this and all other chapters, we specify notebook_url in the first cell of the notebook along with the imports.

Note

In the static HTML rendering of this notebook, the app will not appear. This is because it is running a Python engine for the interactivity, and no such engine is available in the browser. In a later chapter, we will show how you can get interactivity in an HTML rendering of a notebook using JavaScript.

Also, Bokeh apps currently are not supported by Google Colab.

[24]:

if fully_interactive_plots:
    bokeh.io.show(app, notebook_url=notebook_url)
else:
    bokeh.io.show(layout)

Serving an app¶

If you wanted to have a stand-alone interactive plot on its own browser tab, you can put all of the code necessary to generate it in a single .py file and then serve it from the command line. For this example, we could have a file sine_wave.py with the following contents.

import numpy as np
import bokeh.plotting
import bokeh.models

# x-y data for plotting
x = np.linspace(0, 1, 200)
y = np.sin(2 * np.pi * x)

# Place the data in a ColumnDataSource
cds = bokeh.models.ColumnDataSource(dict(x=x, y=y))

# Build the plot
p = bokeh.plotting.figure(
    frame_height=200,
    frame_width=400,
    x_axis_label="x",
    y_axis_label="sin(2πfx)",
    x_range=[0, 1],
    y_range=[-1.1, 1.1]
)
p.line(source=cds, x="x", y="y", line_width=2)

freq_slider = bokeh.models.Slider(
    title="freq", start=0, end=5, step=0.01, value=1, width=150
)

def callback(attr, old, new):
    cds.data["y"] = np.sin(2 * np.pi * freq_slider.value * cds.data["x"])

freq_slider.on_change("value", callback)

layout = bokeh.layouts.row(p, bokeh.models.Spacer(width=30), freq_slider)

def app(doc):
    doc.add_root(layout)

# Build the app in the current doc
app(bokeh.plotting.curdoc())

After saving that file, you can serve it be doing the following on the command line.

bokeh serve --show sine_wave.py

An app for the cascade model¶

We now follow the same steps to build an interactive plot for the response of the cascade model to a step in input concentration \(x_0\).

Step 1: Build the plot¶

To build the plot, we use the same function as above, but we modify it to use a ColumnDataSource, which is also returned along with the plot.

[25]:

def cascade_response_plot(
    log_beta=0,
    log_gamma=0,
    n_x=2,
    n_y=2,
    log_x0=np.log10(2),
    t_max=10,
    normalize=False,
):
    """Generate a plot showing response of the cascade circuit
    to a step in x0.
    """
    # Package parameters into a tuple
    args = (
        10 ** log_beta,
        10 ** log_gamma,
        n_x,
        n_y,
        10 ** log_x0,
    )

    # Integrate ODES
    t = np.linspace(0, t_max, 1000)
    yz = scipy.integrate.odeint(cascade_rhs, yz_0, t, args=args)
    y, z = yz.transpose()

    # Normalize Y and Z responses
    if normalize:
        y /= y.max()
        z /= z.max()

    # Set up the column data source
    cds = bokeh.models.ColumnDataSource(dict(t=t, y=y, z=z))

    # Set up plot
    p = bokeh.plotting.figure(
        frame_width=325,
        frame_height=250,
        x_axis_label="dimensionless time",
        y_axis_label=f"{'normalized ' if normalize else ''}dimensionless y, z",
        x_range=[0, t_max],
    )

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

    # Place the legend
    p.legend.location = "bottom_right"

    return p, cds

Let’s call the function to get our plot and ColumnDataSource.

[26]:

p, cds = cascade_response_plot()

Step 2: Build the widgets¶

Now we will build our sliders, one for each parameter. We are interested in varying the parameters \(\beta\), \(\gamma\), and \(x_0\) on a logarithmic scale, so we will set up the sliders to specify \(\log_{10} \beta\), \(\log_{10}\gamma\) and \(\log_{10} x_0\).

[27]:

log_beta_slider = bokeh.models.Slider(
    title="log₁₀ β", start=-1, end=2, step=0.1, value=0, width=150
)
log_gamma_slider = bokeh.models.Slider(
    title="log₁₀ γ", start=-1, end=2, step=0.1, value=0, width=150
)
nx_slider = bokeh.models.Slider(
    title="nx", start=0.1, end=10, step=0.1, value=2, width=150
)
ny_slider = bokeh.models.Slider(
    title="ny", start=0.1, end=10, step=0.1, value=2, width=150
)
log_x0_slider = bokeh.models.Slider(
    title="log₁₀ x0", start=-1, end=2, step=0.1, value=np.log10(2), width=150
)

We also want to be able to toggle between display of normalized versus unnormalized responses of \(y\) and \(z\). We can make a toggle button for that.

[28]:

normalize_toggle = bokeh.models.Toggle(label='Normalize', active=False, width=50)

Step 3: Build the callbacks¶

We can now build our callback. This callback will be a bit more complicated. For each new slider value, we need to re-integrate the dynamical equations. We also need to check to see if the normalization toggle button is clicked and appropriately scale the data and change the y-axis label. We also will need to recalculate the result if the range of the time axis changes, so we need to read the time points off of the x_range property of the plot.

[29]:

def cascade_callback(attr, old, new):
    # Set up time values, keeping minimum at zero
    t = np.linspace(0, p.x_range.end, 2000)

    # Package slider values
    args = (
        10 ** log_beta_slider.value,
        10 ** log_gamma_slider.value,
        nx_slider.value,
        ny_slider.value,
        10 ** log_x0_slider.value,
    )

    # Integrate ODES
    yz = scipy.integrate.odeint(cascade_rhs, yz_0, t, args=args)

    # Normalize if desired
    if normalize_toggle.active:
        yz[:, 0] /= yz[:, 0].max()
        yz[:, 1] /= yz[:, 1].max()
        p.yaxis.axis_label = "normalized dimensionless y, z"
    else:
        p.yaxis.axis_label = "dimensionless y, z"

    # Update data source
    y, z = yz.transpose()
    cds.data = dict(t=t, y=y, z=z)

Now we link the callback to the sliders, and also to the normalization toggle and the range of the time axis.

[30]:

if fully_interactive_plots:
    log_beta_slider.on_change("value", cascade_callback)
    log_gamma_slider.on_change("value", cascade_callback)
    nx_slider.on_change("value", cascade_callback)
    ny_slider.on_change("value", cascade_callback)
    log_x0_slider.on_change("value", cascade_callback)
    normalize_toggle.on_change("active", cascade_callback)
    p.x_range.on_change("end", cascade_callback)
else:
    log_beta_slider.disabled = True
    log_gamma_slider.disabled = True
    nx_slider.disabled = True
    ny_slider.disabled = True
    log_x0_slider.disabled = True
    normalize_toggle.disabled = True

Step 4: The layout¶

We can now lay things out. I will put the sliders and normalization toggle in a column next to the plot.

[31]:

layout = bokeh.layouts.row(
    p,
    bokeh.layouts.Spacer(width=30),
    bokeh.layouts.column(
        log_beta_slider,
        log_gamma_slider,
        nx_slider,
        ny_slider,
        log_x0_slider,
        normalize_toggle,
    ),
)

Step 5: The app¶

And now for the app! (Note again that this app will not appear in the HTML rendering of this notebook and will not work using Google Colab.)

[32]:

def app(doc):
    doc.add_root(layout)

if fully_interactive_plots:
    bokeh.io.show(app, notebook_url=notebook_url)
else:
    bokeh.io.show(layout)

Importantly, when exploring the plot interactively, we see that when Y cooperatively activates Z (that is, \(n_y\) is large), the delay is longer. The delay is also naturally longer as \(\gamma\) gets small. This makes sense, since \(\gamma\) is the ratio of the decay rate of Z to that of Y. As we saw in the previous chapter, the decay rate sets the speed of the response, and if Z decays more slowly than Y, its dynamics will be slower.

Note that the above conclusions are clearer when we plot the normalized curves.

A challenge¶

In the above analysis of periodic forcing, I made a fresh plot each time of the response to an oscillating input so that the plots would display nicely in the HTML rendered version of this notebook. Can you make the sample plot(s) by making a single interactive plot?

Computing environment¶

[33]:

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

Python implementation: CPython
Python version       : 3.8.8
IPython version      : 7.22.0

numpy     : 1.20.1
scipy     : 1.6.2
bokeh     : 2.3.1
colorcet  : 2.0.6
jupyterlab: 3.0.14