A complete cardiac cycle on a prestressed biventricular mesh#
In this demo we take both ventricles of a UK Biobank atlas mesh through
whole heartbeats, using the five-phase, Alya-style cycle of pulse.cycle.
Each ventricle has a single pressure unknown, held by a
pulse.problem.CavityControl. A pulse.cycle.CycleController decides, step
by step, which constraint that unknown satisfies: a prescribed pressure while
the ventricle is loaded up to end diastole, a fixed volume while both valves
are shut, a three-element Windkessel while the ventricle ejects, and a volume
growing at a set rate while it fills. Switching between them
only changes the values of a few constants, so the problem is built once and
never rebuilt.
The pipeline is:
Geometry. We generate the mesh from the atlas, rotate it so that the base normal points along x, and compute fibres with LDRB, much as in the rotated BiV demo.
Prestress. The mesh is imaged at end diastole, so it is already loaded. We recover the unloaded reference configuration by solving the inverse elasticity problem, as in the BiV prestress demo.
PRELOAD. The first phase of the cycle ramps each cavity pressure from zero back up to its end-diastolic value, which inflates the unloaded mesh back to the imaged shape.
Beats. From there the controller takes both ventricles through contraction, ejection, relaxation and filling, and on into the next beat. The active tension follows the Bestel activation model of the
circulationpackage, which we integrate once withscipy, so this demo needs both of them besidespulse(thedemoextra installs them).
A run like this one takes a while. Restarting a simulation shows how to write
a checkpoint of a CycleController run and carry on from it later.
import logging
import os
import shutil
from pathlib import Path
from mpi4py import MPI
# A sibling module in this directory, not a package.
import animation
import dolfinx
import io4dolfinx
import ldrb
import matplotlib.pyplot as plt
import numpy as np
from circulation import bestel
from scipy.integrate import solve_ivp
import cardiac_geometries
import cardiac_geometries.geometry
import pulse
from pulse import cycle
from pulse.circulation import mL, mmHg
# Setup logging to print only from rank 0
class MPIFilter(logging.Filter):
def __init__(self, comm, *args, **kwargs):
super().__init__(*args, **kwargs)
self.comm = comm
def filter(self, record):
return 1 if self.comm.rank == 0 else 0
outdir = Path("results_biv_complete_cycle")
outdir.mkdir(parents=True, exist_ok=True)
geodir = outdir / "geometry"
comm = MPI.COMM_WORLD
logging.basicConfig(level=logging.INFO)
# The filter sits on the handler, so that it also catches the messages of
# `pulse.cycle`, which logs every phase transition.
mpi_filter = MPIFilter(comm)
for handler in logging.getLogger().handlers:
handler.addFilter(mpi_filter)
logger = logging.getLogger("pulse")
logging.getLogger("scifem").setLevel(logging.WARNING)
logging.getLogger("matplotlib").setLevel(logging.WARNING)
We run two beats of 0.8 s with a time step of 1 ms. Under CI, where this
page is built, the run is cut to two steps. Setting PULSE_MAX_STEPS to a
positive number cuts it to that many steps instead, which is a convenient way
to check the first part of a beat without running the whole thing.
_ci = os.getenv("CI", "").strip().lower()
IN_CI = _ci not in ("", "0", "false", "no", "off")
PERIOD = 0.8 # s
NUM_BEATS = 2
DT = 1e-3 # s
max_steps = 2 if IN_CI else int(round(NUM_BEATS * PERIOD / DT))
MAX_STEPS = int(os.getenv("PULSE_MAX_STEPS", "0"))
if MAX_STEPS > 0:
max_steps = MAX_STEPS
Geometry generation and rotation#
We generate the BiV geometry from the UK Biobank Atlas, rotate it to align the base normal with the x-axis, and generate fiber fields using LDRB. The fibers are based on the fiber orientation angles from [DSIB+19]. The fibers used for the mechanics simulation are in a quadrature space to avoid interpolation errors. We also store fibers in a DG 1 space as additional data, which is useful if we want to compute fiber stress or strain at intermediate points later on.
if not (geodir / "geometry.bp").exists():
logger.info("Generating and processing geometry...")
mode = -1
std = 0
char_length = 10.0
geo = cardiac_geometries.mesh.ukb(
outdir=geodir,
comm=comm,
mode=mode,
std=std,
case="ED",
char_length_max=char_length,
char_length_min=char_length,
clipped=True,
)
# Rotate Mesh (Base Normal -> X-axis)
geo = geo.rotate(target_normal=[1.0, 0.0, 0.0], base_marker="BASE")
fiber_angles = dict(
alpha_endo_lv=60,
alpha_epi_lv=-60,
alpha_endo_rv=90,
alpha_epi_rv=-25,
beta_endo_lv=-20,
beta_epi_lv=20,
beta_endo_rv=0,
beta_epi_rv=20,
)
# Generate Fibers (LDRB)
system = ldrb.dolfinx_ldrb(
mesh=geo.mesh,
ffun=geo.ffun,
markers=cardiac_geometries.mesh.transform_markers(geo.markers, clipped=True),
**fiber_angles,
fiber_space="Quadrature_6",
)
# Additional Vectors for Analysis in DG 1 Space for computing stress/strain later
fiber_space = "DG_1"
system_fibers = ldrb.dolfinx_ldrb(
mesh=geo.mesh,
ffun=geo.ffun,
markers=cardiac_geometries.mesh.transform_markers(geo.markers, clipped=True),
**fiber_angles,
fiber_space=fiber_space,
)
# Save Everything
additional_data = {
"f0_DG_1": system_fibers.f0,
"s0_DG_1": system_fibers.s0,
"n0_DG_1": system_fibers.n0,
}
if (geodir / "geometry.bp").exists():
shutil.rmtree(geodir / "geometry.bp")
cardiac_geometries.geometry.save_geometry(
path=geodir / "geometry.bp",
mesh=geo.mesh,
ffun=geo.ffun,
markers=geo.markers,
info=geo.info,
f0=system.f0,
s0=system.s0,
n0=system.n0,
additional_data=additional_data,
)
INFO:pulse:Generating and processing geometry...
INFO:ukb.atlas:Generating points from /github/home/.ukb/UKBRVLV.h5
INFO:ukb.atlas:Using mode -1 and std 0.0
INFO:ukb.surface:Saved results_biv_complete_cycle/geometry/EPI_ED.stl
INFO:ukb.surface:Saved results_biv_complete_cycle/geometry/MV_ED.stl
INFO:ukb.surface:Saved results_biv_complete_cycle/geometry/AV_ED.stl
INFO:ukb.surface:Saved results_biv_complete_cycle/geometry/TV_ED.stl
INFO:ukb.surface:Saved results_biv_complete_cycle/geometry/PV_ED.stl
INFO:ukb.surface:Saved results_biv_complete_cycle/geometry/LV_ED.stl
INFO:ukb.surface:Saved results_biv_complete_cycle/geometry/RV_ED.stl
INFO:ukb.surface:Saved results_biv_complete_cycle/geometry/RVFW_ED.stl
INFO:ukb.clip:Folder: results_biv_complete_cycle/geometry
INFO:ukb.clip:Case: ED
INFO:ukb.clip:Origin: [-13.612554383622273, 18.55767189380559, 15.135103714006394]
INFO:ukb.clip:Normal: [-0.7160843664428893, 0.544394641424108, 0.4368725838557541]
INFO:ukb.clip:Reading results_biv_complete_cycle/geometry/LV_ED.stl
Warning: PLY writer doesn't support multidimensional point data yet. Skipping Normals.
Warning: PLY doesn't support 64-bit integers. Casting down to 32-bit.
INFO:ukb.clip:Saved results_biv_complete_cycle/geometry/lv_clipped.ply
INFO:ukb.clip:Reading results_biv_complete_cycle/geometry/RV_ED.stl
INFO:ukb.clip:Reading results_biv_complete_cycle/geometry/RVFW_ED.stl
INFO:ukb.clip:Merging RV and RVFW
INFO:ukb.clip:Smoothing RV
INFO:ukb.clip:Saving results_biv_complete_cycle/geometry/rv_clipped.ply
Warning: PLY writer doesn't support multidimensional point data yet. Skipping Normals.
Warning: PLY doesn't support 64-bit integers. Casting down to 32-bit.
INFO:ukb.clip:Reading results_biv_complete_cycle/geometry/EPI_ED.stl
INFO:ukb.clip:Saving results_biv_complete_cycle/geometry/epi_clipped.ply
Warning: PLY writer doesn't support multidimensional point data yet. Skipping Normals.
Warning: PLY doesn't support 64-bit integers. Casting down to 32-bit.
INFO:ukb.mesh:Creating clipped mesh for ED with char_length_max=10.0, char_length_min=10.0
0
INFO:ukb.mesh:Created mesh results_biv_complete_cycle/geometry/ED_clipped.msh
2026-10-10 10:40:49 [debug ] Convert file results_biv_complete_cycle/geometry/ED_clipped.msh to dolfin
INFO:ldrb.ldrb:Calculating scalar fields
INFO:ldrb.ldrb:Compute scalar laplacian solutions with the markers:
lv: [1]
rv: [2]
epi: [3]
base: [4]
INFO:ldrb.ldrb: Num vertices: 708
INFO:ldrb.ldrb: Num cells: 2176
INFO:ldrb.ldrb: Apex coord: (47.18, -18.85, -29.35)
INFO:ldrb.ldrb:
Calculating gradients
Info : Reading 'results_biv_complete_cycle/geometry/ED_clipped.msh'...
Info : 11 entities
Info : 708 nodes
Info : 3496 elements
Info : 3 parametrizations
Info : [ 0%] Processing parametrizations
Info : [ 10%] Processing parametrizations
Info : [ 40%] Processing parametrizations
Info : Done reading 'results_biv_complete_cycle/geometry/ED_clipped.msh'
INFO:ldrb.ldrb:Compute fiber-sheet system
INFO:ldrb.ldrb:Angles:
INFO:ldrb.ldrb:alpha:
endo_lv: 60
epi_lv: -60
endo_septum: 60
epi_septum: -60
endo_rv: 60
epi_rv: -60
INFO:ldrb.ldrb:beta:
endo_lv: 0
epi_lv: 0
endo_septum: 0
epi_septum: 0
endo_rv: 0
epi_rv: 0
2026-10-10 10:40:51 [debug ] Write f0: f0
2026-10-10 10:40:51 [debug ] Write s0: s0
2026-10-10 10:40:51 [debug ] Write n0: n0
2026-10-10 10:40:51 [debug ] Write lv: f
2026-10-10 10:40:51 [debug ] Write rv: f
2026-10-10 10:40:51 [debug ] Write epi: f
2026-10-10 10:40:51 [debug ] Write lv_rv: f
2026-10-10 10:40:51 [debug ] Write apex: f
2026-10-10 10:40:51 [debug ] Write lv_scalar: f
2026-10-10 10:40:51 [debug ] Write rv_scalar: f
2026-10-10 10:40:51 [debug ] Write epi_scalar: f
2026-10-10 10:40:51 [debug ] Write lv_rv_scalar: f
2026-10-10 10:40:51 [debug ] Write lv_gradient: f
2026-10-10 10:40:51 [debug ] Write rv_gradient: f
2026-10-10 10:40:51 [debug ] Write epi_gradient: f
2026-10-10 10:40:51 [debug ] Write apex_gradient: f
2026-10-10 10:40:51 [debug ] Write markers_scalar: f
2026-10-10 10:40:51 [info ] Rotated geometry. Base normal [-0.71608437 0.54439464 0.43687258] aligned to [1.0, 0.0, 0.0]
2026-10-10 10:40:51 [debug ] Rotation matrix:
[[-0.71608437 0.54439464 0.43687258]
[-0.54439464 -0.04385068 -0.83768227]
[-0.43687258 -0.83768227 0.32776631]]
INFO:ldrb.ldrb:Calculating scalar fields
INFO:ldrb.ldrb:Compute scalar laplacian solutions with the markers:
lv: [1]
rv: [2]
epi: [3]
base: [4]
INFO:ldrb.ldrb: Num vertices: 708
INFO:ldrb.ldrb: Num cells: 2176
INFO:ldrb.ldrb: Apex coord: (-56.87, -0.27, -14.44)
INFO:ldrb.ldrb:
Calculating gradients
INFO:ldrb.ldrb:Compute fiber-sheet system
INFO:ldrb.ldrb:Angles:
INFO:ldrb.ldrb:alpha:
endo_lv: 60
epi_lv: -60
endo_septum: 60
epi_septum: -60
endo_rv: 90
epi_rv: -25
INFO:ldrb.ldrb:beta:
endo_lv: -20
epi_lv: 20
endo_septum: -20
epi_septum: 20
endo_rv: -20
epi_rv: 20
INFO:ldrb.ldrb:Calculating scalar fields
INFO:ldrb.ldrb:Compute scalar laplacian solutions with the markers:
lv: [1]
rv: [2]
epi: [3]
base: [4]
INFO:ldrb.ldrb: Num vertices: 708
INFO:ldrb.ldrb: Num cells: 2176
INFO:ldrb.ldrb: Apex coord: (-56.87, -0.27, -14.44)
INFO:ldrb.ldrb:
Calculating gradients
INFO:ldrb.ldrb:Compute fiber-sheet system
INFO:ldrb.ldrb:Angles:
INFO:ldrb.ldrb:alpha:
endo_lv: 60
epi_lv: -60
endo_septum: 60
epi_septum: -60
endo_rv: 90
epi_rv: -25
INFO:ldrb.ldrb:beta:
endo_lv: -20
epi_lv: 20
endo_septum: -20
epi_septum: 20
endo_rv: -20
epi_rv: 20
comm.barrier()
We load the generated geometry
geo = cardiac_geometries.geometry.Geometry.from_folder(comm=comm, folder=geodir)
INFO:cardiac_geometries.geometry:Reading geometry from results_biv_complete_cycle/geometry
and scale it from millimetres to metres. CavityControl works in SI units
throughout (cubic metres, pascals), so a controlled cavity needs
mesh_unit == "m".
scale = 1e-3
geo.mesh.geometry.x[:] *= scale
mesh_unit = "m"
geometry = pulse.HeartGeometry.from_cardiac_geometries(
geo, metadata={"quadrature_degree": 6},
)
The sliding-base condition below holds a single displacement component, the x component, on the base, so we check that the base normal really is the x axis. We keep the normal to stand the mesh upright in the video.
up = animation.base_normal(geometry, "BASE")
if abs(up[0]) < 0.99:
raise RuntimeError(
f"the base normal is {up.round(3)}, not the x axis the sliding-base "
"condition assumes -- the rotation above did not take effect",
)
These are the volumes of the imaged, end-diastolic mesh.
lvv_target = comm.allreduce(geometry.volume("LV"), op=MPI.SUM)
rvv_target = comm.allreduce(geometry.volume("RV"), op=MPI.SUM)
logger.info(
f"ED Volumes: LV={lvv_target / mL:.2f} mL, RV={rvv_target / mL:.2f} mL",
)
INFO:pulse:ED Volumes: LV=109.04 mL, RV=74.92 mL
The cycle#
Each ventricle carries a pulse.cycle.CavityCycle, the state machine that
moves it through five phases. Every phase sets one constraint on the
cavity’s pressure unknown, computed from the state at the start of the step:
Phase |
Constraint |
|---|---|
|
pressure: a linear ramp, in two pieces |
|
volume: \(V\) = the end-diastolic volume |
|
pressure: an implicit three-element Windkessel, affine in \(V\) |
|
volume: \(V\) = the end-systolic volume |
|
volume: \(V = V_n + \text{rate}\,\Delta t\) |
In the first beat, PRELOAD’s ramp goes from 0 to preload_pressure at
t_zero. In every beat it then goes on from preload_pressure to
p_end_diastole at t_end_diastole. \(V_n\) is the volume at the start of the
step and \(\Delta t\) its length.
During ejection the outflow is \(Q = -(V - V_n)/\Delta t\), and the Windkessel advances its compliance (arterial) pressure by backward Euler,
with \(C\) the compliance, \(R_p\) the peripheral resistance and \(R_c\) the characteristic impedance. Written out in terms of the still-unknown volume \(V\), this makes the cavity pressure \(P_v\) an affine function of \(V\), which Newton then enforces together with the mechanics.
The controller moves on from a phase when:
PRELOADreachest_end_diastoleof the current beat;in
ISOVOLUMIC_CONTRACTION, the cavity pressure exceeds the compliance pressure \(P_c\), i.e. the outflow valve opens;in
EJECTION, more thanmin_ejection_duration(10 ms) after the valve opened, either the outflow stops or the cavity pressure drops below \(P_c\);in
ISOVOLUMIC_RELAXATION, the cavity pressure falls belowp_fill;FILLINGreachest_zeroof the next beat, which starts the nextPRELOAD.
Outside ejection the valve is shut, and once a cavity has ejected, its \(P_c\) drains through \(R_p\).
For the left ventricle we use the timings and the Windkessel that pulse’s
own tests run the cycle with. Note that everything in pulse.cycle is SI. We
write the circuit-side values in millilitres and mmHg, and mL and mmHg
convert them.
The filling rate is our own. Filling at a prescribed rate knows nothing of
the next beat, whose PRELOAD starts again from preload_pressure. A faster
rate fills the ventricle past the volume it has at that pressure, and the
pressure then drops as PRELOAD takes over; a slower one leaves it short, and
the pressure jumps. We set each ventricle’s rate so that its pressure at the
end of filling is close to preload_pressure: when the second beat’s
PRELOAD starts, the left ventricle holds 96.1 mL at 3.41 mmHg
(preload_pressure is 500 Pa, 3.75 mmHg), and the right 67.1 mL at
1.53 mmHg (200 Pa, 1.50 mmHg). A constant rate is slower than the wall
recoils early in filling, so the cavity pressure first dips below zero, to
about -2.5 mmHg in the left ventricle and -2.3 mmHg in the right.
lv_params = cycle.CycleParams(
t_zero=0.05,
preload_pressure=500.0,
t_end_diastole=0.12,
p_end_diastole=1000.0,
p_fill=500.0,
period=PERIOD,
windkessel=cycle.Windkessel(
p_init=9000.0,
compliance=1.5 * mL / mmHg,
resistance=1.1 * mmHg / mL,
characteristic_impedance=0.03 * mmHg / mL,
),
filling=cycle.PrescribedInflow(rate=0.104 * mL / 1e-3),
)
The right ventricle pumps into the pulmonary circulation, which has a lower
resistance (0.5 against 1.1 mmHg s/mL) and a higher compliance than the
systemic one, and it fills at a lower pressure. We have tuned these values by
hand on this mesh; they are not taken from a reference. In particular, \(R_p\)
was raised from 0.15 to 0.5 mmHg s/mL to give a right ventricular systolic
pressure of about 24 mmHg. p_init is tuned so that the first beat
opens the pulmonary valve at about the pressure the compliance has drained to
by the second beat (2.50 against 2.33 kPa), so the two beats peak at 24.7 and
23.6 mmHg.
rv_params = cycle.CycleParams(
t_zero=0.05,
preload_pressure=200.0,
t_end_diastole=0.12,
p_end_diastole=400.0,
p_fill=200.0,
period=PERIOD,
windkessel=cycle.Windkessel(
p_init=2500.0,
compliance=4.0 * mL / mmHg,
resistance=0.5 * mmHg / mL,
characteristic_impedance=0.01 * mmHg / mL,
),
filling=cycle.PrescribedInflow(rate=0.057 * mL / 1e-3),
)
params = {"LV": lv_params, "RV": rv_params}
Activation#
The active tension follows the activation model of [BClementS01], the same one that drives LV ellipsoid with time dependent pressure and activation - Benchmark and Monolithic 3D-0D coupling: an LV in a closed circulation. The tension \(\tau\) obeys
where the rate \(a(t)\) switches smoothly from \(a_{\min} = -30\,\text{s}^{-1}\)
to \(a_{\max} = 5\,\text{s}^{-1}\) at t_sys and back at t_dias. In
between, \(\tau\) rises towards \(\sigma_0\) with a time constant of
\(1/a_{\max}\) = 200 ms, and afterwards it decays with one of
\(1/|a_{\min}|\) ≈ 33 ms.
We put t_sys at the end of PRELOAD, so that the tension starts to build
as soon as the ventricles stop being loaded, and t_dias 280 ms later. We
integrate one beat of the model once, here, scale \(\tau\) to a unit peak and
repeat it every PERIOD, as \(T_a(t) = T_{\max}\,\tilde{\tau}(t \bmod
\text{PERIOD})\), with \(\tilde{\tau} = \tau / \max\tau\), the same in every
element of both ventricles. The tension is above a tenth of its peak from
0.143 s to 0.477 s of each beat and peaks at 0.395 s. It has died out long before the next beat starts. We tuned
T_MAX by hand, for a left ventricular peak pressure of about 100 mmHg at an
ejection fraction of about 50 %. As the tension builds slowly, the left
ventricle spends 69 ms of the second beat in isovolumic contraction before
its pressure exceeds the compliance pressure \(P_c\), and it then ejects for
210 ms, until just after the tension peaks.
T_MAX = 100.0 # kPa, uniform over both ventricles
activation_model = bestel.BestelActivation(
parameters={
"t_sys": lv_params.t_end_diastole,
"t_dias": lv_params.t_end_diastole + 0.28,
},
)
beat_times = np.arange(0.0, PERIOD + DT / 2, DT)
# `max_step` keeps the solver from stepping over the switch at `t_sys`, where
# the right-hand side is still zero.
activation_shape = solve_ivp(
activation_model,
[0.0, PERIOD],
[0.0],
t_eval=beat_times,
method="Radau",
max_step=DT,
).y[0]
activation_shape /= activation_shape.max()
above = beat_times[activation_shape > 0.1]
logger.info(
f"Activation peaks at {beat_times[activation_shape.argmax()]:.3f} s and is "
f"above 10% of its peak from {above[0]:.3f} s to {above[-1]:.3f} s",
)
INFO:circulation.bestel:
Bestel activation model
parameters
┏━━━━━━━━━━━┳━━━━━━━━━━┓
┃ Parameter ┃ Value ┃
┡━━━━━━━━━━━╇━━━━━━━━━━┩
│ t_sys │ 0.12 │
│ t_dias │ 0.4 │
│ gamma │ 0.005 │
│ a_max │ 5.0 │
│ a_min │ -30.0 │
│ sigma_0 │ 150000.0 │
└───────────┴──────────┘
INFO:pulse:Activation peaks at 0.395 s and is above 10% of its peak from 0.143 s to 0.477 s
def activation(t: float) -> float:
"""Active tension in kPa at time t, repeating every beat."""
return T_MAX * float(np.interp(t % PERIOD, beat_times, activation_shape))
Crossbridges at every quadrature point in a closed-loop circulation replaces this prescribed shape with a crossbridge model, and a single strength with a different one per region.
The model#
The material is the transversely isotropic Holzapfel-Ogden model, with an active stress along the fibres. The epicardium and the base rest on springs, and the base may slide in its own plane but not leave it.
We will solve the dynamic problem, and we damp it with a viscous term.
Undamped, the wall rings after the switch into the volume constraint of
isovolumic relaxation: the cavity pressure saw-tooths from step to step, and
Newton needs many iterations, or diverges. With Viscous, the pressure falls
smoothly through relaxation. The viscous term acts on the strain rate, which
only the dynamic problem has, so it plays no part in the static prestressing
solve below.
The viscous term resists how fast the wall deforms, but not how fast the
ventricles move as a whole on the elastic springs of the epicardium and the
base. We therefore add a dashpot on the epicardium as well, a Robin
condition with damping=True, as in Monolithic 3D-0D coupling: an LV in a closed circulation.
It resists the normal velocity of the epicardium with 5e3 Pa s/m, much as
the pericardium and the surrounding tissue do. Like the viscous term, it
acts on a velocity, so the static prestressing solve ignores it.
def setup_problem(geometry, f0, s0, material_params):
material = pulse.HolzapfelOgden(f0=f0, s0=s0, **material_params)
Ta = pulse.Variable(
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0)),
"kPa",
)
active_model = pulse.ActiveStress(f0, activation=Ta)
model = pulse.CardiacModel(
material=material,
active=active_model,
compressibility=pulse.Compressible(),
viscoelasticity=pulse.Viscous(),
)
alpha_epi = pulse.Variable(
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1e5)),
"Pa / m",
)
robin_epi = pulse.RobinBC(value=alpha_epi, marker=geometry.markers["EPI"][0])
alpha_base = pulse.Variable(
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1e6)),
"Pa / m",
)
robin_base = pulse.RobinBC(value=alpha_base, marker=geometry.markers["BASE"][0])
beta_epi = pulse.Variable(
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(5e3)),
"Pa s / m",
)
damping_epi = pulse.RobinBC(
value=beta_epi, marker=geometry.markers["EPI"][0], damping=True,
)
robin = [robin_epi, robin_base, damping_epi]
# Dirichlet BC: Sliding Base (ux=0)
def dirichlet_bc(V: dolfinx.fem.FunctionSpace):
facets = geometry.facet_tags.find(geometry.markers["BASE"][0])
dofs = dolfinx.fem.locate_dofs_topological(V.sub(0), 2, facets)
return [dolfinx.fem.dirichletbc(0.0, dofs, V.sub(0))]
return model, robin, dirichlet_bc, Ta
material_params = pulse.HolzapfelOgden.transversely_isotropic_parameters()
model, robin, dirichlet_bc, Ta = setup_problem(
geometry=geometry,
f0=geo.f0,
s0=geo.s0,
material_params=material_params,
)
Prestressing (inverse elasticity)#
The imaged mesh is the shape of the heart at end diastole, under the
end-diastolic pressures. We prestress to the same pressures that PRELOAD
ramps up to, p_end_diastole of each ventricle. TargetPressure takes the
value in the unit of its traction, kilopascals here.
p_LV_ED = lv_params.p_end_diastole / 1e3 # Pa -> kPa
p_RV_ED = rv_params.p_end_diastole / 1e3 # Pa -> kPa
Since we want to apply pressures on both ventricles, we create two Neumann BCs.
pressure_lv = pulse.Variable(dolfinx.fem.Constant(geometry.mesh, 0.0), "kPa")
pressure_rv = pulse.Variable(dolfinx.fem.Constant(geometry.mesh, 0.0), "kPa")
neumann_lv = pulse.NeumannBC(traction=pressure_lv, marker=geometry.markers["LV"][0])
neumann_rv = pulse.NeumannBC(traction=pressure_rv, marker=geometry.markers["RV"][0])
bcs_prestress = pulse.BoundaryConditions(
robin=robin,
dirichlet=(dirichlet_bc,),
neumann=(neumann_lv, neumann_rv),
)
We store the prestressed displacement in a file to avoid recomputing it. The
file name carries the two target pressures, so changing p_end_diastole of
either ventricle prestresses again rather than reusing a reference unloaded
from different pressures.
prestress_fname = outdir / (
f"prestress_biv_inverse_LV{lv_params.p_end_diastole:g}Pa"
f"_RV{rv_params.p_end_diastole:g}Pa.bp"
)
if not prestress_fname.exists():
logger.info(
f"Start prestressing... Targets: p_LV={p_LV_ED:.2f} kPa, p_RV={p_RV_ED:.2f} kPa",
)
prestress_problem = pulse.unloading.PrestressProblem(
geometry=geometry,
model=model,
bcs=bcs_prestress,
parameters={"u_space": "P_2", "mesh_unit": mesh_unit},
targets=[
pulse.unloading.TargetPressure(
traction=pressure_lv, target=p_LV_ED, name="LV",
),
pulse.unloading.TargetPressure(
traction=pressure_rv, target=p_RV_ED, name="RV",
),
],
ramp_steps=20,
)
u_pre = prestress_problem.unload()
io4dolfinx.write_function_on_input_mesh(
prestress_fname, u_pre, time=0.0, name="u_pre",
)
with dolfinx.io.VTXWriter(
comm,
outdir / "prestress_biv_backward.bp",
[u_pre],
engine="BP4",
) as vtx:
vtx.write(0.0)
INFO:pulse:Start prestressing... Targets: p_LV=1.00 kPa, p_RV=0.40 kPa
INFO:pulse.unloading:Ramping LV traction to 0.0000
INFO:pulse.unloading:Ramping RV traction to 0.0000
INFO:pulse.unloading:Ramping LV traction to 0.0526
INFO:pulse.unloading:Ramping RV traction to 0.0211
INFO:pulse.unloading:Ramping LV traction to 0.1053
INFO:pulse.unloading:Ramping RV traction to 0.0421
INFO:pulse.unloading:Ramping LV traction to 0.1579
INFO:pulse.unloading:Ramping RV traction to 0.0632
INFO:pulse.unloading:Ramping LV traction to 0.2105
INFO:pulse.unloading:Ramping RV traction to 0.0842
INFO:pulse.unloading:Ramping LV traction to 0.2632
INFO:pulse.unloading:Ramping RV traction to 0.1053
INFO:pulse.unloading:Ramping LV traction to 0.3158
INFO:pulse.unloading:Ramping RV traction to 0.1263
INFO:pulse.unloading:Ramping LV traction to 0.3684
INFO:pulse.unloading:Ramping RV traction to 0.1474
INFO:pulse.unloading:Ramping LV traction to 0.4211
INFO:pulse.unloading:Ramping RV traction to 0.1684
INFO:pulse.unloading:Ramping LV traction to 0.4737
INFO:pulse.unloading:Ramping RV traction to 0.1895
INFO:pulse.unloading:Ramping LV traction to 0.5263
INFO:pulse.unloading:Ramping RV traction to 0.2105
INFO:pulse.unloading:Ramping LV traction to 0.5789
INFO:pulse.unloading:Ramping RV traction to 0.2316
INFO:pulse.unloading:Ramping LV traction to 0.6316
INFO:pulse.unloading:Ramping RV traction to 0.2526
INFO:pulse.unloading:Ramping LV traction to 0.6842
INFO:pulse.unloading:Ramping RV traction to 0.2737
INFO:pulse.unloading:Ramping LV traction to 0.7368
INFO:pulse.unloading:Ramping RV traction to 0.2947
INFO:pulse.unloading:Ramping LV traction to 0.7895
INFO:pulse.unloading:Ramping RV traction to 0.3158
INFO:pulse.unloading:Ramping LV traction to 0.8421
INFO:pulse.unloading:Ramping RV traction to 0.3368
INFO:pulse.unloading:Ramping LV traction to 0.8947
INFO:pulse.unloading:Ramping RV traction to 0.3579
INFO:pulse.unloading:Ramping LV traction to 0.9474
INFO:pulse.unloading:Ramping RV traction to 0.3789
INFO:pulse.unloading:Ramping LV traction to 1.0000
INFO:pulse.unloading:Ramping RV traction to 0.4000
The unloaded reference configuration#
V = dolfinx.fem.functionspace(geometry.mesh, ("Lagrange", 2, (3,)))
u_pre = dolfinx.fem.Function(V)
io4dolfinx.read_function(prestress_fname, u_pre, time=0.0, name="u_pre")
We use the prestressed displacement to deform the mesh to the reference configuration.
logger.info("Deforming mesh to Reference Configuration...")
geometry.deform(u_pre)
INFO:pulse:Deforming mesh to Reference Configuration...
We now map the fiber fields to the reference configuration. The solve uses the quadrature fibres; the DG 1 fibres are mapped too, for post-processing fibre stress or strain.
logger.info("Mapping fibers to Reference Configuration...")
f0_quad = pulse.utils.map_vector_field(
f=geo.f0, u=u_pre, normalize=True, name="f0_unloaded",
)
s0_quad = pulse.utils.map_vector_field(
f=geo.s0, u=u_pre, normalize=True, name="s0_unloaded",
)
f0_map = pulse.utils.map_vector_field(
geo.additional_data["f0_DG_1"],
u=u_pre,
normalize=True,
name="f0",
)
INFO:pulse:Mapping fibers to Reference Configuration...
Calculate unloaded volumes
lvv_unloaded = comm.allreduce(geometry.volume("LV"), op=MPI.SUM)
rvv_unloaded = comm.allreduce(geometry.volume("RV"), op=MPI.SUM)
logger.info(
f"Unloaded volumes: LV={lvv_unloaded / mL:.2f} mL, RV={rvv_unloaded / mL:.2f} mL",
)
model, robin, dirichlet_bc, Ta = setup_problem(
geometry=geometry,
f0=f0_quad,
s0=s0_quad,
material_params=material_params,
)
bcs_forward = pulse.BoundaryConditions(robin=robin, dirichlet=(dirichlet_bc,))
INFO:pulse:Unloaded volumes: LV=80.90 mL, RV=62.95 mL
The coupled problem#
Each cavity gets a CavityControl rather than a volume. We put no Neumann
pressure on the endocardium: the load on the wall comes from the cavity’s
pressure unknown, whatever constraint it satisfies.
The problem starts from rest in the unloaded configuration, at zero cavity
pressure. initialize reads each cavity’s volume and pressure there and puts
both cavities in PRELOAD. That phase then ramps the pressure back up to
p_end_diastole, which is the pressure we unloaded from, so at
t_end_diastole the wall is back at (close to) the imaged end-diastolic
shape. This is why PRELOAD’s end pressure must match the prestress target:
any other value would start contraction from a different shape. The ramp
also takes the place of an explicit inflation to end diastole, such as the
one in Monolithic 3D-0D coupling: a UK Biobank biventricular mesh in a closed circulation.
problem = pulse.problem.DynamicProblem(
model=model,
geometry=geometry,
bcs=bcs_forward,
cavities=[
pulse.problem.Cavity(marker=name, control=pulse.problem.CavityControl(geometry.mesh))
for name in ("LV", "RV")
],
parameters={
"mesh_unit": mesh_unit,
"u_space": "P_2",
"rho": pulse.Variable(1e3, "kg/m^3"),
"dt": pulse.Variable(DT, "s"),
},
)
controller = cycle.CycleController(problem, params)
controller.initialize(t0=0.0)
Stepping#
The loop is plain Python. At each step we set the active tension, note the
phase each cavity is in, and ask the controller for one step. step sets
each cavity’s constraint, solves, and only then decides on the next phase,
so the phase we note before the call is the one the step is solved under. It
returns False if Newton did not converge, even after one retry, and then
leaves the state as it was before the call.
pulse.cycle logs every phase transition. After each step,
controller.records holds each cavity’s volume V, pressure P, compliance
pressure P_c and outflow Q, all in SI units.
We also keep the moving geometry every few steps so that make_animations.py
can render it afterwards. FrameRecorder keeps each rank’s part of the mesh
separately, so we only record frames in serial runs, and never under CI,
where the run is two steps rather than two beats. The video on the page comes
from a saved run instead.
history: dict[str, list[float]] = {
k: []
for k in (
"time", "V_LV", "V_RV", "p_LV", "p_RV", "Pc_LV", "Pc_RV",
"Q_LV", "Q_RV", "phase_LV", "phase_RV", "Ta_LV",
)
}
recorder = animation.FrameRecorder(
geometry.mesh, every=5, enabled=not IN_CI and comm.size == 1, up=up,
)
vtx = dolfinx.io.VTXWriter(comm, outdir / "displacement.bp", [problem.u], engine="BP4")
logger.info(f"Running {max_steps} steps of {DT * 1e3:.0f} ms...")
t = 0.0
for step in range(1, max_steps + 1):
t = step * DT
Ta.assign(activation(t))
solved_under = {name: controller.cycles[name].phase for name in params}
if not controller.step(t, DT):
raise RuntimeError(f"Step to t={t:.4f} s did not converge (phases {solved_under})")
history["time"].append(t)
history["Ta_LV"].append(activation(t))
for name, record in controller.records.items():
history[f"V_{name}"].append(record.V / mL)
history[f"p_{name}"].append(record.P / mmHg)
history[f"Pc_{name}"].append(record.P_c / mmHg)
history[f"Q_{name}"].append(record.Q / mL)
history[f"phase_{name}"].append(int(solved_under[name]))
recorder.record(problem.u, t, step)
if step % 10 == 0 or step == max_steps:
vtx.write(t)
vtx.close()
INFO:pulse:Running 2 steps of 1 ms...
logger.info("Simulation complete.")
INFO:pulse:Simulation complete.
Results#
We save the traces (volumes in mL, pressures in mmHg, outflows in mL/s and the phase each step was solved under) and the recorded frames, and plot them. The phase intervals of the left ventricle are shaded behind its traces.
saved = recorder.save(outdir / "frames.npz")
if saved is not None:
logger.info(f"Saved {len(recorder.times)} frames of the moving geometry to {saved}")
if comm.rank == 0:
traces = {k: np.asarray(v) for k, v in history.items()}
np.savez(outdir / "traces.npz", **traces)
fig, axes = plt.subplots(2, 2, figsize=(11, 8), layout="constrained")
ax_loop, ax_p, ax_v, ax_c = axes.flat
for name in ("LV", "RV"):
colour = animation.CHAMBER_COLOURS[name]
ax_loop.plot(
traces[f"V_{name}"], traces[f"p_{name}"], color=colour, label=name, linewidth=1.3,
)
ax_p.plot(traces["time"], traces[f"p_{name}"], color=colour, label=f"p {name}")
ax_v.plot(traces["time"], traces[f"V_{name}"], color=colour, label=f"V {name}")
ax_loop.set_xlabel("V [mL]")
ax_loop.set_ylabel("p [mmHg]")
ax_loop.set_title("Pressure-volume loops")
ax_loop.legend(frameon=False)
ax_p.set_xlabel("Time [s]")
ax_p.set_ylabel("p [mmHg]")
ax_p.set_title("Cavity pressures (LV phases shaded)")
ax_v.set_xlabel("Time [s]")
ax_v.set_ylabel("V [mL]")
ax_v.set_title("Cavity volumes (LV phases shaded)")
ax_c.plot(
traces["time"], traces["p_LV"], color=animation.CHAMBER_COLOURS["LV"], label="p LV",
)
ax_c.plot(traces["time"], traces["Pc_LV"], color="0.3", linestyle="--", label="$P_c$ LV")
ax_c.set_xlabel("Time [s]")
ax_c.set_ylabel("p [mmHg]")
ax_c.set_title("LV cavity and compliance pressure")
ax_c.legend(frameon=False)
for axis in (ax_p, ax_v):
animation._shade_phases(axis, traces["time"], traces["phase_LV"])
fig.savefig(outdir / "complete_cycle.png", dpi=140)
plt.show()
if not IN_CI:
# We read the end-diastolic and end-systolic volumes off the phase
# trace of the last complete beat, rather than taking the largest and
# smallest volume: filling at a prescribed rate can carry the volume
# past end diastole before the next beat starts. The first step solved
# under isovolumic contraction holds V at the end-diastolic volume, and
# the first one solved under isovolumic relaxation at the end-systolic
# volume. A run shorter than one beat reports the beat it is in.
time = traces["time"]
complete = int(np.floor(time[-1] / PERIOD + 1e-9))
beat = max(complete - 1, 0)
in_beat = (time > beat * PERIOD + 1e-9) & (time <= (beat + 1) * PERIOD + 1e-9)
label = f"beat {beat + 1}" + ("" if complete > 0 else " (incomplete)")
for name in ("LV", "RV"):
phase = traces[f"phase_{name}"][in_beat]
V_beat = traces[f"V_{name}"][in_beat]
volumes = {}
for key, value in (
("EDV", cycle.Phase.ISOVOLUMIC_CONTRACTION),
("ESV", cycle.Phase.ISOVOLUMIC_RELAXATION),
):
hits = np.flatnonzero(phase == int(value))
if hits.size:
volumes[key] = float(V_beat[hits[0]])
peak = float(traces[f"p_{name}"][in_beat].max())
if len(volumes) < 2:
missing = {"EDV", "ESV"} - set(volumes)
logger.info(
f"{name}, {label}: no {' or '.join(sorted(missing))} "
f"(the beat never reached that phase); peak {peak:.1f} mmHg",
)
continue
EDV, ESV = volumes["EDV"], volumes["ESV"]
logger.info(
f"{name}, {label}: EDV {EDV:.1f} mL, ESV {ESV:.1f} mL, "
f"EF {100 * (1 - ESV / EDV):.1f}%, peak {peak:.1f} mmHg",
)
A whole beat#
The figure and the video come from a full run kept in _static/, rather
than from the two steps this page takes under CI. To regenerate them, run
python3 complete_cycle.py
python3 make_animations.py complete_cycle
Fig. 3 Both ventricles over two beats, the second solid and the first faded. In
the second beat the left ventricle ejects 55.4 mL (EF 50%) over 210 ms
against a peak of 97.6 mmHg, and the right 28.7 mL (EF 40%) against
23.6 mmHg. The first beat peaks a little higher, at 99.2 and 24.7 mmHg,
since its Windkessels start from p_init rather than from the pressures
they have drained to. The first loops close through PRELOAD, which ramps the
pressure back up to end diastole; the second ones stay open because the run
stops at 1.6 s, partway through filling.#
Where to go next#
Here the active tension is a prescribed waveform, the same everywhere, and each ventricle ejects into its own Windkessel. Crossbridges at every quadrature point in a closed-loop circulation drives a biventricular ellipsoid with a crossbridge model and a different strength per region, coupled to a full closed-loop circulation, and Monolithic 3D-0D coupling: a UK Biobank biventricular mesh in a closed circulation solves this mesh and a closed-loop circulation together in a single Newton system.
References#
Julie Bestel, Frédérique Clément, and Michel Sorine. A biomechanical model of muscle contraction. In Medical Image Computing and Computer-Assisted Intervention–MICCAI 2001: 4th International Conference Utrecht, The Netherlands, October 14–17, 2001 Proceedings 4, 1159–1161. Springer, 2001. doi:10.1007/3-540-45468-3_143.
Ruben Doste, David Soto-Iglesias, Gabriel Bernardino, Alejandro Alcaine, Rafael Sebastian, Sophie Giffard-Roisin, Maxime Sermesant, Antonio Berruezo, Damian Sanchez-Quintana, and Oscar Camara. A rule-based method to model myocardial fiber orientation in cardiac biventricular geometries with outflow tracts. International Journal for Numerical Methods in Biomedical Engineering, 35(4):e3185, 2019. doi:10.1002/cnm.3185.