Experimental features¶

Experimental features for modelbase2.

All APIs shown should be considered unstable and may change without notice.

In [1]:

from example_models import get_linear_chain_1v, get_linear_chain_2v
from modelbase2.experimental import (
    generate_model_code_py,
    generate_modelbase_code,
    model_diff,
    to_tex,
)

from example_models import get_linear_chain_1v, get_linear_chain_2v from modelbase2.experimental import ( generate_model_code_py, generate_modelbase_code, model_diff, to_tex, )

code generation¶

Currently, the limitation here is that functions used for reactions etc. cannot call other functions.

modelbase2 can generate own source code from a model.

In [2]:

print(generate_modelbase_code(get_linear_chain_1v()))

print(generate_modelbase_code(get_linear_chain_1v()))

from modelbase2 import Model

def constant(x: Float) -> Float:
    """Constant function."""
    return x

def mass_action_1s(s1: Float, k: Float) -> Float:
    """Irreversible mass action reaction with one substrate."""
    return k * s1

def create_model() -> Model:
    return (
        Model()
        .add_parameters({'k_in': 1.0, 'k_out': 1.0})
        .add_variables({'x': 1.0})
        .add_reaction(
                "v_in",
                fn=constant,
                args=['k_in'],
                stoichiometry={"x": 1},
            )
        .add_reaction(
                "v_out",
                fn=mass_action_1s,
                args=['k_out', 'x'],
                stoichiometry={"x": -1},
            )
    )

modelbase2 can also generate a generic python function from the source code.
The plan here is to generalise this to be able to export models into other programming languages as well.

In [3]:

print(generate_model_code_py(get_linear_chain_2v()))

print(generate_model_code_py(get_linear_chain_2v()))

from collections.abc import Iterable

from modelbase2.types import Float

def model(t: Float, y: Float) -> Iterable[Float]:
    x, y = y
    k1 = 1.0
    k2 = 2.0
    k3 = 1.0
    v1 = k1
    v2 = x * k2
    v3 = y * k3
    dxdt = v1 - v2
    dydt = v2 - v3
    return dxdt, dydt

Diffs¶

modelbase2 can generate diffs between two models to quickly analyse differences between all model elements.

In [4]:

print(
    model_diff(
        get_linear_chain_2v(),
        get_linear_chain_1v(),
    )
)

print( model_diff( get_linear_chain_2v(), get_linear_chain_1v(), ) )

Model Diff
----------
Missing Parameters: k2, k1, k3
Missing Variables: y
Missing Reactions: v2, v3, v1

In [5]:

print(
    model_diff(
        get_linear_chain_2v(),
        get_linear_chain_2v().scale_parameter("k3", 2.0),
    )
)

print( model_diff( get_linear_chain_2v(), get_linear_chain_2v().scale_parameter("k3", 2.0), ) )

Model Diff
----------
Different Parameters:
  k3: 1.0 != 2.0

LaTeX export¶

Currently, the limitation here is that functions used for reactions etc. cannot call other functions.

modelbase2 supports writing out your model as LaTeX.

In [6]:

print(to_tex(get_linear_chain_1v()))

print(to_tex(get_linear_chain_1v()))

\documentclass{article}
\usepackage[english]{babel}
\usepackage[a4paper,top=2cm,bottom=2cm,left=2cm,right=2cm,marginparwidth=1.75cm]{geometry}
\usepackage{amsmath, amssymb, array, booktabs,
            breqn, caption, longtable, mathtools, placeins,
            ragged2e, tabularx, titlesec, titling}
\newcommand{\sectionbreak}{\clearpage}
\setlength{\parindent}{0pt}

\title{Model construction}
\date{} % clear date
\author{modelbase}
\begin{document}
\maketitle
\FloatBarrier\subsubsection*{Variables}
\begin{longtable}{c|c}
    Model name & Initial concentration \\
    \hline
    \endhead
    x & 1.00e+00 \\
    \caption[Model variables]{Model variables}
    \label{table:model-vars}
\end{longtable}

\FloatBarrier\subsubsection*{Parameters}
\begin{longtable}{c|c}
    Parameter name & Parameter value \\
    \hline
    \endhead
    $\mathrm{k\_in}$ & 1.00e+00 \\
    $\mathrm{k\_out}$ & 1.00e+00 \\
    \caption[Model parameters]{Model parameters}
    \label{table:model-pars}
\end{longtable}

\FloatBarrier\subsubsection*{Reactions}
\begin{dmath*}
    v_in = \mathrm{k\_in}
\end{dmath*}

\begin{dmath*}
    v_out = x \cdot \mathrm{k\_out}
\end{dmath*}

\FloatBarrier\subsubsection*{Stoichiometries}
\begin{longtable}{c|c}
    Rate name & Stoichiometry \\
    \hline
    \endhead
    v\_in & $\mathrm{x}$ \\
    v\_out & $-\mathrm{x}$ \\
    \caption[Model stoichiometries.]{Model stoichiometries.}
    \label{table:model-stoichs}
\end{longtable}
\end{document}

Symbolic models & identifiability analysis¶

In [7]:

import sympy

from modelbase2.experimental import strikepy
from modelbase2.experimental.symbolic import SymbolicModel, to_symbolic_model
from modelbase2.model import Model


def check_identifiability(
    sym_model: SymbolicModel, outputs: list[sympy.Symbol]
) -> strikepy.Result:
    strike_model = strikepy.Model(
        states=list(sym_model.variables.values()),
        pars=list(sym_model.parameters.values()),
        eqs=sym_model.eqs,
        outputs=outputs,
    )
    return strikepy.strike_goldd(strike_model)


def infect(s: float, i: float, n: float, beta: float) -> float:
    return beta / n * i * s


def recover(i: float, gamma: float) -> float:
    return gamma * i


def total_population(s: float, i: float, r: float) -> float:
    return s + i + r


def sir() -> Model:
    return (
        Model()
        .add_parameters({"beta": 1.0, "gamma": 0.1})
        .add_variables({"s": 99, "i": 1, "r": 0})
        .add_derived("n", total_population, args=["s", "i", "r"])
        .add_reaction(
            "infect",
            infect,
            args=["s", "i", "n", "beta"],
            stoichiometry={"s": -1, "i": 1},
        )
        .add_reaction(
            "recover", recover, args=["i", "gamma"], stoichiometry={"i": -1, "r": 1}
        )
    )


model = sir()
sym_model = to_symbolic_model(model)
res = check_identifiability(sym_model, [sympy.Symbol("i"), sympy.Symbol("r")])
print(res.summary())

import sympy from modelbase2.experimental import strikepy from modelbase2.experimental.symbolic import SymbolicModel, to_symbolic_model from modelbase2.model import Model def check_identifiability( sym_model: SymbolicModel, outputs: list[sympy.Symbol] ) -> strikepy.Result: strike_model = strikepy.Model( states=list(sym_model.variables.values()), pars=list(sym_model.parameters.values()), eqs=sym_model.eqs, outputs=outputs, ) return strikepy.strike_goldd(strike_model) def infect(s: float, i: float, n: float, beta: float) -> float: return beta / n * i * s def recover(i: float, gamma: float) -> float: return gamma * i def total_population(s: float, i: float, r: float) -> float: return s + i + r def sir() -> Model: return ( Model() .add_parameters({"beta": 1.0, "gamma": 0.1}) .add_variables({"s": 99, "i": 1, "r": 0}) .add_derived("n", total_population, args=["s", "i", "r"]) .add_reaction( "infect", infect, args=["s", "i", "n", "beta"], stoichiometry={"s": -1, "i": 1}, ) .add_reaction( "recover", recover, args=["i", "gamma"], stoichiometry={"i": -1, "r": 1} ) ) model = sir() sym_model = to_symbolic_model(model) res = check_identifiability(sym_model, [sympy.Symbol("i"), sympy.Symbol("r")]) print(res.summary())

Main loop: 0it [00:00, ?it/s]

Main loop: 2it [00:00, 12.01it/s]

Summary
=======
The model is not FISPO.
Identifiable parameters: [beta, gamma]
Unidentifiable parameters: []
Identifiable variables: [s, i, r]
Unidentifiable variables: []
Identifiable inputs: []
Unidentifiable inputs: []