Geometrical and Hemodynamic Variables for Thrombotic Risk Strati cation in Kawasaki Disease


 Objectives

Thrombosis is a major adverse outcome for coronary artery aneurysms (CAA) in Kawasaki disease (KD). We investigated the geometric and hemodynamic abnormalities in patients with CAA and identified the risk factors for thrombosis by computational fluid dynamics (CFD) simulation.
Methods

We retrospectively studied 27 KD patients with 77 CAAs, including 20 CAAs with thrombosis in 12 patients. Patient-specific anatomic models obtained from cardiac magnetic resonance imaging (CMRI) were constructed to perform a CFD simulation. From the simulation results, we produced local hemodynamic parameters comprising of time-averaged wall shear stress (TAWSS), oscillatory shear index (OSI) and relative resident time (RRT). The CAA’s maximum diameter (Dmax) and Z-score were measured on CMRI.
Results

Giant CAAs tended to present with more severe hemodynamic abnormalities. Thrombosed CAAs exhibited lower TAWSS (1.551 ± 1.535 vs. 4.235 ± 4.640dynes/cm2, p = 0.002), higher Dmax (10.905 ± 4.125 vs. 5.791 ± 2.826mm, p = 0.008), Z-score (28.301 ± 13.558 vs. 13.045 ± 8.394, p = 0.002), OSI (0.129 ± 0.132 vs. 0.046 ± 0.080, p = 0.01), and RRT (16.780 ± 11.982s vs. 9.123 ± 11.770s, p = 0.399) than the non-thrombosed group. An ROC analysis for thrombotic risk proved that all of the five parameters had area under the ROC curves (AUC) above 0.7, with Dmax delineating the highest AUC (AUCDmax = 0.871) and a 90% sensitivity, followed by Z-score (AUCZ−score = 0.849).
Conclusions

It is reasonable to combine the geometric index with hemodynamic information to establish a severity classification for KD cases.


Introduction
Kawasaki disease (KD), also called mucocutaneous lymph node syndrome, is a systemic median vasculitis with a predominance for coronary arteries. KD is the leading cause of acquired heart disease in children. Coronary artery aneurysm (CAA) is the most common complication of KD, of which the incidence has been reduced greatly from 20-25% to 3-5% by administering aspirin and intravenous immunoglobulin (IVIG) [1,2]. The American Heart Association (AHA) stated in a scienti c statement that KD cases with associated CAAs were classi ed as at high risk of cardiac events [3]. CAAs may resolve spontaneously in 60% of cases during the subacute and chronic phase, whereas the possibility of thrombosis is the most concerning risk, which can result in myocardial infarction, ischemia-related arrythmias or sudden death. Therapeutic options for CAA in the chronic phase essentially include antiplatelet and anticoagulation therapy. The risk factors for thrombosis have not been totally clari ed and the current recommendation of combining antiplatelet and anticoagulant therapy for CAAs at high risk for cardiovascular events is based mainly on the size of CAA: that is, CAAs with a maximum diameter (Dmax) ≥ 8 mm or Z-score ≥ 10 [4].
Several research articles have suggested that hemodynamic changes in KD patients with CAAs also play a crucial role in the pathogenesis of thrombosis [5][6][7][8]. Dilation of the vascular lumen induces a local ow disturbance which results in injury of the endothelium, accelerated platelet aggregation and hypercoagulation. Computational uid dynamics (CFD) technique can provide multiple hemodynamic parameters including wall shear stress (WSS), relative resident time (RRT), and oscillatory shear index (OSI), which are almost impossible to obtain from usual imaging modalities [9]. The value of CFD in the study of KD has been proven in recent papers, but all the previous research projects were limited by small sample sizes (no more than 10 cases) [6][7][8]. In the present study, we applied patient-speci c CFD in 27 patients with CAAs. We aimed to select variables with high sensitivity and speci city for thrombosis risk strati cation.

Patient population
From May 2016 to February 2019, we enrolled 27 KD patients with CAAs admitted to our medical center. This study was approved by the hospital's Institutional Review Board and written parental consents were obtained. All the data were processed according to the relevant guidelines. A brief clinical history was collected for all the patients, including sex, age, onset of KD, clinical symptoms, treatment, and prognosis.
All the patients received echocardiography as rst screening exam and were recommended to have cardiac magnetic resonance Imaging (CMRI) when the CAAs were identi ed. The size of CAA was categorized principally according to Z-score: small (2.5 ≤ Z-score < 5), medium (5 ≤ Z-score < 10) and giant (Z-score ≥ 10 or Dmax ≥ 8mm) [4]. The thrombi inside the CAAs were also documented on CMRI images. After CMRI, all the patients had digital subtraction angiography (DSA) to identify the existence and location of CAAs and thrombosis. CMRI and DSA was conducted within 3 days after the Page 4/16 echocardiograms. The disease phases were de ned as follows: acute, day 1~14; subacute, day 15~42; and chronic, after day 43. Day 1 was de ned as the rst day of presenting clinical signs [10].

CMRI technique
CMRI was performed using a Siemens 1.5T MRI scanner (Siemens Medical Solutions, Erlangen, Germany). The entire protocol included an ECG-gated steady-state free precession (SSFP) sequence with T2 preparation to suppress the myocardium and fat tissues after a bolus of gadobutrol (0.1-0.2 mmol/kg) at a rate of 1-2 ml/s through intravenous line. A 3D volume of axial images including the whole heart was adapted for simulation [11]. Two-dimensional (2D) steady-state free precession segmented cine images in short axial and longitudinal plans with full coverage of the left and right ventricles were acquired for functional assessment. To assess brosis, late gadolinium enhancement imaging was conducted by an inversion recovery turbo fast low-angle shot sequence. Imaging was carried out under oral sedation, if necessary.

3D anatomical modeling
Patient-speci c modeling and a nite element mesh were established using the SimVascular 2018 (Stanford University, USA) to create a patient speci c model and the main process included import medical image data, create path, segmentations and 3D geometry model. After importing the created model into ICEM CFD (ANSYS, Inc. USA), the model was used to generate tetrahedral grids and set boundary layers. We also tested grid dependence, when the grid size was less than 0.2cm, the maximum velocity of the outlet of aorta tended to converge. Taking all model considerations into account, the nal girds generated were about 150,000 and element numbers 730,000.
The included arteries were the right coronary artery (RCA), the left coronary artery (LCA), the left anterior descending artery (LAD), the left circum ex artery (CCX), as well as the ascending aortic arch.

Simulation methods
Blood ow simulations were conducted with an Ansys CFX 18.0 solver. Blood was assumed as a Newtonian uid with a density of 1.06 g/ml, dynamic viscosity of 0.035 dynes/cm·s, and the vessel wall was assumed to be rigid. We used the velocity-time curve of aorta from echocardiogram to set inlet parameter. The LCA, RCA, LAD and CCX were set as ow boundary conditions according to the crosssectional area of each branch, and the pressure of aorta's outlet was measured by DSA as it was considered as the golden standard. The outlet ow used was the time-average ow and remained the same value throughout the cardiac cycle. The main method was to calculate the total ow according to the inlet velocity, and consider that the total ow of all coronary arteries accounted for 6% of the inlet ow, and then distribute the ow according to the ratio of the outlet cross section area of each coronary artery. We could not measure the time-dependent ow waveform yet. Simulations were operated for 4 cardiac cycles, after which the cycle-to-cycle pressure variations were determined as less than 0.1%. The required simulation time for each model was approximately 20 hours.

Hemodynamic and geometric variables
Simulated pressure and velocity elds were calculated and hemodynamic variables including timeaveraged wall shear stress (TAWSS), OSI and RRT were computed. RRT is the average time a parcel of uid spends in a speci c region based on the pre-computed velocity eld. It is a measure of ow stagnation, which contributes to the form of thrombus. How to calculate it is shown in the formula below [12]: Wss i is the wall shear stress of each grid, and T is the length of a cardiac cycle.
TAWSS, OSI and RRT were temporally averaged over a single cardiac cycle and spatially averaged over each CAA's surface. In addition to hemodynamic variables, we assessed anatomical variables from CMRI, such as location, Dmax and Z-score.

Statistical analysis
Continuous parameters were presented as mean ± standard deviation and inter-quartile range (IQR). Categorized data were expressed as numbers and percentages. Correlations between hemodynamic and geometric data were evaluated using the Pearson's linear correlation coe cient. We used the unpaired Student's t-test and one-way ANOVA to compare anatomical and hemodynamic variables in the different CAA subgroups. The best cutoff values were obtained from receiver operating characteristic (ROC) curves. SPSS version 26 (IBM SPSS Statistics 26, Armonk, NY, USA) was utilized and p values < 0.05 were considered to have statistically signi cant difference.

Results
Patient demographics and clinical data were summarized and illustrated in Tab. 1. All the patients presented relative symptoms like fever, rash and cervical lymphadenopathy at the onset of the disease. All the patients received aspirin and IVIG in the acute stage, and were treated by the combination of antiplatelet (usually aspirin) and anticoagulation (usually warfarin) when giant CAAs or thrombosis were found by echocardiography or CMRI. The dose of anticoagulation was adjusted to the PT-INR target range of 2.0-2.5. All the CAAs were classi ed into three groups based on the Z-score (See Tab.1). Regarding the 25 giant CAAs, 9 occurred at the LCA/LAD branching site, 1 at the LAD/CCX region, while the remaining 15 CAAs were located in the RCA. Moreover, the occurrence of the giant CAAs in the RCA was as follows; 6 were in the proximal segment, 5 in the middle segment, 2 in the remote segment, and 2 within the whole dilated RCA. Figure 1 displays the distribution of hemodynamic variables in some patients with CAAs. In coronary arteries free of aneurysms, the hemodynamic variables were in the normal range. In the dilated RCA, the disturbance was more obvious in the proximal region with abnormal hemodynamics along the dilated artery (Fig. 1A). The distribution of TAWSS, OSI and RRT were uneven, especially in the giant CAAs (Fig. 1B). As listed in Tab. 2, CAAs with a larger diameter had lower TAWSS, along with higher OSI and longer RRT.

Geometrical and hemodynamic variables for thrombotic risk strati cation
We compared geometric and hemodynamic variables in the thrombosed and non-thrombosed CAAs. All 20 thrombosed CAAs were found in 12 KD patients. For the 12 patients with 20 thrombosed CAAs, 4 cases had IVIG resistant who intended to present more severe in ammatory response and more severe coronary artery impairments. In the 20 thrombosed CAAs, 8 aneurysms occurred at the LCA/LAD branching site; 1 at the LAD/CCX region, and 11 in the RCA. 1 case presented long segmental thrombosis in the mid and remote segment of RCA (Fig. 1A, circle). 1 case had CAAs in the RCA and LAD/CCX, with thrombi in the RCA (Fig. 1C, circle); 1 case developed thrombosis in the distal aneurysm within the RCA rather than proximal (Fig. 1D, circle). Out of all the thrombosed CAAs, 4 were of median size and 16 were of giant size. In median and giant CAAs, the incidence of thrombosis was 10.26% and 64%, relatively. In contrast, no thrombosis was detected in small CAAs.
Quantitative analysis revealed signi cant discrepancies between the thrombosed and non-thrombosed groups (Fig. 3 Next, an ROC analysis was conducted based on the geometric and hemodynamic variables, then the associated cutoff value, sensitivities and speci cities were presented in Fig. 3. All of the ve parameters had area under curves (AUC) larger than 0.7, demonstrating their ability for thrombosis risk strati cation. Dmax had the highest AUC (AUC Dmax = 0.871) with 90% sensitivity, followed by Z-score (AUC Z−score = 0.849). OSI (sensitivity = 0.8) and RRT (sensitivity = 0.8) displayed an adequate sensitivity while TAWSS provided a favourable speci city (speci city = 0.8).

Discussion
KD is an idiopathic type of acute vasculitis mainly occurring in infants and children younger than 5 years old. Some potential etiologies like certain genetic markers, viral pathogens or immune responses have been reported to be related with KD [13,14]. It is most common in the Asian race compared to the non-Asian population. In Japan, Korea, and China, the incidence varies from 39 to 250 in 100,000 children < 5 years old, which appears to be on the rise [15,16]. The prevalence is much lower in non-Asian countries: 17-20 per 100 000 children < 5 years old in the USA, 26.2 in Canada, and 8.4-31 in Europe, with boys having a higher incidence than girls [16].
Coronary artery lesions are the main complications of KD, occurring in up to 40% of patients. CAA has a highest incidence of less than 5% following the administration of IVIG [17,18]. Other complications associated with KD include coronary artery dilation, thrombosis, stenosis and occlusion, leading to myocardial infarction, heart failure and sudden death. CAA usually forms 2 weeks after onset and thrombosis may be caused by injury to the endothelium and platelet activation. In the chronic phase, there is reconstruction of the vessel wall leading to progressive stenosis, thrombosis and occlusion [5,19]. Compared with atherosclerosis in adults, thrombus formation in KD is much severe. As instructed by the 2017 AHA and 2014 JCS guidelines, KD necessitates long term management till adulthood.
The size of CAA has always been used as a parameter to distinguish CAAs with higher risks for cardiovascular events. JCS and AHA guidelines of triple-therapy (anticoagulation and dual-antiplatelet therapy) are basically based on the CAA's diameter: Diameter ≥8 mm or z-score ≥10 [4,20]. Cardiac complications have been proven to be more common in giant aneurysms. In a review by Tsuda consisting of 245 patients with giant CAA, death, myocardial infarction, and surgery occurred in 15 (6%), 57 (23%), and 90 patients (37%), respectively [21]. In another study composed of 1006 patients covering 44 participating institutions, the cardiac events free rate was 100%, 96%, and 79% for small, medium, and large CAAs [22]. Over 50% of CAA regress to normal size within 5 years and the likelihood of regression seems to highly correlate with the original size of the lesion [23,24]. In agreement with previous studies, we also con rmed the value of the CAA's size. Dmax ≥ 6.850mm and Z-score ≥ 18.753 may be used as predictors of thrombotic risk.
Some speci c locations were prone to form large CAAs and thrombosis, including the LCA/LAD branching site and the RCA's proximal segment. Risk was increased when extensive dilation or multiple aneurysms were found in one coronary artery (Fig. 1A, 1B). Distal aneurysms may be more prone to thrombosis than proximal ones, even those with same sizes and similar hemodynamic dates (Fig. 1D).
Hemodynamic variables could also be useful in risk strati cation schemes for CAA. In precedent experimental studies, sites exposed to low WSS and high OSI have been demonstrated to be more susceptible to develop intimal thickening and thrombosis [25][26][27]. Local low WSS induced platelets aggregation close to the vascular wall and promoted clot formation [28]. Endothelial cells would respond to the hemodynamic changes of low WSS and high OSI to modulate intracellular signaling, stimulating structural remodeling, thrombogenesis, in ammation, and atherosclerosis [5,25]. Dibyendu discovered that aneurysms with thrombosis presented longer RRT (mean 7.8±2.8 vs. 4.0 ± 2.0) and 108% lower WSS compared to aneurysms without thrombosis [29]. Grande studied 10 KD patients with CAAs and depicted that hemodynamic parameters were sensitive to thrombosis [7].
In Fig. 2, we also noticed a large variation in TAWSS, OSI and RRT in CAAs with similar sizes. The correlation was medium between geometric and hemodynamic variables, in accordance with other papers' reports [7].
Regarding those CAAs with more obvious hemodynamic abnormalities, the risk of thrombosis was higher. For example, compared with larger CAA in LCA/LAD, the dilated RCA exhibited longer RRT along the entire artery, inducing extensive thrombosis (Fig. 1A). Another patient had two CAAs of similar dimensions located in the RCA and LCA/LAD junction, however the hemodynamic abnormalities were more severe in the RCA's CAA (lower TAWSS, higher OSI and RRT), coupled with thrombosis (Fig. 1C). In Fig. 1B and Fig. 1D, RCA had 2 CAAs and had more severe decreased TAWSS. Apart from the decreased TAWSS, the larger CAA demonstrated higher OSI and RRT which prompted the formation of thrombi (Fig. 1B). In Fig. 1D, the thrombosed CAA was in the distal segment of RCA, which had lower blood ow velocity and higher RRT. Our results suggested that hemodynamic parameters can be reliable supplements to the aneurysm size for better predicting thrombotic risk.
Nowadays, CFD studies have been coupled with coronary artery disease, including atherosclerosis, coronary artery bypass grafting surgery, and myocardial bridging [30][31][32]. Vortices, WSS, OSI, and RRT have been identi ed as key hemodynamic quantities having the potential to be conventionally used in clinical work. Other new parameters like aneurysm shape index, aneurysm spherality, and aneurysm surface area have also demonstrated the potential to describe the complexity of aneurysm geometry, but only used in single centers and still require further research to support their accuracy [7,29,33]. Based on these facts, we adopted the classical parameters for analysis in our research instead of other new indices.
The results of our study differed somewhat from the smaller study by Grande Gutierrez et al [7]. They applied uid-structure interactions (FSI) in the modeling which was not adapted in our study. FSI is a technique to combine the change of ow eld with the mechanical change of blood vascular wall.
However, in order to make the results more reliable, it is necessary to know the accurate parameter of blood vascular wall when calculating uid-solid coupling. At present, these parameters are di cult to obtain in the pediatric population. There are some papers showing that the WSS varies between FSI and rigid-wall models but the difference sometimes is small to affect the result and clinical use [34,35]. Moreover, FSI takes a long time and considerable computational power to run, so we chose to set the vessels as rigid walls in this study. In the future, with the progress of technology, FSI can be further used for the study of hemodynamics in our institution.
In our research, we employed a 3D SSFP sequence for patient-speci c modeling. Compared to CT angiography, CMRI does not involve any dose of radiologic exposure, which is particularly essential for children in development and in need of serial scanning for future follow-up. Furthermore, CMRI can also be used to observe thrombi and intimal hyperplasia of the coronary artery. The AHA and JCS both uphold the use of CMRI for myocardial injury evaluation and observation of the complete coronary tree during the follow-up of KD patients [36].

Limitations
Despite the fact that the number of patients included in this study was modest, we could still identify some limitations. First and foremost, being a retrospective study, there was a large variation in patient age, which may affect the patient characteristics and assumptions. Secondly, the images obtained from CMRI have lower image resolution for very young patients, which could potentially affect the accuracy of modelling. These differences may be improved by using the advanced CMRI scanning technique (higher eld strength and spatial resolution). Last but not least, the CFD technique still requires a long period of time for model construction and computing. Hopefully, the increasingly e cient computational methods will enhance the process of simulation to ultimately enable their routine clinical use.

Conclusion
In our study, we discovered several points about geometrical and hemodynamic variables correlating with clinical outcomes. The LCA/LAD branching site and the RCA's proximal segment were prone to the formation of giant CAAs and subsequent thrombosis. Moreover, giant CAAs tend to have more severe hemodynamic abnormalities and higher risk of thrombosis. Aneurysm size (Dmax and Z-score) has been proved to be reliable to determine the thrombotic risk. As the most frequently used hemodynamic parameters based on the CFD technique, TAWSS, OSI and RRT could also be implemented to identify patients at high risk for thrombosis. In future clinical practice, it would be reasonable to combine geometric indices with hemodynamic information for severity classi cation in KD patients.  Unpaired Student's t-test revealed signi cant discrepancies of hemodynamic (TAWSS, OSI and RRT) and geometrical variables (Dmax, Z-score) between thrombosed and non-thrombosed groups. ROC analysis showed that all of the ve parameters had AUC larger than 0.7, demonstrating their ability for thrombosis risk strati cation. Dmax = maximum diameter; TAWSS = time-averaged wall shear stress; OSI = oscillatory shear index; RRT = relative resident time; 0 = non-thrombosed group; 1 = thrombosed group; ROC = receiver operating characteristic; AUC = area under the ROC curve; sen = sensitivity; spe = speci city;