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.

---

import numpy as np
import scipy.integrate

import bokeh.io
import bokeh.plotting

import panel as pn

import colorcet

bokeh.io.output_notebook()
pn.extension()

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 lesson, we will put these tools to use analyzing a genetic circuit.

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} \to \mathrm{Y} \to \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 more quantitatively 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 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. We can do this by defining dimensionless parameters as follows.

\begin{align} \tilde{t} &= \gamma_y t, \\[2mm] \tilde{x} &= x/k_x, \\[2mm] \tilde{y} &= y/k_y, \\[2mm] \tilde{z} &= \gamma_z\,z/\beta_z, \\[2mm] \gamma &= \gamma_z/\gamma_y, \\[2mm] \beta &= \frac{\beta_y}{\gamma_y k_y}, \\[2mm] \end{align}

Inserting these into the dynamical equations, our dimensionless ODEs are

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

For notational convenience, and since we will always be working in dimensionless units, we will drop the tildes.

\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\). Note that this is two fewer parameter than the original equations. Nondimensionalization typically results in a reduction of the number of parameters that vary independently.

The scipy.intergrate module

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. 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\), we suddenly have a concentration of X of \(x_0\). So, we need six parameters for the right hand side of our ODEs, \(\beta_y\), \(\beta_z\), \(\gamma\), \(n_x\), \(n_y\), and \(x_0\).

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

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.

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

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

yz.shape

(400, 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 homeworks. (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.

# 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

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,
):
    # 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, 400)
    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.

p = cascade_response_plot()

# Show plot
bokeh.io.show(p)

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

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.

For our interactivity, we will use Panel. (Note that we imported Panel as pn at the beginning of this notebook and executed pn.extension() to enable Panel to work.) Panel allows convenient functionality for making interactive graphics. As you will see throughout the course, interactive graphics make it easier to explore how circuits work; you can interact directly with plots instead of having to re-run calculations and manually re-plot.

As a first step in creating a graphic, we create sliders, which enable us to graphically change values of parameters. I’ll first build an example slider in Panel.

pn.widgets.FloatSlider(name='example', start=0, end=1, step=0.01, value=0.5)

Data type cannot be displayed:

We would like to build such sliders for all of the parameter of our model.

log_beta_slider = pn.widgets.FloatSlider(
    name="log₁₀ β", start=-1, end=2, step=0.1, value=0
)
log_gamma_slider = pn.widgets.FloatSlider(
    name="log₁₀ γ", start=-1, end=2, step=0.1, value=0
)
nx_slider = pn.widgets.FloatSlider(
    name="nx", start=0.1, end=10, step=0.1, value=2
)
ny_slider = pn.widgets.FloatSlider(
    name="ny", start=0.1, end=10, step=0.1, value=2
)
log_x0_slider = pn.widgets.FloatSlider(
    name="log₁₀ x0", start=-1, end=2, step=0.1, value=np.log10(2)
)
t_max_slider = pn.widgets.FloatSlider(
    start=1, end=20, step=0.1, value=10, name="t_max",
)

The normalize parameter is a binary choice, so we cannot use a slider for it. Instead, we can use a checkbox.

normalize_checkbox = pn.widgets.Checkbox(name='Normalize', value=False)

Once the sliders and checkbox are built, adding the interactivity to the plot is as simple as adding a @pn.depends() decorator to the plotting function. Specifically, we set things up so that the parameters of the function that generates the plot depend on the values of the sliders. We simply add the decorator to the function we already wrote to do this.

@pn.depends(
    log_beta_slider.param.value,
    log_gamma_slider.param.value,
    nx_slider.param.value,
    ny_slider.param.value,
    log_x0_slider.param.value,
    t_max_slider.param.value,
    normalize_checkbox.param.value,
)
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,
):
    # 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, 400)
    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, we need to lay things out. We will construct the widgets (the sliders) in a column using the pn.Column() class and then build the plot next to the widgets using the pn.Row() class. I’ll add a little spacer to keep things well-spaced.

Note: The interactivity from Panel will not be available in the HTML version of this notebook because Python needs to be running to perform the updates. All subsequent interactive plots will appear static in the HTML rendering of this notebook.

widgets = pn.Column(
    pn.Spacer(height=10),
    log_beta_slider,
    log_gamma_slider,
    nx_slider,
    ny_slider,
    t_max_slider,
    log_x0_slider,
    width=150,
)

left_column = pn.Column(
    pn.Row(pn.Spacer(width=30), normalize_checkbox), cascade_response_plot,
)

# Final layout
pn.Row(left_column, pn.Spacer(width=20), widgets)

Data type cannot be displayed:

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 lesson, 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.

Speeding up interactions and Hill functions

The process of building the interactive plot may have seemed like a lot of steps. However, if we look back at what we did to build it, we did the following.

1. Wrote a function to make a plot.

2. Built widgets to control the parameters of the plot.

3. Make a layout positioning the plot(s) and widgets.

You pretty much do step 1 whenever you make any plot. Panel makes the interactivity easy by providing widgets and function decorators.

When we built the plots above, we wrote a function that constructs a new plot every time we change a parameter value. This works well, but it can become slow, since we are re-generating the entire graphic with each update. To have faster responses, we could instead update only the data that is plotted.

To illustrate how this work, and to get a feel for Hill functions in the meantime, we can look at how tuning the Hill coefficient and Hill \(k\) values affect the shape of a Hill function.

We start by writing a function to make a plot of an increasing (in blue) and decreasing (in orange) Hill function. We assume that \(k\) and the concentration \(x\) are in the same (unspecified) units. We will not need any Panel decorator here, as will become clear momentarily.

def hill_plot(k=1, n=1):
    # Set up plot
    p = bokeh.plotting.figure(
        frame_width=450,
        frame_height=250,
        x_axis_label="x",
        y_axis_label="f(x;k,n)",
        x_range=[0, 10],
    )

    # Set up empty data source
    x = np.linspace(0, 10, 400)
    y_decreasing = 1 / (1 + (x / k) ** n)
    y_increasing = 1 - y_decreasing
    source = bokeh.models.ColumnDataSource(
        dict(x=x, y_increasing=y_increasing, y_decreasing=y_decreasing)
    )

    # Populate glyphs
    p.line("x", "y_increasing", source=source, line_width=2, color=colors[0])
    p.line("x", "y_decreasing", source=source, line_width=2, color=colors[1])

    return p

Next, we can build sliders for \(k\) and \(n\).

k_slider = pn.widgets.FloatSlider(
    name="k", start=0.1, end=10, step=0.1, value=1, width=150
)
n_slider = pn.widgets.FloatSlider(
    name="n", start=0.1, end=10, step=0.1, value=1, width=150
)

Next, we need to make sure the Bokeh plot is in a Panel pane, so that we can link it to the slider, which we need to do explicitly since we’re updating the data source and not just replotting.

p_pane = pn.pane.Bokeh(hill_plot())

Next, we need a callback. This is a function that gets executed every time the slider changes. In our previous example, by default, the callback was to regenerate the entire plot. Here, we update the data source in our callback. We have to get into the guts of the Bokeh figure, pulling out the glyph renderer and adjusting its properties, including its data source.

def k_callback(target, event):
    # New value of k
    k = event.new

    # Fetch value of n
    n = n_slider.value

    # Extract data source
    gr = target.object.renderers[0]
    source = gr.data_source

    # Update data source
    x = source.data["x"]
    source.data["y_decreasing"] = 1 / (1 + (x / k) ** n)
    source.data["y_increasing"] = 1 - source.data["y_decreasing"]


def n_callback(target, event):
    # Fetch value of k
    k = k_slider.value

    # New value of n
    n = event.new

    # Extract data source
    gr = target.object.renderers[0]
    source = gr.data_source

    # Update data source
    x = source.data["x"]
    source.data["y_decreasing"] = 1 / (1 + (x / k) ** n)
    source.data["y_increasing"] = 1 - source.data["y_decreasing"]

Finally, we need to link the sliders with the plot.

k_slider.link(p_pane, callbacks={'value': k_callback})
n_slider.link(p_pane, callbacks={'value': n_callback})

Watcher(inst=FloatSlider(end=10, name='n', start=0.1, value=1, width=150), cls=<class 'panel.widgets.slider.FloatSlider'>, fn=<function Reactive.link.<locals>.link at 0x1525a59950>, mode='args', onlychanged=True, parameter_names=('value',), what='value', queued=False)

The result that was printed to the screen simply says that we have set up a watcher so that the plot will get updated whenever the slider is changed. Now, we can look at our result!

pn.Row(
    p_pane,
    pn.Spacer(width=15),
    pn.Column(pn.Spacer(height=40), k_slider, pn.Spacer(height=15), n_slider),
)

Data type cannot be displayed:

The response of sliders constructed this way can be much quicker than in the simpler method we did before, which regenerates the plot. However, as you can see from this example, you have to write more code to implement the faster-responding sliders. Going forward, we will almost always use the simpler construction, where we simple write a function that generates a plot and include that in our Panel layout.

Duration of input

Coming back to our cascade circuit, now imagine that the input is a pulse of duration \(\tau\). We can write a function for this and plot it.

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)

# Show plot
bokeh.io.show(p)

If we want to solve the ODEs for a pulsed input, we need to have a way to pass this function 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.

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\).

# 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\).

# 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\).

# 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, since the Y signal never went up. 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.

# 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 advantages to use scipy.integrate.odeint() due to its low overhead.

Periodic input

Finally, 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\).

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.

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

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

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

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

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

CPython 3.7.7
IPython 7.13.0

numpy 1.18.1
scipy 1.4.1
bokeh 2.0.1
panel 0.9.3
colorcet 2.0.2
jupyterlab 1.2.6