MATLAB · chapter 4 of 13

Building models

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

model = grn.buildModel(S)                                 % defaults
model = grn.buildModel(S, Logic = "multiplicative", Seed = 3)
model = grn.buildModel("repressilator", Stages = 2)       % motif + preset
model = grn.buildModel(S, Params = struct("hill", 2, "gamma", 1))

How parameters are sampled

With rng(Seed) set, 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 option
`gamma_i`	uniform	`GammaRange = [0.5 1.5]`
`vmax_i`	uniform	`VmaxRange = [1 5]`
`gammaP_i`	uniform (used only for two stages)	`GammaPRange = [0.2 1]`
`beta_i`	`gammaP_i` times a uniform ratio	`BetaRatioRange = [1 1]`
`hill_ij`	uniform, on edges only	`HillRange = [1 4]`
`K_ij`	`xchar_j` times a uniform factor	`KRelRange = [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	`BasalFraction = [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 struct whose fields override the sampled values. Motif presets are applied first, then Params.

Field	Accepts	Notes
`k`, `K`, `hill`	scalar or n x n	A scalar is written onto every edge
`basal`, `gamma`, `vmax`, `beta`, `gammaP`	scalar or n x 1
`ic`	scalar or dim x 1	`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 grn:badParameter; unknown fields 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.buildModel("toggle", Params = struct("hill", 3, "K", 1.5));

% Heterogeneous degradation, one value per gene
model = grn.buildModel(S, Params = struct("gamma", [0.5; 1; 2; 1]));

% A hand-built K matrix (only the entries on edges matter)
Kmat = ones(4); Kmat(2, 1) = 0.4; Kmat(3, 2) = 2.0;
model = grn.buildModel(S, Params = struct("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 fields of a model, printed at the prompt:

model = grn.buildModel("repressilator")

Output

            name: "repressilator"
            type: "grn"
               n: 3
               S: [3×3 double]
           logic: "additive"
          stages: 1
               k: [3×3 double]
               K: [3×3 double]
            hill: [3×3 double]
           basal: [3×1 double]
           gamma: [3×1 double]
            vmax: [3×1 double]
            beta: []
          gammaP: []
              ic: [3×1 double]
          labels: ["lacI"    "tetR"    "cI"]
            seed: 0
    perturbation: [1×1 struct]
         options: [1×1 struct]
          hasReg: [3×1 double]
        indegree: [3×1 double]
       outdegree: [3×1 double]
             dim: 3
       stateGene: [1 2 3]
         species: ["protein"    "protein"    "protein"]
     stateLabels: ["lacI"    "tetR"    "cI"]
           decay: [3×1 double]
           scale: [3×1 double]
        stimulus: [0×0 struct]
               f: @(t,x)grn.internal.rhs(x,p)
             jac: @(t,x)grn.internal.jacobian(x,p)
      production: @(t,x)birthOnly(x,p)
          breaks: [1×0 double]

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.buildModel(S, Stages = 2, GammaPRange = [0.1, 0.3]);

% proteins two to five times more abundant than their mRNA
model = grn.buildModel(S, Stages = 2, BetaRatioRange = [2, 5]);

% explicit per-gene rates
model = grn.buildModel(S, Stages = 2, Params = struct("beta", 2, "gammaP", 0.5));

Because thresholds K are drawn relative to the regulator's protein level, changing BetaRatioRange 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 BetaRatioRange = [1 1], or set K yourself.

model = grn.buildModel(grn.motif("cascade"), Stages = 2, Seed = 1, GammaPRange = [0.1, 0.3]);
on  = grn.stimulus("G1", "step", Baseline = 0, Amplitude = 1, Start = 2);
sol = grn.simulate(model, Stimulus = on, InitialCondition = zeros(8, 1), 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 plain struct. To change a parameter after building:

model.gamma(3) = 0.2;
model.K(2, 1) = 0.8;
model = grn.internal.attachHandles(model);   % rebuild f, jac, decay, scale

attachHandles recomputes every derived field (dim, stateLabels, decay, scale, hasReg, degrees) from the parameter fields, so it is safe to call at any time.

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 and export metadata; motif names are used automatically.

← Networks: motifs and random topologies Simulating expression data →