A Novel Scoring System for Rupture Risk Stratification of Intracranial Aneurysms: A Hemodynamic and Morphological Study

Objective: The aim of the present study is to investigate the potential morphological and hemodynamic risk factors related to intracranial aneurysms (IAs) rupture and establish a system to stratify the risk of IAs rupture to help the clinical decision-making. Methods: Patients admitted to our hospital for single-IAs were selected from January 2012 and January 2018. A propensity score matching was conducted to match patients. The morphological parameters were obtained from high solution CTA images, and the hemodynamic parameters were obtained in accordance with the outcomes of computational fluid dynamics (CFDs) simulation. Differences in the morphologic and hemodynamic parameters were compared. The significant parameters were selected to establish a novel scoring system (Intracranial Aneurysm Rupture Score, IARS). The comparison was drawn between the discriminating accuracy of IARS and the Rupture Resemblance Score (RRS) system to verify the value of IARS. Then, a group of patients with unruptured IAs was stratified into the high risk and low risk groups by IARS and RRS system separately and was followed up for 18–27 months to verify the value of IARS. The outcome of different stratifications was compared. Results: The matching process yielded 167 patients in each group. Differences of statistical significance were found in aneurysm length (p = 0.001), perpendicular height (H) (p < 0.001), aspect ratio (AR) (p < 0.001), size ratio (SR) (p < 0.001), deviated angle (DA) (p < 0.001), normalized average wall shear stress (NWSSa) (p < 0.001), wall shear stress gradient (WSSG) (p < 0.001), low shear area ratio (LSAR) (p = 0.01), and oscillatory shear index (OSI) (p = 0.01). Logistic regression analysis further demonstrated that SR, DA, NWSSa, LSAR, and OSI were the independent risk factors of IAs rupture. SR, DA, LSAR, and OSI were finally selected to establish the IARS. Our present IARS showed a higher discriminating value (AUC 0.81 vs. 0.77) in comparison with the RRS (SR, NWSSa, and OSI). After follow-up, seven patients were subject to IAs rupture. 5/26 in high risk group stratified by IARS, yet 7/57 in high risk group stratified by RRS. The accuracy of IARS was further verified (19.2% vs. 12.3%, AUC for the IARS and the RRS was 0.723 and 0.673, respectively). Conclusion: SR, DA, NWSSa, LSAR, and OSI were considered the independent risk factors of IAs rupture. Our novel IARS showed higher accuracy in discriminating IA rupture in comparison with RRS.


INTRODUCTION
Intracranial aneurysms threaten as much as 2-3.5% of the population worldwide and are the leading cause of non-traumatic SAH (Müller et al., 2013). Thus, IAs could cause a devastating outcome with high mortality (50%) (Cross et al., 2003). Accurate discrimination of IAs prone to rupture and optimal treatment decision-making now still pose a great challenge to the clinical practicers (Thompson et al., 2015). The CFDs analysis has been demonstrated be a useful method to assist clinical strategy of cerebral vascular diseases in recent years. Several hemodynamic and morphological parameters have been demonstrated to be the independent risk factors for IAs rupture (Xiang et al., 2011Meng et al., 2012;Paliwal et al., 2017). Until now, several stratification models have tried to establish risk of IAs rupture based on morphological and hemodynamic models (Dhar et al., 2008;Xiang et al., 2011). An effective system called rupture resemblance score (RRS) system [SR, normalized wall shear stress (NWSS), and OSI were included] was proposed recently by Meng et al. through the combination of the hemodynamic and morphological features of ruptured IAs (Xiang et al., 2011(Xiang et al., , 2015(Xiang et al., , 2016). Yet results from recent studies suggested that WSS might be strongly collinearity with the morphological characteristics (including AR). Besides, some other parameters including LSAR and WSSG were also correlated with IAs rupture closely (Fan et al., 2015;Zhang et al., 2016a,b;Levitt et al., 2017;Sano et al., 2017;Skodvin et al., 2017a;Varble et al., 2017). Given these, we assumed that higher discriminating accuracy could be achieved in discriminating ruptured from unruptured IAs.
The aim of this study is to identify the independent risk factors related to IAs rupture and establish a more comprehensive scoring system for rupture risk stratification of IAs using hemodynamic-morphological analysis, and to verify the effectiveness of our novel scoring system by comparing the value with previous scoring system.

Patient Selected and Inclusion/Exclusion Standards
This study retrospectively reviewed the patients admitted into our hospital for ruptured IAs between January 2012 and January Abbreviations: AR, aspect ratio; AUC, area under the curve; CFDs, computational fluid dynamics; CTA, compute tomography angiography; DA, deviated angle; IAs, intracranial aneurysms; LSAR, low shear area ratio; NPa, normalized average pressure; NPm, normalized maximal pressure; NWSSa, normalized average wall shear stress; NWSSm, normalized maximal wall shear stress; OSI, oscillatory shear index; ROC, receiver operating characteristic; RRT, relative resident time; SAH, subarachnoid hemorrhage; SR, size ratio; WSS, wall shear stress; WSSG, wall shear stress gradient.
2018. The patients were selected in accordance with the standards as follows: Inclusion criteria (1) just had one IA, (2) had complete clinical records and follow-up data, (3) intracranial bleeding was identified by CT scan within 24 h after bleeding, and (4) nearby CTA was conducted 20-120 days before bleeding. Exclusion criteria: (1) patients were related to other intracranial tumor, angiostenosis and angio-malformation, including arteriovenous malformation, cavernous malformation, etc. and (2) CTA data was not suitable for morphological analysis or hemodynamic analysis.

Propensity Score Matching
To exclude known clinical risk factors and balance the baseline, a propensity score matching (PSM) was conducted to match patients using STATA (12SE, Stata corporation, United States) based on database of patients with unruptured IAs in our institution. The propensity scores were calculated by using a logistic regression model consisting of the input variables: gender, age, hypertension, atherosclerosis, ever smoker, and IAs' location. The matching rate was 1:1 for ruptured IAs to unruptured IAs. The final matched group was confirmed as appropriately matching using χ 2 testing to validate equivalence of individual variables between each group.

Vascular Modeling
The dicom data of last CTA, conducted before bleeding, were collected from the high solution CTA work station (Siemens, Berlin, Germany) and converted into slice dicom data (about 0.5 mm per slice). The dicom data were introduced into Mimics 17.0 (Mimics Research 17.0, Materialize, Belgium) and reconstructed for further study.

Radiological Measuring and Morphological Parameter
Radiological measuring was performed by two experience neurosurgeons (PJ and JW) using high solution CTA. Maximum length (L), maximum diameter of body (D), diameter of neck (d), perpendicular height (H), diameter of parent artery (P), and volume are measured from CTA as shown in Figures 1a-e. H is the maximum perpendicular distance of the dome from the neck plane. L is the maximum distance of the dome from the neck plane. d is the average diameter of the neck. p is average diameter of the parent artery. D is the maximum diameter of the body. Those parameters were measured twice, and the average was taken. AR, SR, undulation index (UI), ellipticity index (EI), and nonsphericity index (NSI) were calculated according to previous study (Dhar et al., 2008). FIGURE 1 | Morphological parameters were measured from CTA image. D is the maximum diameter of the body (a). L is the maximum distance of the dome from the neck plane (b). d is the average diameter of the neck (c). H is the maximum perpendicular distance of the dome from the neck plane (d). The method of measurement is shown in (e). p is average diameter of the parent artery. Normal vector was combined. The deviated angle (DA), which was between co-velocity and normal vector, was measured (f).

CFD Simulating and Hemodynamic Parameter
To create 4-5 million units of finite tetrahedral and prism elements, each aneurysm model was meshed by STAR-CCM (STAR-CCM+ 12, Siemens, Germany). Then, the models were introduced into STAR-CCM fluid workstation (STAR-CCM+ 12, Siemens, Germany). In line with previous study (Xiang et al., 2011), incompressible Navier-Stokes equation served as the solver under pulsatile blood condition. When the patient enters the operating room, we use ultrasound to record the pulsatile waveform of the internal carotid artery. The pulsatile waveform was obtained by using transcranial Doppler ultrasound device on a representative patient with its magnitude scaled to the desired mean flow rate and was plotted to pulsatile curve. Blood was assumed as a Newtonian fluid with density ρ = 1056 kg/m 3 and viscosity µ = 0.0035 Poise. Pulsatile curve served as the velocity inlet boundary condition, and free boundary condition was implemented at outlet. Four pulsatile cycles were simulated. The last cycle was yielded for the further study.
For each model, the normal vector of the vascular wall (the direction of flow was positive direction) was first reconstructed. The diverse-directed velocities were combined into a co-velocity as shown in Figure 1f. The angle between the normal vector and co-velocity was measured and defined as DA. Pressure maximum (Pm), pressure average (Pa), WSS maximum (WSSm), and WSS average (WSSa) were obtained from IA region, and parent Pa (pPa) and parent WSS average (pWSSa) were obtained from parent artery region. The WSS average and pressure average at peak of the systolic phase were applied for the further study. Low WSS was defined as less than 10% of WSS of parent artery (Xiang et al., 2011). The NPa, NWSSa, NPm, and NWSSm were calculated in line with the equations (1, 2, 3, and 4), respectively. OSI and RRT were calculated according to previous study (Dhar et al., 2008). OSI and RRT were the average over the dome area.

Patients Follow-Up and Validation
To verify the accuracy of our stratification system, a group of other patients with unruptured IAs was followed up, whose aneurysms were found unintentionally. The patients who visited in the outpatient from January 2016 to October 2016 were enrolled. Those patients were followed up for 18-27 months by clinical visiting. To identify the change of IAs, CTA was conducted for every 4-6 months for patients with aneurysm under high rupture risk, and CTA was conducted for every 5-6 months for patients with aneurysm under low rupture risk. Another PSM executed. Propensity scores were calculated for ruptured and unruptured IAs. A group of patients were selected in line with gender, age, hypertension, atherosclerosis, and ever smoker, and then stratified by IARS and RRS system, respectively. We matched the ruptured IAs to the unruptured IAs using a 1:9.

Statistical Analysis
Measurement variables were compared by performing the chi-square test or Fisher's exact test. Continuous variables were firstly assessed visually by performing the P-P plots and the Shapiro-Wilk test and then compared by performing the independent samples t-test. The ROC curve analysis was conducted for significant parameters on univariate analyses. Highest Youden index was employed to find the appropriate FIGURE 2 | The selecting flow chart was presented here. A total of 941 patients were enrolled into our study. 167 patients with ruptured IA met our standards in our series. Another 167 patients with unruptured IA were then selected using the propensity score matching. The matching process yielded 167 patients in each group.
cut-off value. The results are presented in 95% confidence intervals. Intercorrelations between parameters were examined using Person correlation. To clarify the independent risk factors related to the rupture statues, multivariate regression analyses were conducted. ROC analyses were performed for the followup patients based on the IARS and the RRS, and the AUC were compared. The ruptured rate in high risk group and the unruptured rate in low risk group were also compared. A pvalue of <0.05 was assumed to be of statistical significance. The result was expressed in 95% confidence interval. All statistical analysis was performed using SPSS 22.0 (IBM, New York, United States).

Patient and Aneurysm Characteristics
A total of 941 patients were enrolled into our study. 167 patients with ruptured IA met our standards in our series. Another 167 patients with unruptured IA were then selected using the PSM (selecting flow chart was given as Figure 2). Finally, a casecontrol group, which included 334 patients (167 patients with rupture IA and 167 patients with unruptured IA) was created. The female rate was 55.1% in patients with rupture aneurysms and 55.7% for patients with unruptured aneurysms, respectively. Mean age was 52.3 for ruptured group and 53.9 for unruptured group, respectively. The difference in hypertension history, atherosclerosis history, and smoking history were equal between ruptured and unruptured IAs. No significant difference was found in family history of non-traumatic subarachnoid history (p = 0.247). The demographic information was summarized in Table 1.
Large area with low WSS could be found in the dome of aneurysm (b.6).

Further Validation
Two hundred and twenty one patients visited in the outpatient for IA in total. At the end of following up, 7 patients (the rupture rate was 1.6% per year in our series) were subject to IAs rupture. The matching process yielded 70 appropriate patients (included 7 ruptured IAs and 63 unruptured IAs, the hemodynamic-morphological characteristics of the ruptured IAs were given in Table 5). Demography, morphological, and hemodynamic information for each patient were given in Supplementary Table 1. No statistically significant difference was found in gender, age, hypertension, atherosclerosis, now or ever smoker, and aneurysm location. Stratified by IARS, 5/26 patients in high risk group (19.2%) undergoing SAH or intracranial hemorrhage, caused by IAs, and only 2/44 patients in low risk group (4.55%) were subject to that during the followup period. Yet stratified by RRS, 7/57 patients in high risk group (12.3%) were subject to SAH or intracranial hemorrhage, and 0/13 patients in low risk group (0.0%) were subject to that (a hemodynamic-morphological analysis for a ruptured IA was shown in Figure 5). The discriminating accuracy is significantly different between two score system (19.2% vs. 12.3%). As revealed by ROC analysis, the IARS has higher AUC (AUC for the IARS and the RRS was 0.723 and 0.673, respectively, Figures 4d-f). The result of ROC and rupture rate was summarized in Table 5.

DISCUSSION
Some clinical trials have tried to establish a method to assess rupture risk of IAs, such as the International Study of Unruptured Intracranial Aneurysms study and the PHASES study (Greving et al., 2014). Yet models established in those studies included only clinical factors, which were difficult to be explained by the internal mechanisms of IAs rupture. Recent studies suggested that the hemodynamic and morphological characteristics of aneurysm were also important to assess rupture risk. In this study, by using a PSM matched IA data base of 334 aneurysms from our institution, we attempted to find out hemodynamic and morphological characteristics of ruptured IAs and establish a useful score system to help clinical work.
Following the previous studies, we balanced the baseline and excluded the known risk factors by using the PSM. In the International Study of Unruptured Intracranial Aneurysms study, author reported that different location had different rupture rate (Wiebers, 2006). Previous Varble et al. (2017) found that internal carotid aneurysms are the least rupture-prone. Thus, location is an important risk factor to assess rupture risk. By using pool analysis, the PHASES study demonstrated population, hypertension, age, SAH history, and site of aneurysm as independent risk factors for aneurysm rupture (Greving et al., 2014). Here, we balanced those risk factors between ruptured and unruptured IAs by using PSM. This statistical method could exclude the confounding factors.
SR is relative to the parent vessel, has been independently associated with rupture risk in many studies. High SR, combined with abnormal hemodynamic conditions, may result in an aneurysm rupture. Previous CFD studies have found that with increasing SR, aneurysm rupture risk increased. However, previous Meng et al. reported that a SR > 1.75 could increase the rupture risk (Xiang et al., 2011). In another study, Dhar et al. (2008) reported that a SR > 2.05 was FIGURE 4 | The ROC curves of key parameters (DA, SR, LSAR, OSI, and NWSS) were presented here (a,b). The cut-off values for SR, DA, NWSSa, LSAR, OSI, and RRT were 2.3, 35, 0.24, 0.30, 0.008, and 5.3, respectively. The ROC curves of the IARS-SR, IARS-NWSS, and RRS suggested that the AUC of IARS-SR (AUC 0.81) was larger than IARS-NWSS (AUC 0.79) and RRS (AUC 0.77), which suggested that the IARS-SR had higher accuracy to discriminate the ruputre statue (c). For the follow-up group, the distributions of the ruptured aneurysm and unruptured aneurysm in each stratification system were presented (d,e). The IARS has higher AUC (AUC for the IARS and the RRS was 0.723 and 0.673, respectively. a risk factor for IA rupture. Based on our data, we found that cut-off value between ruptured and unruptured IAs was 2.3, which was larger than previous study. The finding is not surprising because previous many studies calculated this parameter based on the image data after rupture. However, a study analyzed the morphological change before and after aneurysm rupture and gave the conclusion that postrupture morphology should not be considered to assess the rupture risk (Skodvin et al., 2017b). It is essential to take this factor into consideration. Low WSS has been previously relative to IA rupture. Chien et al. (2009) found that low WSS was associated IA rupture in poster circulation aneurysms and AcomA aneurysms. Miura et al. (2013) demonstrated low WSS as an independent risk factor for IA rupture. By a PHASES-based study, Varble et al. (2017) found the IAs other than the ICA are subjected to a low WSS hemodynamic condition that may lead to IA rupture. Thus, WSS, the frictional force between blood and aneurysm wall, could use as a parameter to assess the rupture risk. The other important parameter is LSAR. The friction between blood and endothelial cells is essential for normal arterial wall proliferation. Low or stagnant flow can result in an inflammatory response in the vascular wall . This phenomenon was found in the ruptured aneurysm wall as well as apoptosis and degeneration of the matrix (Cebral et al., 2017). Thus, larger low WSS area suggested severe injury for aneurysm wall. In 36 aneurysms, previous study found LSAR was associated with aneurysm rupture by using case-control study (Varble et al., 2017). Our data reflect a similar result, that aneurysm with low WSS and high LSAR have a high risk of incurring an aneurysm rupture. However, there is a difference comparing to result from previous study. Previous Meng et al. has reported that NWSSa in rupture aneurysm was 0.33 ± 0.28 and suggested that a cut-off FIGURE 5 | A female patient with middle cerebral artery aneurysm was followed up by clinical visiting. The patient had hypertension history. Based on the first CTA, SR, DA, LSAR, and OSI were 4.6, 42 • , 0.74 and 0.22, respectively (c.1,3,4,5,6). The aneurysm was scored by the IARS as 4, which was stratified as high-risk group. 3.5 months later, the aneurysm was rupture. SAH and hemotoma were confirmed by CT (c.2).
value of NWSSa was 0.39, which were quite higher than our data (Xiang et al., 2011). In another study, NWSSa was 0.327 ± 0.181 (Qiu et al., 2017). Our data suggested a similar result, that NWSSa in ruptured aneurysm is 0.21 ± 0.14 and the cut-off value of NWSSa is 0.24. Meanwhile Meng et al. also reported LSAR in ruptured aneurysm was 0.38 ± 0.31 and cut-off value was 0.21 (Xiang et al., 2011). However, in a case-control study, the author found a LSAR in ruptured aneurysm was 0.15 ± 0.16 (Skodvin et al., 2017a). In our series, LSAR in ruptured aneurysm was 0.39 ± 0.29 which was similar to Meng et al., but we suggested a cut-off value of LSAR was 0.30 which was much larger than before. For this, we thought a similar explanation of SR may apply to NWSSa and LSAR. However, the real threshold may need a further exploration. OSI was also considered in many previous studies. High OSI is known to upregulate endothelial surface adhesion molecules, which can cause dysfunction of nitrous oxide and inflammatory cell infiltration (Malek et al., 1999).
Many studies demonstrated it as an independent risk factor for aneurysm rupture (Xiang et al., 2011;Qiu et al., 2017). Previous study found that OSI in ruptured aneurysms was 0.016 ± 0.031 (Xiang et al., 2011). Our data have similar result, and it is useful to assess the rupture risk. A novel parameter called DA was further identified, which was defined as the angle between the normal vector and co-velocity. The deviation of blood stream from the vascular center and the possibility of impact effect of vascular wall could be reflected by DA. There have been some studies that have confirmed the role of the impact effect in occurrence and rupture of aneurysms (Omodaka et al., 2012;Zhang et al., 2014;Cebral et al., 2015a). Cebral et al. (2017) reported flow impaction could cause increased stress in the aneurysmal wall where the flow is typically lower, which may lead to the endothelia injury. Those areas are usually presented as low WSS. Malek et al. (1999) found biological responses of the endothelium to low WSS. Also, several studies have confirmed that endothelial cells are damaged by blood flow with low WSS and high OSI, which activates the inflammation or results in wall remolding and atherosclerosis (Malek et al., 1999;Hashimoto et al., 2006). Thus, the deviated blood stream, which could reflect by DA, could damage the vascular endothelia and lead to catastrophic outcome. By considering this, the interaction between blood stream and vascular wall can be understood more thoroughly.
Yet results from many recent studies found that there were many other risk factors contribute significantly to IAs rupture, which were not involved in RRS system. The main reason why RRS did not include these risk factors was that the aneurysm studied was the ruptured aneurysm. After the rupture, the morphology and hemodynamic characteristics of the aneurysm will change dramatically (Skodvin et al., 2017b). Thus, although RRS can already effectively discriminate ruptured and unruptured aneurysms, the accuracy of it could be further increased. As described above, those significant parameters were not totally included in our model because NWSSa suggested a strong collinearity with SR, which should be considered. A study found that NWSSa and OSI are easy to resolve and have good convergence (Evju et al., 2017). Thus, they are suitable to use as discriminating parameters. Accordingly, we established two models in accordance with SR and NWSS, respectively, and confirmed that the IARS-SR has the highest accuracy. Till now, many studies supposed the idea that combining hemodynamic and morphological characteristics is closely related to IAs rupture (Omodaka et al., 2012;Zhang et al., 2014;Cebral et al., 2015a,b). Thus, considering both hemodynamic and morphological factors can thoroughly understand the mechanism of IAs rupture. That is the reason why RRS can effectively discriminate IAs rupture. Yet recent many studies propose other parameters including LSAR, which RSS was not involved, could also suggest the abnormal characteristics (Zhang et al., 2016a,b). In line with recent studies, SR, NWSSa, and OSI, in which RSS is considered, cannot just reflect the possible area of vascular injury and impact effect to vascular wall (Castro, 2013;Chen et al., 2016;Paliwal et al., 2017). In this study, by following up a group of unruptured patients, the unruptured ratio in low risk group was 100% stratified by the RRS but 95.5% stratified by the IARS. However, the ruptured rate in high risk group was 19.2% stratified by the IARS but 12.3% stratified by the RRS. Further ROC analysis was also suggested that the IARS had higher accuracy in discriminating rupture IAs here (AUC for the IARS and the RRS was 0.723 and 0.673, respectively). Meanwhile, if a tool can effectively help clinical work, it should be easy to understand and simple operation. Our study developed a score system to simplify and visualize the morphological-hemodynamic analysis of aneurysm. The proposed model is meant to help physicians who are contemplating timely intervention in face of IA patients.

LIMITATIONS
There are some limitations in our study. First and foremost, as this study is a retrospective study, small sample and short-time follow-up, the conclusion is limited and needs further study to provide more evidences to confirm value of the IARS. Second, although we used the image before rupture to avoid the effect from the event of the rupture, which might not completely reflect the rupture condition. However, the considered statue was nearly close to rupture condition. It may need further prospective study to explore that. Third, we just balanced the risk factors based on previous studies here. However, there may be potential risk factors which may affect the hemodynamic condition in aneurysm. Thus, it was difficult to avoid selection bias because of high selection (we selected 167 from 941 patients with ruptured aneurysm). Fourth, all morphological parameters were calculated from 3D CTA which is suitable for clinical visiting and can exclude the effect from thrombus, but do not give essentially image like 3D angiography. Meanwhile, all hemodynamic models were established from 3D CTA that may affect the accuracy of our simulation. Correlation with angiogram can improve the result, but we did not routinely perform angiogram in clinical visiting. Fifth, our threshold of some parameters was different with previous study. However, all aneurysms in our series came from a single center. The real threshold may require pooling of data from multiple centers. Sixth, some clinical characteristics were not considered, such as multi-IAs and IAs with vascular stenosis.
There are also several limitations related to our CFD analysis. We used the similar modeling and simulating method from Meng et al. The limitations that Meng have given will also exist in our study, such as the traction-free boundary condition at the outlet and definition of low WSS. In our study, differently from Meng, we used the velocity inlet boundary from a representative patient. However, the velocity boundary does not match with the CTA which was scanned before aneurysm rupture. Thus, the pulsatile wave may not be a good representation of the ruptured condition. Considering the CFD is sensitive to the velocity condition, it may affect our result. In the other hand, the effects from magnitude of wave form and so on will also limit our result. In the future study, we may invite reliable ultrasonic expert to join our study. The parameters in the IARS (SR, NWSSa, and OSI) are relatively easy to resolve. However, strict convergence in traditional sense is difficult to reach. Our effort was done to reach the better result. All those questions will be solved in the further study. However, multiple research groups are trying to use CFD analysis to assess aneurysm rupture risk. We believe that CFD may help physicians to make decision-making in the future. CONCLUSION SR, DA, NWSS, LSAR, and OSI were independent risk factors related to IAs rupture. In comparison with RRS system (including SR, NWSS, and OSI), our present IARS system (including SR, DA, LSAR, and OSI) had higher accuracy in discriminating the rupture status of IAs.

ETHICS STATEMENT
The study was approved by the institutional review board. Written informed consents were obtained from all participants or their legally authorized representatives and privacy of patients was effectively protected.

AUTHOR CONTRIBUTIONS
SW and YC were in charge of supervising the whole study. PJ and QL were responsible for the design of the article and drafting. QL was responsible for revising, responding to reviewer and statistical analysis. PJ and JW were also responsible for statistical analysis and confirmed the final article. ML, ZL, SY, and RG were responsible for collecting data. BG and YC provided a reliable technical support.