Reaction carousel
In [1]:
Copied!
import numpy as np
from mxlpy import Model, Simulator, fit, fns, plot
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
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]:
Copied!
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]
100%|██████████| 4/4 [00:00<00:00, 28.84it/s]
Fitting¶
In [3]:
Copied!
data = (
Simulator(get_sir().update_parameters({"beta": 0.3, "gamma": 0.15}))
.simulate(100, steps=11)
.get_result()
.unwrap_or_err()
).variables
data.head()
data = (
Simulator(get_sir().update_parameters({"beta": 0.3, "gamma": 0.15}))
.simulate(100, steps=11)
.get_result()
.unwrap_or_err()
).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]:
Copied!
res = fit.carousel_time_course(
carousel,
p0={
"beta": 0.1,
"gamma": 0.1,
# specific to reaction templates
# "km_bi": 1.0,
},
data=data,
minimizer=fit.LocalScipyMinimizer(),
)
best = res.get_best_fit().model
fig, ax = plot.one_axes()
plot.lines(
Simulator(best).simulate(100).get_result().unwrap_or_err().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,
minimizer=fit.LocalScipyMinimizer(),
)
best = res.get_best_fit().model
fig, ax = plot.one_axes()
plot.lines(
Simulator(best).simulate(100).get_result().unwrap_or_err().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]
lsoda-- warning..internal t (=r1) and h (=r2) are
such that in the machine, t + h = t on the next step
(h = step size). solver will continue anyway
in above, r1 = 0.4411060481182D+02 r2 = 0.1030860755355D-14
lsoda-- warning..internal t (=r1) and h (=r2) are
such that in the machine, t + h = t on the next step
(h = step size). solver will continue anyway
in above, r1 = 0.4411060481182D+02 r2 = 0.1030860755355D-14
lsoda-- warning..internal t (=r1) and h (=r2) are
such that in the machine, t + h = t on the next step
(h = step size). solver will continue anyway
in above, r1 = 0.4411060481182D+02 r2 = 0.1758599387501D-14
lsoda-- warning..internal t (=r1) and h (=r2) are
such that in the machine, t + h = t on the next step
(h = step size). solver will continue anyway
in above, r1 = 0.4411060481182D+02 r2 = 0.4396498468752D-15
lsoda-- warning..internal t (=r1) and h (=r2) are
such that in the machine, t + h = t on the next step
(h = step size). solver will continue anyway
in above, r1 = 0.4411060481182D+02 r2 = 0.4396498468752D-15
lsoda-- warning..internal t (=r1) and h (=r2) are
such that in the machine, t + h = t on the next step
(h = step size). solver will continue anyway
in above, r1 = 0.4411060481182D+02 r2 = 0.2146695402867D-15
lsoda-- warning..internal t (=r1) and h (=r2) are
such that in the machine, t + h = t on the next step
(h = step size). solver will continue anyway
in above, r1 = 0.4411060481182D+02 r2 = 0.2146695402867D-14
lsoda-- warning..internal t (=r1) and h (=r2) are
such that in the machine, t + h = t on the next step
(h = step size). solver will continue anyway
in above, r1 = 0.4411060481182D+02 r2 = 0.1341684626792D-15
lsoda-- warning..internal t (=r1) and h (=r2) are
such that in the machine, t + h = t on the next step
(h = step size). solver will continue anyway
in above, r1 = 0.4411060481182D+02 r2 = 0.1341684626792D-15
lsoda-- warning..internal t (=r1) and h (=r2) are
such that in the machine, t + h = t on the next step
(h = step size). solver will continue anyway
in above, r1 = 0.4411060481182D+02 r2 = 0.2683369253584D-15
lsoda-- above warning has been issued i1 times.
it will not be issued again for this problem
in above message, i1 = 10
25%|██▌ | 1/4 [00:00<00:02, 1.46it/s]
50%|█████ | 2/4 [00:01<00:01, 1.99it/s]
/home/runner/work/MxlPy/MxlPy/.venv/lib/python3.12/site-packages/scipy/integrate/_ivp/lsoda.py:161: UserWarning: lsoda: Repeated convergence failures (perhaps bad Jacobian or tolerances). solver._y, solver.t = integrator.run(
75%|███████▌ | 3/4 [00:25<00:11, 11.22s/it]
100%|██████████| 4/4 [00:25<00:00, 6.27s/it]
lsoda-- at t (=r1) and step size h (=r2), the
corrector convergence failed repeatedly
or with abs(h) = hmin l
in above, r1 = 0.4411060481182D+02 r2 = 0.2336252633388D-21
In [5]:
Copied!
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.3108356251760811), 'gamma': np.float64(0.04033498463024931)}
['mass_action_2s', 'michaelis_menten_1s']
In [ ]:
Copied!