Quantitative CT analysis of honeycombing area predicts mortality in idiopathic pulmonary fibrosis with definite usual interstitial pneumonia pattern: A retrospective cohort study

Background Honeycombing on high-resolution computed tomography (HRCT) images is a key finding in idiopathic pulmonary fibrosis (IPF). In IPF, honeycombing area determined by quantitative CT analysis is correlated with pulmonary function test findings. We hypothesized that quantitative CT-derived honeycombing area (HA) might predict mortality in patients with IPF. Materials and methods Chest HRCT images of 52 IPF patients with definite usual interstitial pneumonia (UIP) pattern were retrospectively evaluated. Mortality data up to July 31, 2016, were recorded. Using a computer-aided system, HA and percentage of HA (%HA) were measured quantitatively. Predictors of 3-year mortality were evaluated using logistic regression models. Results The median %HA, %predicted forced vital capacity (FVC) and composite physiologic index (CPI) were 3.8%, 83.6%, and 33.6, respectively. According to GAP (gender, age, and physiology) stage, 20, 14, and 5 patients were classified under stages I-II-III, respectively. Percentage of HA was significantly correlated with %FVC, CPI, and GAP stage (all, p < 0.001). In univariate analysis, %HA, %FVC, and CPI were statistically significant predictors of mortality. In multivariate analysis using the stepwise regression method, only %HA (odds ratio [OR], 1.27; p = 0.011) was a significant independent predictors of mortality. Patients with %HA ≥ 4.8% had significantly lower survival rates than those with lesser %HA (median survival time, 1.3 vs 5.0 years; log-rank test; p < 0.001). Conclusion Quantitative CT-derived HA might be an important and independent predictor of mortality in IPF patients with definite UIP pattern.


Introduction
Idiopathic pulmonary fibrosis (IPF) is a specific form of chronic and progressive fibrosing interstitial pneumonia of unknown cause [1][2][3]. Honeycombing is one of the key findings on high-resolution computed tomography (HRCT) images in IPF with definite usual interstitial pneumonia (UIP) pattern. The extent of honeycombing area on HRCT images is an independent predictor of mortality in IPF [1,4,5]. Therefore, quantifying the extent of honeycombing on HRCT images allows better management of patients with IPF. Experts have evaluated the extent of honeycombing using CT visual scores [6][7][8][9]. However, this evaluation process is not uniform, yielding only fair to moderate agreement regarding honeycombing area among expert radiologists and pulmonologists [10]. The 2011 international statement also mentions the lack of uniformity in approach towards quantification of honeycombing area [4]. Recently, there are some reports about computer-based quantitative CT analyses [11][12][13][14][15][16][17][18]. These methods may be the useful and reliable for objective assessment of HRCT findings including honeycombing. However, the algorithms of their methods remain unclear, and not everyone can use these methods.
We recently developed a computer-aided system for quantitative CT analysis of honeycombing area in IPF [19]. This system has no interobserver errors and could, therefore, address the issues associated with current quantification methods for honeycombing area. Furthermore, everyone can use this system through the public-domain computer program, ImageJ.
Pulmonary function tests (PFTs) are a standard approach for objective monitoring of patients with IPF. Both forced vital capacity (FVC) and diffusing capacity of the lungs for carbon monoxide (DL CO ) are strongly correlated with disease extent determined histologically or estimated visually using CT images [1]. However, in patients with IPF, interpretation of PFT findings is confounded by coexisting emphysema [20,21]. Therefore, the composite physiologic index (CPI) was developed to improve on previous prognostic measures in IPF by adjusting for emphysema [20]. There exists another composite index, the GAP stage, which uses gender (G), age (A), and physiology (P) (FVC and DL CO ) variables for predicting mortality in IPF [22]. Although PFTs are very important in both these indexes, one should be aware that patients with severe IPF may not be able to perform PFTs.
The aim of this study was to assess the relationship between honeycombing area determined by quantitative CT analysis and mortality in patients with IPF. We hypothesized that quantitative CT-derived honeycombing area (HA) might predict survival in IPF with definite UIP pattern.

Study design and patients
All IPF patients with definite UIP pattern-diagnosed in accordance with the 2011 IPF guidelines-who visited the outpatient clinic of our hospital between April 2011 and March 2014 were enrolled in this study. Patients with secondary interstitial pneumonitis (such as collagen vascular disease, pneumonia caused by occupational or environmental exposure, chronic hypersensitivity pneumonitis, and drug-induced pneumonia) were excluded. Mortality and survival data up to July 31, 2016, were recorded. In order to evaluate 3-year mortality, patients with follow-up period < 3 years were excluded. For comparison of clinical features, patients with IPF were classified as "alive" or "deceased" on the basis of survival status at the end of the 3-year follow-up period. The study protocol conformed to the Declaration of Helsinki and was approved by the ethics committee of Shiga University of Medical Science (approval number 24-182; May 1, 2013), with a waiver of the need for informed consent because of the retrospective study design.
Chest HRCT images at the first visit were acquired using Toshiba Aquillion ONE (Toshiba Medical Systems Corp., Otawara, Tochigi, Japan) or SOMATOM Sensation Cardiac (Siemens Japan K.K., Shinagawa, Tokyo, Japan). The technical parameter for Toshiba Aquillion ONE was as follows: slice thickness 0.5-1.0 mm, x-ray voltage 120kVp, x-ray median current 200mA (range 80-350mA), rotation speed 0.5 sec, and reconstruction algorithm FC52. The technical parameter for SOMATOM Sensation Cardiac: slice thickness 0.5-1.0 mm, x-ray voltage 120kVp, x-ray median current 288mA (range 226-481mA), rotation speed 0.375 sec, and reconstruction algorithm B70f. The images used for the analysis were selected with 10-mm intervals. Among the enrolled patients, only PFT data acquired within 3 months of the selected CT images were analyzed. For these patients, CPI and GAP stage were determined as previously described [20,22]. The PFTs had been performed in accordance with the Japanese Respiratory Society guidelines [23].
A computer-aided method for quantitative CT analysis of honeycombing area was used to automatically measure honeycombing areas on whole HRCT images, as previously described [19]. Briefly, quantitative CT analysis of honeycombing area involved the following process: (1) low-attenuation areas enclosed by thick walls were extracted as honeycombing areas; (2) the sum of honeycombing areas in all slices was calculated as the HA. The values of HA and percentage of HA to total lung area (%HA) were measured for each patient (Fig 1).

Statistical analysis
Data are expressed as median and interquartile range (IQR). The correlation of HA with PFT findings and patient characteristics was evaluated by Spearman's correlation coefficient analysis.
For assessment of 3-year mortality, the clinical features of alive and deceased patients were compared using the Mann-Whitney U test. The effect of each variable on 3-year mortality was evaluated by univariate logistic regression analysis adjusted by age, sex, body mass index (BMI), and pack-years. The model with the greatest predictive utility was determined by multivariable logistic regression analysis adjusted by age, sex, BMI, and pack-years using the stepwise method. The cutoff value for %HA was determined on the basis of the area under the receiver operating characteristic curve and Youden's index [24]. Survival analysis was performed in accordance with the methods of Kaplan-Meier, with the endpoint being death. Differences in mortality were assessed by the log-rank test. All statistical analyses were performed using JMP version 9.0.2 (SAS Institute, Cary, NC). For all tests, p values < 0.05 were considered statistically significant.

Patient characteristics
Of the 72 consecutive IPF patients with definite UIP pattern who visited the outpatient clinic of our hospital between April 2011 and March 2014, one was excluded for the lack of HRCT data at the first visit; additionally, 19 patients with follow-up period < 3 years were excluded. Thus, 52 patients with IPF were finally enrolled in this study. Among the 52 enrolled patients, 45 had corresponding PFT data acquired within 3 months of the selected CT images. The demographic, clinical, and physiologic characteristics of these 52 patients are summarized in Table 1.

3-year mortality and baseline parameters
Of the 52 patients, 23 were deceased and 29 alive at the end of the 3-year follow-up period. Deceased group had the lower BMI and %FVC (p = 0.028, and 0.035, respectively; Table 3) and the higher GAP stage and %HA (p = 0.005, and 0.004, respectively; Table 3). Table 4 presents the results of univariate logistic regression analysis of correlation between baseline variables and 3-year mortality. Among the evaluated clinical, physiologic, and radiographic characteristics, %HA, %FVC, CPI, and GAP stage were statistically significant predictors of

Prediction of mortality by %HA
From the results of receiver operating characteristic analysis, the Youden's index and accuracy of %HA were estimated to be 4.8% and 73.1%, respectively (sensitivity, 86.2%; specificity, 56.5%; S1 Fig). The Kaplan-Meier survival curve demonstrated a significant difference in mortality, which were defined by a cutoff point of 4.8% for %HA (Fig 2). The median (IQR) survival times were 1.3 (0.8-2.8) and 5.0 (2.1-5.9) years, respectively, among patients with %HA greater and lesser than 4.8% (log-rank test; p < 0.001).

Discussion
The 2011 international statement mentions that extent of honeycombing on HRCT images is an independent predictor of mortality in IPF, although there is no uniformity in approach to quantitation [4]. We hypothesized that honeycombing area measured by quantitative CT analysis might be associated with an increased risk of mortality in IPF. The present results demonstrated that %HA was significantly correlated with 3-year mortality (Table 4). Furthermore, % HA was an independent predictor of 3-year mortality as the result of multivariate analysis using the stepwise regression method. These findings show that %HA might be a quantitative prognostic factor in IPF. Both FVC and DL CO are considered the gold standards for assessment of IPF. However, it is important to consider that, in patients with IPF, FVC and DL CO findings are confounded by coexisting emphysema [20,21]. Emphysema is present in approximately one-third of patients with IPF [25]. FVC measurements may not be appropriate for monitoring disease progression in patients with IPF and severe emphysema [26]. In patients with IPF with coexistent emphysema, CPI is thought to be a more accurate prognostic determinant than any individual parameter of PFTs [20]. Additionally, in IPF, GAP stage is thought to be a better prognostic determinant than any individual parameter of PFTs [22,27]. In the present study, %HA was significantly correlated with not only FVC and DL CO but also CPI and GAP stage ( Table 2). This is the retrospective cohort study of the consecutive IPF patients with HRCT scan data. Unfortunately, 7 patients did not have PFT within 3 months of the selected HRCT. It is inaccurate to discuss the relationship between selected HRCT and PFT performed at different time. Therefore, PFT of 7 patients was treated as missing data.
In order to identify better prognostic factors in IPF, we performed logistic regression analysis of correlation of 3-year mortality with %HA and other variables, including PFT findings. We selected the 3-year time point for mortality because several retrospective longitudinal studies have reported that the median survival time from the time of IPF diagnosis is between 2 and 3 years [4,[27][28][29][30][31]. Among the evaluated variables, only %HA was found to be significantly correlated with 3-year mortality upon multivariate logistic regression analysis using the stepwise method. This result showed that %HA was an independent and useful variable for prediction of 3-year mortality in the present study population. Among FVC, forced expiratory volume in 1 s (FEV 1 ), DL CO , CPI, GAP stage, and %HA, all variables exhibited strong correlations with each other. We, therefore, performed bivariate logistic regression analysis for the relationship of FVC, FEV 1 , DL CO , CPI, and GAP stage with %HA (S1 Table) and found that % HA was a significant predictor of mortality in combination with each of the variables except CPI and GAP stage. On the other hand, in combination with each other, neither %HA nor CPI or GAP stage was a significant predictor of mortality. These results seem to be derived from the strong correlation between CPI and %HA, GAP stage and %HA. These results show that in addition to CPI and GAP stage, %HA was a useful variable for prediction of 3-year mortality.
In evaluating the mortality in IPF, the analysis for time to death is also helpful. In univariate Cox proportional hazards models, %HA, %FVC, %FEV 1 , and CPI were statistically significant predictors of mortality (S2 Table). In multivariate Cox proportional hazards models including all variables using the stepwise regression method, only %HA (hazard ratio [HR], 1.13; 95% CI, 1.04-1.21; p = 0.006) was significant predictor of mortality. The results of bivariate Cox proportional hazards models (S3 Table) were almost similar to those of bivariate logistic regression analysis.
Because we wanted to evaluate honeycombing area by quantitative CT analysis, the present study population only included patients with IPF with definite usual interstitial pneumonia pattern on HRCT images. Although our study population was not representative of the entire spectrum of IPF, the relationship between mortality and GAP stage in our population was almost similar to that reported in previous studies [22,32] on IPF (S4 Table). In evaluation of prognosis, we believe that the present results might be comparable to the trend among the general IPF population.
According to the review of computer-based quantitative CT image analyses in IPF, there were 3 types of quantitative CT analyses; the density histogram analysis, the density mask analysis like our method, and texture classification analysis [33]. Including semi-quantitative analysis by visual scoring, each analysis has its merits that can be used in clinical settings. The texture classification analysis may be the most useful and reliable method for objective management of IPF. However, the specific CT index and the algorithm remain unclear [33]. It is desirable to establish a method which is transparent and which is freely available for widespread use, such as the low attenuation area (LAA) for chronic obstructive pulmonary disease (COPD) [34,35].
The present method of quantitative CT analysis of percentage of honeycombing area has important advantages relative to previously reported methods. First, since this method does not require PFT data (especially DL CO ), honeycombing area can be quantified in almost all cases of IPF. Second, this method does not involve interobserver or intraobserver errors and can be used by anybody as a public-domain computer program. Third, it allows quantitative assessment of prognosis in IPF. Our method using 10-mm interval images is not perfect. However, we think it is not inferior to the past methods using 3-4 CT images of fixed position such as aortic arch, carina, or right inferior pulmonary vein.
There are several limitations to the present study. First, this was a retrospective study with a small sample size, conducted at a single institution. Second, our patients did not receive uniform treatment for IPF. The rates of anti-fibrotic drugs usage were not different between the alive group and deceased group (Table 3). Third, we analyzed only IPF patients with definite UIP pattern to investigate whether the extent of honeycombing area was a prognostic predictor or not. In other words, we analyzed only the IPF patients with honeycombing area. We did not analyzed the patients without honeycombing area. A prospective multicenter study on quantitative analysis of honeycombing area is needed to confirm the usefulness of %HA in management of patients with IPF.

Conclusions
Relative extension of honeycombing in comparison to lung volume measured by quantitative CT analysis might be an important and independent predictor of mortality in IPF patients with definite UIP pattern. Furthermore, this analysis might provide new insights into disease progression in the future.
Supporting information S1