Pulmonary vein flow split effects in patient-specific simulations of left atrial flow

Disruptions to left atrial (LA) blood flow, such as those caused by atrial fibrillation (AF), can lead to thrombosis in the left atrial appendage (LAA) and an increased risk of systemic embolism. LA hemodynamics are influenced by various factors, including LA anatomy and function, and pulmonary vein (PV) inflow conditions. In particular, the PV flow split can vary significantly among and within patients depending on multiple factors. In this study, we investigated how changes in PV flow split affect LA flow transport, focusing for the first time on blood stasis in the LAA, using a high-fidelity patient-specific computational fluid dynamics (CFD) model. We use an Immersed Boundary Method, simulating the flow in a fixed, uniform Cartesian mesh and imposing the movement of the LA walls with a moving Lagrangian mesh generated from 4D Computerized Tomography images. We analyzed LA anatomies from eight patients with varying atrial function, including three with AF and either a LAA thrombus or a history of Transient Ischemic Attacks (TIAs). Using four different flow splits (60/40% and 55/45% through right and left PVs, even flow rate, and same velocity through each PV), we found that flow patterns are sensitive to PV flow split variations, particularly in planes parallel to the mitral valve. Changes in PV flow split also had a significant impact on blood stasis and could contribute to increased risk for thrombosis inside the LAA, particularly in patients with AF and previous LAA thrombus or a history of TIAs. Our study highlights the importance of considering patient-specific PV flow split variations when assessing LA hemodynamics and identifying patients at increased risk for thrombosis and stroke. This knowledge is relevant to planning clinical procedures such as AF ablation or the implementation of LAA occluders.


Introduction
Cardiovascular diseases are the leading cause of mortality world-wide.Atrial fibrillation (AF) affects around 35 million people world-wide [1], and is the most common arrhythmia.It is estimated that 20%-25% of ischemic strokes, affecting over 18 million people each year, are caused by thrombi generated in the left atrium (LA) in patients with AF.Furthermore, up to an additional 30% of all ischemic strokes are suspected to be atriogenic in patients with subclinical AF or normal cardiac rhythm [2].Inside the LA, most of these thrombi are formed in the left atrial appendage (LAA), a cavity-shaped protuberance whose morphology and hemodynamics vary significantly among patients [3].However, current medical procedures to estimate ischemic stroke risk are based on demographic and clinical factors and do not consider patient-specific information about LA hemodynamics, a crucial causal thrombosis factor [4,5].
Computational fluid dynamics (CFD) analysis based on patient-specific medical images is a powerful tool for investigating LA thrombosis.In the past decade, advances in CFD algorithms, numerical algorithms and computational resources have supported increasingly complex cardiac flow simulation studies.These works cover four-chamber simulations of the whole heart [6], two-chamber simulations of the left heart including the atrium and ventricle [7][8][9][10], and single-chamber simulations of the LA or LV.These single-chamber simulations have focused mostly on LV hemodynamics [11][12][13][14][15][16], although there is a growing body of literature exploring different aspects of LA hemodynamics.One of the most studied aspects is the effect of LA wall motion on the flow and LAA stasis [17][18][19][20][21], motivated by the fact that AF causes weak irregular LA wall motion and LAA stasis is associated with increased thrombosis risk [1,[22][23][24].Motivated by the clinical association between LAA morphology and stroke risk in AF patients [25][26][27], several groups have used CFD to investigate how LAA shape influences surrogate hemodynamic metrics of thrombosis risk [28][29][30][31].In parallel, recent work [10,[32][33][34][35][36][37] is introducing advanced CFD, fluid-structure interaction, and multi-physics models, and integrating them with multi-modality imaging to realistically simulate the electrophysiological and biomechanical mechanisms underlying LAA stasis and thrombosis.
The aforementioned simulation studies have provided useful guidance for setting up models and boundary conditions that reproduce each patient's hemodynamics representative of a specific time or physiological state.However, both errors in model approximation and the natural variation of physiological state during a patient's daily life (e.g., rest vs. exercise) constitute a source of uncertainty.These errors must be understood, evaluated, and accounted for before model predictions of LAA thrombosis based on hemodynamic computations can be used to support clinical decisions.The PV flow split, i.e., how the flowrate of blood entering the LA from the lungs is divided between the left and right PVs, is a particularly challenging parameter affecting the inflow boundary conditions.Its quantification via transthoracic Doppler echocardiography (TTE) is complicated by the difficulty of imaging all (usually four) pulmonary veins [38] and, while it can be accessed by other modalities, such as 4D flow MRI, these techniques are more cumbersome and time-consuming than TTE.Consequently, the PV flow split is typically assumed to be even for simplicity, even if the difference in size between the left and right lungs makes this split be closer to 47%/53% (left/right) [39], which can vary widely from patient to patient (a range of 43-48/57%-52% for N=206 patients was reported in [39]).Moreover, the PV flow split can experience significant intra-subject variability due to interactions between gravity and body position, with reported differences between lying on the left and right sides of the body being 59/41% and 37/63%, respectively [40].These significant inter-and intra-patient variabilities in PV flow split warrant investigating its effect on LA hemodynamics and LAA blood stasis.
There are few works investigating the influence of PV inflow in CFD simulations of atrial blood flow, and most of them do not directly address the left/right flow split.A systematic parametric study was carried out in [41], including a thorough sweep over a wide range of different PV flow splits and validation vs. in vivo 4D flow MRI measurements.Their results demonstrated that atrial flow can be sensitive to the PV flow split, especially for large departures from normal values.However, they only considered N=3 patients and did not address LAA flow or residence time.Previously, the effect of PV inflow was investigated in [28] as part of a broader sensitivity analysis of the effect of LA and LAA morphological features on atrial flow dynamics.By varying PV cross sectional area, these authors effectively evaluated the effect of total PV flow rate on LA hemodynamics, showing significant effects.However, they kept an even 50/50 split, so the study did not provide information about the influence of the flow split.Subsequently, patient-specific simulations were performed for a large cohort of N=52 patients in [24], concluding that the position and orientation of the right PVs can significantly influence LA hemodynamics.Nevertheless, neither of these previous works [24,28] addressed the effect of the PV flow split.More recently, some works [42,43] have studied LA hemodynamics before and after surgical removal of left upper lung lobe, concluding that the resulting left superior pulmonary vein recession disturbs blood flow all throughout the LA.In summary, there is a paucity of data about how the PV flow split affects LAA blood flow and stasis, and about its differential effects on normal and impaired/pro-thrombotic atria.
Motivated by this paucity of data, this study's goal was to evaluate how the inter-and intra-subject variability in PV flow split affects LA blood flow and stasis.To that end, we built patient-specific models based on 4D computed tomography (CT) images, including total flow rate through the PVs.For each model, we considered four PV flow split variations encompassing the even split commonly used in CFD analysis and more physiological, asymmetric splits with higher flow through the right PVs.Particular attention was given to PV flow split induced variability in LAA blood stasis.
The paper is organized as follows.Section 2 presents the methodology, including the extraction of LA geometries from CT images and the details of the CFD analysis.Section 3 describes the main results of the study, showing the effect of varying PV flow split on LA and LAA flow focusing on LAA stasis.Section 4 critically discusses these results and the study limitations, and outlines future work.Finally, the manuscript concludes in Section 5, which concisely summarizes the main results and contribution of the study.

Imaging and generation of patient-specific anatomical models
We studied N = 8 cases through 3D, time-resolved (a.k.a.4D), patient-specific segmentations of human left atria obtained from computed tomography (CT) images.Six of these segmentations are the same as in [23,36].Two additional subjects were imaged at Hospital General Universitario Gregorio Marañón (HGUGM), Madrid, Spain, following standard clinical protocols using a Philips Brilliance 64 scanner.The doses and injection rates of contrast were chosen based on the patient's weight using standard clinical protocols.The images were reconstructed using CT scanner manufacturers' algorithms, yielding DICOM files with pixel dimensions between 0.375 mm and 0.42 mm in the x-y plane and 0.45 mm in the z-direction.For the other six cases, pixel dimension was 0.32 mm to 0.48 mm in the x-y plane and 0.5 mm to 1 mm in the z-direction.Time-resolved images were obtained at regularly spaced instants across the cardiac cycle, ranging between 5% and 10% of the R-R interval.
The computational LA meshes were generated in four steps using ITK-SNAP [44] and custom-written scripts in MATLAB.The first step comprised segmenting the 3D LA anatomy from CT images and identifying the PV inlets, mitral annulus, and LAA.For each 3D segmentation, a triangular surface mesh was created and then resampled to match the computational fluid dynamics (CFD) solver's resolution [45].The resulting triangular meshes were registered across the cardiac cycle and their positions, yielding a coherent triangle vertex and centroid cloud [46,47].Finally, the positions of these points were expressed as Fourier temporal series to provide interpolated boundary conditions to the CFD solver at time points not contained in the 4D CT sequence.Table 1 summarizes the main anatomical and functional features of the study subjects' atria.Fig. 1 shows the N = 8 patient specific anatomies at the beginning of the R-R interval.Additional details on image acquisition and reconstruction and mesh generation can be found elsewhere [23].
These 4D CT images were also used to measure the LA and LV volumes along the cardiac cycle, and to compute patient-specific flow rates through the mitral valve and pulmonary veins (as described in Section 2.3).

Computational fluid dynamics analysis
The CFD simulations were performed with the in-house code TUCAN, as described in [23].
Here, we summarize the main features of these simulations, emphasizing the aspects that are particularly relevant to the present study of PV flow split variability.
TUCAN solves the Navier-Stokes equations for the incompressible where v and p are the velocity and pressure fields, and ρ and v are the density and the kinematic viscosity of the fluid.The equations are solved using second-order centered finite differences and a three-stage, semi-implicit, low-storage Runge-Kutta scheme.In the present study, we considered a Newtonian fluid of constant kinematic viscosity v = 0.04 cm 2 /s [23].
The flow in the LA was driven by each subject's left heart wall motion, obtained from patient-specific 4D CT images as described above.This motion was imposed using the immersed boundary method (IBM) [48,49], where the flow is solved in a fixed, Cartesian grid and the position and motion of surfaces are imposed by a localized volumetric force.Each simulation was run in a cubic uniform grid of 256 3 points with grid spacing Δx = 0.051 cm in all directions and with a constant time step of t = 5 ⋅ 10 −5 s, using free-slip boundary conditions at the cube boundaries.This time step is chosen to ensure that the Courant-Friedrich-Lewis number is always CF L < 0.3.The numerical error of the time integrator is proportional to a power of the CF L number [50], so it must be small to ensure numerical stability and accuracy in the results.The simulations were initialized from a static flow field and, before running in the nominal conditions, some cardiac cycles were run with a lower resolution of x = 0.09 cm in all directions and a constant time step of t = 10 −4 s.
In total, we run 10 heartbeats at low resolution to accelerate convergence to periodic flow, followed by 5 more heartbeats at the nominal spatial and time resolution.
Similar to [23], the mitral valve was modeled as a flat impermeable surface during LV systole and as an open boundary during LV diastole.

Pulmonary vein boundary conditions
To enforce flow rates through the PVs, a buffer region was defined upstream of each PV inlet plane, and the IBM was applied in this region to bring the flow velocity to the desired vector value.The target inlet velocity at each PV inlet plane was defined as where A i and n are the instantaneous cross-sectional area and the vector normal to each inlet's plane, respectively, and Q i i = 1, …, 4 are the instantaneous flow rates through each PV.
The total flow entering the LA through the four pulmonary veins, Q P V t , was determined from mass conservation (Fig. 2) [23]: where V LA t is the time-dependent LA volume obtained from each patient's 4D CT scan and Q MV t is the flow rate exiting the LA through the mitral valve, calculated from the time rate of change of the LV volume assuming no mitral regurgitation.This approach allows for setting patient-specific PV inflow conditions without the need for measuring flow velocities and should account for AF-mediated changes in PV flow waveforms.
To analyze how the PV flow split affects atrial hemodynamics and blood stasis, we considered four flow split scenarios: a. R60: In this scenario, 60% (40%) of the flow enters the LA through the right (left) PVs, where A RP V and A LP V are the sum of the areas of the right and left PVs respectively.Usually, both the left and right PVs have two inlets into the LA body (Fig. 2).The same velocity value is imposed at each side's inlet.

b. R55:
This scenario is similar to the R60 but with a 55/45%R/L flow split.

c. R50*:
The flow rate is split evenly among the four inlets formed by the right and left PVs, so each inlet's velocity depends on its sectional area.

d.
eqV: Same velocity through each PV, Where A P V is the sum of the areas of all PVs.In this scenario, the flow rate through each PV is not prescribed as it depends on the corresponding sectional area.The eqV flow splits for the 8 cases are shown in Table 1.

Flow analysis
This section describes the post-processing of CFD-computed blood velocity fields.First, we describe a framework to divide the LA in different sub-domains to study the regional impact of PV flow split variations.Then, we define correlation factors that quantify the similarity of LA hemodynamics for different PV flow splits.Finally, we describe the computation of residence time to quantify blood stasis.
To analyze the flow inside the different regions of the LA (see [51], for reference), we divided each anatomical model into three domains (Fig. 3A): the posterior region receiving the PVs (denoted LAPV), the appendage (denoted LAA), and the remaining part of the LA body that includes the septal portion and the vestibule (denoted LAb).While the LAA has a relatively well-defined opening, the other parts do not have clear anatomical demarcations.Thus, we devised an objective procedure to separate the LAPV and LAb domains.We defined an intrinsic coordinate ξ that defines planes parallel to the mean MV crossing the LA body (Fig. 3B).Instantaneous MV planes are computed by least-squares fitting the points that define the MV in each subject's anatomical model.The coordinate ξ is normalized to range from ξ = 0 at the centroid of the PV inlets to ξ = 1 at the centroid of the MV at each instant of time.This coordinate allows for establishing an objective criterion to dynamically demarcate the LAPV and LAb (Fig. 3B).The separation from the LA body and the LAA is defined by a plane that crosses the LAA ostium as part of the CT image segmentation.
To quantify the dependence of LA hemodynamics on the PV flow split, we define correlation factors [52] for the PV inflows and the velocity fields in the LA.The correlation factor of PV flow rates is defined as where is the flowrate waveform corresponding to the kth PV with the flow split condition indicated by i = R60, R55, R50 * , eqV .Likewise, the correlation factor of the 4D velocity field is defined as where u i x is the velocity field at point x of the LA at the flow split condition indicated by i.
To measure blood stasis, we define the residence time T R as the time spent by blood particles inside the LA chamber.This quantity was calculated together with the velocity field by solving the forced transport equation [53], where T R is set to zero at the start of the simulation and at the PV inlets.The spatial discretization of this equation employed a thirdorder WENO scheme [23] to ensure accuracy while avoiding spurious fluctuations near steep gradients.We note that, up to numerical errors, this procedure is equivalent to tracking Lagrangian flow tracers and measuring the time elapsed along each tracer's trajectories [53].

Study population
The LA patient-specific anatomical models and functional characteristics of the study group are indicated in Fig. 1 and Table 1.The cohort's median age was 65 (range 48-92), and four (50%) were female.Subjects 4 and 5 were imaged at HGUGM, whereas the rest of the patients were the same reported in [23].Subjects 1 -5 were imaged in sinus rhythm and had no LAA thrombus.These subjects had normal LA function, as can be deduced from chamber volumes and PV flows derived from CT images.Subjects 6-8 had permanent atrial fibrillation.Subjects 6 and 7 were imaged in atrial fibrillation and had an LAA thrombus, which was digitally removed before the segmentation of the LA geometry.Subject 8 was imaged in atrial fibrillation and did not have an LAA thrombus; however, this patient had a previous history of transient ischemic attacks (TIAs).Subjects 6 -8 had an enlarged LA with impaired global function.While subject 6 had decreased reservoir function but relatively normal booster function, subjects 7 and 8 had impaired reservoir and booster function.In particular, subject 7 had severely impaired function, operating primarily as a conduit between the PVs and the LV.With the exception of patient 8, all subjects had normal left ventricular function as assessed by this chamber's ejection fraction.Based on this information, we categorized subjects 1 -5 as LAA-thrombus/TIA negative (LAAT/ TIA-neg), and subjects 6 -8 as LAAT/TIA-pos.

Effect of PV flow split on left atrial flow patterns
Fig. 4 shows 3D velocity vector maps representative of normal atrial function (subject 2) for the different flow split scenarios considered in this study.Flow maps in a MV-parallel plane of the LA body and the LAA are represented at LA diastole and LV early filling, together with total PV flow rate waveforms.During atrial diastole, the mitral valve is closed while the left atrium and its appendage are filled with blood entering through the PVs.In this phase of the cardiac cycle, the flow patterns in the MV-parallel plane exhibit appreciable differences between flow splits, especially in the LA body (Fig. 4A).These differences are noteworthy considering that the PV inflow waveforms do not experience a striking dependence on the PV flow split (Fig. 4C).In this normal subject, the LAA flow patterns only exhibit modest differences when the PV flow split is varied.During LV relaxation, some flow is still entering through the PVs, however, the MV is now open and LV suction accelerates the blood towards the mitral annulus.In this phase, the flow organization in both the LA body and the LAA seems to be less sensitive to the PV flow split, probably under the influence of the strong atrial emptying jet (Fig. 4B).The 3D velocity vector maps indicative of impaired, pro-thrombotic atria (subject 6) are represented in Fig. 5. Similar to the normal case, the PV flow split has a noticeable effect on the flow patterns contained in the MV-parallel plane.In addition, the LAA flow patterns appear to be more sensitive to the PV flow split than in the normal case.Below, we evaluate these qualitative observations more quantitatively.
Figs. 6 and 7 show the blood kinetic energy KE and residence time T R in a plane section of the LA body and the LAA for the same two subjects of Figs. 4 and 5.During atrial diastole, the KE reaches moderate values at the PV inflow jets and is significantly lower elsewhere.During LV relaxation, the PV region still presents moderate KE values while its peak values are associated with the emptying jet.This overall pattern is conserved for the four flow-split scenarios considered; however, the KE maps show appreciable differences when comparing different flow splits and subjects.These differences are obvious in areas of high KE values corresponding to inflow and outflow jets.However, there are also important relative differences in the flow KE inside the LAA, as evidenced when representing this quantity using a ~ 10-fold lower range (Figs.6C, 7C).
After 13 cardiac cycles, blood residence time inside the left atrium showed stable distributions with small, quasi-periodic variations along the cardiac cycle.The lowest T R values appear near the PV inlets, where blood enters the chamber.The incoming blood travels towards the mitral valve, washing out the LA body.This clearing process is particularly effective in subjects with normal atrial function, keeping T R around 1 cycle inside the LA body.But even in subjects with impaired atrial function, the LA body ultimately acts as a conduit, and T R rarely exceeds 2 cycles in this region of the chamber.Consistent with previous works [23,36], the T R is significantly longer inside the LAA than in the LA body.Moreover, and also in contrast to the LA body, the LAA experiences significant variability in these two hemodynamic variables with respect to the PV flow split.
For both the normal and impaired, pro-thrombotic atrium, T R is elevated in the distal LAA regardless of PV flow split; however, the proximal T R distributions exhibit more variability, especially in the pro-thrombotic case.This variability is captured by computing the spatial average of T R , which is plotted vs. time together with KE in Figs.6-7C.These plots confirm that changing the PV flow split alters the T R inside the LAA all along the cardiac cycle and that, as noted above, these effects are more considerable in the impaired, pro-thrombotic atrium (Fig. 7) than in the normal one (Fig. 6).The next sections quantify this variability in more detail across our entire study cohort.

Flow sensitivity to changes in PV flow split
To study in more detail how PV inflow affects blood flow in different regions of the atrium, we analyzed the correlations between flows obtained with different PV flow splits in different LA regions.We first computed the correlation factor of the PV inflow flowrates along the cardiac cycle (ρ P V , defined in Eq. ( 6), summarized in Table 2), observing small departures from perfect correlation ρ P V = 1 for all cases and low subject-to-subject variability.To investigate the downstream effects of these differences in atrial inflow, we computed correlation factors for the flow velocity inside the atrial body (Eq.( 7)) in the directions parallel ρ LAb, ∥ MV and perpendicular ρ LAb, ⊥ MV to the MV plane.In this calculation, we excluded the mesh points inside the LAA and the region near the PVs, as indicated in Fig. 3.These data are summarized in Table 2 and plotted in Fig. 8A-B as functions of ρ P V .
The results suggest that both ρ LAb, ∥ MV and ρ LAb, ⊥ MV tend to decrease with ρ P V , indicating that stronger variations in PV flow split cause stronger variations downstream in the LA body.
We also found that ρ LAb, ⊥ MV has higher values and a shallower dependence on ρ P V when compared to ρ LAb, ∥ MV across all the simulations.These data indicate that variations in PV flowrates have a weaker effect on the main axial transit of blood from the PVs to the MV than on secondary, transverse flow patterns.In the discussion section, we argue that the different sensitivities to PV flow split are related to global chamber mass conservation.Of note, the thrombus/TIA-positive cases exhibit a steeper decrease of ρ LAb, ∥ MV than the normal cases, suggesting that secondary flow patterns of pro-thrombogenic atria with impaired function (subjects 6-8) are more sensitive to PV inflow flowrates, consistent with the vector velocity, KE, and T R maps presented above.The differences in the sensitivity of the hemodynamics between normal and impaired, prothrombotic atria are particularly noticeable in the LAA, as indicated by the velocity correlation factor ρ LAA vs. ρ P V inside the appendage (Fig. 8C).

Kinetic and residence time statistics in LA body and LAA
Atrial kinetic energy and residence time have been proposed as hemodynamic biomarkers of thrombosis risk [23,36].To study how PV flow split variability affects these variables, we plotted their probability box plots using data from three cardiac cycles after converging the simulations (12 s < t < 15 s).Fig. 9 displays these boxplots for the LA body, showing that the flow in the direction perpendicular to the MV plane (KE ⊥ , Fig. 9A) is more energetic than the secondary flows parallel to that plane (KE ∥ , Fig. 9B).Overall, normal atria experienced higher KE than impaired, pro-thrombotic ones.
To quantify the level of variability across PV flow splits, we computed z-scores for each patient as z = mean μ j − μ i / s i 2 + s j 2 1/2 , where μ i is the mean and s i is the standard deviation of the physical quantity for flow split i, and included these data in the figure.Consistent with the correlation factor data reported above, KE ⊥ was less sensitive to the PV flow split than KE ∥ .Specifically, the z-score medians (range) of KE ⊥ over the normal and pro-thrombotic atria were respectively 0.03(0.01-0.03)and 0.01(0.01-0.03),while the corresponding values for KE ∥ were 0.05(0.05-0.14)and 0.04(0.03-0.05).Interestingly, these z-scores suggest that the flow KE of normal and pro-thrombotic atria have similar sensitivities to the PV flow split, which contrast with the results from the correlation factors shown in Fig. 8.Our interpretation of this seeming discrepancy is that it is the orientation and not the overall magnitude of the velocity field which is more sensitive to PV flow split in the pro-thrombotic atria.It is also worth noting that many of the simulations produce the highest atrial body kinetic energies for the R50* flow split.This peculiarity, which is particularly noticeable for subject 4, can be explained by the differences in PV inlet areas for each subject.The R50* condition splits the PV flowrate evenly among the four PVs, often creating a high inflow velocity through the narrowest PV inlet.
Fig. 9C displays barplots of the residence time averaged over 12 s < t < 15 s inside the LA body.Overall, normal atria experienced lower T R than impaired, pro-thrombotic ones.A recent work [36] suggested that global chamber LA mass conservation and conduit function make T R inside the atrial body to be approximated by the simple model T R = LAV/ α LVSV where LAV and LVSV are respectively the mean LA volume and LV stroke volume, and α > 1 is a parameter that accounts for atrial washout occurring mostly during early LV filling.Fig. 9C confirms the model, which helps understand the low variability of T R with respect to flow split changes for all the cases.Also, the z-score medians (range) of T R over the normal and pro-thrombotic atria were respectively 0.02(0.01-0.03)and 0.05 (0.04-0.07), suggesting that T R in the body of pro-thrombogenic atria is more sensitive to PV flow split when compared to normal atria.
Fig. 10 shows probability box plots of kinetic energy and residence time inside the LAA using data from the same three cardiac cycles (12 s < t < 15 s), analogous to Fig. 9. Overall, the magnitude of KE in the LAA is significantly lower than the total KE in the LA body, although it is comparable to that of atrial body secondary flow patterns (i.e., KE ∥ ).In normal atria, we observe slight differences in LAA KE when the flow split is varied within each subject.These differences become more evident in impaired, pro-thrombotic atria, which, in general, also present smaller LAA KE values.Consonantly, the z-scores for LAA KE have a median (range) of 0.05(0.01-0.08) in normal atria vs. 0.1(0.08-0.14) in pro-thrombotic atria.
Residence time inside the LAA is significantly higher than inside the LA body, specially its maximum values, which often correspond with stagnant blood pools in the LAA apex (see, e.g., Figs. 6 and 7 and refs [23,36]).In normal atrial, the median (range) of the LAA T R z-scores is 0.03(0.01-0.12),and two of the five cases have z-scores > 0.1.In pro-thrombotic atria, the variability induced by PV flow split in LAA T R was more pronounced, with z-score statistics being 0.17 (0.13-0.25).
The variability in LAA residence time induced by the PV flow split could be a significant source of uncertainty when using T R to identify patients at risk of atrial thrombosis.To evaluate this possibility, we plotted the mean LAA residence time (12 s < t < 15 s) vs. PV flow split condition for this study's normal (Fig. 11A) and impaired, pro-thrombotic (Fig. 11B) atria.These plots reveal that the normal atria with highest T R values overlap with the envelope of the pro-thrombotic cases (red shaded patch in Fig. 11A), while the pro-thrombotic atrium with lowest T R overlaps with the envelope of the normal cases (blue shaded patch in Fig. 11B).These overlaps are primarily caused by the significant sensitivity to PV flow split of the residence time in the pro-thrombotic LAA.

Discussion
The past decade has witnessed significant growth in the number of computational fluid dynamics (CFD) analyses of the left atrium [17,18,23,24,29,37,54].These works shed light onto the anatomical and physiological determinants of atrial flow and offered proofof-principle validation for blood stasis mapping as a thrombosis predictor.In parallel, echocardiographic LV stasis mapping has proven predictive value for LV thrombosis and brain embolism in pre-clinical and clinical pilot studies [55][56][57].However, blood flow in the LA is highly three-dimensional and unsteady, and CFD-simulated LA hemodynamics can exhibit significant intra-patient variability depending on modeling parameters [17,21,35,54].
A better understanding of this variability is required to develop reliable CFD-derived biomarkers of atriogenic stroke risk.

Pulmonary vein inflow, atrial hemodynamics, and LAA stasis
The LA contributes to cardiac pumping via three distinct functions.When the mitral valve is closed, it receives flow kinetic energy from the lungs and transforms it into myocardial elastic energy that is released back into the flow during early LV filling (reservoir function).
It actively contracts to drive late LV filling (booster function).And, while the mitral valve is open, it forms a tube through which the pulmonary veins continuously discharge into the LV (conduit function) [58].These three functions involve flow interactions between the atrium and the pulmonary veins and can overlap during some phases of the cardiac cycle, creating complex PV inflow profiles with several forward waves during atrial and early LV diastole, and alternating forward and backward waves during atrial systole.There are usually two pairs of PV inlets in the LA, the right pair opposes the left pair creating two opposing inflow jets.Furthermore, the precise number and anatomical configuration of PV inlets can vary significantly among patients [24].
These factors give rise to intricate atrial flow patterns, first visualized in 3D approximately two decades ago [59,60].However, it has not been until recently that the influence of PV inflow in LA hemodynamics began to be investigated in detail using CFD simulations [24,28,41].In particular, 4D flow CMR measurements and CT-based CFD simulations of the whole left heart, varying the PV flow split, were performed in [41].A strength of this study was its design of experiments (DOE) approach to interrogate the three-dimensional parameter space defined by varying the PV flow split of the four pulmonary veins.Following this systematic approach, they ran 20 simulations per patient in N=3 patients, which provided unprecedented information on the PV inflow dependence of LA and LV flow.They reported a large variability for the LA body KE, although this may be caused by the fact that some of the PV flow splits considered in the DOE analysis were unphysiological.On the other hand, their 25-25-25-25 simulations, which are equivalent to our R50* condition, yielded LA KE profiles in the middle of the DOE envelope and were in fairly good agreement with the 4D flow MRI patient-specific measurements.Unfortunately, they did not compute residence time or report flow data inside the LAA, the most frequent site of intracardiac thrombosis [61].
More recently, an exhaustive investigation of the influence of atrial anatomical parameters on LA hemodynamics and LAA stasis was performed in [24].Their approach was complementary to that of [41]; they considered a large cohort of atrial fibrillation patients (N = 52,25 of them thrombotic), performing one simulation per patient using the inflow/outflow profiles for all simulations and computing a comprehensive set of hemodynamic descriptors in the LAA.Then, they interrogated the database for correlations between anatomical features and hemodynamic descriptors across different patients.Of particular interest for the present study, the orientation of the PV inlets was analyzed in [24], finding that there is significantly more variability in the right-side PVs than in the left-side PVs, and that this variability significantly affected atrial hemodynamics.Nevertheless, they did not find a clear association between right PV orientations and LAA thrombosis.Likewise, they found no direct link between LAA thrombosis and the number of PVs, emphasizing the multi-factorial nature of LAA stasis Thus, the existing literature suggests that the effects of the PV inflow profiles on atrial hemodynamics are significant and may be difficult to disentangle from anatomical factors.Furthermore, despite these recent works, there is a paucity of data on how the PV flow split influences LAA stasis, a major determinant of atrial thrombosis associated with stroke risk.This paucity creates a significant problem because the PV flow split is a crucial parameter to define inflow boundary conditions and its patient-specific measurement is challenging [38].To shed light on this question, we performed a CFD investigation on a moderate-size cohort (N=8) while also considering four different PV flow split scenarios per patient-specific anatomical model.While our study did not perform an exhaustive DOE parametric analysis [41] or considered a large cohort [24], it provides new insights that have not been previously explored in the literature.In particular, this study balances, for the first time, inter-subject and intra-subject variability by including normal and thrombotic atria, focusing on physiologically relevant departures from the 25-25-25-25 flow split.We conducted high-fidelity CFD simulations in terms of spatial and temporal resolution (see, e.g., ref [62]), the number of cycles run, and computed high-cost hemodynamic variables like blood residence time, which are tightly linked with blood stasis [63][64][65][66], and have been clinically validated as thrombosis biomarkers in other cardiac chambers [55][56][57].
Long-term analysis has not been addressed in this study, however, additionally to the effects on transport and thrombosis phenomena, we can speculate that potential variations on secondary flows could involve changes in endocardium stresses and mechanotransduction under AF conditions, as suggested in [67].
In agreement with the previous studies, instantaneous visualizations of our CFD simulations show that the PV flow split modifies the flow patterns in the LA and LAA.While this effect is complex and multifactorial, our analysis suggests a few general trends.First, the velocity along the PV-MV axis delineating the atrial conduit route (denoted as the ⊥ component) is relatively insensitive to the PV flow split.A possible explanation for this behavior is that the flow rate through the MV is dictated by global mass conservation of the LA chamber independent of the flow split, and this flow rate, in turn, dominates fluid motion along the PV-MV axis.Consistent with this hypothesis, a previous work [41] reported CFD-derived flow visualizations in the MV plane that were relatively insensitive to PV flow split and agreed well with 4D flow MRI measurements.Also consistent with this hypothesis, our simulations show that fluid motions parallel to the MV plane (denoted as the ∥ component), which are not constrained by global LA mass conservation, are more sensitive to the PV flow split.Since the flow velocities along the PV-MV axis are stronger, this effect is most apparent when separately computing the flow kinetic energy in the | and ⊥ components.Blood residence time in the atrial body, which is primarily dictated by the ⊥ component, is relatively insensitive to the PV flow split.The sensitivity of LAA hemodynamics to the PV flow split is more involved and is discussed below.

The hemodynamic sensitivity of the left atrial appendage
The LAA is a small tubular sac protruding from the atrial body in a direction roughly perpendicular to the main transit path of LA flow between the PVs and the MV.
Consequently, blood residence time inside the LAA significantly exceeds that in the atrial body [17].Given its anatomical isolation from the conduit that joins the PVs and the MV, one could intuit that LAA hemodynamics should be relatively insensitive to the PV flow split.However, our results indicate that blood flow and, in particular blood residence time, are more sensitive to the PV flow split in the LAA than in the atrial body.As we argued above, secondary flow patterns perpendicular to the PV-MV axis (i.e., the | component) are not significantly constrained by global chamber mass conservation and, thus, could vary appreciably with anatomical parameters, boundary conditions (including the PV flow split), blood rheology (e.g., non-Newtonian effects [36]), etc.Our data suggest that changes in secondary flow patterns near the LAA ostium can affect flow inside the LAA and, as a result, the LAA residence time.This idea is supported by Garcia-Isla et al.'s data showing that perturbations to the PV diameter that increase blood velocity at the ostium also decrease the endothelial cell activation potential [28].This mechanism is highly dependent on each patient's LAA anatomical features and its specific orientation with respect to the PV inlets and the PV-MV axis [24], which would explain why some of our subjects exhibit small variations in mean LAA T R with PV flow split (i.e., cases 1-3), while others exhibit large enough variations to create overlap between the thrombus/TIA group and the normal group.
Unlike the atrial body, the appendage is a closed chamber with no conduit function that relies primarily on the expansion and contraction of its walls for blood clearance.This mechanism is quantified by the LAA emptying fraction, (LAA EF, Table 1), which reached 0.53+/-0.08(mean +/-std) in this study's normal atria and 0.21 +/-0.01 in the impaired, prothrombotic atria.The significant difference in LAA EF between groups is congruent with the higher TR found in the impaired prothrombotic atria and, relevant to the present study, could also explain the higher variability of this group's LAA T R with the PV flow split.This may be because secondary swirling motions with zero net volume flux play a more important role in washing out LAAs of lower emptying fraction, and these fluid motions are more sensitive to anatomical factors and PV flow profiles [36].Consistent with this hypothesis, we recently showed that the LAA T R of normal atria is less sensitive to non-Newtonian effects than the LAA TR of impaired prothrombotic atria [36].Finally, we note that we cannot discard the alternate hypothesis that LA EF and not LAA EF drives the variability of hemodynamic variables, since the two EFs are highly correlated in our cohort (correlation coefficient = 0.90, p = 0.002).

Study limitations and future work
The cohort used in this study is small (N = 8).This is partly due to the computational cost of running high-resolution CFD simulations for 15 cycles at CF L < 0.3 and the need to run four simulations per subject to study the effect of the PV flow split.However, high temporal resolution is crucial to obtain convergence in hemodynamic metrics obtained from CFD analysis of LA hemodynamics [62].To help compensate for the low N, we selected a diverse cohort in terms of LA volume and function (see Table 1).While this strategy allowed us to cover a wide range of patient-specific conditions, it might bias our observations regarding inter-patient vs. intra-patient variability.In addition, all of the 8 patients in this study had normal pulmonary vein anatomy with 2 right and 2 left pulmonary veins.Whether the results extrapolate to anatomic variants was not evaluated.
The small sample size and retrospective nature of our study makes it difficult to ascertain if the observed overlap in LAA stasis between AF and normal patients is significant, as this overlap depends on additional factors, e.g., LAA stasis could vary for each given patient based on how atrial function and flow evolve after the onset of AF or post-ablation [68], and this study did not record such clinical variables or is statistically powered to control for them.
We consider PV flow splits in the ranges (left/right) 40-50%/50%−60% to compare differences in LA hemodynamics between the even split [7,17,23,36,69,70] and equal velocity inflow conditions [19,21,33,71,72] often used in simulation studies, and the more physiological uneven splits caused by the right lung being larger than the left lung.However, multiple factors can alter the PV flow split outside of the range we studied, including common ones like changes in body position [40] and less common ones like abnormalities in the pulmonary artery or veins [73,74], or lung resection surgeries [75].Additionally, body position can also affect superior/inferior PVs flow split, whose effect would be interesting to study in future work.
All our simulations were run at a constant heart rate of 60 min −1 corresponding to resting conditions and constant (i.e., Newtonian) viscosity.This choice is customary [21,23,29,62,70,76] and justified based on the rationale of reducing the number of independent parameters to the PV flow split.The flow split has been reported to change with heart rate, approaching 50%−50% during exercise [77], and heart rate could alter blood viscosity inside the LA via non-Newtonian effects [36].In turn, non-Newtonian effects could be interlinked with flow-split-induced variability in the LAA, where long residence times contribute to increased blood viscosity [36].
To close the manuscript we include some ideas for future work.To increase the statistical relevance of the results, it is necessary to consider a larger cohort of patients, including new cases with both normal and impaired/pro-thrombotic atria and a wider range of PV flow split.Given the anatomical variability in the general population, it would be interesting to study how the results are influenced by variations in the geometry of PVs.This study would analyze the flow from LA geometries where a longer section of the PVs is included and compare it with the results of the present study.In this sense, studying LA geometries with more or less than four PVs [78] would also be pertinent.Additionally, many other relevant problems can be addressed, such as analyzing the effect of superior-inferior PV flow split on LA and LAA flow or studying LAA stasis in patients in particular situations, e.g., after AF ablation or a lung lobectomy surgery [42].

Conclusion
The pulmonary vein flow split is an important parameter governing boundary conditions in patient-specific left atrial blood flow simulations.This work investigated, for the first time, how the PV flow split influences left atrial flow and LAA stasis by performing CFD simulations on subjects with normal and impaired pro-thrombotic atria.We found that secondary flow patterns in planes parallel to the mitral valve are sensitive to variations in PV flow split.Flow inside the LAA, particularly blood residence time, is also affected by the changes in PV flow split.This sensitivity is highest in patients with reduced LAA emptying fraction, highlighting the need for accounting for uncertainty in CFD-derived risk scores of atrial thrombosis.Finally, several ideas for future work have been discussed, including the simulation of a larger cohort of patients and a wider range of PV flow split conditions, studying the influence of PV anatomy, or analyzing the effect of varying superior-inferior PV flow split.Geometry and flow distribution in the left atrium.Q P V is the inflow through the pulmonary veins (RS Right Superior, RI Right Inferior, LS Left Superior, LI Left Inferior).V LA is the time-dependent volume of the left atrium.Q MV is the outflow through the mitral valve.LAA is the left atrial appendage.

Fig. 1 .
Fig. 1.Anatomical geometry of patient-specific left atrium subjects.3D mesh extracted from Computerized Tomography (CT) images of the LA walls and PVs inlets (left PVs in green, right PVs in red) and MV outlet surfaces (blue).The images correspond to an instant at the beginning of the R-R interval.

Fig. 3 .
Fig. 3.(A) Division of the left atrium geometry into three regions: pulmonary veins and the portion of the left atrium adjacent to their ostia (blue), main left atrial body (green), and left atrial appendage (red).(B) Definition of parallel planes (dashed line) dividing the left atrium.Coordinate ξ = 0 is defined with the plane that crosses the centroid of the ensemble of pulmonary veins inlets (CM P V s ) at a given instant of the cardiac cycle, and coordinate ξ = 1 with the plane that crosses the centroid of the mitral valve surface (CM MV ).

Fig. 4 .
Fig. 4. Visualization of the flow in the left atrium for the different flow split conditions during the atrial diastole (A) and during the ventricle diastole (B) in subject 2. For both instants, we represent the flow contained in a plane that crosses the atrial body and the flow inside the left atrial appendage.Panel (C) shows the temporal evolution of the flow rate through the mitral valve (red) and the pulmonary veins (blue), and indicates the position of the instants of (A) and (B) in the cardiac cycle.The right part of panel (C) shows the temporal evolution of the mean velocity through each pulmonary vein inlet for the different flow splits (R60, violet; R55, yellow; R50*, blue; eqV, red).

Fig. 5 .
Fig. 5. Visualization of the flow in the left atrium for the different flow split conditions during atrial diastole (A) and ventricular diastole (B) in subject 6.For both instants, we represent the flow contained in a plane that crosses the atrial body and the flow inside the left atrial appendage.Panel (C) shows the temporal evolution of the flow rate through the mitral valve (red) and the pulmonary veins (blue), and indicates the position of the instants of (A) and (B) in the cardiac cycle.The right part of panel (C) shows the temporal evolution of the mean velocity through each pulmonary vein inlet for the different flow splits (R60, violet; R55, yellow; R50*, blue; eqV, red).

Fig. 6 .
Fig. 6.Subject 2. Visualization of the instantaneous kinetic energy and residence time in an oblique plane section of the left atrium for the same two instants (A) and (B) of Fig. 4. Panel (C) shows the kinetic energy inside the LAA during the ventricle diastole.Panel (D) shows the temporal evolution of the mean kinetic energy and residence time inside the left atrial appendage along five cardiac cycles for the four flow splits (R60, violet; R55, yellow; R50*, blue; eqV, red).

Fig. 7 .
Fig. 7. Subject 6. Visualization of the instantaneous kinetic energy and residence time in an oblique plane section of the left atrium for the same two instants (A) and (B) of Fig. 4. Panel (C) shows the kinetic energy inside the LAA during the ventricle diastole.Panel (D) shows the temporal evolution of the mean kinetic energy and residence time inside the left atrial appendage along five cardiac cycles for the four flow splits (R60, violet; R55, yellow; R50*, blue; eqV, red).

Fig. 8 .
Fig. 8. Scatter plot with the correlation factor of the flow through the pulmonary veins versus the correlation factor of the flow inside the left atrium body in perpendicular (A) and parallel (B) directions to the mitral valve and inside the left atrial appendage (C).Symbols represent different cases (from case 1 to 8: (∘, ⎕, Δ. +, ×, ∇, ◇, *) The straight lines represent the regression lines for thrombus negative (blue) and thrombus positive (red) cases.

Fig. 9 .
Fig. 9. Probability distributions of kinetic energy, in perpendicular (A) and parallel (B) directions to the MV, and residence time (C) inside the left atrium body along three cardiac cycles (12 s < t < 15 s) Boxes and whiskers represent 9, 25, 50, 75 and 91 percentiles for each subject and each flow split condition (from left to right: R60, R55, R50* and eqV).Dashed line represents the residence time prediction given by T R = LAV/ α LVSV , with α = 1.4 and LAV and LVSV are respectively the mean LA volume and LA stroke volume.Mean intra flow split z-scores are represented on the top of each case.

Fig. 10 .
Fig. 10.Probability distributions of kinetic energy (A) and residence time (B) inside the left atrial appendage along three cardiac cycles (12 s < t < 15 s) .Boxes and whiskers represent 9, 25, 50, 75 and 91 percentiles for each subject and flow split condition (from left to right: R60, R55, R50* and eqV).Mean intra flow split z-scores are represented on the top of each case.

Fig. 11 .
Fig. 11.Mean value of the residence time in the left atrial appendage along three cardiac cycles (12 s < t < 15 s) for the four flow splits.Data in (A) correspond to the thrombus-negative cases, and the shaded patch in red corresponds to the envelope of the thrombus-positive cases.Data in (B) correspond to the thrombus-positive cases, and the shaded patch in blue corresponds to the thrombus-negative cases.The colors of the subjects are the same as in Figs. 8 and 9.
is the area of the MV, the peak velocity during the E-wave

Table 2
Median and Inter-Quartile Range of the correlation factors of the flow inside the LA body in perpendicular and parallel to the MV plane directions, and the flow inside the LAA for thrombus/TIA negative and positive cases.
Comput Biol Med.Author manuscript; available in PMC 2023 September 27.