Reaction carousel

In [1]:

import numpy as np

from mxlpy import Model, Simulator, fit, fns, plot, unwrap
from mxlpy.carousel import Carousel, ReactionTemplate


def get_sir() -> Model:
    """Create a simple SIR model."""
    return (
        Model()
        .add_variables({"s": 0.9, "i": 0.1, "r": 0.0})
        .add_parameters({"beta": 0.2, "gamma": 0.1})
        .add_reaction(
            "infection",
            fns.mass_action_2s,
            args=["s", "i", "beta"],
            stoichiometry={"s": -1, "i": 1},
        )
        .add_reaction(
            "recovery",
            fns.mass_action_1s,
            args=["i", "gamma"],
            stoichiometry={"i": -1, "r": 1},
        )
    )


carousel = Carousel(
    get_sir(),
    {
        "infection": [
            ReactionTemplate(fn=fns.mass_action_2s, args=["s", "i", "beta"]),
            ReactionTemplate(
                fn=fns.michaelis_menten_2s,
                args=["s", "i", "beta", "km_bs", "km_bi"],
                additional_parameters={"km_bs": 0.1, "km_bi": 1.0},
            ),
        ],
        "recovery": [
            ReactionTemplate(fn=fns.mass_action_1s, args=["i", "gamma"]),
            ReactionTemplate(
                fn=fns.michaelis_menten_1s,
                args=["i", "gamma", "km_gi"],
                additional_parameters={"km_gi": 0.1},
            ),
        ],
    },
)

import numpy as np from mxlpy import Model, Simulator, fit, fns, plot, unwrap from mxlpy.carousel import Carousel, ReactionTemplate def get_sir() -> Model: """Create a simple SIR model.""" return ( Model() .add_variables({"s": 0.9, "i": 0.1, "r": 0.0}) .add_parameters({"beta": 0.2, "gamma": 0.1}) .add_reaction( "infection", fns.mass_action_2s, args=["s", "i", "beta"], stoichiometry={"s": -1, "i": 1}, ) .add_reaction( "recovery", fns.mass_action_1s, args=["i", "gamma"], stoichiometry={"i": -1, "r": 1}, ) ) carousel = Carousel( get_sir(), { "infection": [ ReactionTemplate(fn=fns.mass_action_2s, args=["s", "i", "beta"]), ReactionTemplate( fn=fns.michaelis_menten_2s, args=["s", "i", "beta", "km_bs", "km_bi"], additional_parameters={"km_bs": 0.1, "km_bi": 1.0}, ), ], "recovery": [ ReactionTemplate(fn=fns.mass_action_1s, args=["i", "gamma"]), ReactionTemplate( fn=fns.michaelis_menten_1s, args=["i", "gamma", "km_gi"], additional_parameters={"km_gi": 0.1}, ), ], }, )

Simulate carousel ensemble¶

In [2]:

carousel_time_course = carousel.time_course(np.linspace(0, 100, 101))
variables_by_model = carousel_time_course.get_variables_by_model()

fig, ax = plot.one_axes()
plot.line_mean_std(variables_by_model["s"].unstack().T, label="s", ax=ax)
plot.line_mean_std(variables_by_model["i"].unstack().T, label="i", ax=ax)
plot.line_mean_std(variables_by_model["r"].unstack().T, label="r", ax=ax)
ax.legend()
plot.show()

carousel_time_course = carousel.time_course(np.linspace(0, 100, 101)) variables_by_model = carousel_time_course.get_variables_by_model() fig, ax = plot.one_axes() plot.line_mean_std(variables_by_model["s"].unstack().T, label="s", ax=ax) plot.line_mean_std(variables_by_model["i"].unstack().T, label="i", ax=ax) plot.line_mean_std(variables_by_model["r"].unstack().T, label="r", ax=ax) ax.legend() plot.show()

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

 75%|███████▌  | 3/4 [00:00<00:00, 29.48it/s]

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

Fitting¶

In [3]:

data = unwrap(
    Simulator(get_sir().update_parameters({"beta": 0.3, "gamma": 0.15}))
    .simulate(100, steps=11)
    .get_result()
).variables

data.head()

data = unwrap( Simulator(get_sir().update_parameters({"beta": 0.3, "gamma": 0.15})) .simulate(100, steps=11) .get_result() ).variables data.head()

Out[3]:

	s	i	r
0.000000	0.900000	0.100000	0.000000
9.090909	0.590585	0.198774	0.210641
18.181818	0.344478	0.175340	0.480182
27.272727	0.238034	0.096976	0.664990
36.363636	0.197740	0.044539	0.757722

In [4]:

res = fit.carousel_time_course(
    carousel,
    p0={
        "beta": 0.1,
        "gamma": 0.1,
        # specific to reaction templates
        # "km_bi": 1.0,
    },
    data=data,
)

best = res.get_best_fit().model

fig, ax = plot.one_axes()
plot.lines(
    unwrap(Simulator(best).simulate(100).get_result()).variables,
    ax=ax,
)
plot.reset_prop_cycle(ax=ax)
plot.lines(data, linestyle="dashed", ax=ax, legend=False)
plot.show()

res = fit.carousel_time_course( carousel, p0={ "beta": 0.1, "gamma": 0.1, # specific to reaction templates # "km_bi": 1.0, }, data=data, ) best = res.get_best_fit().model fig, ax = plot.one_axes() plot.lines( unwrap(Simulator(best).simulate(100).get_result()).variables, ax=ax, ) plot.reset_prop_cycle(ax=ax) plot.lines(data, linestyle="dashed", ax=ax, legend=False) plot.show()

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

 25%|██▌       | 1/4 [00:06<00:20,  6.77s/it]

WARNING:mxlpy.fit:Minimisation failed.

100%|██████████| 4/4 [00:06<00:00,  1.74s/it]

In [5]:

best_fit = res.get_best_fit()

print(best_fit.best_pars)
print([rxn.fn.__name__ for rxn in best_fit.model.get_raw_reactions().values()])

best_fit = res.get_best_fit() print(best_fit.best_pars) print([rxn.fn.__name__ for rxn in best_fit.model.get_raw_reactions().values()])

{'beta': np.float64(0.30825358774362493), 'gamma': np.float64(0.03983034715118903)}
['mass_action_2s', 'michaelis_menten_1s']

In [ ]: