Python · chapter 6 of 13

Perturbations, steady states, and knockout screens

Genetic perturbations

grn.perturb scales the production of chosen genes: the basal rate, the maximal rate vmax, every incoming weight k[i, :], and the gene's entries in model.ic. Genes may be given as 0-based indices or labels. It returns a modified copy; the original model is untouched.

Argument	Factor	Effect
`knockout`	0	No production; the gene decays to zero and stays there
`knockdown`	`knockdown_factor` (0.5)	Reduced production
`overexpress`	`overexpress_factor` (2)	Increased production

ko = grn.perturb(model, knockout=[1, 4])
kd = grn.perturb(model, knockdown="G7", knockdown_factor=0.25)
sol = grn.simulate(model, knockout=3, overexpress="G7")   # same arguments

The perturbation is recorded in model.perturbation["gene"] and ["factor"]. Perturbations accumulate if you call grn.perturb on an already perturbed model; grn.simulate, grn.steady_state, and grn.knockout_screen always start from the model you pass them.

Knockouts act on production, not on the network: the ground truth S is unchanged, which is what a benchmark needs. In two-stage models the transcription of the gene is scaled, so both its mRNA and its protein go to zero.

model = grn.build_model("cascade", seed=1)
wt = grn.simulate(model, initial_condition=np.zeros(4), duration=20)
ko = grn.simulate(model, initial_condition=np.zeros(4), duration=20, knockout="G2")

A four-gene activation cascade switched on from zero, wild type and with G2 knocked out. G2 stays at zero and everything downstream is silenced, while the constitutive G1 is unaffected.

Steady states

xss = grn.steady_state(model, ...) finds a fixed point in two phases: it integrates with BDF (using the analytic Jacobian) in windows of 20 / min(decay) until the relative residual norm(f) / (1 + norm(x)) drops below tol (1e-8) or max_time (2000) is exceeded, then polishes with Newton steps until the residual is at machine precision. xss, info = grn.steady_state(model, ..., info=True) also returns:

Key	Meaning
`converged`	Whether a fixed point was reached
`stable`	All Jacobian eigenvalues at `xss` have negative real part
`residual`, `time`	Final `norm(f)` and integration time used

If the model never settles (an oscillator), converged is False and xss is the time average over the last window.

Arguments: initial_condition (the starting point, which selects the attractor of a multistable model), max_time, tol, the perturbation arguments, and stimulus with time (stimuli are frozen at their value at time, default 0).

model = grn.build_model("toggle")
hi, info = grn.steady_state(model, initial_condition=[5, 0], info=True)   # u high
lo = grn.steady_state(model, initial_condition=[0, 5])                    # v high
_, info = grn.steady_state("repressilator", info=True)   # info["converged"] is False

Output

hi = [4.7913 0.2087]
info = {'converged': True, 'residual': 2.220446049250313e-16, 'time': 40.0, 'stable': True, 'time0': 0.0}

Knockout screens

screen = grn.knockout_screen(model, ...) produces the steady-state dataset used by the DREAM challenges: the wild type followed by one row per single-gene perturbation.

Argument	Default	Meaning
`type`	`"knockout"`	`"knockout"`, `"knockdown"`, or `"overexpress"`
`factor`	0, 0.5, or 2 by type	Production scaling
`genes`	all	Which genes to perturb, indices or labels
`noise_sigma`	0	Gaussian measurement noise on `screen.X` (clean copy in `screen.x_clean`)
`seed`	0	Seed for the measurement noise
`stimulus`, `time`		Inputs held at their value at `time`
`max_time`, `tol`		As in `grn.steady_state`

Each perturbed steady state is computed from the wild-type steady state with the perturbed gene's states scaled by the factor, which keeps the screen on the same attractor branch as the wild type in multistable networks.

Attribute	Meaning
`X`	`(len(genes) + 1, dim)`; row 0 wild type, row i + 1 the perturbation of `genes[i]`
`rows`	Row names: `"wildtype"` then the gene labels
`converged`	Boolean per row
`mrna`, `protein`	The two halves of `X` for two-stage models
`S`, `model`, `labels`, `state_labels`, `options`	Ground truth and settings

ko = grn.knockout_screen(model, noise_sigma=0.05)
kd = grn.knockout_screen(model, type="knockdown", genes=["G1", "G4"])
wt = ko.X[0]
lfc = np.log2((ko.X[1:] + 1e-3) / (wt + 1e-3))     # log fold change per knockout
grn.export(ko, "screen.csv")

A 12-gene network with two autoregulating genes beside its single-knockout screen, shown as log2 fold change against wild type with the knocked-out gene itself masked. Knocking out a repressor raises its targets (blue); knocking out an activator lowers them (red). Every resolvable off-diagonal effect has the sign the network predicts.

The first rows and columns of the screen:

ko.rows[:5], ko.labels[:6], np.round(ko.X[:5, :6], 3)

Output

rows: ('wildtype', 'G1', 'G2', 'G3', 'G4')
columns: ('G1', 'G2', 'G3', 'G4', 'G5', 'G6')
[[1.108 1.216 2.424 3.384 3.329 0.281]
 [0.    1.213 2.423 3.407 3.316 0.161]
 [0.988 0.    2.458 2.766 3.671 0.256]
 [1.25  1.152 0.    4.036 2.943 0.315]
 [0.797 4.739 2.504 0.    4.301 0.221]]

A useful sanity check on any screen: knocking out an activator lowers its targets and knocking out a repressor raises them. The audit in Accuracy and validation verifies this on a random network.

Double perturbations and dose series

Combine the arguments directly:

sol = grn.simulate(model, knockout=0, knockdown=2)
X = np.array([grn.steady_state(model, knockdown=1, knockdown_factor=f)
              for f in (0, 0.25, 0.5, 0.75, 1)])

← Simulating expression data Time-varying stimuli →