Quantitative vessel tortuosity: A potential CT imaging biomarker for distinguishing lung granulomas from adenocarcinomas

Adenocarcinomas and active granulomas can both have a spiculated appearance on computed tomography (CT) and both are often fluorodeoxyglucose (FDG) avid on positron emission tomography (PET) scan, making them difficult to distinguish. Consequently, patients with benign granulomas are often subjected to invasive surgical biopsies or resections. In this study, quantitative vessel tortuosity (QVT), a novel CT imaging biomarker to distinguish between benign granulomas and adenocarcinomas on routine non-contrast lung CT scans is introduced. Our study comprised of CT scans of 290 patients from two different institutions, one cohort for training (N = 145) and the other (N = 145) for independent validation. In conjunction with a machine learning classifier, the top informative and stable QVT features yielded an area under receiver operating characteristic curve (ROC AUC) of 0.85 in the independent validation set. On the same cohort, the corresponding AUCs for two human experts including a radiologist and a pulmonologist were found to be 0.61 and 0.60, respectively. QVT features also outperformed well known shape and textural radiomic features which had a maximum AUC of 0.73 (p-value = 0.002), as well as features learned using a convolutional neural network AUC = 0.76 (p-value = 0.028). Our results suggest that QVT features could potentially serve as a non-invasive imaging biomarker to distinguish granulomas from adenocarcinomas on non-contrast CT scans.


A. Extraction of QVT features:
We extracted 35 QVT features from the nodule vasculature in each CT scan. The extracted QVT features are based on three important properties of the vessels: (I) tortuosity, (II) curvature and (III) branching statistics and the volume measurements of the vasculature. The nodule vasculature is assumed to be comprised of several vessel branches where each branch in turn is comprised of several points in a 3D space. The QVT features were derived from measurements relating to points, vessel branches and the entire vasculature.
The curvature of a point which is located on the center line of a vessel in 3D space is defined as the inverse of the radius of an osculating circle fitted to that point with respect to its immediate neighbors. The statistical moments of the curvature measurements of the points associated with the branches are computed to describe the curvature of a branch and similarly the curvature of the entire vasculature.
The torsion of a branch is defined as the ratio of the Euclidean distance between the starting and end points of a branch to the length of the branch. The torsion of a whole vasculature, is then calculated by computing the first order statistics of the torsion measurements associated with the branches of a vasculature. In addition to torsion and curvature, the branching statistics of the vasculature was computed as well. These measurements refer to the number of small vessel branches in the vasculature.
Each nodule vasculature is comprised of several vessel branches as illustrated in Fig.S1(a) and Fig.S1(b), where each vessel branch in turn is comprised of several points in 3D space ( Fig.S2(a)). Notations:

Notation Definition V
Nodule vasculature which is comprised of several branches. = ⋃ =1 A curve corresponding to the center line of a vessel branch in 3D. = ⋃ =1 A 3D point on a with coordinates of ( , , ) | −1 − | Euclidean distance between −1 and Length of a Euclidean distance between the first ( 1 ) and the last ( ) points of a Torsion of Curvature of a point Curvature of a vessel branch Formulations: Assume that the vasculature (V) for each nodule is comprised of branches which were converted to its corresponding skeleton in the form of a line ( ) in 3D space. We also assume that a 3D line ( ) is in turn comprised of several points ( ). We denote the Euclidean distance between two immediate close points ( −1 , ) on by | −1 − | and hence the length of a vessel branch ( ) is then approximated as follows, where is the number of points on a vessel branch and | −1 − | is computed as follows: Where( −1 , −1 , −1 ) and ( , , ) correspond to the coordinates of −1 and respectively.
We also compute the direct distance between 1 and which corresponds to the length of in Fig.S2(a) as follows: = | 1 − |.
Hence, the torsion of a branch is defined as: The torsion of nodule vasculature (V), is then calculated by computing the first order statistics of the torsion measurements ( ) associated with the branches of a vasculature.
To obtain the curvature of a vessel branch, first, we calculate the curvature of its constituent points (see Fig2.b). Let ⃗( ) = ( ( ), ( ), ( )) be a vector-valued function of parameter that models a smooth vessel branch curve and let be the point on at . The function ⃗ produces a point P on by taking an input value . The curvature at is proportional to the inverse of the radius of an osculating circle at that point. The osculating circle at is the circle that best approximates at . The radius ( ) of the osculating circle at is defined as: where is time parameter and ( ) is the curvature at point and is obtained from the following formula: where and 2 2 denote the first and second order derivative with respect to the parameter t and ( ), ( ) and ( ) are coordinate functions of the parameter t.

For
at each point ∈ , we calculate Min, Max, Mean and Standard deviation (Std) of values to describe the curvature of a vessel branch ( ). Then, the curvature of the whole vasculature (V) is acquired by computing the first order statistics of the aforementioned measurements associated with the branches of a vasculature. Table.S1 illustrates the description of QVT features.     (n=145). AUC values were generated for each of the human readers as well as the machine classifier (SVM). A hard decision was obtained for each case to be a granuloma by using a threshold of <3 and >50% respectively for the human readers and the machine classifier. Each green cell represents a decision either by the machine classifier or the human reader when the hard decision matched the true pathologic diagnosis for that nodule and pink when it did not match.  D. Inclusion and exclusion criteria of the data:

E. Comparison of feature selection strategies:
The AUC performance of the QVT features selected by mRMR, LASSO and PC as well as joint features selected by LASSO and mRMR are presented in Table.S5. Among the 12 features identified by mRMR, 8 features were also selected by LASSO. We believe that mRMR and LASSO are potentially more appropriate for the features considered in this study compared to PCA. This is potentially because PCA transforms the original features into a new feature space and uses a linear embedding of the high dimensional features which may actually lie on a non-linear manifold. Hence the linearity assumptions involved in the use of PCA may result in low dimensional embedding which may not be an accurate reflection of the structure of the original high dimensional feature space.
F. The impact of nodules position to the output of the method:

F.1. Association of the nodule location and diagnostic class:
To determine the association of a nodule's location with its corresponding diagnostic class, we captured the location information of each of the nodules. This included whether a nodule was located in the apical, mid or basal lung regions in 3D (see Fig.S4(b)). In the 2D transverse plane, we also identified whether a nodule was located within the central, upper or lower lung region (see Fig.S4(a)). Figure S4. Upper, central and lower control regions on transverse (a) and sagittal (b) planes. The transverse plane was used to determine whether a nodule was close to the upper pleura, sub pleura or whether it was located in the central lung region. The sagittal plane was used to determine whether the nodule was located in the apical, mid or basal regions along z-axis. Fig.S5 shows a frequency plot for the nodule position in both 3D and 2D for both adenocarcinomas and granulomas. Also we performed 2 test between nodule position and its class to discover possible dependency between these parameters. We have defined control positions and assigned 1, 2, and 3 respectively for nodules located in the upper, central and lower regions in 2D and 3D. The 2 test revealed that variables corresponding to both 2D position of the nodule in lung (Fig.S4-a) and 3D position ( Fig.S4-b) were independent from nodules' diagnostic class. In other words, no significant association was identified between the position of the nodules and its corresponding diagnostic class (p>0.01). These results can be seen in Table  S6.  F.2. The impact of spatial location of a nodule on its corresponding QVT feature measurements: In this section we evaluated to what extent the spatial location of a nodule impacts the corresponding QVT measurements of the nodules from the different diagnostic classes. First, we defined the following control regions to determine whether a nodule was located in the a) apical or lower lung portions, b) Gravity dependent (the lowest part of the lung in relation to gravity) or non-dependent portions of the lung and c) peripheral or central regions. Fig.S6 shows the control regions. Next, the QVT features for the 145 cases of the training set where extracted. Then for both adenocarcinomas and granulomas we determined whether the lung region impacted the QVT measurements. In this regard, an unpaired t-test was applied to the 10 top ranked QVT features of the nodules corresponding to the apical vs lower, dependent vs non-dependent and central vs peripheral lung regions. Fig.S7 shows the QVT features of the apical and lower lung adenocarcinomas. Similarly, Fig.S8 shows the QVT features corresponding to the apical and lower lung granulomas. As may be appreciated from figures S7-S8, QVT features were not found to be significantly different between apical and lower lung regions. As presented in figures S9 and S10, 1 QVT feature was found to be significantly different for the nodules located in the dependent and non-dependent lung regions. The statistically significantly different features were 1: Mean torsion of the vessel branches for adenocarcinomas, and 10: the ratio of vasculature volume to its bounding box for granulomas.
Finally, 5 QVT features (5: Branching count of the vasculature, 6: histogram of curvature, 7-10: torsion histogram) were significantly different between adenocarcinomas of central and peripheral regions. For the granulomas within the same regions, 3 QVT features (1,2: Mean, Std values of the torsion of the vessel branches, 10: histogram of torsion) were found to be significantly different (see Figs.S11-12). a b c Figure S6. Control regions to determine whether a nodule belong to a) apical or lower lung portions, b) dependent or non-dependent lung and c) central or peripheral regions. Figure S7. No significant differences were found between QVT features of the adenocarcinomas located in the apical and lower lung regions. Figure S8. No significant differences were found between QVT features of the granulomas located in the apical and lower lung regions. Figure S9. 1 QVT feature (1: Mean of the vessel branches' torsion) was significantly different between adenocarcinomas corresponding to regions from within the gravity dependent and non-dependent lung regions. The lowest part of the lung in relation to gravity was considered as the dependent region and the upper part was considered as non-dependent region.
Figure S10. 1 QVT feature (10: the ratio of vasculature volume to its bounding box) was significantly different between granulomas corresponding to the regions from with gravity dependent and non-dependent lung. The lowest part of the lung in relation to gravity was considered as the dependent region and the upper part was considered as non-dependent region.

G. Sensitivity of QVT features to errors in vasculature segmentation:
This section presents the sensitivity of QVT features to minor changes/ errors of the vasculature segmentation. We identified the parameter S which controls the growing of the vasculature volume during the segmentation process. This parameter corresponds to the mean difference value in Hounsfield units between a growing volume and the candidate voxels in the neighborhood of a vasculature. The values of S were identified empirically and belong to the following set: S ∈ {−150, −100, −50,0, 50} . For each of the values in S, we generated multiple segmentations for each nodule and its surrounding vasculature on the CT scan. Fig.S13 (b,c,d,e,f) illustrates 5 segmentations of a nodule and its corresponding vasculature shown in Fig.S13 (a). Segmentations were generated based on the values of S. ( 1 , 2 ) . Fig.S14 shows the 4 bar plots in which Fig.S13 (a,b, . As it may be appreciated from the bar plots, some features (13,14,15,16) were highly stable and correlated across changes in vessel segmentation. However, some features (17,18,19,20) were found to show more abrupt changes on account of minor perturbations in vessel segmentation. Consequently, as shown in the panels, these features were not highly correlated between different segmentations. . As it may be appreciated from the bar plots, some features (13,14,15,16) were highly stable and correlated across changes in vessel segmentation. Specially, in panel (c) which corresponds to ( − , − ) most of the QVT features have shown high or moderate correlation. However, some features (17,18,19,20) were found to show more abrupt changes on account of minor perturbations in vessel segmentation. Consequently, as shown in the panels, these features were not highly correlated between different segmentations.
Additionally, we also computed classification performance of the classifier as a function of . Fig.S15 and Table.S7 show the AUC values of the classifier for each , S ∈ {− , − , − , , } . H. The effect of scanner variability and voxel size and smoking history on the QVT features: The effect of scanner variability and voxel size on the QVT features. Fig.S16 presents the distribution of voxel size in both training and validation sets. The statistical significance test result of voxel size and disease outcome is presented in Table S8. Additionally, we divided validation set into two sub-groups based on the voxel size (VS). There were 73 cases with VS>0.7 mm and 0.72 cases with VS<=0.7 mm. The AUC performance of the tortuosity based classifier in these was found to be 0.7 for VS<=0.7 and 0.81 for VS>0.7. Table S9 below provides details regarding the vendor and voxel size effect on the QVT classifier. To show the impact of smoking history on the vessel architecture, we used a stratified sub set including 21 never smoker patients and 22 patients who had had at least 20 pack years. This stratified set include 8 never smokers and 9 smokers in Adenocarcinoma class. The Granuloma class included 13 never smokers and 13 smokers. QVT features were extracted for this subset and the QVT classifier validated on it. Table S10 shows details of the stratified sub set and AUC performance of the QVT based classifier on this dataset. Additionally, the p-values of the null hypothesis that the difference of 5 selected QVT features between "Never smoked" and "Smoker" groups were found to be 0.61, 0.17, 0.03, 0.41, 0.65. In other words, no statistically significant difference were found for 4 out of 5 selected QVT features among "Never smoked" and "Smoker" groups .
I. The performance metric of the deep CNN across the training runs: Figure S17. Training and Validation metrics across 100 epochs for the CNN. This plot shows that the model loss decreases over the epochs, and approaches zero. Similarly, model accuracy approaches 1 as number of epochs increases, and plateaus over 80 epochs.
J. Multi-fold cross validation on the entire dataset: In this section, multi-fold cross validation was performed on the entire dataset, rather than just on the training set, to determine whether there was a difference in the selected features and the performance of the classifier. In this regard, we employed a 3-fold cross validation scheme (repeated 200 times) on the entire dataset which yielded an AUC=0.82±0.05. The corresponding ROC is shown in Fig.S18. Feature selection was performed within each iteration of the cross validation process and during each of the 600 iterations, we recorded the number of times that a feature was ranked among top three. For each QVT, Fig.S19 shows the number of times a feature was identified as being in the top 3 across multiple folds of cross-validation. As illustrated in Fig.S19, the features which were selected more than 100 times within the top ranked features, correspond to the top 12 features which were selected by our method and reported in the main manuscript. ( there were 7 joint features between the feature selection strategies which was performed on training set only or on the entire dataset). Figure S18. 3-fold cross validation across the entire dataset yielded an AUC=0.82±0.05 Figure S19.A frequency plot showing the number of times a QVT feature was identified as being within the top 3 ranked features across multiple runs of the cross-validation procedure.