Validation of a 3D-0D closed-loop model of the heart and circulation -- Modeling the experimental assessment of diastolic and systolic ventricular properties
Christoph M. Augustin, Matthias A. F. Gsell, Elias Karabelas, Erik Willemen, Frits Prinzen, Joost Lumens, Edward J. Vigmond, Gernot Plank
VValidation of a 3D–0D closed-loop model of the heart and circulation –Modeling the experimental assessment of diastolic and systolicventricular properties
Christoph M. Augustin a , Matthias A. F. Gsell a , Elias Karabelas a , Erik Willemen c , Frits Prinzen c ,Joost Lumens c , Edward J. Vigmond d , Gernot Plank a,b, ∗ a Gottfried Schatz Research Center: Division of Biophysics, Medical University of Graz, Graz, Austria b BioTechMed-Graz, Graz, Austria c Department of Biomedical Engineering, CARIM School for Cardiovascular Diseases, Maastricht University,Maastricht, Netherlands d IHU Liryc, Electrophysiology and Heart Modeling Institute, fondation Bordeaux Université, Pessac-Bordeaux,France
Abstract
Computer models of cardiac electro-mechanics (EM) and cardiovascular hemodynamics show promiseas an effective means for the quantitative analysis of clinical data and, potentially, provide capa-bilities to predict the response to therapies. To realize such modeling applications methodologicalkey challenges must be addressed. To endow models with predictive capabilities they must com-prehensively encapsulate relevant cardiovascular mechanism, achieve computational efficiency androbustness to enable simulations of prolonged observation periods under a broad range of conditions,and, thus, facilitate model personalization within tractable time frames.Here, these challenges are addressed by introducing a flexible method for coupling a 3D modelof bi-ventricular EM to the 0D
CircAdapt model representing atrial EM and closed-loop circulation.A detailed mathematical description is given and efficiency, robustness and accuracy of numericalscheme and solver implementation are evaluated. After parameterization and stabilization to alimit cycle under baseline conditions, the model is shown to replicate physiological behaviors bysimulating the transient response to alterations in loading conditions and contractility, induced byprotocols for assessing of systolic and diastolic ventricular properties.
Keywords: ventricular pressure-volume relation, Frank-Starling mechanism, ventricular load
1. Introduction
Cardiovascular diseases (CVDs) are the primary cause of mortality and morbidity in indus-trialized nations, posing a significant burden on health care systems worldwide [1, 2, 3]. Despitecontinuous diagnostic and therapeutic advances, their optimal treatment remains a challenge [4].In no small part, this is due to the complex multiphysics nature of cardiovascular function – theheart is an electrically controlled mechanical pump driving blood through the circulatory system. ∗ Address correspondence to Gernot Plank, Gottfried Schatz Research Center: Division of Biophysics, MedicalUniversity of Graz, Neue Stiftingtalstrasse 6/IV, Graz 8010, Austria. Email: [email protected]
Preprint submitted to Elsevier December 2, 2020 a r X i v : . [ q - b i o . T O ] D ec dvanced clinical modalities provide a wealth of disparate data, but effective tools allowing theircomprehensive quantitative analysis are lacking. Computer models able to capture mechanisticrelations between clinical observations quantitatively show promise to fill this void. In recent singlephysics cardiac electrophysiology (EP) studies, the added value of models in improving therapystratification [5] and planning [6, 7] has been demonstrated already.Models of cardiovascular EM are even more demanding and challenging to apply in a clinicalcontext. Their utility depends on the ability to comprehensively represent mechanisms underlyinga broader range of physiological function, and to tailor these to approximate – with acceptablefidelity – anatomy and cardiovascular function of a given patient. Such models are complex asthe interplay between all major components governing a heart beat must be taken into account.These comprise models of cardiac EP producing electrical activation and repolarization patternsthat drive EM coupling to models of contractile function, cardiac mechanics describing deformationand stresses under given mechanical boundary and hemodynamic loading conditions imposed bythe intra-thoracic embedding of the heart and the circulatory system.Pumping function is regulated through a bidirectional interaction between circulatory systemand the heart. The circulatory system as an extracardiac factor imposes a pressure and volume loadupon the heart and, vice versa , pressure and flow in the circulatory system are determined by themechanical state of cardiac cavities. Optimal function depends on matching the coupling betweenthese two systems [8]. From a physics point of view, coupling poses a fluid-structure interaction(FSI) problem, with pressure and blood flow velocity fields as coupling variables [9, 10]. These arerelevant for investigating flow patterns or wall shear stresses, but are less suitable for systems levelinvestigations. Simpler, computationally less costly 0D and 1D lumped models have been preferredto provide appropriate hemodynamic loading conditions to the heart.Most EM modeling studies consider ventricular afterload only represented by lumped 0D Wind-kessel type models comprising 2-, 3- or 4-elements [11, 12, 13, 14, 15, 16, 17], or, less common,by 1D models derived from Navier–Stokes equations [18, 19, 20, 21, 22]. The latter also accountfor pulse wave transmission and reflection, but identifying parameters is more challenging than for0D models [23]. A fundamental limitation of isolated models of pre- and afterload is the lack ofregulatory loops which respond to altered loading or contractility in one chamber by balancingpreload conditions in all chambers until a new stable limit cycle with common compatible strokevolumes is reached. Isolated afterload models are thus suitable only for approximating the immedi-ate responses in a single beat [24], but less so for predicting transient behaviors over multiple beats.Closed-loop circulatory systems [25, 26, 27, 28, 29] take into account these feedback mechanismsand ascertain the conservation of blood volume throughout the cardiovascular system.Achieving a flexible, robust and efficient coupling of 3D EM models of cardiac chambers to a 0Dclosed-loop model of the circulatory system remains challenging. Hydrostatic pressure, p , in cavitiesand blood flow, q , between cavities and circulatory system serve as coupling variables that act aspressure boundary condition and impose volume constraints on the 3D cavity models. Previousstudies addressed 3D solid-0D fluid coupling problems using simpler partitioned [30, 31, 32] ormore advanced strongly coupled monolithic approaches [33, 34, 35, 36, 23, 37]. Yet, reports oncoupling of a closed-loop 0D to 3D solid models are sparse. These mostly used simplified circuitmodels [35, 31, 36, 38], that were used for simulations of a single heart beat. Model validation todemonstrate agreement with known physiological behaviors has been limited.Based on previous work on cardiac EM models [39, 40, 41] we report on the development ofa monolithic 3D solid - 0D fluid coupling approach. Feasibility is demonstrated by building a3D canine bi-ventricular EM PDE model coupled to the state-of-the art CircAdapt model [25,22] – a non-linear 0D closed-loop model of the cardiovascular system that implements dynamicadaptation processes based on known physiological principles – to represent atrial chambers aswell as systemic and pulmonary circulation. A detailed description of numerical underpinnings isgiven, including a complete mathematical description of the
CircAdapt model that has been lackingso far. Efficiency, robustness and accuracy of numerical scheme and solver implementation areevaluated. The coupled model is first parameterized and stabilized to a limit cycle representingbaseline conditions, and, then, validated by demonstrating its ability to replicate physiologicalbehaviors under various protocols altering loading conditions and contractility that are used for theexperimental assessment of systolic and diastolic ventricular properties. Transient responses underthese protocols are simulated over prolonged observation periods, covering up to 25 beats.
2. Methodology
In a previous study mongrel dog data were acquired to show the maintenance cardiac contractilecoordination, pump function and efficiency of left ventricular (LV) septal and LV apical pacing, see[43]. The animals were handled according to the Dutch Law on Animal Experimentation (WOD)and the European Directive for the Protection of Vertebrate Animals Used for Experimental andOther Scientific Purposes (86/609/EU). The protocol was approved by the Maastricht UniversityExperimental Animal Committee. Anatomical Magnetic Resonance Images (MRI) were acquiredon a Philips Gyroscan . (NT, Philips Medical Systems, Best, the Netherlands) using a stan-dard synergy receiver coil for thorax examinations. Images of seven short-axis cross-sections ofslice thickness with inter-slice distance were obtained to capture the whole heart. LVpressure and volume were determined using the conductance catheter technique (CD-Leycom, TheNetherlands), see [44], and the signals were digitized at . Multilabel segmentations of right ventricular (RV) (tag ) and LV blood pool (tag ) andof the LV myocardium (tag ) were generated from seven MRI short axis slices using Seg3D [45].Each slice was first segmented semi-automatically using thresholding techniques with manual cor-rection. Segmentations were upsampled to isotropic resolution, followed by an automated iterativeerosion and dilation smoothing scheme implemented in Meshtool [46]. The RV wall (tag ) andlids representing the atrio-ventricular valves were automatically generated by dilation of the adja-cent blood pool (tags and − ◦ at theepicardium to +60 ◦ at the endocardium [48]. Universal ventricular coordinates were computed [49]to support the flexible definition of stimulation sites and mechancal boundary conditions. Twomeshes of different resolution were generated, a coarse mesh to reduce computational expenses andto facilitate the fast exploration of experimental protocols over prolonged observation periods, anda higher resolution mesh for investigating potential inaccuracies introduced by the coarser spatialresolution. For the coarse mesh, average edge lengths of ∼ . and ∼ . , were chosen inLV and RV, respectively, to ascertain that at least two elements were generated transmurally acrossthe myocardial walls, as illustrated in Fig. 1. For the finer mesh, average edge lengths of ∼ . and ∼ . , were chosen in LV and RV, respectively.3 igure 1: coarse (A) and fine (B) resolution meshes with domain labels and corresponding fiber fields. Note thedifference in fiber angles due to spatial resolution. Cardiac tissue is mechanically characterized as a hyperelastic, nearly incom-pressible, orthotropic material with a nonlinear stress-strain relationship. The deformation gradi-ent F describes the deformation u of a body from the reference configuration Ω ( X ) to the currentconfiguration Ω t ( x ) , F ij = ∂x i ∂X j , i, j = 1 , , . (1)By convention, we denote J = det F > and introduce the right Cauchy–Green tensor C = F (cid:62) F .The nearly incompressible behavior is modeled by a multiplicative decomposition of the deformationgradient, see, e.g., [50], of the form F = J / F , C = J / C , with det F = det C = 1 . (2)Mechanical deformation is described by the equilibrium equations given as ρ ¨ u ( t, X ) − Div [ FS ( u , X )] = for X ∈ Ω × (0 , T ) , (3)with initial conditions u ( X ,
0) = , ˙ u ( X ,
0) = . Here, ρ is the density in reference configuration; ¨ u are nodal accelerations; ˙ u are nodal velocities; S ( u , X ) is the second Piola–Kirchhoff stress tensor; and Div denotes the divergence operator in thereference configuration.The boundary of the bi-ventricular models was decomposed in three parts, ∂Ω = Γ endo , ∪ Γ epi , ∪ Γ base , , with Γ endo , the endocardium, Γ epi , the epicardium, and Γ base , the base of the4entricles.Normal stress boundary conditions were imposed on the endocardium FS ( u , X ) n out0 ( X ) = − p ( t ) J F −(cid:62) n out0 ( X ) on Γ endo , × × (0 , T ) (4)with n out0 the outer normal vector; omni-directional spring type boundary conditions constrained theventricles at the basal cut plane Γ base , [51]; and to simulate the mechanical constrains imposed bythe pericardium spatially varying normal Robin boundary conditions were applied at the epicardium Γ epi , [7].Apart from external loads the deformation of cardiac tissue is in particular governed by ac-tive stresses intrinsically generated during contraction. To simulate both the active and passiveproperties of the tissue, the total stress S is additively decomposed according to S = S p + S a , (5)where S p and S a refer to the passive and active stresses, respectively. Passive Stress.
Passive stresses are modeled based on the constitutive equation S p = 2 ∂Ψ ( C ) ∂ C , (6)where Ψ is an invariant-based strain-energy function to model the orthotropic behavior of cardiactissue. The prevailing orientation of myocytes, referred to as fiber orientation, is denoted as f .Individual myocytes are surrounded and interconnected by collagen, forming sheets, which is de-scribed by the sheet orientation s , perpendicular to f . Together with the sheet-normal axis n ,orthogonal to the sheet and the fiber orientations, this forms a right-handed orthonormal set ofbasis vectors.Following Usyk et al. [52] the orthotropic constitutive relation is defined as Ψ ( C ) = κ J ) + a Q ) − , (7)where the first term is the volumetric energy with the bulk modulus κ (cid:29) which penalizeslocal volume changes to enforce near incompressible behavior of the tissue and the term in theexponent is Q = b ff E + b ss E + b nn E + b fs (cid:16) E + E (cid:17) + b fn (cid:16) E + E (cid:17) + b ns (cid:16) E + E (cid:17) . (8)Here, the directional strains read E ff = f · Ef , E ss = s · Es , E nn = n · En , E fs = f · Es , E fn = f · En ,E ns = n · Es , E sf = s · Ef , E nf = n · Ef , E sn = s · En , with E = ( C − I ) the modified isochoric Green–Lagrange strain tensor and parameters givenin Tab. 1. Active Stress.
Stresses due to active contraction are assumed to be orthotropic with full contractileforce along the myocyte fiber orientation f and
40 % contractile force along the sheet orientation5 [53, 54]. Thus, the active stress tensor is defined as S a = S a ( f · Cf ) − f ⊗ f + 0 . S a ( s · Cs ) − s ⊗ s , (9)where S a is the scalar active stress describing the contractile force. A simplified phenomenologicalcontractile model was used to represent active stress generation [24]. Owing to its small numberof parameters and its direct relation to clinically measurable quantities such as peak pressure, andthe maximum rate of rise of pressure this model is fairly easy to fit and thus very suitable for beingused in clinical EM modeling studies. Briefly, the active stress transient is given by S a ( t, λ ) = S peak φ ( λ ) tanh (cid:18) t s τ c (cid:19) tanh (cid:18) t dur − t s τ r (cid:19) , for < t s < t dur , (10)with φ = tanh(ld( λ − λ )) , τ c = τ c + ld up (1 − φ ) , t s = t − t a − t emd (11)and t s is the onset of contraction; φ ( λ ) is a non-linear length-dependent function in which λ isthe fiber stretch and λ is the lower limit of fiber stretch below which no further active tensionis generated; t a is the local activation time from Eq. (12), defined when the local transmembranepotential passes the threshold voltage V m , thresh ; t emd is the EM delay between the onsets of electricaldepolarization and active stress generation; S peak is the peak isometric tension; t dur is the durationof active stress transient; τ c is time constant of contraction; τ c is the baseline time constant ofcontraction; ld up is the length-dependence of τ c ; τ r is the time constant of relaxation; and ld is thedegree of length dependence. For the parameter values used in the simulations see Tab. 1. Notethat active stresses in this simplified model are only length-dependent, but dependence on fibervelocity, ˙ λ , is ignored. Electrophysiology.
A recently developed reaction-eikonal (R-E) model [55] was employed to gener-ate electrical activation sequences which serve as a trigger for active stress generation in cardiactissue. The hybrid R-E model combines a standard reaction-diffusion (R-D) model based on themonodomain equation with an eikonal model. Briefly, the eikonal equation is given as (cid:26) (cid:112) ∇ X t (cid:62) a V ∇ X t a = 1 in Ω ,t a = t on Γ ∗ , (12)where ( ∇ X ) is the gradient with respect to the end-diastolic reference configuration Ω ; t a is apositive function describing the wavefront arrival time at location X ∈ Ω ; and t are initialactivations at locations Γ ∗ ⊆ Γ N , . The symmetric positive definite × tensor V ( X ) holds thesquared velocities ( v f ( X ) , v s ( X ) , v n ( X )) associated to the tissue’s eigenaxes f , s , and n . Thearrival time function t a ( X ) was subsequently used in a modified monodomain R-D model given as βC m ∂V m ∂t = ∇ X · σ i ∇ X V m − βI ion + I foot , (13)with β the membrane surface-to-volume ratio; C m the membrane capacitance; V m the unknowntransmembrane voltage; σ i the intracellular conductivity tensor which holds the scalar conductiv-ities ( g f ( X ) , g s ( X ) , g n ( X ) ) and is coupled to V ( X ) proportionally [56]; and I ion the membraneionic current density. Additionally, an arrival time dependent foot current, I foot ( t a ) , was addedwhich is designed to mimic subthreshold electrotonic currents to produce a physiological foot of the6ction potential. The key advantage of the R-E model is its ability to compute activation sequencesat much coarser spatial resolutions that are not afflicted by the spatial undersampling artifactsleading to conduction slowing or even numerical conduction block as it is observed in standard R-Dmodels. Ventricular EP was represented by the tenTusscher–Noble–Noble–Panfilov model of thehuman ventricular myocyte [57]. Computation of Volumes.
To compute the flow into the cardiovascular system, the volume of eachcavity — i.e., left ventricle (lv), left atrium (la), right ventricle (rv), right atrium (ra) — thatis modeled as a 3D elastic body has to be tracked as a function of time: V PDE c ( x , t ) for c ∈{ lv , rv , la , ra } . A reduction in cavitary volume ∂V PDE c ( x , t ) ∂t < , drives a positive flow into the cavitary system. In a purely EM simulation context where thefluid domain is not modeled explicitly, the cavitary volume is not discretized, only the surface Γ c enclosing the cavitary volume is known. Assuming that the entire surface of the cavitary volumeis available, that is, also the faces representing the valves are explicitly discretized, the enclosedvolume V PDE c can be computed from this surface using the divergence theorem V PDE c ( u , t ) = V PDE c ( x , t ) = 13 (cid:90) Γ c,t x · n d Γ c,t . (14)Using this approach, the volume V PDE c ( x , t ) can be computed for each state of deformation at time t and the flow can be derived by differentiating with respect to t . CircAdapt modelCircAdapt [25], as shown schematically in Fig. 2, is a lumped 0D model of the human heart andcirculation. It enables real-time simulation of cardiovascular system dynamics in a wide variety ofphysiological and pathophysiological situations. The entire cardiovascular system is modeled as aconcatenation of modules: a tube module representing the systemic and pulmonary arteries andveins (Appendix A.1); a chamber module modeling actively contracting chambers, i.e., left andright atria and left and right ventricles (Appendix A.3) where myofiber mechanics and contractionis described by a sarcomere module (Appendix A.2); following Lumens et al. [58] this also includesinter-ventricular mechanical interaction through the inter-ventricular septum (Appendix A.4); avalve module representing the aortic, mitral, pulmonary, and tricuspid valves (Appendix A.8); amodule representing systemic and pulmonary peripheral microvasculatures (Appendix A.6); and amodule accounting for effects of the pericardium (Appendix A.5). The modules are connected byflows over valves and venous-atrial inlets (Appendix A.7) and the whole lumped model consists of26 ordinary differential equations (ODEs). In this work we use an adaptive Runge–Kutta–Fehlbergmethod (RKF45) to solve the system of ODEs (Appendix A.9).In Appendix A the mathematical underpinnings of the
CircAdapt model are outlined. Briefly,cavity pressures and cavity volumes are interconnected as follows: volumes regulate cavity wallareas, which in turn determine strain of the myofibers in the wall. Strain is used to calculate my-ofiber stress, (A.9, A.10), which drives wall tension in each cardiac wall (A.22). Using Laplace’s law,transmural pressure is calculated from wall tension and curvature for each wall (A.23). Cavity pres-sures are found by adding the transmural pressures to the intra-pericardial pressure surrounding the7
UBE
Appendix A.1compute A t , p t , Z t (A.1, A.2, A.3) SARCOMERE
Appendix A.2estimate E fib c , σ fib c , κ fib c (A.5, A.11, A.12) CAVITY
Appendix A.3update E fib c (A.17) and C mid c , A mid c , T mid c (A.14, A.16, A.22) TRISEG
Appendix A.4update ˙ y mid , ˙ V midsv (A.24–A.26) CAVITY
Appendix A.3compute A c , Z c , p c (A.18, A.19, A.23) SARCOMERE
Appendix A.2update σ fib c (A.11)compute ˙ L cont c , ˙ C c (A.7,A.8) PERICARD
Appendix A.5compute p peri and update p c (A.30, A.31) VALVE
Appendix A.8compute ˙ q v (A.43) PERIPHERY
Appendix A.6compute q py (A.32) CONNECT
Appendix A.7update p t , p c ,compute ˙ V c , ˙ V t (A.33–A.35) SOLVE
Appendix A.9update V t , V c , C c , L cont c , q v , y mid , V midsv Figure 2:
CircAdapt circulatory system, based on [42], which connects tubes ( t ), cavities ( c ), valves ( v ), and pul-monary and systemic periphery ( py ). In each timestep the ODE system is solved using a Runge–Kutta–Fehlbergmethod, see Appendix A.9, to update the ODE variables (in green), i.e., volumes of tubes ( V t ) and cavities ( V c );sarcomere contractility ( C c ) and sarcomere length ( L cont c ) for each of the cavities and the septum; flow over valves( q v ); and septal midwall volume ( V midsv ) and radius ( y mid ), see Fig. A.9c. In the following steps the updated variablesare used to compute current pressures ( p c , p t ), cross sectional areas ( A c , A t ), and impedances ( Z c , Z t ) for tubes andcavities; fiber strain ( E fib c ), fiber stiffness ( κ fib c ), and fiber stress ( σ fib c ) for the sarcomeres of each cavity and the sep-tum; midwall curvature ( C mid c ), midwall area ( A mid c ), and midwall tension ( T mid c ) for each cavity and the septum;pericardial pressure p peri ; and flow over the systemic and pulmonary periphery q py . myocardial walls (A.31). Consecutively, cavity pressures are used to update flow over valves (A.43)and thus intra-cavitary volumes (A.34).A significant advantage of the modular setup of the model is that a simple 0D module canbe straightforwardly replaced by the more complicated FE model in Section 2.2. In this setup CircAdapt provides realistic boundary conditions to the FE problem, see Section 2.5.The version of the
CircAdapt model used for all simulations has been published previously [42]and can also be downloaded from the
CircAdapt website ( ). Coupling between PDE and ODE models can be achieved in various ways, fundamentally, theproblem is to find the new state of deformation u n +1 as a function of the pressure p n +1 in a givencavity at time n + 1 . The pressure p n +1 is applied as a Neumann boundary condition at thecavitary surface. This pressure is not known and has to be determined in a way which depends onthe current state of the cavity. Basically, two scenarios have to be considered: (i) when all valvesare closed, the cavity is in an isovolumetric state. That is, the muscle enclosing the cavity maydeform, but the volume has to remain constant. Therefore if active stresses vary over time duringan isovolumetric phase, the pressure p n +1 in the cavity has to vary as well to keep the cavitary8 LV SystemicCirculationPulmonaryCirculation q PDE
RV RA q in p in out p out p out PDE-based 3D solid FE model ODE-based 0D closed-loop model q in p in p out q out PDE
LV AoLA q out p out p RV p LV in out p in p q inin qq PA V RV Figure 3: Schematic of the electrical equivalent circuit showing the coupling of the lumped ODE model to the 3DPDE model. Each of the ventricles (LV, RV) and each of the atria (LA, RA) can be modeled either as a PDE, i.e.,as a 3D elastic body, or as a lumped cavity in the
CircAdapt model. In this case the ventricles are modeled as PDEswhile atria are modeled as lumped cavities. Volume changes of the cavities ˙ V • are driven by flow q • of blood overvalves and outlets. Red colors indicate oxygenated and blue colors de-oxygenated blood. volume constant; (ii) when at least one valve is open or regurging, the cavitary volume is changing.In this case the pressure p n +1 is influenced by the state of the circulatory system. Thus p n +1 hasto be determined in a way that matches mechanical deformation and state of the system. Pressure p n +1 in the cardiovascular system depends on flow and flow rate which are governed by cardiacdeformation and as such the two models are tightly bidirectionally coupled.The simplest approach for the PDE-ODE coupling is to prescribe p n +1 explicitly, e.g., [30, 32].During ejection phases this is achieved by updating cavity volumes and flow based on the currentprediction on the change in the state of deformation under the currently predicted pressure p n +1 . Inthis scenario the pressure boundary condition in each non-linear solver step is modified within eachNewton iteration k . The new prediction p k +1 n +1 is then prescribed explicitly as a Neumann boundarycondition. While this explicit approach is easy to implement and may be incorporated into anexisting FE solver package without difficulty, it may introduce some inaccuracies during ejectionphases and it tends to be unstable during isovolumetric phases. Instabilities stem from the problemof estimating the change in pressure necessary to maintain the volume. Inherently, this requiresto know the pressure-volume ( p - V ) relation of the cavity at this given point in time. However,this knowledge on cavitary elastance is not available and thus iterative estimates are necessary togradually inflate or deflate a cavity to its prescribed volume. As the elastance properties of thecavities are highly non-linear, an overestimation may induce oscillations and an underestimationmay lead to very slow convergence and punitively large numbers of Newton iterations.A more elaborate approach is to treat p n +1 as an unknown [59, 31, 60]. In addition to theequilibrium equations (3–4) this requires one further equation for each cavity c that is modeled as a3D elastic body, see Fig. 3. Hence, we get N cav , the number of PDE cavities, additional equationsof the form V PDE c ( u , t ) − V ODE c ( p c , t ) = 0 , (15)where V PDE c ( u , t ) is the cavity volume computed as the integral over the current surface Γ c,t ,9ee Eq. (14), and V ODE c ( p c , t ) is the cavity volumes as predicted by the CircAdapt model for theintra-cavitary pressure p c , see Section 2.4 and Appendix A.We write p c = [ p c ] c ∈{ lv , rv , la , ra } for the vector of up to ≤ N cav ≤ pressure unknowns.Then, linearization of the variational problem, see Appendix B.2, a Galerkin FE discretization,see Appendix B.3, and a time integration using a generalized- α scheme, see Appendix C, resultin solving the block system to find δu ∈ R N and δp c ∈ R N cav such that K (cid:48) ( u k , p k c ) (cid:18) δuδp c (cid:19) = − K ( u k , p k c ) , K ( u k , p k c ) := (cid:18) R α ( u k , p k c ) R p ( u k , p k c ) (cid:19) , (16)with the updates u k +1 = u k + δu, (17) p k +1 c = p k c + δp c . (18)Here, u k ∈ R N and p k c ∈ R N cav are the solution vectors at the k -th Newton step. The block tangentstiffness matrix K (cid:48) is assembled according to Eqs. (B.20)–(B.24) and Eq. (C.12) and the right handside vector R α according to Eqs. (B.25) and (B.26) and Eq. (C.9). The residual R p which measuresthe accuracy of the current coupling is the discrete version of (15), i.e., R p ( u k , p k c ) := V PDE ( u k ) − V ODE ( p k c ) . (19)The whole procedure to perform the PDE-ODE coupling is summarized in Algorithm 1. Temporal synchronization of chamber contraction.
In the lumped
CircAdapt model contractionin individual chambers is controlled by prescribed trigger events. Based on the measured heartrate, HR, of 103 beats per minute contraction of the left atrium (LA) was triggered at intervalscorresponding to a basic cycle length of / HR = 0 .
585 s . In all other chambers contraction wastriggered by prescribed delays relative to the instant of contraction of the LA. In a hybrid coupledmodel contraction times used in 3D EM (10,11) and in the lumped
CircAdapt model (A.8) must besynchronized accordingly. For this sake a finite-state machine (FSM) was used to control activationcycles in both PDE and 0D chamber models. The FSM provides trigger events to the atrialcavities implemented in the 0D model to initiate contraction, and to the fascicular locations in theventricular cavities of the PDE model by prescribing fascicular activation times t a to the Eikonalequation. Mechanical contraction of the ventricles was initiated then within a prescribed EM delay, t emd . The FSM allows for a fine-grained control of the activation sequence in the combined model,providing trigger events for LA, right atrium (RA), atrio-ventricular node entry and exit by a givenatrio-ventricular (AV) delay, left and right His bundle as well antero-septal, septal and posteriorfascicle in the LV, as well as septal and moderator band fascicle in the RV. The triggers providedby the FSM can be flexibly linked to entities within both 3D and 0D model. The overall conceptfor synchronizing contraction in the coupled model is illustrated in Fig.4. Following [40], initial parameters of passive biomechanics, characterized by the material modelgiven in (8), were taken from [61]. The model’s scaling parameter a was determined then by fittingthe LV model to an empirical Klotz relation [62], using p ed and V ed as input. Active stress model10 lgorithm 1 Coupling of the lumped ODE model to the 3D PDE model Initialize time n = 0 Initial displacement u = 0 Initial cavity pressures p c , = [ p c ] c ∈{ lv , rv , la , ra } at time n = 0 Initialize final time point n max and maximal number of Newton iterations k max Initialize Newton tolerance (cid:15) = 10 − Compute initial cavity volumes V PDE ( u k ) = [ V PDE c ( u k )] c ∈{ lv , rv , la , ra } Run
CircAdapt
ODE system, see Fig. 2, until steady-state is found and get V ODE ( p c , ) =[ V ODE c ( p c )] c ∈{ lv , rv , la , ra } while n < n max do Initialize Newton iterator: k = 0 Initial guesses for Newton: u = u n , p c = p c ,n while k < k max do Assemble block matrix (cid:46)
Eqs. (B.20)–(B.24) and Eq. (C.12)and right hand side. (cid:46)
Eqs. (B.25)–(B.28) and (C.9)–(C.11)
Solve linearized system for δu and δp c (cid:46) Eq. (B.17)
Update displacement u k +1 = u k + δu and cavity pressures p k +1 c = p k c + δp c Update cavity volumes V PDE ( u k +1 ) (cid:46) Eq. (14)
Update ODE system and get V ODE ( p k +1 c ) (cid:46) Fig. 2
Convergence test: if (cid:107) R α ( u k +1 , p k +1 c ) (cid:107) L < (cid:15) and (cid:107) R p ( u k +1 , p k +1 c ) (cid:107) ∞ < (cid:15) then Solution at time n + 1 : u n +1 = u k +1 , p c ,n +1 = p k +1 c break (cid:46) Newton converged else k = k + 1 (cid:46) Next Newton step end if end while n = n + 1 (cid:46) Next time step end while igure 4: EM activation of the coupled 3D–0D model is steered by a FSM that provides triggers for electricalactivation of the 3D EM model and for mechanical activation of the 0D lumped atrial cavities. The sino-atrial nodeclock (SA) activates the LA at a prescribed cycle length, SA CL , starting at time SA t . The RA initiates contractionwith a delay of AA delay after the LA. The atrial entrance into the AV node activates at AV a and the ventricular exitof the AV node connected to the His bundle at AV v . Fascicles in the LV are activated then relative to the LV triggerto initiate electrical propagation in the EP model. Similarly, the RV is activated with an interventricular delay ofVV d before (VV d <0) or after VV d > the LV. parameters τ c , S peak , τ r and duration of the force transient, t dur were determined as describedpreviously [23]. Initial values and parameters for the CircAdapt model were chosen following [63].EM activity was simulated over 20 beats at the given heart rate. Execution of the semi-implicit3D-0D model was sped up by limiting the number of Newton steps to k max = 1 for the initialseries of heart beats that were simulated to stabilize the coupled 3D-0D model to a limit cycle.Finally, after arriving at a stable limit cycle two further beats were simulated using a fully con-verging Newton method with k max = 20 and an relative (cid:96) norm error reduction of the residual of (cid:15) = 10 − . To replicate the observed stroke volume (SV) and LV peak pressure ( ˆ p LV ) in the leftventricle, input parameters of the active stress model as well as CircAdapt model parameters wereiteratively adjusted. The final parameterized model beating at 103 bpm produced a cardiac outputof ≈ . / min with a SV of
21 mL . The coupled 3D-0D model was subjected to a thorough physiological validation by evaluatingits transient response to alterations in loading conditions and contractile state. The model underbaseline condition was used as a reference working point relative to which the effect of perturbations12n loading and contractility was compared. In absence of specific experimental data under standardprotocols for assessing of systolic and diastolic properties of the ventricles based on pV analysis [64],the validation steps performed were of qualitative nature where the model’s ability to consistentlyreplicate known cardiovascular physiology was gauged. For all perturbations in preload, afterloador contractility, two points in time were considered, the immediate acute response after perturbingthe system and the new approximate limit cycle reached after 8 beats. For baseline and each limitcycle, the ESPVR of the LV was interrogated by imposing additional step changes in afterload. Forthis sake, data points p es and V es were determined in the pV loops at the instant of end-systole, asdetermined by the cessation of flow out of the LV. Linear regression was used then to determineend-systolic elastance, E es , as the slope of the regression curve, and the volume intercept, V d , of theESPVR. Preload, afterload and contractility step changes were implemented by varying the cross-sectional area of the pulmonary veins, the systemic vascular resistance, R sys , and the active peakstress S peak generated by the myofilament model, respectively. Overall pump function was assessedby constructing a pump function graph from data gathered under all protocols, including afterloadvariations over a wider range, between E a ≈ by setting the system resistance R sys ≈ and E a ≈ ∞ , by closing the aortic valve, to obtain data points under extreme conditions correspondingto the LV beating in absence of external loading and under isovolumetric conditions, respectively. After discretization, at each Newton–Raphson step the block system (16) has to be solved. Inthat regard, we applied a Schur complement approach, see Appendix E, to cast the problem in apure displacement formulation and thus can reuse solver methods established previously [39]. Inbrief, we used the generalized minimal residual method (GMRES) with an relative error reductionof (cid:15) = 10 − . Efficient preconditioning was based on the library PETSc [65] and the incorporatedsolver suite hypre/BoomerAMG [66]. For the time integration we used a generalized- α scheme,see Appendix C, with spectral radius ρ ∞ = 0 and damping parameters β mass = 0 . − , β stiff =0 . .We implemented the coupling scheme in the FE framework Cardiac Arrhythmia Research Pack-age (CARPentry) [67, 55], built upon the openCARP EP framework ( ).Based on the MATLAB code available on the CircAdapt website ( )a C ++ circulatory system module was implemented into CARPentry to achieve a computationallyefficient and strongly scaling numerical scheme that allows fast simulation cycles.
3. Results
The coupled 3D-0D model was fit to approximate the available experimental data on ˆ p and strokevolume under baseline conditions. Electrical activation was driven by a tri- and bi-fascicular modelin LV and RV, respectively (see Fig. 5(A)). Conduction velocities were chosen for the given activationpattern to obtain a total ventricular activation time of . , compatible with the observedQRS duration of the ECG. Mechanical boundary conditions were set to limit radial contraction ofthe model, thus leading to a heart beat where ejection was mediated largely by atrio-ventricularplane displacement, i.e. long-axis shortening of the ventricles, and myocardial wall thickening. Theresulting end-disastolic and end-systolic configuration of the model is shown in Fig. 5(B). Trainsof 20 heart beats were simulated in the coupled model to arrive at an approximate stable limitcycle, as verified by inspecting the slope of the envelope of key hemodynamic state variables (see13ig. 5(C)). Corresponding hemodynamic data on pressure, volume and flows for all four chambersalong with the corresponding pV loops over the last two beats are shown in Fig. 5(B). Parametervalues of the baseline model are given in tables Tabs. A.3 and A.4 and Tab. 1 for the 0D CircAdapt and 3D PDE model, respectively. (A) P r e ss u r e ( mm H g ) V o l u m e ( m l ) . . . . . . Time (s) − . . . F l o w (cid:0) l / s (cid:1) (C)
20 30 40 50 60 70 80
Volume (ml) P r e ss u r e ( mm H g ) SVESPVR EDPVRˆ p LV LA (ODE)RA (ODE)LV (PDE)RV (PDE)
Time (s) P r e ss u r e ( mm H g ) (B) Figure 5: Model parameterization under baseline conditions. (A) Shown are ventricular sinus activation sequenceinduced by three LV and two RV fascicles (top panels) and the mechanical end-diastolic (red) and end-systolic (blue)configuration. Note the minor change in epicardial shape due to the pericardial boundary conditions. (B) Simulatedpressure traces in LV (green) and RV (blue) are shown for the entire pacing protocol using a train of 20 beats..Envelopes (dotted traces) indicate that an approximate limit cycle was reached after 3 beats. (C) Left panels showtime traces of pressure p , flow q and volume V in lumped 0D atrial cavities and PDE-based ventricular cavities for thelast two beats of the limit cycle pacing protocol. All variables traverse the state space along limit cycle trajectories.Right panel shows pressure-volume loops in all four cavities. For PDE-based ventricular cavities EDPVR and ESPVRare indicated. Experimental data on peak LV pressure and stroke volume used for fitting are indicated (dashed lines). The impact of the relatively coarse spatial resolution ( ∼ . and ∼ . for LV andRV, respectively) used was evaluated first by repeating the baseline limit cycle protocol using ahigher resolution mesh ( ∼ . and ∼ . for LV and RV, respectively). Both simulationsused the exact same set of parameters and initial state vectors. With regard to pV behavior thatis governed by the global motion and deformation of the ventricles ventricular end-diastolic andend-systolic configurations were compared(see Fig. 6). Observed discrepancies between coarse andhigher resolution model were marginal and clearly within the limits of experimental data uncertainty,suggesting that the computationally efficient coarse model is suitable for performing a physiologicalvalidation study. The resonse of the coupled 3D-0D system to changes in afterload was probed by altering R sys in the range of ± around its nominal value of − ms − These maneuvers alter14he slope of the arterial elastance curve, E a , pivoting E a around the point V ed and p = 0 in the pV diagramme. Initial response immediately after step changes in afterload and the new limit cycleare shown in Fig. 7(A) and (D), respectively.The response to altering LV preload was then probed by stepwise reducing blood flow from thelungs into the LA. Under such a walk down protocol the E a curve is shifted to the left towardssmaller end-diastolic volumes V ed , without altering its slope, i.e. E a = p es / SV ≈ const, leadingto lower p es . Stroke volumes under this protocol are assumed to gradually reduce due to theFrank-Starling mechanism, mediated by the length-dependence of active stress generation, S a ( λ ) .Transient and limit cycle response under the preload perturbation protocol are shown in Fig. 7(B)and (E), respectively. As contractile properties remained unchanged, the same slope E es of theESPVR was obtained as before under step changes in afterload (compare estimated E es betweenFig. 7(D) and (E)).The effect of step changes in contractility was probed by altering peak active stresses S peak in the LV by ±
20 % around the LV nominal value of
100 kPa
This maneuver steepened/flattenedESPVR (see Fig. 7(C) and (F)). In the initial response to a step change in contractile force V es and p es were affected, with V ed remaining constant, this led to a change in stroke volume and anapparent change in arterial elastance estimated by E a ≈ p es / SV, with ≈ − . / + 36 . relativeto baseline. However, in the limit cycle response after readjustment of preload in all chambers, all E a curves had the same slope and were only shifted according to the new working V ed for the givencontractile state. Transient pV loops under this protocol are shown in Fig. 7(C) and (F).A PFG was constructed by combining data from all tested protocols with additional samplingof data within more extremal ranges of ventriculo-arterial coupling. The models’ PFG is in agree-ment with known cardiovascular physiology, as shown in Fig. 8. Keeping contactility and preloadconstant, the PFG with mean ventricular pressure (MVP) as a function of flow or stroke vol-ume can be approximated by a quadratic function, MVP ( q ) = ˆ P iso (cid:0) − ( q/q mx ) (cid:1) , with ˆ P iso and q mx being maximum mean ventricular pressure and maximum flow under isometric and unloadedconditions, respectively. Increasing preload shifts the PFG towards higher flows and pressures. In-creasing/decreasing contractility pivots the PFG, leading to a steeper/flatter slope ∆ MVP /∆q ofthe PFG.
Computational times for a single heart beat of the lower and higher resolution baseline modelsare given in Tab. 2. Simulations were performed on the Vienna Scientific Cluster (VSC4) and (A) (B)
50 60 70
Volume (ml) P r e ss u r e ( mm H g ) (C) Figure 6: Differences between the coarse (wireframe, red solid line) and higher resolution (solid, blue solid line)model are shown for (A) end-diasolic and (B) end-systolic configuration. (C) Dynamic behavior over the limit cycleprotocol was comparable with minor difference in V ed , V es and ˆ p LV . arameter Value Unit Description Passive biomechanics ρ . / m tissue density κ kPa bulk modulus a kPa stiffness scaling b ff b ss b nn b fs b fn b ns Active biomechanics λ ms minimum fiber stretch V m , Thresh -60.0 mV membrane potential threshold t emd ms EM delay S peak
100 (lv), 80 (rv) kPa peak isometric tension t dur ms duration of active contraction τ c ms baseline time constant of contraction ld ld up ms length-dependence of upstroke time τ r ms time constant of relaxation Electrophysiology t cycle .
585 s cycle time ( = 1 / heartrate )AA delay 20.0 ms inter-atrial conduction delayAV delay 100.0 ms atrioventricular conduction delayVV delay 0.0 ms inter-ventricular conduction delay ( v f , v s , v n ) (0.6, 0.4, 0.2) m / s conduction velocities ( v f , v s , v n ) (1.2, 0.8, 0.4) m / s fast conducting layer velocities ( g f , g s , g n ) (0.44, 0.54, 0.54) m / s conductivities in LV and RV ( g f , g s , g n ) (0.65, 6.27, 6.27) m / s conductivities in the fast conducting layers β cm − membrane surface-to-volume ratio C m F / cm membrane capacitanceTable 1: Input parameters for the 3D PDE model of the left (lv) and right (rv) ventricle. Adjusted to matchsubject-specific data. Volume (ml) P r e ss u r e ( mm H g ) V = 35 . E es = 2 . mmHgml V d = 13 . E a E a , bsl = 3 . E a , inc = 4 . E a , dec = 3 . ↑ afterload ↓ (A)
20 30 40 50 60 70 80
Volume (ml) P r e ss u r e ( mm H g ) V = 35 . E es = 2 . mmHgml V d = 13 . E a , bsl = 3 . E a , inc = 3 . E a , dec = 4 . ↑ preload ↓ (B)
35 40 45 50 55 60 65
Volume (ml) P r e ss u r e ( mm H g ) baselinecontractility ↑ contractility ↓ (C)
20 30 40 50 60 70
Volume (ml) P r e ss u r e ( mm H g ) V = 35 . E es = 2 . mmHgml V d = 13 . E a E a , bsl = 3 . E a , inc = 6 . E a , dec = 2 . ↑ afterload ↓ (D)
10 20 30 40 50 60 70 80
Volume (ml) P r e ss u r e ( mm H g ) V = 35 . E es = 2 . mmHgml V d = 12 . E a , bsl = 3 . E a , inc = 3 . E a , dec = 3 . ↑ preload ↓ (E)
35 40 45 50 55 60 65 70 75
Volume (ml) P r e ss u r e ( mm H g ) baselinecontractility ↑ contractility ↓ (F) Figure 7: Initial response (A-C) and 4 cycles (D-F) after applying a step change in loading conditions and contractility. (A)
Altering afterload by increasing/decreasing the systemic vascular resistance, R sys pivots arterial elastance E a curve. Endsystolic elastance, E es and intercept V d characterizing the ESPVR was determined by linear regression ofend-systolic data points V es and p es , marked by solid circles. (B) Increasing/decreasing preload shifts E a curve andincreases/decreases stroke volume via the Starling mechanism, mediated by the length-dependence of the active stressmodel. Determination of ESPVR was consistent with afterload protocol. (C) Increasing/decreasing contractilityincreases/decreases stroke volume, LV peak pressure and p es . For each contractile state afterload was also perturbedto determine end-systolic elastance E es and V d . we distinguish between solver-time t s which is the accummulated GMRES solver time over allloading/time steps; and assembly-time t a which is the time needed for setup of boundary conditionsand assembling of matrices and vectors of the linearized system (16). In total, for a full simulationwith loading, 18 initialization beats, and 2 final beats with a fully converging Newton method thecomputational costs were .
38 s for the lower resolution model on 24 cores and .
73 s for thehigher resolution model on 256 cores. Here, in addition to GMRES solver and assembly times, alsoEP, postprocessing, ODE, and input-output times are taken into account. It is worth noting thatthe ODE solver alone is very efficient: for a simulation of 103 beats (i.e. 1 minute with a cyclelength of .
585 s ) computational costs were approximately on one core of a desktop computer.Hence, ODE times carry almost no weight in the coupled 3D-0D model.17
Volume (ml) P r e ss u r e ( mm H g ) baselineafterload ↑ afterload ↓ (A)
10 15 20 25 30 35 40
SV (ml) M V P ( mm H g ) baselinecontractilitypreload (B) Figure 8: PFG reconstructed from loading protocols. (A)
Starting from baseline conditions, afterload was variedbetween unloaded conditions, E a ≈ , and isometric conditions, E a ≈ ∞ . (B) PFG constructed from afterloadvariations with constant preload and contractility (solid circles), with increased preload (solid squares) shifting thePFG up and left towards higher flow and pressure. For both cases contractility was also perturbed which pivots thePFG (empty circles and squares, respectively), leading to a steeper/flatter slope MVP /∆SV for increased/decreasedcontractility.
4. Discussion
In this study, we report on the development of a flexible monolithic 3D solid - 0D fluid couplingmethod that allows to mix 3D PDE and non-linear 0D EM representations of cardiac cavities. Thehybrid 3D-0D model of the heart was coupled to a 0D closed-loop model of the cardiovascular sys-tem, where all 0D components were based on the
CircAdapt model [42] The combined model canbe set up to represent the LV only, both ventricles, both atria or all four cavities as 3D PDE modeland all other elements of the cardiovascular system as 0D models based on
CircAdapt . In this studyfeasibility of this approach is demonstrated by coupling a 3D PDE model of bi-ventricular EM to0D model representations of atrial EM function and circulation based on the
CircAdapt model.The combined model was parameterized under baseline conditions and subjected to rigorous phys-
Model cores ¯ h lv / ¯ h rv Elems Nodes t s , /t a , t s , c /t a , c t s , ld /t a , ld [-] [ mm ] [-] [-] [ s ] [ s ] [ s ]coarse 24 . / . . / . . / . . / . fine 256 . / . . / . . / . . / . Table 2: Summary of numerical metrics for coarse and fine model. Given are the number of compute cores usedon VSC4; the average spatial resolution in LV and RV, ¯ h lv and ¯ h rv ; the number of elements and nodes spanningthe mesh; as well as solver, assembly, and total times for a single Newton iteration ( t s , , t a , ) and a fully convergedNewton solution ( t s , c , t a , c ). Timings refer to a single heart beat lasting .
585 s at a time step size of . In additionthe cummulated solver ( t s , ld ) and assembly ( t a , ld ) times for the loading phase using 32 load steps are presented. pV analyses of the hemodynamic model output showed close agreement with established knowledgeon cardiovascular physiology. The underlying numerical scheme is represented in detail, includinga comprehensive mathematical representation of the CircAdapt model. Robustness – in terms ofstability and convergence properties – and computational efficiency – in terms of execution times– are demonstrated. These features combined render advanced EM modeling applications feasible.The model facilitates the efficient and robust exploration of parameter spaces over prolonged obser-vation periods which is pivotal for personalizing models to closely match observations. Moreover,the model can be trusted to provide predictions of the acute transient response to interventions ortherapies altering loading conditions and contractility that are valid within a commonly acceptedphysiological reference frame.
Predictive modeling applications critically rely on the ability of models to encapsulate the mostrelevant mechanisms governing the cardiovascular response to a given intervention that alters load-ing conditions or contractility. In a closed-loop cardiovascular system as represented by
CircAdapt ,isolated changes to a single parameter entail transient adaption processes in the system as a whole.These were studied first in the standalone 0D
CircAdapt model and then replicated, using the sameprotocols, in the 3D-0D coupled model. Qualitatively, the same trends were observed in the coupledmodel (data not shown), with minor quantitative deviations attributable to differences in the activestress generation between the 0D cavity model used in
CircAdapt and the 3D PDE models replacingthese in the coupled 3D-0D system.Validation aimed at replicating, overall, known well established behaviors and not at a 1:1validation against experimental data. For this sake, experimental standard protocols for alter-ing afterload, preload and contractility were applied to the stabilized baseline model to study itsresponse. Two scenarios were analyzed, the initial acute response to a step change in a single pa-rameter – afterload, preload or contractility – where effects on the other unaltered parameters wereminimal (see Fig. 7(A)-(C)), and, the limit cycle response observed after a number of beats wheretransients have largely subsided, but indirect effects led to alteration of other parameters due tothe systemic inter-dependencies between these (see Fig. 7(D)-(F)).Afterload was altered by varying the systemic resistance R sys and, to mimick more extremeconditions closer to isometric contraction, the resistance of the aortic valve. Increasing/decreasingafterload is reflected in pivoting the arterial elastance curve E a around the point of end-diastolicvolume and zero pressure. This behavior is illustrated in Fig. 7 and for the more extensive protocolsused in constructing the PFG in Fig. 8(A). E a was only marginally affected in the initial response,but more significant changes were witnessed after stabilization to a new limit cycle. As expected,changes in slope of E a were proportional to changes in R sys since, in absence of any regulatorymechanisms, heart rate remained unaltered. Thus, E a ≈ MAP / ( SV · HR ) ≈ R sys holds. Thephenomenological active stress model accounted for length-dependent tension, but not for velocitydependence. Thus, afterload effects on the velocity of fiber shortening were ignored.Altering V ed by changing the cross section of the pulmonary vein orificies to increase/reducepreload increased/reduced SV. This was due to the Frank-Starling mechanism, as represented bythe length-tension relationship of the active stress model. In the immediate response to a stepchange in preload the LV emptied to almost the same V es with only minor deviations due to19hanges in arterial pressure induced by the change in SV. After transients subsided, ESPVR andarterial elastance were the same as under baseline conditions, with E a being shifted according to thechanged V ed . In the new limit cycle p es and MAP were increased, but changes in SV were rathersmall. This could be attributed to a rather flat slope of the Frank-Starling curve SV ( p ed ) (notshown) that is related to the significant heterogeneity in fiber stretch in end-diastolic state. Suchheterogeneity is inevitable in biventricular EM models that do not account for residual strains in theunpressurized configuration. As fiber stretch in the reference configuation with p = 0 is assumed λ = 1 , increasing the filling pressure to p ed leads to a significant spread in λ around its mean.Thus, while mean λ increased with p ed the increasing spread of λ led to low fiber stretch λ < λ invarious regions, particularly in antero-septal and postero-septal segments of the LV. These regionscontributed increasingly less to contraction with increasing p ed due to the length-dependence of S a ( λ ) ≈ , thus, leading to an overall attenuation of the Frank-Starling effect.Altering contractility by increasing/reducing the peak active stress S a led to an increase/decreasein SV and systolic pressures in terms of ˆ p and p es . Correspondingly, the ESPVR, as sampled byperturbing LV afterload, was steepened/flattened as expected (see Fig. 7(C)). In the initial responsethe slope of E a was affected, but after stabilization E a was the same for all contractile states. Forthe given parameterization ventriculo-arterial coupling E es /E a fell outside the optimal range of1–2, where external work is maximized around E es /E a ≈ and optimal efficiency is achieved at E es /E a ≈ . Since ˆ p and SV and as such E a ≈ p es / SV were used to parameterize the model theresulting slope of the ESPVR was to flat for the LV to operate within the optimal range. Reasonsare multifactorial and include the absence of velocity-dependence or, potentially, also the role ofmechanical boundary conditions. The main culprit to blame for the limited slope of E es is theimpairment of the Frank-Starling mechanisms due to fiber stretch heterogeneity while a uniformfiber-stretch is assumed to be in place at end-diastolic state [68, 69, 70, 71].The constructed PFG was qualitatively in keeping with physiological expectations (see Fig. 8).Data points obtained from varying afterload between almost unloaded and isometric conditionsagreed well with an assumed quadratic relation between MVP and flow or SV. Increasing preloadshifted the PFG up and left towards higher flows and pressures, with the opposite trend beingobserved for decreasing preload. Altering contractility rotated the PFG. The computational cost imposed by high resolution EM models requires efficient numericalsolvers. Strong scaling characteristics of our numerical framework was reported in detail previ-ously [39, 72]. Also for our coupled 3D-0D model we achieved favorable strong scaling propertiesup to cores depending on the chosen resolution.The compute times reported in Section 3.4 indicate that setup and assembly time were thedominating factors during the initial passive inflation (loading) phase while for the subsequentcoupled 3D-0D EM simulations of a heart beat solver time was the predominant part of total CPUtime. This is due to the Schur complement, see Appendix E, that needs to be solved during the3D-0D active EM phase. This involved – in our case of a bi-ventricular model comprising two PDEcavities – three applications of the GMRES solver while matrices were only assembled once perNewton step.Further, significant savings in compute time could be achieved using a semi-implicit approachuntil a limit cycle was reached. A heartbeat using a fully converging Newton was about five timesmore expensive compared to a heartbeat in the initialization phase. As the deviations of the semi-implicit method from the implicit scheme were negligible and the final two beats of the limit cycle20rotocol were then computed using a fully converging Newton method this gain in performance hadno quantiative impact on the outcome of the simulation.Overall, the total simulation time incurring for 20 heart beats of the coarse model was well below
90 min on 24 cores. Similar performance was achieved on a standard desktop computer (AMD RyzenThreadripper 2990X), demonstrating that realistic multi-beat simulations of the presented 3D-0Dcardiac models deliver sufficient performance for advanced physiological simulation scenarios, evenon a small number of cores. This is of paramount importance for future parameterization studieswhere numerous simulations have to be carried out to personalize the model to patient data. Witharound 2 and a half hours on 256 cores for 20 beats also the higher resolution model could beexecuted within a tractable time frame.
Efforts to parameterize the combined model were limited, only a small subset of available datawere used for model fitting. More comprehensive and efficient procedures are required to furtherenhance compatibility of the high dimensional combined 3D-0D model with all available observa-tions. Building on strategies for fitting the standalone
CircAdapt model to hemodynamic data [], ahybrid parameterization approach appears a pragmatic solution where the 0D model is fitted first toobserved data and in a subsequent parameterization step the 3D cavities are fit to the pV character-istics of the corresponding 0D cavities. However, for the sake of demonstrating overall compatibilitywith known cardiovascular physiology a high fidelity match with experimental data, while desirable,is not crucial. Future more advanced applications that attempt to predict therapeutic responsesin a patient-specific manner will critically depend, beyond a comprehensive representation of themost relevant therapy mechanisms, on the fidelity of personalization.While the framework used in this study is, in principle, complete important factors such as veloc-ity dependence of active force generation or intrinsic physiological details of excitation-contractioncoupling remained unaccounted for. Active force generation in our model is not driven by Calciumtransients of the cellular dynamics model. Rather, this is triggered within a prescribed fixed EMdelay by the upstroke of the action potential and, as such, does not take into account the spa-tial variability of EM delays [73]. Further, length-dependence of the rate of contraction as wellas the peak stresses generated were taken into account by the active stress model. This providedthe basis for representing the Frank-Starling mechanism at the organ scale, as demonstrated bythe protocol for altering preload and the corresponding response in SV. However, inflation of theunpressurized configuation to p ed inevitable introduces a significant spatial heterogeneity in fiberstretch, thus limiting the expected Frank-Starling effect at the organ scale. This is in contrast tothe common assumption of homogeneous fiber stretch in the end-diastolic state [68, 69, 70, 71].This was reflected in a rather flat slope of the Starling curve SV ( p ed ) (not shown).Finally, for the sake of saving computational costs, a coarse spatial resolution was used fordiscretizing the bi-ventricular model. Thus, models were lightweight enough to carry out the largernumber of simulations that were needed for parameterization, the determination of a limit cycleas well as the fine grained testing of physiological maneuvers. The use of such coarse spatialresolutions introduce inaccuracies with regard to fibers and sheet arrangements which are definedon a per element basis. As parts of the model such as the RV wall were composed only of two elementlayers, transmural fiber rotation was esentially reduced to two fiber families. These were mostlyaligned with the endocardial and epicardial fiber orientations as prescribed on a per rule basis (seeFig. 1). Thus, the model’s predictions of motion, strains and stresses deviated quantitatively fromhigher resolution models. However, comparing between the models of different resolution revealed21hat quantitative differences were minor and, qualitatively, models showed essentially the exactsame behavior. Most differences stemmed from differences in stiffness of the simple P1–P0 elementtypes used which tended to be more compliant for higher resolutions. Reference simulations carriedout at higher spatial resolutions of ≈ indicated that the errors in simulated hemodynamicoutcomes were small enough to be considered negligible when put in context to the observationaland residual uncertainty clinical or experimental data are afflicted with.
5. Conclusions
This study reports on a flexible monolithic 3D solid - 0D fluid coupling method for integratedmodels of cardiac EM and cardiovascular hemodynamics. Feasibility of the approach is demon-strated by coupling a 3D PDE model of bi-ventricular EM to 0D model representations of atrialEM and circulation based on the
CircAdapt model. The coupled 3D-0D model is shown to berobust, computational efficient and able to correctly replicate known physiological behaviors in re-sponse to experimental protocols for assessing systolic and diastolic ventricular properties based on pV analysis. These features combined render advanced EM modeling applications feasible. Themodel facilitates the efficient and robust exploration of parameter spaces over prolonged observationperiods which is pivotal for personalizing models to closely match observations. Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relation-ships that could have appeared to influence the work reported in this paper.
Acknowledgements
This project has received funding from the European Union’s Horizon 2020 research and inno-vation programme under the Marie Skłodowska–Curie Action H2020-MSCA-IF-2016 InsiliCardio,GA No. 750835 and under the ERA-NET co-fund action No. 680969 (ERA-CVD SICVALVES)funded by the Austrian Science Fund (FWF), Grant I 4652-B to CMA. Additionally, the researchwas supported by the Grants F3210-N18 and I2760-B30 from the Austrian Science Fund (FWF)and a BioTechMed Graz flagship award “ILearnHeart” to GP.22 ppendix A.
CircAdapt equations summary
Appendix A.1.
CircAdapt
Tube Module
The tube module represents the entrance of a compliant blood vessel capable of propagating apressure-flow wave component added to a constant flow, see [25]. Vessels directly attached to theheart, aorta (ao), arteria pulmonalis (ap), venae cavae (vc), and venae pulmonales (vp) are modeledin a similar fashion in
CircAdapt and for the whole section we define the iterator t ∈ { ao , ap , vc , vp } .The current lumen cross sectional area is computed by A t = V t l t , (A.1)with V t the cavity volume and l t the length of the vessel segment.Using a model of a tube with a fibrous wall, see [74], gives the average extension, λ t , of thefibers in the wall by λ t = (cid:0) V t /V wall t (cid:1) / = (cid:0) A t /A wall t (cid:1) / , with A wall t the wall cross-sectional area and V wall t = A wall t l t the wall volume. Cavity pressuredepends on λ t p t = σ t ( λ t ) λ − t , with the mean Cauchy fiber stress σ t that is modeled by the constitutive equation σ t ( λ t ) = σ ref t · (cid:0) λ t /λ ref t (cid:1) k t , see Arts et al. [74]. Here, k t a stiffness exponent; λ ref t = (cid:0) V ref t /V wall t (cid:1) / and σ ref t = p ref t (cid:0) λ ref t (cid:1) are the fiber extension and fiber state at normal physiological reference state, respectively; and p ref t is the reference tube pressure. By combining above results the current tube pressure is computedas p t = p ref t (cid:18) λλ ref t (cid:19) k λ − t (cid:0) λ ref t (cid:1) = p ref t (cid:18) λ t λ ref t (cid:19) k − = p ref t (cid:18) V wall t + 2 V t V wall t + 2 V ref t (cid:19) k − = p ref t (cid:18) A wall t + 2 A t A wall t + 2 A ref t (cid:19) k − , (A.2)with A ref t the initial cross sectional area and V ref t = A ref t l t the initial vessel volume. The complianceis C t = d p t d V t = dd V t p ref t (cid:18) V wall t + 2 V t V wall t + 2 V ref t (cid:19) k − = p ref t k − (cid:18) V wall t + 2 V t V wall t + 2 V ref t (cid:19) ( k − − ) 2 V wall t + 2 V ref t = 2 p t ( k − (cid:0) V wall t + 2 V t (cid:1) = 2 p t ( k − l t (cid:0) A wall t + 2 A t (cid:1) . Z t compliance C t , and inertance I t Z t = I t C t = ρ b l t A t C t (A.3)yields Z t = (cid:115) ρ b p t l t ( k − V t (cid:0) V wall t + 2 V t (cid:1) = (cid:115) ρ b p t ( k − A t (cid:0) A wall t + 2 A t (cid:1) (A.4)with ρ b the blood density. Appendix A.2. Sarcomere mechanics
In the following a sarcomere contraction model is described that is based on a modified Hillmodel, see [58, 42], for all tissue patches in the wall of the cavity c ∈ { lv , rv , sv , la , ra } . Naturalstrain E fib c of the myofiber is estimated as E fib c = ln (cid:18) L s c L s , ref (cid:19) (A.5)and from this the total sarcomere length L s c can be computed as L s c = L s , ref exp (cid:0) E fib c (cid:1) , (A.6)with L s , ref c a constant describing the reference sarcomere length. The sarcomere is supposed tobe made up of a contractile element of length L cont c in series with an elastic element of length L elast c = L s c − L cont c . Sarcomere active stress
Sarcomere contracting length L cont c varies over time according to ˙ L cont c := d L cont c d t = v max (cid:20) L s c − L cont c L elast , iso − (cid:21) , (A.7)where the constant L elast , iso is the length of the series elastic element during isovolumetric contrac-tion and the constant v max is the maximum velocity of contraction.The governing equation for the contractility C c , describing the density of cross bridge formationwithin the fibers of the patch, is ˙ C c := d C c d t = f rise c ( t ) C s c (cid:0) L cont c (cid:1) − f decay c ( t ) C c . (A.8)In Eq. (A.8) sarcomere contractility rises according to the function f rise c , a phenomenological rep-resentation of the rate of cross bridge formation within the patch, f rise c ( t ) = 1 t rise c . x c (8 − x c ) exp( − x c ) ,x c ( t ) = min (cid:18) , max (cid:18) , t − t act c t rise c (cid:19)(cid:19) , t act c , and the rising time t rise c = 0 . τ R t act , ref c . Here, τ R a constant and t act , ref c is the reference duration of contraction for initial fiber length.Sarcomere contractility in (A.8) decays according to the function f decay c f decay c ( t ) = 12 t decay c (cid:104) (cid:16) sign( y c ) min (cid:16) π , | y c | (cid:17)(cid:17)(cid:105) ,y c ( t ) = t − t act c − t act , dur c t decay c , depending on the decay time t decay c = 0 . τ D t act , ref c , with τ D a constant and t act , dur c is the duration of contraction of the fiber that lengthens withsarcomere length t act , dur c = (0 .
65 + 1 . L norm c ) t act , ref c . Here, L norm c is the normalized sarcomere length for active contraction L norm c = max (cid:0) . , L cont c /L act0 , ref − (cid:1) , where L act0 , ref is the zero active stress sarcomere length. C s c in (A.8) describes the increase in cross bridge formation with intrinsic sarcomere length dueto an increase in available binding sites, C s c (cid:0) L cont c (cid:1) = tanh (cid:16) . ∗ . L norm c ) (cid:17) . Contractility C c (A.8) and sarcomere contracting length L cont c (A.7) are used to compute the activelygenerated fiber stress σ fib , act c = L act0 , ref σ act , max C c L norm c L s c − L cont c L elast , iso , (A.9)with constants L act0 , ref , σ act , max , L elast , iso , see Tab. A.4. Sarcomere passive stress
Passive stress σ fib , pas c is considered to contain two components, σ fib , pas c = σ fib , tit c + σ fib , ecm c , (A.10)first the stress arising from cellular structures such as titin, a highly abundant structural protein ofthe sarcomere, anchoring to the Z-disc, σ fib , tit c , and second the stress arising from the extracellular25atrix (ECM), σ fib , ecm c . Both depend on the passive fiber stretch which is computed as λ pas c = L s c L pas0 , ref , where L pas0 , ref is sarcomere length with zero passive stress and L s c the total sarcomere length, seeabove. Using that we compute σ fib , tit c = 0 . σ act , max (cid:16) [ λ pas c ] k tit − (cid:17) , with σ act , max the maximal isometric stress and the constant exponent k tit = 2 L s , ref d L s , pas . The ECM is modeled as being stiffer than the myocyte contribution using σ fib , ecm c = 0 . σ pas , max (cid:16) ( λ pas c ) − (cid:17) , where σ pas , max is an empirical parameter. Sarcomere total stress
Total myofiber stress σ fib c is the sum of an active (A.9) and a passive (A.10) stress component σ fib c = σ fib , act c + σ fib , pas c . (A.11)Sarcomere stiffness κ fib c is now computed as the derivative of total fiber stress (A.11) with respectto fiber strain (A.5) κ fib c = ∂σ fib c ∂E fib c = ∂σ fib , act c ∂E fib c + ∂σ fib , pas c ∂E fib c , (A.12)with ∂σ fib , act c ∂E fib c = L act0 , ref σ act , max C c L norm c L s c L elast , iso ,∂σ fib , pas c ∂E fib c = 0 . k tit σ act , max c ( λ pas c ) k tit + 0 . ∗ σ pas , max ( λ pas c ) . Appendix A.3.
CircAdapt
Chamber Module
An actively contracting chamber c ∈ { lv , rv , la , ra } is modeled using the state variables volume V c , length of the contractile element of the sarcomere L cont c (A.7), and contractility C c (A.8).Volume changes driven by inflow and outflow of blood induce changes in midwall volume V mid c andarea A mid c . Sphere mechanics
Note that in
CircAdapt ventricles are usually modeled using the TriSeg formulation, see Ap-pendix A.4. If TriSeg is turned on, the calculations in this chapter are only used for the atria whileventricular values are computed as in Appendix A.4.26idwall volume V mid c is estimated as V mid c = V c + 12 V wall c , (A.13)where V wall c is constant wall volume. If not set to a specific value the wall volume is estimated byextruding the sphere enclosing the cavity volume V c by a constant wall thickness h wall c , see Tab. A.3.Chambers are modeled as closed spheres, thus, the following equations result from volume andsurface formulas for spheres C mid c = (cid:18) π V mid c (cid:19) / , (A.14) A mid , tot c = 4 π ( C mid c ) , (A.15) A mid c = A mid , tot c − A mid , dead c , (A.16)where C mid c is midwall curvature, i.e., the inverse of radius; and A mid , dead c is non-contractile area,i.e., valve openings and orifices. Update fiber strain
Natural fiber strain E fib c is calculated by E fib c = 12 ln (cid:18) A mid c A mid , ref c (cid:19) (A.17)with A mid , ref c the surface area in the reference state, see [42]. Note that this updated fiber strain isused in the place of (A.5) to update values in the sarcomere module Appendix A.2.Cross-sectional area A c of chambers is estimated as A c = V c + 0 . V wall c l c , (A.18) l c = 2 (cid:0) V mid c (cid:1) / , with l c the long-axis length of the cavity.The characteristic wave impedance Z c is approximated according to (A.3), see also [25], and byapplying the chain rule Z c = 15 A c (cid:113) ρ b l c | κ mid c | , (A.19)with the sheet stiffness κ mid c = ∂T mid c ∂A mid c = V wall c A mid c ) (cid:18) ∂σ fib c ∂E fib c − σ fib c (cid:19) = V wall c A mid c ) (cid:0) κ fib c − σ fib c (cid:1) (A.20)and the updated fiber stiffness κ fib c , see Eq. (A.12).27 onservation of energyCircAdapt connects midwall tension T mid c and midwall area A mid c to fiber stress σ fib c and strain E fib c through the law of conservation of energy. With the law of Laplace we get T mid c d A mid c = σ fib c V wall c d E fib c (A.21)and with (A.17) we get for the midwall tension T mid c = σ fib c V wall c A mid c . (A.22)Transmural pressure p trans c is finally computed as follows p trans c = 2 T mid c C mid c . (A.23)Since at the moment external pressures are assumed to be zero, the transmural pressure coincideswith the internal pressure of the contracting chamber p c = p trans c . Appendix A.4. TriSeg model of ventricular interaction
In case that one ODE and one PDE cavity is included in the model, ventricles are modeledas atria above. Otherwise, ventricular and septal midwall volumes are modeled as a ventricularcomposite [58] which is defined by the common radius y mid of the wall junction and the enclosedmidwall cap volumes, see Fig. A.9c. Midwall cap volumes of the right and the left ventricle arecomputed as V midlv = − V lv + V midsv − (cid:0) V walllv + V wallsv (cid:1) ,V midrv = V rv + V midsv + 12 (cid:0) V wallrv + V wallsv (cid:1) . Here, the wall volumes of the left, V walllv , and right, V wallrv , ventricle are constants. The blood poolvolumes of the left, V lv , and right, V rv , ventricle are ODE variables as well as the radius y mid andthe septal midwall volume V midsv . Note that the sign of midwall volume V mid c is positive if wallcurvature is convex to the positive x-direction and negative otherwise.Distance x mid c , see Fig. A.9c, is then computed by the relation V mid c = π x mid c (cid:16) ( x mid c ) + 3( y mid ) (cid:17) , for c ∈ { lv , rv , sv } , hence x mid c = q c − ( y mid ) q c , with q c = (cid:118)(cid:117)(cid:117)(cid:116)(cid:115)(cid:18) π V mid c (cid:19) + ( y mid ) + 3 π V mid c . A mid c = π (cid:16) ( x mid c ) + ( y mid ) (cid:17) , for c ∈ { lv , rv , sv } ,C mid c = 2 x mid c ( x mid c ) + ( y mid ) , for c ∈ { lv , rv , sv } , and used to calculate midwall tension T mid c (A.22).Axial T x c and radial T y c tension components are computed using laws of trigonometry T x c = T mid c sin α, with sin α = 2 x mid c y mid ( x mid c ) + ( y mid ) , for c ∈ { lv , rv , sv } ,T y c = T mid c cos α, with cos α = − ( x mid c ) + ( y mid ) ( x mid c ) + ( y mid ) , for c ∈ { lv , rv , sv } . It is required that that the total midwall tension at junctions is zero, i.e., f ( y mid , V midsv ) := (cid:18) T xlv + T xrv + T xsv T ylv + T yrv + T ysv (cid:19) ! = 0 . (A.24)Eq. (A.24) is solved by an iterative Newton scheme f (cid:48) ( y mid k , V mid k, sv ) (cid:0) ∆y mid k , ∆V mid k, sv (cid:1) (cid:62) = − f ( y mid k , V mid k, sv ) , k = 1 , , . . . (A.25)and increments ∆y mid k and ∆V mid k, sv are added to y mid k and V mid k, sv . The solution of (A.25) in the firststep, i.e., for k = 0 is used to define the ODE updates for the septum ˙ V midsv = 1 τ sv ∆V mid0 , sv , ˙ y mid = 1 τ sv ∆y mid0 , (A.26)where τ sv is a time constant.Consequently, the values for the tensions discussed above are updated and the scheme is iterateduntil convergence. Midwall volumes are updated by V mid c = V c + 12 (cid:0) V wall c + V wallsv (cid:1) , long-axis length l c and cross-sectional area A c and of the cavity are computed by l c = 2 (cid:18) V mid c + 12 (cid:0) V wall c + V wallsv (cid:1)(cid:19) / , (A.27) A c = V mid c + (cid:0) V wall c + V wallsv (cid:1) l c . (A.28)Finally, wave impedance Z c is computed according to Eq. (A.19) and transmural pressure p trans c is29 a) (b)(c) ResistanceValveWall SegmentTube LV WallSV WallRV Wall y -direction y mid V walllv V wallsv V wallrv axis ofrotationalsymmetryintersectionjuction circlepositive x -direction A mid c αV mid c x mid c y mid C mid c T mid c T y c T x c TriSeg Model
LVSVRV V lv p lv V rv p rv V wallrv V wallsv V walllv Syst PulmRightAtrium LeftAtriumSystVeins PulmVeinsPulmArteries SystArteriesAorticValvePulmValveTricuspidValve MitralValve
Figure A.9: TriSeg model of ventricular mechanics. (a) The TriSeg model (gray shading) incorporated in themodular
CircAdapt model of the systemic (Syst) and pulmonary (Pulm) circulations. (b) Cross-section of theventricular composite. (c) Cross-section of a single wall segment through the axis of rotational symmetry. Adaptedwith permission from [58]. computed as total axial force p trans c = 2 T x c y mid , for c ∈ { lv , rv , sv } . Assuming the pressure surrounding the ventricular composite to be zero, internal chamber pressurefor the ventricles is now found as p lv = − p translv ,p rv = p transrv . Appendix A.5. Pericardial mechanics
The four cardiac chambers are supposed to have an additional pressure component due to thepericardium. Pressure p peri exerted by the pericardial sack on atria and ventricles was computedas a non-linear function of pericardial volume V peri , computed as the sum of blood pool and wall30olumes of the four cardiac chambers: V peri = V lv + V rv + V la + V ra + V walllv + V wallrv + V wallla + V wallra (A.29) p peri = p refperi (cid:32) V peri V refperi (cid:33) k peri , (A.30)where p refperi and V refperi are constant reference pressure and volume, respectively, and k peri defines thedegree of non-linearity of the pressure-volume relation.Cavity pressures are updated according to p c = p c + p peri , for c ∈ { lv , rv , la , ra } . (A.31) Appendix A.6. Periphery
Pulmonary (pulm) and systemic (sys) periphery are modeled as resistances. The current pressuredrop ∆p py , for py ∈ { pulm , sys } , is computed as the difference of the pressures in the inflow artery p prox t and the outflow vein p dist t : ∆p py = p prox t − p dist t . Using this, the current flow over the periphery is q py = q ref py (cid:18) r py ∆p py ∆p ref py (cid:19) k py , (A.32)where ∆p ref py is the reference arteriovenous pressure drop; q ref py is the reference flow over the periphery; r py is a scaling factor of the ateriovenous resistances; and k py is a factor that accounts for thenonlinearity of the arteriovenous resistances, see Tabs. A.3 and A.4. Appendix A.7. Connect modules
Volume change of inflow arteries ˙ V prox t and outflow veins ˙ V dist t is now updated by ˙ V dist t += q py ˙ V prox t += q py (A.33)Computation of time derivative of flow across valves and venous-atrial inlet requires in input thecross-sectional area of proximal and distal elements to the channel. ˙ V dist c,t += q v ˙ V prox c,t += q v (A.34) p prox c,t += ˙ V prox c,t Z prox c,t p dist c,t += ˙ V dist c,t Z dist c,t (A.35)31 ppendix A.8. Valve dynamics The pressure drop ( ∆p v ) across a valve is the sum of the effects of inertia due to accelerationin time and the Bernoulli effect, see [75] ∆p v = ρ b l v A v ˙ q v + ρ b (cid:16) ( v out v ) − ( v in v ) (cid:17) , (A.36)where ρ b is the density of blood, A v is the current cross-sectional area of the valve, and l v is thelength of the channel with inertia. If not mentioned otherwise this values is estimated as l v = (cid:112) A open v , with A open v the given cross-sectional area of the open valve, see Tab. A.3. For q v ≥ , v in v is thevelocity proximal to the valve v prox v . v out v is the maximum of the blood velocities in the valve region v max v = max( v dist v , v v , v prox v ) . For q v < which indicates that the valve is leaking v in v is the velocitydistal to the valve v dist v and the outflow velocity is the maximum of the blood velocities in the valveregion v out v = v max v . Using v v = q v /A v we can write ∆p v = p prox v − p dist v = α v ˙ q v + β v q v (A.37)with α v = ρ b l v A v (A.38)the inertia of the channel. The open/closed status of the valve is a function of pressure drop andflow. Valves are clearly open/closed if both pressure drop and flow point in the same direction.With forward pressure drop, the valve opens immediately. With backward pressure and forwardflow, the valve is closing softly by a continuous function A closing v = (cid:114) x v x v + ∆p v (cid:0) A open v − A leak v (cid:1) + A leak v (A.39) x v = 40 ρ b q v | q v | ( A open v ) , where A leak v is the given valve cross-sectional areas of the closed (regurging) valve. Using this thecurrent cross sectional area of the valve is A v = A open v for ∆p v > ,A leak v for ∆p v < and q v < ,A closing v for ∆p v < and q v > . (A.40)We define A min v = min (cid:0) A prox v , A v , A dist v (cid:1) , (A.41)32 a) (b) A prox v A v A dist v A prox v A v A dist v Figure A.10: Schematic of the (a) open and (b) regurging valve, based on [76]. with A prox v and A dist v the cross-sectional area of the proximal and distal cavities or tubes respectively,see (A.1,A.18,A.28) and Fig. A.10. Using this β v is given as β v = ρ b (cid:20)(cid:16) A min v (cid:17) − (cid:16) A prox v (cid:17) (cid:21) for q v ≥ , ρ b (cid:20)(cid:16) A dist v (cid:17) − (cid:16) A min v (cid:17) (cid:21) for q v < . (A.42)Flow over the valve is finally updated using (A.37) by ˙ q v = ∆p v − β v q v α v . (A.43) Appendix A.9. Solve ODE system
A Runge–Kutta–Fehlberg method (RKF45), see Appendix D, is used to solve the system of 26ordinary differential equations (ODEs):8 ODEs: for each of the four tubes and the four cavities we get an ODE to update the volumeusing Eqs. (A.33) and (A.34).2 ODEs: for the septum we update midwall volume and the radius according to (A.26).10 ODEs: for the sarcomeres of each cavity and the septum we update sarcomere contractinglength and contractility using (A.7–A.8).6 ODEs: for each of the four valves and the two outlets we update flow by (A.43).
Appendix B. Finite Element Formulation
Appendix B.1. Variational Formulation
We first ignore the acceleration term in
Eq. (3) and look at the stationary version of the boundaryvalue problem (3–4) and (15). For the full nonlinear elastodynamics problem see Appendix C. Thestationary boundary value problem is formally equivalent to the equations (cid:104)A ( u ) , v (cid:105) Ω − (cid:104)F ( u , p c ) , v (cid:105) Ω = , (B.1) (cid:104) V PDE c ( u ) , q (cid:105) Ω − (cid:104) V ODE c ( p c ) , q (cid:105) Ω = 0 , (B.2)33 arameter Value Unit Description General ρ b . / m blood density t cycle .
585 s cycle time ( = 1 / heartrate ) Tubes: aorta (ao), arteria pulmonalis (ap), venae cavae (vc), and venae pulmonales (vp) A wallt
274 (ao), 141 (ap), 58 (vc), 85 (vp) mm cross-sectional wall area l t
500 (ao), 400 (vc), 200 (ap, vp) mm length of vessel A reft adjacent valve area mm initial cross sectional area k t Chambers: left (lv) and right (rv) ventricle; left (la) and right atrium (ra) V c mL cavity volume h wallc mm constant wall thickness ∆t actc s delays of onset of activation in each beatstarting at t bt ; t act c = t bt + ∆t actc Valves: aortic (av), pulmonary (pv), mitral (mv), and tricuspid (tv) valve;pulmonary (po) and systemic (so) outlet A open v
500 (mv, tv) 400 (av, pv, so, po) mm valve cross-sectional area A leak v A open v (so, po) mm cross-sectional area of closed/regurging valve Periphery: systemic (sys) and pulmonary (pulm) circulation ∆p ref py kPa blood pressure drop in pulmonary/systemiccirculation q ref py (pulm, sys) mL / s reference pulmonary/systemic flow r py CircAdapt model. Adjusted to match patient-specific data. which is valid for all smooth enough vector fields v vanishing on the Dirichlet boundary Γ , D ,testfunctions q that are for the cavity c and otherwise, the duality pairing (cid:104)· , ·(cid:105) Ω , and cavities c ∈ { lv , rv , la , ra } . The second term on the left hand side of the variational equation (B.1) has thephysical interpretation of the rate of internal mechanical work and is given by (cid:104)A ( u ) , v (cid:105) Ω := (cid:90) Ω S ( u ) : Σ ( u , v ) d X , (B.3)with the second Piola–Kirchhoff stress tensor S , see (5), and the directional derivative of theGreen–Lagrange strain tensor Σ ( u , v ) , see [39, 77]. The weak form of the contribution of pressureloads (B.1), right term, is computed using (4) (cid:104)F ( u , p c ) , v (cid:105) Ω = − p c (cid:90) Γ , N J F −(cid:62) ( u ) n out0 · v d s X . (B.4)The first term of the coupling equation (B.2) is computed from (14) using Nanson’s formula and x = X + u by (cid:104) V PDE c ( u ) , q (cid:105) Ω = 13 (cid:90) Γ , N ( X + u ) · J F −(cid:62) n out0 q d s X , (B.5)The second term of (B.2) is computed using the lumped CircAdapt model, see 2.4, for c ∈{ lv , rv , la , ra } . 34 arameter Value Unit Description Tubes p ref t kPa reference tube pressure Sarcomeres L s , ref m reference sarcomere length L elast , iso m length of isometrically stressed serieselastic element v max m / s reference shortening velocity τ R t act , ref c τ D t act , ref c L act0 , ref m contractile element length with zeroactive stress L pas0 , ref m sarcomere length with zero passivestress σ act , max
120 (lv, rv) 84 (la, ra) kPa maximal isometric stress σ pas , max
22 (lv, rv) 50 (la, ra) kPa maximal passive stress d L s , pas m TriSeg Module τ sv Pericardium p refperi V refperi Periphery k py CircAdapt model, fitted to general experimental data in [42]. ppendix B.2. Consistent Linearization To solve the nonlinear variational equations (B.1)–(B.2), with a FE approach we first applya Newton–Raphson scheme, see [78]. Given a nonlinear and continuously differentiable operator F : X → Y a solution to F ( x ) = 0 can be approximated by x k +1 = x k + ∆x,∂F∂x (cid:12)(cid:12)(cid:12)(cid:12) x = x k ∆x = − F ( x k ) , which is looped until convergence. In our case, we have X = (cid:2) H ( Ω , Γ , D ) (cid:3) × R , Y = R , ∆x = ( ∆ u , ∆p c ) (cid:62) , x k = ( u k , p kc ) (cid:62) , and F = ( R u , R p ) (cid:62) . We obtain the following linearized saddle-point problem for each ( u k , p kc ) ∈ (cid:2) H ( Ω , Γ , D ) (cid:3) × R , find ( ∆ u , ∆p c ) ∈ (cid:2) H ( Ω ) (cid:3) × R suchthat (cid:104) ∆ u , A (cid:48) ( u k ) v (cid:105) Ω + (cid:104) ∆ u , F (cid:48) ( u k , p kc ) v (cid:105) Ω + (cid:104) ∆p c , F (cid:48) ( u k , p kc ) v (cid:105) Ω = −(cid:104) R u ( u k , p kc ) , v (cid:105) Ω , (B.6) (cid:104) ∆ u , V PDE c ( u k ) q (cid:105) Ω − (cid:104) ∆p c , V ODE c ( p kc ) q (cid:105) Ω = −(cid:104) R p ( u k , p kc ) , q (cid:105) Ω , (B.7)with the updates u k +1 = u k + ∆ u , (B.8) p k +1 c = p kc + ∆p c , (B.9)and the particular terms are introduced below. The Gâteaux derivative of (B.1) with respect tothe displacement change update ∆ u yields the first (cid:104) ∆ u , A (cid:48) ( u k ) v (cid:105) Ω : = D ∆ u (cid:104)A ( u ) , v (cid:105) Ω | u = u k = (cid:90) Ω S k : Σ ( ∆ u , v ) d X + (cid:90) Ω Σ ( u k , ∆ u ) : C k : Σ ( u k , v ) d X , (B.10)and second term of (B.6) (cid:104) ∆ u , F (cid:48) ( u k , p kc ) v (cid:105) Ω : = D ∆ u (cid:104)F ( u , p c ) , v (cid:105) Ω | u = u k ,p c = p kc = p kc (cid:90) Γ , N J k F −(cid:62) k Grad (cid:62) ∆ u F −(cid:62) k n out0 · v d s X − p kc (cid:90) Γ , N J k ( F −(cid:62) k : Grad ∆ u ) F −(cid:62) k n out0 · v d s X , (B.11)with abbreviations F k := F ( u k ) , J k := det( F k ) , S k := S | u = u k , C k := C | u = u k . ∆p c yields the thirdterm of (B.6) (cid:104) ∆p c , F (cid:48) ( u k , p kc ) v (cid:105) Ω : = D ∆p c (cid:104)F ( u , p c ) , v (cid:105) Ω | u = u k ,p c = p kc = − ∆p c (cid:90) Γ , N J k F −(cid:62) k n out0 · v d s X . (B.12)The residual R u , i.e., the right hand side of (B.6), is computed as (cid:104) R u ( u k , p kc ) , v (cid:105) Ω := (cid:104) A ( u k ) , v (cid:105) Ω − (cid:104)F ( u k , p kc ) , v (cid:105) Ω . (B.13)From (B.5), using the known relations, see, e.g., [77], ∂J∂ F : Grad ∆ u = J F −(cid:62) : Grad ∆ u ∂ F −(cid:62) ∂ F : Grad ∆ u = − F −(cid:62) (Grad ∆ u ) (cid:62) F −(cid:62) we can calculate the first term of (B.7) as the Gâteaux derivative with respect to the update ∆ u (cid:104) ∆ u , V PDE c ( u k ) q (cid:105) Ω : = D ∆ u (cid:104) V PDE c ( u ) , q (cid:105) Ω (cid:12)(cid:12) u = u k = D ∆ u (cid:90) Γ , N (cid:0) X + u k (cid:1) · J k F −(cid:62) k n out0 q d s X = 13 (cid:90) Γ , N J k ( F −(cid:62) k : Grad ∆ u ) x · F −(cid:62) k n out0 q d s X − (cid:90) Γ , N J k x · F −(cid:62) k (Grad ∆ u ) (cid:62) F −(cid:62) k n out0 q d s X + 13 (cid:90) Γ , N J k ∆ u · F −(cid:62) k n out0 q d s X , (B.14)with q a testfunction that is for the surface of cavity c , Γ ,c , and otherwise.The second term of (B.7) is computed as a numerical derivative (cid:104) ∆p c , V ODE c ( p kc ) q (cid:105) Ω : = D ∆p c (cid:104) V ODE c ( p c ) q (cid:105) Ω (cid:12)(cid:12) p c = p kc = 1 (cid:15) (cid:0) V ODE c ( p kc + (cid:15) ) − V ODE c ( p kc ) (cid:1) q, (B.15)where (cid:15) = p kc √ (cid:15) m is chosen according to [79, Chapter 5.7] with (cid:15) m = 2 − ≈ . ∗ − the machineaccuracy.Finally, the residual R p , i.e., the right hand side of (B.7), is computed as (cid:104) R p ( u k , p kc ) , q (cid:105) Ω := (cid:104) V PDE c ( u ) , q (cid:105) Ω − (cid:104) V ODE c ( p c ) , q (cid:105) Ω . (B.16) Appendix B.3. Assembling of the block matrices
To apply the finite element method (FEM) we consider an admissible decomposition of thecomputational domain Ω ⊂ R into M tetrahedral elements τ j and introduce a conformal finite37lement space X h ⊂ H ( Ω ) , N = dim X h of piecewise polynomial continuous basis functions ϕ i . The linearized variational problem (B.6)–(B.7) and a Galerkin FE discretization result in solving the block system to find δu ∈ R N and δp c ∈ R N cav such that K (cid:48) ( u k , p k c ) (cid:18) δuδp c (cid:19) = − K ( u k , p k c ) , K ( u k , p k c ) := − (cid:18) R u ( u k , p k c ) R p ( u k , p k c ) (cid:19) , i.e., (cid:18) ( A (cid:48) − M (cid:48) )( u k , p k c ) B (cid:48) p ( u k ) B (cid:48) u ( u k ) C (cid:48) ( p k c ) (cid:19) (cid:18) δuδp c (cid:19) = − (cid:18) A ( u k ) − B p ( u k , p k c ) V PDE c ( u k ) − V ODE c ( p k c ) (cid:19) , (B.17) u k +1 = u k + δu, (B.18) p k +1 c = p k c + δp c (B.19)with the solution vectors u k ∈ R N and p k c ∈ R N cav at the k -th Newton step. The tangent stiffnessmatrix A (cid:48) ∈ R N × N is calculated from (B.10) according to A (cid:48) ( u k )[ j, i ] := (cid:104) ϕ i , A (cid:48) ( u k ) ϕ j (cid:105) Ω (B.20)and the mass matrix M (cid:48) ∈ R N × N is calculated from (B.11) according to M (cid:48) ( u k , p k c )[ j, i ] := (cid:104) ϕ i , F (cid:48) ( u k , p kc ) ϕ j (cid:105) Ω , (B.21)see also [39, 77].The off-diagonal matrices B (cid:48) u ∈ R N × N cav and B (cid:48) p ∈ R N cav × N in (B.17) are assembled us-ing (B.14) B (cid:48) u ( u k , p k c )[ i, j ] = (cid:104) ϕ j , V PDE c ( u k ) ˆ ϕ i (cid:105) Ω , i = 1 , . . . , N cav (B.22)and using (B.12) B (cid:48) p ( u k , p k c )[ i, j ] = (cid:104) ˆ ϕ j , F (cid:48) ( u k , p kc ) ϕ i (cid:105) Ω , j = 1 , . . . , N cav , (B.23)with the constant shape function ˆ ϕ j = 1 if τ j ∈ Γ ,c and ˆ ϕ j = 0 if τ j / ∈ Γ ,c for c ∈ { lv , rv , la , ra } .Using a technique as described in [60, Sect. 4.2] this assembling procedure can be simplified forclosed cavities such that B (cid:48) p ( u k , p k c ) = (cid:104) B (cid:48) u ( u k , p k c ) (cid:105) (cid:62) . The circulatory compliance matrix C (cid:48) ( p k c ) ∈ R N cav × N cav is computed from (B.15) as C (cid:48) ( p k c )[ i, j ] = (cid:104) ˆ ϕ j , V ODE c ( p kc ) ˆ ϕ i (cid:105) Ω , i, j = 1 , . . . , N cav , (B.24)with the constant shape function ˆ ϕ i , ˆ ϕ j = 1 for cavity c and 0 otherwise, leading to a diagonalmatrix. 38he terms on the upper right hand side A ∈ R N , B p ∈ R N are constructed using (B.13)resulting in R u ( u k , p k c ) = A ( u k ) − B p ( u k , p k c ) with A ( u k )[ i ] := (cid:104)A ( u k ) , ϕ i (cid:105) Ω (B.25)and B p ( u k , p k c )[ i ] := (cid:104)F ( u k , p kc ) , ϕ i (cid:105) Ω . (B.26)Finally, the lower right hand side in (B.17), R p ( u k , p k c ) = V PDE ( u k ) − V ODE ( p k c ) ∈ R N cav , isassembled from (B.16) with V PDE ( u k )[ i ] = (cid:104) V PDE c ( u ) , ˆ ϕ i (cid:105) Ω , i = 1 , . . . , N cav , (B.27)and V ODE ( p k c )[ i ] = (cid:104) V ODE c ( p c ) , ˆ ϕ i (cid:105) Ω , i = 1 , . . . , N cav . (B.28) Appendix C. Generalized- α Scheme
After standard discretization we rewrite Eq. (3) using Eqs. (B.25) and (B.26) as a nonlinearODE reading ρ M α ¨ u ( t ) + R u ( u, t ) = 0 , (C.1)with the mass matrix M α [ i, j ] := (cid:90) Ω ϕ i ( X ) · ϕ j ( X ) d X . Following [80] we reformulate Eq. (C.1) as a first order ODE system by introducing the velocity vρ M α ˙ v ( t ) + R u ( u, t ) = 0 , (C.2) M α ˙ u ( t ) − M α v ( t ) = 0 (C.3)and apply a generalized- α approach [81]. To this end we define three parameters α f := 11 + ρ ∞ , α m := 3 − ρ ∞ ρ ∞ ) , γ := 12 + α m − α f , where the spectral radius ρ ∞ is a parameter between and . With this we introduce ˙ v n + α m := α m ˙ v n +1 + (1 − α m ) ˙ v n , ˙ u n + α m := α m ˙ u n +1 + (1 − α m ) ˙ u n ,v n + α f := α f v n +1 + (1 − α f ) v n ,u n + α f := α f u n +1 + (1 − α f ) u n , and reformulate Eq. (C.2) as ρ M α ˙ v n + α m + R u ( u n + α f ) = 0 , (C.4) M α ˙ u n + α m − M α v n + α f = 0 . (C.5)39ere, the second equation gives us ˙ u n + α m = v n + α f and we get for the velocity update v n +1 = α m α f γ∆t (cid:0) u n +1 − u n (cid:1) + γ − α m γα f ˙ u n + α f − α f v n . (C.6)From this and the relationship by Newmark [82] u n +1 = u n + ∆t (cid:0) γ ˙ u n +1 + (1 − γ ) ˙ u n (cid:1) ,v n +1 = v n + ∆t (cid:0) γ ˙ v n +1 + (1 − γ ) ˙ v n (cid:1) , we obtain ˙ v n +1 = α m α f γ ∆t (cid:0) u n +1 − u n (cid:1) − α f γ∆t v n + γ − γ ˙ v n + γ − α m α f γ ∆t . (C.7)Hence, we can rewrite the whole first order system only dependent on the unknowns u n +1 . Newton’s method for the generalized- α scheme. For the implementation of Newton’s method wecompute ∂ ˙ v n + α m ∂u n +1 = α α f γ ∆t , ∂v n + α f ∂u n +1 = α f α m α f γ∆t = α m γ∆t , ∂u n + α f ∂u n +1 = α f . (C.8)To calculate the solution at the current timestep we assume that we know u n , ˙ u n , v n and ˙ v n fromthe previous time step n and get from Eq. (C.4) for the residual R α ( u kn +1 ) := − ρ M α ˙ v kn + α m − R u ( u kn + α f ) , (C.9)with ˙ v kn + α m := ˙ v n + α m ( u kn +1 ) and u kn + α f := u n + α f ( u kn +1 ) . To increase stability we consider Rayleighdamping by adding the two matrices D mass ( u kn +1 ) = ρ β mass M α v kn + α f , (C.10) D stiff ( u kn +1 ) = β stiff A (cid:48) ( u n ) v kn + α f , (C.11)to the residual (C.9) with v kn + α f := v n + α f ( u kn +1 ) and Rayleigh damping parameters β mass ≥ − , β stiff ≥ .The tangent stiffness matrix is now calculated using (C.8) as A (cid:48) α ( u kn +1 , p k c ,n +1 ) := ρ ∂ ˙ v n + α m ∂u n +1 M α + ∂u n + α f ∂u n +1 (cid:16) A (cid:48) ( u kn + α f ) − M (cid:48) ( u kn + α f , p k c ,n +1 ) (cid:17) = ρ α α f γ ∆t M α + α f (cid:16) A (cid:48) ( u kn + α f ) − M (cid:48) ( u kn + α f , p k c ,n +1 ) (cid:17) , (C.12)with A (cid:48) , and M (cid:48) being the known tangent stiffness matrices from the quasi-stationary elasticitycase, see Eqs. (B.20) and (B.21). When using a coupling with the circulatory system we compute40he off diagonal matrices and lower right hand side, see Eqs. (B.22), (B.23) and (B.27), in terms of u kn + α f . Appendix D. Runge–Kutta–Fehlberg method
This methods is of order O ( h ) but, by performing one extra calculation, the error estimator isof order O ( h ) . This allows for an adaptive stepsize that is determined automatically. y n +1 = y n + h (cid:88) i =1 b i k i ,y ∗ n +1 = y n + h (cid:88) i =1 b i k i ,e n +1 = y n +1 − y ∗ n +1 = h (cid:88) i =1 ( b i − b i ) k i , where k = f ( t n , y n ) ,k = f ( t n + c h, y n + h ( a k )) ,k = f ( t n + c h, y n + h ( a k + a k )) , ... k = f ( t n + c h, y n + h ( a k + a k + · · · + a k )) . and coefficients given in Tab. D.5. i c i a i a i a i a i a i / / / /
32 9 / /
13 1932 / − / / / − / − / / − /
27 2 − / / − / j b j b j b j b j b j b j /
135 0 6656 / / − /
50 2 / Table D.5: Butcher tableau for the Runge–Kutta–Fehlberg method. ppendix E. Direct Schur Complement Solver for a Small Number of Constraints Given the block system A ∈ R n × n , D ∈ R m × m (cid:18) A BC D (cid:19) (cid:18) xy (cid:19) = − (cid:18) fg (cid:19) with B = (cid:0) b · · · b m (cid:1) ∈ R n × m , C = (cid:0) c · · · c m (cid:1) (cid:62) ∈ R m × n , we can write the Schur complement system as ( CA − B − D ) y = g − CA − fx = A − f − A − B y. With r = A − f , S = A − B = (cid:0) s · · · s m (cid:1) ∈ R n × m , s i = A − b i , i = 1 , . . . , m (E.1)we get ( CS − D ) y = g − C rx = r − S y. (E.2)The realization of (E.2) involves m + 1 solves and the inversion of an m × m matrix. Since m isgenerally small this can be done symbolically. [ CS ] ij = c i · s j , for i, j = 1 , . . . , m. eferences [1] L. J. Laslett, P. Alagona, B. A. Clark, J. P. Drozda, F. Saldivar, S. R. Wilson, C. Poe,M. Hart, The worldwide environment of cardiovascular disease: prevalence, diagnosis, therapy,and policy issues: a report from the american college of cardiology, Journal of the AmericanCollege of Cardiology 60 (2012) S1–S49.[2] A. Timmis, N. Townsend, C. P. Gale, A. Torbica, M. Lettino, S. E. Petersen, E. A. Mossialos,A. P. Maggioni, D. Kazakiewicz, H. T. May, et al., European society of cardiology: cardiovas-cular disease statistics 2019, European Heart Journal 41 (2020) 12–85.[3] E. Wilkins, L. Wilson, K. Wickramasinghe, P. Bhatnagar, J. Leal, R. Luengo-Fernandez,R. Burns, M. Rayner, N. Townsend, European cardiovascular disease statistics 2017, EuropeanHeart Network (2017).[4] N. P. Smith, A. de Vecchi, M. McCormick, D. A. Nordsletten, O. Camara, a. F. Frangi,H. Delingette, M. Sermesant, J. Relan, N. Ayache, M. W. Krueger, W. H. W. Schulze, R. Hose,I. Valverde, P. Beerbaum, C. Staicu, M. Siebes, J. Spaan, P. J. Hunter, J. Weese, H. Lehmann,D. Chapelle, R. Rezavi, euHeart: personalized and integrated cardiac care using patient-specificcardiovascular modelling, Interface Focus 1 (2011) 349–364.[5] H. J. Arevalo, F. Vadakkumpadan, E. Guallar, A. Jebb, P. Malamas, K. C. Wu, N. A.Trayanova, Arrhythmia risk stratification of patients after myocardial infarction using per-sonalized heart models, Nature Communications (2016).[6] A. Prakosa, H. J. Arevalo, D. Deng, P. M. Boyle, P. P. Nikolov, H. Ashikaga, J. J. Blauer,E. Ghafoori, C. J. Park, R. C. Blake, F. T. Han, R. S. MacLeod, H. R. Halperin, D. J. Callans,R. Ranjan, J. Chrispin, S. Nazarian, N. A. Trayanova, Personalized virtual-heart technology forguiding the ablation of infarct-related ventricular tachycardia, Nature Biomedical Engineering(2018).[7] M. Strocchi, M. A. Gsell, C. M. Augustin, O. Razeghi, C. H. Roney, A. J. Prassl, E. J.Vigmond, J. M. Behar, J. S. Gould, C. A. Rinaldi, M. J. Bishop, G. Plank, S. A. Niederer,Simulating ventricular systolic motion in a four-chamber heart model with spatially varyingrobin boundary conditions to model the effect of the pericardium, Journal of Biomechanics101 (2020) 109645.[8] G. P. Toorop, G. J. van den Horn, G. Elzinga, N. Westerhof, Matching between feline left ven-tricle and arterial load: Optimal external power or efficiency, American Journal of Physiology- Heart and Circulatory Physiology (1988).[9] D. Nordsletten, M. Mccormick, P. J. Kilner, P. Hunter, D. Kay, N. P. Smith, Fluid-solidcoupling for the investigation of diastolic and systolic human left ventricular function, Inter-national Journal for Numerical Methods in Biomedical Engineering 27 (2011) 1017–1039.[10] E. Karabelas, G. Haase, G. Plank, C. M. Augustin, Versatile stabilized finite element formula-tions for nearly and fully incompressible solid mechanics, Computational Mechanics 65 (2020)193–215.[11] G. Elzinga, N. Westerhof, Pressure and flow generated by the left ventricle against differentimpedances, Circulation Research 32 (1973) 178–186.4312] H. Liu, F. Liang, J. Wong, T. Fujiwara, W. Ye, K.-i. Tsubota, M. Sugawara, Multi-scalemodeling of hemodynamics in the cardiovascular system, Acta Mechanica Sinica 31 (2015)446–464.[13] P. Segers, E. Rietzschel, M. De Buyzere, N. Stergiopulos, N. Westerhof, L. Van Bortel, T. Gille-bert, P. Verdonck, Three-and four-element windkessel models: assessment of their fitting per-formance in a large cohort of healthy middle-aged individuals, Proceedings of the Institutionof Mechanical Engineers, Part H: Journal of Engineering in Medicine 222 (2008) 417–428.[14] N. Stergiopulos, B. E. Westerhof, N. Westerhof, Total arterial inertance as the fourth elementof the windkessel model, American Journal of Physiology-Heart and Circulatory Physiology276 (1999) H81–H88.[15] J.-J. Wang, A. B. O’Brien, N. G. Shrive, K. H. Parker, J. V. Tyberg, Time-domain represen-tation of ventricular-arterial coupling as a windkessel and wave system, American Journal ofPhysiology-Heart and Circulatory Physiology 284 (2003) H1358–H1368.[16] N. Westerhof, G. Elzinga, Normalized input impedance and arterial decay time over heartperiod are independent of animal size, American Journal of Physiology-Regulatory, Integrativeand Comparative Physiology 261 (1991) R126–R133.[17] N. Westerhof, N. Stergiopulos, Models of the arterial tree., Studies in health technology andinformatics 71 (2000) 65–77.[18] J. Alastruey, K. H. Parker, S. J. Sherwin, Arterial pulse wave haemodynamics, BHR Group -11th International Conferences on Pressure Surges (2012) 401–442.[19] P. J. Blanco, S. M. Watanabe, E. A. Dari, M. A. R. Passos, R. A. Feijóo, Blood flow distributionin an anatomically detailed arterial network model: criteria and algorithms, Biomechanics andModeling in Mechanobiology 13 (2014) 1303–1330.[20] L. Formaggia, D. Lamponi, A. Quarteroni, One-dimensional models for blood flow in arteries,Journal of Engineering Mathematics 47 (2003) 251–276.[21] J. P. Mynard, P. Nithiarasu, A 1D arterial blood flow model incorporating ventricular pressure,aortic vaive ana regional coronary flow using the locally conservative Galerkin (LCG) method,Communications in Numerical Methods in Engineering (2008).[22] L. O. Müller, E. F. Toro, A global multiscale mathematical model for the human circula-tion with emphasis on the venous system, International Journal for Numerical Methods inBiomedical Engineering 30 (2014) 681–725.[23] L. Marx, M. A. F. Gsell, A. Rund, F. Caforio, A. J. Prassl, G. Toth-Gayor, T. Kuehne, C. M.Augustin, G. Plank, Personalization of electro-mechanical models of the pressure-overloadedleft ventricle: fitting of Windkessel-type afterload models, Philosophical Transactions of theRoyal Society A: Mathematical, Physical and Engineering Sciences 378 (2020) 20190342.[24] S. A. Niederer, G. Plank, P. Chinchapatnam, M. Ginks, P. Lamata, K. S. Rhode, C. A. Rinaldi,R. Razavi, N. P. Smith, Length-dependent tension in the failing heart and the efficacy of cardiacresynchronization therapy, Cardiovascular Research 89 (2011) 336.4425] T. Arts, T. Delhaas, P. Bovendeerd, X. Verbeek, F. Prinzen, Adaptation to mechanical loaddetermines shape and properties of heart and circulation: the CircAdapt model, Am. J. Physiol.Heart Circ. Physiol. 288 (2005) 1943–1954.[26] P. J. Blanco, R. A. Feijóo, et al., A 3D-1D-0D Computational model for the entire cardio-vascular system, Computational Mechanics, eds. E. Dvorking, M. Goldschmit, M. Storti 29(2010) 5887–5911.[27] G. Guidoboni, L. Sala, M. Enayati, R. Sacco, M. Szopos, J. M. Keller, M. Popescu, L. Despins,V. H. Huxley, M. Skubic, Cardiovascular function and ballistocardiogram: a relationshipinterpreted via mathematical modeling, IEEE Transactions on Biomedical Engineering 66(2019) 2906–2917.[28] M. L. Neal, J. B. Bassingthwaighte, Subject-specific model estimation of cardiac output andblood volume during hemorrhage, Cardiovascular Engineering 7 (2007) 97–120.[29] S. Paeme, K. T. Moorhead, J. G. Chase, B. Lambermont, P. Kolh, V. D’orio, L. Pierard,M. Moonen, P. Lancellotti, P. C. Dauby, et al., Mathematical multi-scale model of the cardio-vascular system including mitral valve dynamics. application to ischemic mitral insufficiency,Biomedical engineering online 10 (2011) 86.[30] T. S. E. Eriksson, A. J. Prassl, G. Plank, G. A. Holzapfel, Influence of myocardial fiber/sheetorientations on left ventricular mechanical contraction, Math Mech Solids 18 (2013) 592–606.[31] R. C. Kerckhoffs, M. L. Neal, Q. Gu, J. B. Bassingthwaighte, J. H. Omens, A. D. McCulloch,Coupling of a 3D finite element model of cardiac ventricular mechanics to lumped systemsmodels of the systemic and pulmonic circulation, Annals of biomedical engineering 35 (2007)1–18.[32] T. P. Usyk, I. J. LeGrice, A. D. McCulloch, Computational model of three-dimensional cardiacelectromechanics, Computing and Visualization in Science 4 (2002) 249–257.[33] T. Fritz, C. Wieners, G. Seemann, H. Steen, O. Dössel, Simulation of the contraction ofthe ventricles in a human heart model including atria and pericardium, Biomechanics andModeling in Mechanobiology 13 (2014) 627–641.[34] V. Gurev, T. Lee, J. Constantino, H. J. Arevalo, N. A. Trayanova, Models of cardiac electrome-chanics based on individual hearts imaging data, Biomechanics and Modeling in Mechanobi-ology 10 (2011) 295–306.[35] V. Gurev, P. Pathmanathan, J.-L. Fattebert, H.-F. Wen, J. Magerlein, R. a. Gray, D. F.Richards, J. J. Rice, A high-resolution computational model of the deforming human heart,Biomechanics and Modeling in Mechanobiology 14 (2015) 829–849.[36] M. Hirschvogel, M. Bassilious, L. Jagschies, S. M. Wildhirt, M. W. Gee, A monolithic 3D-0Dcoupled closed-loop model of the heart and the vascular system: Experiment-based parameterestimation for patient-specific cardiac mechanics, International Journal for Numerical Methodsin Biomedical Engineering 33 (2017).[37] J. Sainte-Marie, D. Chapelle, R. Cimrman, M. Sorine, Modeling and estimation of the cardiacelectromechanical activity, Computers and Structures 84 (2006) 1743–1759.4538] K. L. Sack, E. Aliotta, D. B. Ennis, J. S. Choy, G. S. Kassab, J. M. Guccione, T. Franz,Construction and validation of subject-specific biventricular finite-element models of healthyand failing swine hearts from high-resolution DT-MRI, Frontiers in Physiology 9 (2018) 1–19.[39] C. M. Augustin, A. Neic, M. Liebmann, A. J. Prassl, S. A. Niederer, G. Haase, G. Plank,Anatomically accurate high resolution modeling of human whole heart electromechanics: Astrongly scalable algebraic multigrid solver method for nonlinear deformation, Journal ofComputational Physics 305 (2016) 622–646.[40] C. M. Augustin, A. Crozier, A. Neic, A. J. Prassl, E. Karabelas, T. Ferreira da Silva, J. F. Fer-nandes, F. Campos, T. Kuehne, G. Plank, T. da Silva, J. F. Fernandes, F. Campos, T. Kuehne,G. Plank, Patient-specific modeling of left ventricular electromechanics as a driver for haemo-dynamic analysis, Europace 18 (2016) iv121–iv129.[41] A. Crozier, C. M. Augustin, A. Neic, A. J. Prassl, M. Holler, T. E. Fastl, A. Hennemuth,K. Bredies, T. Kuehne, M. J. Bishop, S. A. Niederer, G. Plank, Image-Based Personalization ofCardiac Anatomy for Coupled Electromechanical Modeling, Annals of Biomedical Engineering44 (2016) 58–70.[42] J. Walmsley, T. Arts, N. Derval, P. Bordachar, H. Cochet, S. Ploux, F. W. Prinzen, T. Delhaas,J. Lumens, Fast Simulation of Mechanical Heterogeneity in the Electrically AsynchronousHeart Using the MultiPatch Module, PLOS Computational Biology 11 (2015) e1004284.[43] F. W. Prinzen, R. W. Mills, R. N. Cornelussen, L. J. Mulligan, M. Strik, L. M. Rademakers,N. D. Skadsberg, A. Van Hunnik, M. Kuiper, A. Lampert, T. Delhaas, Left ventricular septaland left ventricular apical pacing chronically maintain Cardiac contractile coordination, pumpfunction and efficiency, Circulation: Arrhythmia and Electrophysiology 2 (2009) 571–579.[44] X. A. Verbeek, K. Vernooy, M. Peschar, R. N. Cornelussen, F. W. Prinzen, Intra-ventricularresynchronization for optimal left ventricular function during pacing in experimental left bundlebranch block, Journal of the American College of Cardiology 42 (2003) 558–567.[45] CIBC, 2016. URL: , seg3D: Volumetric Image Segmentation and Vi-sualization. Scientific Computing and Imaging.[46] A. Neic, M. A. F. Gsell, E. Karabelas, A. J. Prassl, G. Plank, Automating image-based meshgeneration and manipulation tasks in cardiac modeling workflows using Meshtool, SoftwareX11 (2020) 100454.[47] J. D. Bayer, R. C. Blake, G. Plank, N. A. Trayanova, A novel rule-based algorithm for assigningmyocardial fiber orientation to computational heart models, Annals of Biomedical Engineering40 (2012) 2243–2254.[48] D. D. Streeter, H. M. Spotnitz, D. P. Patel, J. Ross, E. H. Sonnenblick, Fiber orientation inthe canine left ventricle during diastole and systole., Circulation research (1969).[49] J. Bayer, A. J. Prassl, A. Pashaei, J. F. Gomez, A. Frontera, A. Neic, G. Plank, E. J. Vigmond,Universal ventricular coordinates: A generic framework for describing position within the heartand transferring data, Medical Image Analysis 45 (2018) 83–93.4650] P. J. Flory, Thermodynamic relations for high elastic materials, Trans Faraday Soc 57 (1961)829–838.[51] S. Land, S. A. Niederer, Influence of atrial contraction dynamics on cardiac function, Inter-national journal for numerical methods in biomedical engineering 34 (2018) e2931.[52] T. P. Usyk, R. Mazhari, A. D. McCulloch, Effect of laminar orthotropic myofiber architectureon regional stress and strain in the canine left ventricle, J. Elast. 61 (2000) 143–164.[53] M. Genet, L. C. Lee, R. Nguyen, H. Haraldsson, G. Acevedo-Bolton, Z. Zhang, L. Ge, K. Or-dovas, S. Kozerke, J. M. Guccione, Distribution of normal human left ventricular myofiberstress at end diastole and end systole: a target for in silico design of heart failure treatments,Journal of Applied Physiology 117 (2014) 142–152.[54] J. C. Walker, M. B. Ratcliffe, P. Zhang, A. W. Wallace, B. Fata, E. W. Hsu, D. Saloner, J. M.Guccione, MRI-based finite-element analysis of left ventricular aneurysm, American Journalof Physiology-Heart and Circulatory Physiology 289 (2005) H692–H700.[55] A. Neic, F. O. Campos, A. J. Prassl, S. A. Niederer, M. J. Bishop, E. J. Vigmond, G. Plank,Efficient computation of electrograms and ECGs in human whole heart simulations using areaction-eikonal model, Journal of Computational Physics 346 (2017) 191–211.[56] C. M. Costa, E. Hoetzl, B. M. Rocha, A. J. Prassl, G. Plank, Automatic ParameterizationStrategy for Cardiac Electrophysiology Simulations., Computing in cardiology 40 (2013) 373–376.[57] K. H. W. J. ten Tusscher, D. Noble, P. J. Noble, A. V. Panfilov, A model for human ventriculartissue., American Journal of Physiology. Heart and Circulatory Physiology 286 (2004) H1573–H1589.[58] J. Lumens, T. Delhaas, B. Kirn, T. Arts, Three-Wall Segment (TriSeg) Model DescribingMechanics and Hemodynamics of Ventricular Interaction, Annals of Biomedical Engineering37 (2009) 2234–2255.[59] V. Gurev, T. Lee, J. Constantino, H. Arevalo, N. A. Trayanova, Models of cardiac electrome-chanics based on individual hearts imaging data, Biomechanics and modeling in mechanobiol-ogy 10 (2011) 295–306.[60] T. Rumpel, K. Schweizerhof, Volume-dependent pressure loading and its influence on thestability of structures, International Journal for Numerical Methods in Engineering 56 (2003)211–238.[61] S. Sugiura, T. Washio, A. Hatano, J. Okada, H. Watanabe, T. Hisada, Multi-scale simulationsof cardiac electrophysiology and mechanics using the University of Tokyo heart simulator,Progress in Biophysics and Molecular Biology 110 (2012) 380–9.[62] S. Klotz, M. L. Dickstein, D. Burkhoff, A computational method of prediction of the end-diastolic pressure–volume relationship by single beat, Nature Protocols 2 (2007) 2152–2158.4763] E. Willemen, R. Schreurs, P. R. Huntjens, M. Strik, G. Plank, E. Vigmond, J. Walmsley,K. Vernooy, T. Delhaas, F. W. Prinzen, J. Lumens, The Left and Right Ventricles RespondDifferently to Variation of Pacing Delays in Cardiac Resynchronization Therapy: A CombinedExperimental- Computational Approach, Front. Physiol. 10 (2019) 1–13.[64] D. Burkhoff, I. Mirsky, H. Suga, Assessment of systolic and diastolic ventricular propertiesvia pressure-volume analysis: a guide for clinical, translational, and basic researchers., Am. J.Physiol. Heart Circ. Physiol. 289 (2005) H501–12.[65] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener,V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, a. Dave A. M, L. C. McInnes, R. T.Mills, T. Munson, K. Rupp, P. Sanan, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, PETScUsers Manual, Technical Report ANL-95/11 - Revision 3.10, Argonne National Laboratory,2018. URL: .[66] V. E. Henson, U. M. Yang, BoomerAMG: A parallel algebraic multigrid solver and precondi-tioner, in: Applied Numerical Mathematics, 2002, pp. 155–177.[67] E. Vigmond, R. Weber dos Santos, A. Prassl, M. Deo, G. Plank, Solvers for the cardiacbidomain equations, Prog Biophys Mol Biol 96 (2008) 3–18.[68] E. D. Carruth, A. D. McCulloch, J. H. Omens, Transmural gradients of myocardial structureand mechanics: Implications for fiber stress and strain in pressure overload, Progress inBiophysics and Molecular Biology 122 (2016) 215–226.[69] Y.-c. Fung, Biomechanics: circulation, Springer Science & Business Media, 2013.[70] J. H. Omens, Y.-C. Fung, Residual strain in rat left ventricle, Circulation Research 66 (1990)37–45.[71] E. K. Rodriguez, J. H. Omens, L. K. Waldman, A. D. McCulloch, Effect of residual stress ontransmural sarcomere length distributions in rat left ventricle, American Journal of Physiology- Heart and Circulatory Physiology (1993).[72] E. Karabelas, M. A. F. Gsell, C. M. Augustin, L. Marx, A. Neic, A. J. Prassl, L. Goubergrits,T. Kuehne, G. Plank, Towards a Computational Framework for Modeling the Impact of AorticCoarctations Upon Left Ventricular Load, Frontiers in Physiology 9 (2018) 1–20.[73] V. Gurev, J. Constantino, J. J. Rice, N. a. Trayanova, Distribution of electromechanical delayin the heart: insights from a three-dimensional electromechanical model., Biophys. J. 99 (2010)745–54.[74] T. J. Arts, P. H. M. Bovendeerd, F. W. Prinzen, R. S. Reneman, Relation between leftventricular cavity pressure and volume and systolic fibre stress and strain in the wall, Biophys.J. 59 (1991) 93–102.[75] M. S. Firstenberg, N. Greenberg, N. Smedira, P. McCarthy, M. J. Garcia, J. D. Thomas,Noninvasive assessment of mitral inertness: clinical results with numerical model validation,Comput. Cardiol. (2001) 613–616. 4876] C. M. Otto, Valvular aortic stenosis: disease severity and timing of intervention., Journal ofthe American College of Cardiology 47 (2006) 2141–51.[77] G. A. Holzapfel, Nonlinear Solid Mechanics. A Continuum Approach for Engineering, JohnWiley & Sons Ltd, Chichester, 2000.[78] P. Deuflhard, Newton Methods for Nonlinear Problems: Affine Invariance and Adaptive Algo-rithms, volume 35, Springer Science & Business Media, 2011.[79] W. H. Press, S. A. Teukolsky, W. T. Vettering, B. P. Flannery, Numerical Recipies: The Artof Scientific Computing, 3rd ed., Cambridge University Press, 2007.[80] C. Kadapa, W. G. Dettmer, D. Perić, On the advantages of using the first-order generalised-alpha scheme for structural dynamic problems, Computers and Structures 193 (2017) 226–238.[81] J. Chung, G. M. Hulbert, A Time Integration Algorithm for Structural Dynamics With Im-proved Numerical Dissipation: The Generalized- αα