GRNLab (Python)

A data generator for gene regulatory networks (GRNs) with Hill-function kinetics. Every dataset comes with the ground truth that produced it: the signed regulatory network, every kinetic parameter, the initial conditions, the seeds, and the solver settings. GRNLab generates data; it does not infer, recover, or score networks.

This is the Python implementation of the GRNLab MATLAB toolbox. The two share the same functions, options, conventions, kinetics, integrators, accuracy audit, and figure style. Seeds are reproducible within each implementation but not across the two (different random-number streams).

Installation

pip install grnlab

Requires Python 3.9 or newer, NumPy, SciPy, and Matplotlib.

Ten lines to a dataset

import grnlab as grn

grn.catalog()                                          # motifs and random topologies

S     = grn.random_network(30, topology="scalefree", seed=2)          # ground truth
model = grn.build_model(S, logic="multiplicative", seed=2)            # kinetics
sol   = grn.simulate(model, cells=50, random_ic=True,
                     noise="langevin", noise_strength=0.2, seed=2)    # data

sol.x.shape            # (501, 30, 50): samples x genes x cells
sol.S                  # the signed adjacency that generated the data
grn.plot_network(model)            # arrowheads = activation, bars = repression
grn.plot_expression(sol, cell=2)
grn.export(sol, "grn30.csv")       # + grn30_truth.csv, grn30_adjacency.csv

Conventions

• S[i, j] = +1 means gene j activates gene i, -1 that it represses it, 0 no direct regulation. Row i lists the regulators of gene i.
• Gene indices are 0-based. Labels ("G1", "G2", ... by default) may be used instead of indices wherever a gene is named: perturbations, stimuli, screens.
• Kinetics: dx_i/dt = production_i(x) - gamma_i x_i with Hill functions h+(x_j) = x_j^h / (K^h + x_j^h) for activators and 1 - h+ for repressors, combined additively (basal + sum k_ij R_ij, OR-like) or multiplicatively (basal + vmax_i prod_j R_ij, AND-like). Genes without regulators are constitutive at vmax_i.
• stages=2 gives every gene an mRNA and a protein state; regulators act through their proteins. The state is [mRNA, protein]; solutions expose sol.mrna and sol.protein.
• Every random draw is governed by a seed (default 0). Cell k of a multi-cell simulation uses seed + k. Pass seed=None to leave the generator unseeded.

Motifs by name

sol = grn.simulate("repressilator")                  # oscillates (preset applied)
model = grn.build_model("toggle")                    # bistable
hi = grn.steady_state(model, initial_condition=[5, 0])
lo = grn.steady_state(model, initial_condition=[0, 5])
grn.plot_motifs()                                    # gallery of all ten

Registered motifs: toggle, repressilator, nar, par, ffl_coherent, ffl_incoherent, cascade, sim, bifan, goodwin.

Perturbations, steady states, screens

model = grn.build_model(grn.random_network(10, seed=3), seed=3)
sol = grn.simulate(model, knockout=3, overexpress="G7")
ko  = grn.knockout_screen(model, noise_sigma=0.05)   # ko.X: (n+1) x n, row 0 = wild type
kd  = grn.knockout_screen(model, type="knockdown")
xss, info = grn.steady_state(model, info=True)       # info["converged"], info["stable"]

Time-varying stimuli

on   = grn.stimulus("X", "step", baseline=0, amplitude=1, start=5)   # switch X on at t = 5
sol  = grn.simulate("ffl_incoherent", stimulus=on, duration=30)
sol.u                                                                # realised input, samples x n

pulse = grn.stimulus(0, "step", amplitude=3, start=5, stop=10)
ramp  = grn.stimulus(0, "ramp", baseline=0, amplitude=3, start=0, stop=40)
drive = grn.stimulus("G2", "sine", amplitude=3, period=12)
xss   = grn.steady_state(model, stimulus=ramp, time=20)              # stimulus held at t = 20

Stimuli scale (or, with mode="add", add to) a gene's production rate. The ODE integration restarts at every stimulus edge, so steps are exact.

Two-stage mRNA/protein kinetics

model = grn.build_model("repressilator", stages=2)
sol = grn.simulate(model, duration=60)
sol.x.shape, sol.state_labels        # (601, 6), ('lacI_mRNA', ..., 'cI_protein')
grn.plot_expression(sol, species="mrna")

Figures

fig = grn.plot_network(model, layout="layered")     # circle, force, or layered
grn.export_figure(fig, "network.pdf")               # vector; .png at 300 dpi by default
fig = grn.plot_expression(sol)
grn.export_figure(fig, "expression.png", resolution=600)

Export

grn.export(data, "name.csv") writes the data table plus name_truth.csv (edge list: source, target, sign) and name_adjacency.csv. .npz stores the arrays and a pickled copy of the whole object under data; .pkl stores the object.

Documentation

The full user guide (13 chapters with worked examples, their output, and figures) is at almomaa.github.io/GRNLab/python.

Function reference

Every function has a full docstring (help(grn.simulate)).

Accuracy

python -m grnlab.audit

runs 19 checks of the public API against closed-form results (exponential relaxation, Hill function values, analytic steady states, the repressilator's Hopf threshold, the toggle switch's bistability threshold, Euler–Maruyama convergence, Ornstein–Uhlenbeck and chemical-Langevin statistics, two-stage limits, exact step stimuli, knockout sign agreement, seeding) and writes AUDIT.md.

Development

pip install -e ".[test]"
pytest
python -m build            # sdist and wheel in dist/

License

MIT. Author: Abd AlRahman AlMomani, Embry-Riddle Aeronautical University.