Numerical prediction of thrombosis risk in left atrium under atrial fibrillation

1.   Introduction

Atrial fibrillation (AF) is a common cardiac arrhythmia, during which the regular and orderly electrical activity in the atrium is replaced by rapid and chaotic heartbeats. The ineffective contraction of the atrium in the post-diastolic phase during AF can induce thrombosis, while the shed thrombus can cause ischemic stroke when flowing to the brain. The presence of AF is independently associated with a 5-fold increased risk of stroke ^[1]. The left atrial appendage (LAA) is a finger- or stump-like extension of the left atrium (LA) with lobes that may harbor up to 90% of thrombi that occur in patients with AF ^[2]. The remodeling process associated with AF causes the LAA to function as a static pouch, predisposing to stagnation and thrombosis. It also should be noticed that the proportion of thrombus in LAA as the origin of cardiogenic emboli is as high as 100% for non-vascular AF patients ^[2].

Prevention strategies of ischemic stroke caused by AF include oral anticoagulation therapy (i.e., warfarin) and minimally invasive surgeries to exclude LAA ^[3,4,5]. Before anticoagulation therapy, the risk of stroke is stratified according to the CHA₂DS₂-VASc score which refers to age, sex, stroke history and so on. A higher score means a higher risk of stroke ^[6]. CHA₂DS₂-VASc score ≡2 means a patient at high risk and needs anticoagulation therapy. However, the scoring system still has some limitations. For example, a young male with AF whose CHA₂DS₂-VASc score is 0, still could develop ischemic stroke. Noticed that some clinical studies have shown that LA and LAA hemodynamic information which is difficult to obtain through current medical technology, can improve stroke risk stratification ^[6], thus computational fluid dynamics (CFD) is widely applied to capture blood flow characteristics ^[7,8,9]. The study of the LA wall motion indicates that LAA is ineffective and acts as a stasis blood reservoir during AF ^[9]. The blood evacuation rate in LAA decreases when AF occurs which might increase thrombosis risk ^[8]. Comparing the lack of active atrial contraction, the occurrence of high-frequency fibrillation might have a larger impact on the blood flow stagnation ^[10]. We find that only a few experimental and numerical research study the flow patterns in LA ^[11] and hemodynamic changes induced by different LAA morphology or AF ^[7,8,9,12]. At the same time, the direct simulation of thrombosis in the LA, especially LAA, has yet to be produced.

The purpose of this study was to simulate the thrombosis in LA using advanced computational modeling analyses. To this end, the patient-specific geometry was reconstructed and a thrombus model was combined with constant inlet pressure and realistic outflow. To the author*s knowledge, the present study is the first to use the CFD method to directly model thrombosis in patient-specific LA. The magnitudes of thrombosis variables, blood flow distribution and wall shear stress (WSS) related indices were analyzed to explore the impact of AF on the hemodynamics.

2.   Materials and method

2.1.   Geometric model

Computed tomography (CT) images of LA and four pulmonary veins (PVs) were obtained from an adult patient with non-valvular AF at Sir Run Run Shaw Hospital, Zhejiang University School of Medicine (Hangzhou, China). The study was conducted in accordance with the ethical standards of the Local Institution Review Board and with the 1964 Helsinki declaration and its later amendments. This project was approved by the Local Bioethics Committee. All patients have given their written informed consent before the study. The patient-specific LA geometry was constructed from the CT images by using Mimics 19.0 (Materialise, Belgium) (Figure 1-A). All the geometric boundaries of inlets and outlet were cropped to get a flat surface in GeoMagic Studio (GeoMagic Inc, USA) (Figure 1-B). The fluid domain was meshed using commercial software ANSYS-ICEM (ANSYS Inc, Canonsburg, USA), while five prism layers of unstructured tetrahedral meshes were created close to the atrial wall. The computational domain had approximately 2,420,000 elements, and the unstructured tetrahedral mesh had a maximum size of 7 mm. Finer meshes with more than 4,900,000 elements were created to assess the mesh sensitivity. The differences in peak wall shear stress and averaged wall shear stress between meshes were both less than 2%.

Figure 1.  Geometry models and blood flow curves. A: Preoperative and postoperative computed tomography images; B: Three-dimensional reconstruction of the left atrial geometries; C: Blood flow curve for atrial fibrillation cases and sinus rhythm cases; D: Geometry for backward-facing step model.

DownLoad: Full-Size Img PowerPoint

2.2.   Numerical model and computational details

The blood was considered to be an incompressible and non-Newtonian fluid with a constant density of 1080 kg/m^{3 [13]}. The LA wall was assumed to be rigid with no-slip conditions. The deformation of the LA should be considered for precise results in most heart flow simulations. However, the focus of the present study is on the application of a thrombus model in LA to simulate thrombosis. In addition, the effect of the atrial kick absence, a feature of the AF, on thrombosis is also explored. To build a fully developed blood flow environment before thrombus model was employed ^[13], Quemada model was used to simulate the shear-thinning characteristics of the blood ^[14] for two complete cardiac cycles (each cardiac cycle of 0.8 s duration including 37.5% ventricular diastole and 62.5% ventricular systole). The thrombus model was developed based on a hemodynamics-based model presented by Menichini, et al. ^[15], which simulated thrombosis by solving transport equations of flow residence time (RT), activated platelets (AP), resting platelets (RP), bound platelets (BP) and coagulation (C). The thrombus model considers the effect of the thrombus on the blood viscosity as well. Thrombus is identified through the local BP concentration. Details of the mathematical equations can be found in the study of Menichini, et al. ^[15], but there are some important changes described below.

The momentum of the flow is controlled by the simplified Navier-Stokes equation. We ignore the frictional force of thrombus on the fluid. For the transport of AP and RP, a non-normalized source item S_i is applied which is the summary of $r_{1}$ and $r_{2}$. BP and C have been normalized by their initial values for avoiding numerical ill-conditioning. Besides, the judgment threshold value of the local time average shear strain rate is 1 s^-1, which determines that the initial concentration of C is 100 or 0 nmol/L.

$k_{1}$ and $k_{2}$ are the kinetic constants ^[16]. AP and RP are the concentrations of the platelets respectively. RRT is the standardized RT which is normalized by cardiac cycle time.

The difficulty in obtaining patient-specific flowrate data in the outlet, mitral valve (MV), would be revealed when using medical measurement. The flowrate curves across the MV for the sinus rhythm (SR) cases were derived from the international regulation, which included E-wave and A-wave representing the first and the second phases, respectively ^[17]. The transmitral Doppler echocardiographic analysis of AF patients shows that A wave which is caused by the atrial kick in late diastole is absent in the flowrate curve ^[10,18]. The flowrate curves for AF cases were differed by dismissing the A wave ^[7]. For the four pulmonary veins inlets, constant pressure of 10 mmHg was adopted. Flowrate curves for SR and AF cases are shown in Figure 1-C, the peaks of two curves are nearly 46.9 ml/s. The maximum Reynolds number of all cases was found lower than 3000 when estimated on the basis of the peak flow rate and the MV diameter. Therefore, the blood flow was assumed as a laminar flow. Simulations were performed on ANSYS Workbench 16.1 (ANSYS Inc, USA). The time step in this transient study was 10 ms and all the cases were simulated for 20 cardiac cycles to get feasible results and the twentieth cycle data was post-processed.

3.   Results

3.1.   Verification of thrombus model

Before applied to LA geometry, the thrombus model was tested in a backward-facing step (BFS) model (Figure 1-D). Quemada model with a 30% hematocrit was selected to simulate experimental blood ^[19], while in the LA model, the hematocrit of human blood was determined to be nearly 45% ^[14]. The boundary conditions were adopted from the experimental test, with 0.76 L/min at the inlets and zero pressure at the outlet ^[19]. Except for the bulk shear rate threshold, all other parameters were the same as in the LA model. The time step in this verified study was 0.05 s and the calculation time was 50 s. Verification results were compared with the experimental study ^[19] and the simulation results ^[15].

Figure 2 shows the thrombosis prediction results in the BFS model and they are compared with the published results. The thrombus initially appears in the downstream of the step then grows continuously. When the height of the thrombus reaches the step height (2.5 mm) and the maximum width is close to that of the step, thrombus growth almost stops in the radial direction. The thrombus growth turns slow after 14 s and almost stops growing when the total thrombus length reaches about 8.8 times the height of the step, with a maximum length of 22 mm. These results are consistent with the simulation data of Menichini, et al. ^[15] and match well with the experimental results ^[19] as we can use it for further study.

Figure 2.  Thrombosis in the backward-facing step model (red for thrombus). A/B/C: Simulation results of our study for 14 s/28 s/50 s respectively; D/E/F: Simulation results for 14 s/28 s/50 s respectively. Reprint from Ref ^[15] with open access; G: The height and length of thrombus at t = 50 s.

DownLoad: Full-Size Img PowerPoint

3.2.   Blood flow and velocity characteristics

The velocity streamlines for two cases are illustrated in Figure 3. The initial streamlines seeds are averagely launched from PVs. At the first peak systole (0.4 s), the streamline distributions in the two cases are similar, the velocities in the inlets (PVs) and MV outlet are higher than in other areas of the LA. During the later peak systole (0.7 s), the blood flow gets highly disordered in the recirculating area, and the maximum velocity appears nearby the left inferior pulmonary vein (LIPV) rather than the PVs and MV in the AF case. Besides, there is little blood flowing through the MV in the AF case while blood can in the SR case. Mentioned that there is little blood flowing out from the LAA as well. Table 1 shows the proportions of blood that pass through each inlet boundaries during the last cardiac cycle. During the later peak systole for AF case, the flow ratios of the four inlets were not calculated because part of the inlet blood returned. In the SR case, the blood flow crossing the LIPV accounts for nearly 30% of the outflow while the minority of the blood (below 20%) goes through the right inferior pulmonary vein.

Figure 3.  The distributions of velocity streamlines in the left atrium for sinus rhythm case (A/B) and atrial fibrillation case (C/D) at two systole peaks (A/C: t = 0.4 s, B/D: t = 0.7 s). The minimum and maximum velocities are also recorded.

DownLoad: Full-Size Img PowerPoint

Table 1.  The flow ratio of each inlet boundaries to the outlet flow for the last cardiac cycle.

Inlet boundaries	0.4 s		0.7 s
	AF	SR	AF	SR
Left Superior Pulmonary Vein	26.86%	27.88%	-	24.49%
Left Inferior Pulmonary Vein	30.27%	30.56%	-	29.02%
Right Superior Pulmonary Vein	26.14%	25.35%	-	27.20%
Right Inferior Pulmonary Vein	16.74%	16.22%	-	19.29%
AF: atrial fibrillation; SR: sinus rhythm.

---

 | Show Table

DownLoad: CSV

3.3.   Wall shear stress

The frictional force on the LA wall exerted by the blood flow is quite difficult to be obtained through current measurement technology but can be studied by the WSS. The distributions of WSS for two cases in peak systole phases are shown in Figure 4. Similar to the flow waveform, the distributions of WSS in two cases have little difference at the earlier peak systole. At the later peak systole, the maximum WSS values of two cases decrease greatly, especially in the AF case. Besides, the region around LIPV shows a bigger WSS than the other three inlets in the AF case (Figure 4-D), while areas around the four PVs have a similar distribution in the other three maps. The LAA has a small WSS value in all maps.

Figure 4.  The distributions of WSS in the left atrium for sinus rhythm case (A/B) and atrial fibrillation case (C/D) at two systole peaks (A/C: t = 0.4 s, B/D: t = 0.7 s). The minimum and maximum WSS are also listed. (WSS: wall shear stress).

DownLoad: Full-Size Img PowerPoint

3.4.   Thrombosis prediction

We choose BP as a thrombosis sign with a threshold of 200 nmol/L ^[20]. Figure 5 shows the thrombosis in the maximum cross-section of the LAA from 5.6 s with 1.6 s interval for AF case. The simulation data illustrates that thrombus forms since the seventh cardiac cycle. The thrombus first appears on the wall of the LAA, then grows inward and fills most space in the LAA at the end of the nineteenth cardiac cycle. The distributions and volume proportions (comparing with the total LA volume) of formed thrombus in the seventh and last cycles of two cases are shown in Figure 6. Thrombus only arises in the LAA with its proportions have increased tremendously for two cases from the eighth to the twentieth cycle. The difference of thrombosis ratio is clear between two cases at 16 s while AF case has a greater severe thrombosis with a 4.93% thrombus volume proportion.

Figure 5.  The thrombosis in the maximum cross-section of the left atrial appendage from 5.6 s with 1.6 s interval in the atrial fibrillation case.

DownLoad: Full-Size Img PowerPoint

Figure 6.  The distributions and thrombus volume proportions comparing with the total left atrium volume. A/B: at the end of the eighth cycle for SR/AF cases; C/D: at the end of the last cycle for SR/AF cases; E: The distributions of bound platelets in two sections of the left atrial appendage. (AF: atrial fibrillation; SR: sinus rhythm).

DownLoad: Full-Size Img PowerPoint

The formation of BP is directly influenced by C, whose initial value is depending on the distribution of the average shear rate on the wall after the first two cardiac cycles. If the average shear rate is smaller than the judgment threshold, the initial value of C is 200 nmol/L ^[15]. In view of the importance of C initial value for outcomes, thrombosis in AF cases with judgment thresholds of 0.5 s^-1 and 1.5 s^-1 were calculated as well (Figure 7). For all cases with different thresholds, the thrombus is visible since the seventh cardiac cycle. In the end, the volume fraction of thrombus reaches 5.68% in the case with a threshold of 1.5 s^-1, which is nearly twice that of the case with the judgment threshold of 0.5 s^-1. Other than that, the relative thrombus fraction difference between cases with thresholds of 1.5 s^-1 and 1 s^-1 is 15.2%, which is smaller than that value between cases with thresholds of 1 s^-1 and 0.5 s^-1. Figure 7-D investigates the time-varying laws of thrombus proportions for three cases. The growth of thrombosis accelerates with the increase of judgment threshold and decelerates with the change of time. Furthermore, the difference of thrombus growth rates between cases with judgment thresholds of 1.5 s^-1 and 1 s^-1 is smaller than that between the latter and the case with a threshold of 0.5 s^-1.

Figure 7.  The distributions and thrombus volume proportions for atrial fibrillation cases comparing with the total left atrium volume. A1/A2: at the end of the eighth/last cycle with judgment thresholds of 0.5 s^-1; B1/B2: at the end of the eighth/last cycle with judgment thresholds of 1 s^-1; C1/C2: at the end of the eighth/last cycle with judgment thresholds of 1.5 s^-1; D: Time variation of thrombus proportion increase in different judgment thresholds.

DownLoad: Full-Size Img PowerPoint

4.   Discussion

There are millions of patients suffering from non-valvular AF in China, with a 20% fatality rate and a 60% disability rate ^[21]. The management of AF in China is unproductive and accompanied by a high risk of mortality due to patients* low awareness and lack of proper treatment ^[22]. Clinical studies have shown that the usage of warfarin or new anticoagulants in AF patients with a high stroke risk can significantly reduce stroke events and improve patient outcomes ^[22]. Therefore, the risk assessment of thrombosis under AF conditions is crucial in the prevention of AF-related stroke. Hemodynamic characteristics of LA can improve stroke risk stratification ^[6], but some of them are difficult to perform high-precision measurement through non-invasive means in clinical work ^[7,23]. Thus we established a thrombus model to study the AF effect on the hemodynamics in the LA by numerical simulation and to build a fundament for further application of the CFD method to study thrombosis risk prediction ^[8].

In this study, we utilized the thrombus model to predict thrombosis in LA with the consideration of non-Newtonian properties of blood for investigating the hemodynamics difference between AF and SR cases to the same LA geometry while most previous studies just simplify blood as Newtonian flow ^{[2,8,9,10,11,12]}. The thrombus model was calibrated in the BFS model before applied to the LA model because of the absence of the thrombosis experiment in a spherical tube ^[19]. The results of the benchmark study were well matched with the experiment results.

Observing the streamlines in the second peak systole, in the AF case, there is little blood flowing through the MV due to the lack of atrial kick but toward the LIPV vice versa, while in the SR case, blood can pass through the MV. It will increase blood*s impact force on the wall which is consistent with the WSS distribution. The blood flow gets highly disordered in the recirculating area which means the formed thrombus is easier to separate from the wall and then move with the blood flow. Comparing to the other areas, the velocity and WSS are small in the LAA dramatically, indicating that the LAA region is susceptible to thrombosis and may be worthy of attention in the thrombosis risk assessment as previous studies have demonstrated that low WSS distributions identify regions more prone to thrombus formation ^[24].

There are many variables in the thrombus model, the distribution of BP predicts the thrombosis area directly while RT, C, AP, and RP are intermediate variables. C has the closest tie with the BP when compared to RT, AP and RP, as the formation of C decides the initiation of BP, in addition, BP will promote the formation of C as well. The impact of C on the thrombosis was also studied in this paper. Observation of BP above the judgment threshold shows that thrombosis mainly happens in the LAA which is consistent with the clinical statistics ^[2]. Thrombosis begins in the wall because of the tissue factor exposure and interaction of cells and fibrinogen there which are enhanced with a low shear rate ^[25]. The simulation results show that the thrombus distributions are distinguished between the different cardiac rhythms. In the AF case, thrombosis is more serious than the SR case with a ratio difference close to 18%. It also means that the impact of the occurrence of high-frequency fibrillation cannot be ignored in the simulation of AF ^[10]. C is an important variable that influences the initiation of thrombosis, and its initial value depends on the underlined threshold of the time-averaged shear stress rate from previous simulation results. From the sensitivity analysis of the threshold, we have found that a higher threshold means more LAA area has non zero initial C and a faster thrombosis growth rate. It can be deduced that the threshold we have chosen possessed certain rationality, noticing that the difference between thresholds of 1.5 s^-1 and 1 s^-1 is much smaller than that of cases with thresholds of 1 s^-1 and 0.5 s^-1.

Treating the LA wall as a rigid wall may be a key limitation of our simulation in this paper. In this study, we investigated the impact of the atrial kick absence by changing the outlet condition as previous researchers did ^[7]. For AF, there are more characters that will influence the hemodynamics in the LA, such as the high-frequency fibrillation of the LA wall. To the author*s knowledge, there are two ways to rebuild the deformation of LA, one is to set the atrial wall as the fibre-reinforced hyperplastic material and consider the fluid-structure interaction effect ^[9,26], while the other way is to load the deformation laws of the LA and the LAA, which can be derived from the MR or CT slice-images of the volunteers to the geometry directly ^[8,10,12,27]. We will consider the deformation of the LA and the LAA for the thrombosis prediction by using these two ways in our future study.

In addition, the boundary conditions are not strictly patient-specific due to the absence of clinical measurement data. The aim of the present study was the application of thrombus model to simulate the thrombosis in LA rather than the accurate thrombosis prediction in special patients. Our further study will focus on the latter strategy and PC-MRI measurement data would be used as boundary conditions for further improvement of the simulation.

LAA has stronger active contractility than LA, due to its thick pectinate muscle ^[28]. For healthy people, the contraction of LAA will empty the inside blood while AF patients vice versa ^[3,28]. The neglect of the influence of the LAA contraction may lead to thrombosis in the SR patient cases. More studies will be displayed on the LAA contraction and geometry in our future research to evaluate their impacts on thrombosis.

Thrombosis is a complex process involves the conversion of many substances ^[29]. In this paper, we only consider the effects of four substances. There is no doubt that other substances also have impacts on the viscosity property of the blood flow and thrombosis which should be further studied.

5.   Conclusion

This is the first study that directly investigates the thrombosis in the LA with the discussion of the hemodynamics changes. The thrombus model used in this paper was calibrated in the BFS model and in good agreement with the experimental data. The hemodynamics of the AF cases and SR cases are compared, showing that the hemodynamics in LA is greatly affected by AF. Thrombosis observed in this study only occurs in the LAA which is consistent with clinical findings, thrombus first forms on the LAA wall then grows inward until filling up the inner space of LAA. AF can promote the degree of thrombosis. In the future, numerical simulation of the thrombosis and clinical statistics can be combined to improve the risk assessment mechanism of thrombosis and stroke.

Conflict of interest

All authors declare no conflicts of interest in this paper.