9. Blinking bacteria: The repressilator enables self-sustaining oscillations

---

Design principle

• Ultrasensitive negative feedback can generate periodic oscillations in cells.

• Oscillator phase can be used as a growth timer.

Techniques

• Composition of functions

• Linear stability analysis

• Linear stability diagrams

• Numerical calculation of a scalar fixed point

• Synthetic biology

---

# Colab setup ------------------
import os, sys, subprocess
if "google.colab" in sys.modules:
    cmd = "pip install --upgrade biocircuits 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 scipy.optimize

import biocircuits

import bokeh.plotting
import bokeh.io

import colorcet

# We will use Matplotlib to make a 3D plot
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline
%config InlineBackend.figure_format = "retina"

# Set to True to have fully interactive plots with Python;
# Set to False to use pre-built JavaScript-based plots
interactive_python_plots = False
notebook_url = "localhost:8888"

bokeh.io.output_notebook()

Loading BokehJS ...

---

How to design a biological clock

In this chapter we will explore ways to build a clock in a living cell. More specifically, we will discuss a synthetic genetic clock circuit called the Repressilator (Elowitz and Leibler, 2000). We first discuss the roles of oscillators and clocks in natural biological systems, and then ask how one might go about designing a synthetic clock circuit that can operate in a cell. To do this, we will use linear stability analysis, a broadly useful approach for analyzing diverse systems. We will also introduce the concept of a limit cycle attractor.

Clocks are necessary for life

Our lives depend on accurate time-keeping. This is enabled by generations of improvement in clocks. A major turning point in global navigation was the development of portable clocks, which enabled sailors to measure longitude. John Harrison engineered a series of increasingly precise, and beautiful, clocks, or “marine chronometers.” A big challenge was making them function accurately despite variations in temperature, humidity, and other conditions. To do this, he invented internal compensation mechanisms that made critical parameters robust to these conditions. Even today, navigation remains dependent on clocks. The Global Positioning System (GPS), which we use to navigate our cities, relies on precise atomic clocks installed in orbiting satellites.

Time-keeping is as fundamental to our inner lives as it is to our outer lives. The cell division cycle drives the cell through a series of states, including periods of growth, chromosome replication, separation of daughter cells, and so on. When cells are rapidly proliferating, this circuit acts as a kind of clock. Organisms also have circadian clocks that allow them to track and anticipate the daily environmental changes that occur during the 24 hour day. Plants have mechanisms to keep track of time on longer timescales, allowing them to anticipate seasonal changes. Hormones can cycle over timescales of weeks or months. On shorter timescales, some factors in the cell oscillate periodically in and out of the nucleus over timescales of minutes to hours. Countless processes in biology involve periodic cycles rather than steady states.

Circadian clocks function in single cells

Circadian clocks control our 24 hour cycles of sleep, hunger, and activity. Humans confined to environments with constant light and temperature still maintain circadian activity cycles with nearly 24-hour periods (Czeisler, et al., 1999). This experiment, along with the more familiar experience of jet lag, shows that our body has an internal clock that enforces temporal organization on our bodies and minds. The cell cycle represents another kind of oscillator that takes cells through repetitive cycles of growth and division. Hormones cycle on a range of timescales from hours to weeks. Plants contain circadian and seasonal clocks that control their movement and flowering, in response to time as well as light, temperature and other inputs.

How do biological clocks work? In 1971, Ronald Konopka and Seymour Benzer identified mutations that altered the circadian behavioral rhythms of fruit flies (Konompa and Benzer, 1971). Over the next decades, biologists discovered key molecular components that enable these clocks to function, including transcription factors, light sensors, and other components, and worked out many aspects of the molecular mechanism of circadian oscillations. This work led to the 2017 nobel prize for Jeffrey Hall, Michael Rosbash, and Michael Young “for their discoveries about how internal clocks and biological rhythms govern human life.”

In multicellular organisms, you might wonder whether the clock emerges from interactions among cells, or is controlled from within each cell. The answer seems to be both. Clocks synchronize between cells and organs. However, clocks are not solely a multicellular phenomenon: Individual mammalian fibroblasts exhibit free-running circadian oscillations (Welsh, et al., 1995), as one can see in this movie (Welsh, et al., 2004). Circadian behavior here is revealed through expression of a luminescent reporter gene knocked into the mPer2 circadian clock gene.

Time traces of individual cells from this movie show periodic 24 hour oscillations:

This work revealed a complex clock containing many interacting protein components.

In cyanobacteria, a three-protein phosphorylation circuit generates oscillations.

But clocks need not be complex. Even single cell cyanobacteria have precise, cell-autonomous circadian clocks (Mihalcescu, et al., 2004). To learn how they work, Kondo and colleagues genetically identified a single locus containing three genes required for the cyanobacterial clock. In 2005 they showed, amazingly, that these three proteins, plus ATP, were sufficient to reconstitute 24 hour oscillations of phosphorylation in a test tube (Nakajima, et al., 2005).

The mechanism for these oscillations was subsequently explained by Rust, et al. in terms of ordered phosphorylation and feedback (Rust, et al., 2007). One protein, KaiC, has two phosphorylation sites. During each cycle, one site is phosphorylated, then the second is phosphorylated, then the first is dephosphorylated, and finally the second is dephosphorylated. The use of two phosphorylation sites ensures the dynamics are at least two dimensional. A one dimensional dynamical system cannot oscillate (see Strogatz, 2015).

Limit cycles are ideal dynamical behaviors for clocks

So far in this course, we have discussed only one kind of “attractor”—the stable fixed point. When a system sits at a stable fixed point, it does not change over time. (While true for a continuous system, biological systems are made of discrete molecules whose dynamics are subject to stochastic fluctuations, or noise, and therefore they cannot really remain precisely at a single state. We will learn more about noise soon.) By contrast, in a functioning clock circuit, the state of the system constantly changes in a periodic fashion, progressing through a cyclic sequence of “phases” that returns it to its starting point. This cycling occurs on its own, without any external input. Furthermore, if such an ideal clock circuit is perturbed in some way—whether by environmental fluctuations or internal noise—it returns to its oscillatory trajectory. It “has to” oscillate.

The behavior we are describing is more formally known as a limit cycle. Stable limit cycles are defined by Strogatz as “isolated closed orbits”, meaning that the system goes around the limit cycle, and that neighboring points ultimately feed into the limit cycle. As a result, if the system is perturbed a little bit away from the limit cycle, it will tend to return back to it, as shown here:

Note that one can also have unstable limit cycles.

As Strogatz notes, linear systems such as a frictionless pendulum can produce a family of orbits but not a limit cycle, because multiplying any solution of

\begin{align} \frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t} =\mathsf{A} \cdot \mathbf{x} \end{align}

by a constant produces another solution. Limit cycle oscillators are thus inherently non-linear systems.

For our purposes, the key point is that limit cycle dynamics are ideal for natural or synthetic clock circuits because they ensure periodic oscillations that do not damp out over time and can resume even if temporarily perturbed.

Build to understand: Designing a synthetic clock circuit

In biology, clocks have long been metaphorically identified with the perplexing mystery of how the undirected process of evolution could generate precise behaviors from seemingly “messy” molecular components. Dawkins’s book, The Blind Watchmaker (Dawkins, 1986), personifies evolution as a watchmaker, who creates devices of astonishing precision (tissues and organisms) without being able to “see,” much less plan, what she is doing. How “hard” is it for evolution—or anyone, really—to make a biological clock?

We will now discuss the design of a synthetic clock circuit in bacteria. By seeing whether we can design a clock out of biological components, we can find out “how hard” it is to make a clock, what circuits are sufficient to produce self-sustaining limit cycle oscillations, and whether the sorts of models we write down can describe biological circuits operating in the more complex, and largely uncharacterized, intracellular milieu. Beyond their value in teaching us about circuit design, synthetic clock circuits can also provide useful functions, as we will see below.

Clock components should ideally be composable, orthogonal, and well-characterized

The first question we must ask is what components we can build a clock from. Among the best characterized regulators in biology are prokaryotic (bacterial) transcriptional repressors and their cognate target promoters. These components are modular, orthogonal, and composable. By modular, we mean that they can be taken out of their natural context and used to generate a new regulatory circuit. By orthogonal, we mean that a variety of variants exist that operate similarly, but independently. Thus, different repressors bind to distinct DNA binding sites. Composability is a stronger form of modularity in which a set of components can regulate each other in the same way. LEGOs are a familiar example of composability: each of the standard bricks can be stuck onto a similar type of brick from below and above. Transcription factors are composable because any one can be engineered to regulate any other simply by combining corresponding target promoter sequences with open reading frames for the transcription factors.

(Transcriptional activators, which we have already encountered, are also excellent components for synthetic design, and have been used to build oscillators. At the time this work was done, there were generally fewer examples that were as well-understood as repressors. Therefore, we will focus below on a circuit design built exclusively from repressors).

Design strategy

Based on the considerations above, we will try to design a biological circuit that generates limit cycle oscillations across a broad range of biochemical parameter values using a set of well-characterized, composable, and orthogonal transcriptional repressors.

One of the first designs one can imagine building with repressors is a “rock-scissors-paper” feedback loop composed of three repressors, each of which represses the next one, in a cycle:

This diagram refers to three specific repressors, TetR, λ cI, and LacI. We will discuss the rationale for choosing these below, after we work out the design. For now, the names of the repressors are unimportant, so we will refer to them as repressors 1, 2, and 3. Repressor 1 represses production of repressor 2, which in turn represses production of repressor 3. Finally, repressor 3 represses production or repressor 1, completing the loop.

This design is a three-component negative feedback loop (analogous to a “three ring oscillator” in electronics). If one were to turn up the level of the first protein in this system, it would lead to a decrease in the second, which would cause an increase in the third, and finally a decrease in the first. Thus, one can see, intuitively, that this system produces a negative feedback that tends to push back in the opposite direction to any perturbation, after the delay required to propagate the perturbation around the loop.

We can try to work out the dynamics of this system by intuitive reasoning. We might achieve a limit cycle oscillation: Say that repressor 1 is present at a high protein concentration, while repressors 2 and 3 are low. The high concentration of repressor 1 will limit expression of repressor 2, keeping its level low. This means that repressor 3 is free to be expressed. As its concentration grows, it will start to repress repressor 1. As repressor 1 goes down, repressor 2 is expressed in higher numbers. The increased repressor 2 concentration leads to less repressor 3. Then, repressor 1 comes back up again. So, we see a cycle, where repressor 1 is high, then repressor 3, and finally repressor 2.

However, this behavior is by no means guaranteed. We might equally well just get a stable steady state, where all three repressors evolve to intermediate values, each sufficient to keep its target repressor at the appropriate level to maintain its target at its steady-state level. This behavior would be much more boring.

In fact, both behaviors are possible.

So, our questions are now:

1. What kinds of behaviors would this circuit be expected to produce?

2. How can we engineer the circuit to favor, or better yet, guarantee, limit cycle oscillations?

Dynamical equations for the repressilator

To analyze the repressilator, we will, as usual, write down a set of differential equations for each of the proteins. For simplicity, we will assume perfect symmetry among the species, with each one having identical biochemical properties (this will not be true in the real system). While we initially consider only protein dynamics here, later on, we will ask how including mRNA modifies the conclusions.

\begin{align} \frac{\mathrm{d}x_1}{\mathrm{d}t} &= \frac{\beta}{1 + (x_3/k)^n} - \gamma x_1, \\[1em] \frac{\mathrm{d}x_2}{\mathrm{d}t} &= \frac{\beta}{1 + (x_1/k)^n} - \gamma x_2, \\[1em] \frac{\mathrm{d}x_3}{\mathrm{d}t} &= \frac{\beta}{1 + (x_2/k)^n} - \gamma x_3. \end{align}

In dimensionless units, we measure protein concentrations in units of \(k\), the relevant protein concentration scale, and time in units of \(\gamma^{-1}\), the only timescale in the system. With this, the equations become:

\begin{align} \frac{\mathrm{d}x_i}{\mathrm{d}t} &= \frac{\beta}{1 + x_j^n} - x_i, \quad \text{ with } (i,j) \text{ pairs } (1,3), (2,1), (3,2). \end{align}

Note that \(\beta\) has been redefined as \(\beta \leftarrow \beta/k\gamma\).

Temporal dynamics

We can integrate the dynamical equations to follow the concentrations of three proteins over time. Interactive plotting further allows us to see how these dynamics depend on the parameters \(\beta\) and \(n\). Oscillations are clear when the curves go up and down periodically in time.

Alternatively, it is also instructive to plot the trajectory of the system as a projection in the \(x_2\)-\(x_1\) plane (or in either of the other two planes this three-dimensional system can be projected onto). When the fixed point is stable, the trajectory in the \(x_2\)-\(x_1\) plane spirals into the fixed point. When it is unstable, the trajectory spirals away from it, eventually cycling around the fixed point to join a limit cycle, corresponding to oscillations.

The layout below shows both plots. For sufficiently large \(\beta\) and \(n\), we see beautiful 3-phase oscillations. A few things to notice:

• If \(n<2\), oscillations diminish over time.

• If \(n>2\), sustained oscillations occur, but only for large enough \(\beta\)

• If \(n=2\), see what happens…

def protein_repressilator_rhs(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])

# Number of points to use in plots
n_points = 1000

# Widgets for controlling parameters
beta_slider_protein = bokeh.models.Slider(title="β", start=0, end=100, step=0.1, value=10)
n_slider_protein = bokeh.models.Slider(title="n", start=1, end=5, step=0.1, value=3)

# Solve for species concentrations
def _solve_protein_repressilator(beta, n, t_max):
    t = np.linspace(0, t_max, n_points)
    x = scipy.integrate.odeint(protein_repressilator_rhs, x0, t, args=(beta, n))

    return t, x.transpose()


# Obtain solution for plot
t, x = _solve_protein_repressilator(beta_slider_protein.value, n_slider_protein.value, 40.0)

# Build the plot
colors = colorcet.b_glasbey_category10[:3]

p_rep = bokeh.plotting.figure(
    frame_width=550, frame_height=200, x_axis_label="t", x_range=[0, 40.0]
)

cds = bokeh.models.ColumnDataSource(data=dict(t=t, x1=x[0], x2=x[1], x3=x[2]))
labels = dict(x1="x₁", x2="x₂", x3="x₃")
for color, x_val in zip(colors, labels):
    p_rep.line(
        source=cds,
        x="t",
        y=x_val,
        color=color,
        legend_label=labels[x_val],
        line_width=2,
    )

p_rep.legend.location = "top_left"


# Set up plot
p_phase = bokeh.plotting.figure(
    frame_width=200, frame_height=200, x_axis_label="x₁", y_axis_label="x₂",
)

p_phase.line(source=cds, x="x1", y="x2", line_width=2)

# Set up callbacks
def _callback(attr, old, new):
    t, x = _solve_protein_repressilator(beta_slider_protein.value, n_slider_protein.value, p_rep.x_range.end)
    cds.data = dict(t=t, x1=x[0], x2=x[1], x3=x[2])


beta_slider_protein.on_change("value", _callback)
n_slider_protein.on_change("value", _callback)
p_rep.x_range.on_change("end", _callback)

# Build layout
protein_repressilator_layout = bokeh.layouts.column(
    p_rep,
    bokeh.layouts.Spacer(height=10),
    bokeh.layouts.row(
        p_phase,
        bokeh.layouts.Spacer(width=70),
        bokeh.layouts.column(beta_slider_protein, n_slider_protein,width=150),
    ),
)

# Build the app
def protein_repressilator_app(doc):
    doc.add_root(protein_repressilator_layout)


# Display
if interactive_python_plots:
    bokeh.io.show(protein_repressilator_app, notebook_url=notebook_url)
else:
    bokeh.io.show(biocircuits.jsplots.protein_repressilator())

Finally, here is a simple three-dimensional plot of the limit cycle in the space of the three protein concentrations.

# We can do with with Bokeh a viz.js, but for now, we use Matplotlib

# Resolve problem for β = 20 and n = 3
t = np.linspace(0, 50, 1000)
x = scipy.integrate.odeint(protein_repressilator_rhs, x0, t, args=(20, 3))

# Generate the plot
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.view_init(30, 30)
ax.plot(x[:, 0], x[:, 1], x[:, 2])
ax.set_xlabel("x₁")
ax.set_ylabel("x₂")
ax.set_zlabel("x₃");

What conditions are necessary and sufficient for oscillations?

Why does the system oscillate for some values of \(n\) and \(\beta\) but not others? To find out, we need to use linear stability analysis to see how the values of the key biochemical parameters, \(\beta\) and \(n\) control oscillations.

As a first step, we must identify fixed points by solving \(\mathrm{d}x_i/\mathrm{d}t = 0\) for each \(x_i\). This is accomplished by appealing to composite functions in Technical Appendix 9a. We find that there is a single fixed point, with \(x_1 = x_2 = x_3 \equiv x_0\), which satisfies

\begin{align} \beta = x_0(1+x_0^n). \end{align}

If this fixed point is stable, then the system will sit there, constant in time. On the other hand, if it is unstable, then there is the potential for limit cycle oscillations. The question is: Under what conditions is this fixed point unstable?

To answer this question, we assess the stability of the unique fixed point using linear stability analysis. The technique and the calculation for this system are described in Technical Appendix 9b. The result is that the fixed point is unstable for \(n > 2\) and

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

We can display the regions in parameter space, that is in the \(n\text{-}\beta\) plane, for which the system is stable and unstable in a linear stability diagram, shown below.

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

# Build the plot
p = bokeh.plotting.figure(
    height=300,
    width=400,
    x_axis_label="n",
    y_axis_label="β",
    y_axis_type="log",
    x_range=[2, 5],
    y_range=[1, 2000],
)
p.patch(
    np.append(n, n[-1]), np.append(beta, beta[0]), color="lightgray", alpha=0.7
)
p.line(n, beta, line_width=4, color="black")
p.text(x=2.1, y=2, text=["stable"])
p.text(x=2.5, y=100, text=["unstable (limit cycle oscillations)"])

bokeh.io.show(p)

The linear stability diagram clearly shows that, from a design point of view, it is desirable to make both \(n\) and \(\beta\) as high as possible.

Intuition from the protein-only model

This analysis shows two conditions that favor oscillations:

• Strong ultrasensitivity (large Hill coefficients)

• Strong promoters

In fact, these results can be understood intuitively: oscillations occur when the feedback “overshoots.” The sharper and stronger the response as one goes around the complete feedback loop, the longer and higher a pulse in one factor can grow before it is, inevitably, yanked back down by the feedback. Consistent with this view, there is a tradeoff between the length of the cycle (number of repressors in the loop) and the minimum Hill coefficient required (Elowitz, 1999).

Including mRNA in the model provides additional insights

In the above analysis, we only considered the three proteins themselves, and we neglected the mRNA dynamics. However, it would be of interest to understand how mRNA properties like stability and translation rate affect whether and how the circuit oscillations. Therefore, we will add additional equations for each mRNA species. We will continue to assume symmetry among species, with all mRNAs sharing the same parameter values. With this assumption, the dynamical equations become:

\begin{align} \frac{\mathrm{d}m_i}{\mathrm{d}t} &= \alpha + \frac{\beta_m}{1 + (x_j/k)^n} - \gamma_m m_i,\\[1em] \frac{\mathrm{d}x_i}{\mathrm{d}t} &= \beta_p m_i - \gamma_p x_i, \\[1em] \end{align}

with \(i,j\) pairs \((1,3), (2,1), (3,2)\). Here, we have introduced \(\rho\) to allow for leaky transcription. In dimensionless units, these equations are

\begin{align} \frac{\mathrm{d}m_i}{\mathrm{d}t} &= \beta\left(\rho + \frac{1}{1 + x_j^n}\right) - m_i, \\[1em] \gamma^{-1}\,\frac{\mathrm{d}x_i}{\mathrm{d}t} &= m_i - x_i, \end{align}

Here,

• \(\gamma \equiv \gamma_p/\gamma_m\) is the ratio of the two timescales in the system–the protein and mRNA degradation/decay rates.

• \(\beta = \beta_m\beta_p/\gamma_m\gamma_p k\) is a dimensionless promoter strength.

• \(\rho = \alpha/\beta_m\) is the relative strength of leaky versus regulated expression.

As above, to start, we will solve this system numerically and explore its dynamics for different parameter values.

# Sliders
beta_slider = bokeh.models.Slider(
    title="β",
    start=0,
    end=4,
    step=0.1,
    value=1,
    width=125,
    format=bokeh.models.CustomJSTickFormatter(code="return Math.pow(10, tick).toFixed(2)"),
)
gamma_slider = bokeh.models.Slider(
    title="γ",
    start=-3,
    end=0,
    step=0.1,
    value=0,
    width=125,
    format=bokeh.models.CustomJSTickFormatter(code="return Math.pow(10, tick).toFixed(3)"),
)
rho_slider = bokeh.models.Slider(
    title="ρ",
    start=-6,
    end=0,
    step=0.1,
    value=-3,
    width=125,
    format=bokeh.models.CustomJSTickFormatter(code="return Math.pow(10, tick).toFixed(6)"),
)
n_slider = bokeh.models.Slider(title="n", start=1, end=5, step=0.1, value=3, width=125)

def repressilator_rhs(mx, t, beta, gamma, rho, n):
    """
    Returns 6-array of (dm_1/dt, dm_2/dt, dm_3/dt, dx_1/dt, dx_2/dt, dx_3/dt)
    """
    m_1, m_2, m_3, x_1, x_2, x_3 = mx
    return np.array(
        [
            beta * (rho + 1 / (1 + x_3 ** n)) - m_1,
            beta * (rho + 1 / (1 + x_1 ** n)) - m_2,
            beta * (rho + 1 / (1 + x_2 ** n)) - m_3,
            gamma * (m_1 - x_1),
            gamma * (m_2 - x_2),
            gamma * (m_3 - x_3),
        ]
    )


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

# Number of points to use in plots
n_points = 1000


# Solve for species concentrations
def _solve_repressilator(log_beta, log_gamma, log_rho, n, t_max):
    beta = 10 ** log_beta
    gamma = 10 ** log_gamma
    rho = 10 ** log_rho
    t = np.linspace(0, t_max, n_points)
    x = scipy.integrate.odeint(repressilator_rhs, x0, t, args=(beta, gamma, rho, n))
    m1, m2, m3, x1, x2, x3 = x.transpose()
    return t, m1, m2, m3, x1, x2, x3


t, m1, m2, m3, x1, x2, x3 = _solve_repressilator(
    beta_slider.value,
    gamma_slider.value,
    rho_slider.value,
    n_slider.value,
    40.0,
)

cds = bokeh.models.ColumnDataSource(
    dict(t=t, m1=m1, m2=m2, m3=m3, x1=x1, x2=x2, x3=x3)
)

p_rep = bokeh.plotting.figure(
    frame_width=500,
    frame_height=200,
    x_axis_label="t",
    x_range=[0, 40.0],
)

colors = bokeh.palettes.d3["Category20"][6]
m1_line = p_rep.line(source=cds, x="t", y="m1", line_width=2, color=colors[1])
x1_line = p_rep.line(source=cds, x="t", y="x1", line_width=2, color=colors[0])
m2_line = p_rep.line(source=cds, x="t", y="m2", line_width=2, color=colors[3])
x2_line = p_rep.line(source=cds, x="t", y="x2", line_width=2, color=colors[2])
m3_line = p_rep.line(source=cds, x="t", y="m3", line_width=2, color=colors[5])
x3_line = p_rep.line(source=cds, x="t", y="x3", line_width=2, color=colors[4])

legend_items = [
    ("m₁", [m1_line]),
    ("x₁", [x1_line]),
    ("m₂", [m2_line]),
    ("x₂", [x2_line]),
    ("m₃", [m3_line]),
    ("x₃", [x3_line]),
]
legend = bokeh.models.Legend(items=legend_items)

p_rep.add_layout(legend, "right")

# Build the layout
layout = bokeh.layouts.column(
    bokeh.layouts.row(
        beta_slider,
        gamma_slider,
        rho_slider,
        n_slider,
        width=575,
    ),
    bokeh.layouts.Spacer(height=10),
    p_rep,
)

# Set up callbacks
def _callback(attr, old, new):
    t, m1, m2, m3, x1, x2, x3 = _solve_repressilator(
        beta_slider.value,
        gamma_slider.value,
        rho_slider.value,
        n_slider.value,
        p_rep.x_range.end,
    )
    cds.data = dict(t=t, m1=m1, m2=m2, m3=m3, x1=x1, x2=x2, x3=x3)

beta_slider.on_change("value", _callback)
gamma_slider.on_change("value", _callback)
rho_slider.on_change("value", _callback)
n_slider.on_change("value", _callback)
p_rep.x_range.on_change("end", _callback)

# Build the app
def repressilator_app(doc):
    doc.add_root(layout)


# Display
if interactive_python_plots:
    bokeh.io.show(repressilator_app, notebook_url=notebook_url)
else:
    bokeh.io.show(biocircuits.jsplots.repressilator())

Using a similar technique as for the simplified three-component system, we can show that there is a unique fixed point with \(m_i = x_i = x_0\) for all \(i\) with

\begin{align} (x_0 - \beta\rho)(1+x_0^n) = \beta. \end{align}

We can perform linear stability for this fixed point. The linear stability matrix is now 6 by 6 rather than 3 by 3. Nonetheless, we can derive analytically that when we get an eigenvalue with a positive real part, it also has a nonzero imaginary part, which means that the instability is oscillatory. We get (not derived here; you can try it in the problems) an eigenvalue with positive real part when

\begin{align} \left(\sqrt{\gamma} + \sqrt{\gamma^{-1}}\right)^2 < \frac{3f_0^2}{4+2f_0}, \end{align}

where

\begin{align} f_0 = \frac{\beta n x_0^{n-1}}{(1+x_0^n)^2}. \end{align}

Note something interesting here: \(x_0\) is independent of \(\gamma\), while the \(\left(\sqrt{\gamma} + \sqrt{\gamma^{-1}}\right)^2\) term is invariant to exchanging \(\gamma \leftrightarrow \gamma^{-1}\). This is telling us that the magnitude of the ratio matters, but not whether protein or mRNA is more stable.

We can compute the phase boundary to be the line for which the inequality above becomes an equality. To compute this, we first need to find the fixed point \(x_0\) as a function of \(\beta\), \(\rho\), and \(n\) for the six-component system. This cannot be done analytically as in the case of the simplified three-component system we considered before. We therefore introduce a numerical technique to compute a fixed point in Technical Appendix 3c.

The linear stability diagram

Now that we can numerically compute the fixed point, we can proceed to construct the linear stability diagram showing the phase boundary (bifurcation line) for various n.

We write a function to compute the values of \(\gamma\) along the bifurcation line for various values of \(\beta\). It is useful to rewrite the expression for the bifurcation.

\begin{align} \gamma -\xi\sqrt{\gamma} + 1 = 0, \end{align}

where

\begin{align} \xi = \sqrt{\frac{3f_0^2}{4+2f_0}}. \end{align}

Then, we have

\begin{align} \gamma = \frac{1}{4}\left(\xi \pm \sqrt{\xi^2-4}\right)^2. \end{align}

Evidently, if \(\xi < 2\), there is no value of \(\gamma\) that satisfies this relation. There is no problem with this; it just says that there are regions in parameter space where oscillations cannot occur.

When we make a plot of the bifurcation, we will only show the plot for \(\gamma > 1\) because of the symmetry of \(\gamma\) and \(1/\gamma\) we have described.

For each bifurcation line, the oscillatory region is below and to the right of the boundary. We will start with \(\rho = 0\) and will vary \(n\). In the following plot, we see something interesting: when \(n>2\), the boundary approaches a vertical asymptote. This means that for strong enough promoters, one can get oscillations for any decay rates. By contrast when \(n<2\) the boundary approaches a horizontal asymptote. In this regime, it is essential the mRNA and protein decay rates must be similar enough to one another (\(\gamma\) must be below the line) or no oscillations are possible at any promoter strength. In between these regimes is the critical case of \(n=2\), where either larger \(\beta\) and \(\gamma\) closer to 1 both favor oscillations.

# Our beta range
beta = np.logspace(0, 5, 2000)

# Color according to n
colors = bokeh.palettes.Blues9[1:-1]

def fixed_point(beta, n, rho):
    """Function derived in technical appendix 9c for finding fixed point."""
    def _root_function(x):
        return beta - (x - beta * rho) * (1 + x ** n)

    return scipy.optimize.brentq(_root_function, 0, beta * (1 + rho))


def gamma_bifurcation(beta, n, rho):
    # Initialize gamma
    gamma = np.empty_like(beta)

    for i, b in enumerate(beta):
        x0 = fixed_point(b, n, rho)

        f0 = -b * n * x0 ** (n - 1) / (1 + x0 ** n) ** 2

        if f0 < -2:
            gamma[i] = np.nan
        else:
            xi = np.sqrt(3 * f0 ** 2 / (4 + 2 * f0))

            if xi < 2:
                gamma[i] = np.nan
            else:
                gamma[i] = (xi + np.sqrt(xi ** 2 - 4)) ** 2 / 4

    return gamma

p = bokeh.plotting.figure(
    width=400,
    height=300,
    x_axis_label="β",
    y_axis_label="γ",
    x_axis_type="log",
    y_axis_type="log",
    y_range=[1.3, 1e3],
    x_range=[1, 1e5],
)

# Make the plots
for n, color in zip([1.5, 1.75, 1.95, 2, 2.05, 2.5, 3][::-1], colors):
    gamma = gamma_bifurcation(beta, n, 0)
    p.line(beta, gamma, line_width=2, color=color, legend_label=f"n = {n}")

bokeh.io.show(p)

We can also see how leaky expression will affect the phase diagram. We will set \(n = 2\) and make a plot of the phase boundary for different values of \(\rho\). As you can see in this plot, leaky transcription kills the oscillations roughly when the amount of leaky protein production becomes sufficient to shut off expression of the next repressor.

p = bokeh.plotting.figure(
    width=400,
    height=300,
    x_axis_label="β",
    y_axis_label="γ",
    x_axis_type="log",
    y_axis_type="log",
    y_range=[1, 1e4],
    x_range=[1, 1e5],
)

for rho, color in zip([1e-2, 1e-3, 1e-4, 0], colors[::2]):
    gamma = gamma_bifurcation(beta, 2, rho)
    p.line(beta, gamma, line_width=2, color=color, legend_label=f"ρ = {rho}")

p.legend.location = "top_left"

bokeh.io.show(p)

The repressilator design objectives

Compared to the protein-only analysis, including the mRNAs explicitly provided additional insight. We learned that when the protein and mRNA decay times are comparable, the delay around the loop is effectively longer, reducing the minimum Hill coefficient and promoter strength required for oscillations. From this analysis we can extract several design objectives to maximize the chances of achieving self-sustaining oscillations in an experimental realization of the circuit:

1. Low “leakiness:” \(\rho \ll 1\) ⟶ Use tight, artificial promoters that can be fully repressed.

2. Strong promoters: \(\beta \gg 1\) ⟶ Can be achieved using strong promoters derived from phages that produce high protein levels

3. Similar protein and mRNA decay rates \(\gamma\approx 1\) ⟶ Destabilize repressors to increase their decay rates to be more comparable to those of mRNA. This can be done by adding destabilizing C-terminal tags based on the ssrA protein degradation system (Karzai, et al., 2000).

4. Ultrasensitive repression curves, ideally \(n > 1.5\) or \(2\), or as large as possible ⟶ Use intrinsically cooperative repression mechanisms, such as those from phage λ, or those that incorporate multiple binding sites, such as those in the TetR system. The phage λ \(P_R\) promoter architecture provides a high regulatory range, and can be adapted to work with binding sites for LacI and TetR (Lutz and Bujard, 1997).

In addition, there are also biological design rules as well:

5. To minimize toxicity from overexpressing repressors, put the circuit on a low copy plasmid (pSC101)…

6. …But to maximize the readout, put a fluorescent reporter gene on a higher copy number plasmid (ColE1).

7. Destabilize the fluorescent protein so that it can track the circuit activity

8. Avoid “read through” from one operon to the next ⟶ add transcriptional terminators between promoter-repressor units.

Based on these considerations, we designed the repressilator as a two plasmid system to be used in an E. coli strain deleted for the natural lac operon.

Here, the repressilator consists of three repressors on the low copy pSC101 plasmid, with TetR additionally repressing a green fluorescent protein reporter on the higher copy ColE1 plasmid. The lite suffix on the repressors signifies that they have a destruction tag to decrease their stability. The aav suffix on the GFP indicates that it incorporates a less active degradation tag variant.

Will it work?

It took a long time to design and then construct the repressilator. In addition to the molecular cloning, there were many steps to characterize the repressors and promoters and make sure signals could propagate through sequential repressor cascades.

After all of this, there was more than a little uncertainty as to whether the system would, indeed, exhibit self-sustaining oscillations.

To see how the circuit behaved, we used time-lapse imaging to record movies of individual repressilator-containing cells growing into microcolonies. (Lacking automated autofocus systems, one author slept near the microscope with an alarm clock to refocus the microscope every hour, all night long, exemplifying another important role for clocks in science and technology.)

In these movies, the changing fluorescence intensity in each cell provided a glimpse into the state of the oscillator over time:

This movie shows both clear oscillations in individual cells, as well as variability among cells in the amplitude, phase, and duration of each pulse. Analyzing these movies was done by manually tracking each cell backwards in time. This would now be done in a more automated fashion.

This procedure revealed clear oscillations in most cells, such as this:

Analysis of many cells showed a typical repressilator period of 160 ± 40 min (SD, n = 63), with a cell division time of ≈50-60 min at 30°C. Sibling cells desynchronized with one another over about two cell cycles (95 ± 10 min).

Evidently, a simple three element negative feedback loop can indeed generate oscillations. But those oscillations are variable. Do the dynamics have to be so variable? Next, we will see that oscillations can be improved dramatically, and even used as a timer to record bacterial growth in the mammalian gut.

Improving the repressilator

In Potvin-Trottier et al, Nature, 2016, Johan Paulsson and colleagues asked what accounted for the repressilator’s variability and whether it could be further improved (see Gao & Elowitz, Nature, 2016 for a summary of this work).

• Observation method. Instead of growth on agarose pads, where waste products can build up and influence cell growth and behavior, they switched to a microfluidic device developed by Suckjoon Jun termed “the mother machine,” (Wang, et al., 2010)which allows continuous observation of single cells over hundreds of generations by trapping it at the end of a channel (for more information, see Suckjoon Jun’s website). This revealed that, despite its variability, the original repressilator exhibited self-sustaining oscillations that never terminate.

• Integration of reporter. Now able to analyze the dynamics in more constant conditions, they found that much of this variation could be attributed to the reporter plasmid itself. Integrating the reporter into the repressilator plasmid reduced this variability, as seen in the movie and traces below. (In this movie, a cyan fluorescent protein, shown in the blue channel, is expressed at a constant level to allow image analysis, while the mVenus fluorescent protein, shown in the green channel, reports on the state of the repressilator.)

• The problem with TetR: It’s too good. A nice property of TetR is that it binds extremely tightly to its operator site. But this tightness of binding creates a problem as well, when one considers the discreteness of molecules and the stochasticity of their removal. For a DNA binding protein with a weaker affinity for its site (higher \(K\)), de-repression occurs at a relatively high concentration of the repressor, where discreteness can be ignored, and concentration dynamics are approximately continuous, as we have assumed. But when binding is very tight, the timing of de-repression by TetR depends sensitively on when the final molecules of TetR are degraded or diluted from the cell, as shown in the following figure from Potvin-Trottier et al. This leads to variability in the overall period. The authors showed that this effect can be mitigated by inclusion of a “DNA sponge”–a plasmid containing extra TetR binding sites. The presence of these extra, competing sites, pushes the effective threshold for de-repression to higher concentrations. (In the original repressilator design the reporter plasmid fortuitously played a similar role.) We will analyze stochasticity in more detailed in coming chapters.

• Eliminating enzymatic degradation makes the repressilator an astonishingly precise clock. Finally, an additional source of variability stemmed from fluctuations in the enzymatic protein degradation machinery. This source of variation could be circumvented by reducing the amounts of repressor made, or more dramatically by eliminating their degradation altogether, allowing dilution to dominate protein removal. Without protein degradation, oscillation periods stretched out to as much as 14 cell cycles. However, their precision increased dramatically, with a standard deviation in period duration of only about 14%, and an incredibly regular pulse shape and amplitude, as you can see in the image below. In fact, the repressilator became so precise that it would take on average 180 cell cycles to accumulate even a half-period of phase drift!

Here, with three colors to allow simultaneous observation of all three repressors, is one of the final repressilator designs as a movie and a typical trace.

In fact, this is so accurate that you can see it in a test tube. Here, the cells have no means of synchronization, yet they stay in phase for multiple oscillations.

The dilution-based repressilator circuit is so precise that even though cells are not synchronized with one another, they stay in sync over many generations of growth. When the cells grow as a colony on a plate, most growth occurs at the leading edge of the colony. Behind that front, cells stop growing, effectively leaving a record of the state of their repressilators at each point in time, somewhat like tree rings, as shown in this image (red, green, and blue channels show fluorescence from three reporters, one expressed with each repressor):

These improvements are summarized in this image:

This is a great example of a “less is more” principle of circuit design. One might have imagined that increasing the precision of the repressilator would require additional regulatory circuitry, but this does not seem to be the case, at least under these well-controlled conditions. As we wrote in an accompanying piece, “Evidently, precision does not necessarily demand circuit complexity and, in this case, even seems to benefit from minimalism.” (Gao and Elowitz, 2016) At the same time, achieving a constant period irrespective of growth conditions will likely require compensation mechanisms, much as early navigational clocks required additional systems to compensate for the effects of varying temperature on clock period (Love, 2016).

Precision repressilators can record bacterial growth in the mammalian gut

Riglar and coworkers took advantage of the incredible precision offered by this circuit to interrogate bacterial growth dynamics in the mammalian gut, which we describe in this section. (Riglar, et al., 2019).

The approach of Riglar and coworkers relied on the ability to read out repressilator phase from colony ring phases. More specifically, they found that the phase of the repressilator in the initial cell that seeded the colony controls the radial phasing of the rings in the colony. Cells that are plated at the peak of reporter fluorescence produce colonies with a high intensity “dot” at the center, whereas cells plated at the trough of the oscillations have rings shifted by 180 degrees. This is shown schematically in part c of the image below. They also took advantage of the ability to modulate the phase of the repressilator directly with two small molecule drugs, anhydrotetracycline (aTc) and IPTG, which respectively inhibit TetR and LacI proteins (a, below). After using these drugs to synchronize a population of cells at one phase, they could observe a beautiful linear relationship between the repressilator phase (as determined by the phasing of colony rings) and the elapsed growth of the cells (bottom-right plot, below). Further, they verified that it was really growth and not time per se that correlated with phase (see Figure 3 of the paper).

Colony phasing indicates cumulative cell growth. (a) The repressilator can be perturbed by adding the drugs aTc and IPTG, which inhibit TetR and LacI, respectively. This allows one to put the system into a defined phase. (b) Colony showing the rings of fluorescent protein expression that indicate phase. (c) Schematic indicating how repressilator phase correlates with colony ring phase. (lower left) Examples of colonies generated by cells exposed to either aTc or IPTG. Note the distinct phasing. (lower right) Colony ring phase correlates with number of generations of growth.

Taken together, this means that the phase of the oscillator can be used as a measure of bacterial growth.

They wondered whether this might allow one to recover the history of cell growth in an environment that would otherwise be inaccessible. E. coli normally spend part of their life cycle within the mammalian gut, an environment that is difficult to observe directly. Ideally, one would like ot know how much bacterial cells typically grow within the gut, how much cell-cell variability there is in that growth, and how the amount of growth depends on conditions in the gut, including which other bacteria are present, or whether it has been recently depleted of bacteria by antibiotics. Using the repressilator as a growth recorder, they identified conditions, such as inflammation, that led to greater variability in growth. They also saw that growth varied between different spatial areas or niches in the gut.

Conclusions

The process of designing and building a synthetic oscillator in bacteria taught us many things about oscillatory dynamics, and synthetic circuit design:

• Limit cycles are stable orbits in phase space. Circuits that generate limit cycle dynamics are ideal for self-sustaining oscillations in living cells.

• The repressilator uses a cycle of three repressors to create a delayed negative feedback loop.

• Linear stability analysis allows us to determine the stability of a fixed point.

• The simplified protein-only repressilator model has a single fixed point, and generates limit cycle oscillations when this point becomes unstable.

• High Hill coefficients and strong promoters favor oscillations.

• Including mRNA dynamics in the model revealed that comparable mRNA and protein decay rates favors oscillations and allows them to occur with less ultrasensitivity. This suggests destabilizing proteins to make their decay rates more similar to those of mRNA.

• A repressilator designed to meet these conditions shows sustained oscillations in individual E. coli cells, with significant variability.

Improving the repressilator revealed additional conclusions:

• Even a relatively simple synthetic circuit of three genes can operate with incredible precision in a living cell.

• The three-repressor structure of the repressilator, and the way its dynamics progress sequentially from expression of one repressor to the next, means that one can define the phase of the oscillator. Phase is a critical aspect of many natural biological oscillators. For example, the cell cycle is controlled by an oscillator that advances a cell from one phase to the next, in a sequential cycle, with each phase performing unique activities. These phases begin with growth during G1 phase, followed by DNA replication during S phase, a second G phase (G2), and finally an M phase corresponding to mitosis.

• The accurate “2.0” repressilator can provide insight into bacterial growth dynamics within the gut.

One could imagine further exploring how engineered probiotic cells behave when introduced into animals.

Taken together, this work shows the complementary power of thoughtful rational design and engineering in creating dynamic cellular circuits. Can our own rationally designed circuits eventually compete those produced by evolution, the “blind watchmaker” itself? Only time will tell.

---

References

• Czeisler, C. A., et al., Stability, precision, and near-24-hour period of the human circadian pacemaker, Science, 284, 2177–2181, 1999. (link)

• Dawkins, R., The Blind Watchmaker, W. W. Norton & Company, 1986 (link)

• Elowitz, M. B. and Leibler, S., A synthetic oscillatory network of transcriptional regulators, Nature, 403, 335–338, 2000. (link)

• Elowitz, M. B., Transport, assembly, and dynamics in systems of interacting proteins, Ph.D. thesis, Princeton University, 1999. (link)

• Gao, X. J. and Elowitz, M. B., Precision timing in a cell, Nature, 538, 462–463, 2016. (link)

• Karzai, A. W., Roche, E. D., Sauer, R. T., The SsrA–SmpB system for protein tagging, directed degradation and ribosome rescue, Nat. Struct. Biol., 7, 449–455, 2000. (link)

• Konopka, R. J. and Benzer, S., Clock mutants of Drosphila melanogaster, _Proc. Natl. Acad. Sci. USA, 68, 2112–2116, 1971. (link)

• Love, S., Building an impossible clock, The Atlantic, January 19, 2016. (link)

• Lutz, R. and Bujard, H., Independent and tight regulation of transcriptional units in Escherichia coli via the LacR/O, the TetR/O and AraC/I₁-I₂ regulatory elements, Nucl. Acid Res., 25, 1203–1210, 1997. (link)

• Mihalcescu, I., Hsing, W., Leiber, S., Resilient circadian oscillator revealed in individual cyanobacteria, Nature, 430, 81–85, 2004. (link)

• Nakajima, et al., Reconstitution of circadian oscillation of cyanobacterial KaiC phosphorylation in vitro, Science, 308, 414–415, 2005.(link)

• Potvin-Trottier, L., et al., Synchronous long-term oscillations in a synthetic gene circuit, Nature, 538, 514–517, 2016 (link)

• Riglar, D. T., Bacterial variability in the mammalian gut captured by a single-cell synthetic oscillator, Nature Comm., 10, 4665, 2019. (link)

• Rust, M. J., et al., Ordered phosphorylation governs oscillation of a three-protein circadian clock, Science, 318, 809–812, 2007. (link)

• Strogatz, S. H., Nonlinear Dynamics and Chaos With Applications to Physics, Biology, Chemistry, and Engineering, 2nd Ed., CRC Press, 2015. (link)

• Wang, P., et al., Robust growth of Escherichia coli, Curr. Biol., 20, 1099–1103, 2010. (link)

• Welsh, D. K., et al., Individual neurons dissociated from rat suprachiasmatic nucleus express independently phased circadian firing rhythms, Neuron, 14, 697–706, 1995. (link)

• Welsh, D. K., Bioluminescence imaging of individual fibroblasts reveals persistent, independently phased circadian rhythms of clock gene expression, Curr. Biol., 14, 2289–2295, 2004. (link)

---

Technical appendices

• 9a. Fixed points and composite functions
• 9b. Linear stability analysis
• 9c. Numerical one-dimensional bounded root finding

---

Problems

• 9.1: Coupled repressilators
• 9.2: The KaiABC clock
• 9.3: Linear stability analysis of the repressilator with mRNA

Computing environment

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

Python implementation: CPython
Python version       : 3.10.10
IPython version      : 8.12.0

numpy      : 1.23.5
scipy      : 1.10.1
bokeh      : 3.1.0
matplotlib : 3.7.1
colorcet   : 3.0.1
biocircuits: 0.1.11
jupyterlab : 3.5.3