Python · chapter 10 of 13

Recipes

Complete, copy-and-run workflows. Each one ends with data and its ground truth on disk.

1. A benchmark time-series dataset

Fifty cells of a 30-gene scale-free network with intrinsic and measurement noise.

import grnlab as grn

S     = grn.random_network(30, topology="scalefree", mean_degree=2,
                           fraction_repressors=0.4, seed=2)
model = grn.build_model(S, logic="multiplicative", seed=2, name="sf30")
sol   = grn.simulate(model, cells=50, random_ic=True, duration=20, sample_dt=0.25,
                     noise="langevin", noise_strength=0.2, noise_sigma=0.05, seed=2)
grn.export(sol, "sf30_timeseries.csv")
grn.export(sol, "sf30_timeseries.pkl")
grn.export_figure(grn.plot_network(model), "sf30_network.pdf")
grn.export_figure(grn.plot_expression(sol, cell=0), "sf30_expression.png")

Recipe 1 output: one cell of the benchmark dataset with the other cells behind it.

2. A DREAM-style package: knockouts, knockdowns, and wild-type time series

model = grn.build_model(grn.random_network(10, seed=3), seed=3, name="net10")
ko = grn.knockout_screen(model, noise_sigma=0.02, seed=3)
kd = grn.knockout_screen(model, type="knockdown", noise_sigma=0.02, seed=4)
ts = grn.simulate(model, cells=10, random_ic=True, duration=20, noise_sigma=0.02, seed=5)
grn.export(ko, "net10_knockouts.csv")
grn.export(kd, "net10_knockdowns.csv")
grn.export(ts, "net10_timeseries.csv")
assert ko.converged.all() and kd.converged.all()     # check before shipping

3. A dose-response experiment

Steady-state expression of every gene as a source gene is driven from 0 to 3 times its normal rate.

import numpy as np
model = grn.build_model(grn.motif("sim", size=5), seed=1)   # keeps labels TF, T1..T4
doses = np.linspace(0, 3, 13)
X = np.array([grn.steady_state(model, stimulus=grn.stimulus("TF", "step", baseline=d, amplitude=d))
              for d in doses])
import matplotlib.pyplot as plt
plt.plot(doses, X[:, 1:]); plt.xlabel("input to TF"); plt.ylabel("target expression")

A ramp gives the same information dynamically:

ramp = grn.stimulus("TF", "ramp", baseline=0, amplitude=3, start=0, stop=60)
sol = grn.simulate(model, stimulus=ramp, duration=80)
grn.plot_expression(sol)

4. Single-cell-like snapshots

Many cells, one sample each at a random time after a stimulus.

model = grn.build_model(grn.random_network(20, seed=6), stages=2, seed=6)
on = grn.stimulus(0, "step", amplitude=4, start=0)
sol = grn.simulate(model, stimulus=on, cells=200, random_ic=True,
                   noise="langevin", noise_strength=0.3, duration=10, sample_dt=0.1, seed=6)
rng = np.random.default_rng(6)
t_idx = rng.integers(sol.t.size, size=sol.cells)              # one time per cell
cells_x = np.array([sol.mrna[t_idx[c], :, c] for c in range(sol.cells)])   # mRNA, like scRNA-seq
pseudotime = sol.t[t_idx]

Recipe 4 output. Left: two genes across 200 cells, coloured by sampling time. Right: the snapshot matrix with cells ordered by their (hidden) sampling time.

5. Comparing motif variants

Coherent and incoherent feed-forward loops under the same step input.

on = grn.stimulus("X", "step", baseline=0, amplitude=1, start=2)
for name in ("ffl_coherent", "ffl_incoherent"):
    model = grn.build_model(name, logic="multiplicative", seed=1,
                            params=dict(gamma=[5, 1, 3], vmax=5, basal=[5, 0, 0]))
    sol = grn.simulate(model, stimulus=on, initial_condition=np.zeros(3), duration=20)
    grn.export_figure(grn.plot_expression(sol), f"{name}.pdf")

Recipe 5 output. The coherent loop delays Z and then holds it; the incoherent loop produces a pulse.

Ten networks that share a topology but differ in kinetics.

S = grn.random_network(15, seed=1)
for k in range(10):
    model = grn.build_model(S, seed=k, name=f"net15_{k}")
    sol = grn.simulate(model, cells=5, random_ic=True, seed=k)
    grn.export(sol, f"net15_{k:02d}.pkl")

7. Sweeping a parameter of a motif

Where does the repressilator start to oscillate? (The analytic answer for Hill 4 is alpha = 2.)

alphas = np.arange(1.5, 3.01, 0.1)
amp = []
for a in alphas:
    m = grn.build_model("repressilator", params=dict(k=a, vmax=a))
    sol = grn.simulate(m, duration=100, transient=300)
    amp.append(np.ptp(sol.x[:, 0]))
plt.plot(alphas, amp, "o-"); plt.xlabel("alpha"); plt.ylabel("amplitude")

Recipe 7 output (left) and the toggle switch's bifurcation diagram (right). Both thresholds land on the analytic value alpha = 2.

← Exporting data and reproducing it Accuracy and validation →