Regazzoni 2020 with scipy#
In this example we show how to use scipy.integrate.solve_ivp to solve the ODEs of the Regazzoni 2020 model.
We also compare the results with the circulation package’s built-in solver (which uses a forward euler scheme).
from circulation.log import setup_logging
from circulation.regazzoni2020 import Regazzoni2020
import matplotlib.pyplot as plt
setup_logging()
circulation = Regazzoni2020()
[12/09/25 10:36:51] INFO INFO:circulation.base: base.py:134 Circulation model parameters (Regazzoni2020) ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ┃ Parameter ┃ Value ┃ ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ │ HR │ 1.25 hertz │ │ chambers.LA.EA │ 0.07 millimeter_Hg / milliliter │ │ chambers.LA.EB │ 0.18 millimeter_Hg / milliliter │ │ chambers.LA.TC │ 0.17 second │ │ chambers.LA.TR │ 0.17 second │ │ chambers.LA.tC │ 0.9 second │ │ chambers.LA.V0 │ 4.0 milliliter │ │ chambers.LV.EA │ 4.482 millimeter_Hg / milliliter │ │ chambers.LV.EB │ 0.17 millimeter_Hg / milliliter │ │ chambers.LV.TC │ 0.25 second │ │ chambers.LV.TR │ 0.4 second │ │ chambers.LV.tC │ 0.1 second │ │ chambers.LV.V0 │ 42.0 milliliter │ │ chambers.RA.EA │ 0.06 millimeter_Hg / milliliter │ │ chambers.RA.EB │ 0.07 millimeter_Hg / milliliter │ │ chambers.RA.TC │ 0.17 second │ │ chambers.RA.TR │ 0.17 second │ │ chambers.RA.tC │ 0.9 second │ │ chambers.RA.V0 │ 4.0 milliliter │ │ chambers.RV.EA │ 0.2 millimeter_Hg / milliliter │ │ chambers.RV.EB │ 0.029 millimeter_Hg / milliliter │ │ chambers.RV.TC │ 0.25 second │ │ chambers.RV.TR │ 0.4 second │ │ chambers.RV.tC │ 0.1 second │ │ chambers.RV.V0 │ 16.0 milliliter │ │ valves.MV.Rmin │ 0.0075 millimeter_Hg * second / milliliter │ │ valves.MV.Rmax │ 75006.2 millimeter_Hg * second / milliliter │ │ valves.AV.Rmin │ 0.0075 millimeter_Hg * second / milliliter │ │ valves.AV.Rmax │ 75006.2 millimeter_Hg * second / milliliter │ │ valves.TV.Rmin │ 0.0075 millimeter_Hg * second / milliliter │ │ valves.TV.Rmax │ 75006.2 millimeter_Hg * second / milliliter │ │ valves.PV.Rmin │ 0.0075 millimeter_Hg * second / milliliter │ │ valves.PV.Rmax │ 75006.2 millimeter_Hg * second / milliliter │ │ circulation.SYS.R_AR │ 0.733 millimeter_Hg * second / milliliter │ │ circulation.SYS.C_AR │ 1.372 milliliter / millimeter_Hg │ │ circulation.SYS.R_VEN │ 0.32 millimeter_Hg * second / milliliter │ │ circulation.SYS.C_VEN │ 11.363 milliliter / millimeter_Hg │ │ circulation.SYS.L_AR │ 0.005 millimeter_Hg * second ** 2 / milliliter │ │ circulation.SYS.L_VEN │ 0.0005 millimeter_Hg * second ** 2 / milliliter │ │ circulation.PUL.R_AR │ 0.046 millimeter_Hg * second / milliliter │ │ circulation.PUL.C_AR │ 20.0 milliliter / millimeter_Hg │ │ circulation.PUL.R_VEN │ 0.0015 millimeter_Hg * second / milliliter │ │ circulation.PUL.C_VEN │ 16.0 milliliter / millimeter_Hg │ │ circulation.PUL.L_AR │ 0.0005 millimeter_Hg * second ** 2 / milliliter │ │ circulation.PUL.L_VEN │ 0.0005 millimeter_Hg * second ** 2 / milliliter │ │ circulation.external.start_withdrawal │ 0.0 second │ │ circulation.external.end_withdrawal │ 0.0 second │ │ circulation.external.start_infusion │ 0.0 second │ │ circulation.external.end_infusion │ 0.0 second │ │ circulation.external.flow_withdrawal │ 0.0 milliliter / second │ │ circulation.external.flow_infusion │ 0.0 milliliter / second │ └───────────────────────────────────────┴─────────────────────────────────────────────────┘
INFO INFO:circulation.base: base.py:141 Circulation model initial states (Regazzoni2020) ┏━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ┃ State ┃ Value ┃ ┡━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ │ V_LA │ 87.183 milliliter │ │ V_LV │ 118.52 milliliter │ │ V_RA │ 86.833 milliliter │ │ V_RV │ 166.177 milliliter │ │ p_AR_SYS │ 87.675 millimeter_Hg │ │ p_VEN_SYS │ 35.898 millimeter_Hg │ │ p_AR_PUL │ 19.545 millimeter_Hg │ │ p_VEN_PUL │ 15.004 millimeter_Hg │ │ Q_AR_SYS │ 71.104 milliliter / second │ │ Q_VEN_SYS │ 94.039 milliliter / second │ │ Q_AR_PUL │ 94.084 milliliter / second │ │ Q_VEN_PUL │ 473.279 milliliter / second │ └───────────┴─────────────────────────────┘
from time import perf_counter
from scipy.integrate import solve_ivp
import numpy as np
time = np.linspace(0, 10, 1000)
t0 = perf_counter()
res = solve_ivp(
circulation.rhs,
[0, 10],
circulation.state,
t_eval=time,
method="RK45",
)
t1 = perf_counter()
fig, ax = plt.subplots(2, 1, sharex=True, sharey=True, figsize=(10, 5))
state_names = circulation.state_names()
var_names = circulation.var_names()
vars = circulation.update_static_variables(time, res.y)
ax[0].plot(time, vars[var_names.index("p_LV"), :], label="p_LV (numpy)")
ax[0].plot(time, vars[var_names.index("p_LA"), :], label="p_LA (numpy)")
ax[0].plot(time, res.y[state_names.index("p_AR_SYS"), :], label="p_AR_SYS (numpy)")
[<matplotlib.lines.Line2D at 0x7fadfaff0290>]
ax[1].plot(time, res.y[state_names.index("V_LA"), :], label="V_LA (numpy)")
ax[1].plot(time, res.y[state_names.index("V_LV"), :], label="V_LV (numpy)")
[<matplotlib.lines.Line2D at 0x7fadfa858c50>]
t2 = perf_counter()
history = circulation.solve(num_beats=10)
t3 = perf_counter()
circulation.print_info()
INFO INFO:circulation.base: base.py:508 Volumes ┏━━━━━━━━┳━━━━━━━━━┳━━━━━━━━┳━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━┳━━━━━━━━━┳━━━━━━━━━┳━━━━━━━━━━┓ ┃ V_LA ┃ V_LV ┃ V_RA ┃ V_RV ┃ V_AR_SYS ┃ V_VEN_SYS ┃ V_AR_PUL ┃ V_VEN_PUL ┃ Heart ┃ SYS ┃ PUL ┃ Total ┃ ┡━━━━━━━━╇━━━━━━━━━╇━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━┩ │ 87.204 │ 118.659 │ 86.768 │ 166.139 │ 120.380 │ 408.167 │ 390.770 │ 239.789 │ 458.770 │ 528.547 │ 630.559 │ 1617.876 │ └────────┴─────────┴────────┴─────────┴──────────┴───────────┴──────────┴───────────┴─────────┴─────────┴─────────┴──────────┘ Pressures ┏━━━━━━━━┳━━━━━━━━┳━━━━━━━┳━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┓ ┃ p_LA ┃ p_LV ┃ p_RA ┃ p_RV ┃ p_AR_SYS ┃ p_VEN_SYS ┃ p_AR_PUL ┃ p_VEN_PUL ┃ ┡━━━━━━━━╇━━━━━━━━╇━━━━━━━╇━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━┩ │ 14.938 │ 12.988 │ 5.801 │ 4.348 │ 87.740 │ 35.921 │ 19.538 │ 14.987 │ └────────┴────────┴───────┴───────┴──────────┴───────────┴──────────┴───────────┘ Flows ┏━━━━━━━━━┳━━━━━━━━┳━━━━━━━━━┳━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┓ ┃ Q_MV ┃ Q_AV ┃ Q_TV ┃ Q_PV ┃ Q_AR_SYS ┃ Q_VEN_SYS ┃ Q_AR_PUL ┃ Q_VEN_PUL ┃ ┡━━━━━━━━━╇━━━━━━━━╇━━━━━━━━━╇━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━┩ │ 257.845 │ -0.001 │ 191.446 │ -0.000 │ 71.163 │ 94.123 │ 94.382 │ 470.060 │ └─────────┴────────┴─────────┴────────┴──────────┴───────────┴──────────┴───────────┘
ax[0].plot(history["time"], history["p_LV"], linestyle="--", label="p_LV (orig)")
ax[0].plot(history["time"], history["p_LA"], linestyle="--", label="p_LA (orig)")
ax[0].plot(history["time"], history["p_AR_SYS"], linestyle="--", label="p_AR_SYS (orig)")
ax[0].legend()
ax[1].plot(history["time"], history["V_LV"], linestyle="--", label="V_LV (orig)")
ax[1].plot(history["time"], history["V_LA"], linestyle="--", label="V_LA (orig)")
ax[1].legend()
<matplotlib.legend.Legend at 0x7fadfa82b9e0>
fig.savefig("regazzoni2020_comp.png", dpi=300, bbox_inches="tight")
print("IVP solve time: ", t1 - t0)
print("Circulation solve time: ", t3 - t2)
IVP solve time: 4.469457665999926
Circulation solve time: 2.2226567620000424