A computationally efficient physiologically comprehensive 3D–0D closed-loop model of the heart and circulation

Computer models of cardiac electro-mechanics (EM) show promise as an effective means for the quantitative analysis of clinical data and, potentially, for predicting therapeutic responses. To realize such advanced applications methodological key challenges must be addressed. Enhanced computational efficiency and robustness is crucial to facilitate, within tractable time frames, model personalization, the simulation of prolonged observation periods under a broad range of conditions, and physiological completeness encompassing therapy-relevant mechanisms is needed to endow models with predictive capabilities beyond the mere replication of observations. Here, we introduce a universal feature-complete cardiac EM modeling framework that builds on a flexible method for coupling a 3D model of bi-ventricular EM to the physiologically comprehensive 0D CircAdapt model representing atrial mechanics and closed-loop circulation. A detailed mathematical description is given and efficiency, robustness, and accuracy of numerical scheme and solver implementation are evaluated. After parameterization and stabilization of the coupled 3D–0D model to a limit cycle under baseline conditions, the model’s ability to replicate physiological behaviors is demonstrated, by simulating the transient response to alterations in loading conditions and contractility, as induced by experimental protocols used for assessing systolic and diastolic ventricular properties. Mechanistic completeness and computational efficiency of this novel model render advanced applications geared towards predicting acute outcomes of EM therapies feasible.


Introduction
Cardiovascular diseases (CVDs) are the primary cause of mortality and morbidity in industrialized nations, posing a significant burden on health care systems worldwide [1,2,3]. Despite continuous 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 -the heart is an electrically controlled mechanical pump driving blood through the circulatory system. Advanced clinical modalities provide a wealth of disparate data, but effective tools allowing their comprehensive quantitative analysis are lacking. Computer models able to capture mechanistic relations between clinical observations quantitatively show promise to fill this void. In recent single physics cardiac electrophysiology (EP) studies, the added value of models in improving therapy stratification [5] and planning [6,7] has been demonstrated already.
Multiphysics models of cardiovascular EM are even more challenging to apply in a clinical context. Their utility depends on the ability to comprehensively represent mechanisms underlying a broader range of physiological function, and to tailor these to approximate -with acceptable fidelity -anatomy and cardiovascular function of a given patient. Such models are complex as all major mechanisms governing a heart beat bidirectionally interact with each other and, thus, must be taken into account. These comprise models of cardiac EP producing electrical activation and repolarization patterns that drive EM coupling to models of contractile function, cardiac mechanics describing deformation and stresses under given mechanical boundary and hemodynamic loading conditions imposed by the intra-thoracic embedding of the heart and the circulatory system, respectively.
Pumping function is regulated through a bidirectional interaction between the heart and both the systemic and pulmonary vascular systems. The circulatory system as an extracardiac factor imposes a pressure and volume load upon the heart and, vice versa, pressure and flow in the circulatory system are determined by the mechanical state of cardiac cavities. Optimal function depends on matching the coupling between these 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 are relevant for investigating flow patterns or wall shear stresses, but are less suitable for systems level investigations. Simpler, computationally less costly 0D and 1D lumped models have been preferred to provide appropriate hemodynamic loading conditions to the heart [11,12].
Most EM modeling studies consider ventricular afterload only represented by lumped 0D Windkessel type models comprising 2-, 3-or 4-elements [13,14,15,16,17,18,19], or, less common, by 1D models derived from Navier-Stokes equations [20,21,22,23,24]. The latter also account for pulse wave transmission and reflection, but identifying parameters is more challenging than for 0D models [25]. A fundamental limitation of isolated models of pre-and afterload is the lack of regulatory loops which respond to altered loading or contractility in one chamber by balancing preload conditions in all chambers until a new stable limit cycle with common compatible stroke volumes is reached. Isolated afterload models are thus best suited for approximating the immediate responses in a single beat [26], but less so for predicting transient behaviors over multiple beats. Closed-loop circulatory systems [27,28,29,30,31] take into account these feedback mechanisms and 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 0D closed-loop model of the circulatory system remains challenging. Hydrostatic pressure, p, in cavities and blood flow, q, between cavities and circulatory system serve as coupling variables that act as pressure boundary condition and impose volume constraints on the 3D cavity models. Previous studies addressed 3D solid-0D fluid coupling problems using simpler partitioned [32,33,34] or more advanced strongly coupled monolithic approaches [35,36,37,38,25,39]. Yet, reports on coupling of a closed-loop 0D to 3D solid models are sparse. Mostly simplified circuit models [37,33,38,40], were used for simulating a single heart beat, where fixed compliances and 0D chambers based on time-varying elastance models are used are used that do not account for pressure-volume relations or the Frank-Starling effect, respectively. Thus, attempts to demonstrate agreement with known physiological principles -fundamental to cardiac pump function -under experimental protocols requiring multibeat simulations have been limited.
Based on previous work on cardiac EM models [41,42,43] we report on the development of a monolithic 3D solid -0D fluid coupling approach. Feasibility is demonstrated by building a 3D canine bi-ventricular EM PDE model coupled to the state-of-the art CircAdapt model [27,44] -a non-linear 0D closed-loop ODE model of the cardiovascular system that implements dynamic adaptation processes based on physiological principles -to represent physiologically realistic atrial EM as well as systemic and pulmonary circulation. A detailed description of numerical underpinnings is given, including a complete mathematical description of the CircAdapt model in a single manuscript that has been lacking so far. Efficiency, robustness, and accuracy of numerical scheme and solver implementation are evaluated. The coupled model is first parameterized and stabilized to a limit cycle representing baseline conditions, and then rigorously tested by demonstrating its ability to predict physiological behaviors under experimental standard protocols altering loading conditions and contractility that are used for the experimental assessment of systolic and diastolic ventricular properties. Transient responses under these protocols are simulated over prolonged observation periods, covering up to 25 beats. The presented framework can be considered a first feature-complete realization of an universal cardiac EM simulator that can be applied, given appropriate parameterization and initialization, under a much broader range of protocols and conditions as any previously reported model.

Experimental data acquisition
In a previous study, see [45], mongrel dog data were acquired to investigate the influence of different pacing protocols on cardiac mechanics, pump function and efficiency. 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 and Other Scientific Purposes (86/609/EU). The protocol was approved by the Maastricht University Experimental Animal Committee. Anatomical Magnetic Resonance Images (MRI) were acquired on a Philips Gyroscan 1.5 T (NT, Philips Medical Systems, Best, the Netherlands) using a standard synergy receiver coil for thorax examinations. Images of seven short-axis cross-sections of slice thickness 8 mm with 0 mm inter-slice distance were obtained to capture the whole heart. LV pressure and volume were determined using the conductance catheter technique (CD-Leycom, The Netherlands), see [46], and the signals were digitized at 1 kHz.

Biventricular Finite Element Models
Multilabel segmentations of right ventricular (RV) (tag 36) and LV blood pool (tag 31) and of the LV myocardium (tag 1), see Fig. 1, were generated from seven MRI short axis slices using Seg3D [47]. Each slice was first segmented semi-automatically using thresholding techniques with manual correction. Segmentations were upsampled to isotropic resolution, followed by an automated iterative erosion and dilation smoothing scheme implemented in Meshtool [48]. The RV wall (tag 6) and lids representing the atrio-ventricular valves were automatically generated by dilation of the adjacent blood pool (tags 41 and 46). Biventricular multilabel meshes were created then from labeled segmentations [43] using the Computational Geometry Algorithms Library, CGAL (www.cgal.org) and subsequently smoothed with Meshtool [48]. A rule-based method according to [49] was applied to define fiber and sheet architecture, with fiber angles changing linearly from −60 • at the epicardium to +60 • at the endocardium [50]. Universal ventricular coordinates were computed [51] to support the flexible definition of stimulation sites and mechanical boundary conditions. Two meshes of different resolution were generated, a coarse mesh to reduce computational expenses and to facilitate the fast exploration of experimental protocols over prolonged observation periods, and a higher resolution mesh for investigating potential inaccuracies introduced by the coarser spatial resolution. For the coarse mesh, average edge lengths of ∼ 3.4 mm and ∼ 2.4 mm, were chosen in LV and RV, respectively, to ascertain that at least two elements were generated transmurally across the myocardial walls, as illustrated in Fig. 1. For the finer mesh, average edge lengths of ∼ 1.3 mm and ∼ 1.2 mm, were chosen in LV and RV, respectively.

Electromechanical PDE Model
Tissue Mechanics. Cardiac tissue is mechanically characterized as a hyperelastic, nearly incompressible, orthotropic material with a nonlinear stress-strain relationship. The deformation gradient F describes the deformation u of a body from the reference configuration Ω 0 (X) to the current configuration Ω t (x), By convention, we denote J = det F > 0 and introduce the right Cauchy-Green tensor C = F F. The nearly incompressible behavior is modeled by a multiplicative decomposition of the deformation gradient [52] of the form Mechanical deformation is described by Cauchy's first equation of motion given as with initial conditions u(X, 0) = 0,u(X, 0) = 0.
Here, ρ 0 is the density in reference configuration;ü are nodal accelerations;u are nodal velocities; S(u, X) is the second Piola-Kirchhoff stress tensor; and Div denotes the divergence operator in the reference configuration. The boundary of the bi-ventricular models was decomposed in three parts, ∂Ω 0 = Γ endo,0 ∪ Γ epi,0 ∪ Γ base,0 , with Γ endo,0 the endocardium, Γ epi,0 the epicardium, and Γ base,0 the base of the ventricles.
Normal stress boundary conditions were imposed on the endocardium with p(t) the pressure and n out 0 the outer normal vector; omni-directional spring type boundary conditions constrained the ventricles at the basal cut plane Γ base,0 [53]; and to simulate the mechanical constrains imposed by the pericardium spatially varying normal Robin boundary conditions were applied at the epicardium Γ epi,0 [7].
Apart from external loads the deformation of cardiac tissue is in particular governed by active stresses intrinsically generated during contraction. To simulate both the active and passive properties of the tissue, the total stress S is additively decomposed according to 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 where Ψ is a strain-energy function to model the orthotropic behavior of cardiac tissue. The prevailing orientation of myocytes, referred to as fiber orientation, is denoted as f 0 . Individual myocytes are surrounded and interconnected by collagen, forming sheets, which is described by the sheet orientation s 0 , perpendicular to f 0 . Together with the sheet-normal axis n 0 , orthogonal to the sheet and the fiber orientations, this forms a right-handed orthonormal set of basis vectors. Following Usyk et al. [54] the orthotropic constitutive relation is defined as where the first term is the volumetric energy with the bulk modulus κ 0 kPa which penalizes local volume changes to enforce near incompressible behavior of the tissue, parameter a is a stress-like scaling parameter, and the term in the exponent is Here, b • are dimensionless parameters and the directional strains read with E = 1 2 (C − I) the modified isochoric Green-Lagrange strain tensor. All passive material parameters are given in Tab. 2.
Active Stress. Stresses due to active contraction are assumed to be orthotropic with full contractile force along the myocyte fiber orientation f 0 and 40 % contractile force along the sheet orientation s 0 [55,56]. Thus, the active stress tensor is defined as where S a is the scalar active stress describing the contractile force. A simplified phenomenological contractile model was used to represent active stress generation [26]. Owing to its small number of parameters and its direct relation to clinically measurable quantities such as peak pressure, and the maximum rate of rise of pressure this model is fairly easy to fit and thus very suitable for being used in clinical EM modeling studies. Briefly, the active stress transient is given by and t s is the onset of contraction; φ(λ) is a non-linear length-dependent function in which λ is the fiber stretch and λ 0 is the lower limit of fiber stretch below which no further active tension is generated; t a is the local activation time from Eq. (12), defined when the local transmembrane potential passes the threshold voltage V m,thresh ; t emd is the EM delay between the onsets of electrical depolarization and active stress generation; S peak is the peak isometric tension; t dur is the duration of active stress transient; τ c is time constant of contraction; τ c0 is the baseline time constant of contraction; ld up is the length-dependence of τ c ; τ r is the time constant of relaxation; and ld is the degree of length dependence. For the parameter values used in the simulations see Tab. 2. Note that active stresses in this simplified model are only length-dependent, but dependence on fiber velocity,λ, is ignored.

Electrophysiology.
A recently developed reaction-eikonal (R-E) model [57] was employed to generate electrical activation sequences which serve as a trigger for active stress generation in cardiac tissue. The hybrid R-E model combines a standard reaction-diffusion (R-D) model based on the monodomain equation with an eikonal model. Briefly, the eikonal equation is given as where (∇ X ) is the gradient with respect to the end-diastolic reference configuration Ω 0 ; t a is a positive function describing the wavefront arrival time at location X ∈ Ω 0 ; and t 0 are initial activations at locations Γ * 0 ⊆ Γ N,0 . The symmetric positive definite 3 × 3 tensor V(X) holds the squared velocities (v f (X), v s (X), v n (X)) associated to the tissue's eigenaxes f 0 , s 0 , and n 0 . The arrival time function t a (X) was subsequently used in a modified monodomain R-D model given as with β the membrane surface-to-volume ratio; C m the membrane capacitance; V m the unknown transmembrane voltage; σ m the monodomain conductivity tensor which holds the scalar conductivities (g f (X), g s (X), g n (X)) and is coupled to V(X) proportionally [58]; and I ion the membrane ionic current density. Additionally, an arrival time dependent foot current, I foot (t a ), was added which is designed to mimic subthreshold electrotonic currents to produce a physiological foot of the action potential. The key advantage of the R-E model is its ability to compute activation sequences at coarser spatial resolutions that are not afflicted by the spatial undersampling artifacts leading to conduction slowing or even numerical conduction block, as it is observed in standard R-D models [59]. Ventricular EP was represented by the tenTusscher-Noble-Noble-Panfilov model of the human ventricular myocyte [60].
Computation of Volumes. To compute the flow across the interface between 3D cavities and the 0D cardiovascular system, the cavitary volume of each chamber that is described as a 3D PDE model has to be tracked as a function of time: drives a positive flow into the circulatory system. In a pure EM simulation context where the fluid domain is not modeled explicitly, the cavitary blood pool volume is not discretized, only the surface Γ enclosing the volume is known. Assuming that the entire surface of the cavitary volume is available, that is, also the faces representing the valves are explicitly discretized, the enclosed volume V PDE can be computed from this surface using the divergence theorem Using this approach, the volume V PDE (x, t) can be computed for each state of deformation at time t and the flow can be derived by a numerical approximation using a difference quotient.
2.4. Lumped ODE model of the circulatory system: the CircAdapt model CircAdapt [27], as shown schematically in Fig. 2, is a lumped 0D model of heart and circulation. It enables real-time simulation of cardiovascular system dynamics under a wide variety of physiological and pathophysiological situations. The entire cardiovascular system is modeled as a concatenation of modules: a tube module representing the systemic and pulmonary arteries and veins (Appendix A.1); a chamber module modeling actively contracting chambers, i.e., left and right atria and ventricles (Appendix A.3), respectively, where myofiber mechanics and contraction is described by a sarcomere module (Appendix A.2); following Lumens et al. [61] this also includes inter-ventricular mechanical interaction through the inter-ventricular septum (Appendix A.4); a valve module representing the aortic, mitral, pulmonary, and tricuspid valves (Appendix A.8); a module representing systemic and pulmonary peripheral microvasculatures (Appendix A.6); and a module accounting for effects of the pericardium (Appendix A.5). The modules are connected by  flows over valves and venous-atrial inlets (Appendix A.7). The whole lumped model consists of 26 ordinary differential equations (ODEs) which are solved using an adaptive Runge-Kutta-Fehlberg method (RKF45) (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 wall areas, which in turn determine strain of the myofibers in the wall. Strain is used to calculate myofiber 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 pressures are found by adding the transmural pressures to the intra-pericardial pressure surrounding the 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 can be straightforwardly replaced by the more complicated finite element (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 [44] and can also be downloaded from the CircAdapt website (http://www.circadapt.org).

PDE-ODE Coupling
We introduce the set of cavities C = {LV, RV, LA, RA}, with the left ventricle (LV), the left atrium (LA), the right ventricle (RV), and the right atrium (RA); the set of cavities C PDE ⊆ C that are modeled as a 3D PDE model; and the set of cavities C ODE = C \ C PDE that are modeled as a 0D ODE model. Coupling between PDE and ODE models can be achieved in various ways. Fundamentally, the problem is to find the new state of deformation u n+1 as a function of the pressure p n+1 in a given cavity at time n+1. The pressure p n+1 is applied as a Neumann boundary condition at the cavitary surface, see Eq. (4). This pressure is not known and has to be determined in a way which depends on the current state of the cavity. Basically, two scenarios have to be considered: (i) when all valves are closed, the cavity is in an isovolumetric state. That is, the muscle enclosing the cavity may deform, but the volume has to remain constant. Therefore if active stresses vary over time during an isovolumetric phase, the pressure p n+1 in the cavity has to vary as well to keep the cavitary 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 or of a connected cavity. Thus p n+1 has to 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 cardiac deformation and as such the two models are tightly bidirectionally coupled.
The simplest approach for the PDE-ODE coupling is to determine p n+1 using a partitioned scheme [32,33,34]. During ejection phases this is achieved by updating cavity volumes and flow based on the current prediction on the change in the state of deformation under the currently predicted pressure p n+1 . In this scenario the pressure boundary condition in each non-linear solver step is modified within each Newton iteration k. The new prediction p k+1 n+1 is then prescribed explicitly as a Neumann boundary condition. While this partitioned approach is easy to implement and may be incorporated into an existing FE solver package without difficulty, it may introduce inaccuracies during ejection phases and its convergence may deteriorate during isovolumetric phases [38]. Instabilities are related to the so-called balloon dilemma [62] and stem from the problem of estimating the change in pressure necessary to maintain the volume. Inherently, this requires to know the pressure-volume (pV ) relation of the cavity at this given point in time. However, this knowledge on chamber elastance is not available and thus iterative estimates are necessary to gradually inflate or deflate a cavity to its prescribed volume. As the elastance properties of the cavities are highly non-linear, an overestimation may induce oscillations and an underestimation may lead to very slow convergence and a punitively large numbers of Newton iterations.
A more elaborate approach is to treat p n+1 as an additional unknown in a monolithic scheme [33,63,38,64]. In addition to the equilibrium equations (3)(4) this requires one further equation for each cavity c ∈ C PDE . Using this approach, we get N cav = C PDE , the number of PDE cavities, additional equations of the form where V PDE c (u, t) is the cavity volume computed as the integral over the current surface Γ c,t , see Eq. (14), and V ODE c (p c , t) is the cavity volumes as predicted by the CircAdapt model for the intra-cavitary pressure p c , see Section 2.4 and Appendix A.
We write p C = [p c ] c∈C for the vector of up to 1 ≤ N cav ≤ 4 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, result in solving the  block system to find δu ∈ R 3N and δp C ∈ R Ncav such that with the updates Here, u k ∈ R 3N and p k C ∈ R Ncav are the solution vectors at the k-th Newton step. The block tangent stiffness matrix K is assembled according to Eqs. (B.20)-(B.24) and Eq. (C.12) and the right hand side vector R α according to Eqs. (B.25) and (B.26) and Eq. (C.9). The residual R p which measures the accuracy of the current coupling is the discrete version of (15), i.e., The whole procedure to perform the PDE-ODE coupling is given in Algorithm 1. In short, volume changes of the 3D cavities are driven by flow of blood over valves and outlets computed by the 0D model. In turn, updated pressures in the 3D cavities are used as an input to the lumped model in the next time step. Note that Eq. (16) is a block system with δp C holding at most four unknowns. Hence, we can apply a Schur complement approach for a small number of constraints, as described in Appendix D, to simplify the numerical solution of this linearized system, see Section 2.8. the cavitary volume of each chamber that is described as a 3D FE model has to be tracked as a function of time: V PDE c (x, t) for c ∈ C, While the described approach works for any combination of 3D PDE chambers and 0D ODE chambers, we consider biventricular FE models for our numerical examples in Section 3. This is, the ventricles are modeled as 3D PDEs as in Section 2.2 and the atria are modeled as 0D ODEs described by the CircAdapt model as in Appendix A.3. See Fig. 3 for a schematic of this 3D solid-0D fluid coupling.
Temporal synchronization of chamber contraction. In the lumped CircAdapt model contraction in individual chambers is controlled by prescribed trigger events. Based on the measured heart rate (HR) of 103 beats per minute contraction of the RA was triggered at intervals corresponding to a basic cycle length of 1/HR = 0.585 s. In all other chambers contraction was triggered by prescribed delays relative to the instant of contraction of the RA. In a hybrid coupled model contraction times used in 3D EM (10,11) and in the lumped CircAdapt model (A.8) must be synchronized accordingly. For this sake an interconnected event-driven finite-state machine (FSM) was used to control activation cycles in both 3D and 0D chamber models. Two types of FSMs were used, an autorhythmic FSM to generated triggers at a prescribed cycle length independently of any input, and a reactive excitable FSM of two possible states, excitable or non-excitable. The excitable FSM reacts to external trigger input. If the machine is in excitable state a transition is initiated to the non-excitable state, otherwise, if in non-excitable state, the FSM does not accept the input and remains in a non-excitable state. The FSM returns to its excitable state automatically after a prescribed effective refractory period. A transition from an excitable to a non-excitable state sends out a trigger event to all interconnected FSMs. These interconnections are implemented as delays representing the travel time needed for depolarization wave fronts to propagate to all neighboring interconnected FSMs.
The triggers provided by the FSM can be flexibly linked to entities within both 3D and 0D model. Specifically, the sino-atrial node is represented by an autorhythmic FSM that is directly interconnected to the RA. Both atrial cavities are implemented as a 0D model that initiate contraction in RA and LA based on FSM trigger events. The RA FSM connects to the atrial entrance to the atrio-ventricular node that transduces excitation through the atrio-ventricular (AV) node with a given delay. The ventricular exit of the AV node is connected to the left and right His bundle that trigger electrical activation of LV and RV. In the LV excitation is initiated by antero-septal, septal and posterior fascicle, and in the RV by a septal and a moderator band fascicle in the RV. As the timing of all fascicles was synchronous between all fascicles of a given chamber, fascicular timings were lumped together under RV and LV (see Fig. 4). RV and LV triggers prescribe fascicular activation times t a to the Eikonal equation that governs electrical activation of the ventricular cavities implemented as PDE model. Mechanical contraction of the ventricles was initiated then within a prescribed EM delay, t emd . The overall concept for synchronizing contraction in the coupled model is illustrated in Fig. 4 and FSM input parameters are given in Tab Table 1: FSM input parameters used for synchronizing electro-mechanical activity comprise heart rate (HR) or cycle length (CL), right atrium (RA), left atrium (LA), intra-atrial (AA), atrio-ventricular (AV) and inter-ventricular (VV) delays and effective refractory period (ERP).

Parameterization of the baseline model
For the sake of physiological validation the available experimental data were used to calibrate Update displacement u k+1 = u k + δu and cavity pressures p k+1

18:
Convergence test: Solution at time n + 1: u n+1 = u k+1 , p C,n+1 = p k+1  The sino-atrial node clock (SA) activates the RA at a prescribed cycle length, SA CL , starting at time SAt 0 . The LA initiates contraction with a delay of AA delay after the RA. The atrial entrance into the AV node activates at AVa which triggers, after the AV delay elapsed, the ventricular exit of the AV node that is connected to the His bundle at AVv. Fascicles in the LV are activated then relative to the LV trigger to initiate electrical propagation in the EP model. Similarly, the RV is activated with an interventricular delay of VV d before (VV d <0) or after VV d > 0 the LV.
the model in terms of stroke volume (SV) and peak systolic pressure in the LV (p LV ). Following [42], initial parameters of passive biomechanics, characterized by the material model given in (8), were taken from [65]; the model was unloaded using a backward displacement algorithm [66] and the material law's scaling parameter a was determined then by fitting the LV model to an empirical Klotz relation [67], using end-diastolic pressure (p ed ) and volume (V ed ) as input. Active stress model parameters τ c , S peak , τ r and duration of the force transient, t dur were determined as described previously [25]. Initial values and parameters for the CircAdapt model were chosen following [68].
To replicate the observed SV andp LV in the left ventricle, input parameters of the active stress model as well as CircAdapt model parameters were iteratively adjusted. The final parameterized model beating at 103 bpm produced a cardiac output of ≈2.1 L/ min with a SV of 21 mL, in keeping with the experimental data.

Physiological testing
The coupled 3D-0D model was subjected to thorough physiological testing by evaluating its transient response to alterations in loading conditions and contractile state. The model under baseline condition was used as a reference working point relative to which the effect of perturbations in loading and contractility was compared. Standard protocols for assessing of systolic and diastolic properties of the ventricles based on pV analysis [69] were implemented to qualitatively gauge the model's ability to consistently predict known cardiovascular physiology. For all perturbations in preload, afterload, or contractility, two points in time were considered, the immediate acute response after perturbing the system and the new approximate limit cycle reached after 8 beats. For baseline and each limit cycle, the end-systolic pressure volume relation (ESPVR) of the LV was interrogated by imposing additional step changes in afterload. For this sake, end-systolic pressure (p es ) and volume (V es ) were determined in the pV loops at the instant of end-systole, as determined by the cessation of flow out of the LV. Linear regression was used then to determine end-systolic elastance, E es , as the slope of the regression curve, and the volume intercept, V d , of the ESPVR.
Step changes in preload, afterload, and contractility were implemented by varying the cross-sectional area of the pulmonary veins, the systemic vascular resistance, R sys , and the active peak stress S peak generated by the myofilament model, respectively. Overall pump function was also assessed. Following [70], the heart as a pump can be described by the pump function graph (PFG), the relation between mean ventricular pressure, i.e., the ventricular pressure averaged over the entire cardiac cycle, and Cardiac Output. A PFG comprehensively describes cardiac pump function similar to the characterization of industrial pumps and ventricular assist devices. To construct a PFG, data were gathered under all protocols, including additional afterload variations over a wider range, between E a ≈ 0 by setting the system resistance R sys ≈ 0 and E a ≈ ∞, by closing the aortic valve, to obtain data points under extreme conditions corresponding to the LV beating in absence of external loading and under isovolumetric conditions, respectively.

Numerical framework
After discretization, at each Newton-Raphson step the block system (16) has to be solved. For this sake, we applied a Schur complement approach, see Appendix D, to cast the problem in a pure displacement formulation, to be able to reuse previously established solver methods [41]. In brief, we used the generalized minimal residual method (GMRES) with an relative error reduction of = 10 −8 . Efficient preconditioning was based on PETSc [71] and the incorporated solver suite hypre/BoomerAMG [72].
The dynamic version of the mechanics equations (3) was also used in other recent studies on cardiac EM [73,74,38] and showed advantages in performance of the linear solver -compared to the more common quasi-static approach -due to a more diagonal-dominant tangent stiffness matrix [75]. For the time integration we used a generalized-α scheme, see Appendix C, with spectral radius ρ ∞ = 0 and damping parameters β mass = 0.1 ms −1 , β stiff = 0.1 ms.
We implemented the coupling scheme in the FE framework Cardiac Arrhythmia Research Package (CARPentry) [76,57], built upon extensions of the openCARP EP framework (http: //www.opencarp.org). Based on the MATLAB code presented in [44], which is available on the CircAdapt website (http://www.circadapt.org), a C++ circulatory system module was implemented into CARPentry to achieve a computationally efficient and strongly scalable numerical scheme that allows fast simulation cycles.
Execution of the 3D-0D model was sped up by limiting the number of Newton steps to k max = 1 for the initial series of heart beats that were simulated to stabilize the coupled 3D-0D model to a limit cycle. This corresponds to a semi-implicit (linearly-implicit) discretization method [77] which worked very well in combination with the generalized-α scheme. Finally, after arriving at a stable limit cycle two further beats were simulated using a fully converging Newton method with k max = 20 and an relative 2 norm error reduction of the residual of = 10 −6 .

Parameterization of the baseline model
The coupled 3D-0D model was fit to approximate the experimental observed data on peak pressurep LV and stroke volume in the LV under baseline conditions. Electrical activation was driven by a tri-and bi-fascicular model in LV and RV, respectively, see Fig. 5(A). Conduction velocities were chosen for the given activation pattern to obtain a total ventricular activation time of ∼75 ms, compatible with the observed QRS duration of the ECG. Mechanical boundary conditions were set to limit radial contraction of the model, thus leading to a heart beat where ejection was mediated largely by atrio-ventricular plane displacement, i.e. long-axis shortening of the ventricles, and myocardial wall thickening. The resulting end-diastolic and end-systolic configuration of the model is shown in Fig. 5(A). Trains of 20 heart beats were simulated to arrive at an approximate stable limit cycle, as verified by inspecting the slope of the envelope of key hemodynamic state variables, see

Effect of spatial resolution
The impact of the relatively coarse spatial resolution (∼ 3.4 mm and ∼ 2.4 mm for LV and RV, respectively) used was evaluated first by repeating the baseline limit cycle protocol using a higher resolution mesh (∼ 1.3 mm and ∼ 1.2 mm for LV and RV, respectively). Both simulations used the exact same set of parameters and initial state vectors. With regard to pV behavior that is governed by global deformation of the ventricles, end-diastolic and end-systolic configurations were compared, see Fig. 6. Observed discrepancies between coarse and higher resolution model were marginal and well below the limits of experimental data uncertainty. Differences inp LV and SV were less than 5.4% and 6.9%, respectively, suggesting that the computationally efficient coarse model is suitable for performing a physiological validation study.

Physiological testing
The response of the coupled 3D-0D system to changes in afterload was probed by altering R sys in the range of ±65% around its nominal value of 6350 mmHg mL −1 ms −1 These maneuvers alter the slope of the arterial elastance curve, E a , pivoting E a around the point V ed and p = 0 in the pV diagramme. The initial response immediately after step changes in afterload and the new limit cycle are shown in Fig. 7(A) and (D), respectively.
The response to altering LV preload was then probed by stepwise reducing blood flow from the lungs into the LA by varying the cross sectional area of the pulmonary veins. Under such a walk down protocol the E a curve is shifted to the left towards smaller end-diastolic volumes V ed , without altering its slope, i.e. E a = p es /SV ≈ const, leading to lower p es . Stroke volumes under this protocol are assumed to gradually reduce due to the Frank-Starling mechanism, mediated by the length-dependence of active stress generation, S a (λ). Initial and limit cycle response under this preload perturbation protocol are shown in Fig. 7(B) and (E), respectively. As contractile  properties remained unchanged, the same slope E es of the ESPVR was obtained as before under step changes in afterload; compare estimated E es between Fig. 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/flattened the ESPVR, see Fig. 7(C) and (F). In the initial response V es and p es were affected, with V ed remaining constant, this led to a change in stroke volume and an apparent change in arterial elastance estimated by E a ≈ p es /SV, with ≈ −13.09/ + 36.81% relative to 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 given contractile state. Transient pV loops under this protocol are shown in Fig. 7(C) and (F). Increasing/decreasing preload shifts Ea curve and increases/decreases stroke volume via the Starling mechanism, mediated by the length-dependence of the active stress model. Determination of ESPVR was consistent with afterload protocol. (C) Increasing/decreasing contractility increases/decreases stroke volume, LV peak pressure and pes. For each contractile state afterload was also perturbed to determine end-systolic elastance Ees and V d .
of data within more extreme ranges of ventriculo-arterial coupling. The models' PFG is in agreement with known cardiovascular physiology, see Fig. 8. Keeping contactility and preload constant, the PFG with mean ventricular pressure (MVP), as a function of flow or stroke volume can be approximated by a quadratic function, MVP(q) =P iso 1 − (q/q mx ) 2 , withP iso and q mx being the maximum MVP and maximum flow under isometric and unloaded conditions, respectively. The MVP as a function of time is given by the integral with a constant cycle length t cycle . Increasing preload shifts the PFG towards higher flows and pressures. Increasing/decreasing contractility pivots the PFG, leading to a steeper/flatter slope ∆MVP/∆q of the PFG.

Numerical performance
Computational times for a single heart beat of the lower and higher resolution baseline models are given in Tab. 3. Simulations were performed on the Vienna Scientific Cluster (VSC4) and we distinguish between solver-time, t s , which is the accummulated GMRES solver time over all loading/time steps; and assembly-time, t a , which is the time spent on the setup of boundary conditions and on the assembly of matrices and vectors of the linearized system (16). In total, for a full simulation with loading, 18 initialization beats, and 2 final beats with a fully converging    Table 3: Summary of numerical metrics for coarse and fine model. Given are the number of compute cores used on VSC4; the average spatial resolution in LV and RV,h LV andh RV ; the number of elements and nodes spanning the mesh; as well as solver, assembly, and total times for a single Newton iteration (t s,1 , t a,1 ) and a fully converged Newton solution (ts,c, ta,c), and the total simulation time per heart beat for single iteration and fully converged Newton scenarios, T b,1 and T b,c , respectively. Timings refer to a single heart beat lasting 0.585 s at a time step size of 1 ms. In addition the cumulated solver (t s,ld ) and assembly (t a,ld ) times for the loading phase using 32 load steps are presented.

Discussion
In this study, we report on the development of a monolithic 3D solid -0D fluid coupling method that allows to flexibly combine 3D PDE and non-linear 0D EM representations of cardiac cavities. The hybrid 3D-0D model of the heart was coupled to a 0D closed-loop model of the cardiovascular system, where all 0D components were based on the CircAdapt model [44]. The combined model can be set up to represent one, two, three, or all four cavities as 3D PDE model and all other elements of the cardiovascular system as 0D models based on CircAdapt. In this study feasibility of this approach is demonstrated by coupling a 3D PDE model of bi-ventricular EM to 0D model representations of atrial EM function and circulation based on the CircAdapt model. The combined model was parameterized under baseline conditions and subjected to comprehensive physiological testing to demonstrate the model's ability to correctly predict known physiological behaviors. A broad range of experimental protocols for altering loading conditions and contractility were simulated to interrogate the models' transient responses to these maneuvers [69]. Overall, pV analyses of the hemodynamic model output showed close agreement with established knowledge on cardiovascular physiology. The underlying numerical scheme is also represented in detail, including a comprehensive mathematical representation of the CircAdapt model. Robustness -in terms of stability 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 observation 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 or therapies altering loading conditions and contractility that are valid within a commonly accepted physiological reference frame.

Physiological validation
Predictive modeling applications critically rely on the ability of models to encapsulate the most relevant mechanisms governing the cardiovascular response to a given intervention that alters loading 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.
Validation aimed at replicating, overall, known well established behaviors and not at a 1:1 validation against experimental data. For this sake, experimental standard protocols for altering afterload, preload and, contractility were applied to the stabilized baseline model to study its response. Two scenarios were analyzed, the initial acute response to a step change in a single parameter -afterload, preload, or contractility -where effects on other unaltered parameters were minimal, see Fig. 7(A)-(C), and, the limit cycle response observed after a number of beats where transients have largely subsided, but indirect effects led to alteration of other parameters due to the systemic inter-dependencies between these, see Fig. 7

(D)-(F).
Afterload was altered by varying the systemic resistance R sys and, to mimick more extreme conditions closer to isometric contraction, the resistance of the aortic valve. Increasing/decreasing afterload is reflected in pivoting the arterial elastance curve E a around the point of end-diastolic volume and zero pressure. This behavior is illustrated in Fig. 7 and for the more extensive protocols used 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 regulatory mechanisms, heart rate remained unaltered. Thus, E a ≈ MAP/(SV · HR) ≈ R sys holds. The phenomenological active stress model as given in Eqs. (10)-(11) accounted for length-dependent tension, but not for velocity dependence. Thus, afterload effects on the velocity of fiber shortening were ignored.
Altering V ed by changing the cross section of the pulmonary vein orifices to increase/reduce preload increased/reduced SV. This was due to the Frank-Starling mechanism, as represented by the length-tension relationship of the active stress model in Eq. (11). In the immediate response to a step change in preload the LV emptied to almost the same V es with only minor deviations due to changes in arterial pressure induced by the change in SV. After transients subsided, ESPVR and arterial elastance were the same as under baseline conditions, with E a being shifted according to the changed V ed . In the new limit cycle p es and MAP were increased, but changes in SV were rather small. This could be attributed to a rather flat slope of the Frank-Starling curve SV(p ed ) (not shown) that is related to the significant heterogeneity in fiber stretch in end-diastolic state. Such heterogeneity inevitably arises in biventricular EM models that do not account for residual strains in the unpressurized configuration. As fiber stretch in the reference configuration 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 λ in our simulations increased with p ed the increasing spread of λ around its mean led to low fiber stretch λ < λ 0 in various regions, particularly in antero-septal and postero-septal segments of the LV. These regions contributed increasingly less to contraction with increasing p ed due to the length-dependence of S a (λ) ≈ 0, 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/decrease in SV and systolic pressures in terms ofp and p es . Correspondingly, the ESPVR, as sampled by perturbing LV afterload, was steepened/flattened as expected, see Fig. 7(C). In the initial response the slope of E a was affected, but after stabilization E a was the same for all contractile states. For the given parameterization ventriculo-arterial coupling E es /E a fell outside the optimal range of 1-2, where external work is maximized around E es /E a ≈ 1 and optimal efficiency is achieved at E es /E a ≈ 2. Sincep and SV and as such E a ≈ p es /SV were used to parameterize the model the resulting slope of the ESPVR was to flat for the LV to operate within this optimal range. Reasons are multifactorial and include the absence of velocity-dependence and, potentially, also the role of mechanical boundary conditions. The main culprit to blame for the limited slope of E es is the impairment of the Frank-Starling mechanisms due to fiber stretch heterogeneity. Our model deviates here from the general physiology-based assumption of uniform fiber-stretch to be in place in an end-diastolic state [78,79,80,81].
The constructed PFG was qualitatively in keeping with physiological expectations, see Fig. 8. Data points obtained from varying afterload between close to unloaded and isometric conditions agreed well with an assumed quadratic relation between MVP and flow or SV. Increasing preload shifted the PFG up and left towards higher flows and pressures, with the opposite trend being observed for decreasing preload. Altering contractility rotated the PFG.

Numerical aspects
The computational cost imposed by higher resolution EM models requires efficient numerical solvers. Strong scaling characteristics of our numerical framework was reported in detail previously [41,82]. The compute times reported in Section 3.4 indicate that setup and assembly time were the dominating factors during the initial passive inflation (loading) phase while for the subsequent coupled 3D-0D EM simulations of a heart beat solver time was the predominant part of total CPU time. This is due to the Schur complement, see Appendix D, that needs to be solved during the 3D-0D active EM phase. This involved -in our case of a bi-ventricular model comprising two PDE cavities -three applications of the GMRES solver while matrices were only assembled once per Newton step.
Further, significant savings in compute time could be achieved using a semi-implicit approach, see Section 2.6, until a limit cycle was reached. A heartbeat using a fully converging Newton was about five times more 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 cycle protocol were then computed using a fully converging Newton method this gain in performance had no quantitative impact on any of the primary simulation outcomes.
Reducing spatial resolution down to 2.4 millimeter and 3.4 millimeter in RV and LV, respectively, allowed for multibeat simulations within tractable time frames using a desktop computer. The impact of this reduction on the hemodynamic output variables was very minor, particularly when viewed in context of the significant observational uncertainties the type of measured data used for model calibration are afflicted with. While probing only two resolutions by far cannot be considered a rigorous convergence study, our results suggest that the use of relatively coarse meshes, with resolutions in the range between 2 to 4 mm are adequate. Indeed, this appears to fall within the range of resolutions used in other recent cardiac mechanics simulation studies. For instance, meshes comprising 208 561 [83] and 167 323 [74] tetrahedral elements were used to represent human-sized hearts with all four chambers whereas our coarse model used 45 686 of the same elements only for the ventricles of a much smaller sized canine heart. While our study cannot offer any conclusive recommendations on choosing an appropriate spatial resolution our results suggest that for solving cardiac mechanics problems significantly less spatial resolution is necessary to achieve acceptable accuracy [84] relative to solving EP problems where spatial resolution is critical to resolve the steep propagating depolarization wavefronts [59]. The reaction-Eikonal model we used to represent electrical activation and repolarization is not sensitive to spatial resolution, as shown previously [57], and yields accurate activation patterns on coarse meshes, allowing to use the same grid for both EP and mechanics, without the need for projection of data between the physical grids.
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 Ryzen Threadripper 2990X), demonstrating that realistic multi-beat simulations of the presented 3D-0D cardiac models deliver sufficient performance for advanced physiological simulation scenarios, even on a small number of cores. This is of paramount importance for future parameterization studies where numerous simulations have to be carried out to personalize models to patient data. With around 2 and a half hours on 256 cores for 20 beats also the higher resolution model could be executed within a tractable time frame.

Relation to previous work
The holistic framework described in this study constitutes a major step towards a universal cardiac electro-mechanics simulation engine that can be applied, in principle and after appropriate parameterization, to a very broad range of applications. Our study builds on and further advances various concepts that have been reported previously in a number of excellent studies [35,33,37,38]. While coupling 3D PDE models to a closed loop circulatory system is important to ensure consistency, by allowing blood to redistribute between compliances in the system, only a few have been reported [37,33,38,40]. However, these were limited in some of the following regards which restricted their universal applicability. For instance, models were discretized with a small number of cubic Hermite elements which led to anatomically stylized representations of the ventricles [85,33,26,86], with an artificially chopped base and a hole or collapsed elements [87] in the apex, owing to their limited ability to accurately accommodate more complex anatomical shapes without greatly increasing computational times [88]. Often, artificial boundary conditions were used that fixed the motion of the base [37,89] and, thus, enforced a zero atrio-ventricular plane displacement. These studies, with only a few exceptions [35,74,90], were unable to replicate a physiological kinematics characterized by the reciprocal filling properties of the heart by maintaining a constant pericardial shape. In other studies EP was not modeled at all, assuming that contraction in the ventricles is initiated simultaneously [38,40,74], or non-physiological activation sequences were used to trigger contraction, both of which impair length-dependent tension mechanisms [37]. Computational cost of numerical methods is not addressed in most previous works; notable exceptions include [37,90,38,74,89] where compute times for one heart beat range between 1.8 and 24 hours. This is considerably less efficient compared to compute times presented in this study in Section 3.4. Mostly simplified circuit models were used [37,33,38,40] for simulating a single heart beat or the simulated PV loops featured unphysiological edgy morphologies due to the usage of (i) simple diode valve models not accounting for inertia, Bernoulli effects and flow resistance between compartments and (ii) time-varying elastance models for compliances such as the atria that yield fixed pressure-volume relations and cannot contract in a load-dependent manner as the sarcomere-based contractile 0D chambers used in the CircAdapt model. None of these methodological limitations apply to the modeling framework presented in here.

Limitations
Computational modeling relies on assumptions and approximations, especially multiphysics simulations on organ scale level as presented in this study. While the accuracy of most individual model components was already assessed in previous studies, all these components have limitations and were chosen (i) to achieve great computational performance while preserving almost the full biophysical details of possibly more accurate, computationally costly models; and (ii) to ease parameterization, i.e., to achieve physiological results even with a smaller amount of parameters compared to more accurate but also more complex models. Dependent on the clinical application it might be required to consider other individual components than those chosen in this study.
Efforts to parameterize the combined model were limited, only a small subset of available data were used for model fitting. A generic five-fascicular representation of the cardiac conduction system was used to drive ventricular activation, without attempting to match the recorded ECG [91]. More comprehensive and efficient procedures are required to further enhance compatibility of the high dimensional combined 3D-0D model with all available observations. Building on strategies for fitting the standalone CircAdapt model to hemodynamic data [92,93], a hybrid parameterization approach appears a pragmatic solution where the 0D model is fitted first to observed data and in a subsequent parameterization step the 3D cavities are fit to the pV characteristics of the corresponding 0D cavities. However, for the sake of demonstrating overall compatibility with known cardiovascular physiology a high fidelity match with experimental data, while desirable, is not crucial. Future more advanced applications that attempt to predict therapeutic responses in a patient-specific manner will critically depend, beyond a comprehensive representation of the most relevant therapy mechanisms, on the fidelity of personalization. Given the large number of model parameters this is not a trivial task that will require the development of dedicated parameter identification strategies. In this regard the outstanding computational performance of the model is essential to facilitate a detailed model personalization which requires a large amount of forward simulations.
While the framework used in this study can be considered universal and feature complete, two important aspects remained unaccounted for. First, active stresses generated by the phenomenological model featured length-but no velocity dependence. However, a velocity-dependent term could be incorporated, as methods that avoid numerical instabilities due to velocity dependence have been reported [73], or a biophysically highly detailed model of excitation-contraction coupling could be used instead as in our previous work [41]. Secondly, residual stresses in the unpressurized ventricles were ignored which impairs, to some extent, the length-dependent Starling effect.
Inflation of the unpressurized configuration to p ed without considering residual stresses, inevitably introduces a significant spatial heterogeneity in fiber stretch. This is in contrast to the common assumption of homogeneous fiber stretch in the end-diastolic state [78,79,80,81]. This was reflected in a rather flat slope of the Starling curve SV(p ed ). As shown in Fig. 8C, for the given inotropic state and afterload, the Starling curve was close to linear, with a slope of ≈ 1.3 mL/mmHg. In humans, slopes in the range between ≈ 2.7 − 5.5mL/mmHg were measured for non-athletes and athletes, respectively [94].
For the sake of saving computational costs, a coarse spatial resolution was used for discretizing the bi-ventricular model. Thus, models were lightweight enough to carry out the larger number of simulations that were needed for parameterization, the determination of a limit cycle as well as the fine grained testing of physiological maneuvers. The use of such coarse spatial resolutions introduces inaccuracies with regard to fibers and sheet arrangements which are defined on a per element basis. As parts of the model such as the RV wall were composed only of two element layers, transmural fiber rotation was essentially reduced to two fiber families. These were mostly aligned with the endocardial and epicardial fiber orientations as prescribed on a per rule basis, see Fig. 1. Thus, the model's predictions of motion, strains and stresses may deviate quantitatively from models discretized at higher resolution. However, comparing between two models of different resolution revealed that quantitative differences were minor and, qualitatively, models showed essentially the same behavior. Most differences stemmed from differences in stiffness of the simple P1-P0 element types used which tended to be more compliant for higher resolutions. Nonetheless, errors in simulated hemodynamic outcomes were small enough to be considered negligible when put in context to observational uncertainties clinical or experimental data are afflicted with.
Finally, the presented model is not validated rigorously against experimental data. Typically, an independent validation is performed by the comparison of displacement [95,74,96] and/or strain [55,97,40] to observations from cine MRI or 3D tagged MRI data. However, an accurate cardiac motion and deformation field can be obtained by tuning boundary conditions and in vivo MRI strain measurements have major caveats [98,99]. Further, clinical validation can be conducted by comparing simulated ECG, pressure, and volume traces and derived quantities, e.g., SV or ejection fraction, against clinically measured data [25,97,90]. However, commonly, these measurements are all used for model calibration to optimize goodness of fit of simulated outputs to the data and, thus, these cannot be used for the purpose of model validation. Overall, difficulty of proper model validation remains a point of concern and in many cardiac modeling studies validation against actual patient data is limited [100]. In this work, we focused on replicating known, well-established behaviors to show and prove a physiological predictive power under a broad range of experimental protocols. Together with a rigorous, independent validation against image data in future studies, this will set a new standard in cardiac modeling.

Conclusions
This study reports on a flexible monolithic 3D solid -0D fluid coupling method for integrated models of cardiac EM and cardiovascular hemodynamics. Feasibility of the approach is demonstrated by coupling a 3D PDE model of bi-ventricular EM to 0D model representations of atrial EM and circulation based on the CircAdapt model. The coupled 3D-0D model is shown to be robust, computational efficient and able to correctly replicate known physiological behaviors in response to experimental protocols for assessing systolic and diastolic ventricular properties based on pV analysis. These features combined render advanced EM modeling applications feasible. The model facilitates the exploration of parameter spaces over prolonged observation periods 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 relationships that could have appeared to influence the work reported in this paper.

Appendix A. CircAdapt equations summary
Appendix A. 1

. CircAdapt Tube Module
The tube module represents the entrance of a compliant blood vessel capable of propagating a pressure-flow wave component added to a constant flow, see [27]. Vessels directly attached to the heart, aorta (AO), arteria pulmonalis (AP), venae cavae (VC), and venae pulmonales (VP) are modeled in 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 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 [101], gives the average extension, λ t , of the fibers in the wall by the wall cross-sectional area and V wall t = A wall t l t the wall volume. Cavity pressure depends on λ t p t = σ t (λ t ) λ −3 t , with the mean Cauchy fiber stress σ t that is modeled by the constitutive equation see Arts et al. [101]. Here, k t a stiffness exponent; 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 computed as Finally, the basic relation between characteristic wave impedance Z t compliance C t , and inertance I t 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 Hill model, see [61,44] where the constant L elast,iso is the length of the series elastic element during isovolumetric contraction 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 formation within the fibers of the patch, iṡ In Eq. (A.8) sarcomere contractility rises according to the function f rise c , a phenomenological representation of the rate of cross bridge formation within the patch, where L pas0,ref is sarcomere length with zero passive stress and L s c the total sarcomere length, see above. Using that we compute σ fib,tit with σ act,max the maximal isometric stress and the constant exponent The ECM is modeled as being stiffer than the myocyte contribution using where σ pas,max is an empirical parameter. A mid,tot

Sarcomere total stress
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 with A mid,ref c the surface area in the reference state, see [44]. Note that this updated fiber strain is used 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 with l c the long-axis length of the cavity. The characteristic wave impedance Z c is approximated according to (A.3), see also [27], and by applying the chain rule with the sheet stiffness In case that one ODE and one PDE cavity is included in the model, ventricles are modeled as atria above. Otherwise, ventricular and septal midwall volumes are modeled as a ventricular composite [61] which is defined by the common radius y mid of the wall junction and the enclosed midwall cap volumes, see Here Midwall area and curvature are consequently computed for c ∈ {LV, RV, Sep}, 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 It is required that that the total midwall tension at junctions is zero, i.e., Eq. (A.24) is solved by an iterative Newton scheme where τ Sep is a time constant. Consequently, the values for the tensions discussed above are updated and the scheme is iterated until convergence. Midwall volumes are updated by long-axis length l c and cross-sectional area A c and of the cavity are computed by Assuming the pressure surrounding the ventricular composite to be zero, internal chamber pressure for the ventricles is now found as Appendix A.5. Pericardial mechanics The four cardiac chambers are supposed to have an additional pressure component due to the pericardium. Pressure p peri exerted by the pericardial sack on atria and ventricles was computed as a non-linear function of pericardial volume V peri , computed as the sum of blood pool and wall volumes of the four cardiac chambers: where p ref peri and V ref peri are constant reference pressure and volume, respectively, and k peri defines the degree of non-linearity of the pressure-volume relation.
Cavity pressures are updated according to 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 the nonlinearity of the arteriovenous resistances, see Tabs Computation of time derivative of flow across valves and venous-atrial inlet requires in input the cross-sectional area of proximal and distal elements to the channel.
The pressure drop (∆p v ) across a valve is the sum of the effects of inertia due to acceleration in time and the Bernoulli effect, see [102] ∆p where ρ b is the density of blood, A v is the current cross-sectional area of the valve, and l v is the length of the channel with inertia. If not mentioned otherwise this values is estimated as . For q v < 0 which indicates that the valve is leaking v in v is the velocity distal to the valve v dist v and the outflow velocity is the maximum of the blood velocities in the valve With forward pressure drop, the valve opens immediately. With backward pressure and forward flow, the valve is closing softly by a continuous function where A leak v is the given valve cross-sectional areas of the closed (regurging) valve. Using this the current cross sectional area of the valve is Figure Flow over the valve is finally updated using (A.37) bẏ Appendix A.9. Solve ODE system A Runge-Kutta-Fehlberg method (RKF45), see, e.g., [104], is used to solve the system of 26 ordinary differential equations (ODEs): 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 contracting length 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 boundary value problem (3)(4) and (15). For the full nonlinear elastodynamics problem see Appendix C. The stationary boundary value problem is formally equivalent to the equations   with the second Piola-Kirchhoff stress tensor S, see (5), and the directional derivative of the Green-Lagrange strain tensor Σ(u, v), see [41,105]. The weak form of the contribution of pressure loads (B.1), right term, is computed using (4) The first term of the coupling equation (B.2) is computed from (14) using Nanson's formula and The second term of (B.

Appendix B.2. Consistent Linearization
To solve the nonlinear variational equations (B.1)-(B.2), with a FE approach we first apply a Newton-Raphson scheme, see [77]. Given a nonlinear and continuously differentiable operator F : X → Y a solution to F (x) = 0 can be approximated by which is looped until convergence. In our case, we have X = H 1 (Ω 0 , Γ 0,D ) 3 × R, Y = R 2 , ∆x = (∆u, ∆p c ) , x k = (u k , p k c ) , and F = (R u , R p ) . We obtain the following linearized saddlepoint problem for each (u k , p k c ) ∈ H 1 (Ω 0 , Γ 0,D ) with the updates u k+1 = u k + ∆u, The Gâteaux derivative of (B.1) with respect to the pressure change update ∆p c yields the third term of (B.6) The residual R u , i.e., the right hand side of (B.6), is computed as Grad ∆u)x · F − k n out 0 q ds X − 1 3 Γ0,N J k x · F − k (Grad ∆u) F − k n out 0 q ds X + 1 3 Γ0,N J k ∆u · F − k n out 0 q ds X , (B.14) with q a testfunction that is 1 for the surface of cavity c, Γ 0,c , and 0 otherwise. The second term of (B.7) is computed as a numerical derivative ∆p c , V ODE element space 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 3N and δp C ∈ R Ncav such that i.e., , (B.17) with the solution vectors u k ∈ R 3N and p k C ∈ R Ncav at the k-th Newton step. The tangent stiffness matrix A ∈ R 3N ×3N is calculated from (B.10) according to see also [41,105]. The off-diagonal matrices B u ∈ R 3N ×Ncav and B p ∈ R Ncav×3N in (B.17) are assembled using (B.14) with the constant shape functionφ j = 1 if τ j ∈ Γ 0,c andφ j = 0 if τ j / ∈ Γ 0,c for c ∈ {LV, RV, LA, RA}. Using a technique as described in [64,Sect. 4.2] this assembling procedure can be simplified for closed cavities such that B p (u k , p k C ) = B u (u k , p k C ) .
The circulatory compliance matrix C (p k C ) ∈ R Ncav×Ncav is computed from (B.15) as with the constant shape functionφ i ,φ j = 1 for cavity c and 0 otherwise, leading to a diagonal matrix.
The terms on the upper right hand side A ∈ R 3N , B p ∈ R 3N 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 From this and the relationship by Newmark [109] u n+1 = u n + ∆t γu n+1 + (1 − γ)u n , v n+1 = v n + ∆t γv n+1 + (1 − γ)v n , 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 we compute To calculate the solution at the current timestep we assume that we know u n ,u n , v n andv n from the previous time step n and get from Eq. (C.4) for the residual R α (u k n+1 ) := −ρ 0 M αv k n+αm − R u (u k n+α f ), (C.9) withv k n+αm :=v n+αm (u k n+1 ) and u k n+α f := u n+α f (u k n+1 ). To increase stability we consider Rayleigh damping by adding the two matrices D mass (u k n+1 ) = ρ 0 β mass M α v k n+α f , (C.10) to the residual (C.9) with v k n+α f := v n+α f (u k n+1 ) and Rayleigh damping parameters β mass ≥ 0 ms −1 , β stiff ≥ 0 ms.
The tangent stiffness matrix is now calculated using (C.8) as The realization of (D.2) involves m + 1 solves and the inversion of an m × m matrix. Since m is generally small this can be done symbolically.