LV ellipsoid with time dependent pressure and activation

In this example we will solve a time dependent mechanics problem for the left ventricle ellipsoid geometry.

The equations of motion can be expressed using the following notation in the reference configuration

(1)\[\begin{split} \begin{aligned} \rho \ddot{u} - \nabla\cdot P &= 0, \mbox{ in } \Omega ,\\ PN &= pJF^{-T}N, \mbox{ on } \Gamma_{\rm endo}, \\ PN\cdot N + \alpha_{\rm epi}u\cdot N + \beta_{\rm epi}\dot{u}\cdot N &= 0, \mbox{ on } \Gamma_{\rm epi}, \\ PN\times N &=0, \mbox{ on } \Gamma_{\rm epi}, \\ PN + \alpha_{\rm top}u + \beta_{\rm top}\dot{u} &= 0, \mbox{ on } \Gamma_{\rm top}. \end{aligned} \end{split}\]

If we introduce the notation \(a= \ddot{u}, v=\dot{u}\) for the acceleration and velocity, respectively, a weak form of (1) can be written as

(2)\[\begin{split} \begin{aligned} \int_{\Omega} \rho a \cdot w \, \mathop{}\!\mathrm{d}{}X+ \int_{\Omega} P:Grad(w) \, \mathop{}\!\mathrm{d}{}X-\int_{\Gamma_{\rm endo}} p I JF^{-T}N \cdot w \, \mathop{}\!\mathrm{d}{}S\\ +\int_{\Gamma_{\rm epi}} \big(\alpha_{\rm epi} u \cdot N + \beta_{\rm epi} v \cdot N \big) w \cdot N \, \mathop{}\!\mathrm{d}{}S\\ +\int_{\Gamma_{\rm top}} \alpha_{\rm top} u \cdot w + \beta_{\rm top} v \cdot w \, \mathop{}\!\mathrm{d}{}S= 0 \quad \forall w \in H^1(\Omega). \end{aligned} \end{split}\]

In order to integrate (2) in time, we need to express \(a\) and \(v\) in terms of the displacement \(u\). For this we use the Generalized \(\alpha\)-method

The generalized \(\alpha\) method

The generalized \(\alpha\) or G-\(\alpha\) methods introduce additional parameters \(\alpha_f\) and \(\alpha_m\) to approximate \(v\) and \(a\) by evaluating the terms of the weak form at times \(t_{i+1}-\Delta t\alpha_m\) and \(t_{i+1}-\Delta t\alpha_f\). Specifically, the inertia term is evaluated at \(t_{i+1}-\Delta t\alpha_m\), and the other terms at \(t_{i+1}-\Delta t\alpha_f\). The weak form becomes

(3)\[\begin{split} \begin{aligned} \int_{\Omega} \rho a_{i+1-\alpha_m} \cdot w \, \mathrm{d}X + \int_{\Omega} P_{i+1-\alpha_f}:Grad(w) \, \mathrm{d}X -\int_{\Gamma_{\rm endo}} p I JF_{i+1-\alpha_f}^{-T}N \cdot w \, dS \\ +\int_{\Gamma_{\rm epi}} \big(\alpha_{\rm epi} u_{i+1-\alpha_f} \cdot N + \beta_{\rm epi} v_{i+1-\alpha_f} \cdot N \big) w \cdot N \, dS \\ +\int_{\Gamma_{\rm top}} \alpha_{\rm top} u_{i+1-\alpha_f} \cdot w + \beta_{\rm top} v_{i+1-\alpha_f} \cdot w \, dS = 0 \quad \forall w \in H^1(\Omega), \end{aligned} \end{split}\]

with

\[\begin{split} \begin{align*} u_{i+1-\alpha_f} &= (1-\alpha_f)u_{i+1}-\alpha_f u_i, \\ v_{i+1-\alpha_f} &= (1-\alpha_f)v_{i+1}-\alpha_f v_i, \\ a_{i+1-\alpha_m} &= (1-\alpha_m)a_{i+1}-\alpha_m a_i, \end{align*} \end{split}\]

\(v_{i+1},a_{i+1}\) given by

(4)\[ v_{i+1} = v_i + (1-\gamma) \Delta t ~ a_i + \gamma \Delta t ~ a_{i+1} \]

(5)\[ a_{i+1} = \frac{u_{i+1} - (u_i + \Delta t ~ v_i + (0.5 - \beta) \Delta t^2 ~ a_i)}{\beta \Delta t^2} \]

and

\[\begin{split} \begin{align*} F_{i+1-\alpha_f} &= I + \nabla u_{i+1-\alpha_f}, \\ P_{i+1-\alpha_f} &= P(u_{i+1-\alpha_f}, v_{i+1-\alpha_f}). \end{align*} \end{split}\]

Different choices of the four parameters \(\alpha_m, \alpha_f, \beta, \gamma\) yield methods with different accuracy and stability properties. Tables 1–3 in [EBB02] provides an overview of parameter choices for methods in the literature, as well as conditions for stability and convergence. We have used the choice \(\alpha_m =0.2, \alpha_f=0.4\), and

\[\begin{split} \begin{align*} \gamma &= 1/2 + \alpha_f-\alpha_m ,\\ \beta &= \frac{(\gamma + 1/2)^2}{4} . \end{align*} \end{split}\]

For this choice the solver converges through the time interval of interest, and the convergence is second order.

Pressure and activation model

We will use a pressure and activation model from [BClementS01] to drive the simulation. We consider a time-dependent pressure derived from the Bestel model. The solution \(p = p(t)\) is characterized as solution to the evolution equation

\[ \dot{p}(t) = -|b(t)|p(t) + \sigma_{\mathrm{mid}}|b(t)|_+ + \sigma_{\mathrm{pre}}|g_{\mathrm{pre}}(t)| \]

with \(b(\cdot)\) being the activation function described below:

\[\begin{split} b(t) =& a_{\mathrm{pre}}(t) + \alpha_{\mathrm{pre}}g_{\mathrm{pre}}(t) + \alpha_{\mathrm{mid}} \\ a_{\mathrm{pre}}(t) :=& \alpha_{\mathrm{max}} \cdot f_{\mathrm{pre}}(t) + \alpha_{\mathrm{min}} \cdot (1 - f_{\mathrm{pre}}(t)) \\ f_{\mathrm{pre}}(t) =& S^+(t - t_{\mathrm{sys}-\mathrm{pre}}) \cdot S^-(t t_{\mathrm{dias} - \mathrm{pre}}) \\ g_{\mathrm{pre}}(t) =& S^-(t - t_{\mathrm{dias} - \mathrm{pre}}) \end{split}\]

with \(S^{\pm}\) given by

\[ S^{\pm}(\Delta t) = \frac{1}{2}(1 \pm \mathrm{tanh}(\frac{\Delta t}{\gamma})) \]

Similarly, the active stress is characterized through a time-dependent stress function \(\tau\) solution to the evolution equation

\[ \dot{\tau}(t) = -|a(t)|\tau(t) + \sigma_0|a(t)|_+ \]

with \(a(\cdot)\) being the activation function and \sigma_0 contractility, where each remaining term is described below:

\[\begin{split} |a(t)|_+ =& \mathrm{max}\{a(t), 0\} \\ a(t) :=& \alpha_{\mathrm{max}} \cdot f(t) + \alpha_{\mathrm{min}} \cdot (1 - f(t)) \\ f(t) =& S^+(t - t_{\mathrm{sys}}) \cdot S^-(t - t_{\mathrm{dias}}) \\ S^{\pm}(\Delta t) =& \frac{1}{2}(1 \pm \mathrm{tanh}(\frac{\Delta t}{\gamma})) \end{split}\]

Constitutive model

We will use a nearly incompressible, orthotropic and viscoelastic version of the Holzapfel Ogden model. The material parameters are taken from [HO09]. The anistropic material strain energy function is given by

\[ \Psi_{\text{aniso}} = \frac{a}{2 b} \left( e^{ b (I_1 - 3)} -1 \right) + \frac{a_f}{2 b_f} \mathcal{H}(I_{4\mathbf{f}_0} - 1) \left( e^{ b_f (I_{4\mathbf{f}_0} - 1)_+^2} -1 \right) + \frac{a_s}{2 b_s} \mathcal{H}(I_{4\mathbf{s}_0} - 1) \left( e^{ b_s (I_{4\mathbf{s}_0} - 1)_+^2} -1 \right) + \frac{a_{fs}}{2 b_{fs}} \left( e^{ b_{fs} I_{8 \mathbf{f}_0 \mathbf{s}_0}^2} -1 \right) \]

The viscoelastic energy is given by

\[ \Psi_{\text{visco}} = \frac{\eta}{2} \mathrm{tr} ( \dot{\mathbf{E}}^2 ) \]

where \(\dot{\mathbf{E}}\) is temporal derivative the Green-Lagrange strain tensor and \(\eta\) is the viscosity parameter. Here the temporal derivative is computed by first computing the temporal deformation gradient

\[ \dot{\mathbf{F}} = \dot{\mathbf{F}} = \dot{\mathbf{I} + \nabla \mathbf{u}} = \nabla \dot{\mathbf{u}} = \nabla \mathbf{v} \]

and then

\[\begin{split} \mathbf{l} = \dot{\mathbf{F}} \mathbf{F}^{-1} \\ \mathbf{d} = \frac{1}{2} ( \mathbf{l} + \mathbf{l}^T ) \\ \dot{\mathbf{E}} = \mathbf{F}^T \mathbf{d} \mathbf{F} \end{split}\]

see [Hol02] Equation 2.139 and 2.163 for more details.

For compressibility we use the following strain energy function

\[ \Psi_{\text{comp}} = \frac{\kappa}{4} \left( J^2 - 1 - 2 \mathrm{ln}(J) \right) \]

where \(J\) is the determinant of the deformation gradient.

To solve the problem we first import the necessary packages

from pathlib import Path
import logging
import math
import os
from mpi4py import MPI
import dolfinx
from dolfinx import log
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import circulation.bestel
import cardiac_geometries
import cardiac_geometries.geometry
import fenicsx_pulse

Next we set up the logging and the MPI communicator

logging.basicConfig(level=logging.INFO)
comm = MPI.COMM_WORLD

and create an output directory

outdir = Path("time-dependent-bestel-lv")
outdir.mkdir(parents=True, exist_ok=True)

Next we create the geometry, which is an ellipsoid with a fiber field varying from -60 to 60 degrees in the endocardium and epicardium and interpolated into a Quadrature_6 function space

geodir = outdir / "geometry"
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,
    )

2025-03-05 19:48:31 [debug    ] Convert file time-dependent-bestel-lv/geometry/lv_ellipsoid.msh to dolfin

INFO:cardiac_geometries.geometry:Reading geometry from time-dependent-bestel-lv/geometry

Info    : Reading 'time-dependent-bestel-lv/geometry/lv_ellipsoid.msh'...
Info    : 54 entities
Info    : 191 nodes
Info    : 990 elements
Info    : Done reading 'time-dependent-bestel-lv/geometry/lv_ellipsoid.msh'

We load the geometry and convert it to the format used by fenicsx-pulse.

geo = cardiac_geometries.geometry.Geometry.from_folder(
    comm=comm,
    folder=geodir,
)
geometry = fenicsx_pulse.HeartGeometry.from_cardiac_geometries(geo, metadata={"quadrature_degree": 6})

INFO:cardiac_geometries.geometry:Reading geometry from time-dependent-bestel-lv/geometry

Next we create the material object using the class Holzapfel Ogden model

material_params = fenicsx_pulse.HolzapfelOgden.orthotropic_parameters()
material = fenicsx_pulse.HolzapfelOgden(f0=geo.f0, s0=geo.s0, **material_params)  # type: ignore
print(material)

a/2b (exp(b(I1 - 3)) - 1) + af/2bf H(I4f - 1) (exp(bf (I4f - 1)_+^2) - 1) + as/2bs H(I4s - 1) (exp(bs (I4s - 1)_+^2) - 1) + afs/2bfs (exp(bfs I8fs^2) - 1)

We set up the active stress model

Ta = fenicsx_pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0)), "Pa")
active_model = fenicsx_pulse.ActiveStress(geo.f0, activation=Ta)
print(active_model)

Ta (I4f - 1 + η ((I1 - 3) - (I4f - 1)))

and define the compressibility

comp_model = fenicsx_pulse.compressibility.Compressible2()
print(comp_model)

κ (J ** 2 - 1 - 2 ln(J))

and viscoelasticity model

viscoelastic_model = fenicsx_pulse.viscoelasticity.Viscous()
print(viscoelastic_model)

0.5η tr (E_dot E_dot)

Finally we assembles the CardiacModel

model = fenicsx_pulse.CardiacModel(
    material=material,
    active=active_model,
    compressibility=comp_model,
    viscoelasticity=viscoelastic_model,
)

Next we define the boundary conditions. First we define a traction on the endocardium

traction = fenicsx_pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0)), "Pa")
neumann = fenicsx_pulse.NeumannBC(traction=traction, marker=geometry.markers["ENDO"][0])

and Robin boundary conditions on the epicardium and base

alpha_epi = fenicsx_pulse.Variable(
        dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1e8)), "Pa / m",
)
robin_epi_u = fenicsx_pulse.RobinBC(value=alpha_epi, marker=geometry.markers["EPI"][0])
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)

alpha_base = fenicsx_pulse.Variable(
    dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1e5)), "Pa / m",
)
robin_base_u = fenicsx_pulse.RobinBC(value=alpha_base, marker=geometry.markers["BASE"][0])
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)

We assemble the boundary conditions

bcs = fenicsx_pulse.BoundaryConditions(robin=(robin_epi_u, robin_epi_v, robin_base_u, robin_base_v), neumann=(neumann,))

Finally we create a DynamicProblem

problem = fenicsx_pulse.problem.DynamicProblem(model=model, geometry=geometry, bcs=bcs, parameters={"base_bc": fenicsx_pulse.problem.BaseBC.free})

Note that we also specify that the base is free to move, meaning that there will be no Dirichlet boundary conditions on the base. Now we can do an initial solve the problem

log.set_log_level(log.LogLevel.INFO)
problem.solve()

INFO:scifem.solvers:Newton iteration 1: r (abs) = 1.0145062213561107e-05 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 2: r (abs) = 4.5704973235210315e-08 (tol=1e-06), r (rel) = 0.004505144697300679 (tol=1e-10)

2

The next step is to get the activation and pressure as a function of time. For this we use the time step from the problem parameters

dt = problem.parameters["dt"].to_base_units()
times = np.arange(0.0, 1.0, dt)

We solve the Bestel model for the pressure and activation which is already implemented in the circulation package

pressure_model = circulation.bestel.BestelPressure()
res = solve_ivp(
    pressure_model,
    [0.0, 1.0],
    [0.0],
    t_eval=times,
    method="Radau",
)
# Convert the pressure from Pa to kPa
pressure = res.y[0]

INFO:circulation.bestel:
 Bestel pressure model  
       parameters       
┏━━━━━━━━━━━━┳━━━━━━━━━┓
┃ Parameter  ┃ Value   ┃
┡━━━━━━━━━━━━╇━━━━━━━━━┩
│ t_sys_pre  │ 0.17    │
│ t_dias_pre │ 0.484   │
│ gamma      │ 0.005   │
│ a_max      │ 5.0     │
│ a_min      │ -30.0   │
│ alpha_pre  │ 5.0     │
│ alpha_mid  │ 1.0     │
│ sigma_pre  │ 7000.0  │
│ sigma_mid  │ 16000.0 │
└────────────┴─────────┘

activation_model = circulation.bestel.BestelActivation()
res = solve_ivp(
    activation_model,
    [0.0, 1.0],
    [0.0],
    t_eval=times,
    method="Radau",
)
# Convert the pressure from Pa to kPa
activation = res.y[0]

INFO:circulation.bestel:
Bestel activation model 
       parameters       
┏━━━━━━━━━━━┳━━━━━━━━━━┓
┃ Parameter ┃ Value    ┃
┡━━━━━━━━━━━╇━━━━━━━━━━┩
│ t_sys     │ 0.16     │
│ t_dias    │ 0.484    │
│ gamma     │ 0.005    │
│ a_max     │ 5.0      │
│ a_min     │ -30.0    │
│ sigma_0   │ 150000.0 │
└───────────┴──────────┘

Let us also plot the profiles of the pressure and activation

fig, ax = plt.subplots(2, 1, sharex=True, figsize=(10, 10))
ax[0].plot(times, pressure)
ax[0].set_title("Pressure")
ax[1].plot(times, activation)
ax[1].set_title("Activation")
fig.savefig(outdir / "pressure_activation.png")

New let us write the displacement to a file using the VTXWriter from dolfinx

vtx = dolfinx.io.VTXWriter(geometry.mesh.comm, outdir / "displacement.bp", [problem.u], engine="BP4")
vtx.write(0.0)

We can also calculate the volume of the left ventricle at each time step, so let us first define the volume form

volume_form = dolfinx.fem.form(geometry.volume_form(u=problem.u) * geometry.ds(geometry.markers["ENDO"][0]))
initial_volume = geo.mesh.comm.allreduce(dolfinx.fem.assemble_scalar(volume_form))
print(f"Initial volume: {initial_volume}")

Initial volume: 0.00013552928610385148[2025-03-05 19:48:39.251] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-03-05 19:48:39.251] [info] Requesting connectivity (3, 0) - (2, 0)

[2025-03-05 19:48:39.251] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-03-05 19:48:39.251] [info] Requesting connectivity (3, 0) - (2, 0)

and then loop over the time steps and solve the problem for each time step

volumes = []
for i, (tai, pi, ti) in enumerate(zip(activation, pressure, times)):
    print(f"Solving for time {ti}, activation {tai}, pressure {pi}")
    traction.assign(pi)
    Ta.assign(tai)
    problem.solve()
    vtx.write((i + 1) * dt)

    volumes.append(geo.mesh.comm.allreduce(dolfinx.fem.assemble_scalar(volume_form)))

    if geo.mesh.comm.rank == 0:
        # Plot data as we go
        fig, ax = plt.subplots(4, 1, figsize=(10, 10))
        ax[0].plot(times[:i + 1], pressure[:i + 1])
        ax[0].set_title("Pressure")
        ax[1].plot(times[:i + 1], activation[:i + 1])
        ax[1].set_title("Activation")
        ax[2].plot(times[:i + 1], volumes)
        ax[2].set_title("Volume")
        ax[3].plot(volumes, pressure[:i + 1])
        fig.savefig(outdir / "lv_ellipsoid_time_dependent_bestel.png")
        plt.close(fig)

    if os.getenv("CI") and i > 2:
        # Early stopping for CI
        break

INFO:scifem.solvers:Newton iteration 1: r (abs) = 1.051666490739738e-05 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)

Solving for time 0.0, activation 0.0, pressure 0.0

INFO:scifem.solvers:Newton iteration 2: r (abs) = 7.111185719151402e-09 (tol=1e-06), r (rel) = 0.0006761825903713469 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 1: r (abs) = 4.3792581856585495e-06 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)

Solving for time 0.001, activation 0.0, pressure 6.916668004930999

INFO:scifem.solvers:Newton iteration 2: r (abs) = 6.831272542970362e-09 (tol=1e-06), r (rel) = 0.001559915459047793 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 1: r (abs) = 7.400240303774599e-06 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 2: r (abs) = 7.220696891292544e-09 (tol=1e-06), r (rel) = 0.0009757381645579163 (tol=1e-10)

Solving for time 0.002, activation 0.0, pressure 13.668630239081294

INFO:scifem.solvers:Newton iteration 1: r (abs) = 1.4723190123813685e-05 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)

Solving for time 0.003, activation 0.0, pressure 20.261093815707657

INFO:scifem.solvers:Newton iteration 2: r (abs) = 1.4469828017039121e-08 (tol=1e-06), r (rel) = 0.0009827916297593163 (tol=1e-10)

References

[BClementS01]

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.

[EBB02]

Silvano Erlicher, Luca Bonaventura, and Oreste S Bursi. The analysis of the generalized-α method for non-linear dynamic problems. Computational mechanics, 28(2):83–104, 2002.

[Hol02]

Gerhard A Holzapfel. Nonlinear solid mechanics: a continuum approach for engineering science. 2002.

[HO09]

Gerhard A Holzapfel and Ray W Ogden. Constitutive modelling of passive myocardium: a structurally based framework for material characterization. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 367(1902):3445–3475, 2009.