In [1]:

from __future__ import annotations

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

from mxlpy import (
    OdeModelBuilder,
    Simulator,
    cartesian_product,
    OdeRunner,
    fit,
    make_protocol,
    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 ( OdeModelBuilder, Simulator, cartesian_product, OdeRunner, fit, make_protocol, mca, plot, scan, ) from mxlpy.meta import generate_model_code_mxlweb

ODE models¶

mxlpy has three model builders, which differ in how the right hand side of the system is written down:

Builder	Right hand side	Typical use
KineticModelBuilder (Model)	\(N \cdot v(x)\): reactions with rate laws and a stoichiometric matrix	metabolic & signalling networks
OdeModelBuilder	\(\frac{dx}{dt} = f(x)\): one differential equation per variable	phenomenological models, equations taken from a paper
SteadyStateModelBuilder	no variables, only explicit functions of parameters	algebraic / steady-state models, see Steady-state models

This notebook covers the OdeModelBuilder. Use it when your model is naturally written as differential equations rather than as a reaction network. Since there are no reactions, an OdeModelBuilder has no fluxes and no stoichiometries.

As an example we use logistic growth with a constant harvesting rate \(h\)

\[\frac{dx}{dt} = r x \left(1 - \frac{x}{K}\right) - h x\]

which has the non-trivial steady state \(x^* = K \left(1 - \frac{h}{r}\right)\).

Building the model¶

Instead of add_reaction, every variable gets exactly one differential equation via add_diff_eq, which also sets the variable's initial value. Parameters and derived quantities work just like in the KineticModelBuilder.

As with all builders, use named functions instead of lambdas (they need to be inspectable and pickle-able).

In [2]:

def logistic_harvest(x: float, r: float, capacity: float, h: float) -> float:
    return r * x * (1 - x / capacity) - h * x


def fraction(x: float, capacity: float) -> float:
    return x / capacity


def get_model() -> OdeModelBuilder:
    return (
        OdeModelBuilder()
        .add_parameters({"r": 1.0, "capacity": 10.0, "h": 0.2})
        .add_diff_eq(
            "x",
            fn=logistic_harvest,
            args=["x", "r", "capacity", "h"],
            initial_value=0.5,
        )
        .add_derived("saturation", fn=fraction, args=["x", "capacity"])
    )


model = get_model()
model.diff_eqs

def logistic_harvest(x: float, r: float, capacity: float, h: float) -> float: return r * x * (1 - x / capacity) - h * x def fraction(x: float, capacity: float) -> float: return x / capacity def get_model() -> OdeModelBuilder: return ( OdeModelBuilder() .add_parameters({"r": 1.0, "capacity": 10.0, "h": 0.2}) .add_diff_eq( "x", fn=logistic_harvest, args=["x", "r", "capacity", "h"], initial_value=0.5, ) .add_derived("saturation", fn=fraction, args=["x", "capacity"]) ) model = get_model() model.diff_eqs

Out[2]:

	value	initial	unit
x	\(- h x + r x \left(1.0 - \frac{x}{capacity}\right)\)	0.5

In [3]:

model.derived

model.derived

Out[3]:

	value	unit
saturation	\(\frac{x}{capacity}\)

The right hand side can be evaluated directly, either at the initial conditions or at custom values

In [4]:

model.get_right_hand_side({"x": 5.0})

model.get_right_hand_side({"x": 5.0})

Out[4]:

x    1.5
dtype: float64

Initial values can be changed with update_variable / update_variables

In [5]:

get_model().update_variable("x", 2.0).get_initial_conditions()

get_model().update_variable("x", 2.0).get_initial_conditions()

Out[5]:

{'x': 2.0}

Simulation¶

An OdeModelBuilder is simulated with the same Simulator as a KineticModelBuilder. The result is an OdeSimulation, which offers the same access patterns as a Simulation, minus everything flux related.

Routines like scans and fits configure the simulation with an OdeRunner instead, see the integrator configuration.

In [6]:

res = Simulator(get_model()).simulate(t_end=15).get_result().unwrap_or_err()

fig, (ax1, ax2) = plot.two_axes(figsize=(7, 3))
plot.lines(res.variables, ax=ax1)
plot.lines(res.get_right_hand_side(), ax=ax2)
ax1.set(xlabel="time / a.u.", ylabel="value / a.u.")
ax2.set(xlabel="time / a.u.", ylabel="dx/dt / a.u.")
plt.show()

res = Simulator(get_model()).simulate(t_end=15).get_result().unwrap_or_err() fig, (ax1, ax2) = plot.two_axes(figsize=(7, 3)) plot.lines(res.variables, ax=ax1) plot.lines(res.get_right_hand_side(), ax=ax2) ax1.set(xlabel="time / a.u.", ylabel="value / a.u.") ax2.set(xlabel="time / a.u.", ylabel="dx/dt / a.u.") plt.show()

get_combined returns all variables, derived variables and readouts as a single pandas.DataFrame. get_args additionally allows including parameters and derived parameters.

In [7]:

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

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

Out[7]:

	x	r	capacity	h	saturation
14.141277	7.998535	1.0	10.0	0.2	0.799853
14.372871	7.998782	1.0	10.0	0.2	0.799878
14.604465	7.998988	1.0	10.0	0.2	0.799899
14.836058	7.999159	1.0	10.0	0.2	0.799916
15.000000	7.999263	1.0	10.0	0.2	0.799926

Steady states and protocols work as usual. For the default parameters the steady state is at \(x^* = 10 \cdot (1 - 0.2) = 8\)

In [8]:

Simulator(get_model()).simulate_to_steady_state().get_result().unwrap_or_err().get_new_y0()

Simulator(get_model()).simulate_to_steady_state().get_result().unwrap_or_err().get_new_y0()

Out[8]:

{'x': np.float64(8.0)}

In [9]:

protocol = make_protocol(
    [
        (10, {"h": 0.2}),
        (10, {"h": 0.6}),  # overharvesting
        (10, {"h": 0.2}),
    ]
)
res = (
    Simulator(get_model())
    .simulate_protocol(protocol, time_points_per_step=20)
    .get_result()
    .unwrap_or_err()
)

fig, ax = plot.lines(res.get_variables(include_derived_variables=False))
plot.shade_protocol(protocol["h"], ax=ax, alpha=0.1)
ax.set(xlabel="time / a.u.", ylabel="x / a.u.")
plt.show()

protocol = make_protocol( [ (10, {"h": 0.2}), (10, {"h": 0.6}), # overharvesting (10, {"h": 0.2}), ] ) res = ( Simulator(get_model()) .simulate_protocol(protocol, time_points_per_step=20) .get_result() .unwrap_or_err() ) fig, ax = plot.lines(res.get_variables(include_derived_variables=False)) plot.shade_protocol(protocol["h"], ax=ax, alpha=0.1) ax.set(xlabel="time / a.u.", ylabel="x / a.u.") plt.show()

Parameter scans¶

All functions of the scan module accept an OdeModelBuilder. Their results have the same shape as for a KineticModelBuilder, just without fluxes. The same goes for the Monte Carlo routines in mc.

In [10]:

res = scan.steady_state(
    get_model(),
    to_scan=pd.DataFrame({"h": np.linspace(0, 0.9, 10)}),
    runner=OdeRunner(tolerance=1e-8),
)

fig, ax = plot.lines(res.variables)
ax.set(xlabel="harvesting rate h / a.u.", ylabel="steady state / a.u.")
plt.show()

res = scan.steady_state( get_model(), to_scan=pd.DataFrame({"h": np.linspace(0, 0.9, 10)}), runner=OdeRunner(tolerance=1e-8), ) fig, ax = plot.lines(res.variables) ax.set(xlabel="harvesting rate h / a.u.", ylabel="steady state / a.u.") plt.show()

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

 10%|█         | 1/10 [00:02<00:21,  2.43s/it]

 20%|██        | 2/10 [00:03<00:10,  1.34s/it]

100%|██████████| 10/10 [00:03<00:00,  3.27it/s]

Scanning multiple parameters works via cartesian_product

In [11]:

res = scan.steady_state(
    get_model(),
    to_scan=cartesian_product(
        {
            "h": np.linspace(0, 0.9, 10),
            "capacity": np.linspace(5, 15, 5),
        }
    ),
)

fig, ax, _ = plot.heatmap_from_2d_idx(res.variables, "x")
ax.set(xlabel="h", ylabel="capacity", title="steady state of x")
plt.show()

res = scan.steady_state( get_model(), to_scan=cartesian_product( { "h": np.linspace(0, 0.9, 10), "capacity": np.linspace(5, 15, 5), } ), ) fig, ax, _ = plot.heatmap_from_2d_idx(res.variables, "x") ax.set(xlabel="h", ylabel="capacity", title="steady state of x") plt.show()

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

  2%|▏         | 1/50 [00:04<03:54,  4.78s/it]

 10%|█         | 5/50 [00:05<00:34,  1.32it/s]

 20%|██        | 10/50 [00:05<00:12,  3.12it/s]

 30%|███       | 15/50 [00:05<00:07,  4.97it/s]

 42%|████▏     | 21/50 [00:05<00:03,  7.59it/s]

100%|██████████| 50/50 [00:05<00:00,  8.70it/s]

The to_scan DataFrame can also contain initial conditions

In [12]:

res = scan.time_course(
    get_model(),
    to_scan=pd.DataFrame({"x": [0.1, 1.0, 5.0, 15.0]}),
    time_points=np.linspace(0, 10, 101),
)

fig, ax = plot.lines(res.get_by_name("x"))
ax.set(xlabel="time / a.u.", ylabel="x / a.u.")
plt.show()

res = scan.time_course( get_model(), to_scan=pd.DataFrame({"x": [0.1, 1.0, 5.0, 15.0]}), time_points=np.linspace(0, 10, 101), ) fig, ax = plot.lines(res.get_by_name("x")) ax.set(xlabel="time / a.u.", ylabel="x / a.u.") plt.show()

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

 25%|██▌       | 1/4 [00:01<00:05,  1.95s/it]

100%|██████████| 4/4 [00:01<00:00,  2.05it/s]

Fitting¶

All fitting routines except the reaction carousel accept an OdeModelBuilder. The returned Fit contains a copy of the OdeModelBuilder.

Here we create some synthetic data with known parameters and fit the model to it.

In [13]:

data = (
    Simulator(get_model().update_parameters({"r": 0.7, "h": 0.35}))
    .simulate_time_course(np.linspace(0, 15, 31))
    .get_result()
    .unwrap_or_err()
    .get_variables(
        include_derived_variables=False,
        include_readouts=False,
        include_surrogate_variables=False,
    )
)

fitted = fit.time_course(
    get_model(),
    p0={"r": 1.0, "h": 0.2},
    data=data,
    minimizer=fit.LocalScipyMinimizer(),
).unwrap_or_err()
fitted.best_pars

data = ( Simulator(get_model().update_parameters({"r": 0.7, "h": 0.35})) .simulate_time_course(np.linspace(0, 15, 31)) .get_result() .unwrap_or_err() .get_variables( include_derived_variables=False, include_readouts=False, include_surrogate_variables=False, ) ) fitted = fit.time_course( get_model(), p0={"r": 1.0, "h": 0.2}, data=data, minimizer=fit.LocalScipyMinimizer(), ).unwrap_or_err() fitted.best_pars

Out[13]:

{'r': np.float64(1.275587677108538), 'h': np.float64(0.8088565783084489)}

In [14]:

res = (
    Simulator(fitted.model.update_parameters(fitted.best_pars))
    .simulate_time_course(data.index)
    .get_result()
    .unwrap_or_err()
)

fig, ax = plot.lines(res.get_variables(include_derived_variables=False))
ax.scatter(data.index, data["x"], color="black", label="data")
ax.set(xlabel="time / a.u.", ylabel="x / a.u.")
ax.legend()
plt.show()

res = ( Simulator(fitted.model.update_parameters(fitted.best_pars)) .simulate_time_course(data.index) .get_result() .unwrap_or_err() ) fig, ax = plot.lines(res.get_variables(include_derived_variables=False)) ax.scatter(data.index, data["x"], color="black", label="data") ax.set(xlabel="time / a.u.", ylabel="x / a.u.") ax.legend() plt.show()

Response coefficients¶

Response coefficients describe how much the steady state changes if a parameter is perturbed. For an OdeModelBuilder they are available for the variables. With the default normalisation, the steady state \(x^* = K (1 - h / r)\) responds proportionally to the carrying capacity \(K\).

In [15]:

mca.response_coefficients(get_model()).variables

mca.response_coefficients(get_model()).variables

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

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

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

Out[15]:

	r	capacity	h
x	0.25	1.000000e+00	-0.25
saturation	0.25	1.630640e-09	-0.25

Export¶

OdeModelBuilder models can be saved to the native .mxl.json format with mxlpy.save / mxlpy.load, and exported to MxlWeb, which natively supports differential equation models.

In [16]:

print(generate_model_code_mxlweb(get_model()))

print(generate_model_code_mxlweb(get_model()))

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

export function initModel(): OdeModelBuilder {
    return new OdeModelBuilder()
      .addParameter("r", {
        value: 1.0,
        texName: 'r',
      })
      .addParameter("capacity", {
        value: 10.0,
        texName: 'capacity',
      })
      .addParameter("h", {
        value: 0.2,
        texName: 'h',
      })
      .addVariable("x", {
        value: 0.5,
        texName: 'x',
      })
      .addAssignment("saturation", {
        fn: new Divide([new Name("x"), new Name("capacity")]),
        texName: 'saturation',
      })
      .setDifferential("x", new Add([new Minus([new Mul([new Name("h"), new Name("x")])]), new Mul([new Name("r"), new Name("x"), new Add([new Num(1.0), new Minus([new Divide([new Name("x"), new Name("capacity")])])])])]))
  }

Current limitations¶

Everything about reactions, fluxes or stoichiometries only applies to the KineticModelBuilder. Not supported for the OdeModelBuilder are thus

• flux elasticities and rate characteristics in mca, reaction carousels, model comparisons and label models
• Jacobians (Simulator(..., use_jacobian=True) falls back to finite differences with a warning) and events