The Impact of Standard Ablation Strategies for Atrial Fibrillation on Cardiovascular Performance in a Four-Chamber Heart Model

Purpose Atrial fibrillation is one of the most frequent cardiac arrhythmias in the industrialized world and ablation therapy is the method of choice for many patients. However, ablation scars alter the electrophysiological activation and the mechanical behavior of the affected atria. Different ablation strategies with the aim to terminate atrial fibrillation and prevent its recurrence exist but their impact on the performance of the heart is often neglected. Methods In this work, we present a simulation study analyzing five commonly used ablation scar patterns and their combinations in the left atrium regarding their impact on the pumping function of the heart using an electromechanical whole-heart model. We analyzed how the altered atrial activation and increased stiffness due to the ablation scars affect atrial as well as ventricular contraction and relaxation. Results We found that systolic and diastolic function of the left atrium is impaired by ablation scars and that the reduction of atrial stroke volume of up to 11.43% depends linearly on the amount of inactivated tissue. Consequently, the end-diastolic volume of the left ventricle, and thus stroke volume, was reduced by up to 1.4 and 1.8%, respectively. During ventricular systole, left atrial pressure was increased by up to 20% due to changes in the atrial activation sequence and the stiffening of scar tissue. Conclusion This study provides biomechanical evidence that atrial ablation has acute effects not only on atrial contraction but also on ventricular performance. Therefore, the position and extent of ablation scars is not only important for the termination of arrhythmias but is also determining long-term pumping efficiency. If confirmed in larger cohorts, these results have the potential to help tailoring ablation strategies towards minimal global cardiovascular impairment.


INTRODUCTION
Although atrial fibrillation (AF) can often be managed pharmacologically, catheter ablation is now considered as a first line treatment for paroxysmal AF patients. 28 Unfortunately, the long-term success rates of catheter ablation of AF are not satisfactory yet. One reason could be long term adaptations of the atria to compensate for ablation-induced impairment of atrial function. Thomas et al. 57 analyzed the effect of atrial radiofrequency ablation (RFA) on the atrial mechanical function. They concluded that multiple linear RFA lesions may impair atrial contractility and that the reduced atrial function is partly due to loss of viable myocardial tissue but may also be caused by the altered atrial activation. Different ablation strategies with the aim to prevent AF and its recurrence exist but their respective impacts on the ventricular mechanics and hemodynamics have not been studied systematically. In this context, cardiac modeling opens up the possibility to get a better understanding of how RFA affects the mechanical behavior and the pumping function of the heart by applying different treatment strategies to the same virtual patient under controlled conditions, which is not feasible in a clinical context.
While in the past most attention in the field of cardiac modeling had been paid to the simulation of the ventricles, and still is an active field of research, [4][5][6]27 modeling of the atria has moved into focus in the last years. However, most publications deal with atrial electrophysiology and there are only a few addressing atrial mechanics. 16 Jernigan et al. 30 published a study on the mechanical properties of porcine left atrium, which were assessed with uniaxial tests. In 2011, Di Martino et al. 14 presented a computational model of the porcine left atrium, which allowed to analyze the wall stress due to the mitral valve movement, obtained from raw multi-detector computed tomography data. The mechanical properties were modeled based on biaxial experiments with porcine atrial tissue. In the same year, Di Martino et al. 15 analyzed how ventricular tachypacing affects the spatial and temporal stress distribution of the left atrial wall. Bellini et al. 9 provided a comprehensive characterization of the passive biomechanics of the left human atria based on a Fung-type elastic strain energy potential. These publications all deal with the passive mechanical behavior of the atria. Models for the active contraction of the atria have only been published recently 3,35 and the simulation of whole heart models is becoming more feasible. 22 A recent study by Ho¨rmann et al. 29 modeled atrial systole using standard ablation strategies in a non-pathological and a dilated atrium to study the impact of ablation lesions on the mechanical performance of the atria. Their model was able to detect differences in left atrial contractility and ejection fraction for multiple activation sequences resulting from RFA and atrial fibrillation-induced atrial remodelling represented by a reduction in global conduction velocity. Since they were only interested in atrial systole, the ventricles were only represented by a passive stiffness and they only included a Windkessel to model the hemodynamics of the atria, thus neglecting a change in pre-and afterload over multiple heart beats. Another study that dealt with atrial mechanics was published by Phung et al. 48 Unlike Ho¨rmann et al., 29 they connected the left atrium to a circulation model for more realistic hemodynamic conditions and validated their baseline model against pre-treatment data by means of left atrial wall motion. Phung et al. found that ablating the posterior wall had a smaller impact on atrial function compared to other ablation patterns such as pulmonary vein isolation or wide area circumferential ablation. They explained this with the little movement of this region in the pretreatment model and thus small contribution to atrial function. A limitation of the study by Phung et al. is the missing representation of electrical propagation in the model.
We hypothesize that ablation scars in the atria also have an effect on the ventricles. During ventricular contraction, the atrioventricular plane is pulled towards the apex, the atria are stretched and their volume is increased. This mechanism supports the filling of the atria with blood from the venae cavae and the pulmonary veins. Eventually, this blood volume is available for ventricular filling during ventricular relaxation. An increased stiffness of the atria due to ablation scars may impede atrioventricular plane displacement(AVPD) and consequently also atrial filling. Accordingly, two effects may affect the ventricular filling negatively: on the one hand, a reduced filling capacity of the atria during the ventricular systole and therefore less blood volume available for the ventricular filling. On the other hand, a reduced active contribution of the atrium due to an altered atrial activation induced by the ablation scars. Here, the effect certainly depends on the respective ablation pattern and the amount of ablated tissue. Alhogbani et al. 1 assessed the contribution of the left atrium to the filling of the left ventricle using cardiac magnetic resonance imaging (MRI) for 120 normal subjects. They reported that the contribution of the atria to ventricular filling increases with age and is in the range of 10-40%. Considering that about 70% of the patients with AF are between 65 and 85 years old, 20 the implications of RFA on the atrial contraction should not be disregarded. To overcome the limitations of the studies by Ho¨rmann et al., 29,48 we utilize our previously published four-chamber heart model including a closed-loop circulatory system, electrophysiological wave propagation and myocardial contraction 22 to extend a previous study from our group 10 by analyzing the impact of five commonly used ablation lesions on cardiovascular performance. Furthermore, we analyze how the stiffening of scar tissue due to an accumulation of collagen in the myocardium affects the deformation of the heart and in particular AVPD. The sensitivity of the model towards changes in the stiffness of the scars is evaluated for one of the scar patterns. Finally, we introduce the global effects of AF related fibrosis by reducing the diffusion coefficient.

Whole Heart Finite Element Model
The cardiac anatomy was manually segmented from MRI data (0:72 Â 0:72 Â 1:8 mm 3 ) of a 33 year old male volunteer. The volunteer provided informed consent and the study was approved by the IRB of Heidelberg University Hospital. 19 The MRI data were acquired using a 1.5 T MR tomography system and consist of a static whole heart image stack taken during diastasis as well as time-resolved images in several long and short axis slices. Based on the segmentation, we first labeled the atria (X A ¼ X LA [ X RA ) and the ventricles (X V ¼ X LV [ X RV ). As shown in Ref. [22], we extended the geometry with a simplified volumetric representation of the mitral valve, the tricuspid valve, the aortic valve and the pulmonary valve (X Valves ) to close the cavities. We found that this helps to improve the movement of the valve plane compared to just using surface elements, since volumetric elements can resist the pressure exerted onto them. Additionally, we closed the endo-and epicardial surfaces of the atria and added truncated pulmonary veins, vena cavae as well as the ascending aorta and pulmonary artery (X Vessels ). Compared to the study in Ref. [22], no adipose tissue in the atrioventricular and interventricular grooves was included due to numerical issues caused by locally insufficient element quality. Furthermore, we added a concentric layer of tissue around the entire heart (X Peri ) which phenomenologically represents the influence of the pericardium. 19 Two tetrahedral meshes were created using Gmsh 26 : with 128,976 elements and (2) the electrophysiological reference domain X EP ¼ X V [ X A as a subset of X M with 50,058,295 elements. In case of the mechanical mesh, we made sure that the model featured a minimum of two elements throughout the wall of the heart. We used rule-based methods to assign the myofiber orientation on X EP as a local basis Q ¼ ðf 0 ; s 0 ; n 0 Þ for the atria 60,61 and the ventricles. 8,50 The fiber angle in the ventricles was chosen as 60 and À60 on the endocardial and epicardial surface, respectively. The sheet angle was set to À65 on the endocardium and 25 on the epicardium. 8 For the mechanical mesh, we interpolated the myocyte orientation from the vertices of X EP to the quadrature points of X M .

Electromechanical PDE Model
The electromechanical model used in this study to simulate a four-chamber model of a human heart comprises three main components. First, we have the electrical propagation model driven by cellular electrophysiology. Second, we describe the solid mechanics model that is used to simulate the active and passive forces in the myocardium. Last, the mechanical model is extended by boundary conditions, which originate from physiologically important mechanisms such as the circulatory system or the tissue surrounding the heart. A complete description of the fully coupled electromechanical model was recently given in Gerach et al. 22 Since the focus of this study is on mechanical features and not on the subtleties of mechanoelectric feedback, we only implement excitation-contraction coupling mechanisms and forego the effects of deformation on electrical propagation.

Electrical Propagation Model
The electrical propagation in the myocardium is modeled by the monodomain equation on the reference domain X EP Â ð0; T: with the transmembrane voltage V m , the cellular surface-to-volume ratio b, the membrane capacitance C m , the total ionic current I ion , and an externally applied stimulus current I ext . The anisotropic conductivity tensor depends on the fiber f 0 , sheet s 0 , and sheet-normal n 0 orientation of the myocardium in the reference configuration with generally higher conductivity values in the fiber direction r f ! r s ! r n ! 0. We set the conductivity in the atria such that the total atrial activation takes approximately 100 ms. 44 For the ventricles, we tuned the conductivities according to Mendonca Costa et al. 40 to match conduction velocities of 0.6, 0.4, 0.2 m/s in fiber, sheet, and sheet-normal directions. 5 The transmembrane voltage V m is coupled to a system of ordinary differential equations, where w describes a state vector of different gating mechanisms and c describes the vector of intracellular ion concentrations. To account for the inter-chamber differences in electrical activity, we use the model of Courtemanche et al. 13 in the atria and the model of O'Hara et al. 45 in the ventricles. The O'Hara et al. model was further modified with the changes proposed by Refs. [46] and [18] to ensure robust conduction. For the implementation and parameterization of these models, we refer to the original publications.
The external stimulus I ext ðx; tÞ is applied at the junction of the right atrial appendage and the superior vena cava 37 to initiate the depolarization wave in the atria. In the ventricles, we use consistent biventricular coordinates 51 to define seven positions ( Table 1) that are known to be common sites of earliest activation. 12,17 To account for the conduction delay introduced by the atrioventricular node, the ventricles are stimulated 160 ms after the atria. Additionally, a thin, fast-conducting endocardial layer is defined in the ventricles to mimic the His-Purkinje system. We choose the conduction velocity in this layer to be twice as high as in the ventricular bulk tissue.
Finally, the electrophysiological model is complemented by homogeneous Neumann boundary conditions on the entire surface C EP . All parameters of the electrical propagation model are given in Table 4.
Equation (1) is solved numerically using linear tetrahedral finite elements in space and a first order operator splitting in time with explicit integration for the reaction term and implicit integration for the diffusion term (Dt EP ¼ 10 ls) as described in Ref. [22].

Active and Passive Mechanics
The biomechanical behavior of the heart is described by the balance of linear momentum equation in its total Lagrangian formulation where the deformation tensor F ¼ I þ $ X u describes the deformation from the reference configuration X M ¼ X 0 ðXÞ into the current configuration X t ðxÞ, u is the displacement vector, q 0 the mass density, and S is the second Piola-Kirchhoff stress. To account for the active stress generated during cardiac contraction, we additively decompose S into passive and active components, such that Assuming that the myocardium is an hyperelastic and nearly incompressible material (J :¼ detðFÞ % 1), we can describe the passive stress response with adopting the orthotropic strain energy function WðEÞ from 58 : where l and j ) l are the shear and bulk modulus, respectively. E ij ¼ Ei 0 Á j 0 for i; j 2 ff; s; ng are the entries of the Green-Lagrange strain tensor E ¼ 1 2 ðC À IÞ with the right Cauchy-Green deformation tensor C ¼ F T F. We keep the anisotropic scaling factors b ij fixed using the values b ff ¼ 1, b ss ¼ 0:4, b nn ¼ 0:3, b fs ¼ 0:7, b fn ¼ 0:6, b ns ¼ 0:2 and only adjust l and a. X Valves , X Vessels , and X Peri are purely passive tissue without a preferred direction of orientation. Therefore, we use a Neo-Hookean model for these materials: withĈ ¼ J À2=3 C. All constitutive law parameters are given in Table 2.
The active second Piola-Kirchhoff stress tensor reads where n f , n s , and n n are orthotropic activation parameters in fiber, sheet, and sheet-normal orientation, respectively, and is the stretch ratio in the orientation k ¼ f 0 ; s 0 ; n 0 . Distributing the active stress across these directions can result in a more physiological deformation pattern as shown in Refs. [36] and [24]. In this study however, we choose n f ¼ 1 and n s ¼ n n ¼ 0 since an increase in stress in the sheet- Active tension development is triggered by the local activation time t A ðxÞ, which is defined as the time when the transmembrane voltage V m reaches a threshold of V thresh m ¼ À20 mV. For a description and values of the remaining parameters, we refer to Table 5. Normally, the process of excitation-contraction-coupling is mediated by an increase of intracellular calcium 32 immediately after depolarization. Atrial myocytes typically have shorter action potentials and consequently shorter calcium transients compared to ventricular myocytes due to different gene expressions and ionic channel composition. 42 To account for these differences, we adapted the atrial tension to have a shorter contraction duration such that tension development in the atria and ventricles do not overlap in time. Furthermore, the model was personalized to match MRI derived data of the subject that serves as the control in this study.
After discretizing Eq. (3) in space using linear tetrahedral elements, the balance of linear momentum equation at time n þ 1 reads where M is the mass matrix, K is the stiffness matrix, v and a are the first and second derivative of the displacement u w.r.t time, and f ext are external forces. To represent viscoelastic effects in the material model that arise in the dynamical system (3), we extend the model by a Rayleigh damping term Dv nþ1 with the damping matrix The damping coefficients were set to a 1 ¼ 500 1/s and a 2 ¼ 0:005 s. Time integration is done implicitly using a Newmark method Dt 2 a n þ bDt 2 a nþ1 ; with the parameters b ¼ 0:3, c ¼ 0:6, and Dt ¼ 1 ms. Substituting Eqs. (12), (13), (14) into (10) results in a nonlinear algebraic equation that only depends on the unknown displacements u nþ1 which is solved using Newton's method and LU-decomposition implemented in PETSc. 6

Boundary Conditions
We apply three types of mechanical boundary conditions to the surfaces shown in Fig. 1. Zero displacement Dirichlet boundary conditions are applied to the surfaces defined by C D , which comprise the distal ends of the pulmonary veins, the vena cava superior and inferior, and the outer surface of the pericardial tissue. Additionally, we use the frictionless and permanent contact problem proposed by Ref. [19] that reflects the constraints of the pericardium and the surrounding tissues on the heart: The contact is established between the surface C P on the epicardium and the surface C s , which is defined as the inner surface of the pericardium. Technically, C P and C s are connected by a spring with spring stiffness k epi ¼ 10 MPa through the gap function gðxÞ. We determine gðxÞ by projecting a point x 2 C P onto the surface C s into the direction of the current normal direction n. Depending on the distance between x and the projected point x proj , the gap function reads gðxÞ ¼ x À x proj 2 =2d for x À x proj <d; x with the transition distance d ¼ 0:1 mm and the maximal distance d M ¼ 8 mm below which contact is maintained. We made this adaptation to the original formulation to yield a smooth transition in case the surfaces C P and C s start to overlap, which results in a change of direction of n. Although we use a fixed spring stiffness k epi , the spatially varying effect of the surrounding tissue as described by Ref. [55] is incorporated into the model by assigning a smaller bulk modulus j to the pericardial tissue in atrial and basal regions of the heart. Thus, the pericardium becomes more compliant to volumetric changes in these regions and enables a radial contraction of the ventricles around the base without a loss of contact to the surrounding tissue. Finally, we use a 0D closed-loop model of the human circulatory system to determine a pressure p C for C 2 fLV,RV,LA,RAg in each of the four chambers. This pressure is applied uniformly on the endocardium as described by the Neumann boundary condition The details of this model were recently described in Ref. [22]. Additionally, the model was extended by a valve model based on the description of the instantaneous transvalvular pressure gradient according to Ref. [21] and a description of valve dynamics according to Ref. [41]. For an updated model description, we refer the interested reader to Ref. ([25], Chap. 3.4). Parameters for the circulatory system are given in Table 6. All four heart chambers of the 3D model are coupled to the 0D model using the condition with a tolerance e ¼ 10 À7 mL.

Modeling of Ablation Scars
We performed a semi-automatic ablation procedure in the left atrium of our four chamber heart model using the rule-based approach described by Refs. [34,38]. Five commonly used ablation lesions were placed in the domain X EP : pulmonary vein isolation (PVI), a mitral isthmus line (MIL), an anterior line (AL), a roof line (RL), and a posterior box lesion (BL) which includes the RL. All lesions were applied transmurally with an average width of 5 mm and mapped to the mechanical domain X M (Fig. 2). Since X M and X EP are nested meshes with different resolutions, the mapping error is minimal. In fact, the ablation lesions cover the same fraction of the atrial myocardium in both meshes. For the simulations, we consider the isolated effect of each single ablation lesion as well as all medically feasible combinations. The electrophysiological tissue conductivity in the affected regions was effectively set to zero to reflect perfect ablations. Therefore, the threshold voltage V thresh m is never reached and no active tension is developed. Additionally, we model ablated tissue mechanically as isotropic regions with increased tissue stiffness. Consequently, we adapt the constitutive model in Eq. (6) by setting the anisotropic scaling parameters The increased stiffness of the scar tissue is achieved by choosing a five times higher and l two times higher than in the bulk tissue. 43 This approximates the expected increase in tissue stiffness due to a build-up of collagen in the myocardium as observed in animal studies. 31 The sensitivity of the results regarding the scar stiffness is evaluated in an additional experiment.

Simulation Setup and Initialization
In total, we conduct ten simulations including nine combinations of the five ablation lines shown in Fig. 2 and a control case without ablation.
First, the cellular models are paced at a cycle length of 1.2 s for a total of 1000 cycles. The final values of the state variables fV m ; w; cg are used as initial values on the vertices of X EP . Second, the reference configuration X M was unloaded using a backward displacement method. 39 We used typical diastatic pressure values in the four chambers: Subsequently, the unloaded configuration was inflated with the same pressure values to pre-stress the tissue. Even though the ablation lines have different material parameters than healthy tissue, we used the same unloaded configuration for all ten simulation setups. As long as the global stiffness of the material is similar, this assumption should not affect the overall results in this study. 33 A comparison of the model outputs is only relevant if the system has reached a stable limit cycle. In the four chamber model, this means that the blood volume distribution does not change anymore during subsequent cardiac cycles. Estimating the time it takes for the model to reach this state is difficult, hence we implement an automated stopping criterion based on the stroke volume difference SV diff between the LV and the right ventricle (RV) by solving an additional differential equation in the 0D circulation model: where Q SysArt and Q PulArt is the flow through the aorta and the pulmonary artery, respectively. At integer multiples of the cycle length, SV j j diff is compared to a pre-defined threshold (1 mL) and reset to zero if the criterion is not fulfilled yet.

RESULTS
We investigated the influence of standard clinical ablation patterns on the pumping efficiency of the human heart by simulating a total of ten cases: first, we calibrated the model parameters to match the cine MRI data of the volunteer to simulate a healthy control case. Next, we introduced different combinations of the basic ablation patterns shown in Fig. 2 into the left atrium in nine additional simulations. All ten simulations were analyzed with regards to the activation sequence, blood volume, and correlation between the amount of inactive and ablated tissue with ejection fraction (EF) of the left atrium. Furthermore, we investigated the influence of the ablation lesions on the ventricular function with a focus on AVPD and ventricular stroke volume (SV).  Figure 3 shows the atrial activation time maps for all ten simulations during sinus rhythm. In the healthy control case (Fig. 3a), the activation starts at the sinus node and propagates from the right atrium (RA) to the left atrium (LA) via four pathways: the Bachmann bundle on the anterior side, a middle and upper posterior interatrial connection, and via the coronary sinus. The depolarization in the left atrium starts after 29 ms at the entry site of the Bachmann bundle. The last activation of the atria occurs after 97 ms at the posterior side of the mitral valve annulus. The ablation lesions that have the biggest influence on atrial acti-vation time are PVI + RL and PVI + AL. As a consequence of the RL lesion, the activation of the left atrial roof can only occur from the posterior side resulting in an activation delay of around 30 ms in this area (% 100 ms in PVI + AL + RL compared to % 70 ms in Control). The AL lesion blocks an activation of the left atrial appendage (LAA) through the Bachmann bundle on the anterior side, thus delaying the depolarization of the LAA significantly. Therefore, total atrial activation takes 147 ms in cases with an AL lesion. Compared to PVI + AL, PVI + MIL does not alter the activation sequence significantly, the MIL  lesion coincides with the area of latest activation in the LA in our setup. Circumferential PVI alone does not notably alter the activation sequence of the LA. Although the activation sequence is not drastically changed by PVI + BL, a large amount of tissue stays inactive due to the isolation of the LA posterior wall by this specific scar. Figure 4 shows the end-diastolic pressure-volume relationship (EDPVR) of all four chambers resulting from the inflation of the pressure-free state during the initialization process described in ''Simulation Setup and Initialization'' section. Starting from an identical pressure-free state, we observe that the volume of the LA in the pre-stressed state is reduced by up to 5.4 mL (À7:8%) if stiffer scar tissue is present compared to the healthy case. However, the reduction in volume does not correlate with the amount of ablated tissue in the LA (e.g., PVI + AL has a lower volume and less ablated tissue than PVI + BL). The volume of the other chambers remains almost identical with differences of less than AE0:1%. Figure 5 shows the amount of blood in the LA until peak LA systole. Peak systole was determined as the time when the maximal output of blood volume is reached. We can observe for all cases that it takes a slightly different amount of time to reach this point.
The change in atrial EF due to the different ablation patterns is summarized in Table 3 alongside the amount of ablated and subsequently inactive tissue. A linear regression model shows that the amount of inactive tissue is highly correlated (R 2 ¼ 0:95) with the change in ejection fraction (DEF) compared to the healthy case (Fig. 6). Ablated tissue and DEF show a much weaker correlation (R 2 ¼ 0:70).
The strain maps of the Control and PVI + AL + BL cases are shown in Fig. 7. During atrial systole, two differences can be observed in the PVI + AL + BL case compared to the Control case. First, strain values at the position of the ablation lesions suggest that the fibers remain elongated during atrial contraction (Fig. 7, dashed green circle). However, they shorten a bit compared to the diastolic state. Second, the isolated posterior roof of the left atrium gets stretched past the strain values during diastole since there is no contraction of the fibers in that area (Fig. 7, green circle). During ventricular systole, the differences in the strain maps between the Control and PVI + AL + BL are very small. Only in the scars themselves a 50% reduced fiber strain can be observed (Fig. 7, green circles). The increased stiffness of the scars has a more pronounced effect on the mean stress S ff ¼ f 0 Á Sf 0 of the tissue. In general, mean diastolic stress is increased two to three times in the posterior roof and scar surroundings (3.7 vs. 1.2 kPa). Due to the missing contraction of the posterior wall in the PVI + AL + BL case, mean stress is generally lower during atrial systole. During ventricular systole however, mean stress is increased twofold in scar tissue (42 vs. 18 kPa) while stress in the posterior wall is nearly equal between the two cases.  BIOMEDICAL ENGINEERING SOCIETY Figure 8 shows the time course of pressure and volume in all four chambers of the heart. Besides the decrease in systolic function of the LA, we observe that passive filling of the LA during ventricular contraction is restricted (7 to 10% lower compared to control) in the presence of ablation lesions. At the same time, we observe an increase in LA pressure of up to 12% compared to the control case. Furthermore, there is an additional peak in LA pressure at the beginning of ventricular systole (0.28 s) which is more pronounced in simulations that include the ablation lesions AL and BL, suggesting the changed activation sequence or the increase in inactive tissue as a source. In the case of PVI + AL + BL, the pressure is 20% higher compared to the control case. Except for small deviations in volume (within 1 mL of each other) at the beginning of the heart beat, the RA is unaffected by the ablation lesions in the LA.
The contribution of the atrial contraction to ventricular filling reduces in the presence of scars in the LA. This can be observed in the reduced end-diastolic volume (EDV) of the LV. However, this effect is not very pronounced and only reduces the volume contribution by the atria by up to 1.4% as shown in Fig. 9 (right graph). End-systolic volume (ESV) of the LV lie within 1% of each other while peak systolic pressure is reduced by 6 mmHg. Furthermore, the increased stiffness of the scars in the LA caused a reduction of the AVPD of the mitral valve (Fig. 9, middle graph) from + 3.90 mm in the control case to + 2.76 mm in PVI + AL + BL during atrial contraction. During ventricular contraction, AVPD is reduced from À 12.33 mm in the control case to À 11.87 mm in PVI + MIL + RL and PVI + RL. Similar behavior for AVPD and blood volume can be observed in the RV (Figs. 8 and 9, left graph).

Model Sensitivity Towards Scar Stiffness
Since there is not a lot of definitive data on how the stiffness of atrial scar tissue changes compared to healthy myocardium, a sensitivity analysis was conducted to investigate how scar stiffness relates to pressure and volume changes in the heart. Stiffness is mainly controlled by the parameters l and a in the constitutive model (6). Therefore, different combinations of these parameters were chosen for the scars in the PVI + AL + BL case and the outcome on pres-  sure and volume was evaluated based on the relative difference x diff to a PVI + AL + BL simulation with 2 Á l and 5 Á a: The results for the left atrium are shown in Fig. 10. The biggest difference was observed when scar tissue is modeled using the same parameters as for healthy myocardium (1 Á l and 1 Á a, anisotropic). Both, pressure and volume differences, are largest during ventricular contraction and relaxation. Pressure changed by up to 9% and volume by up to 25%. In general, an increase in stiffness is related to a decrease in volume and an increase in pressure (vice versa with decreased stiffness). Changes in the other chambers of the heart are in the range of 1 to 3%.

Reduction of Conduction Velocity
Atrial myocardium that was exposed to long lasting AF can appear to have a reduced conduction velocity due to the presence of fibrosis as a result of electrical and structural remodeling. Therefore, six additional simulations were performed to investigate how a reduced conduction velocity in the atria affects the rest of the heart. For the Control, PVI + AL, and PVI + RL cases the conductivity in both atria was reduced by 25 and 50%, respectively, to cover the diverse spectrum of atrial remodeling. Figure 11 shows the activation time maps for the three cases in both scenarios. As expected, reduced conductivities result in a reduced conduction velocity and subsequently in a prolongation of atrial activation. In the case of 50% reduction in conductivity, this results in a 19 ms delay of left atrial emptying in the three scenarios. The change in ejected volume due to this delay is  insignificant. The delay in contraction does, however, affect blood pressure in the left atrium. Especially in the case of an AL lesion, since the left atrial appendage takes up to 200 ms to activate which in turn results in a temporal overlap of atrial and ventricular contraction.

DISCUSSION
We presented an electromechanically coupled fourchamber heart model including a closed-loop circulatory system that can be used to predict acute effects of atrial ablation therapy on the heart's performance. We used patient specific measurements of AVPD and LV volume from MRI data to parameterize a healthy control case and performed a total of nine virtual ablations in the LA based on combinations of five standardized ablation lesions. The impact of these scars on deformation and cardiovascular performance was evaluated using common biomarkers such as local activation time, EF, and AVPD.
For the control case, we parameterized the active tension model for the atria and the ventricles such that the simulation results qualitatively matched the LV volume and AVPD evaluated from the volunteer's cine MRI data (see Fig. 9). Due to missing data from our healthy volunteer, the majority of parameters for the electrophysiological model and the circulatory system have been adapted from literature values of healthy subjects. Hence, most clinical biomarkers can be reproduced faithfully. With a root-mean-square error of RMSE LV Vol ¼ 0:02, RMSE LV AVPD ¼ 0:38 mm, RMSE RV AVPD ¼ 0:99 mm for the normalized LV volume, mitral valve displacement, and tricuspid valve displacement, respectively, the simulated data match the experimental data well. Therefore, we are confident that our model can reproduce physiological deformation patterns.  The ablation lesions in this work were modeled as perfectly isolating tissue. Since there is a scarcity of data on elastomechanical properties of ablation scars in the atria, we increased the mechanical stiffness induced by an accumulation of collagen in the myocardium as it was observed in animal studies. 31 Although these data are based on ventricular remodeling after myocardial infarction, we assume that it is a valid assumption, which is supported by the findings of Ref. [56]. They suggest that RFA scars behave elastomechanically like tissue in zones of chronic myocardial infarction. Additionally, ablated tissue was modeled as isotropic with no preferred direction, since Ref. [11] showed that non-reinforced infarcts formed scars with no significant alignment of collagen fibers. This should also be valid for ablations using electroporation or pulse field ablation, since the electrically isolating effect is the same as in RFA. Mechanical changes of the tissue as a result of electroporation are less clear. However, it most likely results in a stiffening of the tissue as well. The increased stiffness resulted in a change of the left atrial EDPVR that did not only correlate with the amount of ablated tissue as one might expect (compare Fig. 4 and Table 3) but rather with a combination of the extent and position of the scars which is in agreement with Ref. [48].
With the model presented in this manuscript being more comprehensive than the one published by Refs. [29] and [48], our simulations confirm their results for the behavior during atrial systole. In particular, our findings confirm a linear correlation between the change in left atrial EF and the amount of inactive tissue (see Fig. 6). This concurs with in vivo observations by Ref. [56] in 40 patients who showed a reduced left atrial EF after RFA procedures (18 AE 11%). Furthermore, a significant change in the activation sequence of the LA due to the ablation lesions AL and RL was observed. Especially the delay in activation of the LAA through the AL lesion can be problematic since a large delay in LAA activation might lead to an overlap with LV contraction. This was observed in our simulations in the form of a sudden increase in LA pressure at the beginning of LV contraction. Such asynchronous contraction patterns are known to cause electrical and mechanical remodeling of the LA. 53,54 The reduction in stroke volume of the LA results in a lower EDV of the LV. However, this only accounts for maximally 1.8% of LV stroke volume, which would likely not make a significant difference in well being for the patient. The volunteer in this study was a young and healthy male and thus the atrial kick only contributed 11-15% to the EDV of the LV. This contribution can be up to 40% especially in older people. 1 Hence, a reduced atrial contraction would likely have a much higher impact on LV performance in these patients.
The increased stiffness of the scar tissue did have a measurable effect during ventricular systole. We observed a reduced AVPD during both atrial and ventricular contraction. The restricted movement of the valve plane resulted in reduced maximal volumes and increased maximal pressures of the LA during LV systole. In the LV, the reduced stroke volume and valve plane movement result in a lower peak pressure during systole. Overall, the right side of the heart is seemingly unaffected by the ablation scars in the LA. Nevertheless, there are slight changes in volume during mid-diastole in the RV due to the adaptations of the closed-loop circulatory system to the changes in preand afterload.

Limitations and Perspectives
The focus of this study was to investigate the impact of ablation on the cardiovascular performance of the heart using a computational framework that accounts for the majority of known physiological mechanisms. In order to do so, we had to make assumptions on some aspects of the model, which may affect the results: (1) Long-lasting AF can lead to electrical, contractile, and structural remodeling of atrial myocytes. 2 In this study, we assumed that all remodeling processes are reversed after a sufficient amount of time has passed. Therefore, the results presented here are more relevant for long-term performance changes due to ablation. (2) We did not account for AF-related atrial fibrosis, which can lead to more complex propagation patterns. For this purpose, a biomechanical model of fibrosis needs to be developed first. (3) Even though it is a valid approach to use linear tetrahedral elements to solve the monodomain equation (1), it is known that low order elements especially in thin walls are prone to different locking phenomena in solid mechanics problems. 7 Nevertheless, a recent study has shown that errors in simulated circulatory system markers are small enough when compared to observational uncertainties in clinical measurements. 4 In terms of motion, strains and stresses, Ref. [4] observed minor quantitative differences but essentially the same qualitative behavior. Hence, we decided to accept small errors due to locking phenomena in exchange for faster computation to make the amount of simulations required for this study feasible. (4) The simulation results in this study are limited to a single anatomy of a healthy human heart. Furthermore, the atria in our model have a larger wall thickness than what is typically observed in healthy individuals (on average 4.5 mm compared to 2:6 AE 0:7 mm 59 ). However, the results with respect to atrial mechanics agree with independent observations made with different anatomical models in Refs. [29] and [48]. We would expect a similar impact on the ventricles and cardiovascular system despite of anatomical differences. Therefore, we are confident that our results still hold true for most patients. Nevertheless, due to the limited amount of anatomies that were investigated in a similar setting, we cannot completely rule out what kind of impact anatomical differences between individual patients would have.

CONCLUSION
This work provides an insight on how standard ablation strategies affect cardiovascular performance in a four-chamber heart model. We confirmed the results by Refs. [29] and [48] for atrial systole on a different anatomical representation of the heart and extend the study based on preliminary work of Ref.
[10] by a physiologically valid closed-loop circulatory system model and an electromechanically coupled simulation framework for the whole heart. With these additions, it is possible to show that ablation scars in the LA lower the amount of blood pumped into the LV depending on the extent of inactivated tissue. Furthermore, we can observe an increase in pressure in the LA during ventricular contraction and a reduction in AVPD due to the increased stiffness in scar tissue. These results support the hypothesis that extensive ablation not only negatively impacts cardiovascular performance but might also lead to further remodeling of the LA due to increased pressure and changes in the activation sequence.

ACKNOWLEDGMENTS
The authors would like to thank Olaf Do¨ssel for valuable discussions.

DATA AVAILABILITY
The geometry and associated data used in this manuscript is publicly available [23].

CODE AVAILABILITY
Fiber generation algorithms and other useful tools can be found at https://github.com/KIT-IBT.

CONFLICT OF INTEREST
Tobias Gerach, Steffen Schuler, Andreas Wachter, and Axel Loewe declare that they have no conflict of interest.

OPEN ACCESS
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://crea tivecommons.org/licenses/by/4.0/.

APPENDIX A: MODEL PARAMETERS
Here, we provide a complete list of model parameters that were used for all simulations in this study. In particular, Table 4 contains the parameters for the electrophysiological model, Table 5 refers to the  parameters of the active stress model, and Tables 6 and  7 contain the parameters for the circulatory system    Refer to Ref. [22] for the system of equations.