An analysis of the regional heterogeneity in tissue elasticity in lung cancer patients with COPD

Purpose: Recent advancements in obtaining image-based biomarkers from CT images have enabled lung function characterization, which could aid in lung interventional planning. However, the regional heterogeneity in these biomarkers has not been well documented, yet it is critical to several procedures for lung cancer and COPD. The purpose of this paper is to analyze the interlobar and intralobar heterogeneity of tissue elasticity and study their relationship with COPD severity.

Methods: We retrospectively analyzed a set of 23 lung cancer patients for this study, 14 of whom had COPD. For each patient, we employed a 5DCT scanning protocol to obtain end-exhalation and end-inhalation images and semi-automatically segmented the lobes. We calculated tissue elasticity using a biomechanical property estimation model. To obtain a measure of lobar elasticity, we calculated the mean of the voxel-wise elasticity values within each lobe. To analyze interlobar heterogeneity, we defined an index that represented the properties of the least elastic lobe as compared to the rest of the lobes, termed the Elasticity Heterogeneity Index (EHI). An index of 0 indicated total homogeneity, and higher indices indicated higher heterogeneity. Additionally, we measured intralobar heterogeneity by calculating the coefficient of variation of elasticity within each lobe.

Results: The mean EHI was 0.223 ± 0.183. The mean coefficient of variation of the elasticity distributions was 51.1% ± 16.6%. For mild COPD patients, the interlobar heterogeneity was low compared to the other categories. For moderate-to-severe COPD patients, the interlobar and intralobar heterogeneities were highest, showing significant differences from the other groups.

Conclusion: We observed a high level of lung tissue heterogeneity to occur between and within the lobes in all COPD severity cases, especially in moderate-to-severe cases. Heterogeneity results demonstrate the value of a regional, function-guided approach like elasticity for procedures such as surgical decision making and treatment planning.

1. Introduction

Respiratory diseases are a major cause of death, including Chronic Obstructive Pulmonary Disease (COPD), which leads to over 3 million deaths annually (6% worldwide), making it the third leading cause. More than 90% of these deaths occur in low-income or middle-income countries (1). Based on a meta-analysis conducted over studies from 1960 to 2010, COPD is associated with an increased risk of lung cancer with a hazard ratio of 2.22 (2). In addition to the heightened risk of cancer with COPD, one study showed that the cancer mortality rate is also greatly increased with increasing COPD severity, with the worst hazard ratio being 3.36 for GOLD stage 4 (3).

Since COPD is a notoriously heterogeneous disease (4–8), several pulmonary procedures for patients with COPD exploit regional lung function differences. One surgical treatment that has been developed for patients with severe COPD is lung volume reduction surgery (LVRS) (9). When a lobe experiences reduced ventilation due to air trapping, LVRS aims to collapse the lobe and release the trapped air, thus allowing the other more functional lobes to compensate (10, 11). An additional surgical technique building off LVRS is bronchoscopic lung volume reduction (BLVR), in which a one-way collapsible coil is placed in the airway (12, 13). LVRS, BLVR, and other lung interventions depend on disease heterogeneity to ensure that lung function can be compensated post-intervention (14, 15). Not only has heterogeneity been specifically identified as a factor in patient selection for these surgeries, but further research into patient selection, physiological testing, and disease characteristics have been called for in the literature (11, 13, 16).

In addition to surgical interventions, radiation therapy is another example of a treatment that could benefit from an understanding of the regional functionality of lung tissue, specifically when sparing organs at risk (OARs). A recent study showed that patients with lung cancer and COPD receive less curative treatments and experience higher mortality rates with a hazard ratio of 1.2 (17). In addition to poorer survival rates, normal tissue complications have also been observed in lung cancer patients with COPD comorbidities (18, 19). These increased risks of mortality and side effects call for a heightened attention to the regional effects of COPD in treatment planning. This information could aid functionally guided organs at risk (OARs) as opposed to the current lung-based contours (20–22). Studies using CT-based or PET-CT-based ventilation for functional tissue sparing have already shown reduction in doses to functional tissue (23–25) as well as reduction in grades 2+ and 3+ pneumonitis (24, 26).

To diagnose regional functionality, ventilation has been previously employed as a key biomarker. Evidence has been shown that ventilation mapping could aid in treatment planning by highlighting areas of high and low ventilation (27). While these studies of functional tissue mapping using ventilation are promising, accurately calculating ventilation from CT is still an ongoing effort (28–31). Many of the developed approaches incorporate transformation-based calculations, or calculations based on a mapping of CT images and comparison of densities, both of which rely on the accuracy of deformable image registration. It has been shown that the even small errors in the image registration can cause much larger errors in ventilation calculations (32, 33). For this reason, it may be beneficial to explore other functional properties to guide surgeries or radiotherapy until image registration can be performed with sufficient accuracy.

We propose elasticity as an additional functional property that can be measured from CT. Tissue elasticity is a biomechanical property that describes tissue stiffness (34). Elasticity can be calculated from medical images in several ways. In our work, it is done so with the anatomy from an end-exhalation image and the deformation vector field (DVF) mapping it to end-inhalation. These inputs are provided to an iterative model that estimates the Young’s modulus (YM) of each voxel to represent the elasticity, as has been validated in previous work (34, 35). Though also based on image registration, elasticity can be reliably calculated on a voxel-by-voxel basis since it is regularized by a physics-based model. Reduced elasticity is indicative of lung disease, so it can be another regional marker for COPD. Previous work showed that elasticity of voxels in the 1–3 kPa range was a better biomarker for COPD than the traditional RA950 (36). CT-based elasticity measurements thus offer an additional way to characterize regional lung function based on CT and a physics-based model.

In this study, we statistically characterized the regional lung function heterogeneity in patients with lung cancer and varying COPD comorbidities using tissue elasticity. Lobar elasticity distributions were calculated using our inverse biomechanical model. Interlobar heterogeneity was analyzed to provide information specifically relevant to surgical interventions. Additionally, intralobar heterogeneity was analyzed to provide insight to support regionally defined treatments such as radiotherapy. We related each heterogeneity measure to COPD severity. By providing evidence of tissue elasticity heterogeneity, we have shown evidence of how it could benefit surgical or radiotherapy planning for patients with all levels of COPD.

2. Materials and methods

We used a dynamic imaging protocol to obtain all patient data in this study. We acquired fast-helical free-breathing CTs (FHFBCT), constructed motion models using the 5DCT approach, and generated end-exhalation and end-inhalation images to serve as input for our tissue elasticity estimation model. Each step is outlined in Figure 1 with the relevant Materials and Methods sections labeled.

FIGURE 1

Figure 1. Flowchart of the study design beginning with image acquisition and resulting in heterogeneity analysis of tissue elasticity.

2.1. Patient data acquisition

We retrospectively employed a set of 23 lung cancer patients for this study. Each patient was identified as having no COPD, mild COPD, or moderate-to-severe COPD as defined by their physicians and based on spirometry. These severities were noted based on a retrospective review of their charts. Five patients had mild COPD and six had either moderate or severe COPD. Three patients were noted to have COPD, though the severity was unknown from their charts. These patients were only included in the analysis when examining all COPD patients.

For each patient, we collected a set of 25 FHFBCTs through an IRB approved study, as in previous investigations (37, 38). The scans were acquired in alternating directions with 120 kVp and 40 mAs. Three similar multi detector row CT scanners were used (Definition Flash, Biograph 64, Definition AS 64; Siemens Healthcare, Erlangen, Germany). Table 1 summarizes the scanning parameters. For all images, the field of view was set to 500 mm. The in-plane pixel resolution was 0.976 × 0.976 mm, and the slice thickness was 1.0 mm. After reconstruction, we resampled all images to 1 mm isotropic voxels.

TABLE 1

Table 1. Summary of scanners used to acquire patient data.

The FHFBCT acquisition process also included a simultaneous acquisition of breathing amplitude and flow rate (time-derivative of the amplitude) signals, which were needed for the 5DCT model construction to generate the end-exhalation and end-inhalation images. The breathing amplitude signal was obtained by using a pneumatic bellows (Lafayette Instrument Company, Lafayette, IN). The bellows was placed around the abdomen since its expansion was observed to provide the best correlation to the diaphragm motion (37). The bellows converted the change in pressure resulting from expansion into a voltage signal. The signal was sampled at 100 Hz, and amplitudes were assigned to each transverse slice. The bellows signals were finally synchronized and drift-corrected to account for measurement-related errors as previously demonstrated in O’Connell et al. (39).

2.2. 5DCT modeling

5DCT is a model-based CT approach that has been well-validated in previous works (38, 40–42). The 5DCT modeling process used in this study is briefly explained as follows.

The model generation process takes as input the 25 FHFBCT scans and the breathing signal amplitudes and flow rates. From the 25 FHFBCT images, we arbitrarily chose the first scan as the reference for the 5DCT model construction. Using an open-source deformable image registration software, deeds (43–45), the other 24 images were deformably registered to the reference scan as previously demonstrated (46, 47). We used the 24 DVFs with the breathing amplitude, $v$ , and rate, $f$ , to determine tissue-specific motion parameters, $\overrightarrow{\alpha }$ and $\overrightarrow{\beta }$ , by solving the relation shown below.

$\begin{array}{ll}\overrightarrow{X}=\overrightarrow{X_0}+\overrightarrow{\alpha }v+\overrightarrow{\beta }f& (1)\end{array}$

In Equation 1, $\overrightarrow{X_0}$ describes the tissue position at zero amplitude and flow, and $\overrightarrow{X}$ describes the tissue position at $v$ and $f$ . The inhalation motion is represented by the product of $\overrightarrow{\alpha }$ and the amplitude. Similarly, the hysteresis motion is represented by the product of $\overrightarrow{\beta }$ and the breathing rate.

To perform lobar HU-based and biomechanical property measurements, end-exhalation and end-inhalation images along with their corresponding DVFs needed to be generated. We selected the 5^th and 85^th percentile amplitudes with zero flow to represent end-exhalation and end-inhalation respiratory phases as shown in previous studies (42). Using Equation 1, the tissue-specific motion parameters were used to deform the reference image to its position in the end-exhalation and end-inhalation breathing phases.

2.3. Lobe segmentation

To obtain lobar elasticity distributions, we generated lobe masks that grouped lung voxels into one of the five lung lobes. We performed lobe segmentations semi-automatically on the FHFBCT reference scans using the open-source software Pulmonary Toolkit. The software first built a lobar approximation, then applied a “fissureness” filter, and finally fit a smooth multi-level B-spline curve through the fissureness and extrapolated to the lung boundaries to create the lobe segmentation (48). In some cases, the automated segmentation results experienced minor errors, so manual corrections were made using the graphical user interface in the Pulmonary Toolkit and verified by medical experts.

To only include lung parenchyma in the analysis, we removed blood vessels and tumors from the lobe masks by excluding voxels with greater than −700 HU. This threshold was chosen based on the HU distribution found in the lungs of clinical CT scans in a published study (49). Though this reference used inspiration CT rather than free-breathing CT, we experimented with higher thresholds and −700 HU was optimal for removing all vessels from the images. This is important because blood vessels do not expand or ventilate during respiration and have very high elasticity, so their values would misrepresent the lung parenchyma distributions. Therefore, our analysis only pertained to the parenchymal tissue with HU less than −700 in the 5th percentile images.

2.4. Tissue elasticity estimation

We used a previously developed and validated biomechanical model to estimate elasticity of the lung parenchyma, evaluated as the YM. The biomechanical model is based on changes in boundary constraints leading to corrective forces on a distribution of finite elements. These corrective forces are a summation of elastic, shear, and dashpot damping forces, which are included in Equations 2–4, respectively. In these equations, $YM$ is Young’s Modulus, $\Delta L_{ab}$ is the change in length between elements $a$ and $b$ , $L_{ab}$ is the resting length between $a$ and $b$ , $S_{ab}$ is the shear moduli (4 kPa), ${\mu }_{ab}$ is the local damping factor, and $\overrightarrow{v}$ is relative velocity.

$\begin{array}{ll}\overrightarrow{f_{E,ab}}={\displaystyle \sum }_b\left(YM\ast \frac{\Delta L_{ab}}{L_{ab}}\right)& (2)\end{array}$

$\begin{array}{ll}\overrightarrow{f_{2S,ab}}={\displaystyle \sum }_b\left(S_{ab}\ast \frac{L_{ab}-\Delta L_{ab}}{L_{ab}}\right)& (3)\end{array}$

$\begin{array}{ll}\overrightarrow{f_{v,ab}}={\displaystyle \sum }_b\left({\mu }_{ab}{}^{\ast }\left(\overrightarrow{v_b}-\overrightarrow{v_a}\right)\right)& (4)\end{array}$

In this study, the model used the DVF pointing from the 5DCT-based end-exhalation image to the end-inhalation image as the ground-truth DVF. Then, an initial elasticity distribution was set based on the HU of the end-exhalation scan. The elasticity distribution was then optimized to minimize the difference between the model-calculated DVF, based on boundary conditions and the calculated deformations from elastic forces, and the ground-truth DVF. In each iteration, the elasticity was updated, the new DVF was calculated from the updated elasticity values, and the DVF was compared to the ground truth. This process is shown in the flowchart in Figure 2. Further details of the finite element approach, governing equations, boundary conditions, and description of the inverse approach to optimize the elasticity can be found in several publications describing the model (34, 35, 50). The lobar elasticity was obtained by calculating the mean elasticity across each lobe distribution.

FIGURE 2

Figure 2. Flowchart of the iterative parameter optimization problem to estimate elasticity.

2.5. Heterogeneity analysis

First, we analyzed the mean elasticity of each lobe for patients in each COPD severity group. This was to investigate how well elasticity represented COPD severity in our cohort, as well as to observe the lobar trends in elasticity.

To measure interlobar heterogeneity, we defined the Elasticity Heterogeneity Index (EHI). According to a previous study, elasticities in the range of 1–3 kPa indicated diseased lung because COPD causes the lungs to poorly respond to expansion or contraction (36). Therefore, we based the EHI on the maximum percent of voxels in the COPD biomarker range among the five lobes, $E_{\max }$ , and the mean of the other four lobe percentages, $\overline{E_{~\max }}$ . This index indicated the least elastic lobe compared to the other four while emphasizing diseased tissue. We compared the mean EHI across COPD severity groups to study the interlobar heterogeneity in elasticity for these patients.

$\begin{array}{ll}EHI=\left|1-\frac{100-E_{\max }}{100-\overline{E_{~\max }}}\right|& (5)\end{array}$

To measure intralobar heterogeneity in elasticity, we calculated the coefficient of variation, or the standard deviation as a percent of the mean, of the elasticity distributions within each lobe. This metric summarized the spread of data within each lobe distribution. We examined histograms of the coefficient of variation across all lobes of patients in each COPD severity group to study differences in interlobar elasticity heterogeneity.

The coefficient of variation offers information about the dispersion of the data. To investigate the range of elasticity in each scenario, we also calculated an index of non-uniformity. This value has been defined slightly differently across the literature, but we have modified a definition from Jadhav, et al. because it offers a good sense of the range of any given metric (51). Our non-uniformity index is defined in Equation 6, where $NU$ is the non-uniformity index, $E_{{95}^{th}}$ is the 95th percentile elasticity, and $\overline{E}$ is the mean elasticity.

$\begin{array}{ll}NU=\frac{E_{{95}^{th}}-\overline{E}}{\overline{E}}& (6)\end{array}$

We performed two sample t-tests to test the differences between different COPD severity groups for each previously mentioned parameter. We used an F-test for variance equality to determine if the t-tests should be tested with or without equal variances. We tested the differences between these groups in elasticity, EHI, and intralobar coefficient of variation at the 5% significance level.

2.6. Tumor presence

Since all of these patients had lung cancer, there was always at least one lobe that contained a tumor. To explore the effect of the tumor on our analysis, we first separated the lobe(s) with the tumor (five of the patients had two lobes containing tumors) from the other non-cancerous lobes for each patient and calculated the mean elasticity of each group of lobes. We calculated the mean of these differences across patients as well as the mean difference relative to the mean elasticity of all five lobes. We also performed a two sample t-test for each patient to see if the lobe with the tumor had a significantly different elasticity than the other lobes.

However, lobe-dependent trends in elasticity may mask this effect and render interpretation of these results difficult. Therefore, we also calculated how often the lobe with the highest elasticity was also the lobe containing the tumor. We also tested across patients while keeping the lobe constant to see if the group with the tumors had a significantly different elasticity than the group without the tumors. For example, testing if the right upper lobes with tumors were significantly different than the right upper lobes without tumors.

3. Results

3.1. Elasticity

Figure 3 shows examples of our results for two patients. Figures 3A,B show coronal slices of the end-exhalation and end-inhalation scans, respectively, for the patient with high interlobar heterogeneity (EHI = 0.560). The pronounced elasticity difference between the left upper and left lower lobes is apparent in the elasticity distribution shown in Figure 3C. Figures 3D,E show the end-exhalation and end-inhalation scans for a patient with lower interlobar heterogeneity (EHI = 0.248). Regions of high and low elasticity are more dispersed throughout the lobes in this case.

FIGURE 3

Figure 3. Example distributions for two patients (A) End-exhalation scan, (B) end-inhalation scan, and (C) elasticity distribution for patient with high interlobar heterogeneity. (D) End-exhalation scan, (E) end-inhalation scan, and (F) elasticity distribution for patient with low interlobar heterogeneity.

Figure 4 shows a bar graph of the mean elasticity across patients in each COPD severity group separated into lobes, with error bars denoting the standard deviations. It is notable that elasticity was greatest in the lower lobes, which is reflective of their larger deformations during inhalation. This figure not only shows the typical elasticity differences between lobes, but also shows that elasticity decreased for patients with moderate-to-severe COPD, though the error bars still overlapped. When comparing the mean elasticity of the entire lungs across severity groups, significant differences were seen between patients with moderate-to-severe COPD and patients without COPD (p < 0.01), patients with mild COPD (p < 0.01), and all patients (p < 0.01). Significant differences were seen between patients with COPD and without COPD as well (p = 0.02), which was most likely heavily weighted by the patients with moderate-to-severe COPD. No statistically significant differences were seen between patients with mild COPD and patients without COPD. All comparisons were found to have equal variances as a result of F-tests. Overall, we found that the effect of COPD on tissue elasticity was most prominent in moderate-to-severe cases.

FIGURE 4

Figure 4. Bar graph of the mean lobar elasticity of patients in each COPD severity group. Error bars indicate standard deviations.

3.2. Heterogeneity

The mean EHI was 0.388 ± 0.162. Therefore, on average, there were moderate levels of interlobar heterogeneity with a relatively high level of variation across the patient cohort. For patients with no COPD, mild COPD, and moderate-to-severe COPD, the mean EHIs were 0.385 ± 0.183, 0.258 ± 0.133, and 0.473 ± 0.128, respectively. The greatest degree of interlobar heterogeneity was seen in patients with moderate-to-severe COPD. All of these patients had EHIs of at least 0.2. Additionally, among patients with mild COPD, EHIs varied greatly. The mean EHI in patients with mild COPD was lower than those without COPD, though the difference was not statistically significant (p = 0.20). However, the difference in EHI between patients with mild COPD and patients with moderate-to-severe COPD was statistically significant (p = 0.02). All comparisons were found to have equal variances as a result of F-tests.

The coefficients of variation of each lobar elasticity distribution are summarized by the histograms in Figure 5. The histogram in Figure 5A includes all lobes to offer a general sense of intralobar heterogeneity regardless of disease state. The mean of this distribution was 51.1% ± 16.6%, showing that there was a high level of intralobar heterogeneity in general. Figure 5B shows the histogram of the coefficient of variation for the lobes of patients without COPD (mean 47.9% ± 17.5%), and Figure 5C shows the histogram for lobes of patients with mild COPD (mean 47.3% ± 11.2%). These two groups exhibited a similar distribution to each other and to the distribution of all patient lobes. On the other hand, the histogram of the moderate-to-severe COPD group coefficients shown in Figure 5D exhibited a higher distribution (mean 58.9% ± 17.7%). The only significant differences were between the coefficient of variation of lobes with moderate-to-severe COPD and each of the other groups (p = 0.01 when comparing to lobes without COPD, p < 0.01 when comparing to lobes with mild COPD, p = 0.03 when comparing to all lobes). All comparisons to patients with mild COPD were tested with unequal variances as determined by F-test results, indicating that patients with mild COPD may have much more variation in their intralobar heterogeneity. All other comparisons were tested with equal variances. This indicates that patients with more severe COPD may experience higher levels of intralobar heterogeneity, which may be potentially useful as a metric for function-preserving interventions.

FIGURE 5

Figure 5. Relative histograms of the coefficients of variation of elasticity distributions for (A) lobes of all patients, (B) lobes of patients without COPD, (C) lobes of patients with mild COPD, and (D) lobes of patients with moderate-to-severe COPD.

Figure 6 shows violin plots of the non-uniformity index across the five lobes for all patients as well as each COPD severity group. Figure 6A offers a sense of the non-uniformity across the entire patient cohort. The mean value was about 1 for each group and lobe, which shows that the highest elasticity was most often twice the mean value. Figure 6 also shows that in all COPD severity groups, the index of non-uniformity was slightly less in the lower lobes. In some cases, the non-uniformity index was as high as 2.5, showing that the elasticity could be 3.5 times the mean in some regions of the lungs. These results show a high degree of heterogeneity expressed as the range of elasticity among all patients.

FIGURE 6

Figure 6. Violin plots of the index of non-uniformity for (A) lobes of all patients, (B) lobes of patients without COPD, (C) lobes of patients with mild COPD, and (D) lobes of patients with moderate-to-severe COPD.

3.3. Tumor presence

Comparing the elasticity of the lobe(s) with the tumor to the mean of the elasticity of the non-cancerous lobes, we found that the mean of the absolute differences between these values across all patients was 1.011 ± 0.470 kPa. The mean of the elasticity differences relative to the mean elasticity of all five lobes was 27.5% ± 15.2%. A two sample t-test to compare the lobe(s) with the tumor against the other lobes showed that two patients had significant differences (both p < 0.01).

4. Discussion

In this paper, we analyzed the lung tissue functional heterogeneity on an interlobar and intralobar basis that could be useful for decision making in pulmonary surgeries or planning in radiation therapy. We constructed 5DCT models to obtain end-exhalation and end-inhalation scans that served as the basis for all calculations. We used the scans and the deformation vectors mapping them to estimate elasticity using a validated biomechanical model. Heterogeneity was analyzed across lobes as well as within lobes. Since heterogeneity was consistently observed in elasticity at both levels and was dependent on COPD in more severe cases, we conclude that the heterogeneity of lung function warrants further investigation to improve surgical decision making and radiotherapy planning.

Elasticity offers a biomechanical property to indicate lung tissue functionality. Though the biomechanical model used in this study takes a DVF from image registration as input, the iterative process calculated a model-based DVF as elasticity is updated. This iteration continues until the DVF converges to the registration-based DVF. However, since the model-based DVF is calculated from the governing equations (Equations 2–4), small, nonphysical errors in the original DVF will be regularized by the iterative approach. This is an advantage over gradient-based ventilation techniques because detectable registration errors will be less likely to impact elasticity. However, a holistic approach to characterizing lung function should include elasticity and ventilation as complements to offer a more comprehensive picture of lung function. Our future work will include a focus on developing a reliable and validated ventilation calculation technique. Once this is realized, we will perform a thorough comparative study between elasticity and ventilation in terms of heterogeneity, impact of COPD, and how to interpret the different physiology represented by each. Once registration techniques improve in accuracy and robustness, and a dedicated comparative analysis is performed, the combination of these metrics will make for a very powerful tool in medicine.

The 5DCT approach used in this study has been well-validated in accurately modeling tissue motion. However, the accuracy of the modeling still has two primary limitations. One is motion blur. Though mitigated by the fast-scanning protocol, some motion blur artifacts may persist. The second limitation is from image registration. The registration technique used has been well-validated and shown to produce very accurate results. Based on a TG-132-based (52) analysis of the registrations used in this study, we calculated the target registration error of 50 manually defined anatomical landmarks per patient to be 1.31 ± 0.87 mm on average. In some instances, some inaccuracies were still found in the inferior lungs possibly due to blur in the images as well as the need to register larger deformations. In the future, the motion blur correction and improvements in image registration techniques should ameliorate these issues, respectively. Finally, more detailed model terms such as a cardiac motion term could be included in future developments to fine tune model accuracy.

Lobe segmentation was performed using open-source software with published validation (48). However, in certain cases, minor corrections were required possibly due to fissure incompleteness. To reduce segmentation time and increase accuracy, a recently published machine learning technique (53) will be implemented in future work. Additionally, the vessels were segmented using a threshold technique. There may be noise in the HU and boundary voxels that could cause some parenchymal voxels to be mistakenly segmented out of the images. Therefore, in the future, we will incorporate blood vessel tracking algorithms to increase the sophistication of the technique and maintain as many parenchymal voxels as possible.

The tissue elasticity estimation using YM has been well-validated in previous studies (34). However, inaccuracies in the images or DVFs could potentially limit the model. Future work will include using a machine learning approach to estimate elasticity from just the end-exhalation scans, which would help to limit the effects of DVF inaccuracies (50). In the future, this will replace the current iterative approach for faster and more accurate results. This technique could greatly reduce the computation time (from several days per patient to the order of a few seconds). This would enable faster data collection and processing as well as a clinical path to relevance for tissue elasticity estimation.

Limitations of the study results are largely derived from the retrospective data collection. For example, patient information like spirometry-derived function data and smoking history could have provided a more holistic explanation for tissue elasticity in these patients, especially patients with no COPD but an extensive smoking history. In the future, with a prospective protocol, we will acquire consistent spirometry function measurements before or after 5DCT acquisition to compare our results to a more quantitative assessment of COPD as well as carefully documented smoking history. Additionally, our sample size is a major limitation in this study. With small numbers of patients in each COPD severity group, our conclusions are interesting, but mostly hypothesis generating, and we will require a larger cohort in our next studies to continue investigating the relationship between COPD and elasticity heterogeneity. Furthermore, our statistical testing is limited as well by the small cohort size. However, if we use a Mann–Whitney U-test instead of the t-test, which offers a nonparametric alternative with fewer assumptions about the distribution (54), our findings remain the same. Even still, future studies will require additional patients in each COPD severity group.

In this study, we chose to investigate elasticity as it related to physician-determined COPD status. Parametric response mapping is a quantitative alternative to assessing COPD status that categorizes voxels based on thresholds in deep inspiration or expiration scans (55, 56). However, in our study, we used 5DCT-based deformations calculated from free-breathing CT and thus did not have the required adjusted HUs to perform this analysis. Moreover, new parametric response mapping category definitions would need to be determined to assess quiet respiration scans. Another planned development of our technique is to work towards incorporating these approaches to capture a more detailed comparison of voxel-specific disease to elasticity.

In conclusion, heterogeneity of tissue elasticity was consistently observed within and across lobes in patients with no COPD, mild COPD, and moderate-to-severe COPD. Increased heterogeneity was observed with patients with moderate-to-severe COPD. Therefore, elasticity measurement on a lobar or sublobar basis could enable the best guidance for decision making during function sparing treatment planning.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics statement

The studies involving human participants were reviewed and approved by Institutional Review Board IRB #19-002033. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.

Author contributions

ML, DL, and AS: research design. DL, DO’C, and AR: 5DCT development and data collection. MM-G: general study design and radiology context. AS and BS: tissue elasticity development and data processing. JG and IB: pulmonology expertise and guidance during heterogeneity analysis. ML and BS: data processing and analysis. ML, BS, and LN: drafting and critically revising the paper. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the Tobacco Related Disease Research Program 27IR-0056, NIH R56 1R56HL139767-01A1, Ken and Wendy Ruby Foundation, and the UCLA Department of Radiation Oncology.

Conflict of interest

DL is the founder of Pulmonum, LLC. The Department of Radiological Sciences at UCLA has a Master Research Agreement with Siemens Healthineers; MM-G is a Grant recipient of Siemens Healthineers.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. World Health Organization . Chronic obstructive pulmonary disease (COPD). (2021). Available at: https://www.who.int/news-room/fact-sheets/detail/chronic-obstructive-pulmonary-disease-(copd) (Accessed 21 June 2021).

Google Scholar

2. Brenner, DR , McLaughlin, JR , and Hung, RJ . Previous lung diseases and lung cancer risk: a systematic review and meta-analysis. PLoS One. (2011) 6:e17479. doi: 10.1371/journal.pone.0017479

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Mattila, T , Vasankari, T , Kanervisto, M , Laitinen, T , Impivaara, O , Rissanen, H, et al. Association between all-cause and cause-specific mortality and the GOLD stages 1–4: a 30-year follow-up among Finnish adults. Respir Med. (2015) 109:1012–8. doi: 10.1016/j.rmed.2015.06.002

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Valipour, A , Shah, PL , Gesierich, W , Eberhardt, R , Snell, G , Strange, C, et al. Patterns of emphysema heterogeneity. Respiration. (2015) 90:402–11. doi: 10.1159/000439544

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Labaki, WW , Martinez, CH , Martinez, FJ , Galbán, CJ , Ross, BD , Washko, GR, et al. The role of chest computed tomography in the evaluation and management of the patient with chronic obstructive pulmonary disease. Am J Respir Crit Care Med. (2017) 196:1372–9. doi: 10.1164/rccm.201703-0451PP

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Cho, MH , Castaldi, PJ , Hersh, CP , Hobbs, BD , Barr, RG , Tal-Singer, R, et al. A genome-wide association study of emphysema and airway quantitative imaging phenotypes. Am J Respir Crit Care Med. (2015) 192:559–69. doi: 10.1164/rccm.201501-0148OC

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Barrecheguren, M , and Miravitlles, M . COPD heterogeneity: implications for management. BioMed Central. (2016) 11:1–2. doi: 10.1186/s40248-016-0053-4

CrossRef Full Text | Google Scholar

8. Ju, J , Li, R , Gu, S , Leader, JK , Wang, X , Chen, Y, et al. Impact of emphysema heterogeneity on pulmonary function. PLoS One. (2014) 9:e113320. doi: 10.1371/journal.pone.0113320

PubMed Abstract | CrossRef Full Text | Google Scholar

9. DeCamp, MM Jr, McKenna, RJ Jr, Deschamps, CC , and Krasna, MJ . Lung volume reduction surgery: technique, operative mortality, and morbidity. Proc Am Thorac Soc. (2008) 5:442–6. doi: 10.1513/pats.200803-023ET

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Cooper, J , Trulock, E , Triantafillou, A , Patterson, G , Pohl, M , Deloney, P, et al. Bilateral pneumectomy (volume reduction) for chronic obstructive pulmonary disease. J Thorac Cardiovasc Surg. (1995) 109:106–19. doi: 10.1016/S0022-5223(95)70426-4

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Weinmann, GG , and Hyatt, R . Evaluation and research in lung volume reduction surgery. Am J Respir Crit Care Med. (1996) 154:1913–8. doi: 10.1164/ajrccm.154.6.8970386

CrossRef Full Text | Google Scholar

12. Criner, GJ , Sue, R , Wright, S , Dransfield, M , Rivas-Perez, H , Wiese, T, et al. A multicenter randomized controlled trial of Zephyr endobronchial valve treatment in heterogeneous emphysema (LIBERATE). Am J Respir Crit Care Med. (2018) 198:1151–64. doi: 10.1164/rccm.201803-0590OC

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Labarca, G , Uribe, JP , Pacheco, C , Folch, E , Kheir, F , Majid, A, et al. Bronchoscopic lung volume reduction with endobronchial zephyr valves for severe emphysema: a systematic review and meta-analysis. Respiration. (2019) 98:268–78. doi: 10.1159/000499508

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Russi, E , Bloch, K , and Weder, W . Functional and morphological heterogeneity of emphysema and its implication for selection of patients for lung volume reduction surgery. Eur Respir J. (1999) 14:230–6. doi: 10.1034/j.1399-3003.1999.14a39.x

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Slebos, D-J , Klooster, K , Ernst, A , Herth, FJ , and Kerstjens, HA . Bronchoscopic lung volume reduction coil treatment of patients with severe heterogeneous emphysema. Chest. (2012) 142:574–82. doi: 10.1378/chest.11-0730

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Shah, PL , Herth, FJ , van Geffen, WH , Deslee, G , and Slebos, D-J . Lung volume reduction for emphysema. Lancet Respir Med. (2017) 5:147–56. doi: 10.1016/S2213-2600(16)30221-1

CrossRef Full Text | Google Scholar

17. Goffin, JR , Corriveau, S , Tang, GH , and Pond, GR . Management and outcomes of patients with chronic obstructive lung disease and lung cancer in a public healthcare system. PLoS One. (2021) 16:e0251886. doi: 10.1371/journal.pone.0251886

CrossRef Full Text | Google Scholar

18. Kimura, T , Matsuura, K , Murakami, Y , Hashimoto, Y , Kenjo, M , Kaneyasu, Y, et al. CT appearance of radiation injury of the lung and clinical symptoms after stereotactic body radiation therapy (SBRT) for lung cancers: are patients with pulmonary emphysema also candidates for SBRT for lung cancers? International journal of radiation oncology* biology*. Physics. (2006) 66:483–91. doi: 10.1016/j.ijrobp.2006.05.008

CrossRef Full Text | Google Scholar

19. Rancati, T , Ceresoli, GL , Gagliardi, G , Schipani, S , and Cattaneo, GM . Factors predicting radiation pneumonitis in lung cancer patients: a retrospective study. Radiother Oncol. (2003) 67:275–83. doi: 10.1016/S0167-8140(03)00119-1

CrossRef Full Text | Google Scholar

20. Kong, FM , Ritter, T , Quint, DJ , Senan, S , Gaspar, LE , Komaki, RU, et al. Consideration of dose limits for organs at risk of thoracic radiotherapy: atlas for lung, proximal bronchial tree, esophagus, spinal cord, ribs, and brachial plexus. International journal of radiation oncology* biology*. Physics. (2011) 81:1442–57. doi: 10.1016/j.ijrobp.2010.07.1977

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Dolz, J , Kirişli, H , Fechter, T , Karnitzki, S , Oehlke, O , Nestle, U, et al. Interactive contour delineation of organs at risk in radiotherapy: clinical evaluation on NSCLC patients. Med Phys. (2016) 43:2569–80. doi: 10.1118/1.4947484

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Giżyńska, MK , Rossi, L , den Toom, W , Milder, MT , de Vries, KC , Nuyttens, J, et al. Largely reduced OAR doses, and planning and delivery times for challenging robotic SBRT cases, obtained with a novel optimizer. J Appl Clin Med Phys. (2021) 22:35–47. doi: 10.1002/acm2.13172

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Yamamoto, T , Kabus, S , Von Berg, J , Lorenz, C , and Keall, PJ . Impact of four-dimensional computed tomography pulmonary ventilation imaging-based functional avoidance for lung cancer radiotherapy. International journal of radiation oncology* biology*. Physics. (2011) 79:279–88. doi: 10.1016/j.ijrobp.2010.02.008

CrossRef Full Text | Google Scholar

24. Vinogradskiy, Y , Rusthoven, CG , Schubert, L , Jones, B , Faught, A , Castillo, R, et al. Interim analysis of a two-institution, prospective clinical trial of 4DCT-ventilation-based functional avoidance radiation therapy. International journal of radiation oncology* biology*. Physics. (2018) 102:1357–65. doi: 10.1016/j.ijrobp.2018.07.186

CrossRef Full Text | Google Scholar

25. Siva, S , Thomas, R , Callahan, J , Hardcastle, N , Pham, D , Kron, T, et al. High-resolution pulmonary ventilation and perfusion PET/CT allows for functionally adapted intensity modulated radiotherapy in lung cancer. Radiother Oncol. (2015) 115:157–62. doi: 10.1016/j.radonc.2015.04.013

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Faught, AM , Miyasaka, Y , Kadoya, N , Castillo, R , Castillo, E , Vinogradskiy, Y, et al. Evaluating the toxicity reduction with computed tomographic ventilation functional avoidance radiation therapy. International journal of radiation oncology* biology*. Physics. (2017) 99:325–33. doi: 10.1016/j.ijrobp.2017.04.024

CrossRef Full Text | Google Scholar

27. Patton, TJ , Gerard, SE , Shao, W , Christensen, GE , Reinhardt, JM , and Bayouth, JE . Quantifying ventilation change due to radiation therapy using 4 DCT Jacobian calculations. Med Phys. (2018) 45:4483–92. doi: 10.1002/mp.13105

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Scharm, SC , Vogel-Claussen, J , Schaefer-Prokop, C , Dettmer, S , Knudsen, L , Jonigk, D, et al. Quantification of dual-energy CT-derived functional parameters as potential imaging markers for progression of idiopathic pulmonary fibrosis. Eur Radiol. (2021) 31:6640–51. doi: 10.1007/s00330-021-07798-w

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Feng, A , Shao, Y , Wang, H , Chen, H , Gu, H , Duan, Y, et al. A novel lung-avoidance planning strategy based on 4DCT ventilation imaging and CT density characteristics for stage III non-small-cell lung cancer patients. Strahlenther Onkol. (2021) 197:1084–92. doi: 10.1007/s00066-021-01821-1

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Low, DA , O'Connell, D , Lauria, M , Stiehl, B , Naumann, L , Lee, P, et al. Ventilation measurements using fast-helical free-breathing CT. Med Phys. (2021) 48:6094–105. doi: 10.1002/mp.15173

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Kipritidis, J , Tahir, BA , Cazoulat, G , Hofman, MS , Siva, S , Callahan, J, et al. The VAMPIRE challenge: a multi-institutional validation study of CT ventilation imaging. Med Phys. (2019) 46:1198–217. doi: 10.1002/mp.13346

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Castillo, E , Castillo, R , Vinogradskiy, Y , Dougherty, M , Solis, D , Myziuk, N, et al. Robust CT ventilation from the integral formulation of the Jacobian. Med Phys. (2019) 46:2115–25. doi: 10.1002/mp.13453

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Castillo, E , Castillo, R , Vinogradskiy, Y , and Guerrero, T . The numerical stability of transformation-based CT ventilation. Int J Comput Assist Radiol Surg. (2017) 12:569–80. doi: 10.1007/s11548-016-1509-x

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Hasse, K , O'Connell, D , Min, Y , Neylon, J , Low, DA , and Santhanam, A . Estimation and validation of patient-specific high-resolution lung elasticity derived from 4DCT. Med Phys. (2018) 45:666–77. doi: 10.1002/mp.12697

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Hasse, K , Neylon, J , and Santhanam, AP . Feasibility and quantitative analysis of a biomechanical model-guided lung elastography for radiotherapy. Biomed Phys Engineer Exp. (2017) 3:025006. doi: 10.1088/2057-1976/aa5d1c

CrossRef Full Text | Google Scholar

36. Hasse, K , Neylon, J , Min, Y , O'Connell, D , Lee, P , Low, DA, et al. Feasibility of deriving a novel imaging biomarker based on patient-specific lung elasticity for characterizing the degree of COPD in lung SBRT patients. Br J Radiol. (2019) 92:20180296. doi: 10.1259/bjr.20180296

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Lauria, M , Navaratna, R , O’Connell, D , Santhanam, A , Lee, P , and Low, DA . Investigating internal–external motion correlation using fast helical CT. Med Phys. (2021) 48:1823–31. doi: 10.1002/mp.14759

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Thomas, D , Lamb, J , White, B , Jani, S , Gaudio, S , Lee, P, et al. A novel fast helical 4D-CT acquisition technique to generate low-noise sorting artifact–free images at user-selected breathing phases. International journal of radiation oncology* biology*. Physics. (2014) 89:191–8. doi: 10.1016/j.ijrobp.2014.01.016

CrossRef Full Text | Google Scholar

39. O’Connell, D , Thomas, DH , Lamb, JM , Lewis, JH , Dou, T , Sieren, JP, et al. Dependence of subject-specific parameters for a fast helical CT respiratory motion model on breathing rate: an animal study. Phys Med Biol. (2018) 63:15. doi: 10.1088/1361-6560/aaaa15

CrossRef Full Text | Google Scholar

40. Low, DA , Parikh, PJ , Lu, W , Dempsey, JF , Wahab, SH , Hubenschmidt, JP, et al. Novel breathing motion model for radiotherapy. International journal of radiation oncology* biology*. Physics. (2005) 63:921–9. doi: 10.1016/j.ijrobp.2005.03.070

CrossRef Full Text | Google Scholar

41. O'Connell, DP , Thomas, DH , Dou, TH , Lamb, JM , Feingold, F , Low, DA, et al. Comparison of breathing gated CT images generated using a 5DCT technique and a commercial clinical protocol in a porcine model. Med Phys. (2015) 42:4033–42. doi: 10.1118/1.4922201

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Low, DA , White, BM , Lee, PP , Thomas, DH , Gaudio, S , Jani, SS, et al. A novel CT acquisition and analysis technique for breathing motion modeling. Phys Med Biol. (2013) 58:L31–6. doi: 10.1088/0031-9155/58/11/L31

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Heinrich, MP , Jenkinson, M , Bhushan, M , Matin, T , Gleeson, FV , Brady, M, et al. MIND: modality independent neighbourhood descriptor for multi-modal deformable registration. Med Image Anal. (2012) 16:1423–35. doi: 10.1016/j.media.2012.05.008

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Heinrich, MP , Jenkinson, M , Brady, M , and Schnabel, JA . MRF-based deformable registration and ventilation estimation of lung CT. IEEE Trans Med Imaging. (2013) 32:1239–48. doi: 10.1109/TMI.2013.2246577

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Heinrich, MP , Jenkinson, M , Papież, BW , Brady, M , and Schnabel, JA . Towards realtime multimodal fusion for image-guided interventions using self-similarities. In: International conference on medical image computing and computer-assisted intervention; (2013): Springer.

Google Scholar

46. Dou, TH , Thomas, DH , O'Connell, DP , Lamb, JM , Lee, P , and Low, DA . A method for assessing ground-truth accuracy of the 5DCT technique. International journal of radiation oncology* biology*. Physics. (2015) 93:925–33. doi: 10.1016/j.ijrobp.2015.07.2272

CrossRef Full Text | Google Scholar

47. O’Connell, D , Thomas, D , Dou, T , Aliotta, E , Lewis, J , Lamb, J, et al. Adaptive weighted median filtering for reduced blurring when fusing co-registered fast helical CT images. Biomed Phys Engineer Exp. (2017) 3:067002. doi: 10.1088/2057-1976/aa889d

CrossRef Full Text | Google Scholar

48. Doel, T , Matin, TN , Gleeson, FV , Gavaghan, DJ , and Grau, V . Pulmonary lobe segmentation from CT images using fissureness, airways, vessels and multilevel B-splines. In: 2012 9th IEEE International Symposium on Biomedical Imaging (ISBI); (2012): IEEE.

Google Scholar

49. Abadi, E , Sanders, J , and Samei, E . Patient-specific quantification of image quality: an automated technique for measuring the distribution of organ Hounsfield units in clinical chest CT images. Med Phys. (2017) 44:4736–46. doi: 10.1002/mp.12438

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Santhanam, AP , Stiehl, B , Lauria, M , Hasse, K , Barjaktarevic, I , Goldin, J, et al. An adversarial machine learning framework and biomechanical model-guided approach for computing 3D lung tissue elasticity from end-expiration 3DCT. Med Phys. (2021) 48:667–75. doi: 10.1002/mp.14252

CrossRef Full Text | Google Scholar

51. Jadhav, S , Tekale, R , and Date, PP , editors. Analysis of strain non-uniformity index (SNI) for DifferentGeometries of drawn sheet metal parts. In: Journal of physics: Conference series; (2016): IOP Publishing.

Google Scholar

52. Brock, KK , Mutic, S , McNutt, TR , Li, H , and Kessler, ML . Use of image registration and fusion algorithms and techniques in radiotherapy: report of the AAPM radiation therapy committee task group no. 132. Med Phys. (2017) 44:e43–76. doi: 10.1002/mp.12256

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Stiehl, B , Lauria, M , Singhrao, K , Goldin, J , Barjaktarevic, I , Low, D, et al. Scalable quorum-based deep neural networks with adversarial learning for automated lung lobe segmentation in fast helical free-breathing CTs. Int J Comput Assist Radiol Surg. (2021) 16:1775–84. doi: 10.1007/s11548-021-02454-6

CrossRef Full Text | Google Scholar

54. McKnight, PE , and Najab, J . Mann-Whitney U test. Corsini Encycloped Psychol. (2010) 1–1. doi: 10.1002/9780470479216.corpsy0524

CrossRef Full Text | Google Scholar

55. Boes, JL , Hoff, BA , Bule, M , Johnson, TD , Rehemtulla, A , Chamberlain, R, et al. Parametric response mapping monitors temporal changes on lung CT scans in the subpopulations and intermediate outcome measures in COPD study (SPIROMICS). Acad Radiol. (2015) 22:186–94. doi: 10.1016/j.acra.2014.08.015

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Pompe, E , van Rikxoort, EM , Schmidt, M , Rühaak, J , Estrella, LG , Vliegenthart, R, et al. Parametric response mapping adds value to current computed tomography biomarkers in diagnosing chronic obstructive pulmonary disease. Am J Respir Crit Care Med. (2015) 191:1084–6. doi: 10.1164/rccm.201411-2105LE

PubMed Abstract | CrossRef Full Text | Google Scholar