Python · chapter 5 of 13

Simulating expression data

sol = grn.simulate(model, ...) integrates a model and returns a Solution. model may be a Model, a motif name, a Motif, or a signed adjacency (the last three are built with grn.build_model using the logic and seed you pass).

Time grid

Argument	Default	Meaning
`duration`	50	Recorded time span
`sample_dt`	0.1	Recording step; `sol.t = np.arange(0, duration + sample_dt, sample_dt)`
`transient`	0	Time integrated and discarded before recording

Recording starts at t = 0. The transient, if any, runs at negative times, which matters for stimuli (their clock is sol.t).

Deterministic integration

solver="ode45" (default, SciPy's RK45) or "ode15s" (SciPy's BDF). Both use rtol=1e-7, atol=1e-9; the solution is clipped at zero. BDF receives the analytic Jacobian and is the better choice for stiff systems: widely different degradation rates, very high Hill coefficients, or two-stage models with slow proteins. The integration restarts at every stimulus edge, so step inputs are handled exactly rather than smeared across a solver step.

Cells and initial conditions

Argument	Default	Meaning
`cells`	1	Number of independent trajectories. `sol.x` becomes `(samples, dim, cells)`
`random_ic`	False	Draw each cell's initial state uniformly in `[0, scale]` per state
`initial_condition`	`model.ic`	Length-`dim` array shared by all cells, or `(dim, cells)`

Cell k uses seed seed + k for its initial condition and its process noise. With cells > 1 and neither random_ic nor noise, every cell is identical; combine cells with random initial conditions, process noise, or both.

Perturbed genes (see Perturbations) have their random initial conditions scaled by the perturbation factor, so a knocked-out gene starts at zero.

sol = grn.simulate("repressilator", cells=20, random_ic=True, duration=30, seed=3)
grn.plot_expression(sol, cell=0)

Twenty cells of the repressilator from random initial conditions. Cell 0 is drawn in full; the other 19 sit faintly behind it. Every cell oscillates with the same period but its own phase.

The same 20 cells: single traces on the left, the population mean on the right. Dephasing flattens the mean within two periods, which is exactly why bulk measurements of oscillating circuits look constant.

Process noise

noise switches from the ODE solvers to a fixed-step Euler–Maruyama integrator with step dt (default sample_dt / 20; sample_dt must be an integer multiple of dt). States are clipped at zero after every step.

`noise`	Diffusion term per state	Use when
`"none"` (default)		Deterministic ODE
`"additive"`	`sigma`	Simple, level-independent fluctuations; the linear case has exact statistics
`"multiplicative"`	`sigma * x`	Fluctuations proportional to expression, never negative
`"langevin"`	`sigma * sqrt(birth + decay * x)`	Chemical Langevin equation: intrinsic noise from production and degradation events. `sigma` plays the role of `1 / sqrt(system size)`; the Fano factor of an unregulated gene equals `sigma**2`

noise_strength is sigma (default 0.1). For two-stage models the Langevin birth rate of an mRNA is its transcription rate and that of a protein is its translation rate.

The three process-noise models on the repressilator, five cells each. Additive noise jitters uniformly, multiplicative noise scales with expression, and the chemical Langevin form is quiet near zero and loudest at the peaks.

Observation noise

noise_sigma adds Gaussian measurement noise after integration and clips at zero. It accepts a scalar, one value per gene (shared by that gene's mRNA and protein), or one value per state. The clean data are kept in sol.x_clean.

Genetic perturbations and stimuli

knockout, knockdown, overexpress, knockdown_factor, overexpress_factor are forwarded to grn.perturb; stimulus is forwarded to grn.stimulate. See the next two chapters.

Reading the result

sol.t                  # (samples,)
sol.x                  # (samples, dim) or (samples, dim, cells)
sol.state_labels       # names of the dim columns
sol.mrna, sol.protein  # two-stage models only
sol.u                  # (samples, n) multiplicative input (all ones without stimuli)
sol.ic                 # (dim, cells) initial conditions actually used
sol.options            # every option, with defaults filled in
sol.model              # the generating model

To stack all cells into one (samples * cells, dim) array (time within cell):

X = sol.x.transpose(0, 2, 1).reshape(-1, sol.dim)

A solution printed at the prompt:

sol = grn.simulate("repressilator", cells=4, random_ic=True, noise_sigma=0.1)
sol

Output

Solution(name='repressilator', x.shape=(501, 3, 4), cells=4, stages=1, noise='none')
t: (501,)  x: (501, 3, 4)  x_clean: (501, 3, 4)  ic: (3, 4)  u: (501, 3)
state_labels: ('lacI', 'tetR', 'cI')
options: {'duration': 50.0, 'sample_dt': 0.1, 'cells': 4, 'noise': 'none', 'seed': 0}

Choosing settings

Relaxation data (responses from random initial conditions): random_ic=True, duration around 5 to 10 mRNA lifetimes, no transient.
Stationary data (fluctuations around a steady state or oscillations): a transient of at least 20 lifetimes, then a long duration with process noise.
Single-cell-like ensembles: many cells, random_ic, noise="langevin", then keep one or a few time points per cell.
Bulk-like data: average sol.x over the cell axis; oscillators dephase and average out, which is the correct behaviour.

Performance notes

Cost is linear in cells and in duration / dt for stochastic runs. The Euler–Maruyama loop is pure Python over time steps, so long noisy runs of large networks are the slowest case; a 30-gene network with 50 cells and the default dt takes on the order of a minute. The multiplicative-logic Jacobian is O(n^3), which is only noticeable with solver="ode15s" above a few hundred genes.

← Building models Perturbations, steady states, and knockout screens →