Energy consistent boundary conditions for cardiac models

[dseyler]

Use Case

When simulating a ventricular model or any other truncated pressurized vessel, a Neumann pressure BC is applied to the endocardial surface to model the hydrostatic pressure of blood. This results in a net force on the body that should be counteracted by an equal and opposite force distributed over the boundary with the truncated portion of the domain.

Several papers [1] [2] have accounted for this pressure imbalance by implementing "energy consistent boundary conditions" that integrate the pressure over the Neumann face and apply a traction of equal and opposite magnitude to another specified face (usually the annuli or base). This traction force can be computed as:

or in the reference configuration:

This boundary condition is consistent with the physics of the system as the pressure of a "closed" chamber should contribute to the internal energy of the system, not the translational momentum of the myocardium or elastic energy of the surrounding tissue.

Alternatively, other papers [3] have simply placed a physical surface (not just RIS) representing a valve to reach pressure equilibrium.

Problem

Applying a Neumann pressure BC to the endocardium of a ventricular model or any other pressurized open surface results in a net force on the body which results in rigid-body translation if not otherwise constrained, especially during systolic contraction. This is inconsistent with the actual physics of the heart where the AV valves and their attachments (during systole) or atria (during diastole) should also experience an equal and opposite force, assuming the process is quasi-static.

We have usually handled this pressure imbalance by encasing the exterior surface of the model in stiff Robin BCs applied in the normal direction that are extremely stiff at the apex and much less stiff at the base (varying by ~4 orders of magnitude). This is analogous to putting a brick wall in front of our model and encasing it in Jello to prevent rigid-body motion, while allowing the model to deform in the regions with less stiff BCs. This is problematic for a number of reasons:

1. Simulated cardiac motion is far more dependent on the stiffness of Robin BCs than it is on the material model, cardiomyocyte orientation, or force of active contraction. We are essentially tuning the stiffness of the Jello that our model is encased in so that it deforms in a physiologically realistic way. This can have a considerable impact on the stress distribution in the model that may not be realistic. Such a large dependence on Robin BC stiffness is obviously not physical as Heart in a Box machines work fairly well without such a finely tuned Jello.

2. As a result of the excessively stiff Robin BCs, we have to apply a force of active contraction that is much higher in simulations than experimental results show for cardiomyocytes to achieve the same motion since most of the energy is spent resisting the Robin BC stiffness rather than compressing the fluid.

3. The force resisting basal motion should be coupled to the endocardial pressure (coming from 0D coupling), not modeled as a spring force that is merely proportional to basal displacement via a Robin BC. Currently, basal motion is highly dependent on the ratio between basal and apical BC stiffness which must be tuned for a given LPN and active stress magnitude. As a result, even if the Robin BCs can be tuned to reproduce physiological motion in one configuration, that tuning may not generalize. If we were to try to study how vascular remodeling (using a different LPN), infarction/cardiomyopathy (changing active stress magnitude), or cardiomyocyte orientation would affect simulations, the original tuned Robin BCs would likely no longer be valid, making any conclusions about simulated kinematics dependent on the inaccurate Robin BCs rather than the physiology of interest.

4. Regardless of whether stiff, spatially varying Robin BCs are the correct way to counteract the pressure imbalance, we don't have a very systematic way to tune them despite how much of an impact they have on cardiac motion.

Solution

I made a prototype for the energy consistent boundary condition feature in this branch. Currently the feature works for steady and unsteady Neumann BCs but is unstable for coupled BCs due to the incorrect tangent contribution.

On the front end, energy consistent BCs are specified in the XML file under a Neumann BC:

   <Add_BC name="endo" > 
      <Type> Neu </Type> 
      <Time_dependence> Unsteady </Time_dependence> 
      <Temporal_values_file_path> pressure.dat </Temporal_values_file_path> 
      <Energy_balance_face> top </Energy_balance_face>    <!-- Specify energy balance surface -->
      <Ramp_function> false </Ramp_function> 
      <Follower_pressure_load> true </Follower_pressure_load> 
   </Add_BC> 

The parser then creates a new traction BC on the specified face and links the Neumann BC to the traction BC via the variable iEnergyBalanceBC.

For each Neumann BC linked to an energy balance face, at each iteration, the solver:

1. integrates the product of pressure and the normal vectors over the Neumann face, F_neu
2. Computes the area of the energy balance face, A_eb
3. Computes the traction vector as t = F_neu/A_eb

The deformed configuration is used if folder pressure is true. Otherwise, the reference configuration is used.

This enforces global force balance between the pressurized surface and the truncated boundary.

See the video below where the simulation with an energy balanced BC (red) on the top surface shows significantly less rigid body motion than the unbalanced simulation (blue) given the same epicardial Robin BCs. To equivalently constrain the unbalanced model in space, the epicardial BCs needed to be ~1000x stiffer.

Alternatives considered

As mentioned above, one alternative is to simply add a slab of material over the open face rather than enforcing global force balance via a computed traction.

Just to weigh some pros and cons for each method:

Energy Consistent BCs

Pros

• Works for any mesh
• No need for extra mesh processing to add slab face
• Results "look nicer" - no slab over open face
• Easy to control where traction is distributed (annulus? papillary muscles? atria?)
• Easy to turn on/off
• Can be used in conjunction with RIS valves

Cons

• Balances force but not necessarily moments (could be implemented though)
• 0D coupling may be difficult to implement (likely similar to coupled Neumann implementation)
• Could become obsolete if we transition to whole heart simulations with 3D valve mechanics.

Slab over open face

Pros

• No changes to the solver
• No need for "3D-0D virtual cap coupling"
• Slab could have realistic valve material properties

Cons

• Can't turn off - slab material still influences mechanics when valve is open
• Slabs make fiber generation, universal ventricular coordinates difficult to compute
• Extra mesh processing step
• Less visually appealing

Would love to hear your opinions on whether this feature would be useful or whether I'm misunderstanding the impact of the unbalanced pressure contribution.

Additional context

No response

Code of Conduct

• I agree to follow this project's Code of Conduct and Contributing Guidelines

[mrp089]

Having pericardial conditions shouldn't require fixating the apex (that's kinda the point). An easy solution is to have stiffer springs at the base to represent the atria, see, e.g., this paper (Table 1, "spring stiffness).