2. Introduction to Python for biological circuits¶
(c) 2019 Justin Bois and Michael Elowitz. 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 lesson was generated from a Jupyter notebook. Click the links below for other versions of this lesson.
Prior to executing this lesson, you should be sure to update the biocircuits module by doing the following on the command line.
pip install --upgrade biocircuits import numpy as np
import scipy.integrate
import biocircuits
import bokeh.application
import bokeh.application.handlers
import bokeh.io
import bokeh.models
import bokeh.plotting
bokeh.io.output_notebook()
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} &= z/k_y, \\[2mm] \gamma &= \gamma_z/\gamma_y, \\[2mm] \tilde{\beta}_y &= \frac{\beta_y}{\gamma_y k_y}, \\[2mm] \tilde{\beta}_z &= \frac{\beta_z}{\gamma_z k_y}. \end{align}Inserting these into the dynamical equations, our dimensionless ODEs are
\begin{align} \frac{\mathrm{d}\tilde{y}}{\mathrm{d}\tilde{t}} &= \tilde{\beta}_y\,\frac{\tilde{x}^{n_x}}{1+\tilde{x}^{n_x}} - \tilde{y}, \\[1em] \gamma^{-1}\,\frac{\mathrm{d}\tilde{z}}{\mathrm{d}\tilde{t}} &= \tilde{\beta}_z\,\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_y\,\frac{x^{n_x}}{1+x^{n_x}} - y, \\[1em] \gamma^{-1}\,\frac{\mathrm{d}z}{\mathrm{d}t} &= \beta_z\,\frac{y^{n_y}}{1+y^{n_y}} - z. \end{align}Thus, in addition to the specifics of our input $x(t)$, we have five parameters, $\beta_y$, $\beta_z$, $\gamma$, $n_x$, and $n_y$. Note that this is one fewer parameter than the original equations. Nondimensionalization typically results in a reduction of the number of parameters that vary independently.
The scipy.intergrate 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 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_y, beta_z, 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_y * x**n_x / (1 + x**n_x) - y
# Compute dz/dt
dz_dt = gamma * (beta_z * 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_y = 1.0
beta_z = 1.0
gamma = 1.0
n_x = 2
n_y = 2
x_0 = 2.0
# Package parameters into a tuple
args = (beta_y, beta_z, gamma, n_x, n_y, x_0)
# 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
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.)
# Set up color palette for this notebook
colors = bokeh.palettes.d3['Category10'][10]
# Pluck out y and z
y, z = yz.transpose()
# Set up plot
p = bokeh.plotting.figure(plot_width=500,
plot_height=300,
x_axis_label='dimensionless time',
y_axis_label='dimensionless y, z')
# Populate glyphs
p.line(t, y, line_width=2, color=colors[0], legend='y')
p.line(t, z,line_width=2, color=colors[1], legend='z')
# Place the legend
p.legend.location = 'bottom_right'
#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.
To make an interactive plot in Bokeh, there are three major components.
- 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/Falsestate 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
ColumnDataSourceof the plot. AColumnDataSourceis 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 changes in the range of the axis of the plot due to zooming. - The function handler. The Bokeh object takes a Python function and uses it to make a plot. The function you provide sets up the entire plot, including the callbacks.
We will start by setting up the widgets, starting with sliders. For response of the cascade circuit to the sudden presence of X, there are six parameters. For each slider, we need to specify the following.
| property | description |
|---|---|
title |
Title of slider, essentially the name of the variable |
start |
Lower bound on the slider values |
end |
Upper bound on the slider values |
step |
Step size for the slider, going from start to end |
value |
Initial value of the slider |
It is convenient to store the properties of the parameters as a tuple of objects with attributes. To do this, we set up a generic class, which we will call an AttributeContainer.
# Generic class to hold attributes
class AttributeContainer(object):
def __init__(self, **kw):
self.__dict__ = kw
# Parameters for each slider
slider_params = (AttributeContainer(title='log₁₀ βy', start=-1, end=2, value=0, step=0.1),
AttributeContainer(title='log₁₀ βz', start=-1, end=2, value=0, step=0.1),
AttributeContainer(title='log₁₀ γ', start=-1, end=2, value=0, step=0.1),
AttributeContainer(title='nx', start=1, end=10, value=2, step=0.1),
AttributeContainer(title='ny', start=1, end=10, value=2, step=0.1),
AttributeContainer(title='log₁₀ x₀', start=-1, end=2, value=np.log10(2), step=0.1))
Note that we have chosen to vary the parameters $\beta_y$, $\beta_z$, $\gamma$, and $x_0$ on a logarithmic scale. The ordering of the parameters in the tuple is important, since this is the ordering that is expected by the callback function. We can also define Boolean inputs for the plot that will be included as toggles. For our interactive plot of the cascade dynamics, we will have a single toggle that determines whether or not the concentrations of the respective species are normalized such that their maximum value is unity. We specify this again as a tuple of dictionaries (in this case a 1-tuple). The dictionary contains the name of the toggle and its initial value.
toggle_params = (AttributeContainer(title='normalize', active=False),)
We can now turn our attention to the callback function. The callback function for a slider must have call signature callback_function(attr, old, new). For reasons that will become clear in a moment, we should limit the scope of the callback function to be within the function we will pass to the function handler. This is as simple as defining a function that calls a user-supplied callback function with the necessary arguments. It looks like this:
def _callback(attr, old, new):
supplied_callback(source=source,
x_range=p.x_range,
y_range=p.y_range,
sliders=sliders,
toggles=toggles,
*extra_args)
In our interactive plot, it is important to note that the function supplied_callback function has the call signature shown above. It takes as input a ColumnDataSource, the x_range and y_range of the plot, a tuple of arguments from the sliders, a tuple of arguments from the toggles, and finally any extra arguments it needs. Its main purpose is to update the ColumnDataSource according to the inputs from the sliders, toggles, range of the x-axis, and any extra arguments it needs. For the callback for the cascade system, or for any plot that is the solution to a system of ODEs, the extra_args contain the right-hand side of the ODEs (cascade_rhs), the initial conditions (yz_0), and the number of time points to consider (n). Here is the supplied callback for the cascade system.
extra_args = (cascade_rhs, yz_0, n)
def cascade_callback(source, x_range, y_range, sliders, toggles, cascade_rhs, yz_0, n):
# Set up time values, keeping minimum at zero
t = np.linspace(0, x_range.end, n)
# Convert logarithmic sliders
slider_args = tuple(10**slider.value if 'log' in slider.title else slider.value
for slider in sliders)
# Integrate ODES
yz = scipy.integrate.odeint(cascade_rhs, yz_0, t, args=slider_args)
# Normalize if desired
if toggles[0].active:
yz[:,0] /= yz[:,0].max()
yz[:,1] /= yz[:,1].max()
# Update data source
y, z = yz.transpose()
source.data['t'] = t
source.data['y_vals'] = y
source.data['z_vals'] = z
Now that we have our callback function, we can use it to generate our plot. We need to write a function to do this because when we build the interactive plot, the plot objects we create have to be self-containing in a function we pass to the FunctionHandler. The function must return both a Bokeh figure, and a ColumnDataSource for the figure. Given the function we will write for making the interactive plot, the function must have the call signature shown below.
def cascade_plot(callback, sliders, toggles, extra_args):
# Set up plot
p = bokeh.plotting.figure(plot_width=500,
plot_height=300,
x_axis_label='dimensionless time',
y_axis_label='dimensionless y, z')
# Set up empty data source
source = bokeh.models.ColumnDataSource()
# Populate glyphs
p.line('t', 'y_vals', source=source, line_width=2, color=colors[0], legend='y')
p.line('t', 'z_vals', source=source, line_width=2, color=colors[1], legend='z')
# Place the legend
p.legend.location = 'bottom_right'
# Update data according to callback
callback(source, bokeh.models.Range1d(0, 10), None,
sliders, toggles, *extra_args)
return p, source
The base plot looks good, and now we are ready to make it interactive. One important note in making it interactive: If you want the range of the x-axis (in this case, time) to change upon zooming, you need to set the padding of the x_range and y_range to zero, using p.x_range.range_padding = 0. Otherwise, if the callback function resets the x-values of the plot every time the x-range is changed, Bokeh will then pad the x-range. But this means the x-range was changed, so the callback function will reset the x-values, Bokeh will pad it, and so on ad infinitum.
def interactive_xy_plot(base_plot, callback, slider_params, toggle_params, extra_args):
"""Build an interactive x-y plot app."""
def _plot_app(doc):
# Build the initial plot and data source
p, source = base_plot(callback, slider_params, toggle_params, extra_args)
# Make sure axis ranges have no padding
if type(p.x_range) == bokeh.models.ranges.Range1d:
start, end = p.x_range.start, p.x_range.end
p.x_range = bokeh.models.ranges.DataRange1d(p.x_range)
p.x_range.start = start
p.x_range.end = end
if type(p.y_range) == bokeh.models.ranges.Range1d:
start, end = p.y_range.start, p.y_range.end
p.y_range = bokeh.models.ranges.DataRange1d(p.y_range)
p.y_range.start = start
p.y_range.end = end
p.x_range.range_padding = 0
p.y_range.range_padding = 0
# Callbacks
def _callback(attr, old, new):
callback(source, p.x_range, p.y_range, sliders, toggles,
*extra_args)
# Callback for the toggle with required call signature
def _callback_toggle(new):
_callback(None, None, new)
# Set up sliders
sliders = tuple(bokeh.models.Slider(start=param.start,
end=param.end,
value=param.value,
step=param.step,
title=param.title)
for param in slider_params)
for slider in sliders:
slider.on_change('value', _callback)
# Set up toggles
toggles = tuple(bokeh.models.Toggle(label=param.title) for param in toggle_params)
for toggle in toggles:
toggle.on_click(_callback_toggle)
# Execute callback upon changing axis values
p.x_range.on_change('start', _callback)
p.x_range.on_change('end', _callback)
p.y_range.on_change('start', _callback)
p.y_range.on_change('end', _callback)
# Add the plot to the app
widgets = bokeh.layouts.widgetbox(*sliders, *toggles)
doc.add_root(bokeh.layouts.column(widgets, p))
handler = bokeh.application.handlers.FunctionHandler(_plot_app)
return bokeh.application.Application(handler)
We are now finally ready to build our app and interact with it. When showing the interactive app in JupyterLab, you should specify the URL of the notebook. It is almost always 'localhost:8888', but it may have a different port (for example 'localhost:8889') if you have more than one instance of JupyterLab running.
Note: The plot will not display nor be interactive in the HTML version of this notebook; it will only be interactive in a live JupyterLab session. This is because Python functions are being called to update the plot. If you want interactivity in a notebook rendered into HTML, the interactions must be pure JavaScript.
# Build the interactive plotting app
app = interactive_xy_plot(cascade_plot, cascade_callback,
slider_params, toggle_params, extra_args)
# Show the app
notebook_url = 'localhost:8888'
bokeh.io.show(app, notebook_url=notebook_url)
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.
A short aside on interactive plots 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.
- Specify the sliders and toggles as tuples of
AttributeContainerinstances. - Specify the callback function.
- Specify the function to build the plot.
- Use the
interactive_xy_plot()to build the app. - Show the app.
You pretty much do steps 2, 3, and 4 whenever you make any plot. You just need to be sure to specify a ColumnDataSource in the case of interactive plots. Step 1 is necessary to specify what you want to have control of with your widgets. The only other step is step 4. But the function we wrote is general, and can work for any x-y plot with any callbacks. So, we can reuse the interactive_xy_plot() function whenever we need to make an interactive plot. Now that we have gone through building this function, we can use it whenever we like. It is included in the biocircuits package, as is the AttributeContainer class.
To investigate, we can look at how tuning the Hill coefficient and Hill $k$ values affect the shape of a Hill function.
slider_params = (biocircuits.AttributeContainer(title='k', start=0.1,
end=10, value=1, step=0.05),
biocircuits.AttributeContainer(title='n', start=1,
end=10, value=1, step=0.1))
extra_args = (400,)
def hill_callback(source, x_range, y_range, sliders, toggles, n_points):
# Set up x values, keeping minimum at zero
x = np.linspace(0, x_range.end, n_points)
# Pull out slider values
k = sliders[0].value
n = sliders[1].value
# Compute Hill function
source.data['x'] = x
source.data['y'] = (x/k)**n / (1 + (x/k)**n)
def hill_plot(callback, sliders, toggles, extra_args):
# Set up plot
p = bokeh.plotting.figure(plot_width=500,
plot_height=300,
x_axis_label='x',
y_axis_label='f(x;k,n)')
# Set up empty data source
source = bokeh.models.ColumnDataSource()
# Populate glyphs
p.line('x', 'y', source=source, line_width=2)
# Update data according to callback
callback(source, bokeh.models.Range1d(0, 10), None,
sliders, toggles, *extra_args)
return p, source
hill_app = biocircuits.interactive_xy_plot(hill_plot, hill_callback,
slider_params, (), extra_args)
bokeh.io.show(hill_app, notebook_url=notebook_url)
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(plot_width=500,
plot_height=300,
x_axis_label='dimensionless time',
y_axis_label='dimensionless x')
# 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_y, beta_z, 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_y, beta_z, 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_y, beta_z, 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(plot_width=500,
plot_height=300,
x_axis_label='dimensionless time',
y_axis_label='dimensionless y, z')
# Populate glyphs
p.line(t, y, line_width=2, color=colors[0], legend='y')
p.line(t, z,line_width=2, color=colors[1], legend='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 3$.
# Set up parameters for the pulse (on at t = 1, off at t = 3, x_0 = 2)
x_args = (1.0, 0.7, 2.0)
# Package parameters into a tuple
args = (beta_y, beta_z, 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(plot_width=500,
plot_height=300,
x_axis_label='dimensionless time',
y_axis_label='dimensionless y, z')
# Populate glyphs
p.line(t, y, line_width=2, color=colors[0], legend='y')
p.line(t, z,line_width=2, color=colors[1], legend='z')
# Place the legend
p.legend.location = 'top_right'
#Show plot
bokeh.io.show(p)
We see that Z basically does not respond 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_y, beta_z, 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(plot_width=500,
plot_height=300,
x_axis_label='dimensionless time',
y_axis_label='dimensionless y, z')
# Populate glyphs
p.line(t, y, line_width=2, color=colors[0], legend='y')
p.line(t, z,line_width=2, color=colors[1], legend='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_y, beta_z, 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_y, beta_z, 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_y, beta_z, 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(plot_width=500,
plot_height=300,
x_axis_label='dimensionless time',
y_axis_label='dimensionless y, z')
# Populate glyphs
p.line(t, y, line_width=2, color=colors[0], legend='y')
p.line(t, z,line_width=2, color=colors[1], legend='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.
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(plot_width=500,
plot_height=300,
x_axis_label='dimensionless time',
y_axis_label='dimensionless y, z')
# 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 cooperativity
n_y = 10
args = (beta_y, beta_z, 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(plot_width=500,
plot_height=300,
x_axis_label='dimensionless time',
y_axis_label='dimensionless y, z')
# Populate glyphs
p.line(t, y, line_width=2, color=colors[0], legend='y')
p.line(t, z, line_width=2, color=colors[1], legend='z')
p.line(t, x, line_width=2, color=colors[2], alpha=0.2, legend='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_y, beta_z, 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(plot_width=500,
plot_height=300,
x_axis_label='dimensionless time',
y_axis_label='dimensionless y, z')
# Populate glyphs
p.line(t, y, line_width=2, color=colors[0], legend='y', line_join='bevel')
p.line(t, z, line_width=2, color=colors[1], legend='z', line_join='bevel')
p.line(t, x, line_width=2, color=colors[2], alpha=0.2, legend='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_y, beta_z, 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(plot_width=500,
plot_height=300,
x_axis_label='dimensionless time',
y_axis_label='dimensionless y, z')
# Populate glyphs
p.line(t, y, line_width=2, color=colors[0], legend='y', line_join='bevel')
p.line(t, z, line_width=2, color=colors[1], legend='z', line_join='bevel')
p.line(t, x, line_width=2, color=colors[2], alpha=0.2, legend='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, 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,jupyterlab,biocircuits