Analysis of the repressilator

(c) 2018 Justin Bois. With the exception of pasted graphics, where the source is noted, this work is licensed under a Creative Commons Attribution License CC-BY 4.0. All code contained herein is licensed under an MIT license.

This document was prepared at Caltech with financial support from the Donna and Benjamin M. Rosen Bioengineering Center.

This tutorial was generated from an Jupyter notebook. You can download the notebook here.

As usual, before we begin, I will import the modules necessary to do the analysis.

In [1]:
import warnings

# The workhorses
import numpy as np
import scipy.integrate
import scipy.optimize

# Plotting modules
import bokeh.models
import bokeh.plotting
import bokeh.io
import ipywidgets
bokeh.io.output_notebook()
Loading BokehJS ...

The dynamical equations and fixed point

The dynamical equations for the simplified repressilator are

\begin{align} \frac{\mathrm{d}x_1}{\mathrm{d}t} &= \frac{\beta}{1+x_3^n} - x_1 \\[1em] \frac{\mathrm{d}x_2}{\mathrm{d}t} &= \frac{\beta}{1+x_1^n} - x_2 \\[1em] \frac{\mathrm{d}x_3}{\mathrm{d}t} &= \frac{\beta}{1+x_2^n} - x_3 \end{align}

We derived in lecture that there is a single fixed point, $x_1 = x_2 = x_3 \equiv x_0$ with $\beta = x_0(1 + x_0^n)$. We will first write a function to compute this fixed point.

In [2]:
def f(x, beta, n):
    """
    Equation for nullcline for repressilator.
    """
    return beta / (1 + x**n)

def fixed_point(beta, n):
    """
    Computes the fixed point of the simplified dimensionless repressilator.
    """
    
    def opt_fun(x):
        return x - f(f(f(x, beta, n), beta, n), beta, n)
    
    return scipy.optimize.brentq(opt_fun, 0, 10)

Temporal dynamics

We can integrate the dynamical equations to see the levels of the respective proteins. We use interactive plotting of the result so we can see how the dynamics depend on the parameters $\beta$ and $n$. Note that to have these interactions, you need to be running this Jupyter notebook and have ipywidgets installed; the interactivity is lost in the static HTML rendering.

In [3]:
# Define right hand side of ODEs
def dx_dt(x, t, beta, n):
    """
    Returns 3-array of (dx_1/dt, dx_2/dt, dx_3/dt)
    """
    x_1, x_2, x_3 = x
    return np.array([beta / (1 + x_3**n) - x_1,
                     beta / (1 + x_1**n) - x_2,
                     beta / (1 + x_2**n) - x_3])

# Initial condiations
x0 = np.array([1, 1, 1.2])

# Time points
t = np.linspace(0, 30, 1000)

# Choose parameters
beta = 10.0
n = 3

# Solve it!
x = scipy.integrate.odeint(dx_dt, x0, t, args=(beta, n))

# Plot the solutions
source = bokeh.models.ColumnDataSource(
            data=dict(x1=x[:,0], x2=x[:,1], x3=x[:,2], t=t))
p = bokeh.plotting.figure(width=600, height=300, x_axis_label='t',
                           border_fill_alpha=0, background_fill_alpha=0)
p.line('t', 'x1', source=source, line_width=3, color='dodgerblue', legend='1')
p.line('t', 'x2', source=source, line_width=3, color='tomato', legend='2')
p.line('t', 'x3', source=source, line_width=3, color='slateblue', legend='3')
p.legend.location = 'top_left'
bokeh.io.show(p, notebook_handle=True)

# Set up callbacks
def update(n=3, beta=10):
    # Generate the new curve
    x = scipy.integrate.odeint(dx_dt, x0, t, args=(beta, n))

    # Re-source data
    source.data = dict(x1=x[:,0], x2=x[:,1], x3=x[:,2], t=t)
    bokeh.io.push_notebook()
In [4]:
ipywidgets.interactive(update, n=ipywidgets.FloatSlider(min=1, max=5, value=3), 
                       beta=ipywidgets.FloatSlider(min=1, max=100, value=10))

We clearly see oscillations and can investigate how the parameter values change them using the widgets.

Trajectory in phase space

Making a stream plot is difficult in multiple dimensions. Instead, we will plot the trajectory of $x_1$ and $x_2$ in the $x_1$-$x_2$ plane. We indicate the fixed point with a circle.

In [5]:
# Compute fixed point
x_fp = fixed_point(beta, n)

# Plot trajectory
source = bokeh.models.ColumnDataSource(data=dict(x1=x[:,0], x2=x[:,1]))
source_fp = bokeh.models.ColumnDataSource(data=dict(x1=[x_fp], x2=[x_fp]))
p = bokeh.plotting.figure(height=400, width=400, x_axis_label='x1', y_axis_label='x2')
fp = p.circle('x1', 'x2', source=source_fp, color='tomato', size=12)
p.line('x1', 'x2', source=source, line_width=2)
bokeh.io.show(p, notebook_handle=True)

# Set up callbacks
def update_traj(n=3, beta=10):
    # Generate the new curve
    x = scipy.integrate.odeint(dx_dt, x0, t, args=(beta, n))
    
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        x_fp = fixed_point(beta, n)

    # Re-source data
    source.data = dict(x1=x[:,0], x2=x[:,1])
    source_fp.data = dict(x1=np.array([x_fp]), x2=np.array([x_fp]))
    
    bokeh.io.push_notebook()
In [6]:
ipywidgets.interactive(update_traj, n=ipywidgets.FloatSlider(min=1, max=5, value=3), 
                       beta=ipywidgets.FloatSlider(min=1, max=100, value=10))

Linear stability diagram

It is useful to make a linear stability diagram, which is a map of phase space highlighting stable and unstable regions. We know the bifurcation line is

\begin{align} \beta = \frac{n}{2}\left(\frac{n}{2} - 1\right)^{-\frac{n+1}{n}} \end{align}

We can plot this line and delineate the regions of stability and instability.

In [7]:
# Get bifurcation line
n = np.linspace(2.01, 5, 200)
beta = n / 2 * (n / 2 - 1)**(-(1 + 1/n))

p = bokeh.plotting.figure(height=400, width=500, x_axis_label='n', y_axis_label='β',
                          y_axis_type='log', x_range=[2, 5])
p.patch(np.append(n, n[-1]), np.append(beta, beta[0]), color='seashell', alpha=0.7)
p.line(n, beta, line_width=4, color='black')
p.text(x=2.1, y=2, text=['stable'])
p.text(x=3.5, y=100, text=['unstable'])
bokeh.io.show(p);