Dose reconstruction technique using non-rigid registration to evaluate spatial correspondence between high-dose region and late radiation toxicity: a case of tracheobronchial stenosis after external beam radiotherapy combined with endotracheal brachytherapy for tracheal cancer.

PURPOSE: Small organ subvolume irradiated by a high-dose has been emphasized to be associated with late complication after radiotherapy. Here, we demonstrate a potential use of surface-based, non-rigid registration to investigate how high-dose volume topographically correlates with the location of late radiation morbidity in a case of tracheobronchial radiation stenosis.
MATERIAL AND METHODS: An algorithm of point set registration was implemented to handle non-rigid registration between contour points on the organ surfaces. The framework estimated the global correspondence between the dose distribution and the varying anatomical structure. We applied it to an 80-year-old man who developed tracheobronchial stenosis 2 years after high-dose-rate endobronchial brachytherapy (HDR-EBT) (24 Gy in 6 Gy fractions) and external beam radiotherapy (EBRT) (40 Gy in 2 Gy fractions) for early-stage tracheal cancer. RESULTS AND
CONCLUSIONS: Based on the transformation function computed by the non-rigid registration, irradiated dose distribution was reconstructed on the surface of post-treatment tracheobronchial stenosis. For expressing the equivalent dose in a fractional dose of 2 Gy in HDR-EBT, α/β of linear quadratic model was assumed as 3 Gy for the tracheal bronchus. The tracheobronchial surface irradiated by more than 100 Gyαβ3 tended to develop severe stenosis, which attributed to a more than 50% decrease in the luminal area. The proposed dose reconstruction technique can be a powerful tool to predict late radiation toxicity with spatial consideration.

Purpose

For early-stage tracheal cancer, high-dose-rate endobronchial brachytherapy (HDR-EBT) can be a treatment choice as a boost for external beam radiotherapy (EBRT) or as a definitive therapy [1, 2, 3]. High-dose-rate endobronchial brachytherapy is also used successfully as a palliative treatment for symptomatic airway obstruction by malignant tumors [4, 5]. However, in some cases with longer follow-up, HDR-EBT leads to severe complications, which are sometimes lethal, such as radiation-induced bronchitis, bronchial stenosis, and fetal hemoptysis [6, 7]. Although these complications might reflect the high-dose region resulted from the small distance between the applicator and tracheobronchial mucosa [8], further investigation is needed on how high-dose volume topographically correlates with the location of late radiation injury. However, there exists technical challenge because the location of small high-dose volumes in fractionated radiotherapy can be affected by intra- and inter-fractional organ motion, and the anatomical features of an organ post-treatment will be different from those pre-treatment due to the sequelae of radiation treatment [9].

For the purpose of spatiotemporal dose assessment of HDR-EBT in combination with EBRT, one important technique is non-rigid image registration, which enables the estimation of dose summation based on the corresponding structures across the image data sets. Previously, we have demonstrated an image-processing framework using surface-based, non-rigid registration for assessing the spatially cumulative dose in consecutive brachytherapy treatments [10]. The core of our approach is based on the point set registration using a Gaussian mixture model (GMM-REG), which was originally introduced by Jian and Vemuri [11]. Because GMM-REG demonstrated the biomechanically natural registration with high accuracy and validity [12], we considered that it can be a plausible technique for transforming planned dose distribution to the post-treatment organ.

Here, we present a case of tracheobronchial radiation stenosis after the definitive therapy using HDR-EBT and EBRT for early-stage tracheal cancer. By using the proposed dose reconstruction technique, the spatial association between the dosimetric profile and the site of the radiation morbidity with a severe tracheobronchial stenosis was evaluated.

Material and methods

Patient data

A 78-year-old man developed a second primary tracheal cancer 5 years after a right lower lobectomy for early stage lung carcinoma. Bronchoscopy revealed a tumor localized at the left lateral wall of the middle part of the trachea and histology confirmed a squamous cell carcinoma. The patient was then treated with HDR-EBT (24 Gy in 6 Gy fractions) followed by EBRT (40 Gy in 2 Gy fractions) with a curative intent. Follow-up bronchoscopies revealed a patent airway with granulomatous mucositis for 1 year after the treatment, with no evidence of biopsy-proven malignancy. Twenty-four months after the treatment, at the age of 80, the patient was admitted to our department because of progressive dyspnea. Computed tomography (CT) detected a tracheobronchial stenosis with aspiration pneumonia of the right lung. Massive pneumonia infiltrated the right lung completely due to the stenosis of the right primary bronchus. Bronchoscopy revealed that the surface of the tracheal mucosa was irregular with inflammation, and bled easily on contact. It also revealed a small ulcer above the circumferential tracheal stenosis. Because the tissue examination of the tracheal mucosa was consistently negative for malignancy before the admission, we diagnosed late radiation-induced mucositis with severe circumferential stenosis.

Treatment protocol

Our HDR-EBT treatment protocol has been described in detail elsewhere [6]. In short, a 12 Fr Malecot-type catheter with two wings was placed from the middle part of the trachea to the right primary bronchi during bronchoscopy under local anesthesia. Based on the CT images with 2 mm slice thickness, treatment planning was done on Oncentra (Nucletron, an Elekta company, Elekta AB, Stockholm, Sweden). A dose of 6 Gy was prescribed at 1 cm from the catheter. High-dose-rate endobronchial brachytherapy was performed using the Microselectron 192Ir-HDR (Nucletron). Subsequently, EBRT was delivered using a three-dimensional conformal technique with linear accelerator (Clinac iX, Varian Medical System, Palo Alto, CA, USA). Treatment planning was based on the CT images of 3 mm slice thickness taken with a CT simulator (Aquilion™ LB, Toshiba Medical Systems, Tokyo, Japan). The EBRT fields included the primary tumor with adequate margins of 2 cm.

Image-processing framework

The analysis was performed using software developed in-house, which was written in Python using VTK/ITK library and published modules of GMM-REG [11]. A single observer (KK) contoured the inner surface of the tracheal bronchus on every CT image at the treatment planning (the 1st to 4th fraction of HDR-EBT and EBRT) and at post-treatment evaluation for the radiation stenosis. Delineation, which was bounded by the glottis superiorly and the level of second carina inferiorly, was performed using treatment planning software (Oncentra). Then, the contours were exported as binary masks from digital imaging and communication in medicine, and these binary masks were smoothened via Gaussian filter with a sigma value of 2.0 for reducing artifacts. A marching cubes filter was applied to produce surface mesh, which was represented by points that were connected with each other to form triangles. For refining the triangles which defined the surface, a quadric clustering filter with the division spacing of 2.0 mm was used as the decimate modifier to reduce the number of vertices of the mesh.

Next, non-rigid registration based on GMM-REG was guided by the point set distributed on the surface mesh. The theoretical background has been described in detail in the original article written by Jian and Vemuri [7]. Briefly, a Gaussian mixture is generally defined as

where ω is the weight factor associated with the Gaussian component density φ (x|μi, Σi), μi is the mean vector and Σ is the covariance matrix. The scale factor σ represents the size of the spherical covariance matrix of the model, which effects the magnitude of transformation.

GMM-REG treated the problem of point set registration as that of aligning two Gaussian mixtures by minimizing the L distance between them. Thin-plate spline (TPS) function was selected for the non-rigid registration, which added a penalty term to the final cost function to control the strength of regularization. Therefore, the registration method between the model set M and the scene set S finds the warping function u of the TPS function, which minimizes the following cost function in the presence of the weight parameter λ for the bending energy:

The process was driven by a numerical optimization algorithm called the annealing step to minimize local minima. We defined the total number of steps for annealing as five. The transformation parameters of each step was as follows:

Finally, the transformation function computed vectors that connected the points on one tracheobronchial surface to another tracheobronchial surface, which enabled the estimation of cumulative dose on the post-treatment tracheobronchial stenosis. For expressing the equivalent dose in a fractional dose of 2 Gy (EQD2) in HDR-EBT, α/β of a linear quadratic model was assumed as 3 Gy for tracheobronchial mucosa [13].

Validation for transformation

To measure the geometric accuracy of the transformation, surface distance error (SDE) and surface coverage error (SCE) were used [14]. Surface distance error was defined as the mean distance between the transformed points and the fixed surface. Surface coverage error represented the coverage of the fixed surface by the transformed points and was derived by calculating the standard deviation (SD) of the distances from each vertex of the fixed structure to the nearest transformed point. Transformation with small SDE and small SCE indicated a high geometric accuracy.

Homogeneity index

For accessing the homogeneity of the dose distribution, we calculated the homogeneity index (HI) as follows:

HI = D /D

where D is the maximum point dose to the inner surface of the tracheal bronchus and D is the prescribed dose.

Quantification of tracheobronchial stenosis

To quantify the severity of the stenosis, the luminal area of the tracheobronchial stenosis was evaluated as an intersection between the tracheobronchial surface and the plane, which was orthogonal to the principal axis of the organ. The principal axis was obtained using principal component analysis for the points on the structure surface. Each corresponding plane was computed based on the deformation fields for each tracheal bronchus at the time of treatment planning (the 1st to 4th fraction of HDR-EBT and EBRT), and the mean area of intersection was defined as the pre-treatment luminal area (Figure 1A). The local dose at each plane was calculated to be the average dose irradiated to the points constituting the outer line of intersection of the post-treatment tracheobronchial structure. The relationship between the local dose and the ratio of luminal areas was calculated from the level of the plane 100 mm above the first carina (inferiorly) to the plane just below the glottis (superiorly). Spearman's rank correlation coefficient was calculated for the association between the irradiated dose and the ratio of luminal areas. Statistical significance was set at a two-sided p-value of less than 0.05. Data are presented as the mean ± standard deviation unless otherwise specified.

Fig. 1

The dose-response relationship between irradiated local dose and severity of the stenosis was quantitatively evaluated. A) The luminal area of the tracheobronchial stenosis was computed as the intersectional area between the surface and the plane orthogonal to the principal axis of the organ subvolume. The corresponding intersectional area of each pre-treatment structure (the 1st to 4th fraction of high-dose-rate endobronchial brachytherapy and external beam radiotherapy) was obtained on the basis of deformation fields. B) The relationship between the local dose and the ratio of luminal areas was calculated

Results

Transformation results

An example of the deforming model point set to the scene point set on the basis of GMM-REG was shown in Figure 2. When the structure of post-treatment tracheobronchial stenosis was transformed to each tracheobronchial surface at the treatment planning (the 1st to 4th fraction of HDR-EBT and EBRT), SDE and SCE was 0.18 ± 0.01 mm and 0.39 ± 0.01 mm, respectively.

Fig. 2

An example result of the registration using Gaussian mixture model (GMM-REG). A) Control points on the surface of post-treatment tracheobronchial stenosis were used as model point set, which were transformed to the pre-treatment tracheobronchial surface as scene point set. B) Deformation vectors calculated by the transformation function. C) According to the number of iterations, the model point set was gradually registered to the scene point set with biomechanically natural transformation

Mapping reconstructed dose distribution

We reviewed the planned dose distribution on the organ surface at the time of treatment planning. As shown in Figure 3, the surface dose of EBRT was homogenously distributed (HI = 1.1); on the other hand, that of HDR-EBT varied unevenly from fraction to fraction (HI = 3.7 ± 0.2). Unexpected hot spots irradiated by more than 20 Gy per fraction resulted from the small distance between the catheter and the tracheobronchial inner surface (Figure 4), as the minimum distance between them was 2.4 ± 1.1 mm for the four fractions of HDR-EBT. The dose-surface histogram of the HDR-EBT was calculated (Figure 5) and the area of the tracheobronchial mucosa irradiated by more than 20 Gy was 94.7 ± 82.7 mm2 (0.6 ± 0.5%).

Fig. 3

Surface doses on the pre-treatment tracheobronchial surfaces. For each fraction of high-dose-rate endobronchial brachytherapy (24 Gy in 6 Gy fraction) and the total external beam radiotherapy dose (40 Gy in 2 Gy fraction), irradiated doses were mapped on the inner surface of the tracheal bronchus

Fig. 4

An axial computed tomography image at the time of the second high-dose-rate endobronchial brachytherapy treatment planning is shown as an example of a close distance between the catheter (red point) and the tracheal mucosa (white contour) corresponding to the high-dose region (arrow)

Fig. 5

Dose-surface histogram for each fraction of high-dose-rate endobronchial brachytherapy

According to the non-rigid transformation function, we mapped cumulative dose during the course of the radiation therapy on the surface of the pre-treatment (the structure at the time of the treatment planning of EBRT) (Figure 6) and post-treatment tracheal bronchus (Figure 7). For each fraction of HDR-EBT (24 Gy in 6 Gy fractions) and the total EBRT dose (40 Gy in 2 Gy fractions), the fractional dose was normalized to EQD2 with α/β of 3. On the pre-treatment tracheal bronchus, the area with cumulative dose more than 100 Gyαβ3 and 200 Gyαβ3 was 23.3 × 102 mm2 (18.6%) and 5.7 × 102 cm2 (4.5%), respectively. Figure 7 shows that the subvolume of the tracheobronchial surface irradiated by more than 100 Gyαβ3 developed stenosis, and the severity of the stenosis apparently corresponded with the cumulative local dose. When quantitating the structural deformation of radiation stenosis by the ratio of luminal areas, the subvolumes irradiated by more than 100 Gyαβ3 showed severe stenosis with a more than 50% decrease in the luminal area (Figure 1B). The association between the local dose and the ratio of luminal area was statistically significant (p < 0.0001). In particular, small volumes irradiated by more than 200 Gyαβ3 correlated with bronchoscopically-proven ulceration and the circumferential stenosis leading to a severe decrease in the luminal diameter (Figure 8).

Fig. 6

Cumulatively irradiated dose on the pre-treatment tracheal bronchus at the time of the treatment planning of external beam radiotherapy. Irradiated doses over the course of the radiotherapy were summed and represented by EQD2 using α/β = 3

Fig. 7

Cumulatively irradiated dose on the structured trachea as represented by EQD2 using α/β= 3. The degree of tracheal stenosis was higher in the sites which were irradiated by higher doses

Fig. 8

Small regions irradiated by more than 200 Gyαβ3 were associated with severe clinical complications. Arrows represent the structured trachea at each plane on computed tomography image. Arrowhead at plane A indicates the bronchoscopicallyproven ulceration, and arrowhead at plane B shows the circumferential stenosis

Discussion

In our study, we have presented a novel dose reconstruction technique to investigate the topographic correspondence between high-dose volume and the site of late radiation injury represented by tracheobronchial stenosis. Although the post-treatment tracheobronchial surface showed a large anatomical deformation, the measured low-distance error implies that the proposed approach with GMM-REG successfully registered the structures. As a result, the estimated cumulative local dose on the tracheobronchial surface significantly corresponded with the severity of radiation stenosis. The present study is of importance because the surface-based, non-rigid registration can be a core of the technique to predict the location-by-location difference in the severity of late radiation injury.

While there is variation in dose-fractionation regimens and treatment intentions among previous studies, the crude incidence of symptomatic radiation-induced bronchitis is 4-12% in HDR-EBT with EBRT [7, 15, 16, 17, 18]. Speiser and Spratling proposed a grading system to classify radiation-induced bronchitis, which ranges from mild mucosal inflammatory response with swelling (Grade 1) to major fibrosis with circumferential stenosis (Grade 4) [7]. According to the grading scale, our presented case represented the most severe phase of Grade 4 radiation tracheobronchitis in small regions, which were cumulatively irradiated by more than 200 Gyαβ3. Severe inflammatory response with mild fibrosis of the tracheobronchial wall, which corresponds to Grade 3, was apparently associated with mucosa irradiated by more than 100 Gyαβ3.

It is necessary for radiation oncologists to avoid excessive irradiation to the tracheobronchial mucosa because once radiation-induced tracheobronchitis becomes progressive, it never heals spontaneously, and often develops into fatal hemoptysis. For the present case, we prescribed a dose of 6 Gy at 1 cm from the catheter without enough spatial consideration of the distance from the tracheobronchial surface; thus, this may have been a primary cause of incidental high-dose volume.

There have been few reports on the spatial correspondence between high-dose volume and the site of radiation injury. By using rectosigmoidoscopy and rectal morbidity assessment, Georg et al. reported the locational correlation of dosimetric parameters of gynecological brachytherapy with clinical changes [19]. Although they emphasized the clinical relevance of DVH parameters, such as D2cc, D1cc, and D0.1cc, it might be an alternative approach to focus on the surface dose for reporting when non-rigid registration is used, because non-rigid registration sometimes produces artifacts by condensing or stretching volumes [20]. The present study also implies that the surface dose can be used in the prediction of radiation morbidity in the walls of a hollow organ.

Although the present study proposed a non-rigid registration-based assessment for the dose reconstruction of a previously irradiated dose to the post-treatment organ with large deformation, the technique has so far been applied in only one case of tracheal cancer treated with HDR-EBT and EBRT. General validation for other treatment sites with other organs at risk, different inter-fraction variations, and more inhomogeneous and complex dose distributions is still a big challenge.

Conclusions

The image-processing framework based on GMM-REG enabled to correlate the cumulative surface dose of HDR-EBT and EBRT with the location of radiation injury, by applying it to one case of tracheal cancer. This technique seems to be a promising tool that might be applied to other treatment sites in the future.