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()
[07/01/25 08:23:58] INFO INFO:circulation.base: base.py:132 Circulation model parameters (Regazzoni2020) ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ┃ Parameter ┃ Value ┃ ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ │ HR │ 1.0 hertz │ │ chambers.LA.EA │ 0.07 millimeter_Hg / milliliter │ │ chambers.LA.EB │ 0.09 millimeter_Hg / milliliter │ │ chambers.LA.TC │ 0.17 second │ │ chambers.LA.TR │ 0.17 second │ │ chambers.LA.tC │ 0.8 second │ │ chambers.LA.V0 │ 4.0 milliliter │ │ chambers.LV.EA │ 2.75 millimeter_Hg / milliliter │ │ chambers.LV.EB │ 0.08 millimeter_Hg / milliliter │ │ chambers.LV.TC │ 0.34 second │ │ chambers.LV.TR │ 0.17 second │ │ chambers.LV.tC │ 0.0 second │ │ chambers.LV.V0 │ 5.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.8 second │ │ chambers.RA.V0 │ 4.0 milliliter │ │ chambers.RV.EA │ 0.55 millimeter_Hg / milliliter │ │ chambers.RV.EB │ 0.05 millimeter_Hg / milliliter │ │ chambers.RV.TC │ 0.34 second │ │ chambers.RV.TR │ 0.17 second │ │ chambers.RV.tC │ 0.0 second │ │ chambers.RV.V0 │ 10.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.8 millimeter_Hg * second / milliliter │ │ circulation.SYS.C_AR │ 1.2 milliliter / millimeter_Hg │ │ circulation.SYS.R_VEN │ 0.26 millimeter_Hg * second / milliliter │ │ circulation.SYS.C_VEN │ 130.0 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.1625 millimeter_Hg * second / milliliter │ │ circulation.PUL.C_AR │ 10.0 milliliter / millimeter_Hg │ │ circulation.PUL.R_VEN │ 0.1625 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:138 Circulation model initial states (Regazzoni2020) ┏━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━┓ ┃ State ┃ Value ┃ ┡━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━┩ │ V_LA │ 65.0 milliliter │ │ V_LV │ 120.0 milliliter │ │ V_RA │ 65.0 milliliter │ │ V_RV │ 145.0 milliliter │ │ p_AR_SYS │ 80.0 millimeter_Hg │ │ p_VEN_SYS │ 30.0 millimeter_Hg │ │ p_AR_PUL │ 35.0 millimeter_Hg │ │ p_VEN_PUL │ 24.0 millimeter_Hg │ │ Q_AR_SYS │ 0.0 milliliter / second │ │ Q_VEN_SYS │ 0.0 milliliter / second │ │ Q_AR_PUL │ 0.0 milliliter / second │ │ Q_VEN_PUL │ 0.0 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 0x7fb1d9ec1e20>]
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 0x7fb1d9ee31d0>]
t2 = perf_counter()
history = circulation.solve(num_beats=10)
t3 = perf_counter()
circulation.print_info()
INFO INFO:circulation.base: base.py:423 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 ┃ ┡━━━━━━━━╇━━━━━━━━━╇━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━┩ │ 74.087 │ 135.887 │ 71.246 │ 165.457 │ 92.272 │ 3869.306 │ 347.703 │ 369.042 │ 446.676 │ 3961.578 │ 716.746 │ 5125.000 │ └────────┴─────────┴────────┴─────────┴──────────┴───────────┴──────────┴───────────┴─────────┴──────────┴─────────┴──────────┘ Pressures ┏━━━━━━━━┳━━━━━━━━┳━━━━━━━┳━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┓ ┃ p_LA ┃ p_LV ┃ p_RA ┃ p_RV ┃ p_AR_SYS ┃ p_VEN_SYS ┃ p_AR_PUL ┃ p_VEN_PUL ┃ ┡━━━━━━━━╇━━━━━━━━╇━━━━━━━╇━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━┩ │ 10.866 │ 10.467 │ 8.460 │ 7.768 │ 76.893 │ 29.764 │ 34.770 │ 23.065 │ └────────┴────────┴───────┴───────┴──────────┴───────────┴──────────┴───────────┘ Flows ┏━━━━━━━━┳━━━━━━━━┳━━━━━━━━┳━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┓ ┃ Q_MV ┃ Q_AV ┃ Q_TV ┃ Q_PV ┃ Q_AR_SYS ┃ Q_VEN_SYS ┃ Q_AR_PUL ┃ Q_VEN_PUL ┃ ┡━━━━━━━━╇━━━━━━━━╇━━━━━━━━╇━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━┩ │ 51.125 │ -0.001 │ 90.093 │ -0.000 │ 59.299 │ 81.867 │ 72.166 │ 74.828 │ └────────┴────────┴────────┴────────┴──────────┴───────────┴──────────┴───────────┘
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 0x7fb1e97446e0>
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: 3.4581446219999634
Circulation solve time: 2.7490772400000196