LV ellipsoid coupled to a 0D circulatory model and a 0D cell model#
In this example we will couple a 3D LV ellipsoid model to a 0D circulatory model and a 0D cell model.
The 0D cell model#
The 0D cell model we will use is the one TorORd model [TBOP+19] which is an electrophysiology model developed for human ventricular myocytes. The model is available at jtmff/torord in CellML format (we use the endo cell version of the model). To convert the model to Python code we use the gotranx package. See https://finsberg.github.io/gotranx/docs/cli.html for an example of how to convert a CellML model to Python code. We have also coupled this model to the Land model [LPHS+17] which is an excitation-contraction model for cardiac cells. In order to use this model in the 3D simulations we will first run the model to steady state, and then we will use the active tension (Ta), which is an intermediate variable in the model, as the active stress in the 3D model. We will also scale the single cell active stress by a factor of 5.0 to get a more pronounced effect in the 3D model.
The 0D circulatory model#
The 0D circulatory model we will use is the one from [RSA+22] which is a lumped parameter model of the circulatory system. We use an existing implementation in the circulation package.
Coupling the 0D cell model to the 3D model#
TBW
Now we will start by importing the necessary modules
from pathlib import Path
from mpi4py import MPI
import dolfinx
import logging
import circulation
import os
import math
from functools import lru_cache
from dolfinx import log
from matplotlib.gridspec import GridSpec
import matplotlib.pyplot as plt
import numpy as np
import gotranx
import adios4dolfinx
import fenicsx_pulse
import cardiac_geometries
import cardiac_geometries.geometry
Next we set up the logging and the MPI communicator
circulation.log.setup_logging(logging.INFO)
comm = MPI.COMM_WORLD
Next we create the geometry of the LV ellipsoid
geodir = Path("lv_ellipsoid-time-dependent")
if not geodir.exists():
comm.barrier()
cardiac_geometries.mesh.lv_ellipsoid(
outdir=geodir,
create_fibers=True,
fiber_space="Quadrature_6",
r_short_endo=0.025,
r_short_epi=0.035,
r_long_endo=0.09,
r_long_epi=0.097,
psize_ref=0.03,
mu_apex_endo=-math.pi,
mu_base_endo=-math.acos(5 / 17),
mu_apex_epi=-math.pi,
mu_base_epi=-math.acos(5 / 20),
comm=comm,
fiber_angle_epi=-60,
fiber_angle_endo=60,
)
Info : Reading 'lv_ellipsoid-time-dependent/lv_ellipsoid.msh'...
Info : 54 entities
Info : 191 nodes
Info : 990 elements
Info : Done reading 'lv_ellipsoid-time-dependent/lv_ellipsoid.msh'
[03/05/25 19:49:59] INFO INFO:cardiac_geometries.geometry:Reading geometry from lv_ellipsoid-time-dependent geometry.py:145
If the folder already exist, then we just load the geometry
geo = cardiac_geometries.geometry.Geometry.from_folder(
comm=comm,
folder=geodir,
)
INFO INFO:cardiac_geometries.geometry:Reading geometry from lv_ellipsoid-time-dependent geometry.py:145
Next we transform the geometry to a HeartGeometry object
geometry = fenicsx_pulse.HeartGeometry.from_cardiac_geometries(geo, metadata={"quadrature_degree": 6})
Next we create the material object, and we will use the transversely isotropic version of the Holzapfel Ogden model
material_params = fenicsx_pulse.HolzapfelOgden.transversely_isotropic_parameters()
material = fenicsx_pulse.HolzapfelOgden(f0=geo.f0, s0=geo.s0, **material_params) # type: ignore
We use an active stress approach with 30% transverse active stress (see fenicsx_pulse.active_stress.transversely_active_stress())
Ta = fenicsx_pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0)), "kPa")
active_model = fenicsx_pulse.ActiveStress(geo.f0, activation=Ta)
a compressible material model
comp_model = fenicsx_pulse.compressibility.Compressible2()
and assembles the CardiacModel
model = fenicsx_pulse.CardiacModel(
material=material,
active=active_model,
compressibility=comp_model,
)
Next we set up the boundary conditions. We use a Robin boundary condition on the epicardium and the base of the LV
alpha_epi = fenicsx_pulse.Variable(
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1e8)), "Pa / m",
)
robin_epi = fenicsx_pulse.RobinBC(value=alpha_epi, marker=geometry.markers["EPI"][0])
alpha_base = fenicsx_pulse.Variable(
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1e5)), "Pa / m",
)
robin_base = fenicsx_pulse.RobinBC(value=alpha_base, marker=geometry.markers["BASE"][0])
For the pressure we use a Lagrange multiplier method to enforce a given volume. The resulting Lagrange multiplier will be the pressure in the cavity.
To do this we create a Cavity object with a given volume, and specify which marker to use for the boundary condition.
initial_volume = geometry.volume("ENDO")
Volume = dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(initial_volume))
cavity = fenicsx_pulse.problem.Cavity(marker="ENDO", volume=Volume)
We also specify the parameters for the problem and say that we want the base to move freely and that the units of the mesh is meters
parameters = {"base_bc": fenicsx_pulse.problem.BaseBC.free, "mesh_unit": "m"}
Next we set up the problem. We can choose between a static and a dynamic problem by setting the static variable to True or False. Currently the dynamic problem is not working (when coupled to a circulation model), so we will use the static problem for now.
static = True
if static:
outdir = Path("lv_ellipsoid_time_dependent_circulation_static")
bcs = fenicsx_pulse.BoundaryConditions(robin=(robin_epi, robin_base))
problem = fenicsx_pulse.problem.StaticProblem(model=model, geometry=geometry, bcs=bcs, cavities=[cavity], parameters=parameters)
else:
outdir = Path("lv_ellipsoid_time_dependent_circulation_dynamic")
beta_epi = fenicsx_pulse.Variable(
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(5e3)), "Pa s/ m",
)
robin_epi_v = fenicsx_pulse.RobinBC(value=beta_epi, marker=geometry.markers["EPI"][0], damping=True)
beta_base = fenicsx_pulse.Variable(
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(5e3)), "Pa s/ m",
)
robin_base_v = fenicsx_pulse.RobinBC(value=beta_base, marker=geometry.markers["BASE"][0], damping=True)
bcs = fenicsx_pulse.BoundaryConditions(robin=(robin_epi, robin_epi_v, robin_base, robin_base_v))
problem = fenicsx_pulse.problem.DynamicProblem(model=model, geometry=geometry, bcs=bcs, cavities=[cavity], parameters=parameters)
outdir.mkdir(exist_ok=True)
Now we can solve the problem
log.set_log_level(log.LogLevel.INFO)
problem.solve()
[03/05/25 19:50:10] INFO INFO:scifem.solvers:Newton iteration 1: r (abs) = 17.22683609893559 (tol=1e-06), r (rel) = 1.0 (tol=1e-10) solvers.py:212
INFO INFO:scifem.solvers:Newton iteration 2: r (abs) = 1.9300070449246391 (tol=1e-06), r (rel) = 0.1120349107543834 (tol=1e-10) solvers.py:212
INFO INFO:scifem.solvers:Newton iteration 3: r (abs) = 0.005944450390252497 (tol=1e-06), r (rel) = 0.00034506919065769676 (tol=1e-10) solvers.py:212
INFO INFO:scifem.solvers:Newton iteration 4: r (abs) = 1.0269109795830773e-06 (tol=1e-06), r (rel) = 5.961111916810586e-08 (tol=1e-10) solvers.py:212
INFO INFO:scifem.solvers:Newton iteration 5: r (abs) = 2.3268565718629844e-13 (tol=1e-06), r (rel) = 1.3507161492102174e-14 (tol=1e-10) solvers.py:212
5
We also use the time step from the problem to set the time step for the 0D cell model
if static:
dt = 0.001
else:
dt = problem.parameters["dt"].to_base_units()
times = np.arange(0.0, 1.0, dt)
Next we will load the 0D cell model and run it to steady state. Here we use gotranx to load the ode, remove potential singularities and convert it to Python code. We then run the model for 200 beats and save the state of the model at the end of the simulation. We also plot the results. Note also that we extract the index of the active tension (Ta) which will be used to drive the 3D model.
ode = gotranx.load_ode("TorOrdLand.ode")
ode = ode.remove_singularities()
code = gotranx.cli.gotran2py.get_code(
ode, scheme=[gotranx.schemes.Scheme.generalized_rush_larsen], shape=gotranx.codegen.base.Shape.single,
)
Path("TorOrdLand.py").write_text(code)
import TorOrdLand
2025-03-05 19:50:10 [info ] Load ode TorOrdLand.ode
2025-03-05 19:50:11 [info ] Num states 52
2025-03-05 19:50:11 [info ] Num parameters 140
2025-03-05 19:50:14 [warning ] Cannot apply black, please install 'black'
TorOrdLand_model = TorOrdLand.__dict__
Ta_index = TorOrdLand_model["monitor_index"]("Ta")
y = TorOrdLand_model["init_state_values"]()
# Get initial parameter values
p = TorOrdLand_model["init_parameter_values"]()
import numba
fgr = numba.njit(TorOrdLand_model["generalized_rush_larsen"])
mon = numba.njit(TorOrdLand_model["monitor_values"])
V_index = TorOrdLand_model["state_index"]("v")
Ca_index = TorOrdLand_model["state_index"]("cai")
In the cell model the unit of time is in millisenconds, and we set the time step to 0.1 ms
# Time in milliseconds
dt_cell = 0.1
Now we solve the cell model for 200 beats and save the state at the end of the simulation
state_file = outdir / "state.npy"
if not state_file.is_file():
@numba.jit(nopython=True)
def solve_beat(times, states, dt, p, V_index, Ca_index, Vs, Cais, Tas):
for i, ti in enumerate(times):
states[:] = fgr(states, ti, dt, p)
Vs[i] = states[V_index]
Cais[i] = states[Ca_index]
monitor = mon(ti, states, p)
Tas[i] = monitor[Ta_index]
# Time in milliseconds
nbeats = 200
T = 1000.00
times = np.arange(0, T, dt_cell)
all_times = np.arange(0, T * nbeats, dt_cell)
Vs = np.zeros(len(times) * nbeats)
Cais = np.zeros(len(times) * nbeats)
Tas = np.zeros(len(times) * nbeats)
for beat in range(nbeats):
print(f"Solving beat {beat}")
V_tmp = Vs[beat * len(times) : (beat + 1) * len(times)]
Cai_tmp = Cais[beat * len(times) : (beat + 1) * len(times)]
Ta_tmp = Tas[beat * len(times) : (beat + 1) * len(times)]
solve_beat(times, y, dt_cell, p, V_index, Ca_index, V_tmp, Cai_tmp, Ta_tmp)
fig, ax = plt.subplots(3, 2, sharex="col", sharey="row", figsize=(10, 10))
ax[0, 0].plot(all_times, Vs)
ax[1, 0].plot(all_times, Cais)
ax[2, 0].plot(all_times, Tas)
ax[0, 1].plot(times, Vs[-len(times):])
ax[1, 1].plot(times, Cais[-len(times):])
ax[2, 1].plot(times, Tas[-len(times):])
ax[0, 0].set_ylabel("V")
ax[1, 0].set_ylabel("Cai")
ax[2, 0].set_ylabel("Ta")
ax[2, 0].set_xlabel("Time [ms]")
ax[2, 1].set_xlabel("Time [ms]")
fig.savefig(outdir / "Ta_ORdLand.png")
np.save(state_file, y)
Solving beat 0
Solving beat 1
Solving beat 2
Solving beat 3
Solving beat 4
Solving beat 5
Solving beat 6
Solving beat 7
Solving beat 8
Solving beat 9
Solving beat 10
Solving beat 11
Solving beat 12
Solving beat 13
Solving beat 14
Solving beat 15
Solving beat 16
Solving beat 17
Solving beat 18
Solving beat 19
Solving beat 20
Solving beat 21
Solving beat 22
Solving beat 23
Solving beat 24
Solving beat 25
Solving beat 26
Solving beat 27
Solving beat 28
Solving beat 29
Solving beat 30
Solving beat 31
Solving beat 32
Solving beat 33
Solving beat 34
Solving beat 35
Solving beat 36
Solving beat 37
Solving beat 38
Solving beat 39
Solving beat 40
Solving beat 41
Solving beat 42
Solving beat 43
Solving beat 44
Solving beat 45
Solving beat 46
Solving beat 47
Solving beat 48
Solving beat 49
Solving beat 50
Solving beat 51
Solving beat 52
Solving beat 53
Solving beat 54
Solving beat 55
Solving beat 56
Solving beat 57
Solving beat 58
Solving beat 59
Solving beat 60
Solving beat 61
Solving beat 62
Solving beat 63
Solving beat 64
Solving beat 65
Solving beat 66
Solving beat 67
Solving beat 68
Solving beat 69
Solving beat 70
Solving beat 71
Solving beat 72
Solving beat 73
Solving beat 74
Solving beat 75
Solving beat 76
Solving beat 77
Solving beat 78
Solving beat 79
Solving beat 80
Solving beat 81
Solving beat 82
Solving beat 83
Solving beat 84
Solving beat 85
Solving beat 86
Solving beat 87
Solving beat 88
Solving beat 89
Solving beat 90
Solving beat 91
Solving beat 92
Solving beat 93
Solving beat 94
Solving beat 95
Solving beat 96
Solving beat 97
Solving beat 98
Solving beat 99
Solving beat 100
Solving beat 101
Solving beat 102
Solving beat 103
Solving beat 104
Solving beat 105
Solving beat 106
Solving beat 107
Solving beat 108
Solving beat 109
Solving beat 110
Solving beat 111
Solving beat 112
Solving beat 113
Solving beat 114
Solving beat 115
Solving beat 116
Solving beat 117
Solving beat 118
Solving beat 119
Solving beat 120
Solving beat 121
Solving beat 122
Solving beat 123
Solving beat 124
Solving beat 125
Solving beat 126
Solving beat 127
Solving beat 128
Solving beat 129
Solving beat 130
Solving beat 131
Solving beat 132
Solving beat 133
Solving beat 134
Solving beat 135
Solving beat 136
Solving beat 137
Solving beat 138
Solving beat 139
Solving beat 140
Solving beat 141
Solving beat 142
Solving beat 143
Solving beat 144
Solving beat 145
Solving beat 146
Solving beat 147
Solving beat 148
Solving beat 149
Solving beat 150
Solving beat 151
Solving beat 152
Solving beat 153
Solving beat 154
Solving beat 155
Solving beat 156
Solving beat 157
Solving beat 158
Solving beat 159
Solving beat 160
Solving beat 161
Solving beat 162
Solving beat 163
Solving beat 164
Solving beat 165
Solving beat 166
Solving beat 167
Solving beat 168
Solving beat 169
Solving beat 170
Solving beat 171
Solving beat 172
Solving beat 173
Solving beat 174
Solving beat 175
Solving beat 176
Solving beat 177
Solving beat 178
Solving beat 179
Solving beat 180
Solving beat 181
Solving beat 182
Solving beat 183
Solving beat 184
Solving beat 185
Solving beat 186
Solving beat 187
Solving beat 188
Solving beat 189
Solving beat 190
Solving beat 191
Solving beat 192
Solving beat 193
Solving beat 194
Solving beat 195
Solving beat 196
Solving beat 197
Solving beat 198
Solving beat 199
Next we load the state of the cell model and create an ODEState object. This is a simple wrapper around the cell model that allows us to step the model forward in time and get the active tension at a given time.
y = np.load(state_file)
class ODEState:
def __init__(self, y, dt_cell, p, t=0.0):
self.y = y
self.dt_cell = dt_cell
self.p = p
self.t = t
def forward(self, t):
for t_cell in np.arange(self.t, t, self.dt_cell):
self.y[:] = fgr(self.y, t_cell, self.dt_cell, self.p)
self.t = t
return self.y[:]
def Ta(self, t):
monitor = mon(t, self.y, p)
return monitor[Ta_index]
ode_state = ODEState(y, dt_cell, p)
@lru_cache
def get_activation(t: float):
# Find index modulo 1000
t_cell_next = t * 1000
ode_state.forward(t_cell_next)
return ode_state.Ta(t_cell_next) * 5.0
Next we will save the displacement of the LV to a file
vtx = dolfinx.io.VTXWriter(geometry.mesh.comm, f"{outdir}/displacement.bp", [problem.u], engine="BP4")
vtx.write(0.0)
Here we also use adios4dolfinx to save the displacement over at different time steps. Currently it is not a straight forward way to save functions and to later load them in a time dependent simulation in FEniCSx. However adios4dolfinx allows us to save the function to a file and later load it in a time dependent simulation. We will first need to save the mesh to the same file.
filename = Path("function_checkpoint.bp")
adios4dolfinx.write_mesh(filename, geometry.mesh)
Next we set up the callback function that will be called at each time step. Here we save the displacement of the LV, the pressure volume loop, and the active tension, and we also plot the pressure volume loop at each time step.
def callback(model, t: float, save=True):
model.results["Ta"].append(get_activation(t))
if save:
adios4dolfinx.write_function(filename, problem.u, time=t, name="displacement")
vtx.write(t)
if comm.rank == 0:
fig = plt.figure(layout="constrained")
gs = GridSpec(3, 2, figure=fig)
ax1 = fig.add_subplot(gs[:, 0])
ax2 = fig.add_subplot(gs[0, 1])
ax3 = fig.add_subplot(gs[1, 1])
ax4 = fig.add_subplot(gs[2, 1])
ax1.plot(model.results["V_LV"], model.results["p_LV"])
ax1.set_xlabel("V [mL]")
ax1.set_ylabel("p [mmHg]")
ax2.plot(model.results["time"], model.results["p_LV"])
ax2.set_ylabel("p [mmHg]")
ax3.plot(model.results["time"], model.results["V_LV"])
ax3.set_ylabel("V [mL]")
ax4.plot(model.results["time"], model.results["Ta"])
ax4.set_ylabel("Ta [kPa]")
for axi in [ax2, ax3, ax4]:
axi.set_xlabel("Time [s]")
fig.savefig(outdir / "pv_loop_incremental.png")
plt.close(fig)
fig, ax = plt.subplots(4, 1)
ax[0].plot(model.results["V_LV"], model.results["p_LV"])
ax[0].set_xlabel("V [mL]")
ax[0].set_ylabel("p [mmHg]")
ax[1].plot(model.results["time"], model.results["p_LV"])
ax[2].plot(model.results["time"], model.results["V_LV"])
ax[3].plot(model.results["time"], model.results["Ta"])
fig.savefig(outdir / "pv_loop_incremental.png")
plt.close(fig)
Now we will set ut the function to calculate the pressure in the LV. In this function we will receive a volume and a time from the 0D circulation model, we will set the active tension from the 0D cell model, and solve the 3D model to get the pressure in the LV.
def p_LV_func(V_LV, t):
print("Calculating pressure at time", t)
value = get_activation(t)
print("Time", t, "Activation", value)
Ta.assign(value)
Volume.value = V_LV * 1e-6
problem.solve()
pendo = problem.cavity_pressures[0]
pendo_kPa = pendo.x.array[0] * 1e-3
return circulation.units.kPa_to_mmHg(pendo_kPa)
Finally we create the circulation model and and pass in al the arguements. We also set the initial volume of the LV to be the same as the one computed from the 3D model.
mL = circulation.units.ureg("mL")
add_units = False
surface_area = geometry.surface_area("ENDO")
initial_volume = geo.mesh.comm.allreduce(geometry.volume("ENDO", u=problem.u), op=MPI.SUM) * 1e6
print(f"Initial volume: {initial_volume}")
init_state = {"V_LV": initial_volume * mL}
[2025-03-05 19:51:20.258] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-03-05 19:51:20.258] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-03-05 19:51:20.258] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-03-05 19:51:20.258] [info] Requesting connectivity (3, 0) - (2, 0)
Initial volume: 135.52291214199383
[2025-03-05 19:51:20.261] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-03-05 19:51:20.261] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-03-05 19:51:20.261] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-03-05 19:51:20.261] [info] Requesting connectivity (3, 0) - (2, 0)
circulation_model_3D = circulation.regazzoni2020.Regazzoni2020(
add_units=add_units,
callback=callback,
p_LV_func=p_LV_func,
verbose=True,
comm=comm,
outdir=outdir,
initial_state=init_state,
)
# Set end time for early stopping if running in CI
end_time = 2 * dt if os.getenv("CI") else None
circulation_model_3D.solve(num_cycles=5, initial_state=init_state, dt=dt, T=end_time)
circulation_model_3D.print_info()
[03/05/25 19:51:20] INFO INFO:circulation.base: base.py:131 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:137 Circulation model initial states (Regazzoni2020) ┏━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ┃ State ┃ Value ┃ ┡━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ │ V_LA │ 65.0 milliliter │ │ V_LV │ 135.52291214199383 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 │ └───────────┴───────────────────────────────┘
Calculating pressure at time 0.0
Time 0.0 Activation 0.2849891894884486
INFO INFO:scifem.solvers:Newton iteration 1: r (abs) = 54.85743991019322 (tol=1e-06), r (rel) = 1.0 (tol=1e-10) solvers.py:212
INFO INFO:scifem.solvers:Newton iteration 2: r (abs) = 0.09235387373496964 (tol=1e-06), r (rel) = 0.001683525040289186 (tol=1e-10) solvers.py:212
INFO INFO:scifem.solvers:Newton iteration 3: r (abs) = 0.004979426070979202 (tol=1e-06), r (rel) = 9.077029622838743e-05 (tol=1e-10) solvers.py:212
INFO INFO:scifem.solvers:Newton iteration 4: r (abs) = 1.082025064853668e-06 (tol=1e-06), r (rel) = 1.972430843701501e-08 (tol=1e-10) solvers.py:212
INFO INFO:scifem.solvers:Newton iteration 5: r (abs) = 1.3389410871828277e-12 (tol=1e-06), r (rel) = 2.440764806696776e-14 (tol=1e-10) solvers.py:212
Calculating pressure at time 0.0
Time 0.0 Activation 0.2849891894884486
[03/05/25 19:51:21] INFO INFO:scifem.solvers:Newton iteration 1: r (abs) = 2.3470475880201367e-13 (tol=1e-06), r (rel) = 2.3470475880201368e-07 (tol=1e-10) solvers.py:212
INFO INFO:circulation.base: base.py:321 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 ┃ ┡━━━━━━━━╇━━━━━━━━━╇━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━┩ │ 63.779 │ 136.744 │ 64.881 │ 145.119 │ 96.000 │ 3900.000 │ 350.000 │ 384.000 │ 410.523 │ 3996.000 │ 734.000 │ 5140.523 │ └────────┴─────────┴────────┴─────────┴──────────┴───────────┴──────────┴───────────┴─────────┴──────────┴─────────┴──────────┘ Pressures ┏━━━━━━━┳━━━━━━━┳━━━━━━━┳━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┓ ┃ p_LA ┃ p_LV ┃ p_RA ┃ p_RV ┃ p_AR_SYS ┃ p_VEN_SYS ┃ p_AR_PUL ┃ p_VEN_PUL ┃ ┡━━━━━━━╇━━━━━━━╇━━━━━━━╇━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━┩ │ 9.440 │ 0.267 │ 7.656 │ 6.750 │ 80.000 │ 30.000 │ 35.000 │ 24.000 │ └───────┴───────┴───────┴───────┴──────────┴───────────┴──────────┴───────────┘ Flows ┏━━━━━━━━━━┳━━━━━━━━┳━━━━━━━━━┳━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┓ ┃ Q_MV ┃ Q_AV ┃ Q_TV ┃ Q_PV ┃ Q_AR_SYS ┃ Q_VEN_SYS ┃ Q_AR_PUL ┃ Q_VEN_PUL ┃ ┡━━━━━━━━━━╇━━━━━━━━╇━━━━━━━━━╇━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━┩ │ 1220.919 │ -0.001 │ 118.628 │ -0.000 │ 10.000 │ 44.688 │ 22.000 │ 29.120 │ └──────────┴────────┴─────────┴────────┴──────────┴───────────┴──────────┴───────────┘
Calculating pressure at time 0.001
Time 0.001 Activation 0.28501827278006825
INFO INFO:scifem.solvers:Newton iteration 1: r (abs) = 93.71486838103476 (tol=1e-06), r (rel) = 1.0 (tol=1e-10) solvers.py:212
INFO INFO:scifem.solvers:Newton iteration 2: r (abs) = 2.319147349366899 (tol=1e-06), r (rel) = 0.024746845291800344 (tol=1e-10) solvers.py:212
INFO INFO:scifem.solvers:Newton iteration 3: r (abs) = 0.06051001410276901 (tol=1e-06), r (rel) = 0.0006456821115806478 (tol=1e-10) solvers.py:212
[03/05/25 19:51:22] INFO INFO:scifem.solvers:Newton iteration 4: r (abs) = 9.506435797538463e-06 (tol=1e-06), r (rel) = 1.0143999518717028e-07 (tol=1e-10) solvers.py:212
INFO INFO:scifem.solvers:Newton iteration 5: r (abs) = 6.127782910436119e-11 (tol=1e-06), r (rel) = 6.538752085230703e-13 (tol=1e-10) solvers.py:212
INFO INFO:circulation.base: base.py:321 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 ┃ ┡━━━━━━━━╇━━━━━━━━━╇━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━┩ │ 62.711 │ 137.841 │ 64.813 │ 145.232 │ 95.990 │ 3899.965 │ 349.978 │ 383.993 │ 410.597 │ 3995.955 │ 733.971 │ 5140.523 │ └────────┴─────────┴────────┴─────────┴──────────┴───────────┴──────────┴───────────┴─────────┴──────────┴─────────┴──────────┘ Pressures ┏━━━━━━━┳━━━━━━━┳━━━━━━━┳━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┓ ┃ p_LA ┃ p_LV ┃ p_RA ┃ p_RV ┃ p_AR_SYS ┃ p_VEN_SYS ┃ p_AR_PUL ┃ p_VEN_PUL ┃ ┡━━━━━━━╇━━━━━━━╇━━━━━━━╇━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━┩ │ 9.231 │ 0.988 │ 7.623 │ 6.758 │ 79.992 │ 30.000 │ 34.998 │ 24.000 │ └───────┴───────┴───────┴───────┴──────────┴───────────┴──────────┴───────────┘ Flows ┏━━━━━━━━━━┳━━━━━━━━┳━━━━━━━━━┳━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┓ ┃ Q_MV ┃ Q_AV ┃ Q_TV ┃ Q_PV ┃ Q_AR_SYS ┃ Q_VEN_SYS ┃ Q_AR_PUL ┃ Q_VEN_PUL ┃ ┡━━━━━━━━━━╇━━━━━━━━╇━━━━━━━━━╇━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━┩ │ 1096.870 │ -0.001 │ 113.237 │ -0.000 │ 18.400 │ 66.204 │ 36.850 │ 49.194 │ └──────────┴────────┴─────────┴────────┴──────────┴───────────┴──────────┴───────────┘
INFO INFO:circulation.base: base.py:321 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 ┃ ┡━━━━━━━━╇━━━━━━━━━╇━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━┩ │ 62.711 │ 137.841 │ 64.813 │ 145.232 │ 95.990 │ 3899.965 │ 349.978 │ 383.993 │ 410.597 │ 3995.955 │ 733.971 │ 5140.523 │ └────────┴─────────┴────────┴─────────┴──────────┴───────────┴──────────┴───────────┴─────────┴──────────┴─────────┴──────────┘ Pressures ┏━━━━━━━┳━━━━━━━┳━━━━━━━┳━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┓ ┃ p_LA ┃ p_LV ┃ p_RA ┃ p_RV ┃ p_AR_SYS ┃ p_VEN_SYS ┃ p_AR_PUL ┃ p_VEN_PUL ┃ ┡━━━━━━━╇━━━━━━━╇━━━━━━━╇━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━┩ │ 9.231 │ 0.988 │ 7.623 │ 6.758 │ 79.992 │ 30.000 │ 34.998 │ 24.000 │ └───────┴───────┴───────┴───────┴──────────┴───────────┴──────────┴───────────┘ Flows ┏━━━━━━━━━━┳━━━━━━━━┳━━━━━━━━━┳━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┓ ┃ Q_MV ┃ Q_AV ┃ Q_TV ┃ Q_PV ┃ Q_AR_SYS ┃ Q_VEN_SYS ┃ Q_AR_PUL ┃ Q_VEN_PUL ┃ ┡━━━━━━━━━━╇━━━━━━━━╇━━━━━━━━━╇━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━┩ │ 1096.870 │ -0.001 │ 113.237 │ -0.000 │ 18.400 │ 66.204 │ 36.850 │ 49.194 │ └──────────┴────────┴─────────┴────────┴──────────┴───────────┴──────────┴───────────┘
Fig. 1 Pressure volume loop for the LV.#
References#
[LPHS+17]
Sander Land, So-Jin Park-Holohan, Nicolas P Smith, Cristobal G Dos Remedios, Jonathan C Kentish, and Steven A Niederer. A model of cardiac contraction based on novel measurements of tension development in human cardiomyocytes. Journal of molecular and cellular cardiology, 106:68–83, 2017.
[RSA+22]
Francesco Regazzoni, Matteo Salvador, Pasquale Claudio Africa, Marco Fedele, Luca Dedè, and Alfio Quarteroni. A cardiac electromechanical model coupled with a lumped-parameter model for closed-loop blood circulation. Journal of Computational Physics, 457:111083, 2022.
[TBOP+19]
Jakub Tomek, Alfonso Bueno-Orovio, Elisa Passini, Xin Zhou, Ana Minchole, Oliver Britton, Chiara Bartolucci, Stefano Severi, Alvin Shrier, Laszlo Virag, and others. Development, calibration, and validation of a novel human ventricular myocyte model in health, disease, and drug block. Elife, 8:e48890, 2019.