2.1.

Study Design and Participants

All participants provided informed written consent to a study protocol approved by a local research ethics board and in compliance with the Health Canada approved and registered protocol¹⁹ (clinicaltrials.gov NCT02279329, Institutional Review Board IRB00000940). Inclusion criteria were a history of cigarette smoking $>10$ pack years, age between 50 and 85 years at baseline. Ex-smoker participants were included who had ceased smoking $\ge 1$ year prior to the study visit, with no cut-off in terms of smoking cessation. COPD subjects were classified according to the Global Initiative for Chronic Obstructive Lung Disease (GOLD) grades.²⁰ Participants also completed a longitudinal follow-up visit at $24\pm 6$ months after the baseline visit.¹⁹ Participants that reported a change in smoking status from ex-smokers to current smokers at follow-up visit were included in the analysis. The CONSORT diagram for the Thoracic Imaging Network of Canada (TINCan) study cohort participants is depicted in Fig. 1.

Fig. 1

CONSORT flow diagram. Of the 266 participants enrolled in the TINCan study, 33 were enrolled in a sub-study and 61 either cancelled or did not complete all required tests during visit 1. Of the 172 participants that completed visit 1, 79 participants did not complete a 3-year follow-up visit 2 and 5 had artifacts present in their images, which were excluded from further analysis. Eighty-eight ex-smoker participants with visit 1 and visit 2 data were analyzed in this study, of which 31 had $\mathrm{\Delta }{\mathrm{FEV}}_1\ge 60\mathrm{mL}/\text{year}$ and 57 had $\mathrm{\Delta }{\mathrm{FEV}}_1<60\mathrm{mL}/\text{year}$.

2.2.

Pulmonary Function Tests and Image Acquisition

Spirometry, plethysmography, and the diffusing capacity of the lungs for carbon monoxide (${\mathrm{DL}}_{\mathrm{CO}}$) were measured according to the American Thoracic Society/European Respiratory Society guidelines²¹ using a whole-body plethysmography system (MedGraphics Corporation, St Paul, Minnesota, United States) and attached gas analyzer.²² COPD was defined as post-bronchodilator spirometry according to the GOLD criteria.²⁰ The 6MWD²³ test and St. George’s Respiratory Questionnaire (SGRQ)²⁴ were administered under supervision.

Anatomic proton (${}^1\mathrm{H}$) and hyperpolarized ${}^3\mathrm{He}$ MR images were acquired using a whole-body 3.0 Tesla Discovery MR750 system (GE Healthcare, Milwaukee, Wisconsin, United States), a whole-body radiofrequency coil and a fast gradient recalled echo (FGRE) sequence with a partial echo implementation, with acquisition parameters as previously described.²⁵ Hyperpolarized ${}^3\mathrm{He}$ MRI was acquired using a linear bird-cage transmit/receive chest coil (RAPID Biomedical GmbH, Wuerzburg, Germany). A turn-key system (HeliSpin™, Polarean Inc, Durham, North Carolina, United States) was used to polarize ${}^3\mathrm{He}$ gas to 30% to 40% and doses ($5\mathrm{mL}/\mathrm{kg}$ body weight) diluted with ${\mathrm{N}}_2$ were administered in 1.0 L Tedlar^® bags. Hyperpolarized ${}^3\mathrm{He}$ MRI diffusion-weighted imaging was performed using a 2D multi-slice fast gradient-echo method, as previously described,²⁵ during breath-hold for acquisition of two interleaved images with and without additional diffusion sensitization with $b=1.6\sec /{\mathrm{cm}}^2$ (maximum gradient amplitude [G] = 1.94 G/cm, rise and fall-time = 0.5 ms, gradient duration = 0.46 ms, and diffusion time = 1.46 ms).¹⁹ Pulmonary function data and imaging were acquired during both baseline and follow-up visits.

2.3.

Image Analysis and Proposed Algorithm

Baseline and follow-up visit ${}^1\mathrm{H}$ and ${}^3\mathrm{He}$ MR images were processed by a single observer where the thoracic cavity was segmented from the ${}^1\mathrm{H}$ images using a seeded region-growing algorithm, and the ${}^3\mathrm{He}$ ventilation region was segmented using k-means clustering.¹ The generated maximum entropy mask was then applied to identify the ventilated region-of-interest (ROI) for feature extraction. Diffusion-weighted images were automatically processed to generate apparent diffusion coefficient (ADC) values and images, as previously described.²²

MRI VDP was generated using a semi-automated segmentation approach, as previously described.¹ Ventilation defect cluster percent (VDCP), which is the sum of ventilation-defect cluster volume normalized to the total lung volume, and defect cluster sizes were generated by an automated in-house custom-developed algorithm, as described in the prior proceedings paper.²⁶ Briefly, the proposed cluster algorithm iteratively segments unventilated MRI volumes until the maximum sphere fit (or multiple spheres of the same size) within the unventilated volume is identified. This sphere volume(s) is then removed from the non-ventilated region and the process is iteratively repeated until the non-ventilated region is filled by spheres. This is similar in approach to sphere packing previously investigated for radio-surgical treatment planning.²⁷

The approach was implemented using a naïve greedy algorithm where $S=[b_1,b_2,\dots ,b_n]$ is a set with $n$ elements, where each element $b_n=B_n(r,l)$ is an open sphere of radius $r$ at locations $\overrightarrow{l}$. Determining the required minimum number of spheres of unequal sizes resulted in the following minimization problem

A volume constraint was also imposed such that the total volume of spheres was equal to the volume of the specified unventilated region $R$

The specified regions were filled with spheres, where the minimum sphere diameter was equal to one voxel ($5\times 5\times 5{\mathrm{mm}}^3$) such that the total volume of spheres with diameter equal to one voxel was equal to be the total residual volume not clustered into large (>one voxel) sphere sizes. To further simplify the problem, there were no location constraints on the sphere spatial position.

We used MATLAB R2021a (MathWorks) to solve the minimization problem and generated VDCP to calculate cluster-defect-diameter voxel size one (CDD1) or the number of small ventilation defects, which is the cumulative number of defect clusters of one voxel, as shown in Fig. 2. Low ventilation cluster slopes were also calculated based on the log–log relationship between the cumulative number of spheres and cluster size.^28,29

Fig. 2

Ventilation defect cluster volume output from custom-developed algorithm. A representative 66-year old male ex-smoker participant with COPD: ${\mathrm{FEV}}_1=1.93\mathrm{L}$, VDP = 25%, $\mathrm{\Delta }{\mathrm{FEV}}_1=-0.22\mathrm{L}$, and $\mathrm{\Delta }\mathrm{VDP}=1\%$ between visits. Three-dimensional isotropic ventilation volume shown in cyan and ventilation defects represented by different sphere sizes ranging from small (yellow = sphere diameters of 3 to 5 voxels) to large (red = sphere diameters of 9 to 13 voxels).

Further texture feature extraction was conducted using an open-source PyRadiomics software, detailed in the next section. Unlike the proposed algorithm that analyzes the unventilated region, texture features were extracted from the inhaled hyperpolarized gas distribution within the thoracic cavity.

2.4.

Feature Extraction and Selection Pipeline

The complete pipeline for processing the baseline and follow-up visit images is summarized in Fig. 3. First, we generate a maximum entropy mask by segmenting the ${}^3\mathrm{He}$ ventilation image as previously described.¹ We then use a custom-developed algorithm to calculate the ventilation defect clusters, described above, and the PyRadiomics platform³⁰ for extracting texture features. Texture features were calculated from gray-level histograms and matrices generated from the ROI of the original image.

Fig. 3

Image processing and model generation pipeline. Hyperpolarized noble gas MR images were analyzed using custom-developed algorithms to calculate ventilation defect clusters and texture features. Raw MR images and corresponding binary lung masks were used to generate the regions-of-interest (ROI), which were analyzed via the PyRadiomics platform to calculate shape and first- and higher-order texture features. Feature selection was performed using the training set, combined with PCA and Boruta analysis. Classification learner application was used for model generation, which was trained using a fivefold cross-validation scheme, with 20% of data reserved for the testing set.

We generated first-order texture features from the gray-level histograms and also evaluated the texture features calculated from run-length, gap-length, co-occurrence, size-zone, dependence, and neighborhood gray tone matrices using the PyRadiomics open-sourced platform (version 2.2.0) in Python environment (version 3.7.5).³⁰ Image processing filters were also applied for the extraction of wavelet band-pass filtering texture features. This further quadrupled the number of extracted features due to permutations of high-pass and low-pass filters for wavelet decomposition, resulting in 376 additional MRI texture features. Low-pass filtering in both directions (LL) assesses the lowest frequencies, low-pass filtering followed by high-pass filtering (LH) assesses horizontal edges, high-pass filtering followed by low-pass filtering (HL) assesses vertical edges and high-pass filtering in both directions (HH) assesses diagonal details.³¹

The resulting data and 496 extracted features were randomly assigned into training and testing sets with 80%/20% data split. Participants were randomly divided into five groups and models were trained using the training set and a fivefold cross-validation scheme. Therefore, training was performed during five iterations, whereby for each iteration, the model was trained on four groups and validated on one group. Classification learner application (MATLAB R2021a) was used for generating all the models using the training dataset, with 20% of unseen data reserved as the final testing set to determine model accuracy. Feature selection was performed on the training set using principal component analysis (PCA) and Boruta analysis to rank and determine non-redundant texture features significantly contributing to the predictive power of the machine learning models.

PCA begins by standardizing the data and calculating the covariance matrix, which describes the relationships between features. This matrix is decomposed into eigenvectors and eigenvalues, representing the directions of maximum variance and their magnitudes, respectively. The top eigenvectors, known as principal components, are selected based on their corresponding eigenvalues to capture the most variance in the data. These principal components are used to transform the original data into a lower-dimensional space while preserving as much variance as possible. Boruta analysis is a feature selection technique designed to identify pertinent variables amidst high-dimensional datasets and used a two-step correction for multiple testing. Leveraging a random forest algorithm in conjunction with shadow feature creation, Boruta evaluates the importance of predictor variables iteratively, distinguishing between genuine predictors and random noise. This methodological approach enhances model interpretability and reduces the risk of overfitting. After creating shadow features using the original features, they are concatenated with the original dataset and used for training the random forest classifier. Feature importance is computed for the highest-rated shadow features and all original features that are more important than the corresponding shadow features are kept. This is repeated for multiple iterations, keeping track of the most important features using $z$-score statistics to determine the final subset of features [number of trees in the forest = 200, maximum iterations = 300, maximum tree depth = 10 (branches), percentage of shadow feature threshold = 95%, and alpha-level = 0.05]. We utilized all wavelet band-pass filtering features that were available in the PyRadiomics platform (version 2.2.0), with the detailed mathematical descriptions of all the selected and extracted texture features provided in Table 3 and Table S1 in the Supplementary Material.³⁰

2.5.

Machine learning and Statistical Analysis

Once all the features and parameters were selected, machine learning models were generated based on (1) demographic measurements alone, (2) spirometry measurements alone, (3) imaging and texture measurements alone, and (4) combination of all available measurements. Fivefold cross-validation training was performed using several machine learning algorithms, including single classifiers and ensemble classifiers, to determine the best model for identifying accelerated disease progression. The data were standardized and hyper-parameter optimization was performed through MATLAB R2021a (Classification Learner App) for each model individually. We compared the performance of multiple machine learning algorithms, including variations of: logistic regression, Naïve Bayes,³² support vector machines (SVMs),³³ decision trees,³⁴ K-nearest neighbors (KNNs),³⁵ and four ensemble-classifiers: bagged and boosted trees,³⁶ subspace discriminant,³⁷ subspace K-nearest-neighbors,³⁷ and random under-sampling boosted (RUSBoosted) trees. Model performance was evaluated using the mean cross-validation area under the receiver-operator curve (AUC) from training, as well as sensitivity, specificity, and F1-score using the test set and models’ confusion matrix. The DeLong’s test was used to compare the performance of all machine learning models.³⁸

Statistical analysis was performed using SPSS Statistics v28.0 (IBM Statistics, Armonk, New York, United States). Shapiro–Wilk tests were used to determine the normality of the data and non-parametric tests were performed for data that were not normal. Differences between subject groups were determined using analysis of variance with post-hoc analysis using the Benjamini–Hochberg correction. The relationship between measurements was determined using Pearson and Spearman coefficients for parametric and non-parametric data, respectively. Holm–Bonferroni correction was applied for multiple comparison tests for the selected texture features. Results were considered significant when the probability of two-tailed type I error was less than 5% ($p<0.05$).

3.1.

Participant Demographics

We evaluated 88 ex-smoker participants, of which 49 had spirometry evidence of COPD and 39 with no spirometry evidence of COPD. As shown in Table 1, 57 participants (22 females/35 males, $70\pm 9$ years) demonstrated a smaller ${\mathrm{FEV}}_1$ decline ($<60\mathrm{mL}/\text{year}$) and 31 participants (7 females/24 males, $68\pm 9$ years) demonstrated a rapid decline in ${\mathrm{FEV}}_1$ ($\ge 60\mathrm{mL}/\text{year}$), between baseline and follow-up visit $31\pm 7$ months later. At baseline, there were significant differences only in forced vital capacity (FVC) between subgroups [FVC (L) $p=0.003$ and FVC (${\%}_{\mathrm{pred}}$) $p=0.03$, respectively]. As summarized in Table 2, these subgroups did not have any statistically significant differences at follow-up.

Table 1

Baseline participant demographics and pulmonary function measurements.

Table 2

Participant demographics and pulmonary function measurements at follow-up visit.

3.2.

Imaging Measurements and Texture Features

Table 3 and Fig. S1 in the Supplementary Material summarize the texture feature definitions and descriptions used throughout the text. Table 4 summarizes quantitative MR imaging measurements and ranked texture features after feature selection step. Ex-smokers with accelerated lung function decline had significantly different local homogeneity normalized (Idmn, $p=0.048$) feature, large unventilated region emphasis (LGLZE, $p=0.01$) feature, small unventilated region emphasis (SRLGLE, $p=0.007$) feature, and local distribution of unventilated regions [small dependence low gray level emphasis (SDLGLE), $p<0.001$] feature. As summarized in Table 5, these subgroups did not have any statistically significant differences between imaging measurements at follow-up. Hyperpolarized gas MRI for representative participants in the stable and rapid decliner groups are shown in Fig. S1 in the Supplementary Material.

Table 3

MRI texture feature descriptors for machine-learning modeling.

Table 4

Imaging measurements of participants with stable and accelerated lung function decline.

Table 5

Imaging measurements of participants with stable and accelerated lung function decline at follow-up visit.

3.3.

Machine learning Modeling

As summarized in Table 6, the best performing machine learning model trained on demographic measurements (age, sex, BMI, % ${\mathrm{SpO}}_2$, pack years, and years since quit) achieved 64% prediction accuracy. Best performing spirometry model (FVC, TLC, ${\mathrm{FEV}}_1$, ${\%}_{\mathrm{pred}}$ IC, SVC, and RV/TLC) was cosine KNN algorithm with 68% accuracy, which was not statistically different from the models based on demographics measurements. The texture-based model achieved the highest sensitivity (86%) and an 81% accuracy via the ensemble RUSBoosted trees algorithm, exclusively trained on selected MR image texture features. The combined model achieved the highest accuracy of 82% using the Medium-Gaussian SVM algorithm, trained on FVC, age, longest distance across ventilated lung, local homogeneity, large unventilated region emphasis, and local distribution of unventilated regions.

Table 6

Machine-learning performance at predicting accelerated lung function decline.

As summarized in Table 7, the machine learning models trained exclusively on MRI texture features outperformed the machine learning models trained using participant demographic and spirometry measurements ($p<0.05$). In addition, the combined model trained on all available measurements also outperformed the demographic and spirometry-based models ($p<0.05$); however, the performance of the combined model failed to show a significant difference ($p=0.9$) to the machine learning models trained exclusively on MRI texture features.

Table 7

DeLong’s test for comparing the models for predicting accelerated disease progression in ex-smokers.

The ensemble models outperformed the single machine learning models, indicating the presence of more complex and non-linear relationships of texture features and accelerated lung function decline. Logistic regression models for predicting accelerated lung function decline were generated for individual clinical and imaging texture measurements with the receiver-operator characteristic curve AUC, which is summarized in Fig. 4. The best performing clinical measurements for predicting patients with accelerated ${\mathrm{FEV}}_1$ decline were FVC (AUC = 0.68) and TLC (AUC = 0.65). The overall best predictive measurement was wavelet-based local distribution of unventilated regions (SDLGLE, AUC = 0.77), which also outperformed standard imaging measurements, such as MRI VDP (AUC = 0.63).

Fig. 4

Receiver-operator characteristic curves of clinical and texture measurements. Top panel: logistic regression analysis of individual top-performing demographic and spirometry variables at predicting $\mathrm{\Delta }{\mathrm{FEV}}_1\ge 60\mathrm{mL}/\text{year}$ in ex-smoker participants. Bottom panel: logistic regression of selected top-performing MR imaging and texture features at predicting $\mathrm{\Delta }{\mathrm{FEV}}_1\ge 60\mathrm{mL}/\text{year}$. Individual texture features clearly outperformed established clinical variables available to physicians at predicting accelerated lung function decline. FVC, forced vital capacity; TLC, total lung capacity; ${\mathrm{FEV}}_1$, forced expiratory volume in 1 s; LL, low-low pass filter; SDLGLE, short distance low gray level emphasis; SRLGLE, short run low gray level emphasis; Idmn, inverse difference moment normalized; CDD1, cluster-defect diameter of one voxel; LGLZE, low gray level zone emphasis; and MAL, major axis length.

3.4.

Relationships with Clinical Measurements

Spearman correlations were used to evaluate the relationships between well-established clinical measurements and MRI texture features identified as significant predictors of clinically relevant ${\mathrm{FEV}}_1$ changes. As shown in Fig. 5, the best performing clinical measurements of FVC and TLC correlated with $\mathrm{\Delta }{\mathrm{FEV}}_1$ between visits ($\rho =-0.24$, $p=0.01$; $\rho =-0.23$, and $p=0.03$, respectively). Similarly, texture features from the original unfiltered image representing the number of small ventilation defects and the longest distance across ventilated lung correlated with $\mathrm{\Delta }{\mathrm{FEV}}_1$ ($\rho =-0.20$, $p=0.047$; $\rho =-0.21$, and $p=0.046$, respectively). The best performing wavelet-based local distribution of unventilated regions (SDLGLE) feature exhibited the strongest correlation with $\mathrm{\Delta }{\mathrm{FEV}}_1$ ($\rho =-0.29$, $p=0.006$), and only the longitudinal change in this specific texture correlated with the clinically relevant changes in ${\mathrm{FEV}}_1$ ($\rho =0.27$, $p=0.041$).

Fig. 5

Relationships between selected texture features and change in FEV1.Top panel: Spearman correlation for FVC and TLC with $\mathrm{\Delta }{\mathrm{FEV}}_1$ ($\rho =-0.24$, $p=0.01$; $\rho =-0.23$, and $p=0.03$, respectively) between baseline and follow-up visits. Middle panel: Spearman correlation for custom CDD1 feature and shape-major axis length with $\mathrm{\Delta }{\mathrm{FEV}}_1$ ($\rho =-0.20$, $p=0.047$; $\rho =-0.21$, and $p=0.046$, respectively) between baseline and follow-up visits. Bottom panel: Spearman correlation for wavelet-LL-filtered-GLDM-SDLGLE and $\mathrm{\Delta }\text{wavelet}$-LL-filtered-GLDM-SDLGLE with $\mathrm{\Delta }{\mathrm{FEV}}_1$ ($\rho =-0.29$, $p=0.006$; $\rho =0.27$, and $p=0.041$, respectively) between baseline and follow-up visits. ${\mathrm{FEV}}_1$, forced expiratory volume in 1 s; FVC, forced vital capacity; TLC, total lung capacity; CDD1, cluster-defect diameter of one voxel; LL, low-low pass filter; GLDM, gray level dependence matrix; and SDLGLE, short distance low gray level emphasis.