The importance of the pericardium for cardiac biomechanics: From physiology to computational modeling

10/12/2018 ∙ by Martin R. Pfaller, et al. ∙ Technische Universität München 0

The human heart is enclosed in the pericardial cavity. The pericardium consists of a layered thin sac and is separated from the myocardium by a thin film of fluid. It provides a fixture in space and frictionless sliding of the myocardium. The influence of the pericardium is essential for predictive mechanical simulations of the heart. However, there is no consensus on physiologically correct and computationally tractable pericardial boundary conditions. Here we propose to model the pericardial influence as a parallel spring and dashpot acting in normal direction to the epicardium. Using a four-chamber geometry, we compare a model with pericardial boundary conditions to a model with fixated apex. The influence of pericardial stiffness is demonstrated in a parametric study. Comparing simulation results to measurements from cine magnetic resonance imaging reveals that adding pericardial boundary conditions yields a better approximation with respect to atrioventricular plane displacement, atrial filling, and overall spatial approximation error. We demonstrate that this simple model of pericardial-myocardial interaction can correctly predict the pumping mechanisms of the heart as previously assessed in clinical studies. Utilizing a pericardial model can not only provide much more realistic cardiac mechanics simulations but also allows new insights into pericardial-myocardial interaction which cannot be assessed in clinical measurements yet.



There are no comments yet.


page 2

page 4

page 7

page 12

page 15

page 18

page 23

page 26

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

1.1 Motivation

(a) Dissected mediastinum with cut pericardium and heart surface. Image by G. M. Gruber, Medical University of Vienna, Austria.
(b) Location of the heart with respect to serous and fibrous pericardium. Inspired by iaizzo15 .
(c) Cross-sectional view of transmural layers of heart and pericardium. Inspired by iaizzo15 .
Figure 1: Heart and pericardium.

Cardiac mechanics simulations consist of solving a nonlinear elastodynamic boundary value problem sainte-marie06 . Physiological boundary conditions are essential to achieve predictive results for any clinical purposes. The boundary conditions on the structure field of the myocardium are mainly governed by two physiological aspects: Blood flow within the chambers near the inside surface of the myocardium (endocardium) and the pericardial sac on the outside surface (epicardium), see figure (a)a. There are many applications for simulating heart blood flow doost16 . However, for many relevant questions the exact fluid dynamics of blood within the heart or a resolved fluid-solid interaction simulation are often not needed for simulating the myocardium. Instead, a realistic pressure-flow relationship stemming from the circulatory system is sufficient, which is commonly represented by lumped-parameter fluid models that provide the correct normal pressure to the endocardial wall kerckhoffs07 .

However, there is no consensus on boundary conditions to represent the effects of the pericardial sac. The goal of this work is twofold: (a) to provide a detailed literature review of pericardial biomechanics, hence justifying its modeling using a computationally inexpensive viscoelastic model, and (b) to highlight the relevance of such boundary conditions through a detailed quantitative analysis using a subject-specific cine MRI data set. We employ a four-chamber geometry including parts of the great vessels, as it provides us with additional options to asses the physiological correctness of our boundary condition, e.g. through the interplay between ventricles and atria during ventricular systole. Note, however, that the pericardial boundary condition is independent of the geometry and is meant to be applied to any kind of cardiac mechanics simulation that includes the epicardial surface.

This work is structured as follows. Following a review of the anatomy and physiology of the pericardium in section 1.2, we review pericadial boundary conditions currently used in cardiac mechanics simulations. We propose a model to represent the influence of the pericardium by parallel springs and dashpots acting in normal direction to the epicardium in section 2. Furthermore, we summarize a three-dimensional elastodynamical continuum model for the myocardium which is monolithically coupled to a zero-dimensional reduced-order windkessel model for the circulatory system. We demonstrate the influence of the pericardial boundary condition in section 3 through a detailed quantitative comparison of simulation results to cine MRI. For that we evaluate ventricular volume, atrioventricular-plane-displacement, atrioventricular interaction, and introduce a quantitative error measurement by calculating a distance error at endo- and epicardial surfaces between simulation results and cine MRI. We close this work with a discussion of the results, the limitations of our study, future perspectives, and some conclusions in section 4.

1.2 The pericardium

In the following, we review the anatomy of the pericardium and its physiology, where we focus on the mechanical interaction between the pericardium and the heart. Based on this review, we evaluate variants of pericardial boundary conditions and propose a model for pericardial-myocardial interaction.

1.2.1 Anatomy

As shown in figure (a)a, the pericardium is a sac-like structure with a combined thickness of 1-2 mm that contains the heart and parts of the great vessels holt70 . Figures (b)b and (c)c show a cross-sectional view of the myocardium and the layers of the pericardium. A common analogy for the location of the heart within the pericardium is that of a fist pushed into an inflated balloon martini15 .

The fibrous pericardium consists of a fibrous layer that forms a flask-like sac with a wavy collagenous structure of three interwoven main layers that are oriented to each other standring15 . It has a higher tensile stiffness than the myocardium and is dominated by the viscoelastic behavior of extracellular collagen matrix and elastin fibers lee85 . The fibrous pericardium is fixed in space by a ”three point cardiac seat belt” via the pericardial ligaments to the sternum. Furthermore, it is thoroughly attached to the central tendon of the thoracic diaphragm and additionally supported by the coats of the great vessels spodick96 . The various tissues, the fibrous pericardium is in contact with, can be seen in figure 2.

The fibrous pericardium contains a serous membrane, the serous pericardium, forming a closed sac. The serous pericardium is connected to the myocardium (visceral pericardium) and the fibrous pericardium (parietal pericardium). The composite of fibrous and parietal pericardium is commonly referred to as pericardium, whereas the visceral pericardium in contact with the myocardium is referred to as epicardium spodick96 . The space between the visceral and parietal pericardium is the pericardial cavity, which is filled by a thin film of fluid with an average volume of 20-25 ml holt70 . Beneath the visceral pericardium the heart is covered by a layer of adipose tissue, accumulated especially in the interventricular and atrioventricular grooves and around the coronary vessels, constituting about 20% of the heart weight rabkin07 .

(a) Coronal plane
(b) Transverse plane
(c) Sagittal plane
Figure 2: Position of the pericardium indicated in 3D MRI taken during diastasis. The neighboring tissue is color-coded: lungs (red), diaphragm (orange), sternum and ribs (light blue), aorta (dark blue), esophagus (purple), other (yellow). MRI courtesy of R. Chabiniok, J. Harmer, E. Sammut, King’s College London, UK.

1.2.2 Mechanical physiology

The pericardium serves multiple purposes spodick83 that can be grouped in: (a) membranous, it serves as a barrier against the spread of infection standring15 and (b) mechanical, it secures cardiac stability via its attachments within the thorax shabetai03 , as will be explained in the following. The mechanical properties of the pericardium itself can be found in sacks03 .

There is clear empirical evidence that the pericardium has a direct mechanical impact on the acute and chronic biomechanics of the heart. For example, in glantz78 it was discovered that the correlation of left and right ventricular pressure is higher with intact pericardium than after its complete removal. Maximal cardiac output during exercise can be increased acutely by the complete removal of the pericardium (pericardiectomy) through utilizing the Frank-Starling mechanism hammond92 . However, removing the pericardial restraint chronically promotes eccentric hypertrophy, i.e. an increase in dimension and mass of the heart and a change in shape from elliptical to spherical. The pericardium thus acts as a diastolic constraint for the heart by exercising a radial compression stress. This was confirmed in joebsis07 where it was observed that the opening angle of the myocardium with intact visceral pericardium is much higher than after its removal. The visceral pericardium is thus important for residual stress and passive stiffness.

It is widely accepted that the mechanism of the myocardium-pericardium interaction is through the pericardial fluid. In holt60 it was found that while increasing the volume of fluid within the pericardial cavity, the pericardial liquid pressure remains constant until it suddenly rises sharply. This led to the conclusion that most of the fluid is contained in pericardial sinuses and groves. This mostly empty space forms the so-called pericardial reserve volume, acting as a buffer against increasing pericardial liquid pressure. Only a small portion of the pericardial fluid remains as a thin film on the interface between parietal and visceral pericardium. In santamore90

, a dye was injected into the pericardial cavity near the apex. Fifteen minutes after injection the dye was almost exclusively found in the interventricular and atrioventricular grooves. This suggests that there is no significant fluid movement on the large surface areas of the ventricular free walls, leaving just a very thin film of fluid with an estimated thickness of less than 0.5 mm.

The mechanical constraint of the pericardium on diastolic cardiac function can be quantified by pericardial pressure. Here it is important to distinguish between liquid pressure and contact pressure smiseth85 ; tyberg86 . Liquid pressure describes the hydrostatic pressure inside the pericardial fluid and is measured by an open-ended catheter. However, liquid pressure does not describe the constraining effect of the pericardium on the myocardium. The constraint is assessed by contact pressure, which can be measured by a thin, flat, air-filled balloon catheter. In smiseth85 it was found that liquid pressure is substantially below contact pressure unless the pericardium contains a significant amount of pericardial fluid, which happens e.g. due to pericardial effusion. Furthermore, contact stress and thus ventricular restraint was maintained even though pericardial fluid was completely removed and liquid pressure at the epicardial surface was zero. Pericardial fluid therefore acts as lubrication rather than a load balancing mechanism, providing low-friction sliding between pericardium and epicardium hills85 .

There is less information available on the influence of the pericardium during systole. A pericardial restraining effect during systole would require a tension force to be transmitted by the myocardial-pericardial interface. The restraining effect of the pericardium during systole can be well observed in fish, where the parietal pericardium is almost rigid holt70 . It was observed in sudak65 that pericardial liquid pressure in smooth dogfish is always negative and decreases further during cardiac contraction. In man, sutton77 found that pericardial liquid pressure also drops during ventricular systole but remains positive throughout the cardiac cycle. However, to the best of the authors’ knowledge there is no study on the change of contact pressure during systole. It can be observed from mammal cine MRI that surrounding tissue moves toward the heart during systole, indicating attachment of pericardium and epicardium. We hypothesize that during systole, through the effect of adhesion, the pericardium remains in contact with the epicardium. This is analogous to the simple experiment of “gluing” two glass plates together with a drop of water. The glass plates can hardly be separated in normal direction but can be easily moved relatively to each other in tangential direction.

1.2.3 Current pericardial boundary conditions

For biventricular geometries, the constraining effect of the pericardium in diastole is accounted for in mansi10 ; chabiniok12 ; marchesseau13 , where a no penetration condition is enforced on the epicardium by a unidirectional penalty contact with a rigid pericardial reference surface. However, this neglects any constraining in systole by not allowing the pericardial interface to transmit any tension forces. Recently, it was proposed in santiago18a to completely prohibit movement normal to the epicardial surface, neglecting any elastic effects. The bi-directional elastic constraining effect of the pericardium is accounted for in spring-type boundary conditions, where a spring-dashpot boundary condition is enforced either on the base sainte-marie06 or on apex and valve annuli sermesant12 with homogeneous Neumann conditions applied to the rest of the epicardium. These boundary conditions are analogous to the external tissue support of the aorta in moireau12 ; moireau13 . However, they do not cover the whole epicardial surface thus representing pericardial-myocardial interaction only partially.

Fewer references exist for four chamber geometries. A common combination of boundary conditions for four chamber geometries are homogeneous Dirichlet conditions on vessel cut-offs and a soft material connected to the apex chabiniok12 ; augustin16 , or springs on the outside of great vessels land17 . In those cases, however, homogeneous Neumann conditions are applied on the remaining epicardial surface, neglecting any influence of the pericardium as in the biventricular case. In baillargeon14 ”omni-directional” springs acting in all directions are applied to the epicardium, artificially constraining any sliding movement along the pericardial-epicardial interface. To the authors’ best knowledge, the most detailed and physiologically correct representation of the pericardium so far was implemented in fritz13 . The pericardial-myocardial interaction was here modeled by a frictionless sliding, bi-directional penalty contact interaction in normal direction between the epicardium and a solid pericardial reference body. However, this condition is computationally very expensive as it requires solving an adhesial contact interaction problem. It also requires an additional solid body to be created, representing the surrounding tissue. Furthermore, no boundary conditions could be enforced at the great vessels since they were not included in the geometry. Thus, a fixation of the apex was necessary. All models based on four chamber geometries reviewed here lack a quantitative validation through comparison of simulation results to measurements, e.g. medical images like magnetic resonance imaging (MRI).

2 Models

In this work, we use a patient-specific four-chamber geometry from high-resolution static 3D MRI, including ventricles, atria, adipose tissue, and great vessels. Our cardiac model is formulated in a large displacement, constitutive nonlinear framework with nonlinear boundary conditions. It features high-resolution quadratic tetrahedral finite elements for structural dynamics with implicit time integration. Blood pressure is incorporated through monolithic coupling of the left and right ventricle to windkessel models which include each the atrioventricular and semilunar valves. The reference configuration is prestressed in all four cardiac chambers. The passive myocardial material features a state of the art orthotropic exponential material law proposed in holzapfel09 . Myofiber contraction in atria and ventricles is modeled with an active stress approach. Passive and active material behavior is based on local fiber orientations.

We follow the classic approach of nonlinear large deformation continuum mechanics to model the elastodynamic problem of 3D cardiac contraction. We define the reference configuration and current configuration which are connected by the displacements . We calculate the deformation gradient

, the Green-Lagrange strain tensor

, and the right Cauchy-Green tensor


2.1 Modeling the pericardium

Our aim in this work is to propose and justify a pericardial boundary condition that is both realistic and computationally inexpensive. The pericardial model we propose is based on rabkin75 and is sketched in figure 3. Using our code it was also already applied to a two-chamber geometry in hirschvogel16 . It consists of a spring and a dashpot in parallel acting in normal direction to the epicardial surface. Within the tangential plane we allow frictionless sliding to account for the lubricating effect of the pericardial fluid. A spring stiffness  and dashpot viscosity  contain the combined effects of serous pericardium, fibrous pericardium, and neighboring tissue. Generalizing the effect on the ventricles, spring compression models the pericadium’s constraining effect during passive ventricular filling, whereas spring expansion models the pericadium’s support during ventricular systole.

Figure 3: Serous pericardium, fibrous pericardium, and neighbouring tissue modeled by a spring (stiffness ) and a dashpot (viscosity ) in parallel. Spring and dashpot act in normal direction of the epicardial surface.

Note that only in the limit case of , we would obtain a boundary condition that penalizes and therefore prohibits any movement in normal direction to the epicardium, as it was recently proposed in santiago18a . However, our pericardial boundary condition is meant to be used with finite values for and , more realistically representing the visco-elastic support of the pericardium and its surrounding tissue and permitting movement normal to the epicardial surface. Furthermore, our parametric study in section 3.2 shows that small values of lead to physiological results.

In the following, we will derive a simple mathematical formulation for the pericardial boundary condition depicted in figure 3. This derivation will be carried out in two steps, where different assumptions are introduced in each step. Only the spring component will be considered during the derivation. However, all conclusions hold equivalently for the dashpot component.

Our goal is to preserve the features of the detailed pericardial boundary condition in fritz13 but arrive at a much simpler and cheaper formulation. As reviewed in section 1.2.3, pericardial-myocardial interaction is modeled in fritz13 by adhesial contact between the epicardium and an elastic reference body that is fixed in space and representing the surrounding tissue, see figure (a)a.

(a) Adhesial contact interaction from fritz13 between myocardium and fixed surrounding tissue (blue).
(b) Current normal spring from (3).
(c) Reference normal spring from (4).
Figure 4: Different formulations of interaction between myocardium (red) and pericardium.

In the first step, we replace the elastic body representing the surrounding tissue in fritz13 with springs acting in normal direction to the epicardium. Here, we assume that the elasticity of the surrounding tissue is linear with respect to the small movements of the epicardium in its normal direction. Note that we enforce the boundary conditions on the epicardial side of the myocardial-pericardial interface, as this does not require a representation of the actual pericardial surface. We therefore do not model the pericardium itself but the forces acting on the myocardium because of its presence. The elastic potential of a linear spring distributed on the epicardial surface in current configuration surface is


with spring stiffness , gap , and surface integral in current configuration . The calculation of the gap is illustrated in figure (b)b. We project a point on the current epicardial surface onto the point

on the reference epicardial surface. The distance between both points projected in the direction of the current outward normal vector

yields the gap function


Though reducing algorithmic and computational demands compared to contact interaction, this boundary condition still requires updates of the normal vector and its linearization with respect to the displacements as well as a projection of each evaluation point onto in each Newton iteration at each time step.

In a second step, we introduce two further simplifications. Instead of calculating the spring deformation from a projection, we directly use the spatial displacements . Furthermore, we use the epicardial normal vector in reference configuration (i.e. instead of ), neglecting any change in normal direction throughout the simulation. The formulation of the gap in (3) is then simplified to


The simplifications leading to (4) are valid for small rotations of the epicardium, an assumption that is not valid for all parts of the epicardium. However, the performance of both formulations (3) and (4) is reviewed in appendix 0.A.

We then arrive at the final expression for the pericardial boundary traction acting on the epicardial surface


For the sake of simplicity, we use here constant boundary condition parameters and on the whole epicardial surface. As it will be shown in the numerical examples, this simple approach already leads do greatly improved results. But of course, a regional distribution based on neighboring organs as visualized in figure 2 is also possible.

2.2 Geometrical model

To illustrate the effects of our pericardial model, we employ a four chamber geometry obtained in vivo from a 33 year old healthy female volunteer. The imaging data was acquired at King’s College London, UK using a Philips Achieva 1.5T magnetic resonance imaging (MRI) scanner with a dual-phase whole-heart 3D b-SSFP sequence uribe08 , acquisition matrix , acquired voxel size  mm, repetition time 4.5 ms, echo time 2.2 ms, echo train length 26 and flip angle . The diastolic rest period (diastasis) was used to generate the computational mesh. The geometry was meshed using Gmsh geuzaine09 with a resolution of , yielding nodes and quadratic tetrahedral elements, totalling a

structural degrees of freedom. Additionally, our geometry contains triangular surface elements with no additional degrees of freedom to track the movement of the planes of cardiac valves, allowing us to monitor the volumes of all four cardiac cavities.

All four cardiac cavities are closed with surface elements with no additional degrees of freedom at the valve planes depicted in red in figure (c)c at the left and right atrioventricular plane, respectively. The atria are additionally closed at their respective connections to the vasculature. We can thus monitor the volumes of all four cardiac cavities and track the movement of cardiac valve planes. The meshed geometry is shown in figure (a)a. The different materials are depicted in figure (c)c.

(a) Computational mesh.
(b) Atrial fibers and ventricular fibers at endocardium. Color shows helix angle with respect to long axis.
(c) Materials.
Figure 5: Four chamber patient-specific cardiac geometry.

Additionally, we prescribe the local angles of cardiac myofibers at epi- and endocardium of the ventricles. Using harmonic lifting, the fiber vectors

are interpolated to the interior of the domain by solving a Laplace problem

nagler16 . We interpolate the fiber orientation smoothly at each integration point. It is well known that the fiber orientation has a strong impact on passive and active cardiac mechanics ubbink06 ; wong10 ; gil13 ; eriksson13a ; asner16 ; nikou16 . In order to make a more clear statement about the pericardial boundary conditions independently of the fiber orientation and to show the interplay between boundary conditions and fiber orientations, we compare in this work three different fiber distributions: , , and . The first and second angle describe the fiber helix angle at the endo- and epicardial surface, respectively, with respect to the local circumferential axis. The transverse angle is zero throughout the myocardium. The sheet normal vector is perpendicular to the epi- and endocardial surfaces. The sheet vector is then obtained from . The atrial fiber architecture is obtained using a semi-automatic registration method based on the fiber definition in atlas atria hoermann17 ; hoermann18b . See figure (b)b for atrial fibers and ventricular fibers visualized at the endocardium.

2.3 Modeling cardiac contraction

Balance of momentum, a Neumann windkessel coupling condition with left ventricular pressure  acting on the left endocardium , omni-directional spring-dashpot boundary conditions, and pericardial boundary conditions yield the energy of the boundary value problem


with density , accelerations , the second Piola-Kirchhoff stress tensor , and spring stiffnesses and viscosities for vessel and apical surfaces and , respectively. Furthermore, we have virtual displacements and virtual strains . Omni-directional springs and dashpots are placed on the outsides of the great vessels and apical surface . The energy is complemented in section 3 by the energy of different boundary conditions.

We define different materials with the volumes defined as in figure 5 for adipose tissue (7), ventricular and atrial myocardium (8), and aorta and pulmonary artery (9):


Each material is composed of a hyperelastic and a viscous contribution depending on the rate of Green-Lagrange strains . The viscous behavior of myocardial tissue is modeled with a viscous pseudo-potential. Only the ventricular tissue in (8) has an additional active stress component . The strain energy density functions for exponential orthotropic solid holzapfel09 , Mooney-Rivlin solid sainte-marie06 ; chabiniok12 , Neo-Hookean solid , penalty function , and viscous pseudo-potential chapelle12 are given as


with the Jacobian of the deformation gradient and material parameters , , , and , penalty parameter , viscosity , and invariants

(a) Prescribed indicator functions for atria (blue) and ventricles (red).
(b) Active stress for atria (blue) and ventricles (red) with maximum values and , respectively.
Figure 6: Active stress.

We use the same active stress approach for both atrial and ventricular myocardium, though with different parameters. However, for simplicity of notation, we do not distinguish between atria or ventricles in the following equations. The active stress is computed as


with fiber stress in local reference fiber direction . Based on bestel01 , we model fiber stress by the evolution equation


with activation function

, contractility , and the function . The activation function is modeled by


with maximum and minimum activation rates and , respectively, and functions


with steepness

 s and descending and ascending sigmoid functions

and , respectively. The indicator function indicates systole. The times and model the onset of systole and diastole, respectively. Note that the active stress tensor is the only input of our solid model we prescribe over time. Our structural model can be coupled to a model of electric signal propagation as shown in hoermann18 . However, as the focus in this work is on pericardial boundary conditions, a coupled electro-mechanical model would only introduce unnecessary complexity and variability. Using the parameters in table (a)a, we obtain the active stress curve depicted in figure (b)b. The values and denote the maximum of the active stress for atria and ventricles, respectively.

Using the finite element method, we discretize displacements and virtual displacements arising in the weak form (6) by


with quadratic basis functions and nodal displacements on each tetrahedral element . Assembly of the discretized problem leads to the matrix notation


with mass matrix , force vector , and discrete displacements, velocities, and accelerations , and , respectively. We discretize the boundary value problem in time with Newmark’s method newmark59


with parameters and , and time step size . Additionally, we apply the generalized- method chung93 , yielding quantities at a generalized time step


depending on the weights and for force vector and mass matrix respectively. Newmark’s method in combination with the generalized -scheme is a common technique for implicit one-step time integration for finite elements in nonlinear solid dynamics. Finally, we obtain the time and space discrete structural residual


All parameters used for the elastodynamical model are summarized in table (a)a.

(a) Schematic.
(b) Prescribed left (blue) and right (red) atrial pressure .
Figure 7: Windkessel model.

2.4 Modeling the circulatory system

We use the same windkessel model for each ventricle with different parameters. For simplicity of notation, we drop in this section the index of the ventricle in all windkessel parameters and variables. Note that this cardiovascular model does not represent a closed-loop system since the total blood volume is not conserved, i.e. blood exiting the right ventricle into the lungs does not enter the left atrium. However, using a windkessel model for each ventricle provides us with a reasonable approximation of ventricular pressures.

In this work we use a four element windkessel model based on the ideas in westerhof08 and kerckhoffs07 . A comprehensive review of different windkessel models is given in shi11 . The schematic of our windkessel model is given in figure (a)a using resistances , compliances , and an inertance . Pressures at different parts of the model are denoted by . We distinguish between a proximal (index ) and a distal part (index ) of the outlets, i.e. lung and aorta for the right and left ventricle, respectively. The atrial pressure is prescribed to simulate atrial systole, see figure (b)b. The reference pressure  is kept constant.

We model the atrioventricular and semilunar valves with a smooth diode-like behavior by non-linear resistances and , respectively, depending on the sigmoid function


with the steepness of the sigmoid function and the minimal and maximal valve resistance and , respectively. This yields the set of differential equations


The 0D windkessel model is strongly coupled to the 3D structural model. The 0D model depends on the structural displacements of the 3D model via the change in ventricular volume . On the other hand, the 3D model depends on left and right ventricular pressure from the 0D model. The coupling between both models is described in section 2.6. The vector of primary variables yields , including the flux through the inertance . We discretize the set of windkessel equations (23) in time with the one-step- scheme


This yields the discrete windkessel residual evaluated at time step . The parameters and initial conditions of the cardiovascular model are summarized in table (b)b and table (c)c respectively. Windkessel parameters are motivated by values from literature and adapted for this heart to yield physiological pressures as well as approximately a periodic state of the windkessel systems.

Name Par. Value Unit
All tissues
Tissue density
Viscosity 0.1
Volumetric penalty
Active myocardial tissue
Atrial contractility 9.72 kPa
Ventricular contractility see table (b)b
Activation rate
Deactivation rate
Atrial systole 70
Atrial diastole 140
Ventricular systole see table (b)b
Ventricular diastole 484
Passive myocardial tissue (holzapfel09 table 1, shear, figure 7)
Matrix 0.059
Fiber 18.472
Sheet 2.481
Fiber-sheet 0.216
Great vessels
Mooney-Rivlin 5.0
Mooney-Rivlin 0.04
Spring stiffness
Dashpot viscosity
Adipose tissue
Neo-Hooke 1.0
Pericardial boundary condition
see table (a)a
(a) Parameters of the elastodynamical model.
Name Par. Value Unit Proximal inertance Proximal capacity Distal capacity Proximal resistance Distal resistance Reference pressure 0 Closed valve resistance Open valve resistance Valve steepness (b) Parameters of the reduced order cardiovascular model (identical for left and right ventricle). Name Par. Value Unit Left Right Atrial pressure 6.0 4.0 Ventricular pressure 8.0 6.0 Proximal pressure 61.8 24.0 Distal pressure 59.7 23.2 Proximal flow 38.3 14.9 (c) Initial conditions of the reduced order cardiovascular model. Par. Value Generalized- , , 0.5 0.25 One-step- 1.0 (d) Numerical time integration parameters.
Table 1: Overview of parameters in four-chamber cardiac model.

2.5 Prestress

For our reference configuration we use a patient-specific geometry segmented from static 3D MRI at diastolic rest period (diastasis), see section 2.2. Diastasis is very suitable for the reference configuration, since both ventricular and atrial myofibers are relaxed, the heart is not accelerated, and blood pressures are minimal and constant. This simplifies the task of obtaining the stress state of the reference configuration, which in this case is determined by the static blood pressures within the cardiac cavities. We therefore have to prestress our geometry with the initial ventricular pressure from table (c)c. In this work, we use the Modified Updated Lagrangian Formulation as proposed in gee09a ; gee10 . This method incrementally calculates a deformation gradient with respect to an unknown stress-free reference configuration. From this deformation gradient, a stress field is calculated so that the segmented geometry of the heart is in balance with the prestressed pressure state. This yields a prestress within the myocardium as well as in the pericardial boundary condition in case pericardium. Note that while this technique allows to model prestress, we do not account for the residual stress inherent in myocardial tissue joebsis07 .

2.6 Solving the 0D-3D coupled problem

We solve the coupled 0D-3D model with the structural and windkessel residuals and , respectively, at time step with the Newton-Raphson method


in a monolithic fashion for increments in displacements and windkessel variables and , respectively, at iteration until convergence. We build and solve this coupled system using our in-house code BACI wall10 . The numerical parameters for the time integration of our model are listed in table (d)d.

3 Results

(a) Case apex with omni-directional spring-dashpots on (green).
(b) Case pericardium with normal spring-dashpots on (red).
Figure 8: Surface definitions for boundary conditions with omni-directional spring-dashpots on (blue) and homogeneous Neumann boundary conditions (white).

In order to investigate the influence of pericardial boundary conditions, we compare simulations with and without pericardial boundary condition on the epicardial surface . The simulations will be denoted by apex and pericardium in the following. See table (a)a for an overview of used parameters.

Case apex depicted in figure (a)a yields the boundary-vale problem


adding the energy for omni-directional spring dashpots to the energy (6), where is the apical surface. The apical surface is defined as the epicardial surface within 10 mm of the apex, see figure 8. It resembles homogeneous Neumann boundary conditions on , i.e. the absence of any pericardial boundary conditions as frequently found in literature sermesant12 ; chabiniok12 ; augustin16 .

Case pericardium depicted in figure (b)b yields the boundary-vale problem


adding the energy for the pericardial boundary condition to the energy (6), where is the pericardial surface. The choice of parameters and is detailed in appendix 3.2.

The remainder of this section is structured as follows. We first give an overview of all methods used in this work to quantify the difference of simulation and MRI in section 3.1. Next, we calibrate model parameters for both cases in section 3.2. In the following sections, we validate various outputs of both simulation cases apex and pericardium with measurements from cine MRI. We begin with scalar windkessel outputs and compare the simulation volume curves to MRI in section 3.4. A qualitative evaluation of displacement results is given in section 3.5 by comparing end-systolic simulation results to cine MRI frames at multiple views. We quantify the differences in pumping motion for simulation cases apex and pericardium by comparing the displacements of the left and right atrioventricular plane to MRI in section 3.6. The interplay between ventricles and atria with and without the presence of pericardial boundary conditions is investigated in section 3.7. In section 3.8 we calculate a spatial error for the left and right endocardium to quantify the overall approximation quality. Finally, we evaluate the contact stress of our pericardial boundary condition of case pericardium in section 3.9.

3.1 Assessment of cardiac function

In this section, we briefly describe the various methods we use throughout this work to quantify cardiac function of different simulations.

Cine MRI

We use cine MRI with a temporal resolution of 30 ms in four- (figure (a)a), three-, and two-chamber views and short axis planes with a slice distance of 8 mm (figures (c)c, (d)d, (b)b). All cine MRI data used in this work is rigidly registered to the static 3D image taken during diastasis and used for geometry creation to account for any movement of the subject during image acquisition.

It is important to note that the reference configuration of our simulation is obtained from static 3D imaging with a fine isotropic resolution and acquired in free-breathing, as explained in section 2.2. For the comparison of simulation results to cine MRI however, we have to rely on sparsely distributed images acquired in expiration breath-hold. The used image types rely on different MRI acquisition parameters and pulse sequences. Therefore, it is impossible for our simulation to match the cine MRI data perfectly, even in reference configuration. This error however is usually smaller than the approximation error of the cardiac model.

(a) Four chamber view.
(b) Short axis endocardial contours used for error calculation.
(c) Short axis slice 9.
(d) Short axis slice 6.
Figure 9: Post processing planes for simulations and cine MRI.
Left ventricular volume

We obtain a reference for left ventricular volume by manually segmenting the left endocardial surface obtained from the short axis cine MRI stack at all time steps. We add the sum of areas in each short axis slice multiplied by the slice thickness. We cut the volume at the top and bottom according to the limits of the left ventricle at each time step, as observed in two chamber and four chamber views.

In order to be fair and not introduce a bias towards our more realistic model, for each simulation, we calibrate the contractility . It is a key parameter describing cardiac elastodynamics, resembling the asymptotical active fiber stress in (13). It controls maximum deformation during systole. In order to make simulations comparable, we adapt for each combination of boundary condition and fiber orientation to match the left ventricular volume at end-systole as segmented from cine MRI of . The heart thus yields a stroke volume of 75 ml and an ejection fraction (EF) of 57%.

Atrioventricular plane displacement (AVPD)

The movement of the left or right plane of the valve separating atrium and ventricle in long axis direction during the cardiac cycle is described by AVPD. For left and right ventricle those valves are termed mitral and tricuspital valve, respectively. As a scalar parameter, AVPD at end-systole is an important clinical parameter to describe and predict cardiac vitality carlsson07 ; willenheimer97 .

We evaluate AVPD in this work as it gives us a quantitative measurement of the displacements in long axis direction. We semi-automatically extract the displacement of the left and right atrioventricular plane from two, three, and four chamber cine MRI using the freely available software Segment version 2.0 R5585 heiberg10 . In our simulations, we average the displacements on all nodes on the valve plane (see the red planes in figure (c)c) and project them onto the long axis direction. A positive sign indicates a movement of the base towards the apex.

Spatial error

We validate displacement in long axis direction using AVPD as measurement. To validate radial displacement we compare the movements of cardiac surfaces in simulations to the ones from short axis cine MRI. For comparison, we select the left and right endocardium, as it shows how pericardial boundary conditions prescribed on the epicardium act on the interior of the domain.

For each MRI time step (temporal resolution 30 ms) we select the closest simulation time step (temporal resolution 1 ms). Spatially, we extract the simulations’ displacement results at the same positions where the cine MRI slices were acquired. This is possible since we use the MRI scanner’s global coordinate system for all images and the simulation. This method can be thought of as taking a virtual cine MRI of the simulation. This yields an Eulerian description of motion, as the observer is fixed in space. The difference of simulated displacements to cine MRI data was used previously, e.g. in chabiniok12 to estimate local tissue contractility. Note that this technique does not allow us to track rotations of the left ventricle due to its rotational symmetry.

We manually segment the contours of left and right endocardium from short axis cine MRI for slices 5 to 9 at all MRI time steps, see figure (b)b. These slices are selected because the myocardium is recognizable for all MRI time steps and not disturbed by either apex or AVP. The function converts both MRI and simulated endocardial contours and , respectively, to binary images with a resolution of for every slice . We use the Dice metric to compare both binary images


where denotes the area of the binary image.

Ventricular-atrial interaction

Utilizing a four chamber geometry allows us to investigate the interaction between ventricles and atria. Specifically, we want to study the influence of ventricular contraction on atrial filling. We therefore analyze atrial volumes over time. Furthermore, we segmented left and right atrial volumes at ventricular diastasis and end-systole from isotropic 3D MRI.

Pericardial contact stress

We evaluate the stresses transmitted between the epicardial boundary conditions and the myocardium for both cases apex and pericardium. We use different averaged stresses for both cases to quantify the constraining effect of each boundary condition. In case apex the stresses are concentrated on the small apical area and acting in any direction. We thus integrate the stress vectors of the apical boundary condition over the apical surface and normalize by the apical area to obtain the mean apical stress


In case pericardium the boundary stresses are distributed over the whole epicardial surface and acting only in normal direction. Therefore, we extract the (signed) normal component and integrate it over the epicardial surface to obtain the mean pericardial stress


normalized by the epicardial area.

3.2 Selection of pericardial paramaters

Since in case apex the purpose of the apical boundary condition is fixing the apex throughout cardiac contraction, we chose a high spring stiffness permitting only little motion. For case pericardium, the parameters and describing pericardial stiffness and viscosity, respectively, need to be calibrated. The chosen value for pericardial viscosity has on its own, i.e. without parallel spring, only little influence on cardiac dynamics. However, in combination with the spring, it prevents unphysiological oscillations of the heart. Pericardial stiffness controls the amount of displacement perpendicular to the epicardial surface and thus the radial motion of the myocardium.

We investigate in the following the influence of the parameter on the contraction of the heart. For this study, we limit ourselves to the fiber distribution, as it is commonly used in cardiac simulations, see e.g. chabiniok12 ; lee15a ; hirschvogel16 . We tested the following parameter values:


For each choice of , we calibrated active stress to yield the same end-systolic volume as measured from MRI. All parameters except are kept constant throughout this study. Specifically, we did not adjust the timing of ventricular systole to match the volume curve from MRI. However, since all simulations reach the end-systolic state it will be used in this section for quantitative comparisons.

The results of the calibration are shown in figure 10, where maximum active stress is plotted against pericardial stiffness. For comparison, the result for case apex with fibers is included. Active stress required to yield identical end-systolic volume rises strongly with increasing pericardial stiffness. The temporal maximum of pericardial contact stress averaged over the epicardium also increases strongly with , as shown in figure 11. For high , contact stress has the same order of magnitude as active stress and exceeds maximum left ventricular pressure. For small , it has the same order of magnitude as atrial pressure, which experimentally shown to be a good predictor for pericardial contact stress tyberg86 .

Figure 12 shows the volume within the pericardial cavity, calculated as the combined volume of all tissue inside the pericardium and the volume within the four cardiac cavities. Case pericardium yields a lower volume change than case apex and decreases further with increasing .

The end-systolic state of the simulations is shown in figure 13 compared to MRI. The images contain all simulated variants for , where the color changes continuously from  kPa/mm (blue) to  kPa/mm (red). All MRI views in figure 13 show clearly that pericardial stiffness controls radial displacement of the epicardium. High stiffness values result in less radial inward motion during ventricular systole than visible in cine MRI and vice versa. This is also well observable in figure (a)a for the atria in four-chamber view. The short axis views in figures (b)b and (c)c additionally show that the interventricular septum is stretched and rotated as compared to MRI for high .

The spatial error at left and right ventricular endocardium is shown in figures (a)a and (b)b, respectively. The increasing mismatch between simulations and MRI for increasing as visible in figures (b)b and (c)c is quantified as increasing spatial error.

Left and right AVPD is displayed in figures (c)c and (d)d, respectively. In the left ventricle, i.e. at the mitral valve, AVPD is not very sensitive to the choice of . However, it is higher than in case apex but much lower than in MRI. For the right ventricle, i.e. at the tricuspid valve, AVPD is greatly enhanced by increasing towards the value measured in MRI. An identical trend is observable for left and right atrial volume in figures (e)e and (f)f, respectively.

Figure 10: Maximum ventricular active stress calibrated to yield identical end-systolic volume. Shown for case pericardium with varying pericardial stiffness compared to case apex, both with fiber distributions.
Figure 11: Maximum of mean pericardial contact stress for case pericardium with varying pericardial stiffness and fiber distributions.
Figure 12: Volume change for case pericardium with varying pericardial stiffness compared to case apex, both with fiber distributions.
(a) Four chamber view
(b) Short axis view slice 9
(c) Short axis view slice 6
Figure 13: Cine MRI at end-systole for case pericardium with fiber distribution from (blue) to (red). Views as defined in figure 9. MRI courtesy of R. Chabiniok, J. Harmer, E. Sammut, King’s College London, UK.
(a) Left ventricular endocardial error.
(b) Right ventricular endocardial error.
(c) Left atrioventricular plane displacement.
(d) Right atrioventricular plane displacement.
(e) Left atrial volume.
(f) Right atrial volume.
Left ventricular endocardial error.
Figure 14: Kinematic scalar cardiac quantities at end-systole for case pericardium with varying pericardial stiffness compared to case apex both with fiber distributions and MRI.

To conclude the parametric study for pericardial stiffness, we choose for all following simulations the value  kPa/mm. It offers a low spatial error at the ventricles but has higher atrial volume and AVPD than the simulation with  kPa/mm.

3.3 Model personalization

All simulations are carried out in the following using three different fiber distributions, i.e. , , and . The results for the calibration of are shown in table (b)b. Note that we show here the maximum value of active stress instead of . It can be observed that is larger in case pericardium than in case apex. Furthermore, increases from to fibers for more vertical fiber distributions.

The onset of systole and diastole, and , as well as the myofiber activation and deactivation rates, and , are adapted to the left ventricular volume curve for ventricles and atria. Here, parameters for the atria are fitted from the interval and for the ventricles. The material parameter controlling the viscosity of the tissue is fitted during ventricular diastole, i.e. . Since active stress is zero during this interval, viscosity controls the relaxation speed of the model. A summary of all calibrated model-specific material parameters is given in table 2. For parameters identical in all models see table 1.

apex 0 0
pericardium 0 0 0.2
(a) Spring stiffness and dashpot viscosity on apical and epicardial surface.
[kPa] [ms]
apex 63.5 72.4 91.0 143 155 172
pericardium 79.4 90.7 129 161 170 193
(b) Maximal myocardial active stress and ventricular activation time .
Table 2: Calibrated parameters for simulation cases apex and pericardium and different fiber orientations.

3.4 Scalar windkessel results

Firstly, in figure 15 we compare the scalar outputs volume (left) and pressure (right) of the left ventricle of our windkessel model. As explained in section 3.2, the contractility was calibrated in all simulations to match end-systolic volume as segmented from cine MRI. Therefore, in figures (a)a and (c)c the volumes of MRI and all simulations match at  s. Furthermore, although they result from simulations with very different boundary conditions and fiber orientations the volume curves are very similar. The maximum volume due to atrial contraction and the prescribed atrial pressure in figure (b)b is similar in both cases but lower than in MRI. As for the volume curves, the pressure curves in figures (b)b and (d)d are remarkably similar despite the different simulation settings. Because case pericardium exhibits a faster decay in volume during systole than in case apex, the pressure peak during systole is more pronounced.

(a) Case apex: LV volume
(b) Case apex: LV pressure
(c) Case pericardium: LV volume
(d) Case pericardium: LV pressure
Figure 15: Simulation results for volume (left) and pressure (right) of the left ventricle (LV) for boundary condition cases apex (top) and pericardium (bottom). Volume results are compared to cine MRI.

3.5 Displacements at end-systole

As demonstrated in section 3.4, the results of the scalar output parameters left ventricular volume and pressure are mostly invariant to changes in boundary conditions or fiber orientation. Validating the elastodynamical model of cardiac contraction thus requires a comparison of displacement results to spatially distributed MRI observations, see figure 16. The reference configuration (diastasis) of the simulation is shown in figures (a)a, (b)b, (c)c. We compare the MRI frames at end-systole to our simulation results using the four chamber view, see figures (d)d and (g)g, and two different short axis views, see figures (e)e, (h)h, (f)f, (i)i. The location of the view planes is visualized in figures (a)a, (c)c, (d)d.

For case apex, there is a radial inward movement of the myocardial wall. In figure (d)d, this is especially visible at right atrial free wall and at the left and right epicardial free wall. There is a large mismatch between simulation and MRI at the interatrial septum. Due to the radial contraction motion, atrioventricular plane displacement (AVPD) is lower than in MRI. The fixation of the apex in case apex causes a mismatch between simulations and MRI at the apex, as the apex slightly moves during cardiac contraction. The interventricular septum’s matches well with MRI in figures (e)e and (f)f. However, the mismatch of epicardial contours is clearly visible and sensitive to fiber orientation.

Comparing figures (d)d and (g)g, the influence of the pericardial boundary condition becomes clearly visible. It can be observed for case pericardium in figure (g)g that the epicardial contour matches the MRI much closer than case apex in figure (d)d for any fiber orientation. The movements of the left and right atrioventricular plane also match well with MRI, for both orientation and displacement in normal direction. The displacements at the apical region are also predicted more accurately than in case apex. Comparing the shape of the right ventricle in figures (d)d and (g)g, one can observe that the pumping motion of the right ventricle in case apex is the result of radial movement, whereas in case pericardium it is the result of a downward movement of the atrioventricular plane. The same observation holds for a less visible degree for the left ventricle. Through the constraining effect of the pericardium, the atria are visibly more stretched than in case apex. There is also an influence of the fiber orientation in case pericardium, although it is more bound to the endocardial surfaces. The more vertical the fiber orientation, i.e. from (red) to (blue), the larger the displacements of the atrioventricular planes and the smaller the displacement of the apex in anterior direction. There are some mismatches between simulation and MRI at the interatrial and interventricular septum but less pronounced than in case apex. The deviation at the interventricular septum can be observed for short axis slice 9 in figure (h)h. For short axis slice 6 in figure (i)i there is a good agreement with simulation and MRI at all regions of the left and right myocardium.

(a) Reference configuration (diastasis) four chamber view
(b) Reference configuration (diastasis) short axis view slice 9
(c) Reference configuration (diastasis) short axis view slice 6
(d) Case apex four chamber view
(e) Case apex short axis view slice 9
(f) Case apex short axis view slice 6
(g) Case pericardium four chamber view
(h) Case pericardium short axis view slice 9
(i) Case pericardium short axis view slice 6
 four chamber view
Figure 16: Reference configuration (diastasis) as well as simulation results and cine MRI at end-systole in four chamber view and short axis views as defined in figure 9. MRI courtesy of R. Chabiniok, J. Harmer, E. Sammut, King’s College London, UK

3.6 Atrioventricular plane displacement

The AVPDs of simulations and MRI are compared in figure 17. The left and right AVPD from MRI (black) is zero at the beginning and at the end of the cardiac cycle. During atrial systole, the left and right atriovencular planes (AVP) move away from the apex and reach both their minimal value at atrial end-systole at . Followed by ventricular systole, the AVPs move towards the apex and both reach their maximal value at ventricular end-systole at .

During atrial systole for , negative AVPD, i.e. movement of the AVP towards the atria, is less pronounced and delayed in both cases as compared to MRI. However, extremal AVPD at atrial end-systole is slightly higher in case pericardium than in case apex.

Comparing AVPD cases apex and pericardium in figures (a)a, (c)c, (b)b, (d)d, one can observe that in both cases maximum AVPD depends on fiber orientation: Maximum AVPD increases from horizontal fibers (red) to vertical fibers (blue). The dependence on fiber orientation is more pronounced in case apex than in case pericardium. In general, AVPD is slightly higher in case pericardium than in case apex but still underestimates measurements from MRI.

(a) Case apex: left AVPD
(b) Case apex: right AVPD
(c) Case pericardium: left AVPD
(d) Case pericardium: right AVPD
Figure 17: Simulated atrioventricular plane displacement for left and right ventricle compared to cine MRI.

3.7 Ventricular-atrial interaction

Atrial systole is visible by the drop in atrial volume in both cases. Passive atrial filling is non-existent in case apex, as the volume in figures (a)a and (b)b stay constant during ventricular systole. This is also visible at the end-systolic four-chamber view in figure (d)d. For fibers, the right atrium is even slightly emptied during ventricular systole, as observed in figure (b)b. Atrial filling can be observed for case pericardium in figures (c)c and (d)d. Both atria are visibly filled during ventricular systole, although maximum atrial volume remains smaller than in MRI.

(a) Case apex: Left atrium
(b) Case apex: Right atrium
(c) Case pericardium: Left atrium
(d) Case pericardium: Right atrium
Figure 18: Simulated volume curves for left and right atrium compared to 3D MRI at ventricular diastasis and end-diastole.

3.8 Spatial error

For case apex in figures (a)a and (b)b the error is lowest in both ventricular endocardia during contraction at end-systole at . The error rises during ventricular contraction and relaxation. Errors at the end of the simulation higher than the ones at suggest that the state at the end of the simulation differs from the reference configuration. The overall error is much lower in case pericardium than in case apex.

(a) Case apex: left endocardium spatial error
(b) Case apex: right endocardium spatial error
(c) Case pericardium: left endocardium spatial error
(d) Case pericardium: right endocardium spatial error
Figure 19: Relative spatial error of simulation results and cine MRI at left and right endocardium.

3.9 Boundary stresses

Both scalar boundary stresses and are visualized in figure 20 over time for all fiber orientations. It can be observed that apical stress in case apex is orders of magnitude higher than pericardial stress in case pericardium and more dependent on fiber orientation. Positive values of indicate predominant tensile stresses between epicardium and pericardium. It can be seen that mean pericardial stress in figure (b)b is a compressive stress for most of the cardiac cycle, except at the end of systole and onset of diastole.

(a) Mean apical contact stress for case apex.
(b) Mean pericardial contact stress for case pericardium.
Mean apical contact stress
Figure 20: Boundary condition stress.

Boundary stresses are visualized in figure 21. For case apex, the mean stress vectors for all three fiber distributions are shown in figure (a)a at and scaled according to their magnitude. Fiber orientation has not only a strong influence on the magnitude but also the direction of mean apical stress.

The local distribution of pericardial contact stress with fibers at end-systole is shown in figure (b)b in reference configuration. At end-systole, compressive as well as tensile stresses occur. Stresses are centered around a tensile stress of 20 mmHg. Areas of high compressive stresses are at the left atrium, the anterior and posterior right ventricle, the posterior left ventricle, and the anterior left ventricular apex. Areas of high tensile stresses are the right ventricle close to the anterior part of the right ventricular outflow tract and the left and right ventricular free wall. Overall, pericardial contact stress is evenly distributed around the epicardial surface.

(a) Case apex mean apical stress vectors (scaled by magnitude) for (red), (green), and (blue) fiber orientations at .
(b) Case pericardium pericardial contact stress on epicardial surface with fibers at end-systole .
Figure 21: Visualization of boundary stresses.

4 Discussion

Our objective was to analyze the effects of the pericardial boundary condition proposed in section 2.1 based on the physiology of the pericardium, comparing simulation cases pericardium and apex. We first performed a parametric study to explore the influence of pericardial stiffness. Each simulation case was then personalized and evaluated for the fiber orientations , , and . We then compared scalar left ventricular pressure and volume. The displacements at end-systole were qualitatively compared to multi-view cine MRI. Additionally, we quantified the differences of both simulation cases to MRI by atrioventricular plane displacement (AVPD), passive atrial filling, and spatial approximation error at the left and right ventricular endocardium.

4.1 Pericardial stiffness

The parametric study for pericardial stiffness in case pericardium in section 3.2 revealed that the ventricles are well approximated by the lowest tested stiffness values, e.g.  kPa/mm. Here, the error at left and right ventricular endocardium was minimized and much lower than in case apex.

In contrast, right AVPD and right atrial passive filling matched well with measurements from MRI for high stiffness values, e.g.  kPa/mm. Choosing this value globally for pericardial stiffness lead however to some undesirable consequences, namely unphysiologically high myocardial contractility and pericardial stress as well as bad approximation of the interventricular septum.

In future studies, it might thus be reasonable to select spatially varying pericardial parameters. This hypothesis is supported by the fact that the pericardial tissue is in contact with various organs of different material properties as outlined in section 1.2.2. A starting point could be the estimation of regional pericardial parameters based on the surface definitions in figure 2 with the objective to match MRI measurements in section 3.2.

In case of a biventricular geometry, no atria are present. Thus AVPD is not controlled by the interaction of atria and pericardium. Furthermore, atrial filling is not taken into account. We thus expect that a global value of  kPa/mm for pericardial stiffness yields good results for a biventricular geometry with fibers. This value was also used in hirschvogel16 , although it was not really analyzed there, e.g. with respect to MRI.

4.2 Pumping mechanism

We calibrated cardiac contractility in all simulations in section 3.2 to yield the same end-systolic volume. It was shown that in case pericardium, higher contractilities are required than in case apex. Therefore, for a given contractility, a heart constrained with the pericardial boundary condition yields less output. This result is in agreement with the experimental observation that cardiac output is greatly increased after the removal of the pericardium hammond92 . The result further agrees with the numerical experiments performed in fritz13 . For identical active stress, left ventricular ejection fraction decreased from 71 % to 63 % when including the pericardium.

The main pumping mechanism of the heart is shortening in long axis direction, which is quantified by AVPD arutunyan15 ; arvidsson15 . In maksuti15 , the pumping function of the heart was compared to a piston unit with the AVP as a piston. This mechanism could be observed in section 3.6 for case pericardium, where left and right AVPD is higher than in case apex but still lower than in MRI.

The upper part of the left atrium is fixated by pulmonary veins. Ventricular contraction forces the mitral ring towards the apex and promotes the filling of the left atrium from the pulmonary veins fujii94 . In section 3.7 we compared atrial filling during ventricular systole with and without pericardium. It was observed that the simulations of case pericardium which promoted higher AVPD in section 3.6 contribute more to atrial filling during ventricular systole. Case pericardium predicts maximal atrial volume at ventricular end-systole as segmented from isotropic 3D MRI better than case apex. The simulated values are however still lower than in MRI for the chosen pericardial parameters.

Keeping in mind that all simulations yield the same end-systolic volume, it was shown that the pumping mechanism of the heart is very different for cases apex and pericardium although their pressure and volume curves were similar in section 3.4. Comparison of four-chamber and short axis slices of the left and right ventricle from simulation results to cine MRI in section 3.5 revealed an unphysiological radial pumping motion without pericardial boundary conditions in case apex. In emilsson01 it was found that the outer diameter of the left ventricle shortens only about 2 mm during systole. Furthermore, the total volume enclosed by the pericardium changes only by about 5-8 % during the cardiac cycle bowman03 ; arvidsson15 . We found that for the fiber orientation the total change in pericardial volume is 24 % and 21 % for cases pericardium and apex, respectively. This mismatch is mainly due to the unphysiological change in atrial volume during ventricular contraction.

As demonstrated in the parametric study in section 3.2, AVPD and atrial filling could be increased to the values measured in MRI by increasing the global pericardial stiffness. However, this was shown to lead to a worse approximation of the interventricular septum. This motivates the use of a regionally distributed pericardial stiffness. Another reason for underestimating AVPD and atrial filling might be an atrial material model, which is in our case identical to the ventricular one, that is too stiff.

4.3 Fiber orientation

We compared in this work three fiber orientations within the myocardium, namely , , and . We studied the influence of fiber direction for both boundary conditions cases apex and pericardium. It was shown that fiber orientation has a strong influence on the displacements. In section 3.6 fiber orientations exhibited larger AVPD for both boundary condition cases. This can be attributed to the fact that the fiber orientation is more vertical, i.e. more aligned with the long axis, than and fiber orientations. Since myofiber contraction is prescribed in fiber direction, more vertical fiber orientations inherently apply a greater force pushing apex and AVP together, thus yielding higher AVPD. Since the AVP is also attached to the atria, AVPD is also linked to atrial filling. In section 3.7 it was shown that the more vertical fiber orientation also yielded the highest atrial filling during ventricular systole. Comparing results to short axis cine MRI slices, it was shown that a more horizontal fiber orientation leads to a more radial contraction of the heart. The maximum pericardial stress at end-systole was highest for fibers. This can be explained by the observation in figure (h)h where fibers (red) exhibited the most radial inward movement during systole. Since the myocardial-pericardial interface can only transmit forces in normal direction, a more radial contraction exerts a higher pericardial tensile stress. The overall spatial approximation error was also shown to be dependent on fiber direction. However, the dependence was more pronounced in case apex than in case pericardium.

4.4 Pericardial contact stress

In smiseth85 end-diastolic pericardial contact pressure was measured with a flat balloon catheter at the left ventricular anterolateral epicardial surface with around 15 mmHg. In vivo experiments on humans in tyberg86 showed pericardial pressures on the left lateral surface of the heart between 0 and 15 mmHg. The prestressing procedure in our model not only includes the myocardium but also the pericardial boundary condition in case pericardium. Here, we measure a contact pressure of 20 mmHg at diastasis on the left ventricular epicardial surface, agreeing well with experimental observations.

The stresses exerted by the boundary conditions on the epicardial surface of the heart were found to be over one order of magnitude higher in case apex than in case pericardium. The exact stress values in case apex depend on the choice of apical spring stiffness, which was not calibrated in this study. It is nevertheless evident, that unphysiologically high stresses are concentrated in a very small area of the heart. In case pericardium, all boundary stresses are evenly distributed on the epicardial surface.

The similar maximum values of mean pericardial contact stress for all fiber directions in case pericardium suggest that pericardial constraint is displacement-controlled. Pericardial constraint is determined by the deviation of the heart throughout the cardiac cycle from its end-diastolic state. However, as outlined in section 1.2.2, to the best of the authors’ knowledge there are no measurements of pericadial contact pressure during the cardiac cycle to validate the stresses experienced in our computational study. Pericardial contact stress is thus an output of our computational model which is not yet available in clinical practice. Given that our model was solely calibrated to kinematic data, the pericardial contact stresses predicted by our model should be considered as qualitative results.

4.5 Numerical performance

We ran all simulations on two nodes of our Linux cluster. One node features 64 GB of RAM and two Intel Xeon E5-2680 ”Haswell” processors, each equipped with 12 cores operating at a frequency of 2.5 GHz. The computation time of cases apex and pericardium was almost identical, which was about 18 hours for each of the simulations performed in this work, including prestressing. The pericardial boundary condition requires little effort to evaluate, since the pericardial boundary condition only requires the displacement field, which is computed anyway, and reference surface normals, which are computed once at the initialization of the simulation. Some differences in numerical performance arise since the calculated displacement fields of both cases are different. Taking the fiber distribution, case apex had an average of 7.8 Newton iterations per time step and 28 linear solver iterations per Newton iteration. For case pericardium, these values were 8.3 and 25, respectively.

4.6 Limitations and future perspectives

As mentioned earlier, in this work, we did not account for the propagation of the electrical signal sent from the sinus node. Rather, all myocardial tissue in our simulations was activated simultaneously. We recently demonstrated the ability to couple our mechanical model to an electrophysiological model hoermann18 , which we can include in further studies. However, since the data came from a healthy volunteer, we do no expect relevant variations.

Ex-vivo experiments on myocardial tissue in yin87 ; dokos02 ; sommer15 showed anisotropic tissue characteristics, depending on myocardial fiber and sheet orientation. In our model, we used the anisotropic material model proposed in holzapfel09 for myocardial tissue. Due to the lack of sufficient experimental data, we used identical material properties for left and right myocardium, as well as the atria. However, no studies have been carried out how material parameters obtained from experiments on ex vivo tissue correlate to in vivo material behavior. Furthermore, it should be noted that vastly different material parameters have been estimated in holzapfel09 and gueltekin16 when being fitted to measurements from either biaxial extension tests or shear tests.

Our structural model was coupled to a lumped-parameter windkessel model of hemodynamics of the systemic and pulmonary circulation with prescribed atrial pressures. The interaction between atria and ventricles should be investigated in further studies using a volume-preserving closed-loop model, including both pulmonary and systemic circulation. Furthermore, none of our cardiac simulations are perfectly periodic, i.e. the values at the end of the cardiac cycle are not equal to the initial conditions. In future studies, achieving a periodic state should be incorporated into parameter estimation.

In this work we interpolated the local helix fiber directions at the integration points from three different prescribed constant-per-surface fiber orientations. Results showed that fiber orientation has a large influence on AVPD. However, we have no knowledge of patient-specific fiber orientation and assumed equal distributions in left and right ventricle. Patient-specific cardiac fiber orientations can be estimated from diffusion tensor MRI nagler17 (DTMRI). However, while applicable to in vivo DTMRI (as shown in nagler17 ), to the best of our knowledge, fiber estimation has not been tested and validated with in vivo DTMRI yet. Further quantitative studies of cardiac dynamics require a fine resolution of patient-specific fibers.

We further assumed constant stiffness and viscosity parameters of our pericardial boundary condition over the epicardial surface. Given reliable material parameters for the myocardium, constant pericardial stiffness and viscosity could be estimated from measured AVPD. The choice of constant parameters might however be oversimplified, as the pericardium is in contact with various tissues of different mechanical behaviors, as illustrated in figure 2. For example, the movement of the apex in anterior direction in case pericardium as observed in figure (g)g suggests a higher pericardial stiffness to model the influence of the sternum and the diaphragm. This will however introduce more parameters to the model, which will need to be calibrated to measurements from e.g. cine or 3D tagged MRI. For this study we kept the number of parameters small in order to make evident the general effect of the pericardial boundary condition even by using a simplified modeling approach.

From a machine learning perspective, we split our limited available data from cine MRI into a training set and a test set. The training set data is used during model personalization. The rest of the data can then be used in the test set to check how well the model actually predicts data that was not used during personalization. In our case, we used as training set left ventricular volume and ventricular epicardial contours to tune timing, (de-) activation rates, and contractility for atria and ventricles and global material viscosity and pericardial stiffness. We then used as test set AVPD, atrial volume, and ventricular endocardial contours, each left and right, to quantify the simulations’ approximation error. Many more parameters of our cardiac model could be personalized for this patient-specific study. However, using the metrics in our test set for model calibration would disqualify using them to test model accuracy and limit our abilities to test the model.

We validated our simulation results solely with cine MRI data. Cine MRI can be interpreted as an Eulerian description of cardiac movement, as the imaging planes stay fixed in space throughout the cardiac cycle. This observation however cannot detect any rotational movement with respect to the long axis, as the left ventricle is almost rotationally symmetric. To properly validate any rotational movement of the myocardium, a comparison to data from 3D tagged MRI is necessary, which can be interpreted as a Lagrangian observation of cardiac motion. Furthermore, pressure measurements from within ventricles and atria are required. Pressure values at end-diastole yield initial values for the stress state of the myocardium, which cannot be assessed from imaging alone. Pressure curves over the cardiac cycle would yield a ground truth to validate the outputs of our windkessel model. Figure 18 demonstrates that the model, while using the pericardial constraint, does predict accurately the atrial volume at ventricular end-systole. However, we have no data available at atrial end-systole. In future studies, if detailed cine data of atria are available (e.g. cine stack in trasverse orientation with respect to the body, and using thin slices of 5mm), we will consider a more detailed analysis of atrial contraction.

4.7 Concluding remarks

In this work we gave an overview of the anatomy and mechanical function of the pericardium and motivated to model its influence on the myocardium as a parallel spring and dashpot acting on the epicardial surface. Following a review of pericardial boundary conditions currently used in mechanical simulations of the heart, we proposed to compare two simulation cases, one with and one without pericardial boundary conditions. Following calibration to stroke volume as measured from short axis cine MRI, we compared several physiological key outputs of our model and validated them using multi-view cine MRI. Although exhibiting similar volume and pressure curves, the displacement results of both simulation cases were radically different. The simulations with pericardial boundary conditions matched MRI measurements much closer than without, especially with respect to atrioventricular plane displacement and atrial filling during ventricular systole, quantities which were not included in the calibration of the model. By establishing an overall spatial approximation error at the left and right endocardium, we showed that the introduction of only two global parameters for the pericardial boundary condition already yields a big gain in model accuracy. Our ultimate goal is to obtain more comprehensive data sets, adding 3D tagged MRI and pressure measurements, to further validate our model of pericardial-myocardial interaction.

Measurements of pericardial contact stress at different locations on the epicardium throughout the cardiac cycle would help to test the qualitative predictions of pericardial contact stresses by our model and will probably lead to further model improvements.

Compliance with ethical standards

Conflict of interest The authors declare that they have no conflict of interest.

(a) Ellipsoid model in reference configuration with cross-section (yellow).
(b) Frontal view of epi- and endocardial cross-section at end-systole.
(c) Top-down view of epicardial cross-section at end-systole.
Figure 22: Ellipsoid model in reference configuration (black) with cases free (green), pericardium (blue), and pseudo-contact (orange).

Appendix 0.A Comparison of spring formulations

We show in section 2.1 how the pericardial boundary condition in case pericardium can be derived from adhesive sliding contact by introducing several simplifications. To justify the simplifications made by our pericardial boundary condition, we use a very simple geometry of a hollow half-ellipsoid with fibers, which roughly represents the shape of the left ventricle, see figure (a)a. It is able to show the consequences of each approach while being simple enough to isolate the effects of the boundary condition. The parameters of the ellipsoid model are given in table 3. We use the same active stress model introduced in (8) to mimic cardiac contraction. All three simulations use the same contractility parameter.

As in section 3, case pericardium utilizes the pericardial boundary condition proposed in (5) using the gap (4). Additionally, we introduce case pseudo-contact, which uses the definition of the gap in (3) based on projection and the current normal vector to the epicardium. Case free has homogeneous zero-Neumann boundary conditions on the whole epicardial surface.

The results of the contraction simulation are shown in figures (b)b and (c)c at end-systole. Displayed is the reference configuration and all three boundary condition cases for a cross-section of the ellipsoid. Figure (b)b shows in a frontal view the shortening of the ellipsoid with visible epi- and endocardial contours. While cases pericardium and pseudo-contact are very similar with little differences only in radial direction, case free exhibits much less longitudinal shortening. There is almost no longitudinal shortening but a translational movement of the whole geometry instead.

Figure (b)b shows the epicardial contour of the ellipsoid in a top-down view to observe the twisting motion of the ellipsoid. All three boundary condition cases are very similar. This confirms that the normal springs in cases pericardium and pseudo-contact in fact allow tangential sliding and do not prohibit any rotational movement, as they are very similar to case free. Furthermore, the similarity of cases pericardium and pseudo-contact shows that the simplified spring formulation (4) in case pericardium is sufficient to represent the effects of the pericardium compared to the more detailed formulation (3) in case pseudo-contact.

Name Par. Value Unit
Tissue density
Viscosity 10
Volumetric penalty
Ventricular contractility 185
Mooney-Rivlin 10
Mooney-Rivlin 40
Spring stiffness base 1
Pericardial spring stiffness
Case free 0
Case pericardium 20
Case pseudo-contact 20
(a) Parameters of the elastodynamical model.
Table 3: Overview of material parameters in ellipsoid model. For numerical parameters see table (d)d.


  • (1) Iaizzo P. A. Handbook of cardiac anatomy, physiology, and devices. Springer Science & Business Media; 2015.
  • (2) Sainte-Marie J, Chapelle D, Cimrman R, Sorine M. Modeling and estimation of the cardiac electromechanical activity. Computers & Structures. 2006;84(28):1743–1759.
  • (3) Doost S. N, Ghista D, Su B, Zhong L, Morsi Y. S. Heart blood flow simulation: a perspective review. BioMedical Engineering OnLine. 2016;15(1):101.
  • (4) Kerckhoffs R. C, et al. Coupling of a 3D finite element model of cardiac ventricular mechanics to lumped systems models of the systemic and pulmonic circulation. Annals of Biomedical Engineering. 2007;35(1):1–18.
  • (5) Holt J. P. The normal pericardium. The American Journal of Cardiology. 1970;26(5):455–465.
  • (6) Martini F. H, Timmons M. J, Tallitsch R. B. Human Anatomy. Pearson Education; 2015.
  • (7) Standring S. Gray’s Anatomy: The Anatomical Basis of Clinical Practice. Elsevier Health Sciences; 2015.
  • (8) Lee J. M, Boughner D. R. Mechanical properties of human pericardium. Differences in viscoelastic response when compared with canine pericardium.. Circulation Research. 1985;57(3):475–481.
  • (9) Spodick D. H. The pericardium: a comprehensive textbook. Informa Health Care; 1996.
  • (10) Rabkin S. W. Epicardial fat: properties, function and relationship to obesity. Obesity Reviews. 2007;8(3):253–261.
  • (11) Spodick D. H. The normal and diseased pericardium: Current concepts of pericardial physiology, diagnosis and treatment. Journal of the American College of Cardiology. 1983;1(1):240-251.
  • (12) Shabetai R. The Pericardium. Springer; 2003.
  • (13) Sacks M. S. Incorporation of experimentally-derived fiber orientation into a structural constitutive model for planar collagenous tissues. Journal of Biomechanical Engineering. 2003;125(2):280–287.
  • (14) Glantz S. A, et al. The pericardium substantially affects the left ventricular diastolic pressure-volume relationship in the dog. Circulation Research. 1978;42(3):433-41.
  • (15) Hammond H. K, White F. C, Bhargava V, Shabetai R. Heart size and maximal cardiac output are limited by the pericardium. American Journal of Physiology - Heart and Circulatory Physiology. 1992;263(6):H1675–H1681.
  • (16) Jöbsis P. D, et al. The visceral pericardium: macromolecular structure and contribution to passive mechanical properties of the left ventricle. American Journal of Physiology - Heart and Circulatory Physiology. 2007;293(6):H3379–H3387.
  • (17) Holt J. P, Rhode E. A, Kines H, Ruth H. Pericardial and Ventricular Pressure. Circulation Research. 1960;8(6):1171-1181.
  • (18) Santamore W. P, Constantinescu M. S, Bogen D, Johnston W. E. Nonuniform distribution of normal pericardial fluid. Basic Research in Cardiology. 1990;85(6):541–549.
  • (19) Smiseth O. A, Frais M. A, Kingma I, Smith E. R, Tyberg J. V. Assessment of pericardial constraint in dogs. Circulation. 1985;71(1):158-64.
  • (20) Tyberg J. V, et al. The relationship between pericardial pressure and right atrial pressure: an intraoperative study. Circulation. 1986;73(3):428-32.
  • (21) Hills B. A, Butler B. D. Phospholipids identified on the pericardium and their ability to impart boundary lubrication. Annals of Biomedical Engineering. 1985;13(6):573–586.
  • (22) Sudak F. Intrapericardial and intracardiac pressures and the events of the cardiac cycle in Mustelus canis (Mitchill). Comparative Biochemistry and Physiology. 1965;14(4):689 - 705.
  • (23) Sutton J, Gibson D. G. Measurement of postoperative pericardial pressure in man. British Heart Journal. 1977;39(1):1-6.
  • (24) Mansi T. Image-based physiological and statistical models of the heart: application to tetralogy of Fallot. Phd-thesis, École Nationale Supérieure des Mines de Paris 2010.
  • (25) Chabiniok R, et al. Estimation of tissue contractility from cardiac cine-MRI using a biomechanical heart model. Biomechanics and Modeling in Mechanobiology. 2012;11(5):609-630.
  • (26) Marchesseau S, Delingette H, Sermesant M, Ayache N. Fast parameter calibration of a cardiac electromechanical model from medical images based on the unscented transform. Biomechanics and Modeling in Mechanobiology. 2013;12(4):815–831.
  • (27) Santiago A, et al. Fully coupled fluid-electro-mechanical model of the human heart for supercomputers. International Journal for Numerical Methods in Biomedical Engineering. 2018;.
  • (28) Sermesant M, et al. Patient-specific electromechanical models of the heart for the prediction of pacing acute effects in CRT: A preliminary clinical validation. Medical Image Analysis. 2012;16(1):201 - 215.
  • (29) Moireau P, et al. External tissue support and fluid–structure simulation in blood flows. Biomechanics and Modeling in Mechanobiology. 2012;11(1-2):1–18.
  • (30) Moireau P, et al. Sequential identification of boundary support parameters in a fluid-structure vascular model using patient image data. Biomechanics and Modeling in Mechanobiology. 2013;12(3):475-496.
  • (31) Augustin C. M, et al. Anatomically accurate high resolution modeling of human whole heart electromechanics: A strongly scalable algebraic multigrid solver method for nonlinear deformation. Journal of Computational Physics. 2016;305:622 - 646.
  • (32) Land S, Niederer S. A. Influence of atrial contraction dynamics on cardiac function. International Journal for Numerical Methods in Biomedical Engineering. 2017;DOI:10.1002/cnm.2931.
  • (33) Baillargeon B, Rebelo N, Fox D. D, Taylor R. L, Kuhl E. The Living Heart Project: A robust and integrative simulator for human heart function. European Journal of Mechanics - A/Solids. 2014;48(0):38 - 47.
  • (34) Fritz T, Wieners C, Seemann G, Steen H, Dössel O. Simulation of the contraction of the ventricles in a human heart model including atria and pericardium. Biomechanics and Modeling in Mechanobiology. 2013;13(3):1–15.
  • (35) Holzapfel G. A, Ogden R. W. Constitutive modelling of passive myocardium: a structurally based framework for material characterization. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences. 2009;367(1902):3445–3475.
  • (36) Rabkin S, Hsu P. Mathematical and mechanical modeling of stress-strain relationship of pericardium. American Journal of Physiology – Legacy Content. 1975;229(4):896–900.
  • (37) Hirschvogel M, Bassilious M, Jagschies L, Wildhirt S, Gee M. A monolithic 3D-0D coupled closed-loop model of the heart and the vascular system: Experiment-based parameter estimation for patient-specific cardiac mechanics. International Journal for Numerical Methods in Biomedical Engineering. 2016;33(8):e2842.
  • (38) Uribe S, et al. Volumetric Cardiac Quantification by Using 3D Dual-Phase Whole-Heart MR Imaging. Radiology. 2008;248(2):606–614.
  • (39) Geuzaine C, Remacle J.-F. Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities. International Journal for Numerical Methods in Engineering. 2009;79(11):1309–1331.
  • (40) Nagler A, Bertoglio C, Ortiz M, Wall W. A. A spatially varying mathematical representation of the biventricular cardiac fiber architecture. technical report: Institute for Computational Mechanics, Technische Universität München ; Center for Mathematical Modeling, Universidad de Chile ; 2016.
  • (41) Ubbink S, Bovendeerd P, Delhaas T, Arts T, Vosse F. Towards model-based analysis of cardiac MR tagging data: Relation between left ventricular shear strain and myofiber orientation. Medical Image Analysis. 2006;10(4):632 - 641. Special Issue on Functional Imaging and Modelling of the Heart (FIMH 2005).
  • (42) Wong K. C. L, et al. Cardiac Motion Estimation Using a ProActive Deformable Model: Evaluation and Sensitivity Analysis. In: Springer Berlin Heidelberg; 2010; Berlin, Heidelberg.
  • (43) Gil D, et al. What a Difference in Biomechanics Cardiac Fiber Makes. In: Camara O, Mansi T, Pop M, Rhode K, Sermesant M, Young A, eds. Statistical Atlases and Computational Models of the Heart. Imaging and Modelling Challenges, Lecture Notes in Computer Science, vol. 7746: Springer Berlin Heidelberg 2013 (pp. 253-260).
  • (44) Eriksson T, Prassl A, Plank G, Holzapfel G. Influence of myocardial fiber/sheet orientations on left ventricular mechanical contraction. Mathematics and Mechanics of Solids. 2013;18(6):592-606.
  • (45) Asner L, et al. Estimation of passive and active properties in the human heart using 3D tagged MRI. Biomechanics and Modeling in Mechanobiology. 2016;15(5):1121–1139.
  • (46) Nikou A, Gorman R. C, Wenk J. F. Sensitivity of left ventricular mechanics to myofiber architecture: A finite element study. Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine. 2016;230(6):594-598. PMID: 26975892.
  • (47) Hörmann J. M, et al. Multiphysics Modeling of the Atrial Systole under Standard Ablation Strategies. Cardiovascular Engineering and Technology. 2017;8(2):205–218.
  • (48) Hörmann J. M, Pfaller M. R, Bertoglio C, Avena L, Wall W. A. Automatic Mapping of Atrial Fiber Orientations for Patient-Specific Modeling of Cardiac Electromechanics using Image-Registration. submitted to International Journal for Numerical Methods in Biomedical Engineering. 2018;.
  • (49) Chapelle D, Le Tallec P, Moireau P, Sorine M. An energy-preserving muscle tissue model: formulation and compatible discretizations. International Journal for Multiscale Computational Engineering. 2012;10(2):189–211.
  • (50) Bestel J, Clément F, Sorine M. A Biomechanical Model of Muscle Contraction. In: Niessen W. J, Viergever M. A, eds. Medical Image Computing and Computer-Assisted Intervention – MICCAI 2001, Springer Berlin Heidelberg; 2001; Berlin, Heidelberg.
  • (51) Hörmann J. M, et al. An adaptive Hybridizable Discontinuous Galerkin approach for cardiac electrophysiology. International Journal for Numerical Methods in Biomedical Engineering. 2018;34(5):e2959.
  • (52) Newmark N. M. A method of computation for structural dynamics. Journal of the Engineering Mechanics Division. 1959;85(3):67–94.
  • (53) Chung J, Hulbert G. A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized- method. Journal of Applied Mechanics. 1993;60(2):371–375.
  • (54) Westerhof N, Lankhaar J.-W, Westerhof B. E. The arterial windkessel. Medical & Biological Engineering & Computing. 2008;47(2):131–141.
  • (55) Shi Y, Lawford P, Hose R. Review of zero-D and 1-D models of blood flow in the cardiovascular system. Biomed Eng Online. 2011;10:33.
  • (56) Gee M. W, Reeps C, Eckstein H. H, Wall W. A. Prestressing in finite deformation abdominal aortic aneurysm simulation. Journal of Biomechanics. 2009;42(11):1732 - 1739.
  • (57) Gee M. W, Förster C, Wall W. A. A computational strategy for prestressing patient-specific biomechanical problems under finite deformation. International Journal for Numerical Methods in Biomedical Engineering. 2010;26(1):52–72.
  • (58) Wall W. A, et al . Baci: a parallel multiphysics simulation environment. : Technical report, Institute for Computational Mechanics, Technische Universität München; 2018.
  • (59) Carlsson M, Ugander M, Mosén H, Buhre T, Arheden H. Atrioventricular plane displacement is the major contributor to left ventricular pumping in healthy adults, athletes, and patients with dilated cardiomyopathy. American Journal of Physiology - Heart and Circulatory Physiology. 2007;292(3):H1452–H1459.
  • (60) Willenheimer R, Cline C, Erhardt L, Israelsson B. Left ventricular atrioventricular plane displacement: an echocardiographic technique for rapid assessment of prognosis in heart failure. Heart. 1997;78(3):230–236.
  • (61) Heiberg E, et al. Design and validation of Segment - freely available software for cardiovascular image analysis. BMC Medical Imaging. 2010;10(1):1.
  • (62) Lee L. C, Sundnes J, Genet M, Wenk J. F, Wall S. T. An integrated electromechanical-growth heart model for simulating cardiac therapies. Biomechanics and Modeling in Mechanobiology. 2016;15(4):791–803.
  • (63) Arutunyan A. H. Atrioventricular plane displacement is the sole mechanism of atrial and ventricular refill. American Journal of Physiology-Heart and Circulatory Physiology. 2015;308(11):H1317-H1320. PMID: 25795710.
  • (64) Arvidsson P. M, Carlsson M, Kovács S. J, Arheden H. Letter to the Editor: Atrioventricular plane displacement is not the sole mechanism of atrial and ventricular refill. American Journal of Physiology-Heart and Circulatory Physiology. 2015;309(6):H1094-H1096. PMID: 26374902.
  • (65) Maksuti E, Bjällmark A, Broomé M. Modelling the heart with the atrioventricular plane as a piston unit. Medical Engineering & Physics. 2015;37(1):87 - 92.
  • (66) Fujii K, et al. Effect of left ventricular contractile performance on passive left atrial filling - Clinical study using radionuclide angiography. Clinical Cardiology. 1994;17(5):258–262.
  • (67) Emilsson K, Brudin L, Wandt B. The mode of left ventricular pumping: is there an outer contour change in addition to the atrioventricular plane displacement?. Clinical Physiology. 2001;21(4):437-446.
  • (68) Bowman A. W, Kovács S. J. Assessment and consequences of the constant-volume attribute of the four-chambered heart. American Journal of Physiology-Heart and Circulatory Physiology. 2003;285(5):H2027-H2033. PMID: 12869381.
  • (69) Yin F. C, Strumpf R. K, Chew P. H, Zeger S. L. Quantification of the mechanical properties of noncontracting canine myocardium under simultaneous biaxial loading. Journal of Biomechanics. 1987;20(6):577 - 589.
  • (70) Dokos S, Smaill B. H, Young A. A, LeGrice I. J. Shear properties of passive ventricular myocardium. American Journal of Physiology - Heart and Circulatory Physiology. 2002;283(6):H2650–H2659.
  • (71) Sommer G, et al. Biomechanical properties and microstructure of human ventricular myocardium. Acta Biomaterialia. 2015;24:172–192.
  • (72) Gültekin O, Sommer G, Holzapfel G. A. An orthotropic viscoelastic model for the passive myocardium: continuum basis and numerical treatment. Computer Methods in Biomechanics and Biomedical Engineering. 2016;0(0):1-18. PMID: 27146848.
  • (73) Nagler A, Bertoglio C, Stoeck C. T, Kozerke S, Wall W. A. Maximum likelihood estimation of cardiac fiber bundle orientation from arbitrarily spaced diffusion weighted images. Medical Image Analysis. 2017;39:56 - 77.