modelbase tutorial¶
from __future__ import annotations
from importlib.metadata import version
from typing import Iterable
import matplotlib.pyplot as plt
import numpy as np
from assets.models import upper_glycolysis, phase_plane_model
from assets.utils import print_as_table
from modelbase.ode import DerivedStoichiometry, Model, Simulator, mca
from scipy.integrate import solve_ivp
for pkg in ("modelbase",):
print(f"{pkg:<10} {version(pkg)}")
def kelvin_from_celsius(t: float) -> float:
return t + 273.15
UserWarning: Assimulo not found, disabling sundials support.
modelbase 1.59.0
Motivation¶
Let's say you want to model the following chemical network
$$ \Large \varnothing \xrightarrow{v_0} S \xrightarrow{v_1} P \xrightarrow{v_2} \varnothing $$
To translate that into a set of differential equations we first write out the stoichiometric matrix $N$
| $v_0$ | $v_1$ | $v_2$ | |
|---|---|---|---|
| S | 1 | -1 | 0 |
| P | 0 | 1 | -1 |
which translates into
$$\begin{align*} \frac{dS}{dt} &= v_0 - v_1 \\ \frac{dP}{dt} &= v_1 - v_2 \\ \end{align*} $$
and then choose rate equations for each rate to get the flux vector $v$
$$\begin{align*} v_0 &= k_0 \\ v_1 &= k_1 * S \\ v_2 &= k_2 * P \\ \end{align*}$$
Then the system of ODEs is given by $\frac{d\text{Model}}{dt} = Nv$, which gives us
$$\begin{align*} \frac{dS}{dt} &= k_0 - k_1 * S \\ \frac{dP}{dt} &= k_1 * S - k_2 * P \\ \end{align*}$$
Since the stoichiometric matrix is sparse, writing it out completely would do a lot of unnecessary computations.
The easiest way to do that would be to create your model as a function that returns the derivatives dSdt and dPdt and then to integrate that using an ODE solver.
def model(
_: float, y: Iterable[float], k_in: float, k1: float, k_out: float
) -> Iterable[float]:
# Unpack state vector
s, p = y
# Flux vector
v_in = k_in
v1 = k1 * s
v_out = k_out * p
# Stoichiometric matrix
dSdt = v_in - v1
dPdt = v1 - v_out
return dSdt, dPdt
t_span = (0, 10)
t_eval = np.linspace(*t_span, 100)
integration = solve_ivp(
model,
t_span=t_span,
t_eval=t_eval,
y0=(0, 0),
args=(1, 1, 1),
)
fig, ax = plt.subplots()
ax.plot(integration.t, integration.y.T)
ax.set_xlabel("time / a.u.")
ax.set_ylabel("Concentration / a.u.")
ax.legend(["S", "P"], loc="upper left", bbox_to_anchor=(1.01, 1), borderaxespad=0)
plt.show()
Creating your first model¶
While this works well for small models, changing the model or creating model-variants requires you to re-write that function, probably introducing errors along the way.
What you want is to be able to have a modular interface to building your models.
modelbase supplies you with the Model object, which you can use to iteratively build your models.
Let's begin by defining rate functions.
Note that these should be general and re-usable whenever possible, to make your model clear to people reading it.
Try to give these functions names that are meaningful to your audience, e.g. a rate function k * s could be named proportional or mass-action.
def constant(k: float) -> float:
return k
def proportional(k: float, s: float) -> float:
return k * s
Create the model by
- adding compounds (variables) and parameters (constants).
- adding the three reactions using
add_reaction_from_args
The method add_reaction_from_args has the following arguments:
rate_name: The name you want to give your rate. Has to be uniquefunction: The python function you defined abovestoichiometry: A dictionary defining which compounds will be consumed (negative value) or produced (positive value) by the reactionargs: A list of all arguments that will be passed to the python function. The order is the order in which the function will receive the arguments!
def linear_chain_2cpds() -> Model:
m = Model()
m.add_compounds(["S", "P"])
m.add_parameters({"k_in": 1, "k_1": 1, "k_out": 1})
m.add_reaction_from_args(
rate_name="v0",
function=constant,
stoichiometry={"S": 1},
args=["k_in"],
)
m.add_reaction_from_args(
rate_name="v1",
function=proportional,
stoichiometry={"S": -1, "P": 1},
args=["k_1", "S"],
)
m.add_reaction_from_args(
rate_name="v2",
function=proportional,
stoichiometry={"P": -1},
args=["k_out", "P"],
)
return m
Note that we defined our model in a function.
This isn't strictly necessary, but is highly recommended, as it avoids problems of hidden state.
There are some legacy variants of building the stoichiometric matrix, which you might encounter.
Model.add_reactionModel.add_reaction_from_ratelawModel.add_rateModel.add_ratesModel.add_stoichiometryModel.add_stoichiometry_by_compoundModel.add_stoichiometriesModel.add_stoichiometries_by_compounds
Do note however, that add_reaction_from_args is the preferred way.
Simulating the model¶
Create a Simulator object by
- passing the model into it
- initialising the simulator with a dictionary containing the initial conditions
- simulate the model until time point
t_end = 10 - plot the result
If we just want to plot a single result, we can make use of the fact that initialise and simulate_and both return the Simulator object and we thus make an easy-to-read and beautiful pipeline.
fig, ax = (
Simulator(linear_chain_2cpds())
.initialise({"S": 0, "P": 0})
.simulate_and(t_end=10)
.plot(
xlabel="time / a.u.",
ylabel="concentration / a.u.",
title="Linear chain",
figure_kwargs={"figsize": (5, 3.5)},
)
)
plt.show()
If you require multiple plots I recommend to assign the result of simulate_and to a variable.
s = (
Simulator(linear_chain_2cpds())
.initialise({"S": 0, "P": 0})
.simulate_and(t_end=10) # stop assignment here
)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(7, 3))
s.plot(
xlabel="time / a.u.",
ylabel="concentration / a.u.",
title="Linear chain",
ax=ax1,
)
s.plot_phase_plane(
cpd1="S",
cpd2="P",
title="Linear chain phase plane",
ax=ax2,
)
plt.show()
There are a lot more plot functions to be explored.
Explore them - usually they will be sufficient for most of your needs.
print_as_table([i for i in dir(s) if i.startswith("plot")])
plot | plot_3d_trajectories plot_against_variable | plot_all plot_all_against_variable | plot_derivatives plot_derived | plot_derived_against_variable plot_flux_selection | plot_flux_selection_against_variable plot_fluxes | plot_fluxes_against_variable plot_fluxes_grid | plot_grid plot_log | plot_phase_plane plot_phase_space | plot_producing_and_consuming plot_selection | plot_selection_against_variable plot_semilog | plot_trajectories
Simulate a steady-state¶
t_ss, y_ss = (
Simulator(linear_chain_2cpds())
.initialise({"S": 0, "P": 0})
.simulate_to_steady_state()
)
print(y_ss)
[[1. 1.]]
Parameter scans¶
Very often you want to systematically investigate the steady-state concentrations depending on different values of a parameter.
For this you can use the Simulator.parameter_scan method.
y_scan = (
Simulator(linear_chain_2cpds())
.initialise({"S": 0, "P": 0})
.parameter_scan(parameter_name="k_in", parameter_values=np.linspace(0, 10, 11))
)
y_scan
| S | P | |
|---|---|---|
| 0.0 | 0.0 | 0.0 |
| 1.0 | 1.0 | 1.0 |
| 2.0 | 2.0 | 2.0 |
| 3.0 | 3.0 | 3.0 |
| 4.0 | 4.0 | 4.0 |
| 5.0 | 5.0 | 5.0 |
| 6.0 | 6.0 | 6.0 |
| 7.0 | 7.0 | 7.0 |
| 8.0 | 8.0 | 8.0 |
| 9.0 | 9.0 | 9.0 |
| 10.0 | 10.0 | 10.0 |
If you are also interested in the steady-state fluxes, you can use the Simulator.parameter_scan_with_fluxes method.
y_scan, v_scan = (
Simulator(linear_chain_2cpds())
.initialise({"S": 0, "P": 0})
.parameter_scan_with_fluxes(
parameter_name="k_in", parameter_values=np.linspace(0, 10, 11)
)
)
v_scan
| v0 | v1 | v2 | |
|---|---|---|---|
| 0.0 | 0.0 | 0.0 | 0.0 |
| 1.0 | 1.0 | 1.0 | 1.0 |
| 2.0 | 2.0 | 2.0 | 2.0 |
| 3.0 | 3.0 | 3.0 | 3.0 |
| 4.0 | 4.0 | 4.0 | 4.0 |
| 5.0 | 5.0 | 5.0 | 5.0 |
| 6.0 | 6.0 | 6.0 | 6.0 |
| 7.0 | 7.0 | 7.0 | 7.0 |
| 8.0 | 8.0 | 8.0 | 8.0 |
| 9.0 | 9.0 | 9.0 | 9.0 |
| 10.0 | 10.0 | 10.0 | 10.0 |
Phase plane analysis¶
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
s = Simulator(phase_plane_model()).initialise({"s1": 10, "s2": 0})
s.plot_trajectories(
cpd1="s1",
cpd2="s2",
y0={"s1": 0, "s2": 0},
cpd1_bounds=(0, 2),
cpd2_bounds=(0, 2),
n=20,
ax=ax,
)
for s1 in np.linspace(0, 1, 4):
for s2 in np.linspace(0, 2, 4):
(
s.initialise({"s1": s1, "s2": s2})
.simulate_and(t_end=1.5, steps=1000)
.plot_phase_plane("s1", "s2", ax=ax)
)
fig.tight_layout()
plt.show()
Metabolic control analysis¶
modelbase provides routines to calculate two important measures from metabolic control analysis
- elasticities, which are measured at any given state of the system
- response coefficients, which are measured at steady-state
Elasticities¶
You can get the elasticities of the compounds / substrates via mca.get_compound_elasticities_df and the ones for the parameters using mca.get_parameter_elasticities_df.
There are routines for different data structures as well, but usually DataFrames are the preferred choice.
mca.get_compound_elasticities_df(
upper_glycolysis(),
compounds=["GLC", "F6P"],
y={
"GLC": 0.3,
"G6P": 0.4,
"F6P": 0.5,
"FBP": 0.6,
"ATP": 0.4,
"ADP": 0.6,
},
)
| v1 | v2 | v3 | v4 | v5 | v6 | v7 | |
|---|---|---|---|---|---|---|---|
| GLC | 0.0 | 1.0 | -0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| F6P | 0.0 | 0.0 | 5.0 | 1.0 | 0.0 | 0.0 | 0.0 |
mca.get_parameter_elasticities_df(
upper_glycolysis(),
parameters=["k1", "k2"],
y={
"GLC": 0.3,
"G6P": 0.4,
"F6P": 0.5,
"FBP": 0.6,
"ATP": 0.4,
"ADP": 0.6,
},
)
| v1 | v2 | v3 | v4 | v5 | v6 | v7 | |
|---|---|---|---|---|---|---|---|
| k1 | 1.0 | 0.0 | -0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| k2 | 0.0 | 1.0 | -0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
Response coefficients¶
You can get the control and flux response coefficients using mca.get_response_coefficients_df.
There are routines for different data structures as well, but usually DataFrames are the preferred choice.
crc, frc = mca.get_response_coefficients_df(
upper_glycolysis(),
parameters=["k1", "k2", "k3", "k4", "k5", "k6", "k7"],
y={
"GLC": 0,
"G6P": 0,
"F6P": 0,
"FBP": 0,
"ATP": 0.5,
"ADP": 0.5,
},
)
Plotting coefficients¶
You can plot the response coefficients using mca.plot_coefficient_heatmap.
This function will also work with the elasticities.
_ = mca.plot_coefficient_heatmap(
crc,
title="Concentration response coefficient",
annotate=False,
)
Quite often you might want to plot multiple response coefficients.
For example, that might be both concentration and flux response coefficients.
You can do this using mca.plot_multiple, which can plot an arbitrary amount of heatmaps.
fig, axs = mca.plot_multiple(
[crc, frc],
titles=["Concentration response coefficients", "Flux response coefficients"],
annotate=False,
figsize=(8, 4),
)
Deeper dive¶
modelbase contains a number of fundamental building blocks
- parameters (constants)
- compounds (variables)
- reactions, which consist of
- rates
- stoichiometries
Frequently, one wants to build building blocks that depend on those fundamental blocks.
For example, the product of the gas constant R with the temperature T occurs quite often in thermodynamics.
Thus, one might just want to define R and T and use a derived item RT = R * T.
modelbase offers a number of those derived building blocks, which differ in when they are calculated.
- derived parameters are calculated from other parameters and thus before the model is integrated
- derived compounds (also called algebraic modules) are calculated from other variables and thus before any rates
- derived stoichiometries are calculated from other variables and thus before any rates
Derived parameters¶
Derived parameters depend only on other parameters.
They are defined using Model.add_derived_parameter(name, function, parameters).
They are calculated at every parameter change.
def derived_parameter_example() -> Model:
m = Model()
m.add_parameter("temperature_celsius", 25.0)
m.add_derived_parameter(
"temperature_kelvin",
kelvin_from_celsius,
parameters=["temperature_celsius"],
)
return m
derived_parameter_example().parameters
{'temperature_celsius': 25.0, 'temperature_kelvin': 298.15}
Derived compounds / algebraic modules¶
Algebraic modules can depend on parameters and compounds (dynamic variables).
They are calculated every time the model is called.
Algebraic modules can define multiple derived compounds.
As frequently only a single derived compound is required, there exists the shortcut function add_derived_compound, which doesn't require explicitly setting the derived compounds.
Thus, their function signatures are as follows
Model.add_derived_compound(name, function, args)
Model.add_algebraic_module_from_args(name, function, derived_compounds, args)
def derived_compound_example() -> Model:
m = Model()
m.add_compound("temperature_celsius")
m.add_derived_compound(
"temperature_kelvin",
kelvin_from_celsius,
args=["temperature_celsius"],
)
return m
derived_compound_example().get_derived_variables({"temperature_celsius": 25.0})
temperature_kelvin 298.15 Name: derived_variables, dtype: float64
def distribute(s: float) -> tuple[float, float]:
return s / 3, s * 2 / 3
def algebraic_module_example() -> Model:
m = Model()
m.add_compound("a")
m.add_algebraic_module_from_args(
"distribute",
distribute,
derived_compounds=["a1", "a2"],
args=["a"],
)
return m
algebraic_module_example().get_derived_variables({"a": 1.0})
a1 0.333333 a2 0.666667 Name: derived_variables, dtype: float64
Derived stoichiometries¶
Like derived parameters, derived stoichiometries only depend on parameters.
They are defined along with their respective reactions using the additional
derived_stoichiometry argument.
def diffusion(inside: float, outside: float, k: float) -> float:
return k * (outside - inside)
def div(x: float, y: float) -> float:
return x / y
def derived_stoichiometries_example() -> Model:
m = Model()
m.add_compounds(["x_in"])
m.add_parameters(
{
"x_base_stoichiometry": 1.0,
"cell_size": 1.0,
"x_out": 1.0,
"k_diffusion": 1.0,
}
)
m.add_reaction_from_args(
"rxn",
diffusion,
stoichiometry={},
derived_stoichiometry={
"x_in": DerivedStoichiometry(div, ["x_base_stoichiometry", "cell_size"])
},
args=["x_in", "x_out", "k_diffusion"],
)
return m
print(derived_stoichiometries_example().stoichiometries)
print(
derived_stoichiometries_example()
.update_parameter("cell_size", 2.0) # update cell size
.stoichiometries
)
{'rxn': {'x_in': 1.0}}
{'rxn': {'x_in': 0.5}}
Modularisation and composition¶
The best way to build larger models is to compose them from smaller building blocks.
That makes testing and analysing the building blocks a lot easier.
For larger scale modularisation I recommend our qtbmodels package which builds on top of modelbase.
For small scale modularisation however, you can also quickly build model variants by expanding and modifying an existing model.
To make it easy for others to see exactly what you changed, I recommend building them as functions that take other model functions as arguments.
There are two ways of doing this.
- implicitly, creating the old model inside the new one
- explicitly, passing the old model as an argument to the new one
def create_model_v2() -> Model:
m = create_model_v1()
m.add_parameter("k_1_rev", 0.5)
m.add_reaction_from_args(
"v1_rev", proportional, {"P": -1, "S": 1}, ["k_1_rev", "P"]
)
return m
def create_model_v2(m: Model) -> Model:
m.add_parameter("k_1_rev", 0.5)
m.add_reaction_from_args(
"v1_rev", proportional, {"P": -1, "S": 1}, ["k_1_rev", "P"]
)
return m
I would generally recommed the latter approach, as it allows you to re-use the changes you are making in the second model variant.
def create_model_v1() -> Model:
m = Model()
m.add_compounds(["S", "P"])
m.add_parameters({"k_in": 1, "k_1": 1, "k_out": 1})
m.add_reaction_from_args("v0", constant, {"S": 1}, ["k_in"])
m.add_reaction_from_args("v1", proportional, {"S": -1, "P": 1}, ["k_1", "S"])
m.add_reaction_from_args("v2", proportional, {"P": -1}, ["k_out", "P"])
return m
def create_model_v2(m: Model) -> Model:
m.add_parameter("k_1_rev", 0.5)
m.add_reaction_from_args(
"v1_rev", proportional, {"P": -1, "S": 1}, ["k_1_rev", "P"]
)
return m
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(7, 3), sharey=True)
_ = (
Simulator(create_model_v1())
.initialise({"S": 0, "P": 0})
.simulate_and(t_end=10)
.plot(
xlabel="time / a.u.",
ylabel="concentration / a.u.",
title="Linear chain v1",
ax=ax1,
)
)
_ = (
Simulator(create_model_v2(create_model_v1()))
.initialise({"S": 0, "P": 0})
.simulate_and(t_end=10)
.plot(
xlabel="time / a.u.",
ylabel="concentration / a.u.",
title="Linear chain v2",
ax=ax2,
)
)
plt.show()
Introspection & bug fixing¶
Not all the times model creation goes as smoothly as in the example.
To allow you to quickly find an error, modelbase contains introspection methods, to quickly find out where things went wrong.
We can check the constructed system, by taking a look at the stoichiometric matrix, the fluxes and the right hand side.
Let's start with the stoichiometric matrix, which modelbase returns as a dataframe.
Here you can see which reaction affects which compound by how much.
If a number here is unexpected, you know that you need to change the respective stoichiometry dictionary.
linear_chain_2cpds().get_stoichiometric_df()
| v0 | v1 | v2 | |
|---|---|---|---|
| P | 0.0 | 1.0 | -1.0 |
| S | 1.0 | -1.0 | 0.0 |
The next step is to check the fluxes of the system.
This is especially useful to see whether the direction of all reactions is as expected or to check if a reaction is overly large or small
linear_chain_2cpds().get_fluxes_df({"S": 1.0, "P": 0.5})
| v0 | v1 | v2 | |
|---|---|---|---|
| 0.0 | 1.0 | 1.0 | 0.5 |
You can also check the right hand side of the differential equation.
Since that is just the stoichiometric matrix multiplied by the fluxes this doesn't give any new information, but it is an easy way to check if all reactions are going in the correct direction.
linear_chain_2cpds().get_right_hand_side({"S": 1.0, "P": 0.5})
{'dSdt': np.float64(0.0), 'dPdt': np.float64(0.5)}
Some bugs only appear after a certain amount of time.
The simulator object allows checking the concentrations over time (results) and fluxes over time.
If you need to access the concentrations over time for something other than plotting, you can do that using .get_results_df()
(
Simulator(linear_chain_2cpds())
.initialise({"S": 0, "P": 0})
.simulate_and(t_end=10)
.get_results_df()
.head(5) # type: ignore
)
| S | P | |
|---|---|---|
| 0.00000 | 0.000000 | 0.000000 |
| 0.10101 | 0.096076 | 0.004771 |
| 0.20202 | 0.182922 | 0.017855 |
| 0.30303 | 0.261423 | 0.037612 |
| 0.40404 | 0.332383 | 0.062639 |
A just like with the model you can get the fluxes.
(
Simulator(linear_chain_2cpds())
.initialise({"S": 0, "P": 0})
.simulate_and(t_end=10)
.get_fluxes_df()
.head(5) # type: ignore
)
| v0 | v1 | v2 | |
|---|---|---|---|
| 0.00000 | 1.0 | 0.000000 | 0.000000 |
| 0.10101 | 1.0 | 0.096076 | 0.004771 |
| 0.20202 | 1.0 | 0.182922 | 0.017855 |
| 0.30303 | 1.0 | 0.261423 | 0.037612 |
| 0.40404 | 1.0 | 0.332383 | 0.062639 |
Sharp edges¶
Misc¶
Explicit time dependency¶
If a rate is explicily dependent on time, modelbase gives access to that using the time argument in any algebraic module or reaction.
def time_dependency() -> Model:
m = Model()
m.add_compound("x")
m.add_reaction_from_args(
"v1",
proportional,
{"x": -1},
["time", "x"],
)
return m
fig, ax = (
Simulator(time_dependency())
.initialise({"x": 1})
.simulate_and(t_end=10)
.plot(
xlabel="time / a.u.",
ylabel="amount / a.u.",
title="Time-dependent reaction",
figure_kwargs={"figsize": (5, 3)},
)
)
Readouts¶
Readouts are a small performance enhancement that essentially function like derived compounds, but cannot be used by anything else in the model.
Model().add_readout(name, function, args)
Generating source code¶
modelbase models can generate their own source code using Model().generate_model_source_code().
The source code generated is in canonical form instead of trying to reproduce the exact code used to generate the model.
print(linear_chain_2cpds().generate_model_source_code())
import math
import numpy as np
from modelbase.ode import Model, Simulator
def constant(k: float) -> float:
return k
def proportional(k: float, s: float) -> float:
return k * s
m = Model()
m.add_parameters(parameters={"k_in": 1, "k_1": 1, "k_out": 1})
m.add_compounds(compounds=["S", "P"])
m.add_rate(
rate_name="v0",
function=constant,
substrates=[],
products=["S"],
modifiers=[],
parameters=["k_in"],
reversible=False,
args=["k_in"],
)
m.add_rate(
rate_name="v1",
function=proportional,
substrates=["S"],
products=["P"],
modifiers=[],
parameters=["k_1"],
reversible=False,
args=["k_1", "S"],
)
m.add_rate(
rate_name="v2",
function=proportional,
substrates=["P"],
products=[],
modifiers=[],
parameters=["k_out"],
reversible=False,
args=["k_out", "P"],
)
m.add_stoichiometries(
rate_stoichiometries={"v0": {"S": 1}, "v1": {"S": -1, "P": 1}, "v2": {"P": -1}}
)