The discovery dataset included a total of 760 patients recruited from Fujian Medical University Union Hospital (Fujian, China) from April 2017 to January 2019. The inclusion criteria were as follows: patients diagnosed with DM (American Diabetes Association, 2019), older than 18 years, first visited our clinic, and underwent a baPWV test. The exclusion criteria were as follows: patients with an ankle-brachial index (ABI) less than 0.9 (Ato, 2018); diagnosis of severe arrhythmia, pulmonary, renal, rheumatic diseases, heart valve disease, aortopathy, and myocarditis; and antibiotic and probiotic usage in the past 3 months. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). The study was approved by the Medical Faculty of Fujian Medical University Union Hospital Ethics Committee (NO.: 2020KY031) and individual consent for this retrospective analysis was waived.

The validation dataset from an existing study from the Dryad Digital Repository¹ (Fukuda et al., 2014b) was used to further evaluate the performance of the predictive model. A total of 912 patients from Murakami Memorial Hospital in Japan from March 2004 to December 2012 were recruited in this study. Detailed information about this cohort is described in the original study publication (Fukuda et al., 2014a).

The automatic artery stiffness tester BP203RPE-II (VP-1000; Omron, Kyoto, Japan) was used to measure baPWV, blood pressure, and ABI. The patients were divided into two groups: baPWV ≥ 1,400 cm/s as the elevated arterial stiffness (EAS) group and baPWV < 1,400 cm/s as the non-EAS group (non-EAS).

A standardized questionnaire regarding demographic characteristics, blood test indicators, arterial elasticity indicators, hemodynamic parameters, echocardiographic parameters, and carotid artery ultrasound parameters was administered by the same trained team of interviewers. The body mass index (BMI) was based on the height and weight: BMI (kg/m²) = weight (kg)/height² (m²). The estimated glomerular filtration rate (eGFR) was calculated according to the CKD-EPI formula (Biljak et al., 2017). The ascending aortic diameter (AO) and other parameters were measured via echocardiography according to the American Society of Echocardiography guidelines (Mitchell et al., 2019). The internal diameter of the common carotid artery and other parameters were measured by carotid vascular ultrasound. Alcohol consumption was categorized into two groups: no alcohol and > 30 g/week beginning at least 1 year after drinking (Yang et al., 2010). The smoking status was classified into two groups: non-smoker and current smoker (continuously smoking one or more cigarettes a day for at least 6 months) (Qian et al., 2010). Postmenopausal state was defined as amenorrhea for 12 consecutive months, excluding other pathological or physiological causes (Hamaguchi et al., 2012). Coronary heart disease was diagnosed according to the European Society of Cardiology (ESC) guidelines (Taylor, 2013). Hypertension was diagnosed according to the ESC guidelines (Williams et al., 2018). Diabetes was diagnosed according to the American Diabetes Association guidelines (American Diabetes Association, 2019). Carotid artery plaque was diagnosed according to the European Mannheim consensus (Touboul et al., 2012). Carotid intima-media thickness (CIMT) was evaluated using ultrasound, and CIMT > 1 mm signified the thickening of carotid intima (Guo et al., 2020).

Least absolute shrinkage and selection operator (LASSO) is a compression estimation algorithm, which adds a penalty parameter to least squares regression to compress the estimated variables, thereby improving the prediction accuracy and interpretation of a model (Cecelja et al., 2020; Zhang K. et al., 2021). Thus, we used LASSO to select candidate variables. In addition, a support vector machine-recursive feature elimination (SVM-RFE) analysis was performed for variable selection (Wang L. et al., 2016). Finally, we combined variables from either the LASSO or SVM-RFE algorithm and then selected variables that are easily available in clinical practice for subsequent model development. LASSO and SVM-RFE were performed using glmnet (version 3.0−2) and e1071 (version 1.7−3) R packages, respectively.

Accurate prediction of EAS is important for clinical treatment decisions and can avoid excessive medical treatment caused by false-positive prediction. Thus, in this study, we aimed to achieve a simple and high-accuracy predictive model. Machine learning algorithms, including decision tree (DT), support vector machine (SVM), random forest (RF), and gradient boosting (GB), were used to construct the model, and then their performances were compared to determine the best model.

Decision tree is a tree structure model that consists of a root node and several internal nodes and leaf nodes. The root node contains all samples, each internal node represents a decision point corresponding to a single attribute, and each leaf represents a single class label (Podgorelec et al., 2002; Krzywinski and Altman, 2017). The sample was classified based on the structure of the DT model level by level. Given that DT has a high degree of transparency and is not affected by data scaling, we first used DT to construct the model. Although DT can provide a complete decision-making process for clinical problems, it often suffers from overfitting, which increases the complexity of the model and may result in poor performance on generalization. Thus, the second algorithm, SVM, with excellent generalization capability, is also used to construct an optimal classification hyperplane in an N-dimensional feature space (N: the number of features) to separate the two classes of data points. SVM is a supervised learning method based on the principle of structural risk minimization for classification prediction and non-linear regression (Nedaie and Najafi, 2018). Finally, ensemble learning methods including RF and GB, which aim to reduce the variance in models and further improve the accuracy of predictions by combining multiple models instead of using a single model, were used to develop the models. The RF model is based on the DT method, which parallelly combines a large number of DTs using bootstrap resampling to generate a model with a lower variance and better generalization than a single DT (Friedman, 2001). GB goes one step further, improving performance over iterations rather than averaging predictive results from all DTs in an RF (Friedman, 2001). GB generates a new DT based on previous DTs by reducing prediction errors when blended with previous ones.

To obtain optimal hyperparameters, the area under the receiver operating characteristic curve (AUROC) was evaluated based on a 10-fold cross-validation with different parameters in the discovery cohort. We tuned the complexity parameter for the DT model, the ntree and mtry parameters for the RF model, and multiple parameters (interaction.depth, n.tree, shrinkage, and n.minobsinnode) for the GB model. The SVM, DT, RF, and GB models were constructed using svm (version 1.4.0), rpart (version 4.1−15), randomForest (version 4.6−14), and gbm (version 2.1.5) R packages, respectively. p < 0.05 indicates a statistically significant difference.

The discovery data were randomly split into two groups 100 times: training data (70%) and testing data (30%). Each time, we first developed the four different machine learning models on the training data based on the previous tuning parameters. We then calculated the AUROC and area under the precision-recall curve (AUPRC) of the four machine learning algorithms on the testing data. Finally, we compared the values of AUROC and AUPRC from the four models to determine which model performed best. After selecting the best-performing model as the final model, we constructed the GB model using the full discovery data. Youden’s index was calculated to determine the best cutoff value of the GB model. To further validate the classification capacity of the GB model, we applied our trained GB model on an independent validation cohort. The ROC and PRC were analyzed using pROC (version 1.16.2) and PRROC (version 1.3.1) R packages, respectively.

To develop a web application for EAS assessment that is applicable in daily clinical practice, we designed a web-based tool, an EAS predictor, allowing access to our final trained model. Specifically, we used front-end development technologies, Node.js (v12.14.0), React (v16.13.1), and Ant Design (v4.5.4), to simplify the development process. RestRserve (v0.3.0) back-end development technology was used to load the final trained model. Data required for prediction were received by the model using the TCP/IP method, and then the predictive result was returned. This web tool is hosted on our server, which is freely accessible via http://vascularagingpredictor.top/.

The Kolmogorov–Smirnov (K–S) test was used to assess the normality of data. Continuous data with a normal distribution are presented as mean values ± standard deviation (SD), whereas continuous data with non-normal distribution are presented as median values (quartile). Student’s t-test was used for the comparison of continuous data following a normal distribution, and Mann–Whitney U test was used for the comparison of data with a non-normal distribution. Categorical data are presented as frequency (percentage), and comparisons between two groups were performed using the χ2 test or Fisher’s exact test (if theoretical frequency T < 5). The above statistical analysis was performed using R software 3.6.2².

This study enrolled 760 subjects with a mean age of 56 ± 12 years (60.1% male, 39.9% female). Based on the dividing value of 1,400 cm/s of baPWV, the subjects were divided into two groups: non-EAS and EAS. The complete data of patients include demographic information, chemistry indicators, diseases, hemodynamic parameters, and echocardiographic and carotid artery ultrasound parameters in each group (Table 1 and Supplementary Tables 1–3). A total of 230 patients with a mean age ± SD of 48 ± 13 (66.5% male, 33.5% female) are in the non-EAS group, whereas a total of 530 patients with age of 60 ± 9 (57.4% male, 42.6% female) are in the EAS group. Significant differences in baPWV values were observed between the non-EAS and EAS groups.

Two different algorithms, LASSO and SVM-RFE, were applied to select the most significant features for classifying individuals with normal (<1,400 cm/s) or abnormally elevated baPWV (≥1,400 cm/s). First, all features (a total of 99 variables) were included in the LASSO regression analysis and narrowed down to 15 features with non-zero β coefficients in the LASSO regression model (Figures 1A,B and Supplementary Table 4). Second, SVM-RFE was analyzed to select the top 15 important features (Supplementary Table 5). We combined features from either the LASSO or SVM-RFE algorithm, and then we further selected four variables (age, SBP, DBP, and BMI) that are easily available in clinical practice for subsequent model construction.

Before the model construction using the full discovery dataset, we first tuned the parameters of the model based on a 10-fold cross-validation. We found that when the complexity parameter of the DT model was set as 0, the model achieved the highest AUC value, whereas when the RF ntree = 1,000 and mtry = 2, the model achieved the best performance. For the GB model, the best performance was obtained when interaction.depth = 2, n.trees = 400, shrinkage = 0.02, and n.minobsinnode = 5 (Supplementary Figures 1A–C).

Next, we randomly divided the discovery dataset into two groups 100 times: training data (70%) and testing data (30%). Each time, four cutting-edge machine learning algorithms with the optimized parameters were used to develop models on the training data. Based on the obtained models, we used the testing data to assess the probability of EAS of the testing population, and ROC and PRC analyses were performed, followed by calculating the AUC values on the testing data. We compared the models and observed that DT was associated with significantly lower AUROC and AUPRC values, whereas the GB approach has higher AUROC and AUPRC values (Figure 2). Moreover, the two ensemble learning algorithms (RF and GB), especially GB, have lower variances compared to DT and SVM (Figure 2). Altogether, the GB algorithm outperformed the other machine learning algorithms in terms of the classification capacity of EAS.

Based on the GB algorithm, we finalized our EAS predictive model by training the GB model on the full discovery data with optimized parameters and calculated GB-based risk scores. In addition, we applied the GB model to the validation data from an independent Japanese cohort. Each predictor and other demographic information for the discovery and validation data are shown in Table 2. The AUROC and AUPRC values were assessed in both cohorts. The results showed high AUROC values of 0.928 and 0.821 and AUPRC values of 0.964 and 0.798 in the discovery and validation datasets, respectively, for the classification between non-EAS and EAS (Figures 3A,B).

To determine the best cutoff value of the GB model, Youden’s indexes were calculated in both cohorts. The cutoff value (0.75) built on the discovery cohort was higher than that (0.46) of the validation cohort (Figure 4), suggesting that the cutoff value derived from one cohort might not be ideal for other cohorts from different countries. Given that a cutoff value of 0.46 resulted in a better classification performance in both cohorts relative to a cutoff value of 0.75, which led to more false negative findings due to low sensitivity (0.677; Supplementary Table 6), we, therefore, selected 0.46 as a cutoff value for GB scores. The specificity/sensitivity in the discovery and validation cohorts at this cutoff value for the GB scoring system were 0.813/0.875 and 0.761/0.738, respectively (Supplementary Table 6).

To facilitate further study and use of this GB model for EAS prediction, we built a free and user-friendly online web-based tool (elevated arterial stiffness predictor:³). Figure 5 shows the user interface (UI) of the web tool. To use this web application, one only needs to input values for age, SBP, DBP, weight, and height, followed by clicking the “Predictor” button. Then, the UI will display the BMI and GB risk score value for this subject.