Martin R Pfaller
email: martin.pfaller@tum.de

Julia Hoermann

Martina Weigl

Andreas Nagler

Radomir Chabiniok

Cristóbal Bertoglio

Wolfgang A Wall

Julia M Hörmann

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

Keywords: Cardiac mechanical modeling, Pericardium, Boundary conditions, Finite element simulation

come

Introduction

Motivation

Cardiac mechanics simulations consist of solving a nonlinear elastodynamic boundary value problem [START_REF] Sainte-Marie | Modeling and estimation of the cardiac electromechanical activity[END_REF]. 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 1a. There are many applications for simulating heart blood flow [START_REF] Doost | Heart blood flow simulation: a perspective review[END_REF]. However, for many relevant questions the exact fluid dynamics of blood within the heart or a resolved fluidsolid interaction simulation are often not needed for simulating the myocardium. Instead, a realistic pressureflow relationship stemming from the circulatory system is sufficient, which is commonly represented by lumpedparameter fluid models that provide the correct normal pressure to the endocardial wall [START_REF] Kerckhoffs | Coupling of a 3D finite element model of cardiac ventricular mechanics to lumped sys-tems models of the systemic and pulmonic circulation[END_REF].

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 reducedorder 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-planedisplacement, 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.

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.

Anatomy

As shown in figure 1a, 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 [START_REF] Holt | The normal pericardium[END_REF]. Figures 1b and1c 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 [START_REF] Martini | Human Anatomy[END_REF].

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 120 • to each other [START_REF] Standring | s Anatomy: The Anatomical Basis of Clinical Practice[END_REF]. It has a higher tensile stiffness than the myocardium and is dominated by the viscoelastic behavior of extracellular collagen matrix and elastin fibers [START_REF] Boughner | Mechanical properties of human pericardium. Differences in viscoelastic response when compared with canine pericardium[END_REF]. 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 [START_REF] Spodick | The pericardium: a comprehensive textbook[END_REF]. 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 [START_REF] Spodick | The pericardium: a comprehensive textbook[END_REF]. 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 [START_REF] Holt | The normal pericardium[END_REF]. 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 [START_REF] Rabkin | Epicardial fat: properties, function and relationship to obesity[END_REF].

Mechanical physiology

The pericardium serves multiple purposes [START_REF] Spodick | The normal and diseased pericardium: Current concepts of pericardial physiology, diagnosis and treatment[END_REF] that can be grouped in: (a) membranous, it serves as a barrier against the spread of infection [START_REF] Standring | s Anatomy: The Anatomical Basis of Clinical Practice[END_REF] and (b) mechanical, it secures cardiac stability via its attachments within the thorax [START_REF] Shabetai | The Pericardium[END_REF], as will be explained in the following. The mechanical properties of the pericardium itself can be found in [START_REF] Sacks | Incorporation of experimentally-derived fiber orientation into a structural constitutive model for planar collagenous tissues[END_REF].

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 [START_REF] Glantz | The pericardium substantially affects the left ventricular diastolic pressure-volume relationship in the dog[END_REF] 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 [START_REF] Hammond | Heart size and maximal cardiac output are limited by the pericardium[END_REF]. 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 [START_REF] Jöbsis | The visceral pericardium: macromolecular structure and contribution to passive mechanical properties of the left ventricle[END_REF] 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 [START_REF] Holt | Pericardial and Ventricular Pressure[END_REF] it was found that while increas-ing 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 [START_REF] Santamore | Nonuniform distribution of normal pericardial fluid[END_REF], 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 [START_REF] Smiseth | Assessment of pericardial constraint in dogs[END_REF][START_REF] Tyberg | The relationship between pericardial pressure and right atrial pressure: an intraoperative study[END_REF]. 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 [START_REF] Smiseth | Assessment of pericardial constraint in dogs[END_REF] 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 [START_REF] Hills | Phospholipids identified the pericardium and their ability to impart boundary lubrication[END_REF].

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 [START_REF] Holt | The normal pericardium[END_REF]. It was observed in [START_REF] Sudak | Intrapericardial and intracardiac pressures and the events of the cardiac cycle in Mustelus canis (Mitchill)[END_REF] that pericardial liquid pressure in smooth dogfish is always negative and decreases further during cardiac contraction. In man, [START_REF] Sutton | Measurement of postoperative pericardial pressure in man[END_REF] 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, indi- cating 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.

Current pericardial boundary conditions

For biventricular geometries, the constraining effect of the pericardium in diastole is accounted for in [START_REF] Mansi | Image-based physiological and statistical models of the heart: application to tetralogy of Fallot[END_REF][START_REF] Chabiniok | Estimation of tissue contractility from cardiac cine-MRI using a biomechanical heart model[END_REF][START_REF] Marchesseau | Fast parameter calibration of a cardiac electromechanical model from medical images based on the unscented transform[END_REF], 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 [START_REF] Santiago | Fully coupled fluid-electro-mechanical model of the human heart for supercomputers[END_REF] 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 springtype boundary conditions, where a spring-dashpot boundary condition is enforced either on the base [START_REF] Sainte-Marie | Modeling and estimation of the cardiac electromechanical activity[END_REF] or on apex and valve annuli [START_REF] Sermesant | Patient-specific electromechanical models of the heart for the prediction of pacing acute effects in CRT: A preliminary clinical validation[END_REF] 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 [START_REF] Moireau | External tissue support and fluidstructure simulation in blood flows[END_REF][START_REF] Moireau | Sequential identification of boundary support parameters in a fluid-structure vascular model using patient image data[END_REF]. 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 [START_REF] Chabiniok | Estimation of tissue contractility from cardiac cine-MRI using a biomechanical heart model[END_REF][START_REF]Anatomically accurate high resolution modeling of human whole heart electromechanics: A strongly scalable algebraic multigrid solver method for nonlinear deformation[END_REF], or springs on the outside of great vessels [START_REF] Land | Influence of atrial contraction dynamics on cardiac function[END_REF]. In those cases, however, homogeneous Neumann conditions are applied on the remaining epicar-dial surface, neglecting any influence of the pericardium as in the biventricular case. In [START_REF] Baillargeon | The Living Heart Project: A robust and integrative simulator for human heart function[END_REF] "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 [START_REF] Fritz | Simulation of the contraction of the ventricles in a human heart model including atria and pericardium. Biomechanics and Modeling in Mechanobiology[END_REF]. The pericardial-myocardial interaction was here modeled by a frictionless sliding, bidirectional 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).

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 con- figuration is prestressed in all four cardiac chambers. The passive myocardial material features a state of the art orthotropic exponential material law proposed in [START_REF] Holzapfel | Constitutive modelling of passive myocardium: a structurally based framework for material characterization[END_REF]. 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 X and current configuration x which are connected by the displacements u = x -X. We calculate the deformation gradient F , the Green-Lagrange strain tensor E, and the right Cauchy-Green tensor

C F = ∂x ∂X , E = 1 2 (C -I), C = F T F . (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 [START_REF] Rabkin | Mathematical and mechanical modeling of stress-strain relationship of pericardium[END_REF] and is sketched in figure 3. Using our code it was also already applied to a two-chamber geometry in [START_REF] Hirschvogel | A monolithic 3D-0D coupled closed-loop model of the heart and the vascular system: Experiment-based parameter estimation for patient-specific cardiac mechanics[END_REF]. 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 k and dashpot viscosity c 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. Note that only in the limit case of k → ∞, 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 [START_REF] Santiago | Fully coupled fluid-electro-mechanical model of the human heart for supercomputers[END_REF]. However, our pericardial boundary condition is meant to be used with finite values for k and v, 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 k 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 [START_REF] Fritz | Simulation of the contraction of the ventricles in a human heart model including atria and pericardium. Biomechanics and Modeling in Mechanobiology[END_REF] but arrive at a much simpler and cheaper formulation. As reviewed in section 1.2.3, pericardial-myocardial interaction is modeled in [START_REF] Fritz | Simulation of the contraction of the ventricles in a human heart model including atria and pericardium. Biomechanics and Modeling in Mechanobiology[END_REF] by adhesial contact between the epicardium and an elastic reference body that is fixed in space and representing the surrounding tissue, see figure 4a.

In the first step, we replace the elastic body representing the surrounding tissue in [START_REF] Fritz | Simulation of the contraction of the ventricles in a human heart model including atria and pericardium. Biomechanics and Modeling in Mechanobiology[END_REF] 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 Γ epi in current configuration surface is

W = 1 2 Γ epi kg 2 da (2)
with spring stiffness k, gap g, and surface integral in current configuration da. The calculation of the gap is illustrated in figure 4b. We project a point x ∈ Γ epi on the current epicardial surface onto the point X proj ∈ Γ epi 0 on the reference epicardial surface. The distance between both points projected in the direction of the current outward normal vector n yields the gap function

g = (x -X proj ) • n. (3)
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 Γ epi 0 in each Newton iteration at each time step. In a second step, we introduce two further simplifications. Instead of calculating the spring deformation (a) Adhesial contact interaction from [START_REF] Fritz | Simulation of the contraction of the ventricles in a human heart model including atria and pericardium. Biomechanics and Modeling in Mechanobiology[END_REF] between myocardium and fixed surrounding tissue (blue). Fig. 4: Different formulations of interaction between myocardium (red) and pericardium.

from a projection, we directly use the spatial displacements u. Furthermore, we use the epicardial normal vector in reference configuration (i.e. N instead of n), neglecting any change in normal direction throughout the simulation. The formulation of the gap in (3) is then simplified to

g = u • N (4)
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 A.

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

t epi = N (k p u • N + c p u • N ) . (5)
For the sake of simplicity, we use here constant boundary condition parameters k p and c p 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.

Geometrical model

To Additionally, we prescribe the local angles of cardiac myofibers at epi-and endocardium of the ventricles. Using harmonic lifting, the fiber vectors f 0 are interpolated to the interior of the domain by solving a Laplace problem [START_REF] Nagler | A spatially varying mathematical representation of the biventricular cardiac fiber architecture[END_REF]. 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 [START_REF] Ubbink | Towards model-based analysis of cardiac MR tagging data: Relation between left ventricular shear strain and myofiber orientation[END_REF][START_REF] Wong | Cardiac Motion Estimation Using a ProActive Deformable Model: Evaluation and Sensitivity Analysis[END_REF][START_REF] Gil | What a Difference in Biomechanics Cardiac Fiber Makes[END_REF][START_REF] Eriksson | Influence of myocardial fiber/sheet orientations on left ventricular mechanical contraction[END_REF][START_REF] Asner | Estimation of passive and active properties in the human heart using 3D tagged MRI[END_REF][START_REF] Nikou | Sensitivity of left ventricular mechanics to myofiber architecture: A finite element study[END_REF]. 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: ±50 • , ±60 • , and ±70 • . 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 n 0 is perpendicular to the epi-and endocardial surfaces. The sheet vector s 0 is then obtained from s 0 = n 0 × f 0 . The atrial fiber architecture is obtained using a semi-automatic registration method based on the fiber definition in atlas atria [START_REF] Hörmann | Multiphysics Modeling of the Atrial Systole under Standard Ablation Strategies[END_REF][START_REF] Hörmann | Automatic Mapping of Atrial Fiber Orientations for Patient-Specific Modeling of Cardiac Electromechanics using Image-Registration[END_REF]. See figure 5b for atrial fibers and ventricular ±60 • fibers visualized at the endocardium. (c) Materials.

Fig. 5: Four chamber patient-specific cardiac geometry.

Modeling cardiac contraction

Balance of momentum, a Neumann windkessel coupling condition with left ventricular pressure p v acting on the left endocardium Γ endo 0 , omni-directional springdashpot boundary conditions, and pericardial boundary conditions yield the energy δΠ of the boundary value problem

δΠ = Ω0 ρ 0 ü • δu dV + Ω0 S : δE dV + ν∈{l,r} Γ endo,ν 0 p ν v F -T • N • δu dA + Γ vess [k v u + c v u] • δu dA (6)
with density ρ 0 , accelerations ü, the second Piola-Kirchhoff stress tensor S, and spring stiffnesses k v , k a and viscosities c v , c a for vessel and apical surfaces Γ vess and Γ apex , respectively. Furthermore, we have virtual displacements δu and virtual strains δE. Omni-directional springs and dashpots are placed on the outsides of the great vessels Γ vess and apical surface Γ apex . 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 [START_REF] Standring | s Anatomy: The Anatomical Basis of Clinical Practice[END_REF], ventricular and atrial myocardium [START_REF] Boughner | Mechanical properties of human pericardium. Differences in viscoelastic response when compared with canine pericardium[END_REF], and aorta and pulmonary artery ( 9):

S = ∂ ∂E (ψ NH + ψ vol ) + ∂ ∂ Ė ψ visco , (7)
S = ∂ ∂E (ψ exp + ψ vol )+ ∂ ∂ Ė ψ visco + S act , (8)
S = ∂ ∂E (ψ MR + ψ vol ) + ∂ ∂ Ė ψ visco . ( 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 S act . The strain energy density functions for exponential orthotropic solid ψ exp [START_REF] Holzapfel | Constitutive modelling of passive myocardium: a structurally based framework for material characterization[END_REF], Mooney-Rivlin solid [START_REF] Sainte-Marie | Modeling and estimation of the cardiac electromechanical activity[END_REF][START_REF] Chabiniok | Estimation of tissue contractility from cardiac cine-MRI using a biomechanical heart model[END_REF] ψ MR , Neo-Hookean solid ψ NH , penalty function ψ vol , and viscous pseudo-potential [START_REF] Chapelle | An energypreserving muscle tissue model: formulation and compatible discretizations[END_REF] ψ visco are given as

ψ exp = a 2b e b( Ī1-3) -1 + a f s 2b f s e b f s I 2 8,f s -1 + i∈{f,s} a i 2b i e bi(I4,i-3) -1 , ψ MR = C 1 ( Ī1 -3) + C 2 ( Ī2 -3), ψ NH = µ 2 ( Ī1 -3), ψ vol = κ 2 (1 -J) 2 , ψ visco = η 2 tr Ė2 , (10)
with the Jacobian J = det F of the deformation gradient and material parameters a i , b i , µ, C 1 , and C 2 ,

Ī1 = J -2/3 I 1 , Ī2 = J -4/3 I 2 , I 1 = tr C, I 2 = 1 2 tr 2 (C) -tr C 2 , I 4,f = f 0 • Cf 0 , I 4,s = s 0 • Cs 0 , I 8,f s = f 0 • Cs 0 . (11)
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

S act = τ (t) • f 0 ⊗ f 0 ( 12
)
with fiber stress τ in local reference fiber direction f 0 .

Based on [START_REF] Bestel | A Biomechanical Model of Muscle Contraction[END_REF], we model fiber stress by the evolution equation

τ (t) = -|a(t)|τ (t) + σ 0 |a(t)| + (13)
with activation function a, contractility σ 0 , and the function |a(t)| + = max(a(t); 0). The activation function a(t) is modeled by

a(t) = α max • f (t) + α min • [1 -f (t)] (14)
with maximum and minimum activation rates α max and α min , respectively, and functions

f (t) = S + (t -t sys ) • S -(t -t dias ), (15)
S ± (∆t) = 1 2 1 ± tanh ∆t γ (16)
with steepness γ = 0.005 s and descending and ascending sigmoid functions S + and S -, respectively. The indicator function f ∈ ]0, 1[ indicates systole. The times t sys and t dias model the onset of systole and diastole, respectively. Note that the active stress tensor S act 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 [START_REF] Hörmann | An adaptive Hybridizable Discontinuous Galerkin approach for cardiac electrophysiology[END_REF]. 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 1a, we obtain the active stress curve depicted in figure 6b. The values σ a and σ v denote the maximum of the active stress τ for atria and ventricles, respectively. Using the finite element method, we discretize displacements u and virtual displacements δu arising in the weak form [START_REF] Martini | Human Anatomy[END_REF] by

u (e) = ϕ (e) (X) d (e) δu (e) = ϕ (e) (X) δd (e) (17)
with quadratic basis functions ϕ and nodal displacements d on each tetrahedral element (e). Assembly of the discretized problem leads to the matrix notation

M d + F(d, ḋ, p v ) = 0, (18)
with mass matrix M, force vector F, and discrete displacements, velocities, and accelerations d, ḋ, and d, respectively. We discretize the boundary value problem in time with Newmark's method [START_REF] Newmark | A method of computation for structural dynamics[END_REF]

ḋn+1 = γ β∆t (d n+1 -d n ) - γ -β β ḋn - γ -2β 2β ∆t dn , dn+1 = 1 β∆t 2 (d n+1 -d n ) - 1 β∆t ḋn - 1 -2β 2β dn , (19)
with parameters γ ∈ [0, 1] and β ∈ [0, 0.5], and time step size ∆t = t n+1 -t n . Additionally, we apply the generalized-α method [START_REF] Chung | A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-α method[END_REF], yielding quantities at a gen-

eralized time step n + 1 -α i (•) n+1-αi = (1 -α i )(•) n+1 + α i (•) n , α i ∈ [0, 1], i ∈ {f, m} (20)
depending on the weights α f and α m 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

R S = M dn+1-αm + F n+1-α f . (21)
All parameters used for the elastodynamical model are summarized in table 1a.

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 [START_REF] Westerhof | The arterial windkessel[END_REF] and [START_REF] Kerckhoffs | Coupling of a 3D finite element model of cardiac ventricular mechanics to lumped sys-tems models of the systemic and pulmonic circulation[END_REF]. A comprehensive review of different windkessel models is given in [START_REF] Shi | Review of zero-D and 1-D models of blood flow in the cardiovascular system[END_REF]. The schematic of our windkessel model is given in figure 7a using resistances R, compliances C, and an inertance L p . Pressures at different parts of the model are denoted by p. We distinguish between a proximal (index p) and a distal part (index d) of the outlets, i.e. lung and aorta for the right and left ventricle, respectively. The atrial pressure p at is prescribed to simulate atrial systole, see figure 7b. The reference pressure p ref is kept constant.

We model the atrioventricular and semilunar valves with a smooth diode-like behavior by non-linear resistances R av := R(p v -p at ) and R sl := R(p p -p v ), respectively, depending on the sigmoid function

R(∆p) = R min + (R max -R min ) • S + (∆p) (22)
with the steepness k p of the sigmoid function S + and the minimal and maximal valve resistance R min → 0 and R max → ∞, respectively. This yields the set of differential equations

p v -p at R av + p v -p p R sl + V (u) = 0, q p - p v -p p R sl + C p ṗp = 0, q p + p d -p p R p + L p R p qp = 0, p d -p ref R d -q p + C d ṗd = 0. ( 23
)
The 0D windkessel model is strongly coupled to the 3D structural model. The 0D model depends on the structural displacements u of the 3D model via the change in ventricular volume V . 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 p = [p v , p p , p d , q p ] T , including the flux q p through the inertance L p . We discretize the set of windkessel equations [START_REF] Sutton | Measurement of postoperative pericardial pressure in man[END_REF] in time with the one-step-θ scheme

(•) n+1 = (•) n+1 -(•) n ∆t , θ ∈ [0, 1], (•) n+θ = θ(•) n+1 + (1 -θ)(•) n . (24)
This yields the discrete windkessel residual R 0D evaluated at time step n + θ. Proximal flow q p (0) 38.3 14.9 cm 3 s (c) Initial conditions of the reduced order cardiovascular model.

Par. Value

Generalized-α γ, α f , α m 0.5

β 0.25 One-step-θ θ 1.0 (d) Numerical time integration parame- ters.
Table 1: Overview of parameters in four-chamber cardiac model.

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 1c.

In this work, we use the Modified Updated Lagrangian Formulation as proposed in [START_REF] Gee | Prestressing in finite deformation abdominal aortic aneurysm simulation[END_REF][START_REF] Gee | A computational strategy for prestressing patient-specific biomechanical problems under finite deformation[END_REF]. 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 [START_REF] Jöbsis | The visceral pericardium: macromolecular structure and contribution to passive mechanical properties of the left ventricle[END_REF].

Solving the 0D-3D coupled problem

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

   ∂R S ∂d ∂R S ∂p ∂R 0D ∂d ∂R 0D ∂p    i n+1 •    ∆d ∆p    i+1 n+1 = -    R S R 0D    i n+1 , (25)
in a monolithic fashion for increments in displacements and windkessel variables ∆d n+1 and ∆p n+1 , respectively, at iteration i + 1 until convergence. We build and solve this coupled system using our in-house code BACI [START_REF] Wall | Baci: a parallel multiphysics simulation environment[END_REF]. The numerical parameters for the time integration of our model are listed in table 1d.

Results

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

Case apex depicted in figure 8a yields the boundaryvale problem adding the energy for omni-directional spring dashpots to the energy [START_REF] Martini | Human Anatomy[END_REF], where Γ apex 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 Γ epi 0 \ Γ apex 0 , i.e. the absence of any pericardial boundary conditions as frequently found in literature [START_REF] Sermesant | Patient-specific electromechanical models of the heart for the prediction of pacing acute effects in CRT: A preliminary clinical validation[END_REF][START_REF] Chabiniok | Estimation of tissue contractility from cardiac cine-MRI using a biomechanical heart model[END_REF][START_REF]Anatomically accurate high resolution modeling of human whole heart electromechanics: A strongly scalable algebraic multigrid solver method for nonlinear deformation[END_REF].

0 = δΠ + Γ apex [k a u + c a u] • δu dA ( 26
)
Case pericardium depicted in figure 8b yields the boundary-vale problem

0 = δΠ + Γ epi N [k p u • N + c p u • N ] • δu dA ( 27
)
adding the energy for the pericardial boundary condition to the energy [START_REF] Martini | Human Anatomy[END_REF], where Γ epi is the pericardial surface. The choice of parameters k v and c v 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.

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 9a), three-, and twochamber views and short axis planes with a slice distance of 8 mm (figures 9c, 9d, 9b). 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 freebreathing, 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.

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 σ 0 . It is a key parameter describing cardiac elastodynamics, resembling the asymptotical active fiber stress in [START_REF] Sacks | Incorporation of experimentally-derived fiber orientation into a structural constitutive model for planar collagenous tissues[END_REF]. It controls maximum deformation during systole. In order to make simulations comparable, we adapt σ 0 for each combination of boundary condition and fiber orientation to match the left ventricular volume at end-systole as segmented from cine MRI of V min = 57 ml. 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 [START_REF] Carlsson | Atrioventricular plane displacement is the major contributor to left ventricular pumping in healthy adults, athletes, and patients with dilated cardiomyopathy[END_REF][START_REF] Willenheimer | Left ventricular atrioventricular plane displacement: an echocardiographic technique for rapid assessment of prognosis in heart failure[END_REF]. 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 [START_REF] Heiberg | Design and validation of Segmentfreely available software for cardiovascular image analysis[END_REF].

In our simulations, we average the displacements on all nodes on the valve plane (see the red planes in figure 5c) 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 [START_REF] Chabiniok | Estimation of tissue contractility from cardiac cine-MRI using a biomechanical heart model[END_REF] 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 9b. These slices are selected because the myocardium is recognizable for all MRI time steps and not disturbed by either apex or AVP. The function A converts both MRI and simulated endocardial contours d s MRI and d s , respectively, to binary images with a resolution of 1 × 1 mm 2 for every slice s. We use the Dice metric to compare both binary images = 1 -1 5

9 s=5 2 |A (d s MRI ) ∩ A (d s )| |A (d s MRI )| + |A (d s )| ∈ [0, 1] (28)
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 tapex (t) = Γ apex t apex da

Γ apex 1 da , t apex = k a u + c a u. ( 29
)
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 t epi and integrate it over the epicardial surface to obtain the mean pericardial stress tepi (t) = Γ epi t epi da

Γ epi 1 da , t epi = k p u • N + c p u • N , (30)
normalized by the epicardial area.

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 k p and c v 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 k p on the contraction of the heart. For this study, we limit ourselves to the ±60 • fiber distribution, as it is commonly used in cardiac simulations, see e.g. [START_REF] Chabiniok | Estimation of tissue contractility from cardiac cine-MRI using a biomechanical heart model[END_REF][START_REF] Sundnes | An integrated electromechanical-growth heart model for simulating cardiac therapies[END_REF][START_REF] Hirschvogel | A monolithic 3D-0D coupled closed-loop model of the heart and the vascular system: Experiment-based parameter estimation for patient-specific cardiac mechanics[END_REF]. We tested the following parameter values: k p ∈ {0.1, 0.2, . . . , 1.0, 1.5, . . . , 5.0} kPa mm [START_REF]Anatomically accurate high resolution modeling of human whole heart electromechanics: A strongly scalable algebraic multigrid solver method for nonlinear deformation[END_REF] For each choice of k p , we calibrated active stress to yield the same end-systolic volume as measured from MRI.

All parameters except k p 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 endsystolic 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 ±60 • 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 k p , as shown in figure 11. For high k p , contact stress has the same order of magnitude as active stress and exceeds maximum left ventricular pressure. For small k p , it has the same order of magnitude as atrial pressure, which experimentally shown to be a good predictor for pericardial contact stress [START_REF] Tyberg | The relationship between pericardial pressure and right atrial pressure: an intraoperative study[END_REF].

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 k p .

The end-systolic state of the simulations is shown in figure 13 compared to MRI. The images contain all simulated variants for k p , where the color changes continuously from k = 0.1 kPa/mm (blue) to k = 5.0 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 13a for the atria in four-chamber view. The short axis views in figures 13b and 13c additionally show that the interventricular septum is stretched and rotated as compared to MRI for high k p .

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

Left and right AVPD is displayed in figures 14c and 14d, respectively. In the left ventricle, i.e. at the mitral valve, AVPD is not very sensitive to the choice of k p . 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 k p towards the value measured in MRI. An identical trend is observable for left and right atrial volume in figures 14e and 14f, respectively.

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

Model personalization

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

The onset of systole and diastole, t sys and t dias , as well as the myofiber activation and deactivation rates, α max and α max , are adapted to the left ventricular volume curve for ventricles and atria. Here, parameters for the atria are fitted from the interval t ∈ [0, 0.2 s] and t ∈ [0.2 s, 0.9 s] for the ventricles. The material parameter η controlling the viscosity of the tissue is fitted during ventricular diastole, i.e. t ∈ [0.5 s, 0.9 s]. 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

±50 • ±60 • ±70 • ±50 • ±60 • ±70 • apex 63.

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 σ 0 was calibrated in all simulations to match end-systolic volume as segmented from cine MRI. Therefore, in figures 15a and 15c the volumes of MRI and all simulations match at t = 0.51 s. Furthermore, although axis views, see figures 16e, 16h, 16f, 16i. The location of the view planes is visualized in figures 9a, 9c, 9d. For case apex , there is a radial inward movement of the myocardial wall. In figure 16d, 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 16e and 16f. However, the mismatch of epicardial contours is clearly visible and sensitive to fiber orientation.

Comparing figures 16d and 16g, the influence of the pericardial boundary condition becomes clearly visible. It can be observed for case pericardium in figure 16g that the epicardial contour matches the MRI much closer than case apex in figure 16d 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 16d and 16g, 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 ±50 • (red) to ±70 • (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 16h. For short axis slice 6 in figure 16i there is a good agreement with simulation and MRI at all regions of the left and right myocardium. 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 t = 0.17s. Followed by ventricular systole, the AVPs move towards the apex and both reach their maximal value at ventricular end-systole at t = 0.51s. During atrial systole for t ∈ [0, 0.25s], 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 17a, 17c, 17b, 17d, one can observe that in both cases maximum AVPD depends on fiber orientation: Maximum AVPD increases from horizontal ±50 • fibers (red) to vertical ±70 • 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.

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 18a and 18b stay constant during ventricular systole. This is also visible at the end-systolic four-chamber view in figure 16d. For ±70 • fibers, the right atrium is even slightly emptied during ventricular systole, as observed in figure 18b. Atrial filling can be observed for case pericardium in figures 18c and 18d. Both atria are visibly filled during ventricular systole, although maximum atrial volume remains smaller than in MRI.

Spatial error

For case apex in figures 19a and 19b the error is lowest in both ventricular endocardia during contraction at end-systole at t = 0.51. The error rises during ventricular contraction and relaxation. Errors at the end of the simulation higher than the ones at t = 0 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 .

Boundary stresses

Both scalar boundary stresses tapex = tapex 2 and tepi 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 tepi indicate predominant tensile stresses between epicardium and pericardium. It can be seen that mean pericardial stress in figure 20b is a compressive stress for most of the cardiac cycle, except at the end of systole and onset of diastole.

Boundary stresses are visualized in figure 21. For case apex , the mean stress vectors tapex for all three fiber distributions are shown in figure 21a at t = 0.45 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 ±60 • fibers at end-systole is shown in figure 21b 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.

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 ±50 • , ±60 • , and ±70 • . We then compared scalar left ventricular pressure and volume. The displacements at endsystole 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.

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. k p = 0.1 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. k p = 3.0 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 k p = 0.1 kPa/mm for pericardial stiffness yields good results for a biventricular geometry with ±60 • fibers. This value was also used in [START_REF] Hirschvogel | A monolithic 3D-0D coupled closed-loop model of the heart and the vascular system: Experiment-based parameter estimation for patient-specific cardiac mechanics[END_REF], although it was not really analyzed there, e.g. with respect to MRI.

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 [START_REF] Hammond | Heart size and maximal cardiac output are limited by the pericardium[END_REF]. The result further agrees with the numerical experiments performed in [START_REF] Fritz | Simulation of the contraction of the ventricles in a human heart model including atria and pericardium. Biomechanics and Modeling in Mechanobiology[END_REF]. 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 [ [START_REF] Arutunyan | Atrioventricular plane displacement is the sole mechanism of atrial and ventricular refill[END_REF][START_REF] Arvidsson | Letter to the Editor: Atrioventricular plane displacement is not the sole mechanism of atrial and ventricular refill[END_REF]. In [START_REF] Maksuti | Modelling the heart with the atrioventricular plane as a piston unit[END_REF], 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 [START_REF] Fujii | Effect of left ventricular contractile performance on passive left atrial filling -Clinical study using radionuclide angiography[END_REF]. 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 [START_REF] Emilsson | The mode of left ventricular pumping: is there an outer contour change in addition to the atrioventricular plane displacement?[END_REF] 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 [START_REF] Bowman | Assessment and consequences of the constant-volume attribute of the fourchambered heart[END_REF][START_REF] Arvidsson | Letter to the Editor: Atrioventricular plane displacement is not the sole mechanism of atrial and ventricular refill[END_REF]. We found that for the ±60 • 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.

Fiber orientation

We compared in this work three fiber orientations within the myocardium, namely ±50 • , ±60 • , and ±70 • . 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 ±70 • 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 ±60 • and ±50 • 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 ±70 • fiber orienta- tion also yielded the highest atrial filling during ventricular systole. Comparing results to short axis cine MRI slices, it was shown that a more horizontal ±50 • fiber orientation leads to a more radial contraction of the heart. The maximum pericardial stress at end-systole was highest for ±50 • fibers. This can be explained by the observation in figure 16h where ±50 • 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.

Pericardial contact stress

In [START_REF] Smiseth | Assessment of pericardial constraint in dogs[END_REF] 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 [START_REF] Tyberg | The relationship between pericardial pressure and right atrial pressure: an intraoperative study[END_REF] 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 enddiastolic 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.

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 ±60 • 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.

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 [START_REF] Hörmann | An adaptive Hybridizable Discontinuous Galerkin approach for cardiac electrophysiology[END_REF], 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 [START_REF] Yin | Quantification of the mechanical properties of noncontracting canine myocardium under simultaneous biaxial loading[END_REF][START_REF] Dokos | Shear properties of passive ventricular myocardium[END_REF][START_REF] Sommer | Biomechanical properties and microstructure of human ventricular myocardium[END_REF] showed anisotropic tissue characteristics, depending on myocardial fiber and sheet orientation. In our model, we used the anisotropic material model proposed in [START_REF] Holzapfel | Constitutive modelling of passive myocardium: a structurally based framework for material characterization[END_REF] 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 [START_REF] Holzapfel | Constitutive modelling of passive myocardium: a structurally based framework for material characterization[END_REF] and [START_REF] Gültekin | An orthotropic viscoelastic model for the passive myocardium: continuum basis and numerical treatment[END_REF] 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 patientspecific fiber orientation and assumed equal distributions in left and right ventricle. Patient-specific cardiac fiber orientations can be estimated from diffusion tensor MRI [START_REF] Nagler | Maximum likelihood estimation of cardiac fiber bundle orientation from arbitrarily spaced diffusion weighted images[END_REF] (DTMRI). However, while applicable to in vivo DTMRI (as shown in [START_REF] Nagler | Maximum likelihood estimation of cardiac fiber bundle orientation from arbitrarily spaced diffusion weighted images[END_REF]), 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 16g 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 Eule-rian 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 endsystole. 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.

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 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 ±60 • fibers, which roughly represents the shape of the left ventricle, see figure 22a. 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 22b and 22c at end-systole. Displayed is the reference configuration and all three boundary condition cases for a cross-section of the ellipsoid. Figure 22b 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 22b 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.

  (a) Dissected mediastinum with cut pericardium and heart surface. Image by G. M. Gruber, Medical University of Vienna, Austria. Location of the heart with respect to serous and fibrous pericardium. Inspired by[START_REF] Iaizzo | Handbook of cardiac anatomy, physiology, and devices[END_REF]. Cross-sectional view of transmural layers of heart and pericardium. Inspired by[START_REF] Iaizzo | Handbook of cardiac anatomy, physiology, and devices[END_REF].

Fig. 1 :

 1 Fig. 1: Heart and pericardium.

Fig. 2 :

 2 Fig. 2: Position of the pericardium indicated in 3D MRI taken during diastasis. The neighboring tissue is colorcoded: 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.

Fig. 3 :

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

  Reference normal spring from (4).

  (a) Computational mesh. (b) Atrial fibers and ventricular ±60 • fibers at endocardium. Color shows helix angle with respect to long axis.

  Prescribed indicator functions f (t) for atria (blue) and ventricles (Time [s] Active stress τ [Pa] (b) Active stress τ (t) for atria (blue) and ventricles (red) with maximum values σ a and σ v , respectively.

Fig. 6 :

 6 Fig. 6: Active stress.

  Atrial pressure [mmHg] (b) Prescribed left (blue) and right (red) atrial pressure p at (t).

Fig. 7 :

 7 Fig. 7: Windkessel model.

( a )

 a Case apex with omni-directional springdashpots on Γ apex (green). (b) Case pericardium with normal spring-dashpots on Γ epi (red).

Fig. 8 :

 8 Fig. 8: Surface definitions for boundary conditions with omni-directional spring-dashpots on Γ vess (blue) and homogeneous Neumann boundary conditions (white).

  (a) Four chamber view. (b) Short axis endocardial contours used for error calculation. (c) Short axis slice 9. (d) Short axis slice 6.

Fig. 9 :

 9 Fig. 9: Post processing planes for simulations and cine MRI.

Fig. 10 :Fig. 11 :

 1011 Fig.10: Maximum ventricular active stress σ v calibrated to yield identical end-systolic volume. Shown for case pericardium with varying pericardial stiffness compared to case apex , both with ±60 • fiber distributions.

Fig. 12 :

 12 Fig.12: Volume change for case pericardium with varying pericardial stiffness compared to case apex , both with ±60 • fiber distributions.

Fig. 13 :

 13 Fig.13: Cine MRI at end-systole for case pericardium with ±60 • fiber distribution from k p = 0.1 (blue) to k p = 5.0 (red). Views as defined in figure9. MRI courtesy of R. Chabiniok, J. Harmer, E. Sammut, King's College London, UK.

  Spring stiffness and dashpot viscosity on apical and epicardial surface. σv [kPa] tsys [ms]

Fig. 14 : 3 . 5 Fig. 15 :

 143515 Fig.14: Kinematic scalar cardiac quantities at end-systole for case pericardium with varying pericardial stiffness k p ∈ [0.1, 5.0] compared to case apex both with ±60 • fiber distributions and MRI.

Fig. 16 :

 16 Fig. 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

Fig. 17 :

 17 Fig. 17: Simulated atrioventricular plane displacement for left and right ventricle compared to cine MRI.

Fig. 18 :

 18 Fig. 18: Simulated volume curves for left and right atrium compared to 3D MRI at ventricular diastasis and end-diastole.

Fig. 19 :Fig. 20 :

 1920 Fig. 19: Relative spatial error of simulation results and cine MRI at left and right endocardium.

( a )

 a Case apex mean apical stress vectors tapex (scaled by magnitude) for ±50 • (red), ±60 • (green), and ±70 • (blue) fiber orientations at t = 0.45. (b) Case pericardium pericardial contact stress t epi on epicardial surface with ±60 • fibers at end-systole t = 0.51.

Fig. 21 :

 21 Fig. 21: Visualization of boundary stresses.

  (a) Ellipsoid model in reference configuration with cross-section (yellow). (b) Frontal view of epi-and endocardial cross-section at endsystole. (c) Top-down view of epicardial cross-section at end-systole.

Fig. 22 :

 22 Fig. 22: Ellipsoid model in reference configuration (black) with cases free (green), pericardium (blue), and pseudocontact (orange).

  Parameters of the elastodynamical model.

Table 2 :

 2 Calibrated parameters for simulation cases apex and pericardium and different fiber orientations.

given in table 2. For parameters identical in all models see table 1.

Table 3 :

 3 Overview of material parameters in ellipsoid model. For numerical parameters see table 1d.