Fitting

In [1]:

from __future__ import annotations

import matplotlib.pyplot as plt
import numpy as np

from example_models import get_linear_chain_2v
from mxlpy import Simulator, fit, make_protocol, plot, unwrap

from __future__ import annotations import matplotlib.pyplot as plt import numpy as np from example_models import get_linear_chain_2v from mxlpy import Simulator, fit, make_protocol, plot, unwrap

Fitting¶

Almost every model at some point needs to be fitted to experimental data to be validated.

mxlpy offers highly customisable local and global routines for fitting either time series or steady-states.

For this tutorial we are going to use the fit module to optimise our parameter values and the plot module to plot some results.

Let's get started!

Creating synthetic data¶

Normally, you would fit your model to experimental data.
Here, for the sake of simplicity, we will generate some synthetic data.

Checkout the basics tutorial if you need a refresher on building and simulating models.

In [2]:

# As a small trick, let's define a variable for the model function
# That way, we can re-use it all over the file and easily replace
# it with another model
model_fn = get_linear_chain_2v

res = unwrap(
    Simulator(model_fn())
    .update_parameters({"k1": 1.0, "k2": 2.0, "k3": 1.0})
    .simulate_time_course(np.linspace(0, 10, 101))
    .get_result()
).get_combined()

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

# As a small trick, let's define a variable for the model function # That way, we can re-use it all over the file and easily replace # it with another model model_fn = get_linear_chain_2v res = unwrap( Simulator(model_fn()) .update_parameters({"k1": 1.0, "k2": 2.0, "k3": 1.0}) .simulate_time_course(np.linspace(0, 10, 101)) .get_result() ).get_combined() fig, ax = plot.lines(res) ax.set(xlabel="time / a.u.", ylabel="Conc. & Flux / a.u.") plt.show()

Steady-states¶

For the steady-state fit we need two inputs:

1. the steady state data, which we supply as a pandas.Series
2. an initial parameter guess

The fitting routine will compare all data contained in that series to the model output.

Note that the data both contains concentrations and fluxes!

In [3]:

data = res.iloc[-1]
data.head()

data = res.iloc[-1] data.head()

Out[3]:

x     0.500000
y     1.000045
v1    1.000000
v2    1.000000
v3    1.000045
Name: 10.0, dtype: float64

In [4]:

fit_result = unwrap(
    fit.steady_state(
        model_fn(),
        p0={"k1": 1.038, "k2": 1.87, "k3": 1.093},
        data=res.iloc[-1],
        minimizer=fit.LocalScipyMinimizer(),
    )
)

fit_result.best_pars

fit_result = unwrap( fit.steady_state( model_fn(), p0={"k1": 1.038, "k2": 1.87, "k3": 1.093}, data=res.iloc[-1], minimizer=fit.LocalScipyMinimizer(), ) ) fit_result.best_pars

Out[4]:

{'k1': np.float64(1.0000155163152737),
 'k2': np.float64(2.000024696808468),
 'k3': np.float64(0.9999684416227076)}

If only some of the data is required, you can use a subset of it.
The fitting routine will only try to fit concentrations and fluxes contained in that series.

In [5]:

fit_result = unwrap(
    fit.steady_state(
        model_fn(),
        p0={"k1": 1.038, "k2": 1.87, "k3": 1.093},
        data=data.loc[["x", "y"]],
        minimizer=fit.LocalScipyMinimizer(),
    )
)
fit_result.best_pars

fit_result = unwrap( fit.steady_state( model_fn(), p0={"k1": 1.038, "k2": 1.87, "k3": 1.093}, data=data.loc[["x", "y"]], minimizer=fit.LocalScipyMinimizer(), ) ) fit_result.best_pars

Out[5]:

{'k1': np.float64(0.982944694802088),
 'k2': np.float64(1.965884849179406),
 'k3': np.float64(0.9829008077692202)}

Time course¶

For the time course fit we need again need two inputs

1. the time course data, which we supply as a pandas.DataFrame
2. an initial parameter guess

The fitting routine will create data at every time points specified in the DataFrame and compare all of them.

Other than that, the same rules of the steady-state fitting apply.

In [6]:

fit_result = unwrap(
    fit.time_course(
        model_fn(),
        p0={"k1": 1.038, "k2": 1.87, "k3": 1.093},
        data=res,
        minimizer=fit.LocalScipyMinimizer(),
    )
)

fit_result.best_pars

fit_result = unwrap( fit.time_course( model_fn(), p0={"k1": 1.038, "k2": 1.87, "k3": 1.093}, data=res, minimizer=fit.LocalScipyMinimizer(), ) ) fit_result.best_pars

Out[6]:

{'k1': np.float64(1.0000001257062412),
 'k2': np.float64(2.0000059119300104),
 'k3': np.float64(1.0000008172941233)}

Protocol time courses¶

Normally, you would fit your model to experimental data.
Here, again, for the sake of simplicity, we will generate some synthetic data.

In [7]:

protocol = make_protocol(
    [
        (1, {"k1": 1.0}),
        (1, {"k1": 2.0}),
        (1, {"k1": 1.0}),
    ]
)

res_protocol = unwrap(
    Simulator(model_fn())
    .update_parameters({"k1": 1.0, "k2": 2.0, "k3": 1.0})
    .simulate_protocol(
        protocol,
        time_points_per_step=10,
    )
    .get_result()
).get_combined()

fig, ax = plot.lines(res_protocol)
ax.set(xlabel="time / a.u.", ylabel="Conc. & Flux / a.u.")
plt.show()

protocol = make_protocol( [ (1, {"k1": 1.0}), (1, {"k1": 2.0}), (1, {"k1": 1.0}), ] ) res_protocol = unwrap( Simulator(model_fn()) .update_parameters({"k1": 1.0, "k2": 2.0, "k3": 1.0}) .simulate_protocol( protocol, time_points_per_step=10, ) .get_result() ).get_combined() fig, ax = plot.lines(res_protocol) ax.set(xlabel="time / a.u.", ylabel="Conc. & Flux / a.u.") plt.show()

For the protocol time course fit we need three inputs

1. an initial parameter guess
2. the time course data, which we supply as a pandas.DataFrame
3. the protocol, which we supply as a pandas.DataFrame

Note that the parameter given by the protocol cannot be fitted anymore

In [8]:

fit_result = unwrap(
    fit.protocol_time_course(
        model_fn(),
        p0={"k2": 1.87, "k3": 1.093},  # note that k1 is given by the protocol
        data=res_protocol,
        protocol=protocol,
        minimizer=fit.LocalScipyMinimizer(),
    )
)

fit_result.best_pars

fit_result = unwrap( fit.protocol_time_course( model_fn(), p0={"k2": 1.87, "k3": 1.093}, # note that k1 is given by the protocol data=res_protocol, protocol=protocol, minimizer=fit.LocalScipyMinimizer(), ) ) fit_result.best_pars

Out[8]:

{'k2': np.float64(2.000001827016672), 'k3': np.float64(0.9999986401044075)}

Ensemble fitting¶

mxlpy supports ensebmle fitting, which is a multi-model single data approach, where shared parameters will be applied to all models at the same time.

Here you supply an iterable of models instead of just one, otherwise the API stays the same.

In [9]:

ensemble_fit = fit.ensemble_steady_state(
    [
        model_fn(),
        model_fn(),
    ],
    data=res.iloc[-1],
    p0={"k1": 1.038, "k2": 1.87, "k3": 1.093},
    minimizer=fit.LocalScipyMinimizer(tol=1e-6),
)

ensemble_fit = fit.ensemble_steady_state( [ model_fn(), model_fn(), ], data=res.iloc[-1], p0={"k1": 1.038, "k2": 1.87, "k3": 1.093}, minimizer=fit.LocalScipyMinimizer(tol=1e-6), )

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

 50%|█████     | 1/2 [00:01<00:01,  1.29s/it]

100%|██████████| 2/2 [00:01<00:00,  1.47it/s]

To get the best fitting model, you can use get_best_fit on the ensemble fit

In [10]:

fit_result = ensemble_fit.get_best_fit()

fit_result = ensemble_fit.get_best_fit()

And you can of course also access all other fits

In [11]:

[i.loss for i in ensemble_fit.fits]

[i.loss for i in ensemble_fit.fits]

Out[11]:

[np.float64(1.6611309860488136e-05), np.float64(1.6611309860488136e-05)]

Time course¶

Time course fits are adjusted just the same

In [12]:

ensemble_fit = fit.ensemble_time_course(
    [
        model_fn(),
        model_fn(),
    ],
    data=res,
    p0={"k1": 1.038, "k2": 1.87, "k3": 1.093},
    minimizer=fit.LocalScipyMinimizer(tol=1e-6),
)

ensemble_fit = fit.ensemble_time_course( [ model_fn(), model_fn(), ], data=res, p0={"k1": 1.038, "k2": 1.87, "k3": 1.093}, minimizer=fit.LocalScipyMinimizer(tol=1e-6), )

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

 50%|█████     | 1/2 [00:03<00:03,  3.08s/it]

100%|██████████| 2/2 [00:03<00:00,  1.36s/it]

100%|██████████| 2/2 [00:03<00:00,  1.68s/it]

Protocol time course¶

As are protocol time courses

In [13]:

ensemble_fit = fit.ensemble_protocol_time_course(
    [
        model_fn(),
        model_fn(),
    ],
    data=res_protocol,
    protocol=protocol,
    p0={"k2": 1.87, "k3": 1.093},  # note that k1 is given by the protocol
    minimizer=fit.LocalScipyMinimizer(tol=1e-6),
)

ensemble_fit = fit.ensemble_protocol_time_course( [ model_fn(), model_fn(), ], data=res_protocol, protocol=protocol, p0={"k2": 1.87, "k3": 1.093}, # note that k1 is given by the protocol minimizer=fit.LocalScipyMinimizer(tol=1e-6), )

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

 50%|█████     | 1/2 [00:01<00:01,  1.92s/it]

100%|██████████| 2/2 [00:02<00:00,  1.02s/it]

Joint fitting¶

Next, we support joint fitting, which is a combined multi-model multi-data approach, where shared parameters will be applied to all models at the same time

In [14]:

fit.joint_steady_state(
    [
        fit.FitSettings(model=model_fn(), data=res.iloc[-1]),
        fit.FitSettings(model=model_fn(), data=res.iloc[-1]),
    ],
    p0={"k1": 1.038, "k2": 1.87, "k3": 1.093},
    minimizer=fit.LocalScipyMinimizer(tol=1e-6),
)

fit.joint_steady_state( [ fit.FitSettings(model=model_fn(), data=res.iloc[-1]), fit.FitSettings(model=model_fn(), data=res.iloc[-1]), ], p0={"k1": 1.038, "k2": 1.87, "k3": 1.093}, minimizer=fit.LocalScipyMinimizer(tol=1e-6), )

Out[14]:

JointFitResult(
  best_pars={
    'k1': np.float64(1.0000149372806875),
    'k2': np.float64(2.000025722679901),
    'k3': np.float64(0.9999698102562694)
  },
  loss=np.float64(3.316840676386697e-05)
)

In [15]:

fit.joint_time_course(
    [
        fit.FitSettings(model=model_fn(), data=res),
        fit.FitSettings(model=model_fn(), data=res),
    ],
    p0={"k1": 1.038, "k2": 1.87, "k3": 1.093},
    minimizer=fit.LocalScipyMinimizer(tol=1e-6),
)

fit.joint_time_course( [ fit.FitSettings(model=model_fn(), data=res), fit.FitSettings(model=model_fn(), data=res), ], p0={"k1": 1.038, "k2": 1.87, "k3": 1.093}, minimizer=fit.LocalScipyMinimizer(tol=1e-6), )

Out[15]:

JointFitResult(
  best_pars={
    'k1': np.float64(0.9999993178316183),
    'k2': np.float64(1.9999961723337583),
    'k3': np.float64(1.0000008116982133)
  },
  loss=np.float64(1.8228566470869358e-06)
)

In [16]:

fit.joint_protocol_time_course(
    [
        fit.FitSettings(model=model_fn(), data=res_protocol, protocol=protocol),
        fit.FitSettings(model=model_fn(), data=res_protocol, protocol=protocol),
    ],
    p0={"k2": 1.87, "k3": 1.093},
    minimizer=fit.LocalScipyMinimizer(tol=1e-6),
)

fit.joint_protocol_time_course( [ fit.FitSettings(model=model_fn(), data=res_protocol, protocol=protocol), fit.FitSettings(model=model_fn(), data=res_protocol, protocol=protocol), ], p0={"k2": 1.87, "k3": 1.093}, minimizer=fit.LocalScipyMinimizer(tol=1e-6), )

Out[16]:

JointFitResult(
  best_pars={'k2': np.float64(2.0000014942851934), 'k3': np.float64(0.9999997359378191)},
  loss=np.float64(6.838379642522272e-07)
)

Mixed joint fitting¶

Lastly, we support mixed-joint fitting, where each analysis takes it's own residual function to allow fitting both time series and steady-state data for multiple models at the same time.

In [17]:

fit.joint_mixed(
    [
        fit.MixedSettings(
            model=model_fn(),
            data=res.iloc[-1],
            residual_fn=fit.steady_state_residual,
        ),
        fit.MixedSettings(
            model=model_fn(),
            data=res,
            residual_fn=fit.time_course_residual,
        ),
        fit.MixedSettings(
            model=model_fn(),
            data=res_protocol,
            protocol=protocol,
            residual_fn=fit.protocol_time_course_residual,
        ),
    ],
    p0={"k2": 1.87, "k3": 1.093},
    minimizer=fit.LocalScipyMinimizer(tol=1e-6),
)

fit.joint_mixed( [ fit.MixedSettings( model=model_fn(), data=res.iloc[-1], residual_fn=fit.steady_state_residual, ), fit.MixedSettings( model=model_fn(), data=res, residual_fn=fit.time_course_residual, ), fit.MixedSettings( model=model_fn(), data=res_protocol, protocol=protocol, residual_fn=fit.protocol_time_course_residual, ), ], p0={"k2": 1.87, "k3": 1.093}, minimizer=fit.LocalScipyMinimizer(tol=1e-6), )

Out[17]:

JointFitResult(
  best_pars={'k2': np.float64(1.9999979281880325), 'k3': np.float64(1.0000003839540792)},
  loss=np.float64(2.9641396939299067e-05)
)

First finish line

With that you now know most of what you will need from a day-to-day basis about fitting in mxlpy.

Congratulations!

Advanced topics / customisation¶

All fitting routines internally are build in a way that they will call a tree of functions.

• minimizer
  • residual_fn
    • integrator
    • loss_fn

You can therefore use dependency injection to overwrite the minimisation function, the loss function, the residual function and the integrator if need be.

In [18]:

from functools import partial
from typing import TYPE_CHECKING, cast

from mxlpy.integrators import Scipy

if TYPE_CHECKING:
    import pandas as pd

    from mxlpy.fit import LossFn, ResidualFn
    from mxlpy.model import Model
    from mxlpy.types import Array, IntegratorType

from functools import partial from typing import TYPE_CHECKING, cast from mxlpy.integrators import Scipy if TYPE_CHECKING: import pandas as pd from mxlpy.fit import LossFn, ResidualFn from mxlpy.model import Model from mxlpy.types import Array, IntegratorType

Parameterising scipy optimise¶

In [19]:

optimizer = fit.LocalScipyMinimizer(tol=1e-6, method="Nelder-Mead")

optimizer = fit.LocalScipyMinimizer(tol=1e-6, method="Nelder-Mead")

Custom loss function¶

You can change the loss function that is being passed to the minimsation function using the loss_fn keyword.
Depending on the use case (time course vs steady state) this function will be passed two pandas DataFrames or Series.

In [20]:

def mean_absolute_error(
    x: pd.DataFrame | pd.Series,
    y: pd.DataFrame | pd.Series,
) -> float:
    """Mean absolute error between two dataframes."""
    return cast(float, np.mean(np.abs(x - y)))


fit_result = unwrap(
    fit.time_course(
        model_fn(),
        p0={"k1": 1.038, "k2": 1.87, "k3": 1.093},
        data=res,
        loss_fn=mean_absolute_error,
        minimizer=fit.LocalScipyMinimizer(),
    )
)

fit_result.best_pars

def mean_absolute_error( x: pd.DataFrame | pd.Series, y: pd.DataFrame | pd.Series, ) -> float: """Mean absolute error between two dataframes.""" return cast(float, np.mean(np.abs(x - y))) fit_result = unwrap( fit.time_course( model_fn(), p0={"k1": 1.038, "k2": 1.87, "k3": 1.093}, data=res, loss_fn=mean_absolute_error, minimizer=fit.LocalScipyMinimizer(), ) ) fit_result.best_pars

Out[20]:

{'k1': np.float64(0.9999992365531291),
 'k2': np.float64(1.999991221520435),
 'k3': np.float64(0.9999812481484508)}

Custom integrator¶

You can change the default integrator to an integrator of your choice by partially application of the class of any of the existing ones.

Here, for example, we choose the Scipy solver suite and set the default relative and absolute tolerances to 1e-6 respectively.

In [21]:

fit_result = unwrap(
    fit.time_course(
        model_fn(),
        p0={"k1": 1.038, "k2": 1.87, "k3": 1.093},
        data=res,
        integrator=partial(Scipy, rtol=1e-6, atol=1e-6),
        minimizer=fit.LocalScipyMinimizer(),
    )
)

fit_result.best_pars

fit_result = unwrap( fit.time_course( model_fn(), p0={"k1": 1.038, "k2": 1.87, "k3": 1.093}, data=res, integrator=partial(Scipy, rtol=1e-6, atol=1e-6), minimizer=fit.LocalScipyMinimizer(), ) ) fit_result.best_pars

Out[21]:

{'k1': np.float64(1.0000001028767629),
 'k2': np.float64(2.0000056657202365),
 'k3': np.float64(1.0000007455902167)}

In [ ]: