Imaging-based clusters in former smokers of the COPD cohort associate with clinical characteristics: the SubPopulations and intermediate outcome measures in COPD study (SPIROMICS)

Background Quantitative computed tomographic (QCT) imaging-based metrics enable to quantify smoking induced disease alterations and to identify imaging-based clusters for current smokers. We aimed to derive clinically meaningful sub-groups of former smokers using dimensional reduction and clustering methods to develop a new way of COPD phenotyping. Methods An imaging-based cluster analysis was performed for 406 former smokers with a comprehensive set of imaging metrics including 75 imaging-based metrics. They consisted of structural and functional variables at 10 segmental and 5 lobar locations. The structural variables included lung shape, branching angle, airway-circularity, airway-wall-thickness, airway diameter; the functional variables included regional ventilation, emphysema percentage, functional small airway disease percentage, Jacobian (volume change), anisotropic deformation index (directional preference in volume change), and tissue fractions at inspiration and expiration. Results We derived four distinct imaging-based clusters as possible phenotypes with the sizes of 100, 80, 141, and 85, respectively. Cluster 1 subjects were asymptomatic and showed relatively normal airway structure and lung function except airway wall thickening and moderate emphysema. Cluster 2 subjects populated with obese females showed an increase of tissue fraction at inspiration, minimal emphysema, and the lowest progression rate of emphysema. Cluster 3 subjects populated with older males showed small airway narrowing and a decreased tissue fraction at expiration, both indicating air-trapping. Cluster 4 subjects populated with lean males were likely to be severe COPD subjects showing the highest progression rate of emphysema. Conclusions QCT imaging-based metrics for former smokers allow for the derivation of statistically stable clusters associated with unique clinical characteristics. This approach helps better categorization of COPD sub-populations; suggesting possible quantitative structural and functional phenotypes. Electronic supplementary material The online version of this article (10.1186/s12931-019-1121-z) contains supplementary material, which is available to authorized users.


Background
Chronic obstructive pulmonary disease (COPD) is the third leading cause of death in the United States [1] and is identified by airflow limitation and/or obstruction. The severity of COPD is assessed by forced expiratory volume in 1 s (FEV 1 %) predicted values at post bronchodilator [2]. The pulmonary function test (PFT)-based FEV 1 and forced vital capacity (FVC) values are highly recommended to assess the global alteration of lung, but they do not correlate well with symptoms [3]. In addition, PFTs do not reveal local structural and functional alterations, which are essential in examining the heterogeneity of COPD phenotypes. Thus, the ability to quantify these alterations at multiple scales during COPD progression is necessary to characterize COPD phenotypes.
A multicenter study of COPD, i.e., Subpopulations and Intermediate Outcomes in COPD Study (SPIRO-MICS) [2] acquired QCT scans at total lung capacity (TLC) and residual volume (RV) [4]. This is an integral part of the multicenter study to find structural and functional phenotypes. A recent advance of quantitative medical imaging and data analysis techniques allows for derivation of QCT imaging-based metrics, leading to identification of statistically stable clusters/ phenotypes. For instance, using only QCT imagingbased variables, Choi et al. [5] derived clinically meaningful asthmatic sub-groups, being potentially useful in developing clusters-specific treatments. Furthermore, Haghighi et al. [6] expanded the QCT imagingbased clustering approach to identify homogenous clusters within current smokers from SPIROMICS. In this study, we hypothesize that QCT-based imaging metrics could be used to identify distinct COPD former smoker sub-groups with clinically meaningful characteristics, subsequently adding insights to the previous study of current smokers [6]. Shaker et al. [7] and Zach et al. [8] reported that former smokers had significantly higher % low-attenuation areas (%LAAs) on inspiration and expiration CT scans (for emphysema and air trapping measures) than current smokers. This is possibly due to parenchymal inflammation in current smokers serving to mask CT-based indices relative to former smokers [6,7]. Therefore, we divided the subjects into former and current smokers to independently assess phenotypes between these two groups and report on the former smokers in this work.
With the aid of machine learning techniques, QCT imaging-based metrics have been used to find homogeneous sub-groups of COPD subjects. As an example, Bodduluri et al. [9] have employed image registrationbased metrics to discriminate COPD subjects from non-COPD subjects. The study demonstrated the potential of registration-based variables to characterize COPD phenotypes, but this study was limited in supervised learning. In regards to unsupervised learning methods, there have been several efforts to identify COPD subgroups, but they employed either clinical data-only or a mix of clinical and CT data together [10][11][12] as we focus on imaging-only parameters to identify clusters. Although it would be possible to add clinical/physiological/biological measures into our cluster analysis, we used only imaging-based features to focus features of airway structure and lung function. Once established, our clusters were evaluated for their clinical, physiological, or biological measures. The associations between imaging-only clusters and the non-imaging phenotypes provide a validation of the ability of imaging metrics to characterize clinically meaningful phenotypes. Choi et al. [5] pioneered the use of unsupervised cluster analysis using CT image data acquired by the Severe Asthma Research Program (SARP) to identify four asthmatic clusters. Their approach accounted for inter-site and intersubject variations, enabling an analysis of large data sets acquired by multiple centers. Furthermore, Choi et al. [13] successfully identified imaging-based structural and functional features that differentiate asthmatics and COPD patients with chronic functional alteration.
In this study, we adopted the approach by Choi et al. [5]. In addition to the existing imaging-based metrics developed for asthma, we introduced several new metrics to account for tissue alterations and emphysematous lung [5,6]. A comprehensive set of imaging-based metrics were transformed to the principal component domain, and a cluster analysis was performed to explore possible COPD phenotypes of former smokers. The former smokers-clusters were then evaluated in association with severity, GOLD stages [14], sex, BMI and biomarkers, such as neutrophil counts, leukocyte (WBC) count and matrix metalloproteinase (MMP-3). We then compared the cluster membership of former smokers in this study with that of current smokers presented in our previous study [6].

Human subject data and QCT imaging
We analyzed a total of 758 SPIROMICS subjects containing an extensive set of biomarkers. In our analysis, we hypothesized that smoking status may have effects on CT measures of former and current smokers [7,8]. The hypothesis was further consolidated by performing a combined analysis and finding that a mix of both groups cannot provide adequate cluster stability. Hence, we excluded current smokers, so that a total of 406 formers smokers remained. The healthy never smokers without COPD were considered as healthy controls and were not included in the clustering analysis. PFTs were performed for all subjects pre-and post-bronchodilator, and CT was performed post-bronchodilator. Table 1 shows the demographic and PFT measures based on each stratum. Former smokers with post-bronchodilator FEV 1 /FVC > 0.7 were grouped in stratum 2, and former smokers in strata 3 and 4 had post-bronchodilator FEV 1 /FVC < 0.7, with FEV 1 > 50% in stratum 3 and FEV 1 < 50% in stratum 4, respectively [2].
Two QCT scans at TLC and RV were acquired by multiple imaging centers in the NIH-funded SPIRO-MICS multicenter research study [4]. The CT imaging protocols were approved by each center's institutional review boards (IRB). All QCT scans were obtained with post-bronchodilator. They were segmented with an automated commercial airway/lung segmentation software (Apollo 2.0, VIDA Diagnostics), and registered with a non-rigid mass-preserving imaging registration technique [15,16].

Derivation of QCT imaging-based metrics
A total of 75 multiscale imaging-based variables were extracted to derive principal components (Fig. 1). The segmental variables included bifurcation angle (θ), airway circularity (Cr), wall thickness (WT) and hydraulic diameter (D h ), where each variable indicated alteration of skeletal structure, alteration of airway shape, wall thickening and luminal narrowing, respectively. The sizes of WT and D h were normalized by tracheal WT and average diameter (D ave ) predicted from healthy subjects [5], being denoted by WT* and D h *, to eliminate inter-subject variability due to age, sex, and height. The four segmental variables were extracted from ten local regions to reflect characteristics of regional alterations. A detailed derivation of the above structural variables can be found in reference [17]. We further derived both strain-based and density-based functional metrics with the aid of image registration that matched two QCT images at TLC and RV. The strainbased variables included fractional air volume change (ΔV air F ), the determinant of Jacobian (Jacobian), and anisotropic deformation index (ADI). These are estimates of regional ventilation, local volume change, and preferential local lung deformation respectively [18,19]. Next, the density-based functional metrics included functional small airway disease percentage (fSAD%) and emphysema percentage (Emph%) to characterize the portions of small airway narrowing/closure and emphysematous lung, respectively. This approach was devised to dissociate emphysematous region from air-trapping region, previously proposed by Galban et al. [20]. In order to eliminate intersite variation, we employed a fraction-based fSAD% and Emph% using 90 and 98.5% air-fraction as the threshold, instead of using the density threshold of − 856 and − 950, respectively [21]. We further added two more imagingbased metrics that measure tissue fraction [13,22] at TLC and RV (β tissue TLC and β tissue RV ). The tissue fractions measure the portion of tissue volume in each voxel. These are supplementary metrics for Emph% and fSAD%, because β tissue TLC decreases if tissue destruction is captured and β tissue RV decreases if air fraction increases due to air-trapping.
In addition, we included global imaging-based metrics such as the ratio of apical-basal distance over ventral-dorsal distance at TLC (lung shape), the ratio of air-volume changes in upper lobes to those in middle and lower lobes between TLC and RV (U/(M + L)|v), fSAD%, Emph%, β tissue TLC and β tissue RV , Jacobian and ADI in the whole lung. Overall, there were 32 local/segmental structural variables, 35 lobar structural variables and 8 global structural variables.

Cluster and statistical analysis
Raw imaging data were scaled with standard scaler, and a principal component analysis was performed to derive linearly uncorrelated variables, so-called principal components (PCs). To obtain an optimal number of PCs, a parallel analysis [23] with random uncorrelated data was adopted. The analysis led to the number of 7 as an optimal choice of PCs (Additional file 1: Figure S1). Using the 7 derived PCs, to find the optimal clustering method and number, we then assessed internal properties including connectivity, average Silhouette width and Dunn indices [24] for three different clustering methods, i.e., hierarchical, K-means, and Gaussian finite mixture modelbased methods. Connectivity, average Silhouette width and Dunn indices measure the inverse of i th nearest neighbors which is not assigned to the same cluster, how tightly grouped all the points in the cluster are, and the ratio between the minimal inter-cluster distance to maximal intra-cluster distance, respectively. Thus smaller connectivity and larger Silhouette width and Dunn index indicate better clustering properties. First, K-means method was found to be a good clustering method for current data based on connectivity, and average Silhouette width (Additional file 2: Figure S2a). Dunn criteria then suggested that the number of 4 is an optimal choice in using K-means. To further test stability of the clustering membership, a nonparametric bootstrap analysis was performed with 200 bootstrapped data sets. The mean of Jaccard similarity coefficients, defined by the size of intersection divided by the size of the union between clusters [25], was computed to find the optimal cluster number and clustering approach (Additional file 2: Figure S2b).
Kruskal-Wallis and chi-square tests were performed to compare differences of continuous and categorical variables, respectively. The reported P values were significant, if any one group is statistically different from one group or more. We then performed association tests of imagingbased clusters with demographic and clinical variables to investigate the clinical relevance of current clusters.

Structural and functional features of imaging-based clusters
Cluster analysis identified four stable [6] imaging-based clusters with the sizes of 100, 80, 141 and 85, respectively ( Table 2, and Fig. 2). Five major variables with higher Wilk's λ values which best describe the four clusters were selected with a stepwise forward variable selection technique using Wilk's λ criterion [26]. Note that the clusters were differentiated predominantly with whole lung (total) parenchymal metrics including β tissue at RV and TLC, Jacobian, Emph% and fSAD%. Overall whole lung Emph% and fSAD% increased with increasing cluster number. It was noted that Emph% and fSAD% in Cluster 2 fell within the similar range with healthy subjects (Fig. 3).
Structural alterations in segmental airways were also captured between clusters (Table 3). Tracheal bifurcation angle (θ) and circularity (Cr) measured in the sLUL were significantly reduced in Cluster 4. Cluster 1 was characterized by airway wall thickening (WT*↑), whereas Clusters 3 and 4 were demonstrated by airway wall thinning (WT*↓) and airway narrowing (D h *↓). As summarized in Fig. 2, clusters were characterized by airway wall thickeningdominance (Cluster 1), increased tissue fraction at TLC with marginally increased emphysema (Cluster 2), proximal and peripheral airway narrowing (Cluster 3), and severe alterations of tracheal bifurcation angle (θ) and airway shape (Cr) on proximal airways as well as peripheral alterations (Cluster 4).

Associations of imaging-based clusters with clinical features
Clusters 1 and 2 were mostly populated in GOLD 0, 1 and 2 along with a lower BODE index, while Cluster 4 was mostly populated with GOLD 3 and 4 (stratum 4) with the highest BODE index (  The smoking pack-years were significantly greater in Clusters 3 and 4 than those of Clusters 1 and 2 (Table 5). Cluster 4 showed higher associations with pulmonary/ vascular condition, and chronic bronchitis, emphysema, and COPD diagnosed at baseline across all clusters. Shortness of breath during sleep was increased in Clusters 2 and 4. Fathers and mothers of subjects in Cluster 4 were likely to have COPD. The WBC counts were increased in Clusters 2-4, with increased neutrophils (Table 6). Lymphocytes were reduced in Cluster 4. The proteolytic enzymes of matrix-metalloproteinases (MMPs) were reduced especially in Cluster 2. Based on the lowest CAT score and exacerbation, Cluster 1 subjects were likely asymptomatic (CAT< 10) former smokers with the lowest exacerbation across all clusters. In contrast, Cluster 4 showed the highest CAT score  Fig. 3 a Percentage of emphysema (Emph%) for four clusters and the healthy control group (green). † P > 0.05 between clusters 1, 2, 3 and the healthy group. P < 0.05 between Cluster 4 and other groups for all pairwise comparisons b Percentage of small airway disease (fSAD%) for four clusters and the healthy control group (green). ‡ P < 0.05 for comparisons between four clusters 2, 3, 4 (red) and the healthy group for all pairwise comparison. P > 0.05 for between Cluster 1 and the healthy group with the lowest 6-min walk distance along with severe oxygen desaturation.
We further associated the clusters with visual diagnostic assessments including COPD subtypes (CLE: Centrilobular; PSE: Paraseptal; PLE: Panlobular emphysema) as well as interstitial lung disease (ILD) by an experienced thoracic radiologist at the University of Iowa (Table 7) because these subtypes might be associated with airway abnormalities [27]. Cluster 4 was less likely related to ILD and had a significant increase of PLE. Subjects with PLE were not observed in Clusters 1 and 2. We analyzed longitudinal data of 169 available subjects among the current cohort of former smokers to quantify change of Emph%, i.e., emphysema progression index (ΔEmph%) between baseline and oneyear follow-up. ΔEmph% is computed as the percentage of  Post-bronchodilator values after six to eight puffs of albuterol. Full names of each variable were described in Abbreviations used. BODE indexes for 24 subjects were not available voxels within the lung less than − 950 HU and assesses the extend of emphysema (ΔEmph% ≥ 1% and ΔEmph% ± 0.5% are considered as rapid-progressors and nonprogressors, respectively) [28]. ΔEmph was marginal in Cluster 2 (Table 7), whereas it was significantly higher in Cluster 4. Furthermore, we compared two different clustersgrouping derived from current smokers [6] and former smokers, respectively (Table 8). Overall CAT score and exacerbation histories of current smokers were greater than those of former smokers. WBC counts were not differentiable in current smokers-derived clusters because all clusters showed large numbers of WBC count. On the other hand, WBC count of former smokersderived Cluster 1 was the smallest and it was increased as increasing the cluster membership of former smokers. On the contrary to the finding of WBC counts, former smokers demonstrated greater Emph and fSAD% than current smokers, based on kernel density estimation (KDE) plots (Fig. 4). The dispersed density distribution of current smokers may indicate the masking effect of CT-based measures of emphysema and small airway disease, compared to former smokers [7]. The Emph% and fSAD% of former smokers (Table 2) were especially increased in Clusters 3 and 4, as compared with counterparts of current smokers [6].

Decision tree analysis
We performed a decision tree analysis to construct a simple predictive model (Additional file 3: Figure S3) to classify former smokers. The data set was shuffled randomly into training (n = 324) and test sets (n = 82) and the accuracy was assessed on the test set. The model comprising 5 discriminant variables resulted in accuracy of 81%. These variables were β tissue RV (Total), Jacobian (Total), β tissue TLC (Total), D h * (RMB) and ADI (Total). We further evaluated an association between current and former smoker clusters by assessing the membership of former smokers in the decision tree of current smokers [6] and vice versa. The classification accuracy for both cases was about 0.62 based on the confusion matrices (Additional file 4: Table S1). It can give an assessment for possible overlap between clusters of these two cohorts.

Discussion
In this study, we applied an unsupervised clustering method with an expanded set of imaging-based variables to former COPD smokers collected in the multicenter study of SPIROMICS. Four homogeneous clusters were derived within a former-smoker population, exhibiting distinct phenotypic characteristics and strong associations with clinically relevant COPD biomarkers. The imaging-based clusters can provide more information than the conventional PFT-based classification of COPD, such as stratum and GOLD criteria, because they explain structural and functional alterations at lobar and segmental levels. We also included parenchymal metrics including Emph%, fSAD%, tissue fractions at TLC and RV as well as segmental-level structural metrics including wall thickness and diameter of airway branches. The imaging and clinical phenotypes based on the clusters could be explained as follows.

Features of respective clusters
The cluster memberships can suggest possible phenotypes with distinct characteristic correlated with relevant clinical/biomarker measures for former COPD smoker.

Cluster 1: asymptomatic resistant smokers with preserved pulmonary function
Cluster 1 showed preserved pulmonary function (FEV1/ FVC = 0.72) at post bronchodilator and was mostly populated in GOLD stages 0 and 1. This cluster had a relatively low Emph% and fSAD% with structural and functional characteristics close to those of healthy controls. BODE index, exacerbation histories and WBC count of this cluster were relatively lower compared with other clusters. These characteristics along with CAT< 10 and the lowest exacerbation among all clusters suggests that Cluster 1 belongs to asymptomatic resistant smokers. Cluster 1 imaging metrics were very close to those of healthy subjects. Airway wall thickening was the only abnormality in this  [29], reported that long-term smoking may contribute to airway wall thickening prior to the development of more severe imaging features of COPD.

Cluster 2: obese female individuals with preserved lung function and marginal emphysema
Cluster 2 with the highest BMI and over-representation of women indicated clinical and epidemiological importance as reported by Castaldi et al. [10] and Martinez et al. [30]. Castaldi et al. [10] derived four clusters with 10,192 subjects from COPDGene using several imaging-based metrics, e.g., Emph%, upper/lower ratio of Emph%, gas trapping, and PFT results acquired by a feature selection method. Note that our Cluster 2 is aligned with Cluster 2 of Castaldi et al. [10] in high BMI, African-American and women-dominance. Cluster 2 showed the preserved pulmonary function (FEV/FVC = 0.71) close to Cluster 1, but the CAT score and exacerbation of this cluster was greater than that of Cluster 1. This group showed a noticeable increase of tissue fraction at TLC, and a decrease of emphysema index among clusters. This cluster included more CLE-only type while showing the lowest ΔEmph% among clusters. This finding is of interest because most studies showed that development of CLE is associated with severe abnormalities of the small airways, e.g. wall thickening. Thus, CLE may be more related to air-borne risk factors that cause airway inflammatory processes [27]. Cluster 2 also showed the lowest value of MMPs among clusters. Ostridge et al. [31] investigated the association between specific pulmonary MMPs and emphysema as these enzymes degrade the extracellular matrix and have been identified as potentially important in the development of emphysema [31].

Cluster 3: older male individuals with increasing fSAD and emphysema
Unlike Clusters 1 and 2, Cluster 3 demonstrated a significant decrease of FEV1/FVC and FEV1% predicted values, but their FVC % predicted value remained in the normal range. This cluster was mostly populated in GOLD stages 2 and 3 with a significant increase in BODE index. From this cluster, Emph% and fSAD% in parenchymal regions were significantly increased, being similar with Cluster 4. Thus, this cluster showing airway  This cluster showed the highest Emph%, fSAD%, BODE index, WBC count and CAT score along with the lowest FEV1/FVC among all clusters. These characteristics along with structural and functional variables indicated that Cluster 4 belongs to severe symptomatic COPD subjects. The pattern of decreasing D h * with increasing fSAD% (non-emphysematous air trapping) indicates severely narrowed status of both proximal and distal airways. In addition to airway narrowing, this group actually contains most of the significant structural and functional alterations. It is especially noted that prominent airway wall thinning and alteration of airway geometry change were only observed in this cluster. Assuming that this cluster is the most severe COPD group, alterations of airway features including airway wall thinning (WT*), elliptic airway shape (Cr), and change of airway geometry (θ) may occur at the end stage of COPD. Dominance of PLE with diffuse destruction in Cluster 4 along with its highest progression index among all clusters might be related to blood-borne mechanism rather than the possible air-borne mechanism in Cluster 2. These finding shows the possibility of two different pathogenetic mechanisms among subjects. In addition, Koo et al. [32] studied WBC count as a biomarker and their associations with the severity of the disease. WBC count in former smokers has an increasing pattern from Cluster 1 to Cluster 4 (Table 6) along with increasing CAT score and decreasing FEV1/FVC.
With previously analyzed current smokers [6], the comparison for important clinical and biomarker measures between former and current smokers are shown in Table 8. Overall, exacerbation has increasing pattern between clusters of former smokers with Cluster 1 and Cluster 4 with the lowest and highest, respectively. Cluster 2 for both current and former smokers has increased exacerbation compared to clusters 1 and 3 and might be related to the highest tissue fraction and possible inflammation in Cluster 2.
WBC count was lower in former smokers possibly due to the effect of smoking on the WBC [6], which was also significantly elevated as increasing cluster membership. This result indicates that WBC count can serve as an important risk factor such as inflammation especially in former smokers. Furthermore, the CAT score and exacerbation histories were significantly higher in current smokers than in former smokers. An increase in inflammatory markers in current smokers relative to former smokers was contradictory to imaging-based features such as Emph% and fSAD% (Fig. 4). The smoking status could affect parenchymal inflammation, leading to an increase of CT density [6,7]. Thus Emph% and fSAD% could be underestimated, if patients are on smoking. This confounding effect prevents from applying a clustering algorithm for former and current smokers due to the low Jaccard index (< 0.7).
To assess a possible overlap between current and former smokers, we used the trained decision tree on current smokers to classify former smokers and vice versa; the classification accuracy for both cases was about 0.62 (the confusion matrices are reported in Additional file 4: Table S1). This result indicates that two clustering analyses between former and current smokers can be further used to investigate the difference in phenotypic characteristics of these cohorts. The impact of smoking status on cluster membership requires further investigation with larger cohorts as well as with longitudinal data to inspect disease progression and membership transition over time.

Conclusions
We performed a cross-sectional study to derive four unique imaging-based clusters in former smokers with COPD. The current cluster analysis can be used in conjunction with our previously reported cluster analyses in current smokers with COPD to assess the differences in smoking status (former vs current) in the COPD population and explore possible different phenotypes between these two groups.

Additional files
Additional file 1: Figure S1. A scree plot: eigenvalues (magnitude of variances) according to the number of principal components for determining the optimal number of components. (DOCX 65 kb) Additional file 2: Figure S2. (a) Internal properties in different clustering methods to find the best clustering approaches as well as the optimal number of clusters; (b) Bootstrapping stability analysis between K-means and hierarchical clustering with 4 or 5 numbers of clusters. (DOCX 58 kb) Additional file 3: Figure S3.