Subject-Specific Calculation of Left Atrial Appendage Blood-Borne Particle Residence Time Distribution in Atrial Fibrillation

Atrial fibrillation (AF) is the most common arrhythmia that leads to thrombus formation, mostly in the left atrial appendage (LAA). The current standard of stratifying stroke risk, based on the CHA₂DS₂-VASc score, does not consider LAA morphology, and the clinically accepted LAA morphology-based classification is highly subjective. The aim of this study was to determine whether LAA blood-borne particle residence time distribution and the proposed quantitative index of LAA 3D geometry can add independent information to the CHA₂DS₂-VASc score. Data were collected from 16 AF subjects. Subject-specific measurements included left atrial (LA) and LAA 3D geometry obtained by cardiac computed tomography, cardiac output, and heart rate. We quantified 3D LAA appearance in terms of a novel LAA appearance complexity index (LAA-ACI). We employed computational fluid dynamics analysis and a systems-based approach to quantify residence time distribution and associated calculated variable (LAA mean residence time, t_m) in each subject. The LAA-ACI captured the subject-specific LAA 3D geometry in terms of a single number. LAA t_m varied significantly within a given LAA morphology as defined by the current subjective method and it was not simply a reflection of LAA geometry/appearance. In addition, LAA-ACI and LAA t_m varied significantly for a given CHA₂DS₂-VASc score, indicating that these two indices of stasis are not simply a reflection of the subjects' clinical status. We conclude that LAA-ACI and LAA t_m add independent information to the CHA₂DS₂-VASc score about stasis risk and thereby can potentially enhance its ability to stratify stroke risk in AF patients.

Introduction

Arial fibrillation (AF) is the most common arrhythmia, affecting three to six million US patients a year. This number is rapidly increasing with 12.1 million AF patients expected by 2030 (Virani et al., 2020). Due to the complex morphology of the left atrial appendage (LAA), as compared with the relatively smooth-walled left atrium (LA), the LAA is a favored location for thrombus formation (Reddy et al., 2013). These thrombi are known to cause stroke in AF patients. The CHA₂DS₂-VASc score is the most commonly used index for making clinical decisions regarding the management of AF patients. While this index is based on clinical data, it does not incorporate the role of LA–LAA geometry or local hemodynamics in the thromboembolic risk assessment.

The hypothesis that there is a correlation between the LAA morphology and stroke risk has been tested in several studies. Many indices have been examined in this context: LAA orifice diameter; number of branches and twigs; degree of coverage with fine structures (Ernst et al., 1995); LAA volume, depth, and number of lobes (Beinart et al., 2011); LAA takeoff from mitral valve (Nedios et al., 2014); and existence of a bend in LAA with an acute angle (Yaghi et al., 2020). Di Biase et al. (2012) categorized LAA shapes into four groups: chicken wing, windsock, cactus, and cauliflower shapes (Figure 1A). They concluded that patients with the chicken wing morphology are less likely to have a stroke. Although these results are promising, there is a large variability in stroke occurrence within a given LAA shape category (Khurram et al., 2013; Nedios et al., 2014; Sanatkhani and Menon, 2017; Yaghi et al., 2018). The subjective nature of LAA shape categorization may contribute to this variability. In an effort to objectify this, a recent study from our group has quantified LAA morphologies using principal component analysis (Sanatkhani and Menon, 2018). This approach uses the entire three-dimensional cardiac computed tomography (CCT) image, as opposed to isolated measurements of LAA dimensions, and therefore is more objective and comprehensive in quantifying LAA appearance.

FIGURE 1

Figure 1. (A) One representative 3D reconstructed left atrial appendage (LAA) surface geometry from our cohort from each of the four categories proposed by Di Biase et al. (2012). (B) Generic flow waveform over one cardiac cycle used to generate subject-specific waveforms based on subject's cardiac output and heart rate (adapted from Smiseth et al., 1999). To avoid tracers exiting left atrial inlets, the waveform is shifted such that the backflow occurs at the end of the cardiac cycle around t = 0.6 s. (C) Sample geometry with LAA filled with tracer (colored red). All inlets (four pulmonary veins), the outlet (mitral valve), and the location of LAA are displayed. (D) A representative left atrium and LAA geometry meshed with tetrahedrons. A section through the appendage region shows a finer mesh inside the appendage. Furthermore, it shows the use of prismatic layers at the wall boundary.

Several surrogates of thrombus-promoting flow patterns have been used to relate blood flow in vascular structures (including LA and LAA) to probability of clot formation: wall shear stress, shear strain rate, time-averaged wall shear stress, oscillatory shear index (Koizumi et al., 2015), time-averaged velocity, mean resident time (Rayz et al., 2010), local residence time (Esmaily-Moghadam et al., 2013), residual virtual contrast agent (Otani et al., 2016; Bosi et al., 2018), and vortex structure (Masci et al., 2019a,b). The most realistic solution to simulate clot formation is to model the transport of blood cells (i.e., platelets, red blood cells, etc.) in a geometry. A Lagrangian approach can be used for this purpose (Bernsdorf et al., 2006), which requires tracking of a large number of particles and a very fine mesh to resolve the flow field for particle tracking, making it computationally expensive (Rayz et al., 2010). A Eulerian approach, which approximates particle tracking, has been used with reasonable success for quantifying indices correlated with thrombus formation (Rayz et al., 2010; Esmaily-Moghadam et al., 2013; Otani et al., 2016; Bosi et al., 2018; Sanatkhani et al., 2018; Masci et al., 2019b). However, the Eulerian approach has not been employed alongside a systems model approach to correlate the patient-specific LAA geometries with the residence time in LAA and consequently thrombus formation probability.

The present study has two aims: (1) to introduce a Eulerian approach for calculating the blood-borne particles residence time distribution (RTD) in patient-specific LAA and (2) to evaluate the correlation between RTD and known physical (e.g., LAA appearance and location) or clinical (e.g., CHA₂DS₂-VASc score) indices. This patient-specific morphology and hemodynamics-based approach of calculating RTD may add novel value to the current methodology of stratifying stroke risk and therefore enhance AF management-related decisions.

Materials and Methods

Data Acquisition

We performed CCT in a cohort of 16 AF subjects prior to an AF catheter ablation procedure. CCT images were obtained using a multidetector 64-row helical system (Brilliance 64, Philips, Netherlands). Image acquisition was ECG-gated when possible with the following parameters: 70–120 kV, 850 mA s, 0.6 mm beam collimation, 0.625–1.25 mm thickness, and 20–30 cm field-of-view. During an end-inspiratory breath-hold of 20 s and following a timing bolus-chase injection (20 ml at 5 ml s⁻¹), 90 ml of an iodinated contrast medium (Ultravist 370, Bayer Vital, Cologne, Germany) was administered. Furthermore, subject-specific cardiac output and heart rate were measured. Heart rate and cardiac output data were used to scale (both the time and flow axes) a template LA inlet flow waveform (pulmonary vein flow) to generate a subject-specific LA inlet flow waveform (Figure 1B).

Image Segmentation

Contrast-enhanced CCT DICOM images were cropped and then smoothed using a median filter with a kernel of 5 × 5 × 5. Images were segmented in ParaView (version 5.9.0, Kitware, Inc., Albuquerque, NM, USA) using the marching cubes method to create an iso-surface, representing the LA surface. The extracted surface included pulmonary veins (PV), the LA and LAA walls, and the mitral valve plane (excluding the valves themselves; Figure 1C). Extracted surfaces were smoothed out for computational fluid dynamics mesh (Figure 1D) by removing spikes and reducing noise (i.e., simplifying polygons) in Geomagic Studio (version 10, Geomagic, Inc., Research Triangle Park, NC, USA) and ANSYS SpaceClaim (version 2020 R2, ANSYS Inc., Canonsburg, PA, USA). A more detailed flowchart of LAA segmentation is presented in Sanatkhani and Menon (2018).

LAA Appearance, LAA Location, and LAA–LA Size

LAA appearance was quantified according to the methodology described in Sanatkhani and Menon (2018). Briefly, principal component analysis was used to generate the eigenvectors of the LAA appearance of each LAA in our study cohort. Next, we reconstructed each subject-specific LAA using a successively increasing number of principal components (PCs) and calculated a normalized residual error in appearance reconstruction for each step, RE(i), as follows (Equation 1):

$\begin{array}{ll}RE(i)=\text{SQRT}\left[\frac{RSS(i)}{TSS}\right]& (1)\end{array}$

where RSS(i) is the residual sum of squares for the i^th step (i.e., i PCs used in the reconstruction) and TSS is the total sum of squares. RE(i) decreases with increasing i because the more PCs we use to reconstruct an image, the more information is available to describe the details of the original image. Using all PCs would result in the original LAA image with zero residual error [RE(16) = 0]. We have defined the area under the curve of RE(i) vs. i as the LAA appearance complexity index (LAA-ACI); a larger area under the curve would correspond to a more complex appearance (e.g., see Figure 2A).

FIGURE 2

Figure 2. (A) Normalized residual error, RE(i), plotted as a function of the number of principal components (i) used to reconstruct the LAA appearance for three subjects. As i increases, RE(i) decreases, reaching a value of zero when i = 16, corresponding to a perfect reconstruction. LAA appearance complexity index (LAA-ACI) is defined as the area the area under the RE(i)–i curve; larger LAA-ACI corresponds to a more complex LAA appearance. (B) Geometrical features of these three LAAs, including LAA ostium, tip of the LAA, LAA lobes (shown by circles), and LAA centerline bend (shown by curved arrows). The rank ordering of these three LAAs based on the LAA-ACI (most complex to simplest) is subject #3, subject #1, and subject #2.

The LAA location with respect to the LA inlet and LA outlet was characterized by the distance from the center of the LAA ostium to the four inlets (i.e., pulmonary veins; d_PV1, d_PV2, d_PV3, and d_PV4) and the outlet (i.e., mitral valve; d_Mitral). LAA volume (V_LAA) and LA volume (V_LA) were calculated using the 3D reconstructed CCT images. Distances are reported in centimeters and volumes in milliliters.

Computational Fluid Dynamics

The prepared geometries were meshed in ANSYS Meshing. The maximum length for the tetrahedron edge was considered 3.5 mm for the whole geometry (including LA and LAA). Then, mesh was refined based on surface curvature to capture the topology. For instance, mesh at the tip of the LAA is finer than at the center of LA. Five prismatic layers at wall boundaries were used to resolve the boundary layer flow (Figure 1D). We used these settings for our course mesh. Next, we incrementally increased the number of elements until the changes in averaged wall shear stress in LAA were <5%. The number of mesh elements was chosen based on our mesh independency study. Based on the size and tortuosity of each subject, the number of mesh elements varied between 300,000 and 500,000 tetrahedrons, which was considered acceptable according to the literature (Otani et al., 2016; Aguado et al., 2019). LA and LAA walls were assumed to be rigid, impermeable, and with no-slip boundary conditions. Furthermore, for simplicity and lowering the computational costs, the mitral valve was assumed to be open throughout the simulation with both gauge pressure and velocity gradient set to zero. To prevent outlet backflow divergence, we extended the outlet in our geometries (Figure 1C) to develop a uniform with zero gradient velocity and pressure at the outlets. A velocity profile was prescribed at PV inlet boundaries. The PV waveform was generated by modifying a template normal waveform (Figure 1B) based on subject-specific cardiac output and heart rate. Blood was treated as an incompressible, Newtonian fluid with a density of ρ = 1, 060 kg m⁻³ and dynamic viscosity of μ = 0.00371 Pa s (Formaggia et al., 2009).

We first ran the fluid dynamics simulations until a hemodynamic steady-state was reached (after 25 cycles), as defined by the steady state of wall shear stress averaged over the LAA surface area of each subject. Thereafter, we performed simulations to analyze the transport of virtual tracer [i.e., passive scalar, representative of blood-borne particles (cells) that are neutrally buoyant in plasma] out of LAA. These tracer transport-related simulations were initialized with the LAA filled with the tracer concentration, C(t), of unity (representing an impulse filling of LAA with the tracer) (Figure 1C). Tracer advection was simulated using fluid dynamic analysis where the tracer concentration of each cell was calculated in the transport equation coupled with the momentum equations. The volumetric average of tracer concentration inside the LAA was recorded as C(t) for 150 s (Figure 3A). Based on the decay characteristics of C(t), we fitted a triple exponential model to C(t) that included an asymptotic term, C_∞ (Equation A1).

FIGURE 3

Figure 3. (A) Spatially averaged LAA tracer concentration, C(t), plotted as a function of time for a representative LAA. Inset: Tracer concentration contours for selected time, illustrating the tracer washout from most of the LAA, except for the tip of the LAA. (B) E(t), the residence time distribution function quantifying the normalized rate of tracer washout across the LAA ostium, as a function of time for three representative subjects. Data for the first 10 s are shown to highlight the early washout. (C) Box plots showing the mean residence time, t_m, for each of the four traditional LAA shape groups. There is large variability in t_m within each traditional LAA shape group.

The dynamics of the tracer clearance from LAA was quantified in terms of the RTD function, E(t) (Fogler, 2016):

$\begin{array}{ll}\text{RTD Function}:E(t)=\frac{M(t)}{M_{total}}& (2)\end{array}$

where M(t) is the outflow of tracer material (amount of tracer material per unit time) from LAA at the LAA ostium and M_total is the total amount of tracer that will leave the LAA over the period 0 to infinity. Thus, E(t), with the unit per second, represents the normalized outflow of tracer material from LAA at time t. As shown in the Appendix, we can rewrite Equation 2 in terms of the LAA tracer concentration, C(t), as follows:

$\begin{array}{ll}E(t)=\frac{\left[C(t)-C(t+\Delta t)\right]}{\Delta t(1-C_{\infty })}& (3)\end{array}$

where Δt and C_∞ are the time increment used in the finite difference-based estimation of M(t) (Equation A3) and the asymptotic LAA concentration remaining in the LAA (Equation A1), respectively. Two measures of the propensity of particles to stay within the LAA were calculated: mean residence time, t_m, which is the first moment of E(t) (Equation A6), and C_∞ [C_∞ = C(t → ∞), Equation A1]. A larger value for either of these two indices is expected to increase the clot formation risk.

The tracer transport was considered to happen with 0 m² s⁻¹ diffusivity, making the transport purely due to advection. We solved the transient transport equation for a scalar (i.e., tracer) using a laminar solver developed from IcoFoam and ScalarTransportFoam solvers in OpenFOAM (version 8, The OpenFOAM Foundation Ltd, Inc., UK). A time step of 500 μs was chosen based on a time-step independence study. We started with a 2 ms time step and decreased the value until the changes in averaged wall shear stress in LAA were <5%. The first-order implicit method was used for time discretization and second-order least-square scheme was used for pressure and velocity gradient discretization. First-order and second-order upwind schemes were used to discretize the divergence terms in the scalar transport equation and the convection term in Navier–Stokes equations, respectively. Pressure, velocity, and concentration tolerances were set to be 10⁻⁷, 10⁻⁸, and 10⁻⁸, respectively. For these simulations, 16 cores of Intel Zeon CPU with 2.7 GHz clock speed and 8 GB of RAM were used and the average execution time for each case was ~35 h. C(t) was extracted in ParaView from LAA and postprocessed in MATLAB^® (version R2020b, MathWorks, Inc., Natick, MA).

Statistical Analysis

Continuous variables are expressed as mean ± standard deviation. Correlations between variables were determined by Spearman rank correlation. The relationship between LAA t_m (dependent variable) and several independent (predictor) variables (e.g., LAA and LA morphological indices) was analyzed using backward multiple linear regression, which yielded the coefficients (β) of significant predictor variables and the respective 95% confidence interval (CI). The threshold for statistical significance, α, was set to 0.05. All statistical analyses were performed in SAS software (version 9.4, SAS Institute, Inc., NC, USA).

Results

Study Subject Characteristics

A total of 16 subjects with symptomatic AF (eight paroxysmal, eight persistent) were studied (four samples per four LAA shapes; 11 males). The average age, heart rate, and cardiac output were 60.3 ± 11.1 years (range: 33–78 years), 71 bpm (range: 42–100 bpm), and 4.2 L min⁻¹ (range 2.7–6.0 L min⁻¹), respectively. The range of CHA₂DS₂-VASc scores was 0 to 4 (mean = 2.2 ± 1.1).

LAA Appearance Complexity Index (LAA-ACI)

The residual error for step i, RE(i), vs. i curves for three representative subjects is shown in Figure 2A, along with the geometry of the three LAAs (Figure 2B). Subject #3 requires significantly more PCs for accurate reconstruction and, therefore, has the largest LAA-ACI (i.e., area under the curve), indicating that this is the most complex appearance. The rank ordering based on the LAA-ACI is subject #3, subject #1, and subject #2. Interestingly, the rank ordering by the LAA mean residence time, t_m, is subject #1, subject #3, and subject #2, suggesting that the LAA-ACI and t_m are not conveying the same information.

Relating LAA RTD Function to Traditional LAA Shape Classification

The C(t) curve for a representative subject is shown in Figure 3A. The tracer washed out from the regions close to the LAA ostium after 2–10 s and tracer concentration continued to be high at the tip of the LAA even at the end of the simulation (t = 150 s).

The RTD function, E(t), for three representative subjects (same as those in Figure 2) is illustrated in Figure 3B, starting with the instant of the tracer introduction (t = 0) and until 10 s later. Subjects #1 and #3 started with lower initial normalized rate of tracer washout, E(t), across the LAA ostium, and this washout continued to remain lower as time progressed. Based on these E(t) curves and associated t_m values, the tracer exited from the LAA of subject #2 the fastest, followed by subjects #3 and subject #1 in that order.

The LAA appearance for each of the 16 subjects in the present study was classified into one of four groups by an electrophysiologist based on the study by Di Biase et al. (2012). Group data for LAA t_m (Figure 3C) indicated that LAA t_m had a large variability within each group, resulting in a significant overlap of this index of RTD function among the four LAA shape groups.

Relating LAA Mean Residence Time to LAA Asymptotic Tracer Concentration

The propensity of particles to stay within the LAA was characterized in terms of two indices: LAA mean residence time, t_m, and LAA asymptotic tracer concentration, C_∞. Spearman rank correlation analysis showed that there was a significant positive correlation between these two indices (R² = 0.78, P = 0.0003; Figure 4A), suggesting that only one of these indices may be sufficient to characterize the propensity of particles to stay within LAA; we choose LAA mean residence time, t_m.

FIGURE 4

Figure 4. (A) Relationship between the two indices of LAA residence time distribution (RTD) function: mean residence time, t_m, and aymptotic LAA concentration, C_∞. The Spearman rank correlation analysis indicates that LAA t_m and C_∞ are highly correlated, suggesting that only one of these indices is sufficent to characterizie the LAA RTD function. The three representative subjects shown in Figures 2, 3 are identified. (B) Relationship between LAA appearance complexity index (LAA-ACI) and LAA t_m. The Spearman rank correlation analysis indicates a weak (although statistically significant) correlation, suggesting that these two variables do not provide the same infoamtion. (C) Relationship between LAA-ACI and CHA₂DS₂-VASc score, analyzed using Spearman rank correlation analysis, showing a weak and insignificant correlation. (D) Relationship between LAA t_m and CHA₂DS₂-VASc score, analyzed using Spearman rank correlation analysis, showing a weak and insignificant correlation.

Relating LAA Mean Residence Time to LAA Morphology and Location

Given that the calculation of LAA t_m is computationally expensive, we wanted to examine whether LAA morphological indices (size, shape, and location) that are easier to measure can provide the same information. The Spearman rank correlation analysis indicated that only 56% of the variation in LAA t_m was explained by LAA-ACI (Figure 4B). We next carried out backward multiple linear regression analysis to investigate the relationship between LAA t_m (dependent variable) and eight independent variables: LAA-ACI, LAA location measures (d_PV1, d_PV2, d_PV3, and d_PV4, and d_Mitral; see LAA Appearance, LAA Location, and LAA–LA Size section), and LA and LAA volumes (V_LA and V_LAA). Only d_PV2 and V_LAA were found to be significant predictor variables for t_m. Although statistically significant, these two predictor variables could explain only 43% of the variation in t_m. Thus, the LAA morphological indices do not provide the same information as t_m.

Relating LAA Mean Residence Time and LAA-ACI to CHA₂DS₂-VASc Score

Spearman rank correlation analysis was performed between LAA-ACI and CHA₂DS₂-VASc score (Figure 4C) as well as between LAA t_m and CHA₂DS₂-VASc score (Figure 4D). Both LAA-ACI (R² = 0.46, P = 0.07; Figure 4C) and t_m (R² = 0.41, P = 0.11; Figure 4D) were positively correlated with CHA₂DS₂-VASc score, but all these correlations were weak and did not reach statistical significance.

Discussion

Lingering of blood cells inside the LAA could result in an elevated risk of thrombus formation and, consequently, stroke. In the present study, we quantify the propensity of blood cell lingering within the LAA in terms of the RTD function, E(t), and associated calculated variables (mean residence time of blood-borne particles in LAA, LAA t_m, and asymptotic concentration remaining inside LAA, C_∞). Both LAA and LA morphological features and spatially distributed hemodynamic milieu determine the LAA t_m and C_∞. The key contributions of the present study are as follows: (1) the quantitation of the overall LAA appearance in terms of a novel index, the LAA-ACI; (2) the development of a Eulerian and systems-based approach for quantifying the LAA RTD function and associated calculated variable (LAA t_m); (3) the observation that LAA-ACI and LAA location and size measures do not fully capture the information contained in LAA t_m; and (4) the observation that the LAA-ACI and LAA t_m can add independent information to the CHA₂DS₂-VASc score and thereby potentially enhance its ability to stratify stroke risk in AF patients.

The idea that the complexity of LAA geometry plays an important role in stroke risk stratification in AF subjects is not new. As discussed in the Introduction section, many indices have been examined in this context, e.g., LAA orifice diameter; number of branches and twigs; degree of coverage with fine structure; LAA volume, depth, and number of lobes; and existence of a bend in LAA (Ernst et al., 1995; Beinart et al., 2011; Di Biase et al., 2012; Khurram et al., 2013; Nedios et al., 2014). While these are isolated features of the complex LAA geometry, our LAA-ACI utilizes the entire 3D dataset, and we believe that this integrated index incorporates the information provided by isolated measures. Considering the LAA-ACI values for three representative subjects illustrated in Figure 2, the LAA-ACI for subject #2 is the lowest, indicating that this subject has the simplest appearance. Subject #2 has a smoother wall and its total length (length of the LAA centerline from the ostium to the tip) is shorter compared with the other two. The appearance of subject #3 is the most complex, having several lobes (circled regions in Figure 2B). Furthermore, subject #3 has the longest length with a bend along its centerline. Although the LAA of subject #1 is smooth, its ACI falls in between the other two subjects because of the large bend along its centerline. Thus, our LAA-ACI is an objective and quantitative metric that characterizes the complexity of LAA appearance in a holistic way. The next question is whether this integrated index is superior from the perspective of improving the stroke risk stratification. We cannot answer this question at the present time; a longitudinal study will be necessary to address this question.

Other studies have used the classification paradigm of Di Biase et al. (2012) to characterize LAA shape. This is a very subjective approach and even experienced cardiologists may not always agree when classifying a LAA into specific shape categories. In addition, there is a large variability in stroke occurrence within a given LAA shape category (Khurram et al., 2013; Nedios et al., 2014; Sanatkhani and Menon, 2017; Yaghi et al., 2018). This is consistent with the large variability of t_m (Figure 3C) and LAA-ACI (data not shown) within each shape category and in LAA values that exist even within a given LAA shape category (Figure 3C). The intercategorical t_m (and LAA-ACI) variability may explain the differences in the stroke risk seen among subjects with similar overall LAA geometry. This variability underscores the importance of considering subject-specific LA and LAA morphologies in constructing a metric for stroke risk stratification in AF based on hemodynamics.

A systems-based approach was used to calculate E(t) in that it is the tracer washout response to an impulse injection of tracer in LAA. LAA E(t) curves are depicted in Figure 3B for the same three subjects as in Figure 2. Subject #2 and subject #1 had the highest and lowest starting points (i.e., value at t = 0), with subject #3 having an intermediate value. The lower values in subjects #1 and #3 are a consequence of the flow entering the LAA ostium that does not go all the way up to the LAA tip, resulting in a stagnant region at the tip. Subject #2 had the LAA t_m. These data would predict that subject #1 has the highest risk of clot formation and subject #2 has the lowest risk.

There are other studies in the literature that have utilized more sophisticated computational fluid dynamics-based analysis to better mimic physiological conditions such as LA and LAA wall motion and more realistic LA outlet boundary conditions (Otani et al., 2016; Masci et al., 2019a). However, due to the very high computational costs associated with more complex computational fluid dynamic analyses, these studies were limited in terms of the number of subjects [five for Masci et al. (2019a) and four for Otani et al. (2016)] and the number of cardiac cycles analyzed (5 cycles). Based on our LAA mean residence time values (Figures 4A,B,D), we believe that it is necessary to perform the fluid dynamics simulation for at least 100 s. Our less sophisticated computational fluid dynamic analysis took about 35 h to complete an 150 s simulation (Intel Zeon CPU, 2.7 GHz, 16 cores). We believe that the more sophisticated analyses of Masci et al. (2019a) and Otani et al. (2016) will take at least several fold longer to perform an 150 s simulation, limiting their clinical application.

We have shown that the LAA t_m and LAA-ACI are weakly correlated (Figure 4B). This suggests that LAA t_m, representing a holistic measure of subject-specific LA–LAA geometry features and hemodynamics, and LAA-ACI have a potential to contribute independent information. The observation that d_Mitral and V_LA were significant predictors of LAA t_m is in line with clinical data from a recent publication (Nedios et al., 2014) showing that a higher LAA takeoff, remote to the mitral valve plane, was associated with an increased thromboembolic risk.

We found a weak and insignificant correlation between the LAA-ACI and CHA₂DS₂-VASc score (Figure 4C) and between LAA t_m and CHA₂DS₂-VASc score (Figure 4D). In addition, LAA t_m varied significantly for a given CHA₂DS₂-VASc score (Figure 4D). Although subject #1 has a low CHA₂DS₂-VASc score (=1), this subject has the highest LAA t_m (Figure 4D). We would suggest that subject #1 has a high risk of stroke, despite the low CHA₂DS₂-VASc score. In contrast, there are two other subjects with CHA₂DS₂-VASc score of 1 and relatively low values of LAA t_m (Figure 4D); we would suggest that these subjects have a very low risk of stroke. These observations suggest that the hemodynamics-based index (i.e., LAA t_m) and appearance indices (i.e., LAA-ACI) can add independent information to the CHA₂DS₂-VASc score.

Limitations

Although we used subject-specific LA and LAA geometries, two concerns can be raised regarding the subject specificity of the fluid dynamics-based analysis. First, for LA inlet boundary conditions, a template (generic) LA inlet flow (pulmonary vein flow) waveform was used instead of a subject-specific waveform because patient-specific waveforms were not available for our cohort. However, because our preliminary study (Sanatkhani et al., 2020) showed that LAA t_m is significantly affected by mean LAA inlet flow (i.e., cardiac output), we scaled the template waveform of each subject to match the subject-specific cardiac output and heart rate. Future parametric studies will examine whether the temporal features of LAA inlet flow affect LAA t_m. Second, for LA outlet boundary condition, the mitral valve was assumed to be open throughout the simulation, with both gauge pressure and velocity gradient set to zero. A better representation of LA outlet boundary condition would be in terms of left ventricular non-linear diastolic compliance and patient-specific left ventricular end-diastolic or end-systolic volume. These left ventricular diastolic compliance and volume data were not available for our cohort. We plan to conduct parametric studies to examine how left ventricular compliance and end-diastolic (or end-systolic) volume affect the calculated index, LAA t_m.

We used the assumption of rigid LA and LAA walls. This assumption is an approximation for highly quivering LA and LAA in persistent AF patients, corresponding to the loss of coordinated contraction and consequent reduction in wall motion. A more realistic simulation with compliant LA and LAA walls would require the fluid–solid interaction approach for hemodynamic simulations, which significantly increases the complexity of the computational fluid dynamics analysis and associated computational cost. Alternatively, one can prescribe LA and LAA wall motions as boundary conditions; however, such wall motion data were not available for our cohort. It is recognized that some studies have used the prescribed wall motion in the computational fluid dynamics analysis of LA–LAA hemodynamics (Otani et al., 2016; Masci et al., 2019a). However, this approach requires additional data acquisition and significantly increases the computational costs, which may limit the clinical applicability. It has been shown that LAA thrombus formation is associated with both poor LAA contraction and LAA dilation (Pollick and Taylor, 1991). Therefore, the complete loss of wall motion may overestimate the risk of thrombus formation inside the LAA.

Furthermore, we assumed the blood to be a Newtonian fluid considering that LA domain is large and shear rate values are mostly in the range that allows for this assumption (mean shear strain rate, $\overline{\dot{\gamma }}=28.3\pm 45.9s^{-1}$). However, to better simulate the stasis regions inside the LAA, one might need to consider a non-Newtonian model. In our ongoing study, we are investigating the sensitivity of LAA t_m to various hematocrit using a non-Newtonian model in comparison with the Newtonian model.

Although our results indicate that LAA t_m and ACI provide additional information, a longitudinal study with a larger number of subjects will be needed to examine whether adding these indices to the CHA₂DS₂-VASc score can indeed enhance stroke prediction in AF.

Summary and Conclusions

We have presented a novel index for quantifying the LAA geometry (LAA-ACI) and an approach for quantifying LAA residence time distribution and associated calculated variables using the subject-specific morphology, cardiac output, and heart rate with a hemodynamic model. Both the appearance index (i.e., LAA-ACI) and the hemodynamics-based index (i.e., LAA t_m) add independent information to the CHA₂DS₂-VASc score about subject-specific stasis risk and thereby can potentially enhance its ability to stratify stroke risk in AF patients.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics Statement

Ethical review and approval was not required for the study on human participants in accordance with the local legislation and institutional requirements. Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.

Author Contributions

SS, PM, and SGS designed the model and the computational framework and analyzed the data. SS performed the computations. SN, AB, and GH acquired the cardiac-CT data. SS, PM, and SGS wrote the manuscript with input from all authors. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by Gerald E. McGinnis Endowed Chair funds (SGS), American Heart Association grant #20PRE35210594 (SS), and Leonard H. Berenfield Graduate Fellowship in Cardiovascular Bioengineering (SS).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors would like to thank the University of Pittsburgh Center for Research Computing for providing resources.

Abbreviations

AF, atrial fibrillation; LAA, left atrial appendage; LA, left atrium; RTD, residence time distribution; CCT, cardiac computed tomography; DICOM, digital imaging and communications in medicine; PV, pulmonary vein; PC, principal component; RE(i), residual error; RSS(i), residual sum of squares; TSS, total sum of squares; ACI, appearance complexity index; d_PV1, d_PV2, d_PV3, and d_PV4, distance from the center of the LAA ostium to the four inlets (i.e., pulmonary veins), respectively; d_Mitral, distance from the center of the LAA ostium to the mitral valve; V_LA, LA volume; V_LAA, LAA volume; C′(t), concentration of tracer while exiting the LAA; C(t), tracer concentration inside the LAA; E(t), residence time distribution function; t_m, mean residence time; C_∞, asymptotic concentration remaining inside LAA; Q(t), volumetric flow exiting the LAA; β, multiple linear regression coefficient; CI, confidence interval; M(t), outflow of tracer material from LAA at the LAA ostium; M_total, total amount of tracer that will leave the LAA over the period 0 to infinity.

References

Aguado, A. M., Olivares, A. L., Yagüe, C., Silva, E., Nuñez-García, M., Fernandez-Quilez, Á., et al. (2019). In silico optimization of left atrial appendage occluder implantation using interactive and modeling tools. Front. Physiol. 10:237. doi: 10.3389/fphys.2019.00237

PubMed Abstract | CrossRef Full Text | Google Scholar

Beinart, R., Heist, E. K., Newell, J. B., Holmvang, G., Ruskin, J. N., and Mansour, M. (2011). Left atrial appendage dimensions predict the risk of stroke/TIA in patients with atrial fibrillation. J. Cardiovasc. Electrophysiol. 22, 10–15. doi: 10.1111/j.1540-8167.2010.01854.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Bernsdorf, J., Harrison, S. E., Smith, S. M., Lawford, P. V., and Hose, D. R. (2006). Numerical simulation of clotting processes: a lattice boltzmann application in medical physics. Math. Comput. Simul. 72, 89–92. doi: 10.1016/j.matcom.2006.05.008

CrossRef Full Text | Google Scholar

Bosi, G. M., Cook, A., Rai, R., Menezes, L. J., Schievano, S., Torii, R., et al. (2018). Computational fluid dynamic analysis of the left atrial appendage to predict thrombosis risk. Front. Cardiovasc. Med. 5:34. doi: 10.3389/fcvm.2018.00034

PubMed Abstract | CrossRef Full Text | Google Scholar

Di Biase, L., Santangeli, P., Anselmino, M., Mohanty, P., Salvetti, I., Gili, S., et al. (2012). Does the left atrial appendage morphology correlate with the risk of stroke in patients with atrial fibrillation?: results from a multicenter study. J. Am. Coll. Cardiol. 60, 531–538. doi: 10.1016/j.jacc.2012.04.032

PubMed Abstract | CrossRef Full Text | Google Scholar

Ernst, G., Stöllberger, C., Abzieher, F., Veit-Dirscherl, W., Bonner, E., Bibus, B., et al. (1995). Morphology of the left atrial appendage. Anat. Rec. 242, 553–561. doi: 10.1002/ar.1092420411

CrossRef Full Text | Google Scholar

Esmaily-Moghadam, M., Hsia, T.-Y., and Marsden, A. L. (2013). A non-discrete method for computation of residence time in fluid mechanics simulations. Phys. Fluids 25:110802-. doi: 10.1063/1.4819142

PubMed Abstract | CrossRef Full Text | Google Scholar

Fogler, H. S. (2016). “Residence time distributions of chemical reactors,” in Elements of Chemical Reaction Engineering, ed H. S. Fogler (Boston, MA: Prentice Hall), 767–806. Available online at: https://learning.oreilly.com/library/view/elements-of-chemical/9780135486252/

Formaggia, L., Quarteroni, A., and Veneziani, A. (2009). “Cardiovascular mathematics: modeling and simulation of the circulatory system,” in Cardiovascular Mathematics, eds L. Formaggia, A. Quarteroni, and A. Veneziani (Milano: Springer Science and Business Media), 22–24.

Google Scholar

Khurram, I. M., Dewire, J., Mager, M., Maqbool, F., Zimmerman, S. L., Zipunnikov, V., et al. (2013). Relationship between left atrial appendage morphology and stroke in patients with atrial fibrillation. Heart Rhythm 10, 1843–1849. doi: 10.1016/j.hrthm.2013.09.065

PubMed Abstract | CrossRef Full Text | Google Scholar

Koizumi, R., Funamoto, K., Hayase, T., Kanke, Y., Shibata, M., Shiraishi, Y., et al. (2015). Numerical analysis of hemodynamic changes in the left atrium due to atrial fibrillation. J. Biomech. 48, 472–478. doi: 10.1016/j.jbiomech.2014.12.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Masci, A., Alessandrini, M., Forti, D., Menghini, F., Dedé, L., Tomasi, C., et al. (2019b). A proof of concept for computational fluid dynamic analysis of the left atrium in atrial fibrillation on a patient-specific basis. J. Biomech. Eng. 142:011002. doi: 10.1115/1.4044583

CrossRef Full Text | Google Scholar

Masci, A., Barone, L., Dedè, L., Fedele, M., Tomasi, C., Quarteroni, A., et al. (2019a). The impact of left atrium appendage morphology on stroke risk assessment in atrial fibrillation: a computational fluid dynamics study. Front. Physiol. 9:1938. doi: 10.3389/fphys.2018.01938

PubMed Abstract | CrossRef Full Text | Google Scholar

Nedios, S., Kornej, J., Koutalas, E., Bertagnolli, L., Kosiuk, J., Rolf, S., et al. (2014). Left atrial appendage morphology and thromboembolic risk after catheter ablation for atrial fibrillation. Heart Rhythm 11, 2239–2246. doi: 10.1016/j.hrthm.2014.08.016

CrossRef Full Text | Google Scholar

Otani, T., Al-Issa, A., Pourmorteza, A., McVeigh, E. R., Wada, S., and Ashikaga, H. (2016). A computational framework for personalized blood flow analysis in the human left atrium. Ann. Biomed. Eng. 44, 3284–3294. doi: 10.1007/s10439-016-1590-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Pollick, C., and Taylor, D. (1991). Assessment of left atrial appendage function by transesophageal echocardiography. Implications for the development of thrombus. Circulation 84, 223–231. doi: 10.1161/01.CIR.84.1.223

PubMed Abstract | CrossRef Full Text | Google Scholar

Rayz, V. L., Boussel, L., Ge, L., Leach, J. R., Martin, A. J., Lawton, M. T., et al. (2010). Flow residence time and regions of intraluminal thrombus deposition in intracranial aneurysms. Ann. Biomed. Eng. 38, 3058–3069. doi: 10.1007/s10439-010-0065-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Reddy, V. Y., Doshi, S. K., Sievert, H., Buchbinder, M., Neuzil, P., Huber, K., et al. (2013). Percutaneous left atrial appendage closure for stroke prophylaxis in patients with atrial fibrillation. Circulation 127, 720–729. doi: 10.1161/CIRCULATIONAHA.112.114389

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanatkhani, S., and Menon, P. G. (2017). “Relating atrial appendage flow stasis risk from computational fluid dynamics to imaging based appearance paradigms for cardioembolic risk,” in Imaging for Patient-Customized Simulations and Systems for Point-of-Care Ultrasound; Lecture Notes in Computer Science, ed M. J. Cardoso (Québec City, QC: Springer), 86–93.

Google Scholar

Sanatkhani, S., and Menon, P. G. (2018). “Generative statistical modeling of left atrial appendage appearance to substantiate clinical paradigms for stroke risk stratification,” in SPIE Medical Imaging (Houston, TX).

Google Scholar

Sanatkhani, S., Menon, P. G., Shroff, S. G., and Sotirios, N. (2018). Abstract 11255: thrombus risk prediction in atrial fibrillation: hemodynamic model of left atrial appendage. AHA Scientific Sessions. Circulation 138:A11255. doi: 10.1161/circ.138.suppl_1.11255

CrossRef Full Text | Google Scholar

Sanatkhani, S., Nedios, S., Jain, S. K., Saba, S., Menon, P. G., and Shroff, S. G. (2020). Abstract 16439: is pulmonary venous flow pulsatility a critical determinant of left atrial appendage blood stasis risk? AHA Scientific Sessions. Circulation 142:A16439. doi: 10.1161/circ.142.suppl_3.16439

CrossRef Full Text | Google Scholar

Smiseth, O. A., Thompson, C. R., Lohavanichbutr, K., Ling, H., Abel, J. G., Miyagishima, R. T., et al. (1999). The pulmonary venous systolic flow pulse–its origin and relationship to left atrial pressure. J. Am. Coll. Cardiol. 34, 802–809. doi: 10.1016/S0735-1097(99)00300-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Virani, S. S., Alonso, A., Benjamin, E. J., Bittencourt, M. S., Callaway, C. W., Carson, A. P., et al. (2020). Heart disease and stroke statistics −2020 update: a report from the American Heart Association. Circulation 141, e139–e596. doi: 10.1161/CIR.0000000000000757

PubMed Abstract | CrossRef Full Text | Google Scholar

Yaghi, S., Chang, A. D., Akiki, R., Collins, S., Novack, T., Hemendinger, M., et al. (2020). The left atrial appendage morphology is associated with embolic stroke subtypes using a simple classification system: a proof of concept study. J. Cardiovasc. Comput. Tomogr. 14, 27–33. doi: 10.1016/j.jcct.2019.04.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Yaghi, S., Chang, A. D., Hung, P., Mac Grory, B., Collins, S., Gupta, A., et al. (2018). Left atrial appendage morphology and embolic stroke of undetermined source: a cross-sectional multicenter pilot study. J. Stroke Cerebrovasc. Dis. 27, 1497–1501. doi: 10.1016/j.jstrokecerebrovasdis.2017.12.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Appendix

The computational fluid dynamics model yields LAA tracer concentration, C(t), every 1 s over the entire simulation period t = 0, 150 s (Figure 3A). The C(t) decay curve seems to have mutiple time constants and a non-zero asymptotic value (Figure 3A). Accordignly, a triple exponential model with non-zero asymptote was used to characterize the C(t) curve:

$\begin{array}{ll}C(t)=a_1{e}^{-b_{1}t}+a_2{e}^{-b_{2}t}+a_3{e}^{-b_{3}t}+C_{\infty }& (\mathrm{A1})\end{array}$

where a₁, a₂, a₃, b₁, b₂, and b₃ are the parameters to be estimated by fitting the triple exponential model (Equation A1) to the calculated C(t) data. Given that C(t = 0) = 1, C_∞ = 1 − a₁ − a₂ − a₃.

The dynamics of the tracer clearance from LAA was quantified in terms of the RTD function, E(t) (Fogler, 2016):

$\begin{array}{ll}\text{RTD Function}:E(t)=\frac{M(t)}{M_{total}}& (\mathrm{A2})\end{array}$

where M(t) is the outflow of tracer material (amount of tracer material per unit time) from LAA at the LAA ostium and M_total is the total amount of tracer that will leave the LAA over the period 0 to infinity. M(t) can be written in terms of C(t) using the finite difference approximation:

$\begin{array}{ll}M(t)=\frac{[C(t)-C(t+t)]V_{LAA}}{t}& (\mathrm{A3})\end{array}$

where Δt and V_LAA are the time increment used in the finite difference-based estimation of M(t) and the LAA volume, respectively. To minimize errors in the finite difference-based estimation of M(t), we chose Δt = 10^–6 s.

Given that C(t = 0) = 1 and C(t = ∞) = C_∞, M_total is given by:

$\begin{array}{ll}{M}_{total}=(1-C_{\infty })V_{LAA}& (\mathrm{A4})\end{array}$

Substituting Equation A3 and Equation A4 in Equation A2, we get:

$\begin{array}{ll}E(t)=\frac{\left[C(t)-C(t+\Delta t)\right]}{\Delta t(1-C_{\infty })}& (\mathrm{A5})\end{array}$

E(t)Δt is the fraction of fluid exiting the LAA that has spent between time t and t + Δt inside the LAA. Thus, the LAA mean residence time, t_m, is given by the the first moment of E(t):

$\begin{array}{l}{t}_m={\displaystyle {\int }_0^{\infty }}tE(t)dt=\frac{1}{\Delta t(1-C_{\infty })}\\ {\displaystyle {\int }_0^{\infty }}\left[tC(t)-tC(t+\Delta t)\right]dt& (\mathrm{A6})\end{array}$

It should be noted that since C(t) is comprised of exponetially decaying terms, ${\int }_0^{\infty }tC(t)dt$ is finite and the result can be written in terms of an analytical expression. For example, ${\int }_0^{\infty }t{e}^{-at}dt=a^{-2}$. A custom program was developed in the MATLAB^® (version R2020b, MathWorks, Inc., Natick, MA, USA) environment to estiamte the six parameters (a₁, a₂, a₃, b₁, b₂, and b₃) using a non-linear iterative optimization algorithm and perform subsequent data processing to calculate LAA mean residence time, t_m.