Biocircuit 101: ODE Models of Biocircuits
Click code cells to expand/collapse.
Biocircuit 101: ODE Models of Biocircuits
For this tutorial, we do not care about the molecular details (for the sake of time), ODE solving techniques, and visualization coding.
AI could help with coding and solving ODEs. In my opinion, it is more important to know:
Workflow for Each Model
1. Describe the circuit: write ODE for each variable in the circuit (nondimensionalization)
2. Solve for the dynamics: show how all variables evolve over time
- use numerical methods if it is hard to solve explicitly
- tune parameters and interpret the result
3. Do Nullcline Analysis: find fixed points/steady state locus
- tune parameters and interpret the result along with its dynamics
4. Do Bifurcation Analysis: track the stability of the circuit in the parameter space
- determine the parameter values to realize the function you want
0. Environment Setup
Set up a Python environment with all the packages we need.Open your terminal and cd to the folder containing environment.yml
conda env create -f environment.yml
Activate it:
conda activate biocircuits
jupyter lab
Run the cell below to confirm all packages are available:
▼[4] Code
# ============================================================
# IMPORTS -- libraries used throughout the notebook
# ============================================================
import numpy as np # numerical arrays & math
from bokeh.plotting import figure, show # Bokeh: create & display plots
from bokeh.layouts import column, row # Bokeh: arrange plots & widgets
from bokeh.models import ( # Bokeh: interactive components
ColumnDataSource, # -- efficient data <-> plot bridge
CustomJS, # -- browser-side callbacks (no Python kernel)
Slider, # -- interactive parameter control
Span, # -- vertical/horizontal marker line
VArea, # -- filled region (bifurcation diagrams)
)
from bokeh.io import output_notebook # Bokeh: render plots inside Jupyter
from scipy.optimize import fsolve # root-finding (fixed-point analysis)
from scipy.integrate import odeint # ODE numerical integration (LSODA adaptive solver)
import plotly.graph_objects as go # Plotly: 3D interactive plots (repressilator)
# Enable interactive Bokeh plots in the notebook
output_notebook()
1. Get the taste of ODE description
The simplest biocircuit consists of a gene that is constitutively expressed (produced at a constant rate) and degraded/diluted.
$$dx/dt = \beta - \gamma x$$
dimensionless form: $$dx/dt = \beta - x$$
It is a linear first-order ODE and can be solved analytically. So no numerical method is needed.
Using separation of variables (or the integrating factor method), with initial condition $x(0) = x_0$:
$$x(t) = \beta + (x_0 - \beta)e^{-t}$$
If we start from $x_0 = 0$ (no protein initially), this simplifies to:
$$x(t) = \beta\left(1 - e^{-t}\right)$$
> Note: For more complex circuits (autoregulation, toggle switch, repressilator), analytical solutions are generally not available, and we must resort to numerical integration. But for this simplest case, the analytical solution gives us everything we need. >
▼[7] Code
# ============================================================
# SIMPLE PRODUCTION-DEGRADATION DEVICE
# ODE: dx/dt = beta - x (analytical solution: x(t) = beta*(1 - e^(-t)))
# ============================================================
def plot_simple_device():
# --- Parameters ---
beta_init = 5.0 # production rate beta
t = np.linspace(0, 10, 100) # time vector for dynamics plot
x_nc = np.linspace(0, 10, 100) # concentration range for nullcline plot
# --- Analytical solution ---
# dx/dt = beta - x --> x(t) = beta*(1 - e^(-t)) with x(0)=0
x_dyn_init = beta_init * (1 - np.exp(-t))
# Production term is constant (= beta); removal term is linear (= x)
prod_init = np.full_like(x_nc, beta_init)
# --- Bokeh data sources ---
src_dyn = ColumnDataSource(data=dict(t=t, x=x_dyn_init))
src_nc = ColumnDataSource(data=dict(x=x_nc, prod=prod_init, deg=x_nc))
# Fixed point: intersection of production and removal --> x* = beta
src_pt = ColumnDataSource(data=dict(x=[beta_init], y=[beta_init]))
# --- Dynamics plot (x vs t) ---
p_dyn = figure(title="Dynamics", x_axis_label="t", y_axis_label="x", width=400, height=300)
p_dyn.line('t', 'x', source=src_dyn, line_width=3, color='#2ca02c')
# --- Nullcline plot (production & removal rates vs x) ---
p_nc = figure(title="Nullcline Analysis", x_axis_label="x", y_axis_label="rate", width=400, height=300)
p_nc.line('x', 'prod', source=src_nc, line_width=3, color='#1f77b4', legend_label="Production")
p_nc.line('x', 'deg', source=src_nc, line_width=3, color='#ff7f0e', legend_label="Removal")
p_nc.scatter('x', 'y', source=src_pt, size=12, color='black', legend_label="Stable Point")
p_nc.legend.location = "top_left"
p_nc.legend.label_text_font_size = "8pt"
# --- Interactive slider for beta ---
slider_beta = Slider(start=0.0, end=10.0, value=beta_init, step=0.1, title="beta")
# Browser-side callback: update all plots when beta changes (uses analytical formula)
callback = CustomJS(args=dict(sd=src_dyn, snc=src_nc, spt=src_pt, b=slider_beta), code="""
const beta = b.value;
const t = sd.data['t'];
const x_d = sd.data['x'];
for(let i=0; i<t.length; i++) {
x_d[i] = beta * (1 - Math.exp(-t[i]));
}
sd.change.emit();
const x_n = snc.data['x'];
const prod = snc.data['prod'];
for(let i=0; i<x_n.length; i++) prod[i] = beta;
snc.change.emit();
spt.data = {x: [beta], y: [beta]};
spt.change.emit();
""")
slider_beta.js_on_change('value', callback)
# Layout: slider on top, dynamics (left) and nullcline (right) side by side
return column(slider_beta, row(p_dyn, p_nc))
show(plot_simple_device())
Interactive: Dynamics & Nullcline
Varying β gives a straight line of steady states. The system is always monostable. No need for bifurcation anaylsis.
▼[9] Code
# ============================================================
def plot_simple_device_advanced(t_switch=5.0, beta1=2.0, beta2=8.0, x0=0.0):
"""
Demonstrate memory-like behavior in a simple production-degradation system.
The system starts with production rate beta1, then switches to beta2 at t_switch.
This shows how the system "remembers" its history and adapts to new parameters.
Parameters:
t_switch — time to switch from beta1 to beta2 (default: 5.0)
beta1 — initial production rate (default: 2.0)
beta2 — production rate after switch (default: 8.0)
x0 — initial concentration (default: 0.0)
"""
from scipy.integrate import odeint
# Total simulation time and time vector
t_total = 10.0
t = np.linspace(0, t_total, 500)
# Define the ODE: dx/dt = beta - x (production - degradation)
def rhs(x, t, beta):
return beta - x
# --- Solve ODE with a parameter switch using scipy.integrate.odeint ---
# Split time vector at the switch point
t1 = t[t <= t_switch]
t2 = t[t >= t_switch]
# Phase 1: integrate from t=0 to t=t_switch with production rate beta1
sol1 = odeint(rhs, x0, t1, args=(beta1,))
# Phase 2: integrate from t=t_switch to t=t_total with production rate beta2.
# The initial condition for phase 2 is the final state of phase 1, so the
# trajectory is continuous across the switch.
x_switch = sol1[-1, 0]
sol2 = odeint(rhs, x_switch, t2, args=(beta2,))
# Concatenate: if t_switch falls exactly on a grid point, t1 and t2 overlap (drop duplicate);
# otherwise they partition t (keep all points)
if t1[-1] == t2[0]:
sol = np.concatenate([sol1[:-1], sol2]).flatten()
else:
sol = np.concatenate([sol1, sol2]).flatten()
# --- Build the interactive Bokeh plot ---
src = ColumnDataSource(data=dict(t=t, x=sol))
# Vertical dashed line marking the parameter switch time
switch_line = Span(location=t_switch, dimension="height",
line_color="red", line_dash="dashed", line_width=2)
p = figure(title="Simple Device: Parameter Switch Dynamics",
x_axis_label="Time (t)", y_axis_label="Concentration (x)",
width=600, height=350)
p.line('t', 'x', source=src, line_width=3, color='#2ca02c')
p.add_layout(switch_line)
# --- Interactive sliders ---
s_t = Slider(start=1.0, end=9.0, value=t_switch, step=0.5, title="Switch Time")
s_b1 = Slider(start=0.0, end=10.0, value=beta1, step=0.1, title="β₁ (before switch)")
s_b2 = Slider(start=0.0, end=10.0, value=beta2, step=0.1, title="β₂ (after switch)")
s_x0 = Slider(start=0.0, end=10.0, value=x0, step=0.1, title="x₀ (initial concentration)")
# Browser-side (JavaScript) solver using RK4 for smooth, accurate interactivity
code = """
const t_switch = s_t.value, beta1 = s_b1.value, beta2 = s_b2.value, x0 = s_x0.value;
const t = src.data['t'], x = src.data['x'];
let x_curr = x0;
for(let i=0; i<t.length; i++) {
x[i] = x_curr;
if(i < t.length - 1) {
let dt = t[i+1] - t[i];
let beta = (t[i] < t_switch) ? beta1 : beta2;
// RK4 integration step
let k1 = beta - x_curr;
let k2 = beta - (x_curr + 0.5 * dt * k1);
let k3 = beta - (x_curr + 0.5 * dt * k2);
let k4 = beta - (x_curr + dt * k3);
x_curr += (dt / 6.0) * (k1 + 2*k2 + 2*k3 + k4);
}
}
src.change.emit();
switch_line.location = t_switch;
"""
cb = CustomJS(args=dict(src=src, s_t=s_t, s_b1=s_b1, s_b2=s_b2, s_x0=s_x0, switch_line=switch_line), code=code)
s_t.js_on_change("value", cb)
s_b1.js_on_change("value", cb)
s_b2.js_on_change("value", cb)
s_x0.js_on_change("value", cb)
# Layout: plot on the left, sliders stacked on the right
return row(p, column(s_t, s_b1, s_b2, s_x0))
show(plot_simple_device_advanced())Interactive: Parameter Switch
2. Autoregulation (Positive Feedback)
Dimensionless ODE: $$\frac{\mathrm{d}x}{\mathrm{d}t} = \beta \frac{x^n}{1 + x^n} - x$$
Dynamics and nullcline anaylsis
▼[12] Code
def plot_autoactivator(x0_low=0.5, x0_high=3):
# ============================================================
# PART A: CORE MATH -- ODE definition & numerical integration
# (this is what you want to focus on for learning)
# ============================================================
#
# Autoactivator dimensionless ODE:
# dx/dt = beta * x^n / (1 + x^n) - x
# +---- production ----+ +- degradation -+
#
# Parameters:
# beta -- maximal expression rate
# n -- Hill coefficient (cooperativity)
#
# User-configurable initial conditions (for dynamics):
# x0_low -- starting value for the "Starts Low" trajectory (default: 0.5)
# x0_high -- starting value for the "Starts High" trajectory (default: 6.0)
#
# Tip: To see bistability clearly, choose x0_low below the unstable
# fixed point (e.g., 0.01) and x0_high above it (e.g., 6.0).
x_max, beta_val, n_val = 8.0, 5.0, 4.0
x_arr = np.linspace(0, x_max, 500)
t_arr = np.linspace(0, 15, 500)
# Production term (for nullcline plotting)
prod_init = beta_val * (x_arr**n_val) / (1.0 + x_arr**n_val)
# --- Fixed-point analysis ---
# Fixed-point condition: dx/dt = 0 --> beta*x^n/(1+x^n) = x
# Stability condition: |f'(x*)| < 1 is stable, > 1 is unstable
# Note: f'(0) = beta*n*x^(n-1) for x->0
# - n < 1: f'(0) -> infinity (unstable)
# - n = 1: f'(0) = beta
# - n > 1: f'(0) = 0 (stable)
#
# Use fsolve with many initial guesses to find ALL fixed points
# (grid-based sign-change search misses points near x=0 when n~1)
def f_prime(xv, beta, n):
"""Derivative of the ODE right-hand side at xv (for stability check)"""
if xv == 0:
if n < 1: return float('inf')
return beta if n == 1 else 0.0
return (beta * n * (xv**(n - 1))) / ((1.0 + xv**n)**2)
def find_fixed_points(beta, n):
"""Find all fixed points using fsolve with many initial guesses"""
f = lambda x: beta * x**n / (1 + x**n) - x
guesses = np.linspace(0.001, x_max - 0.01, 300)
roots = []
for g in guesses:
try:
r = fsolve(f, g, full_output=False)[0]
if 0.001 <= r <= x_max and abs(f(r)) < 1e-6:
if not any(abs(r - existing) < 1e-4 for existing in roots):
roots.append(r)
except:
pass
return sorted(roots)
nonzero_roots = find_fixed_points(beta_val, n_val)
st_x, st_y, ust_x, ust_y = [], [], [], []
# Classify x = 0
if f_prime(0, beta_val, n_val) < 1.0:
st_x.append(0.0); st_y.append(0.0)
else:
ust_x.append(0.0); ust_y.append(0.0)
# Classify non-zero fixed points
for r in nonzero_roots:
if f_prime(r, beta_val, n_val) < 1.0:
st_x.append(r); st_y.append(r)
else:
ust_x.append(r); ust_y.append(r)
# --- Numerical ODE integration (CORE) ---
# Using scipy.integrate.odeint (LSODA adaptive solver)
# - Auto-detects stiffness, switches between Adams and BDF methods
# - Adaptive step size + error control (default rtol=1e-8, atol=1e-12)
# - Far superior to fixed-step Forward Euler
#
# Signature: odeint(func, y0, t, args=...)
# func(y, t, *args) returns dy/dt
def autoactivator_rhs(x, t, beta, n):
"""ODE right-hand side: dx/dt = beta * x^n / (1 + x^n) - x"""
return beta * (x**n) / (1.0 + x**n) - x
# Integrate from two different initial conditions to show bistability
x1_init = odeint(autoactivator_rhs, x0_low, t_arr, args=(beta_val, n_val)).flatten()
x2_init = odeint(autoactivator_rhs, x0_high, t_arr, args=(beta_val, n_val)).flatten()
# ============================================================
# PART B: VISUALIZATION (Bokeh plotting + interactivity, math-independent)
# ============================================================
src_nc = ColumnDataSource(data=dict(x=x_arr, prod=prod_init, deg=x_arr))
src_dyn = ColumnDataSource(data=dict(t=t_arr, x1=x1_init, x2=x2_init))
src_st = ColumnDataSource(data=dict(x=st_x, y=st_y))
src_unst = ColumnDataSource(data=dict(x=ust_x, y=ust_y))
p_dyn = figure(title="Dynamics (Two different initial conditions)", x_axis_label="t", y_axis_label="x", width=400, height=350)
p_dyn.line('t', 'x1', source=src_dyn, line_width=2, color='purple', legend_label=f"Starts at {x0_low}")
p_dyn.line('t', 'x2', source=src_dyn, line_width=2, color='orange', legend_label=f"Starts at {x0_high}")
p_dyn.legend.location = "top_left"
p_dyn.legend.label_text_font_size = "8pt"
p_nc = figure(title="Nullclines", x_axis_label="x", y_axis_label="rate", width=400, height=350)
p_nc.line('x', 'prod', source=src_nc, line_width=3, color='#1f77b4', legend_label="Production")
p_nc.line('x', 'deg', source=src_nc, line_width=3, color='#ff7f0e', legend_label="Removal")
p_nc.scatter('x', 'y', source=src_st, size=10, color='black', legend_label="Stable")
p_nc.scatter('x', 'y', source=src_unst, size=10, fill_color='white', line_color='black', line_width=2, legend_label="Unstable")
p_nc.legend.location = "top_left"
p_nc.legend.label_text_font_size = "8pt"
s_beta = Slider(start=0.0, end=10.0, value=beta_val, step=0.1, title="beta")
s_n = Slider(start=1, end=8.0, value=n_val, step=0.1, title="n")
# --- CustomJS: browser-side real-time update (no Python kernel needed) ---
# Uses Newton's method + bisection to find ALL fixed points (stable + unstable).
# Newton's method converges quickly; bisection is used as a robust fallback
# to find the unstable middle fixed point between two stable ones.
# This approach is adapted from the biocircuits library (jsplots.py).
code = """
const beta = s_beta.value, n = s_n.value;
const x = snc.data['x'], prod = snc.data['prod'], deg = snc.data['deg'];
// Update production term for nullcline plot
for(let i=0; i<x.length; i++) {
prod[i] = beta * Math.pow(x[i], n) / (1 + Math.pow(x[i], n));
}
snc.change.emit();
// --- Root-finding helpers (Newton's method + bisection) ---
// Fixed-point condition: beta*x^n/(1+x^n) - x = 0
const diff = (xv) => (beta * Math.pow(xv, n) / (1 + Math.pow(xv, n))) - xv;
// Derivative of diff: beta*n*x^(n-1)/(1+x^n)^2 - 1
const ddiff = (xv) => {
if (xv === 0) {
if (n < 1) return Infinity;
return (n === 1) ? beta - 1 : -1;
}
const xn = Math.pow(xv, n);
return (beta * n * Math.pow(xv, n - 1)) / Math.pow(1 + xn, 2) - 1.0;
};
// Stability: |f'(x*)| < 1 is stable, > 1 is unstable
const f_prime = (xv) => {
if (xv === 0) {
if (n < 1) return Infinity;
return (n === 1) ? beta : 0;
}
const xn = Math.pow(xv, n);
return (beta * n * Math.pow(xv, n - 1)) / Math.pow(1 + xn, 2);
};
// Newton's method: fast convergence to a root near the starting guess
function newtonSolve(x0, tol, maxIter) {
tol = tol || 1e-8;
maxIter = maxIter || 100;
let x = x0;
for (let i = 0; i < maxIter; i++) {
let y = diff(x);
let yp = ddiff(x);
if (Math.abs(yp) < 1e-14) break;
let xNew = x - y / yp;
if (xNew < 0) xNew = x / 2; // keep positive
if (Math.abs(xNew - x) <= tol) return xNew;
x = xNew;
}
return x;
}
// Bisection method: guaranteed to converge for bracketed roots
function bisectionSolve(lower, upper, tol, maxIter) {
tol = tol || 1e-8;
maxIter = maxIter || 100;
let a = lower, b = upper;
let fa = diff(a), fb = diff(b);
if (fa * fb > 0) return null;
for (let i = 0; i < maxIter; i++) {
let mid = (a + b) / 2.0;
let fmid = diff(mid);
if (Math.abs(fmid) < tol || (b - a) / 2.0 < tol) return mid;
if (fa * fmid < 0) { b = mid; fb = fmid; }
else { a = mid; fa = fmid; }
}
return (a + b) / 2.0;
}
// Find ALL fixed points using Newton + bisection
let st_x = [], st_y = [], ust_x = [], ust_y = [];
const x_max = 8.0;
// Newton from large x -> finds the rightmost (stable) positive fixed point
let xStable = newtonSolve(x_max);
// Collect positive fixed points (verify diff ≈ 0 to discard false convergence)
let roots = [];
if (xStable !== null && xStable > 0.001 && Math.abs(diff(xStable)) < 1e-4) {
roots.push(xStable);
}
// Use bisection to find the unstable middle point (if it exists).
// diff(x) < 0 just after x=0. The unstable point is where diff changes
// sign. We check two sub-intervals: [0.01, xStable/2] and [xStable/2, xStable].
if (roots.length === 1) {
let xMid = roots[0] / 2;
let fMid = diff(xMid);
let xUnstable = null;
if (diff(0.01) * fMid < 0) {
xUnstable = bisectionSolve(0.01, xMid);
} else if (fMid * diff(roots[0] - 1e-4) < 0) {
xUnstable = bisectionSolve(xMid, roots[0] - 1e-4);
}
if (xUnstable !== null && xUnstable > 0.01 && Math.abs(xUnstable - roots[0]) > 0.1) {
roots.push(xUnstable);
roots.sort((a, b) => a - b);
}
}
// Classify x = 0 (fixed point only if n >= 1)
if (n >= 1) {
if (f_prime(0) < 1.0) { st_x.push(0.0); st_y.push(0.0); }
else { ust_x.push(0.0); ust_y.push(0.0); }
}
// Classify non-zero fixed points
for (let r of roots) {
if (f_prime(r) < 1.0) { st_x.push(r); st_y.push(r); }
else { ust_x.push(r); ust_y.push(r); }
}
sst.data = {x: st_x, y: st_y}; sst.change.emit();
sust.data = {x: ust_x, y: ust_y}; sust.change.emit();
// Browser-side ODE integration via RK4 (4th-order accuracy, far better than Euler)
const t = sdyn.data['t'], x1 = sdyn.data['x1'], x2 = sdyn.data['x2'];
const dt = t[1] - t[0];
let curr_x1 = """ + str(x0_low) + """, curr_x2 = """ + str(x0_high) + """;
for(let i=0; i<t.length; i++) {
x1[i] = curr_x1; x2[i] = curr_x2;
let k1a = diff(curr_x1);
let k2a = diff(curr_x1 + 0.5*dt*k1a);
let k3a = diff(curr_x1 + 0.5*dt*k2a);
let k4a = diff(curr_x1 + dt*k3a);
curr_x1 += (dt/6.0) * (k1a + 2*k2a + 2*k3a + k4a);
let k1b = diff(curr_x2);
let k2b = diff(curr_x2 + 0.5*dt*k1b);
let k3b = diff(curr_x2 + 0.5*dt*k2b);
let k4b = diff(curr_x2 + dt*k3b);
curr_x2 += (dt/6.0) * (k1b + 2*k2b + 2*k3b + k4b);
}
sdyn.change.emit();
"""
cb = CustomJS(args=dict(snc=src_nc, sdyn=src_dyn, sst=src_st, sust=src_unst, s_beta=s_beta, s_n=s_n), code=code)
s_beta.js_on_change('value', cb)
s_n.js_on_change('value', cb)
return column(row(s_beta, s_n), row(p_dyn, p_nc))
# Example usage:
# plot_autoactivator() # uses defaults: x0_low=0.5, x0_high=6.0
# plot_autoactivator(x0_low=0.01, x0_high=6.0) # low start below unstable point
show(plot_autoactivator())
Interactive: Dynamics & Nullcline
Stable Unstable
> Note: The nullcline plot may show two stable fixed points (bistability), while the dynamics plot shows both trajectories converging to the same state. This is not a contradiction! Tt reflects that both initial conditions lie in the basin of attraction of the same stable state. The other stable state (e.g., $x=0$ for $n>1$) has a very small basin of attraction near the origin. To observe it, start a trajectory from a very small initial condition (e.g., $x_0 = 0.001$).
Memorize the state!
▼[15] Code
# ============================================================
# MEMORY DEMO: Autoactivator with TWO parameter switches
# ODE: dx/dt = beta * x^n / (1 + x^n) - x
# Uses scipy.integrate.odeint (three-phase integration)
# Demo: bistable -> monostable -> bistable (irreversible)
# ============================================================
def plot_autoactivator_advanced(t_switch1=5.0, t_switch2=10.0,
beta1=5.0, n1=4.0,
beta2=1.0, n2=4,
beta3=10.0, n3=4.0,
x0_low=0.5, x0_high=6.0):
"""
Demonstrate memory (bistability) in the autoactivator system with TWO parameter switches.
The system transitions through three parameter regimes:
- Phase 1 (t < t_switch1): (beta1, n1) -- typically bistable
- Phase 2 (t_switch1 < t < t_switch2): (beta2, n2) -- typically monostable
- Phase 3 (t > t_switch2): (beta3, n3) -- bistable again (but state may not recover)
This shows how a system can be driven from bistable -> monostable -> bistable,
and that the final state may depend on the path (hysteresis / irreversibility).
Parameters:
t_switch1 -- time of first parameter switch (default: 5.0)
t_switch2 -- time of second parameter switch (default: 10.0)
beta1, n1 -- parameters for phase 1 (default: 5.0, 4.0)
beta2, n2 -- parameters for phase 2 (default: 2.0, 1.5)
beta3, n3 -- parameters for phase 3 (default: 5.0, 4.0)
x0_low -- initial condition for low trajectory (default: 0.5)
x0_high -- initial condition for high trajectory (default: 6.0)
"""
# Total simulation time and time vector
t_total = 15.0
t = np.linspace(0, t_total, 1500)
# ODE right-hand side: dx/dt = beta * x^n / (1 + x^n) - x
def rhs(x, t, beta, n):
return beta * x**n / (1 + x**n) - x
# --- Solve ODE with TWO parameter switches using scipy.integrate.odeint ---
# Three phases per trajectory; each phase's initial condition = previous phase's final state
t1 = t[t <= t_switch1]
t2 = t[(t > t_switch1) & (t <= t_switch2)]
t3 = t[t > t_switch2]
def solve_three_phase(x0, beta1, n1, beta2, n2, beta3, n3):
"""Integrate one trajectory through three parameter regimes."""
# Phase 1
sol1 = odeint(rhs, x0, t1, args=(beta1, n1))
# Phase 2
sol2 = odeint(rhs, sol1[-1, 0], t2, args=(beta2, n2))
# Phase 3
sol3 = odeint(rhs, sol2[-1, 0], t3, args=(beta3, n3))
# Concatenate: drop duplicate points at switch boundaries
parts = [sol1]
if len(t2) > 0:
parts.append(sol2[1:] if t1[-1] == t2[0] else sol2)
if len(t3) > 0:
parts.append(sol3[1:] if t2[-1] == t3[0] else sol3)
return np.concatenate(parts).flatten()
# Compute both trajectories (low and high initial conditions)
sol_high = solve_three_phase(x0_high, beta1, n1, beta2, n2, beta3, n3)
sol_low = solve_three_phase(x0_low, beta1, n1, beta2, n2, beta3, n3)
# --- Build the interactive Bokeh plot ---
src = ColumnDataSource(data=dict(t=t, x_high=sol_high, x_low=sol_low))
# Vertical dashed lines marking the two parameter switch times
switch_line1 = Span(location=t_switch1, dimension="height",
line_color="red", line_dash="dashed", line_width=2)
switch_line2 = Span(location=t_switch2, dimension="height",
line_color="blue", line_dash="dashed", line_width=2)
p = figure(title="Autoactivator: Bistable -> Monostable -> Bistable (Irreversible)",
x_axis_label="Time (t)", y_axis_label="Concentration (x)",
width=550, height=400)
p.line('t', 'x_high', source=src, line_width=2, color='purple', legend_label="Started High")
p.line('t', 'x_low', source=src, line_width=2, color='orange', legend_label="Started Low")
p.add_layout(switch_line1)
p.add_layout(switch_line2)
p.legend.location = "top_left"
p.legend.label_text_font_size = "8pt"
# --- Interactive sliders ---
s_t1 = Slider(start=1.0, end=7.0, value=t_switch1, step=0.5, title="Switch Time 1")
s_t2 = Slider(start=8.0, end=14.0, value=t_switch2, step=0.5, title="Switch Time 2")
s_b1 = Slider(start=0.0, end=10.0, value=beta1, step=0.5, title="beta_1 (phase 1)")
s_n1 = Slider(start=0.1, end=8.0, value=n1, step=0.1, title="n_1 (phase 1)")
s_b2 = Slider(start=0.0, end=10.0, value=beta2, step=0.5, title="beta_2 (phase 2)")
s_n2 = Slider(start=0.1, end=8.0, value=n2, step=0.1, title="n_2 (phase 2)")
s_b3 = Slider(start=1.0, end=10.0, value=beta3, step=0.5, title="beta_3 (phase 3)")
s_n3 = Slider(start=1.0, end=8.0, value=n3, step=0.1, title="n_3 (phase 3)")
s_x0l = Slider(start=0.0, end=2.0, value=x0_low, step=0.1, title="x0 (low)")
s_x0h = Slider(start=2.0, end=8.0, value=x0_high, step=0.1, title="x0 (high)")
# Browser-side solver using RK4 for smooth, accurate interactivity
code = """
const t_switch1 = s_t1.value, t_switch2 = s_t2.value;
const beta1 = s_b1.value, n1 = s_n1.value;
const beta2 = s_b2.value, n2 = s_n2.value;
const beta3 = s_b3.value, n3 = s_n3.value;
const x0_low = s_x0l.value, x0_high = s_x0h.value;
const t = src.data['t'], x_high = src.data['x_high'], x_low = src.data['x_low'];
const rhs = (x, beta, n) => beta * Math.pow(x, n) / (1 + Math.pow(x, n)) - x;
// RK4 integration for the high trajectory (three phases)
let x_curr = x0_high;
for(let i=0; i<t.length; i++) {
x_high[i] = x_curr;
if(i < t.length - 1) {
let dt = t[i+1] - t[i];
let beta = (t[i] < t_switch1) ? beta1 : (t[i] < t_switch2) ? beta2 : beta3;
let n = (t[i] < t_switch1) ? n1 : (t[i] < t_switch2) ? n2 : n3;
let k1 = rhs(x_curr, beta, n);
let k2 = rhs(x_curr + 0.5*dt*k1, beta, n);
let k3 = rhs(x_curr + 0.5*dt*k2, beta, n);
let k4 = rhs(x_curr + dt*k3, beta, n);
x_curr += (dt/6.0) * (k1 + 2*k2 + 2*k3 + k4);
}
}
// RK4 integration for the low trajectory (three phases)
x_curr = x0_low;
for(let i=0; i<t.length; i++) {
x_low[i] = x_curr;
if(i < t.length - 1) {
let dt = t[i+1] - t[i];
let beta = (t[i] < t_switch1) ? beta1 : (t[i] < t_switch2) ? beta2 : beta3;
let n = (t[i] < t_switch1) ? n1 : (t[i] < t_switch2) ? n2 : n3;
let k1 = rhs(x_curr, beta, n);
let k2 = rhs(x_curr + 0.5*dt*k1, beta, n);
let k3 = rhs(x_curr + 0.5*dt*k2, beta, n);
let k4 = rhs(x_curr + dt*k3, beta, n);
x_curr += (dt/6.0) * (k1 + 2*k2 + 2*k3 + k4);
}
}
src.change.emit();
switch_line1.location = t_switch1;
switch_line2.location = t_switch2;
"""
cb = CustomJS(args=dict(src=src, s_t1=s_t1, s_t2=s_t2,
s_b1=s_b1, s_n1=s_n1, s_b2=s_b2, s_n2=s_n2,
s_b3=s_b3, s_n3=s_n3,
s_x0l=s_x0l, s_x0h=s_x0h,
switch_line1=switch_line1, switch_line2=switch_line2), code=code)
s_t1.js_on_change("value", cb)
s_t2.js_on_change("value", cb)
s_b1.js_on_change("value", cb)
s_n1.js_on_change("value", cb)
s_b2.js_on_change("value", cb)
s_n2.js_on_change("value", cb)
s_b3.js_on_change("value", cb)
s_n3.js_on_change("value", cb)
s_x0l.js_on_change("value", cb)
s_x0h.js_on_change("value", cb)
# Layout: plot on left, sliders stacked on right
sliders = column(s_t1, s_t2, s_b1, s_n1, s_b2, s_n2, s_b3, s_n3, s_x0l, s_x0h)
return row(p, sliders)
show(plot_autoactivator_advanced())
Interactive: Two Parameter Switches
Bifurcation Analysis (Autoactivator)
By mapping out the stable and unstable regimes across parameters, we generate a Bifurcation Diagram. The boundary between monostability (1 stable point) and bistability (2 stable points) can be analytically derived by setting $\frac{\mathrm{d}x}{\mathrm{d}t} = 0$ and $\frac{\mathrm{d}^2x}{\mathrm{d}t^2} = 0$.The boundary curve is given by: $$\beta = \frac{n}{(n-1)^{\frac{n-1}{n}}}$$ > Note: This bifurcation boundary above is computed from an analytical formula derived from the fixed-point condition. For more complex models where an analytical boundary cannot be derived, a numerical method can be used instead.
▼[17] Code
# ============================================================
# BIFURCATION ANALYSIS: Autoactivator
# Analytical boundary: beta = n / ((n-1)^((n-1)/n))
# Above the curve: bistable; below: monostable
# ============================================================
def plot_autoactivator_bifurcation():
# Range of Hill coefficients to sweep
n_vals = np.linspace(1.001, 10, 500)
# Analytical bistability boundary derived from fixed-point conditions
beta_boundary = n_vals / ((n_vals - 1)**((n_vals - 1)/n_vals))
# --- Bifurcation diagram ---
p = figure(title="Bifurcation Diagram: Autoactivator",
x_axis_label="n", y_axis_label="nondimensional beta",
x_range=(1, 10), y_range=(0, 3), width=600, height=400)
# Italic axis labels for math notation
p.xaxis.axis_label_text_font_style = "italic"
p.yaxis.axis_label_text_font_style = "italic"
# Gray region = bistable (above boundary)
p.varea(x=n_vals, y1=beta_boundary, y2=3, fill_color="gray", fill_alpha=0.7)
# Light gray region = monostable (below boundary)
p.varea(x=n_vals, y1=0, y2=beta_boundary, fill_color="lightgray", fill_alpha=0.5)
# Boundary curve
p.line(n_vals, beta_boundary, color="black", line_width=3)
# Region labels
p.text(x=[5], y=[2.3], text=["bistable"], text_color="black", text_align="center", text_font_size="16pt")
p.text(x=[5], y=[0.9], text=["monostable"], text_color="black", text_align="center", text_font_size="16pt")
return p
show(plot_autoactivator_bifurcation())
3. The Toggle Switch
Dimensionless ODEs: $$\frac{\mathrm{d}x}{\mathrm{d}t} = \frac{\beta_x}{1 + y^n} - x$$ $$\frac{\mathrm{d}y}{\mathrm{d}t} = \frac{\beta_y}{1 + x^n} - y$$
Note: For simplification, we fix $\beta_y$ and only vary $\beta_x$ to observe how changes in x-production affect the system dynamics.
Dynamics and nullcline anaylsis
▼[20] Code
def plot_toggle_switch(beta_x=5.0, n=3.0, x0=0.5, y0=5.0, beta_y=5.0):
# ============================================================
# PART A: CORE MATH -- ODE definition & numerical integration
# (this is what you want to focus on for learning)
# ============================================================
#
# Toggle Switch dimensionless ODEs (asymmetric beta):
# dx/dt = beta_x / (1 + y^n) - x
# dy/dt = beta_y / (1 + x^n) - y
#
# Two genes mutually repressing each other; can produce bistability (1-bit memory)
#
# Parameters:
# beta_x -- production rate of x (adjustable, default: 5.0)
# beta_y -- production rate of y (fixed, default: 5.0)
# n -- Hill coefficient (cooperativity)
#
# User-configurable initial conditions (for dynamics):
# x0 -- starting value for x (default: 0.5)
# y0 -- starting value for y (default: 5.0)
u_arr = np.linspace(0, 8, 200)
t_arr = np.linspace(0, 15, 200)
# --- Fixed-point analysis ---
# Fixed-point condition: x = beta_x / (1 + y^n), y = beta_y / (1 + x^n)
# Jacobian: J = [-1, -df/dy; -dg/dx, -1]
# where df_dy = beta_x * n * y^(n-1) / (1+y^n)^2
# dg_dx = beta_y * n * x^(n-1) / (1+x^n)^2
# Stability: eigenvalues of J have negative real parts
# trace = -2 (always negative), det = 1 - df_dy * dg_dx
# Stable if det > 0
def get_fixed_points(beta_x, beta_y, n):
def equations(vars):
x, y = vars
return [x - beta_x / (1 + y**n), y - beta_y / (1 + x**n)]
# Use 2D grid of initial guesses to find ALL fixed points
guesses = np.linspace(0.1, 7.9, 50)
roots = []
for gx in guesses:
for gy in guesses:
try:
r = fsolve(equations, [gx, gy], full_output=False)
x_star, y_star = r[0], r[1]
if x_star > 0 and y_star > 0:
err = np.abs(equations(r))
if err[0] < 1e-6 and err[1] < 1e-6:
if not any(np.isclose(x_star, ex, atol=1e-3) and np.isclose(y_star, ey, atol=1e-3) for ex, ey in roots):
roots.append((x_star, y_star))
except:
pass
st_pts, ust_pts = [], []
for x_star, y_star in roots:
df_dy = beta_x * n * y_star**(n-1) / (1 + y_star**n)**2
dg_dx = beta_y * n * x_star**(n-1) / (1 + x_star**n)**2
det = 1 - df_dy * dg_dx
if det > 0:
st_pts.append((x_star, y_star))
else:
ust_pts.append((x_star, y_star))
return st_pts, ust_pts
st_pts, ust_pts = get_fixed_points(beta_x, beta_y, n)
st_x, st_y = zip(*st_pts) if st_pts else ([], [])
ust_x, ust_y = zip(*ust_pts) if ust_pts else ([], [])
# --- Numerical ODE integration (CORE) ---
# Using scipy.integrate.odeint (LSODA adaptive solver)
def toggle_rhs(y, t, beta_x, beta_y, n):
"""Toggle switch ODE: dx/dt = beta_x/(1+y^n) - x, dy/dt = beta_y/(1+x^n) - y"""
x, yv = y
dxdt = beta_x / (1 + yv**n) - x
dydt = beta_y / (1 + x**n) - yv
return [dxdt, dydt]
# Integrate from user-specified initial conditions
sol = odeint(toggle_rhs, [x0, y0], t_arr, args=(beta_x, beta_y, n))
u_dyn, v_dyn = sol[:, 0], sol[:, 1]
# ============================================================
# PART B: VISUALIZATION (Bokeh plotting + interactivity, math-independent)
# ============================================================
# Nullclines: x = beta_x / (1 + y^n) and y = beta_y / (1 + x^n)
v_init = beta_x / (1.0 + u_arr**n) # x-nullcline: x as function of y
w_init = beta_y / (1.0 + u_arr**n) # y-nullcline: y as function of x
src_nc = ColumnDataSource(data=dict(u=u_arr, v1=v_init, w1=w_init))
src_st = ColumnDataSource(data=dict(x=st_x, y=st_y))
src_ust = ColumnDataSource(data=dict(x=ust_x, y=ust_y))
src_dyn = ColumnDataSource(data=dict(t=t_arr, u=u_dyn, v=v_dyn))
p_dyn = figure(title="Dynamics over Time", x_axis_label="t", y_axis_label="Concentration", width=400, height=400)
p_dyn.line('t', 'u', source=src_dyn, line_width=3, color='#1f77b4', legend_label="x(t)")
p_dyn.line('t', 'v', source=src_dyn, line_width=3, color='#ff7f0e', legend_label="y(t)")
p_dyn.legend.location = "top_right"
p_dyn.legend.label_text_font_size = "8pt"
p_nc = figure(title="State Space Nullclines", x_axis_label="x", y_axis_label="y", width=400, height=400, x_range=(0, 8), y_range=(0, 8))
p_nc.line('v1', 'u', source=src_nc, line_width=3, color='#1f77b4', legend_label="x-nullcline (dx/dt=0)")
p_nc.line('u', 'w1', source=src_nc, line_width=3, color='#ff7f0e', legend_label="y-nullcline (dy/dt=0)")
p_nc.scatter('x', 'y', source=src_st, size=12, color='black', legend_label="Stable")
p_nc.scatter('x', 'y', source=src_ust, size=12, fill_color='white', line_color='black', line_width=2, legend_label="Unstable")
p_nc.legend.location = "top_right"
p_nc.legend.label_text_font_size = "8pt"
s_bx = Slider(start=0.0, end=10.0, value=beta_x, step=0.5, title="beta_x (adjustable)")
s_n = Slider(start=1, end=5.0, value=n, step=0.1, title="n")
s_x0 = Slider(start=0.0, end=8.0, value=x0, step=0.1, title="x_0")
s_y0 = Slider(start=0.0, end=8.0, value=y0, step=0.1, title="y_0")
# --- CustomJS: browser-side real-time update ---
# Uses Newton's method + Brent's method to find ALL fixed points (stable + unstable).
# This approach is adapted from the biocircuits library (jsplots.py).
# Key insight: Brent's method is bracketed root-finding, so it can find ANY root
# including unstable (saddle) points, unlike simple iteration which only converges
# to stable points.
js_rootfinding = """
function newtonSolve(x0, f, df, args, tol, maxIter, epsilon) {
let x = Infinity;
let solved = false;
tol = tol || 1e-8;
maxIter = maxIter || 200;
epsilon = epsilon || 1e-14;
for (let i = 0; i < maxIter; i++) {
let y = f(x0, ...args);
let yprime = df(x0, ...args);
if (Math.abs(yprime) < epsilon) break;
x = x0 - y / yprime;
if (Math.abs(x - x0) <= tol) { solved = true; break; }
x0 = x;
}
return solved ? x : null;
}
// Bisection method: guaranteed to converge for bracketed roots.
// Used to find the saddle point between two stable fixed points.
function bisectionSolve(f, lower, upper, args, tol, maxIter) {
let a = lower, b = upper;
let fa = f(a, ...args), fb = f(b, ...args);
tol = tol || 1e-8;
maxIter = maxIter || 100;
if (fa * fb > 0) return null;
for (let i = 0; i < maxIter; i++) {
let mid = (a + b) / 2.0;
let fmid = f(mid, ...args);
if (Math.abs(fmid) < tol || (b - a) / 2.0 < tol) return mid;
if (fa * fmid < 0) { b = mid; fb = fmid; }
else { a = mid; fa = fmid; }
}
return (a + b) / 2.0;
}
"""
js_toggle = """
function f(x, beta, n) { return beta / (1.0 + Math.pow(x, n)); }
function ff(x, betax, nx, betay, ny) {
return betax / (1.0 + Math.pow(f(x, betay, ny), nx));
}
function rootFun(x, betax, nx, betay, ny) {
return x - ff(x, betax, nx, betay, ny);
}
function derivff(x, betax, nx, betay, ny) {
let fy = betay / (1.0 + Math.pow(x, ny));
let fynx = Math.pow(fy, nx);
let num = nx * ny * Math.pow(x, ny - 1.0) * betax * fynx * fy;
let denom = Math.pow(betay * (1.0 + fynx), 2);
return num / denom;
}
function derivRootFun(x, betax, nx, betay, ny) {
return 1.0 - derivff(x, betax, nx, betay, ny);
}
function leftRoot(betax, nx, betay, ny) {
return newtonSolve(0.0, rootFun, derivRootFun, [betax, nx, betay, ny]);
}
function rightRoot(betax, nx, betay, ny) {
return newtonSolve(betay, rootFun, derivRootFun, [betax, nx, betay, ny]);
}
function findRoots(betax, nx, betay, ny) {
let x1 = leftRoot(betax, nx, betay, ny);
let x3 = rightRoot(betax, nx, betay, ny);
let args = [betax, nx, betay, ny];
if (x1 === null) {
if (x3 === null) {
return [bisectionSolve(rootFun, 0.0, betay, args)];
}
else return [x3];
}
if (x3 === null) return [x1];
if (Math.abs(x1 - x3) < 2.0 * 1e-4) return [x1];
if (x1 > x3) [x1, x3] = [x3, x1];
let x2 = bisectionSolve(rootFun, x1 + 1e-4, x3 - 1e-4, args);
if (x2 !== null) return [x1, x2, x3];
else return [x1];
}
"""
js_callback = """
const beta_x = s_bx.value, n = s_n.value, x0 = s_x0.value, y0 = s_y0.value;
const beta_y = """ + str(beta_y) + """;
const u = snc.data['u'], v1 = snc.data['v1'], w1 = snc.data['w1'];
// Update nullclines
for(let i=0; i<u.length; i++) {
v1[i] = beta_x / (1 + Math.pow(u[i], n));
w1[i] = beta_y / (1 + Math.pow(u[i], n));
}
snc.change.emit();
// Find ALL fixed points using Newton + Brent methods
let xfp = findRoots(beta_x, n, beta_y, n);
if (xfp === null || xfp[0] === null) {
sst.data = {x: [], y: []};
sust.data = {x: [], y: []};
}
else if (xfp.length === 1) {
sst.data = {x: xfp, y: [f(xfp[0], beta_y, n)]};
sust.data = {x: [], y: []};
}
else {
sst.data = {x: [xfp[0], xfp[2]], y: [f(xfp[0], beta_y, n), f(xfp[2], beta_y, n)]};
sust.data = {x: [xfp[1]], y: [f(xfp[1], beta_y, n)]};
}
sst.change.emit();
sust.change.emit();
// Update dynamics (RK4 integration)
const t = sdyn.data['t'], u_dyn = sdyn.data['u'], v_dyn = sdyn.data['v'];
const dt = t[1] - t[0];
let cu = x0, cv = y0;
for(let i=0; i<t.length; i++) {
u_dyn[i] = cu; v_dyn[i] = cv;
let k1x = beta_x / (1 + Math.pow(cv, n)) - cu;
let k1y = beta_y / (1 + Math.pow(cu, n)) - cv;
let k2x = beta_x / (1 + Math.pow(cv + 0.5*dt*k1y, n)) - (cu + 0.5*dt*k1x);
let k2y = beta_y / (1 + Math.pow(cu + 0.5*dt*k1x, n)) - (cv + 0.5*dt*k1y);
let k3x = beta_x / (1 + Math.pow(cv + 0.5*dt*k2y, n)) - (cu + 0.5*dt*k2x);
let k3y = beta_y / (1 + Math.pow(cu + 0.5*dt*k2x, n)) - (cv + 0.5*dt*k2y);
let k4x = beta_x / (1 + Math.pow(cv + dt*k3y, n)) - (cu + dt*k3x);
let k4y = beta_y / (1 + Math.pow(cu + dt*k3x, n)) - (cv + dt*k3y);
cu += (dt/6) * (k1x + 2*k2x + 2*k3x + k4x);
cv += (dt/6) * (k1y + 2*k2y + 2*k3y + k4y);
}
sdyn.change.emit();
"""
code = js_rootfinding + js_toggle + js_callback
cb = CustomJS(args=dict(snc=src_nc, sdyn=src_dyn, sst=src_st, sust=src_ust,
s_bx=s_bx, s_n=s_n, s_x0=s_x0, s_y0=s_y0), code=code)
s_bx.js_on_change("value", cb)
s_n.js_on_change("value", cb)
s_x0.js_on_change("value", cb)
s_y0.js_on_change("value", cb)
return column(row(s_bx, s_n, s_x0, s_y0), row(p_dyn, p_nc))
# Example usage:
# plot_toggle_switch() # uses defaults: beta_x=5.0, n=3.0
# plot_toggle_switch(beta_x=8.0, n=3.0) # increase production of x
# plot_toggle_switch(beta_x=2.0, n=3.0) # decrease production of x
show(plot_toggle_switch())
Interactive: Dynamics & Nullclines
Stable Unstable
Stable Unstable
Bifurcation Analysis (Toggle Switch)
Symmetric case ($\beta_x = \beta_y = \beta$, $n_x = n_y = n$):
The toggle switch ODEs reduce to a symmetric system. At the symmetric fixed point ($x = y = u$), the Jacobian has $\text{trace} = -2$ (always negative) and $\det = 1 - (df_u)^2$ where $f_u = n u^n / (1+u^n)$. The saddle-node bifurcation occurs when $\det = 0$, i.e., $f_u = 1$. Solving this together with the fixed-point condition $u = \beta/(1+u^n)$ yields the analytical boundary:
$$\beta = \frac{n}{(n-1)^{(n+1)/n}}$$
Above this curve the system is bistable; below it is monostable. Panel below shows this symmetric case.
▼[22] Code
# ============================================================
# BIFURCATION ANALYSIS: Toggle Switch
# Analytical boundary: beta = n / ((n-1)^((n+1)/n))
# ============================================================
def plot_toggle_bifurcation():
# Range of Hill coefficients to sweep
n_vals = np.linspace(1.001, 10, 500)
# Analytical bistability boundary for the toggle switch
beta_boundary = n_vals / ((n_vals - 1)**((n_vals + 1)/n_vals))
# --- Bifurcation diagram ---
p = figure(title="Bifurcation Diagram: Toggle Switch",
x_axis_label="n", y_axis_label="nondimensional beta",
x_range=(1, 10), y_range=(0, 10), width=600, height=400)
# Gray region = bistable (above boundary)
p.varea(x=n_vals, y1=beta_boundary, y2=10, fill_color="gray", fill_alpha=0.7)
# Light gray region = monostable (below boundary)
p.varea(x=n_vals, y1=0, y2=beta_boundary, fill_color="lightgray", fill_alpha=0.5)
# Boundary curve
p.line(n_vals, beta_boundary, color="black", line_width=3)
# Region labels
p.text(x=[5], y=[6], text=["bistable"], text_color="black", text_align="center", text_font_size="16pt")
p.text(x=[3], y=[0.3], text=["monostable"], text_color="black", text_align="center", text_font_size="16pt")
return p
show(plot_toggle_bifurcation())
For you to explore
Try yourself with your AI assistant:
plot_autoactivator_advanced() above, write plot_toggle_switch_advanced() to show the dynamics with parameters changing after certain time points.4. Extended Task: The Repressilator
What happens if we add a third repressor to the loop in a rock-paper-scissors topology? X represses Y, Y represses Z, Z represses X.
This breaks the bistability and instead generates sustained oscillations (a stable limit cycle) under the right parameters.
Dimensionless ODEs (symmetric case): $$\frac{\mathrm{d}x}{\mathrm{d}t} = \frac{\beta}{1 + z^n} - x$$ $$\frac{\mathrm{d}y}{\mathrm{d}t} = \frac{\beta}{1 + x^n} - y$$ $$\frac{\mathrm{d}z}{\mathrm{d}t} = \frac{\beta}{1 + y^n} - z$$
Dynamics
Because the repressilator is a 3D system, we project the dynamics onto the x-y plane to visualize the phase portrait. The limit cycle appears as a closed loop in this 2D projection. The symmetric fixed point (where x=y=z) is marked:
▼[26] Code
def plot_repressilator(x0=1.0, y0=2.0, z0=3.0):
# ============================================================
# PART A: CORE MATH -- ODE definition & numerical integration
# (this is what you want to focus on for learning)
# ============================================================
#
# Repressilator dimensionless ODEs:
# dx/dt = beta / (1 + z^n) - x
# dy/dt = beta / (1 + x^n) - y
# dz/dt = beta / (1 + y^n) - z
#
# Three genes in a ring repression topology; can produce sustained oscillations (limit cycle)
#
# Parameters:
# beta -- maximal expression rate
# n -- Hill coefficient (cooperativity)
#
# User-configurable initial conditions (for dynamics):
# x0 -- starting value for protein X (default: 1.0)
# y0 -- starting value for protein Y (default: 2.0)
# z0 -- starting value for protein Z (default: 3.0)
t_arr = np.linspace(0, 50, 1000)
beta_init, n_init = 50.0, 3.0
# --- Fixed-point & stability analysis ---
# Symmetric fixed point: x = y = z = u, satisfying u = beta/(1+u^n)
# Hopf bifurcation condition: c = n*u^n/(1+u^n) > 2 becomes unstable -> oscillations
u_star = fsolve(lambda u: beta_init / (1.0 + u**n_init) - u, 1.0)[0]
c_val = n_init * (u_star**n_init) / (1.0 + u_star**n_init)
if c_val > 2.0:
st_x, st_y, ust_x, ust_y = [], [], [u_star], [u_star]
else:
st_x, st_y, ust_x, ust_y = [u_star], [u_star], [], []
# --- Numerical ODE integration (CORE) ---
# Using scipy.integrate.odeint (LSODA adaptive solver)
# Signature: odeint(func, y0, t, args=...)
# func(y, t, *args) returns dy/dt (list)
def repressilator_rhs(y, t, beta, n):
"""Repressilator ODE: 3-node ring oscillator"""
x, yv, z = y
return [beta / (1 + z**n) - x, beta / (1 + x**n) - yv, beta / (1 + yv**n) - z]
# Integrate from user-specified initial conditions
sol = odeint(repressilator_rhs, [x0, y0, z0], t_arr, args=(beta_init, n_init))
x_dyn, y_dyn, z_dyn = sol[:, 0], sol[:, 1], sol[:, 2]
# ============================================================
# PART B: VISUALIZATION (Bokeh plotting + interactivity, math-independent)
# ============================================================
src_st = ColumnDataSource(data=dict(x=st_x, y=st_y))
src_unst = ColumnDataSource(data=dict(x=ust_x, y=ust_y))
src_dyn = ColumnDataSource(data=dict(t=t_arr, x=x_dyn, y=y_dyn, z=z_dyn))
# 1. Time dynamics plot
p_dyn = figure(title="Repressilator Time Dynamics", x_axis_label="Time (t)", y_axis_label="Concentration", width=550, height=350)
p_dyn.line('t', 'x', source=src_dyn, line_width=2, color='#1f77b4', legend_label="Protein X")
p_dyn.line('t', 'y', source=src_dyn, line_width=2, color='#ff7f0e', legend_label="Protein Y")
p_dyn.line('t', 'z', source=src_dyn, line_width=2, color='#2ca02c', legend_label="Protein Z")
p_dyn.legend.location = "top_right"
p_dyn.legend.label_text_font_size = "8pt"
# 2. 2D Phase Portrait (x vs y projection)
# Project the 3D limit cycle onto the x-y plane to visualize the nullcline structure
p_phase = figure(title="2D Phase Portrait (x vs y)", x_axis_label="Protein X", y_axis_label="Protein Y", width=400, height=400)
p_phase.line('x', 'y', source=src_dyn, line_width=2, color='purple', alpha=0.8)
p_phase.scatter('x', 'y', source=src_st, size=12, color='black', legend_label="Stable FP")
p_phase.scatter('x', 'y', source=src_unst, size=12, fill_color='white', line_color='black', line_width=2, legend_label="Unstable FP")
p_phase.legend.location = "top_right"
p_phase.legend.label_text_font_size = "8pt"
s_beta = Slider(start=1.0, end=100.0, value=beta_init, step=1.0, title="beta")
s_n = Slider(start=0.1, end=5.0, value=n_init, step=0.1, title="n")
# --- CustomJS: browser-side real-time update (no Python kernel needed) ---
# Uses Newton's method + bisection to find the symmetric fixed point.
code = """
const beta = s_beta.value, n = s_n.value;
// --- Root-finding helpers (Newton's method + bisection) ---
// Symmetric fixed point condition: u = beta/(1+u^n)
// rootFun(u) = u - beta/(1+u^n) = 0
const rootFun = (u) => u - beta / (1 + Math.pow(u, n));
const drootFun = (u) => {
if (u === 0) return 1.0;
const un = Math.pow(u, n);
return 1.0 + beta * n * Math.pow(u, n - 1) / Math.pow(1 + un, 2);
};
function newtonSolve(x0, tol, maxIter) {
tol = tol || 1e-8;
maxIter = maxIter || 100;
let x = x0;
for (let i = 0; i < maxIter; i++) {
let y = rootFun(x);
let yp = drootFun(x);
if (Math.abs(yp) < 1e-14) break;
let xNew = x - y / yp;
if (xNew < 0) xNew = x / 2;
if (Math.abs(xNew - x) <= tol) return xNew;
x = xNew;
}
return x;
}
function bisectionSolve(lower, upper, tol, maxIter) {
tol = tol || 1e-8;
maxIter = maxIter || 100;
let a = lower, b = upper;
let fa = rootFun(a), fb = rootFun(b);
if (fa * fb > 0) return null;
for (let i = 0; i < maxIter; i++) {
let mid = (a + b) / 2.0;
let fmid = rootFun(mid);
if (Math.abs(fmid) < tol || (b - a) / 2.0 < tol) return mid;
if (fa * fmid < 0) { b = mid; fb = fmid; }
else { a = mid; fa = fmid; }
}
return (a + b) / 2.0;
}
let u_star = newtonSolve(beta);
if (u_star === null || Math.abs(rootFun(u_star)) > 1e-4) {
u_star = bisectionSolve(0.0, beta);
}
let c = 0;
if (u_star !== null && u_star > 0) {
const un = Math.pow(u_star, n);
c = n * un / (1 + un);
}
if (u_star !== null && c > 2.0) {
sst.data = {x: [], y: []};
sust.data = {x: [u_star], y: [u_star]};
} else if (u_star !== null) {
sst.data = {x: [u_star], y: [u_star]};
sust.data = {x: [], y: []};
} else {
sst.data = {x: [], y: []};
sust.data = {x: [], y: []};
}
sst.change.emit(); sust.change.emit();
// Browser-side ODE integration via RK4
const t = sdyn.data['t'];
const x = sdyn.data['x'], y = sdyn.data['y'], z = sdyn.data['z'];
const dt = t[1] - t[0];
let cx = """ + str(x0) + """, cy = """ + str(y0) + """, cz = """ + str(z0) + """;
const f = (v) => beta / (1 + Math.pow(v, n));
for(let i=0; i<t.length; i++) {
x[i] = cx; y[i] = cy; z[i] = cz;
let k1x = f(cz) - cx, k1y = f(cx) - cy, k1z = f(cy) - cz;
let k2x = f(cz + 0.5*dt*k1z) - (cx + 0.5*dt*k1x);
let k2y = f(cx + 0.5*dt*k1x) - (cy + 0.5*dt*k1y);
let k2z = f(cy + 0.5*dt*k1y) - (cz + 0.5*dt*k1z);
let k3x = f(cz + 0.5*dt*k2z) - (cx + 0.5*dt*k2x);
let k3y = f(cx + 0.5*dt*k2x) - (cy + 0.5*dt*k2y);
let k3z = f(cy + 0.5*dt*k2y) - (cz + 0.5*dt*k2z);
let k4x = f(cz + dt*k3z) - (cx + dt*k3x);
let k4y = f(cx + dt*k3x) - (cy + dt*k3y);
let k4z = f(cy + dt*k3y) - (cz + dt*k3z);
cx += (dt/6) * (k1x + 2*k2x + 2*k3x + k4x);
cy += (dt/6) * (k1y + 2*k2y + 2*k3y + k4y);
cz += (dt/6) * (k1z + 2*k2z + 2*k3z + k4z);
}
sdyn.change.emit();
"""
cb = CustomJS(args=dict(sdyn=src_dyn, sst=src_st, sust=src_unst, s_beta=s_beta, s_n=s_n), code=code)
s_beta.js_on_change('value', cb)
s_n.js_on_change('value', cb)
return column(row(s_beta, s_n), row(p_dyn, p_phase))
# Example usage:
# plot_repressilator() # uses defaults: x0=1.0, y0=2.0, z0=3.0
# plot_repressilator(x0=5.0, y0=0.1, z0=2.0) # custom initial conditions
show(plot_repressilator())
Interactive: Dynamics & Phase Portrait
Stable FP Unstable FP
▼[27] Code
# ============================================================
# 3D VISUALIZATION: Repressilator limit cycle
# ODEs: dx/dt = beta / (1 + z^n) - x
# dy/dt = beta / (1 + x^n) - y
# dz/dt = beta / (1 + y^n) - z
# Uses scipy.integrate.odeint + Plotly for 3D interactive plot
# ============================================================
def plot_repressilator_3d(beta=50.0, n=3.0):
# Time vector for integration (long enough to see limit cycle)
t_arr = np.linspace(0, 50, 2000)
# Repressilator ODE right-hand side (3-node ring oscillator)
def repressilator_rhs(y, t, beta, n):
x, yv, z = y
return [beta / (1 + z**n) - x, beta / (1 + x**n) - yv, beta / (1 + yv**n) - z]
# --- Numerical ODE integration ---
sol = odeint(repressilator_rhs, [1.0, 2.0, 3.0], t_arr, args=(beta, n))
# --- 3D interactive plot (Plotly) ---
# Color along the trajectory shows time progression (Viridis colormap)
fig = go.Figure(data=[go.Scatter3d(
x=sol[:, 0], y=sol[:, 1], z=sol[:, 2],
mode='lines',
line=dict(color=t_arr, colorscale='Viridis', width=4)
)])
fig.update_layout(
scene=dict(
xaxis_title='Protein X_1',
yaxis_title='Protein X_2',
zaxis_title='Protein X_3'
),
width=700, height=600
)
fig.show()
plot_repressilator_3d()
Plotly 3D interactive plot - open in Jupyter Lab to view
Bifurcation Analysis (Repressilator)
▼[29] Code
# ============================================================
# BIFURCATION ANALYSIS: Repressilator
# Hopf bifurcation boundary: beta = (n/(n-2)) * ((2/(n-2))^(1/n))
# Above the curve: sustained oscillations; below: stable steady state
# Asymptote at n = 2 (oscillations require n > 2)
# ============================================================
def plot_repressilator_bifurcation():
# n must be > 2 for oscillations (Hopf condition)
n_vals = np.linspace(2.001, 10, 1000)
# Analytical Hopf bifurcation boundary
beta_boundary = (n_vals / (n_vals - 2)) * ((2 / (n_vals - 2))**(1/n_vals))
# Clip the upper boundary strictly for visualization filling
beta_boundary_clipped = np.clip(beta_boundary, 0, 200)
# --- Bifurcation diagram ---
p = figure(title="Bifurcation Diagram: Repressilator",
x_axis_label="n", y_axis_label="nondimensional beta",
x_range=(0, 10), y_range=(0, 100), width=600, height=400)
# Gray region = oscillatory regime (unstable steady state, above boundary)
p.varea(x=n_vals, y1=beta_boundary_clipped, y2=200, fill_color="gray", fill_alpha=0.7)
# Light gray region = stable steady state (no oscillations, below boundary)
# For n <= 2, the entire beta range is stable
n_mono = np.linspace(0, 10, 1000)
y2_mono = np.zeros_like(n_mono)
for i, nv in enumerate(n_mono):
if nv <= 2:
y2_mono[i] = 200
else:
y2_mono[i] = min(200, (nv / (nv - 2)) * ((2 / (nv - 2))**(1/nv)))
p.varea(x=n_mono, y1=0, y2=y2_mono, fill_color="lightgray", fill_alpha=0.5)
# Boundary curve (Hopf bifurcation)
p.line(n_vals, beta_boundary, color="black", line_width=3)
# Vertical asymptote at n = 2
p.line([2.0, 2.0], [0, 200], color="black", line_width=3, line_dash="dashed")
# Region labels
p.text(x=[6], y=[60], text=["sustained oscillations"], text_color="black", text_align="center", text_font_size="16pt")
p.text(x=[1.3], y=[5], text=["stable \n steady state"], text_color="black", text_align="center", text_font_size="16pt")
return p
show(plot_repressilator_bifurcation())
