Klotz curve#
Inlate a geometry to a pressure using different material models, and compare with the Klotz curve.
import math
import matplotlib.pyplot as plt
import numpy as np
import pulse
try:
from dolfin_adjoint import Constant, DirichletBC
except ImportError:
from dolfin import Constant, DirichletBC
ED_pressure = 1.6 # kPa
def setup_material(material_model):
"""
Choose parameters based on
Hadjicharalambous, Myrianthi, et al. "Analysis of passive
cardiac constitutive laws for parameter estimation using 3D
tagged MRI." Biomechanics and modeling in mechanobiology 14.4
(2015): 807-828.
These parameters did not really match the Klotz curve here.
Perhaps they did some more tuning?
"""
if material_model == "guccione":
matparams = pulse.Guccione.default_parameters()
matparams["C"] = 0.18 # kPa
matparams["bf"] = 27.75
matparams["bt"] = 5.37
matparams["bfs"] = 2.445
material = pulse.Guccione(
parameters=matparams,
f0=geometry.f0,
s0=geometry.s0,
n0=geometry.n0,
)
elif material_model == "neo_hookean":
matparams = pulse.NeoHookean.default_parameters()
matparams["mu"] = 10.0 # kPa
material = pulse.NeoHookean(parameters=matparams)
elif material_model == "holzapfel_ogden":
matparams = pulse.HolzapfelOgden.default_parameters()
matparams["a"] = 4.0 # kPa
matparams["a_f"] = 10.0 # kPa
matparams["b"] = 5.0
matparams["b_f"] = 5.0
material = pulse.HolzapfelOgden(
parameters=matparams,
f0=geometry.f0,
s0=geometry.s0,
n0=geometry.n0,
)
return material
def klotz_curve(ED_pressure):
"""
EDPVR based on Klotz curve
Klotz, Stefan, et al. "Single-beat estimation of end-diastolic
pressure-volume relationship: a novel method with potential for
noninvasive application." American Journal of Physiology-Heart and
Circulatory Physiology 291.1 (2006): H403-H412.
"""
# Some point at the EDPVR line
Vm = 148.663
Pm = ED_pressure
# Some constants
An = 27.8
Bn = 2.76
# kpa to mmhg
Pm = Pm * 760 / 101.325
V0 = Vm * (0.6 - 0.006 * Pm)
V30 = V0 + (Vm - V0) / (Pm / An) ** (1.0 / Bn)
beta = math.log(Pm / 30.0) / math.log(Vm / V30)
alpha = 30.0 / V30**beta
# Unloaded volume (not used here)
# P_V0 = alpha * V0 ** beta
vs = [V0]
ps = [0.0]
for p in np.linspace(1.0, 12.0):
vi = (p / alpha) ** (1.0 / beta)
vs.append(vi)
ps.append(p * 101.325 / 760) # Convert from mmhg to kPa
return vs, ps
def fix_basal_plane(W):
V = W if W.sub(0).num_sub_spaces() == 0 else W.sub(0)
bc = DirichletBC(
V,
Constant((0.0, 0.0, 0.0)),
geometry.ffun,
geometry.markers["BASE"][0],
)
return bc
geometry = pulse.HeartGeometry.from_file(pulse.mesh_paths["simple_ellipsoid"])
geometry.mesh.coordinates()[:] *= 4.7
# geometry = pulse.geometries.prolate_ellipsoid_geometry(mesh_size_factor=3.0)
# geometry.mesh.coordinates()[:] *= 0.37
2025-06-06 10:58:51,933 [418] INFO pulse.geometry_utils:
Load mesh from h5
dirichlet_bc = [fix_basal_plane]
lvp = Constant(0.0)
lv_marker = geometry.markers["ENDO"][0]
lv_pressure = pulse.NeumannBC(traction=lvp, marker=lv_marker, name="lv")
neumann_bc = [lv_pressure]
bcs = pulse.BoundaryConditions(dirichlet=dirichlet_bc, neumann=neumann_bc)
fig, ax = plt.subplots()
for material_model in ["neo_hookean", "guccione", "holzapfel_ogden"]:
material = setup_material(material_model)
problem = pulse.MechanicsProblem(geometry, material, bcs)
pressures = [0.0]
volumes = [geometry.cavity_volume()]
for p in np.linspace(0, ED_pressure, 10)[1:]:
pulse.iterate.iterate(problem, lvp, p)
pressures.append(p)
volumes.append(geometry.cavity_volume(u=problem.state.split()[0]))
ax.plot(volumes, pressures, label=" ".join(material_model.split("_")))
# Reset pressure
lvp.assign(Constant(0.0))
2025-06-06 10:58:52,063 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 10:58:52,064 [418] INFO pulse.iterate: Current control: f_20 = 0.000
2025-06-06 10:58:52,065 [418] INFO pulse.iterate: Target: 0.178
2025-06-06 10:58:54,230 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 10:58:54,230 [418] INFO pulse.iterate: Current control: f_20 = 0.178
2025-06-06 10:58:54,231 [418] INFO pulse.iterate: Target: 0.356
2025-06-06 10:58:56,083 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 10:58:56,084 [418] INFO pulse.iterate: Current control: f_20 = 0.356
2025-06-06 10:58:56,085 [418] INFO pulse.iterate: Target: 0.533
2025-06-06 10:58:57,936 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 10:58:57,937 [418] INFO pulse.iterate: Current control: f_20 = 0.533
2025-06-06 10:58:57,937 [418] INFO pulse.iterate: Target: 0.711
2025-06-06 10:58:59,790 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 10:58:59,790 [418] INFO pulse.iterate: Current control: f_20 = 0.711
2025-06-06 10:58:59,791 [418] INFO pulse.iterate: Target: 0.889
2025-06-06 10:59:01,640 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 10:59:01,641 [418] INFO pulse.iterate: Current control: f_20 = 0.889
2025-06-06 10:59:01,642 [418] INFO pulse.iterate: Target: 1.067
2025-06-06 10:59:03,494 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 10:59:03,495 [418] INFO pulse.iterate: Current control: f_20 = 1.067
2025-06-06 10:59:03,496 [418] INFO pulse.iterate: Target: 1.244
2025-06-06 10:59:05,350 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 10:59:05,351 [418] INFO pulse.iterate: Current control: f_20 = 1.244
2025-06-06 10:59:05,352 [418] INFO pulse.iterate: Target: 1.422
2025-06-06 10:59:07,254 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 10:59:07,255 [418] INFO pulse.iterate: Current control: f_20 = 1.422
2025-06-06 10:59:07,255 [418] INFO pulse.iterate: Target: 1.600
2025-06-06 10:59:09,179 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 10:59:09,180 [418] INFO pulse.iterate: Current control: f_20 = 0.000
2025-06-06 10:59:09,180 [418] INFO pulse.iterate: Target: 0.178
2025-06-06 10:59:15,202 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 10:59:15,202 [418] INFO pulse.iterate: Current control: f_20 = 0.178
2025-06-06 10:59:15,203 [418] INFO pulse.iterate: Target: 0.356
2025-06-06 10:59:20,260 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 10:59:20,261 [418] INFO pulse.iterate: Current control: f_20 = 0.356
2025-06-06 10:59:20,261 [418] INFO pulse.iterate: Target: 0.533
2025-06-06 10:59:25,314 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 10:59:25,315 [418] INFO pulse.iterate: Current control: f_20 = 0.533
2025-06-06 10:59:25,315 [418] INFO pulse.iterate: Target: 0.711
2025-06-06 10:59:29,870 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 10:59:29,870 [418] INFO pulse.iterate: Current control: f_20 = 0.711
2025-06-06 10:59:29,871 [418] INFO pulse.iterate: Target: 0.889
2025-06-06 10:59:34,421 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 10:59:34,421 [418] INFO pulse.iterate: Current control: f_20 = 0.889
2025-06-06 10:59:34,422 [418] INFO pulse.iterate: Target: 1.067
2025-06-06 10:59:38,480 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 10:59:38,481 [418] INFO pulse.iterate: Current control: f_20 = 1.067
2025-06-06 10:59:38,481 [418] INFO pulse.iterate: Target: 1.244
2025-06-06 10:59:42,536 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 10:59:42,537 [418] INFO pulse.iterate: Current control: f_20 = 1.244
2025-06-06 10:59:42,538 [418] INFO pulse.iterate: Target: 1.422
2025-06-06 10:59:46,616 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 10:59:46,617 [418] INFO pulse.iterate: Current control: f_20 = 1.422
2025-06-06 10:59:46,618 [418] INFO pulse.iterate: Target: 1.600
2025-06-06 10:59:50,760 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 10:59:50,761 [418] INFO pulse.iterate: Current control: f_20 = 0.000
2025-06-06 10:59:50,761 [418] INFO pulse.iterate: Target: 0.178
2025-06-06 10:59:55,055 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 10:59:55,056 [418] INFO pulse.iterate: Current control: f_20 = 0.178
2025-06-06 10:59:55,056 [418] INFO pulse.iterate: Target: 0.356
2025-06-06 10:59:59,420 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 10:59:59,421 [418] INFO pulse.iterate: Current control: f_20 = 0.356
2025-06-06 10:59:59,421 [418] INFO pulse.iterate: Target: 0.533
2025-06-06 11:00:03,302 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 11:00:03,303 [418] INFO pulse.iterate: Current control: f_20 = 0.533
2025-06-06 11:00:03,303 [418] INFO pulse.iterate: Target: 0.711
2025-06-06 11:00:07,224 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 11:00:07,225 [418] INFO pulse.iterate: Current control: f_20 = 0.711
2025-06-06 11:00:07,225 [418] INFO pulse.iterate: Target: 0.889
2025-06-06 11:00:11,621 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 11:00:11,622 [418] INFO pulse.iterate: Current control: f_20 = 0.889
2025-06-06 11:00:11,623 [418] INFO pulse.iterate: Target: 1.067
2025-06-06 11:00:16,001 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 11:00:16,002 [418] INFO pulse.iterate: Current control: f_20 = 1.067
2025-06-06 11:00:16,002 [418] INFO pulse.iterate: Target: 1.244
2025-06-06 11:00:19,432 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 11:00:19,432 [418] INFO pulse.iterate: Current control: f_20 = 1.244
2025-06-06 11:00:19,433 [418] INFO pulse.iterate: Target: 1.422
2025-06-06 11:00:23,826 [418] INFO pulse.iterate: Iterating to target control (f_20)...
2025-06-06 11:00:23,827 [418] INFO pulse.iterate: Current control: f_20 = 1.422
2025-06-06 11:00:23,828 [418] INFO pulse.iterate: Target: 1.600
vs, ps = klotz_curve(ED_pressure)
ax.plot(vs, ps, linestyle="--", label="Klotz curve")
ax.legend(loc="best")
ax.set_xlabel("Volume (ml)")
ax.set_ylabel("Pressure (kPa)")
plt.show()
# plt.savefig("klotz_curve.png")