Python · chapter 4 of 13

Building models

grn.build_model turns a signed adjacency into a simulatable model by drawing kinetic parameters, then applying any preset or override you give it.

model = grn.build_model(S)                                  # defaults
model = grn.build_model(S, logic="multiplicative", seed=3)
model = grn.build_model("repressilator", stages=2)          # motif + preset
model = grn.build_model(S, params=dict(hill=2, gamma=1))

How parameters are sampled

With numpy.random.default_rng(seed), the parameters are drawn in this order (which is why two models built from the same S and seed but different stages share their transcription parameters):

Parameter	Drawn as	Default range argument
`gamma_i`	uniform	`gamma_range=(0.5, 1.5)`
`vmax_i`	uniform	`vmax_range=(1, 5)`
`gammaP_i`	uniform (used only for two stages)	`gammap_range=(0.2, 1)`
`beta_i`	`gammaP_i` times a uniform ratio	`beta_ratio_range=(1, 1)`
`hill_ij`	uniform, on edges only	`hill_range=(1, 4)`
`K_ij`	`xchar_j` times a uniform factor	`k_rel_range=(0.3, 1.2)`
`k_ij`	`vmax_i / indegree_i` times uniform in (0.5, 1.5), on edges
`basal_i`	`vmax_i` for source genes; `vmax_i` times uniform fraction otherwise	`basal_fraction=(0, 0.1)`
`ic`	uniform in `[0, scale]` per state

xchar_j = vmax_j / gamma_j is the characteristic level of regulator j (multiplied by beta_j / gammaP_j for two-stage models, since regulators act through protein). Drawing K_ij relative to it guarantees that every edge is actually "in play": the regulator crosses its threshold within its natural range. Drawing k_ij as a share of vmax_i keeps the additive production of a gene near vmax_i when all its activators are on.

Off-edge entries are sanitised after every override: k = 0, K = 1, hill = 1 wherever S = 0, so they never affect the dynamics.

The sampled parameters of a 400-gene model with the default ranges. The last panel shows the threshold of every edge divided by its regulator's characteristic level: every threshold lies inside the regulator's range.

Overriding parameters

params is a dict whose entries override the sampled values. Motif presets are applied first, then params.

Key	Accepts	Notes
`k`, `K`, `hill`	scalar or (n, n)	A scalar is written onto every edge
`basal`, `gamma`, `vmax`, `beta`, `gammaP`	scalar or (n,)
`ic`	scalar or (dim,)	`dim = n` or `2n`

Rules enforced after overrides: K > 0 and hill >= 1 on every edge, gamma > 0, gammaP > 0, all production parameters nonnegative, ic nonnegative. Violations raise GRNError with identifier grn:badParameter; unknown keys raise grn:unknownParameter.

Source genes follow vmax: if you override vmax and not basal, the constitutive level of every unregulated gene is updated to the new vmax. Override basal explicitly to decouple them.

# Symmetric toggle switch with a specific Hill coefficient and threshold
model = grn.build_model("toggle", params=dict(hill=3, K=1.5))

# Heterogeneous degradation, one value per gene
model = grn.build_model(S, params=dict(gamma=[0.5, 1, 2, 1]))

# A hand-built K matrix (only the entries on edges matter)
Kmat = np.ones((4, 4)); Kmat[1, 0] = 0.4; Kmat[2, 1] = 2.0
model = grn.build_model(S, params=dict(K=Kmat))

The same toggle switch with Hill coefficient 1 and 2. Without cooperativity the nullclines cross once and every trajectory ends in the same state; with cooperativity they cross three times and the two outer crossings are attractors.

The attributes of a model, printed at the prompt:

model = grn.build_model("repressilator")
model, model.k, model.K, model.gamma, model.state_labels

Output

Model(name='repressilator', n=3, stages=1, logic='additive', edges=3, dim=3)
k =
 [[ 0.  0. 10.]
 [10.  0.  0.]
 [ 0. 10.  0.]]
K =
 [[1. 1. 1.]
 [1. 1. 1.]
 [1. 1. 1.]]
gamma = [1. 1. 1.]
state_labels = ('lacI', 'tetR', 'cI')

Two-stage models

stages=2 adds beta and gammaP and doubles the state. Three ways to control the protein dynamics:

# proteins outlive their mRNA 3 to 10 times, protein level = mRNA level
model = grn.build_model(S, stages=2, gammap_range=(0.1, 0.3))

# proteins two to five times more abundant than their mRNA
model = grn.build_model(S, stages=2, beta_ratio_range=(2, 5))

# explicit per-gene rates
model = grn.build_model(S, stages=2, params=dict(beta=2, gammaP=0.5))

Because thresholds K are drawn relative to the regulator's protein level, changing beta_ratio_range does not break the regulation. Changing beta or gammaP through params does not rescale K, so with the motif presets (which assume protein = mRNA) prefer beta_ratio_range=(1, 1), or set K yourself.

model = grn.build_model(grn.motif("cascade"), stages=2, seed=1, gammap_range=(0.1, 0.3))
on  = grn.stimulus("G1", "step", baseline=0, amplitude=1, start=2)
sol = grn.simulate(model, stimulus=on, initial_condition=np.zeros(8), duration=60)

A four-gene cascade with slow proteins after G1 is switched on at t = 2. Each protein lags its mRNA, and each stage of the cascade adds a further delay.

Editing a model by hand

A model is a dataclass with NumPy arrays. To change a parameter after building:

model.gamma[2] = 0.2
model.K[1, 0] = 0.8
model = grn.attach_handles(model)   # rebuild f, jac, decay, scale

attach_handles recomputes every derived attribute (dim, state_labels, decay, scale, hasReg, degrees) from the parameter arrays, so it is safe to call at any time. To keep the original, edit model.copy() instead.

Labels and names

labels names the genes (default "G1", "G2", ...). Labels are used for CSV column names, plots, and everywhere a gene may be given by name (perturbations, stimuli, screens). name sets model.name, which appears in figure titles; motif names are used automatically.

← Networks: motifs and random topologies Simulating expression data →