Quantifying ventilation change due to radiation therapy using 4DCTJacobian calculations

1.A.4DCT‐derived ventilation

Ventilation, along with blood perfusion and gas exchange, are necessary components of lung function. Four‐dimensional computed tomography (4DCT) can be used to compute a surrogate for regional ventilation of lung tissue and provide a spatial map of the local lung tissue expansion and contraction. The term ventilation in this paper refers to this tissue expansion surrogate. The advantage of using 4DCT data to calculate ventilation is that a 4DCT scan is typically acquired as part of the treatment planning process. In addition, 4DCT‐derived ventilation maps have a high‐spatial resolution and show the temporal changes of the lung over different phases of the breathing cycle. Reinhardt et al.8 calculated ventilation from CT image volumes using the Jacobian determinant of the transformation calculated by deformable image registration. Simon9 and Guerrero et al.10 used density changes caused by the increase in air in the lungs to calculate ventilation. Kipritidis et al.11 omits image registration and uses the intensity values of the phase‐averaged 4DCT dataset to calculate ventilation. These ventilation calculations have been compared to clinically accepted measures of pulmonary function: Xenon CT, pulmonary function tests, and Galligas PET imaging, respectively.11, 12, 13 This paper estimates ventilation using the Jacobian determinant of the transformation.

1.B.Dose dependence of lung function

Many studies have calculated functional change in the lungs during or after radiation therapy using different modalities. Some studies14, 15, 16 assessed changes found in pulmonary function tests (PFT), which represent changes in the global function of the lung. Yuan et al. studied changes during RT and found no significant changes in PFT metrics. Vagane et al. and Guerra et al. showed no change in forced vital capacity or forced expiratory volume during the first second divided by vital capacity, but Guerra saw that the diffusing capacity of the lungs for carbon monoxide was reduced in a majority of patients. Non‐significant results in global PFT measures suggest looking at regional lung function metrics, which can be based on single‐photon emission computed tomography (SPECT) or 4DCT data. Yuan et al. showed that during treatment, ventilation/perfusion improved significantly.16 Farr et al. and Boersma et al.17, 18 looked at the change in SPECT data from pre‐RT to 5 and 3 months post‐RT, respectively. Both found that lung function decreased as dose to the lung region increased. Seppenwoolde et al.19 compared the radiation‐induced SPECT changes in non‐small cell lung cancer (NSCLC) patients to a reference group of non‐lung cancers. They found that in well‐perfused lung regions, the NSCLC patients responded similar to the reference group, but in poorly perfused regions, there was less decrease in function due to local reperfusion. Vinogradskiy et al.20 studied 4DCT‐based functional change during radiation therapy and found no apparent relationship between ventilation and dose. Ding et al.21 used 4DCT data to study the function change pre‐RT to post‐RT and found that there was a decrease in function for regions receiving greater than 24 Gy, and had varied responses in regions with less than 24 Gy. King et al.22 analyzed pre‐RT and post‐RT 4DCT datasets and found a significant dose‐dependent decline in function for voxels receiving over 20 Gy. To our knowledge, this is the first study to examine if the radiation dose response is dependent on 4DCT‐derived regional lung function prior to therapy. The purpose of this study was to quantitate the change in regional ventilation after a course of radiation therapy, and to establish its dependence on pre‐RT ventilation as derived from 4DCT datasets along with deformable image registration.

2.A.4DCT acquisition

This institutional review board‐approved prospective study (NCT01039649) evaluated twelve lung cancer subjects who underwent radiation therapy treatment. Exclusion criteria included prior or future planned surgery for treatment of the lung cancer, prior thoracic radiotherapy, uncontrolled intercurrent illness, pregnancy, severe COPD, underlying collagen vascular disease, oxygen dependence, recent abdominal surgery, and being under 18 yr old. The patient information for the study cohort is shown in Table 1. Two patients had early stage disease (I, II) and ten had later stage disease (III, IV). There was one patient treated with hypofractionation and the rest were treated with 1.8–2 Gy fractions. Prior to treatment, each subject had a 3D breath‐hold CT image acquired, which was used for radiation treatment planning. Three 4DCT scans were acquired for each subject: two scans before RT, with a 5–15 min break between scans, and one scan 3 months after completion of RT. The 4DCT datasets were reconstructed into 10 phases. All scans were acquired in helical mode with a 2 mm slice thickness, 0.5 mm slice increment, 0.5 s tube rotation time, 0.1 pitch, 120 kVp, 700 mAs, and 1.2 mm beam collimation on a Siemens Biograph 30‐slice CT scanner with subjects in the supine position. The scanner used an Anzai AZ‐773V with a strain gauge belt as a surrogate for the respiratory pattern. To increase the repeatability of the respiratory pattern, the subjects were guided with audible timing cues, which were played throughout the scan (RESP@RATE, Intercure Ltd., Lod Israel). The subjects in this study were selected from a protocol with a larger cohort; those excluded either withdrew prior to completion of the follow‐up scans or were missing datasets (31), or provided 4DCT scans with major artifacts (7).

2.B.Data preprocessing

Deformable image registration was used to find the pointwise correspondence between inspiration and expiration phase images. The registration algorithm uses a multi‐resolution cubic B‐spline parameterization and a sum of squared tissue volume difference (SSTVD) similarity metric. The Jacobian determinant of this transformation was used to estimate local ventilation. The registration accuracy was validated using anatomical landmarks,23 and the Jacobian ventilation calculation was validated by comparison to Xenon CT scans on ventilated sheep.8 A Jacobian value equal to one represents no local volume change, greater than one represents local expansion, and less than one represents local contraction of the lung tissue. This study used equivalent tidal volumes (ETV) for phase selection, which has been shown to increase the repeatability of the ventilation measurement, to isolate the effect of radiation in longitudinal changes.24 A 4DCT phase image at the beginning of inhalation was selected as the fixed image, and an image at the end of inhalation was selected as the moving image. The specific fixed and moving images were chosen to meet a 100 cc tidal volume range between all three scans with the goal of normalizing a subject's effort longitudinally.

2.C.Overview of Jacobian analysis

An overview of the Jacobian calculations and subsequent analysis of ventilation change for one subject can be seen in Fig. 1. Before RT, an expiratory breath‐hold 3DCT and two 4DCTs were acquired. The dose distribution was planned on the 3D breath‐hold CT since it was a higher quality image and free of motion artifacts. For subjects with both pre‐RT 4DCT images fitting the ETV criteria, the two pre‐RT Jacobian maps were deformably registered to the 3D breath‐hold CT and then geometrically averaged to get a more robust estimate of each subject's pre‐RT ventilation, J_preRT. For subjects with only one pre‐RT 4DCT image satisfying the volume constraints, this single Jacobian map was used as the estimate of their pre‐RT ventilation, J_preRT. Three months after completion of RT, another 4DCT was acquired and the Jacobian map was calculated to estimate the post‐RT ventilation, J_postRT. To visualize the consequences of the delivered dose, J_preRT and J_postRT were each deformed to the common coordinate space of the 3D breath‐hold CT. A Jacobian ratio map was calculated by taking a voxel‐wise ratio of J_postRT to J_preRT. This ratio map, J_postRT/preRT, quantitates the fractional change in ventilation due to RT. A ratio greater than one represents an increase in local ventilation post‐RT, and a ratio less than one represents a reduction in local ventilation. J_postRT/preRT was used in conjunction with the dose information from treatment to assess dose response at the voxel level. The 4DCT datasets, originally 2 × 0.98 × 0.98 mm, and dose distributions, originally 4 × 4 × 4 mm, were resampled to 1 × 1 × 1 mm for voxel‐wise comparison.

Figure 1

Diagram of the process used to measure the Jacobian change due to dose. A breath‐holdCTscan was acquired beforeRT, which was used to plan the dose distribution. Repeat 4DCTdatasets were acquired beforeRTand one dataset was acquired afterRT; these were used to calculate the pre‐RTand post‐RTJacobian, respectively. A Jacobian ratio map was obtained by taking the voxel‐wise ratio of post‐RTJacobian to pre‐RTJacobian. The Jacobian ratio map was analyzed with the dose distribution to calculate the dose dependence of ventilation.24 The 30 Gy isodose line is outlined in black.27

2.D.Determining response to radiation therapy and pre‐RT ventilation dependence

Dose response was assessed by first determining the calculated radiation dose delivered to each voxel of tissue within the lung. Doses were calculated for a heterogeneous medium derived from the CT dataset using a convolution superposition algorithm (Philips Pinnacle,3 Fitchburg, WI, USA), adjusted to perform adaptive sampling during convolution dose calculation. Specifically, the dose algorithm sampling was relaxed to every fourth point in the 4 mm × 4 mm dose grid in areas of flat dose. As described above, J_postRT/preRT provided the fractional change in ventilation following RT. This response map was deformed into the coordinate system of the computed dose and analyzed to establish the dose‐response relationship. Since there are several fractionation schemes, all delivered doses have been converted to equivalent dose in 2 Gy fractions (EQD2). The α/β of radiation pneumonitis is 4 and this value is used for the effects in this study due to the time to onset of pneumonitis corresponding with the 3‐month post‐RT time point.25

J_postRT/preRT and EQD2 were analyzed to quantitate the ventilation change due to dose. The magnitude of the dose response was assessed for all voxels within the lung, excluding those in the gross tumor volume. All voxels were stratified by a dose threshold of 20 Gy, which was defined due to its importance in treatment planning. The average J_postRT/preRT was calculated for the voxels in each of the two dose bins for each subject. The mean J_postRT/preRT per dose bin was then averaged over all subjects to quantitate the decrease in function between low‐dose and high‐dose tissues. The voxels in each dose bin were further stratified into two functional bins defined by a pre‐RT Jacobian threshold, creating four groups of voxels. This functional stratification was performed to determine if the dose response was also dependent on pre‐RT Jacobian. To establish a functional threshold, the complementary cumulative distribution of J_preRT was created for each subject. The median pre‐RT Jacobian for all subjects was 1.1, and was chosen as the threshold between high function and low function. The average Jacobian ratio was calculated for each of the four bins for each subject. The decline in function due to dose was calculated by the difference between the average of the high‐dose and low‐dose voxels for each functional stratification. To establish the significance of pre‐RT function, the Student's t test was used to compare the functional decline of the low‐function voxels to that of the high‐function voxels.

Keeping constant the two functional bins, the dose was further stratified from two dose bins into twelve dose bins of 5 Gy from 0 to 60 Gy to analyze the linear trend of J_postRT/preRT for both high‐ and low‐function lung tissue. Linear regression was performed on the two groups of data. The homogeneity of regression was tested using analysis of covariance to determine if there was a significant difference in the dose response between high‐ and low‐function voxels. Regression slopes were also tested for significance using analysis of variance. Cumulative dose‐volume histograms of the combined ipsilateral and contralateral lungs excluding the gross tumor volume were created to visualize the dose distributions and the relative volumetric extent of potential damage at different dose levels.

Finally, keeping constant the twelve dose bins, the voxels were further divided from two functional bins into six pre‐RT function bins of equal size from 1 to 1.6. This further stratification allowed the analysis of functional decline for more discrete groups of pre‐RT function. Linear regression was performed for each pre‐RT Jacobian bin looking at the response of J_postRT/preRT to increasing dose. The slopes of these lines quantitate the lung tissue's sensitivity to dose and were tabulated for each pre‐RT Jacobian.

3.A.Radiation‐induced changes

The J_postRT/preRT average for all subjects for voxels receiving less than 20 Gy was 1.014 ± 0.018 and 0.981 ± 0.024 for voxels receiving more than 20 Gy, resulting in a significant (P < 0.01) J_postRT/preRT decline of 0.033 for high‐dose voxels. To determine the amount of lung tissue that expands by more than 10%, the complementary cumulative distributions of J_preRT for all subjects were examined (see Fig. 2).

Figure 2

Complementary cumulative distribution, P(X > x), of J_preRTfor all subjects.

The median for all subjects, represented by the vertical red line, was 1.1, and is used to represent the threshold between high function and low function. The response between the high‐dose and low‐dose tissue for all voxels, low‐function voxels, and high‐function voxels is pictured in Fig. 3.

Figure 3

Bar plot of the average J_postRT_/preRTacross subjects for voxels receiving <20 Gy and >20 Gy for all voxels, voxels with low pre‐RTfunction, and voxels with high pre‐RTfunction. Error bars represent one standard deviation above and below the mean across subjects.

The average J_postRT/preRT for high and low pre‐RT function and high‐ and low‐dose groups for each subject are shown in Table 2. On average low‐function voxels receiving low dose had a J_postRT/preRT of 1.030 ± 0.027 while low‐function voxels receiving high dose had an average J_postRT/preRT of 1.007 ± 0.024, resulting in a relative decrease in J_postRT/preRT of 0.024 for voxels receiving more than 20 Gy. Similarly, high‐function voxels receiving low dose versus high dose had an average J_postRT/preRT of 1.001 ± 0.026 and 0.953 ± 0.036, respectively, resulting in a relative decrease of 0.048. The difference between J_postRT/preRT of high‐ and low‐dose voxels was significantly higher for high‐function voxels (P = 0.037), signifying an increased sensitivity to dose for tissue with high pre‐RT function.

3.B.Dose dependence of ventilation

Figure 4 shows the change in J_postRT/preRT as dose increases, categorized into regions of lung with high function and low function. Both the high function and low‐function trends are fairly linear, indicating an increasing reduction in lung tissue elasticity with increasing radiation dose. The low‐function voxels start at 0 Gy with a modest increase in function of 1.02 and end at 60 Gy with a ratio of one. The low‐function voxels have a very slight downward trend as dose increases, with a 0.3% ± 0.2% decrease in J_postRT/preRT per 10 Gy delivered dose. The high‐function lung voxels also started above a ratio of 1 at 0 Gy with an average Jacobian ratio of 1.01 decreasing to 0.93 at 60 Gy. The high‐function voxels showed a larger decrease in function over the range of doses, with a 1.4% ± 0.1% decrease in J_postRT/preRT per 10 Gy delivered dose. The regressions between the Jacobian ratios of high‐function and low‐function voxels are significantly different (P ≪ 0.001) calculated with analysis of covariance. Both slopes were also significantly different from zero (P ≪ 0.001) calculated using analysis of variance.

Figure 4

Average J_postRT_/preRTfrom 0 to 60 Gy for all subjects, stratified by pre‐RTfunction. The pre‐RTJacobian threshold between high function and low function was 1.1. Error bars represent values within one standard deviation of the mean.

Figure 5 shows the cumulative dose‐volume histograms of both ipsilateral and contralateral lung for the 12 subjects. Although there are a minority of voxels in the high‐dose region, where the separation between high‐ and low‐function voxels is most apparent, there are 12% of voxels over 40 Gy, an average of 443 cc per subject, which corresponds to a 5% reduction in function of high‐function voxels compared to low‐function voxels. Figure 6 shows the functional response with continued stratification of pre‐RT function into six discrete bins of width 0.1. The volume of lung in each of these bins used for averaging is shown in Table 3.

Figure 5

Cumulative dose‐volume histograms of combined ipsilateral and contralateral lung excluding theGTVfor all subjects.

Figure 6

Average J_postRT_/preRTfrom 0 to 60 Gy for all subjects, grouped by pre‐RTfunction. The pre‐RTJacobian was stratified into six equally sized bins between 1 and 1.6; the nominal bin centers represent the specific stratification of the line.

At high doses, the separation between the response of the pre‐RT Jacobian bins shows that higher functioning tissue incurs more radiation‐induced damage. Linear regression was performed for each pre‐RT Jacobian bin as a function of J_postRT/preRT and dose. The changes in J_postRT/preRT with respect to dose are derived from the linear regression slopes and are presented in Table 4. The decline in J_postRT/preRT increases monotonically with pre‐RT Jacobian from a decline of 0.3% ± 0.2% per 10 Gy for the lowest functioning tissue to a decline of 4.8% ± 0.4% per 10 Gy for the highest functioning tissue.