Compressibility

Cardiac tissue is often modelled as a nearly incompressible or a fully incompressible material.

A nearly incompressible material can be modeled using a penalty formulation of volume changes. This might be considered the most physiological correct option. It has been demonstrated that the cells, which mostly consist of water, are around 100–1000 times more incompressible than the matrix surrounding them. On the other hand, a fully incompressible model is a common simplification – which has certain numerical advantages. Here, however, we don’t separate between the two subdomains – and we assume a continuous pressure field across the membrane.

In the papers we haven’t compared the two methods – but, as we will see later, it matters quite a bit, especially for the stress values. Let’s start, as usual, by importing the relevant libraries:

import numpy as np
import matplotlib.pyplot as plt

import fenics as f
import emimechanicalmodel as emi_m

f.set_log_level(30)        # less verbatim output from the Newton solver compared to default

… and load the meshes and geometries:

mesh_file = "tile_connected_10p0.h5"
mesh, volumes = emi_m.load_mesh(mesh_file)

Mesh and subdomains loaded successfully.
Number of nodes: 335, number of elements: 1364

Fundamental equations

When we solve the mechanical problem, we want to find a displacement field \(\mathbf{u}\) which gives us an equilibrium of stresses. We are, essentially, solving for a three-dimensional version of Newton’s third law, expressed by Piola-Kirchhoff stress tensor \(\mathbf{P}\):

\[\nabla \cdot \mathbf{P} = 0\]

Earlier we have seen how this can be decomposed into an active and a passive part, and how different material models can be used for the passive part. There is, however, one more term to take into account - a term for the compressibility model.

Using a nearly incompressible model, we calculate the volume change \(J = \mathrm{det} (\mathbf{F})\), and add a penalty term:

(14)\[\begin{equation} \psi* = \psi + \kappa (J \ln J - J + 1) \end{equation}\]

– and when calculating \(\mathbf{P}\) we use \(\psi*\) instead of \(\psi\). \(\kappa\) is a parameter which we can define differently over each subdomain, as \(\kappa_i\) and \(\kappa_e\). It has been demonstrated experimentally that the intracellular subdomain is 100–1000 times more compressible than the extracellular subdomain (citation needed).

Comparing the compressibility models:

We will try out the different compressibility models doing an active contraction experiment:

time = np.linspace(0, 500, 250)       # 500 ms with 250 steps

active_precomputed = emi_m.compute_active_component(time)
active_precomputed *= 2

plt.figure(figsize=(4, 3))
plt.plot(time, active_precomputed)
plt.xlabel("Time (ms)")
plt.ylabel("Active force (kPa); unscaled")
plt.show()

We will next define two instances of the model, one using a nearly incompressible and one using a fully incompressible approach – and track the fiber direction strain and stresses for comparison. Note that, as demonstrated here, you can specify the compressibility parameters for each subdomain, which only works using the nearly incompressible approach.

emi_model_incompressible = emi_m.EMIModel(
    mesh,
    volumes,
    experiment="contraction",
    compressibility_model="incompressible",
    )

emi_model_nearly_incompressible = emi_m.EMIModel(
    mesh,
    volumes,
    experiment="contraction",
    compressibility_model="nearly_incompressible",
    compressibility_parameters = {"kappa_e" : 100, "kappa_i" : 1000},
    )

fig, axes = plt.subplots(1, 2, figsize=(15, 3))

for emimodel in [emi_model_incompressible, emi_model_nearly_incompressible]:
    fiber_dir_strain_i = np.zeros_like(active_precomputed)
    fiber_dir_stress_i = np.zeros_like(active_precomputed) 
    subdomain_id = 1         # intracellular subdomain

    # then run the simulation
    for step, a_str in enumerate(active_precomputed):
        if step%10==0:
            print(f"Solving for step {step}/{len(active_precomputed)}")
        emimodel.update_active_fn(a_str)
        emimodel.solve()
    
        fiber_dir_strain_i[step] = \
            emimodel.evaluate_subdomain_strain_fibre_dir(subdomain_id)
        fiber_dir_stress_i[step] = \
            emimodel.evaluate_subdomain_stress_fibre_dir(subdomain_id)

    axes[0].plot(time, fiber_dir_strain_i)
    axes[1].plot(time, fiber_dir_stress_i)
        
axes[0].set_xlabel("Time (ms)")
axes[1].set_xlabel("Time (ms)")
axes[0].set_ylabel("Strain (-)")
axes[1].set_ylabel("Stress (kPa)")

axes[0].set_title(r"Fiber direction strain, over $\Omega_i$")
axes[1].set_title(r"Fiber direction stress, over $\Omega_i$")

axes[0].legend(["Fully incompressible model", "Nearly incompressible model"])
axes[1].legend(["Fully incompressible model", "Nearly incompressible model"])

plt.show()

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Domain length=102.0, width=20.0, height=24.0
Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Domain length=102.0, width=20.0, height=24.0

Solving for step 0/250
Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Solving for step 10/250

Solving for step 20/250

Solving for step 30/250

Solving for step 40/250

Solving for step 50/250

Solving for step 60/250

Solving for step 70/250

Solving for step 80/250

Solving for step 90/250

Solving for step 100/250

Solving for step 110/250

Solving for step 120/250

Solving for step 130/250

Solving for step 140/250

Solving for step 150/250

Solving for step 160/250

Solving for step 170/250

Solving for step 180/250

Solving for step 190/250

Solving for step 200/250

Solving for step 210/250

Solving for step 220/250

Solving for step 230/250

Solving for step 240/250

Solving for step 0/250
Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Solving for step 10/250

Solving for step 20/250

Solving for step 30/250

Solving for step 40/250

Solving for step 50/250

Solving for step 60/250

Solving for step 70/250

Solving for step 80/250

Solving for step 90/250

Solving for step 100/250

Solving for step 110/250

Solving for step 120/250

Solving for step 130/250

Solving for step 140/250

Solving for step 150/250

Solving for step 160/250

Solving for step 170/250

Solving for step 180/250

Solving for step 190/250

Solving for step 200/250

Solving for step 210/250

Solving for step 220/250

Solving for step 230/250

Solving for step 240/250

In similar manners we can track stresses and strain values in the extracellular subdomain, identified with the number \(0\). We also have in-built functions for tracking sheet and normal direction values, using functions evaluate_subdomain_strain_sheet_dir, evaluate_subdomain_strain_normal_dir, evaluate_subdomain_stress_sheet_dir, and evaluate_subdomain_stress_normal_dir.