User guide#

The guide follows one model through detection and reduction: a precursor x1 is made at a basal rate and on a stimulus u, and is either converted into the observed product x2 or degraded. The system rests before the stimulus, which is switched on at \(t = 0\) with two doses.

The model#

symident.Model takes a reaction network with the species, the stoichiometry (one row per reaction, one column per species) and the rates:

import symident as si

m = si.Model(reactions={"species": ["x1", "x2", "u"],
                        "stoichiometry": [[1, 0, 0], [1, 0, 0], [-1, 1, 0],
                                          [-1, 0, 0], [0, -1, 0]],
                        "rates": ["kb", "kin*u", "k1*x1", "kl*x1", "k2*x2"]},
             observables={"y": "s*x2"},
             inputs=["u"],
             events=[{"var": "u", "time": 0, "value": "dose", "method": "replace"}],
             conditions={"rows": ["low", "high"], "cols": {"dose": [1, 10]}})

The five reactions are basal production of x1 (kb), production on the stimulus (kin*u), conversion into x2 (k1), degradation of x1 (kl) and degradation of x2 (k2). Only x2 is measured, up to the unknown scale s. u takes part in no reaction: it is an input, zero at rest, and the event sets it to dose at \(t = 0\). The two conditions differ in the dose.

The right-hand sides can also be given directly, as a dict from state to expression (^ or ** for powers), followed by the observables:

m = si.Model({"x1": "kb + kin*u - (k1 + kl)*x1", "x2": "k1*x1 - k2*x2", "u": "0"},
             {"y": "s*x2"},
             inputs=["u"],
             events=[{"var": "u", "time": 0, "value": "dose", "method": "replace"}],
             conditions={"rows": ["low", "high"], "cols": {"dose": [1, 10]}})

Arguments of Model:

inputs: input states, zero at rest.
events: a list of dicts with var, time, value and method (replace, add or multiply); root in place of time fires where an expression changes sign (see Switches).
conditions: one row per experimental condition; a numeric cell fixes a symbol in that condition, a string renames it. A dict {"rows": [...], "cols": {name: [values]}} or a pandas data frame.
trafo: a parameter transformation; an entry named like a state is its initial value.
fixed: known symbols.
parameters: further symbols treated as parameters.

Expressions#

Right-hand sides, observables, initial values and event values are built from +, -, *, /, powers (^ or **) and these functions:

Exponential and logarithm	`exp`, `exp10`, `b^x` with a positive rational `b`, `log`, `ln`, `log2`, `log10`, `log(x, b)`
Roots and powers	`sqrt`, rational powers such as `x^(1/2)`, free exponents `x^n` with a parameter `n`
Trigonometric	`sin`, `cos`, `tan`, `cot`, `sec`, `csc`, `asin`, `acos`, `atan`, `atan2`
Hyperbolic	`sinh`, `cosh`, `tanh`, `coth`, `sech`, `csch`, `asinh`, `acosh`
Switches	`piecewise`, `ifelse`, `Heaviside`, `sign`, `abs`, `max`, `min`, comparisons and logical operators (see Switches)

A logarithm, a root or a fractional power needs an argument that is positive on the coordinates declared positive. A power is rational only when written as a fraction: x^(1/2), not x^0.5. time is the clock and may appear anywhere. Any other function is rejected by the default engine with a message that names the term; engine="symbolic" accepts every SymPy function but is meant for small models.

A function not covered here can be requested from the maintainer, Simon Beyer (simon.beyer@fdm.uni-freiburg.de), or as an issue at dModverse/symident#issues.

Switches#

Right-hand sides and observables may switch between regimes: piecewise(v1, c1, v2, c2, ..., otherwise) in SBML order (0 where no condition holds and no otherwise is given), ifelse(c, a, b), the comparisons <, >, <=, >=, the operators &&, || and !, Heaviside(), sign(), abs(), max() and min(). A condition that the coordinates declared positive decide is resolved, so abs(x) is x and x > 0 is true for a positive x. Any other condition makes the model hybrid:

m = si.Model({"x": "piecewise(-k1*x, x > thr, -k2*x)"}, {"y": "c*x"})
print(si.detect(m, reconstruct=True))

Result:  rank 4 / 5  |  1 non-identifiable direction

Generators  X = Σᵢ η(i) ∂ᵢ

Scalings:
  X₁
      η(c)   :  - c
      η(thr) :  + thr
      η(x)   :  + x

Both regimes keep the scaling of x and c, and the threshold moves along with x, so of c, thr and the initial value only c*thr and c*x are determined. With a known threshold, x > 1 or fixed=["thr"], the switch breaks the scaling and the model is identifiable.

A direction is reported only when the observables of every regime stay unchanged along it, every switching surface whose switch shows in the observables is mapped to itself, and every event and the given initial values are kept. Every sign pattern of the switching functions counts as a regime, whether or not a trajectory reaches it; only a pattern that the declared signs contradict is dropped, such as x > c together with x < c - d for positive d. A parameter that acts only in a regime the trajectory never enters is therefore reported identifiable, while a reported direction always exists. Method states the rule and how it is decided.

A switch is visible when the observables differ between the regimes on its two sides on the surface. A switch that acts only on states that never reach the observables is invisible, and its surface may move:

m = si.Model({"x": "-k*x", "z": "piecewise(-a*z, z > c, -b*z)"}, {"y": "s*x"})
si.detect(m).info["switching"]["conditions"][0]["visibility"]
# {'c - z = 0': 'invisible'}: a, b, c and z are free, and so are s and x up to s*x

regimes names the regimes an experiment visits, by condition: each entry is a comparison or a list of them, the first the regime at the start. Then only these regimes and the surfaces between them count, and a parameter of a regime outside them is reported non-identifiable:

m = si.Model({"x": "piecewise(-k1*x, x > c, -k2*x)"}, {"y": "x"}, trafo={"x": "2"},
             regimes=["x > c"])

Here k2 and c are free. A start that the initial values contradict is an error; a start they do not decide is an assumption, listed in info["switching"]. Without regimes the result holds for every initial state; with it, for the named experiment.

A switch on time alone, such as time > tau, is a surface like the others: a symmetry keeps its switching time when the switch is visible. Events may also fire at a surface: an event with root in place of time resets its state each time the root changes sign,

si.Model({"x": "a - k*x"}, {"y": "x"}, trafo={"x": "0"},
         events=[{"var": "x", "root": "x - c", "value": "0", "method": "replace"}])

info["switching"] lists the surfaces, the regimes counted, the sign patterns excluded and, per condition, the regimes, the visibility of each surface and the assumptions.

The answer is certified from both sides. Necessary conditions, the observability rows of every regime from the start of every condition, from generic points of the visible surfaces and after a switch between two regimes, bound the rank from below; directions proven symbolically bound it from above. When the two bounds do not meet, detect() raises an error naming the candidates it could not prove, instead of returning a result. The same happens for x == c (a condition that holds only on a surface), for a switch on a switched expression, and for gauge, reduce_cq and free_initial together with switches. With equilibrate, every resting state of every regime must have a closed form.

Detection#

symident.detect() returns the non-identifiable directions as a Symmetries object. With equilibrate=True every condition starts at its steady state with the inputs at zero, here \(x_1^\ast = k_b/(k_1 + k_l)\) and \(x_2^\ast = k_1 x_1^\ast / k_2\):

res = si.detect(m, equilibrate=True, reduce_cq=True, reconstruct=True)
print(res)

Result:  rank 4 / 6  |  2 non-identifiable directions

Generators  X = Σᵢ η(i) ∂ᵢ

Scalings:
  X₁
      η(kb)  :  + kb
      η(kin) :  + kin
      η(s)   :  - s

General:
  X₂
      η(k1)  :  + k1
      η(kb)  :  - kb
      η(kin) :  - kin
      η(kl)  :  - k1

Of the six parameters only four combinations are identifiable. A generator \(X = \sum_i \eta_i\,\partial_{z_i}\) is a direction along which the observables do not change:

\(X_1\) raises both production rates and lowers the scale s by the same factor. x2 grows in proportion, and s*x2 stays the same: the absolute production rates are not identifiable, only kb*s and kin*s.
\(X_2\) shifts flux between the two fates of x1: it raises k1, lowers kl by the same amount and lowers the production rates. The total outflow k1 + kl is unchanged, so x1 relaxes at the same rate, and so is the flux into x2, k1*kb/(k1 + kl) at rest and k1*kin/(k1 + kl) per dose. x2 sees only these combinations, never the split itself.

Options of detect():

reconstruct: return general directions as exact rational functions (otherwise by their support).
equilibrate: start every condition at its steady state, with the inputs at zero.
reduce_cq: report moiety freedoms on conserved totals instead of single species.
gauge: False reconstructs every direction; None fixes one coordinate per scaling, chosen to keep the general directions narrow; a list ranks the coordinates to fix.
cores: threads and worker processes.

summary() adds how the result was obtained: the Lie order, whether the saturation is certified and the gauge. The raw result is available with to_dict(), each direction as a Symmetry.

Reduction#

symident.reduce() (or res.reduce()) returns a reparametrisation that removes the directions. A scaling is removed by fixing one of its coordinates; a general direction by its invariants, the combinations that stay constant along it, and each invariant becomes a new parameter q_1, q_2, ...:

red = res.reduce()
print(red)

Reduced 2 of 2 directions.

Trafo (non-identity entries)
  k1  = q_1
  kb  = q_2
  kin = q_3
  kl  = 0
  s   = 1

Invariants:
  q_1 = k1 + kl
  q_2 = k1*kb/(k1 + kl)
  q_3 = k1*kin/(k1 + kl)

Sections:
  {X₁}  s = 1
  {X₂}  kl = 0
  summary(detailed=True) explains each section.

The reduced model fixes the scale at s = 1 and moves all outflow of x1 into the conversion, kl = 0. Its parameters are the identifiable combinations found above: q_1 is the total outflow rate of x1, q_2 and q_3 are the basal and the stimulated flux into x2, each in units of the measured signal. k2 is identifiable and stays as it is.

A section is the equation that picks one point on each orbit. red.summary(detailed=True) explains each section.

red.trafo holds the transformation as a dict of strings, ready to be composed with the model’s parameters:

red.trafo
# {'k1': 'q_1', 'k2': 'k2', 'kb': 'q_2', 'kin': 'q_3', 'kl': '0', 's': '1'}

Settings#

symident.reconst_control() collects the settings of the saturation and the reconstruction, among them the degree bounds of the fits, the relevance caps and a time limit; pass the result as control to detect(). Under equilibrate, directions that stay open are last fitted along parameter lines and may then name the resting states; homotopy=False skips that route.