In [1]:

from __future__ import annotations

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from mxlpy import (
    InitialAssignment,
    SteadyStateModelBuilder,
    cartesian_product,
    fit,
    mca,
    plot,
    scan,
)
from mxlpy.meta import generate_model_code_mxlweb

from __future__ import annotations import matplotlib.pyplot as plt import numpy as np import pandas as pd from mxlpy import ( InitialAssignment, SteadyStateModelBuilder, cartesian_product, fit, mca, plot, scan, ) from mxlpy.meta import generate_model_code_mxlweb

Steady-state models¶

A SteadyStateModelBuilder describes a model without any dynamics: there are no variables, only parameters and quantities derived from them. This makes it a good fit for algebraic models, e.g. analytic steady-state solutions, response curves or allometric relationships.

Because there is nothing to integrate, a SteadyStateModelBuilder is never simulated. Instead it is evaluated directly, which is both simpler and much faster.

As an example we use an enzyme following Michaelis-Menten kinetics, which is competitively inhibited by \(I\)

\[ v = \frac{V_{\max} \cdot S}{K_m \left(1 + \frac{I}{K_i}\right) + S} \qquad V_{\max} = k_{cat} \cdot E_{tot}\]

Building the model¶

Parameters and derived quantities work just like in the other builders. Parameters can also be initial assignments, and derived quantities can depend on each other in any order.

In [2]:

def mul(a: float, b: float) -> float:
    return a * b


def apparent_km(km: float, i: float, ki: float) -> float:
    return km * (1 + i / ki)


def michaelis_menten(s: float, vmax: float, km: float) -> float:
    return vmax * s / (km + s)


def ratio(a: float, b: float) -> float:
    return a / b


def get_model() -> SteadyStateModelBuilder:
    return (
        SteadyStateModelBuilder()
        .add_parameters(
            {"kcat": 10.0, "e_tot": 0.5, "km": 2.0, "ki": 1.0, "S": 1.0, "I": 0.0}
        )
        .add_parameter("vmax", InitialAssignment(fn=mul, args=["kcat", "e_tot"]))
        .add_derived("km_app", fn=apparent_km, args=["km", "I", "ki"])
        .add_derived("v", fn=michaelis_menten, args=["S", "vmax", "km_app"])
        .add_derived("saturation", fn=ratio, args=["v", "vmax"])
    )


model = get_model()
model.derived

def mul(a: float, b: float) -> float: return a * b def apparent_km(km: float, i: float, ki: float) -> float: return km * (1 + i / ki) def michaelis_menten(s: float, vmax: float, km: float) -> float: return vmax * s / (km + s) def ratio(a: float, b: float) -> float: return a / b def get_model() -> SteadyStateModelBuilder: return ( SteadyStateModelBuilder() .add_parameters( {"kcat": 10.0, "e_tot": 0.5, "km": 2.0, "ki": 1.0, "S": 1.0, "I": 0.0} ) .add_parameter("vmax", InitialAssignment(fn=mul, args=["kcat", "e_tot"])) .add_derived("km_app", fn=apparent_km, args=["km", "I", "ki"]) .add_derived("v", fn=michaelis_menten, args=["S", "vmax", "km_app"]) .add_derived("saturation", fn=ratio, args=["v", "vmax"]) ) model = get_model() model.derived

Out[2]:

	value	unit
km_app	\(km \left(\frac{I}{ki} + 1.0\right)\)
v	\(\frac{S vmax}{S + km_{app}}\)
saturation	\(\frac{v}{vmax}\)

Evaluation¶

evaluate returns all derived values as a pandas.Series. This replaces the Simulator for this builder.

In [3]:

model.evaluate()

model.evaluate()

Out[3]:

km_app        2.000000
v             1.666667
saturation    0.333333
dtype: float64

Optionally, the parameters (including initial assignments) can be included as well

In [4]:

model.evaluate(include_parameters=True)

model.evaluate(include_parameters=True)

Out[4]:

kcat          10.000000
e_tot          0.500000
km             2.000000
ki             1.000000
S              1.000000
I              0.000000
vmax           5.000000
km_app         2.000000
v              1.666667
saturation     0.333333
dtype: float64

The model is a normal fluent builder, so you can update parameters and evaluate again

In [5]:

model.update_parameter("I", 2.0).evaluate()

model.update_parameter("I", 2.0).evaluate()

Out[5]:

km_app        6.000000
v             0.714286
saturation    0.142857
dtype: float64

Parameter scans¶

scan.steady_state accepts a SteadyStateModelBuilder and evaluates it for every row of to_scan. It returns a SteadyStateScan holding one Evaluation per scanned value. Since there are no variables, the values are accessed via combined or get_args.

Nothing is integrated, so there is no integrator to configure and no y0. Time course and protocol scans are not available for this builder.

In [6]:

res = scan.steady_state(
    get_model(),
    to_scan=pd.DataFrame({"S": np.geomspace(0.1, 100, 21)}),
)

fig, ax = plot.lines(res.combined.loc[:, ["v"]])
ax.set(xscale="log", xlabel="substrate / a.u.", ylabel="rate / a.u.")
plt.show()

res = scan.steady_state( get_model(), to_scan=pd.DataFrame({"S": np.geomspace(0.1, 100, 21)}), ) fig, ax = plot.lines(res.combined.loc[:, ["v"]]) ax.set(xscale="log", xlabel="substrate / a.u.", ylabel="rate / a.u.") plt.show()

  0%|          | 0/21 [00:00<?, ?it/s]

100%|██████████| 21/21 [00:00<00:00, 8367.89it/s]

get_args can additionally include the parameters

In [7]:

res.get_args(include_parameters=True).head()

res.get_args(include_parameters=True).head()

Out[7]:

	kcat	e_tot	km	ki	S	I	vmax	km_app	v	saturation
S
0.100000	10.0	0.5	2.0	1.0	0.100000	0.0	5.0	2.0	0.238095	0.047619
0.141254	10.0	0.5	2.0	1.0	0.141254	0.0	5.0	2.0	0.329839	0.065968
0.199526	10.0	0.5	2.0	1.0	0.199526	0.0	5.0	2.0	0.453566	0.090713
0.281838	10.0	0.5	2.0	1.0	0.281838	0.0	5.0	2.0	0.617569	0.123514
0.398107	10.0	0.5	2.0	1.0	0.398107	0.0	5.0	2.0	0.830045	0.166009

Scanning multiple parameters works via cartesian_product, which yields a pandas.MultiIndex

In [8]:

res = scan.steady_state(
    get_model(),
    to_scan=cartesian_product(
        {
            "S": np.geomspace(0.1, 100, 10),
            "I": np.linspace(0, 5, 6),
        }
    ),
)

fig, ax, _ = plot.heatmap_from_2d_idx(res.combined, "v")
ax.set(xlabel="S", ylabel="I", title="rate v")
plt.show()

res = scan.steady_state( get_model(), to_scan=cartesian_product( { "S": np.geomspace(0.1, 100, 10), "I": np.linspace(0, 5, 6), } ), ) fig, ax, _ = plot.heatmap_from_2d_idx(res.combined, "v") ax.set(xlabel="S", ylabel="I", title="rate v") plt.show()

  0%|          | 0/60 [00:00<?, ?it/s]

100%|██████████| 60/60 [00:00<00:00, 8125.08it/s]

If an evaluation fails with an arithmetic error (e.g. a division by zero), the respective row is filled with NaN instead of failing the whole scan

In [9]:

scan.steady_state(
    get_model(),
    to_scan=pd.DataFrame({"ki": [1.0, 0.0, 2.0]}),
).combined

scan.steady_state( get_model(), to_scan=pd.DataFrame({"ki": [1.0, 0.0, 2.0]}), ).combined

  0%|          | 0/3 [00:00<?, ?it/s]

100%|██████████| 3/3 [00:00<00:00, 3950.68it/s]

Out[9]:

	km_app	v	saturation
ki
1.0	2.0	1.666667	0.333333
0.0	NaN	NaN	NaN
2.0	2.0	1.666667	0.333333

Fitting¶

fit.steady_state accepts a SteadyStateModelBuilder and fits its parameters to a pandas.Series of data. The data can contain any parameter or derived value.

Here we assume we measured the rate \(v\) and the saturation of the enzyme at \(S = 2\) and fit \(K_m\) and \(k_{cat}\).

In [10]:

fitted = fit.steady_state(
    get_model().update_parameter("S", 2.0),
    p0={"km": 1.0, "kcat": 5.0},
    data=pd.Series({"v": 1.25, "saturation": 0.25}),
    minimizer=fit.LocalScipyMinimizer(),
).unwrap_or_err()
fitted.best_pars

fitted = fit.steady_state( get_model().update_parameter("S", 2.0), p0={"km": 1.0, "kcat": 5.0}, data=pd.Series({"v": 1.25, "saturation": 0.25}), minimizer=fit.LocalScipyMinimizer(), ).unwrap_or_err() fitted.best_pars

Out[10]:

{'km': np.float64(6.0000083815204395), 'kcat': np.float64(10.000006243275802)}

The returned Fit contains a copy of the model, which can be evaluated with the best parameters

In [11]:

fitted.model.update_parameters(fitted.best_pars).evaluate()

fitted.model.update_parameters(fitted.best_pars).evaluate()

Out[11]:

km_app        6.000008
v             1.249999
saturation    0.250000
dtype: float64

Response coefficients¶

Response coefficients describe how much the derived values change if a parameter is perturbed. As the model is evaluated instead of simulated, these are just the (normalised) derivatives of evaluate(). For example, the rate \(v\) responds proportionally to \(k_{cat}\).

In [12]:

mca.response_coefficients(get_model(), to_scan=["kcat", "km", "S"]).combined

mca.response_coefficients(get_model(), to_scan=["kcat", "km", "S"]).combined

  0%|          | 0/3 [00:00<?, ?it/s]

 33%|███▎      | 1/3 [00:01<00:03,  1.75s/it]

100%|██████████| 3/3 [00:01<00:00,  1.68it/s]

Out[12]:

	kcat	km	S
km_app	0.0	1.000000	0.000000
v	1.0	-0.666667	0.666667
saturation	0.0	-0.666667	0.666667

Export¶

SteadyStateModelBuilder models can be saved to the native .mxl.json format with mxlpy.save / mxlpy.load, and exported to MxlWeb, which evaluates them in the browser.

In [13]:

print(generate_model_code_mxlweb(get_model()))

print(generate_model_code_mxlweb(get_model()))

import { SteadyStateModelBuilder } from "@computational-biology-aachen/mxlweb-core";
import { Add, Divide, Mul, Name, Num } from "@computational-biology-aachen/mxlweb-core/mathml";

export function initModel(): SteadyStateModelBuilder {
    return new SteadyStateModelBuilder()
      .addParameter("kcat", {
        value: 10.0,
        texName: 'kcat',
      })
      .addParameter("e_tot", {
        value: 0.5,
        texName: 'e\\_tot',
      })
      .addParameter("km", {
        value: 2.0,
        texName: 'km',
      })
      .addParameter("ki", {
        value: 1.0,
        texName: 'ki',
      })
      .addParameter("S", {
        value: 1.0,
        texName: 'S',
      })
      .addParameter("I", {
        value: 0.0,
        texName: 'I',
      })
      .addParameter("vmax", {
        value: new Mul([new Name("e_tot"), new Name("kcat")]),
        texName: 'vmax',
      })
      .addAssignment("km_app", {
        fn: new Mul([new Name("km"), new Add([new Num(1.0), new Divide([new Name("I"), new Name("ki")])])]),
        texName: 'km\\_app',
      })
      .addAssignment("v", {
        fn: new Divide([new Mul([new Name("S"), new Name("vmax")]), new Add([new Name("S"), new Name("km_app")])]),
        texName: 'v',
      })
      .addAssignment("saturation", {
        fn: new Divide([new Name("v"), new Name("vmax")]),
        texName: 'saturation',
      })
  }

Current limitations¶

The SteadyStateModelBuilder supports every routine whose question makes sense without dynamics: steady-state scans, fits (including group, ensemble and joint fits), Monte Carlo methods, response coefficients and sensitivity analysis. Time courses, protocols and everything about reactions or fluxes don't apply to it.