Klotz curve
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
2022-06-05 11:44:29,034 [791] 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))
2022-06-05 11:44:30,097 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:44:30,098 [791] INFO pulse.iterate: Current control: f_20 = 0.000
2022-06-05 11:44:30,099 [791] INFO pulse.iterate: Target: 0.178
2022-06-05 11:46:09,794 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:46:09,795 [791] INFO pulse.iterate: Current control: f_20 = 0.178
2022-06-05 11:46:09,796 [791] INFO pulse.iterate: Target: 0.356
2022-06-05 11:46:12,535 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:46:12,536 [791] INFO pulse.iterate: Current control: f_20 = 0.356
2022-06-05 11:46:12,537 [791] INFO pulse.iterate: Target: 0.533
2022-06-05 11:46:15,290 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:46:15,291 [791] INFO pulse.iterate: Current control: f_20 = 0.533
2022-06-05 11:46:15,291 [791] INFO pulse.iterate: Target: 0.711
2022-06-05 11:46:18,039 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:46:18,039 [791] INFO pulse.iterate: Current control: f_20 = 0.711
2022-06-05 11:46:18,042 [791] INFO pulse.iterate: Target: 0.889
2022-06-05 11:46:20,764 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:46:20,765 [791] INFO pulse.iterate: Current control: f_20 = 0.889
2022-06-05 11:46:20,767 [791] INFO pulse.iterate: Target: 1.067
2022-06-05 11:46:23,523 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:46:23,524 [791] INFO pulse.iterate: Current control: f_20 = 1.067
2022-06-05 11:46:23,526 [791] INFO pulse.iterate: Target: 1.244
2022-06-05 11:46:26,265 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:46:26,266 [791] INFO pulse.iterate: Current control: f_20 = 1.244
2022-06-05 11:46:26,267 [791] INFO pulse.iterate: Target: 1.422
2022-06-05 11:46:29,007 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:46:29,008 [791] INFO pulse.iterate: Current control: f_20 = 1.422
2022-06-05 11:46:29,010 [791] INFO pulse.iterate: Target: 1.600
2022-06-05 11:46:31,840 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:46:31,841 [791] INFO pulse.iterate: Current control: f_20 = 0.000
2022-06-05 11:46:31,841 [791] INFO pulse.iterate: Target: 0.178
2022-06-05 11:49:29,654 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:49:29,655 [791] INFO pulse.iterate: Current control: f_20 = 0.178
2022-06-05 11:49:29,655 [791] INFO pulse.iterate: Target: 0.356
2022-06-05 11:49:36,449 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:49:36,450 [791] INFO pulse.iterate: Current control: f_20 = 0.356
2022-06-05 11:49:36,451 [791] INFO pulse.iterate: Target: 0.533
2022-06-05 11:49:43,255 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:49:43,256 [791] INFO pulse.iterate: Current control: f_20 = 0.533
2022-06-05 11:49:43,258 [791] INFO pulse.iterate: Target: 0.711
2022-06-05 11:49:49,453 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:49:49,454 [791] INFO pulse.iterate: Current control: f_20 = 0.711
2022-06-05 11:49:49,454 [791] INFO pulse.iterate: Target: 0.889
2022-06-05 11:49:55,643 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:49:55,644 [791] INFO pulse.iterate: Current control: f_20 = 0.889
2022-06-05 11:49:55,645 [791] INFO pulse.iterate: Target: 1.067
2022-06-05 11:50:00,778 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:50:00,778 [791] INFO pulse.iterate: Current control: f_20 = 1.067
2022-06-05 11:50:00,779 [791] INFO pulse.iterate: Target: 1.244
2022-06-05 11:50:06,224 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:50:06,225 [791] INFO pulse.iterate: Current control: f_20 = 1.244
2022-06-05 11:50:06,227 [791] INFO pulse.iterate: Target: 1.422
2022-06-05 11:50:11,378 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:50:11,379 [791] INFO pulse.iterate: Current control: f_20 = 1.422
2022-06-05 11:50:11,379 [791] INFO pulse.iterate: Target: 1.600
2022-06-05 11:50:16,945 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:50:16,945 [791] INFO pulse.iterate: Current control: f_20 = 0.000
2022-06-05 11:50:16,946 [791] INFO pulse.iterate: Target: 0.178
2022-06-05 11:52:24,492 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:52:24,493 [791] INFO pulse.iterate: Current control: f_20 = 0.178
2022-06-05 11:52:24,494 [791] INFO pulse.iterate: Target: 0.356
2022-06-05 11:52:30,262 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:52:30,262 [791] INFO pulse.iterate: Current control: f_20 = 0.356
2022-06-05 11:52:30,264 [791] INFO pulse.iterate: Target: 0.533
2022-06-05 11:52:35,404 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:52:35,405 [791] INFO pulse.iterate: Current control: f_20 = 0.533
2022-06-05 11:52:35,407 [791] INFO pulse.iterate: Target: 0.711
2022-06-05 11:52:40,546 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:52:40,547 [791] INFO pulse.iterate: Current control: f_20 = 0.711
2022-06-05 11:52:40,550 [791] INFO pulse.iterate: Target: 0.889
2022-06-05 11:52:46,407 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:52:46,408 [791] INFO pulse.iterate: Current control: f_20 = 0.889
2022-06-05 11:52:46,409 [791] INFO pulse.iterate: Target: 1.067
2022-06-05 11:52:51,524 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:52:51,525 [791] INFO pulse.iterate: Current control: f_20 = 1.067
2022-06-05 11:52:51,527 [791] INFO pulse.iterate: Target: 1.244
2022-06-05 11:52:56,638 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:52:56,639 [791] INFO pulse.iterate: Current control: f_20 = 1.244
2022-06-05 11:52:56,639 [791] INFO pulse.iterate: Target: 1.422
2022-06-05 11:53:02,369 [791] INFO pulse.iterate: Iterating to target control (f_20)...
2022-06-05 11:53:02,370 [791] INFO pulse.iterate: Current control: f_20 = 1.422
2022-06-05 11:53:02,372 [791] 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")