Use of deformable image registration techniques to estimate dose to organs at risk following prostate external beam radiation therapy and high-dose-rate brachytherapy

Purpose

Prostate cancer is the most commonly reported cancer in men globally and the second most commonly diagnosed cancer in men in Australia [1,2]. For intermediate to high-risk prostate cancer patients, external beam radiation therapy (EBRT) combined with a high-dose-rate brachytherapy (HDR-BT) boost dose has shown better biochemical relapse-free survival, local tumor control, decreased acute rectal and gastrointestinal toxicity, and distant metastasis-free survival compared with EBRT treatment alone [3,4,5]. During prostate cancer treatment planning and radiation treatment, it is important to have accurate estimates of total doses to organs at risk, such as the urethra and rectum, as excessive rectal doses have been shown to correlate with increased gastrointestinal toxicity [6,7,8]. Consensus guidelines provide dose constraints to minimize the risk of excess toxicity [6,8]. According to the Groupe Européen de Curiethérapie of European Society for Radiotherapy and Oncology (GEC/ESTRO), when determining the total dose delivered during combined EBRT and HDR-BT treatments, and to account for a variation in dose delivered per fraction, the total dose from each modality should first be converted to an equivalent dose in 2 Gy fractions (EQD2), and the combined modality dose assessed [6]. According to Hoskin et al., the recommended dose constraint to the rectum for combined modality treatment is D2cc ≤ 75 Gy EQD2 [6].
TRUS-based treatment planning for HDR-BT uses real-time treatment planning and eliminates the need to move the patient to a different room for imaging, since the same TRUS used for catheter placement guidance can also be applied for treatment planning [8]. Furthermore, treatment planning is performed in real-time between catheter placement and treatment, reducing time for edema to occur after catheter placement [9,10]. While the TRUS-based planning for HDR-BT has the above-mentioned advantages, there is a new challenge caused by the placement of transducer probe in the rectum. Inclusion of the probe in the rectum causes significant rectal deformation, forcing the anterior rectal wall to straighten and move closer to the prostate [10,11]. This large rectal deformation becomes challenging to estimate rectal dose for combined EBRT and HDR-BT treatment. In addition, there are further challenges in co-registering different imaging modalities (CT and TRUS), and in that these images do not include the same volumes (i.e., only anterior rectal wall included in TRUS images compared with full cross-section in EBRT).
One dose estimation approach called “parameter addition” includes full parameter addition (FPA) and partial parameter addition (PPA) [12,13]. FPA, also known as the “overlapping high-dose method”, adds together the maximum dose from the dose-volume histogram from each fraction in the dose accumulation approximations [12,13,14,15]. PPA is similar to FPA, except for the dose to target volume, which is considered homogeneously distributed [10,13]. Another dose estimation approach is to use image registration techniques, either rigid image registration (RIR) or deformable image registration (DIR) [14,16]. RIR comprises translation and rotation between moving and stationary images while keeping the same distances between corresponding points of the images [17]. DIR, on the other hand, allows for each voxel in the deforming image to be adapted to match the stationary image [17].
There have been many studies in the field of brachytherapy in recent years that have investigated how DIR dose estimation compares to parameter addition methods. For example, DIR has been applied to dose accumulation for combined EBRT and BT treatments of the cervix [13,18], and the prostate [12,19,20]. However, the afore-mentioned studies dealt with image registrations between the same modality or computed tomography-magnetic resonance imaging (CT-MRI), which have the problem of anatomical differences such as rectal filling [12,13,18,19], rather than the unique challenge faced by deformation of TRUS probe placed in the rectum. While previous investigations have shown a general trend that DIR does not produce significantly different results compared to parameter addition methods in EBRT and BT combined treatments, these studies have not investigated the effect of registering EBRT planning CT and real-time HDR-BT planning based on TRUS. The present study provides the first step in determining whether the use of DIR has the potential to improve dose estimation to organs at risk, given the substantial deformation of rectum caused by placement of the TRUS probe. The differences in dose estimates calculated via parameter addition methods and image registration methods were considered.

Material and methods

Patient treatment characteristics and dose constraints

This retrospective study included 10 patients treated with HDR-BT prescription dose of 15 Gy delivered in 1 fraction, followed 2 weeks later by EBRT 46 Gy in twenty-three 2 Gy fractions. The EBRT component of the treatment was planned with Pinnacle (v.9.10, Philips, Fitchburg, WI, USA) treatment planning system with a volumetric modulated arc therapy (VMAT) technique, and delivered through Truebeam (v.2.7, Varian Medical Systems, Palo Alto, CA, USA) linear accelerator. For the EBRT component of treatment, the PTV was the prostate contour, plus the seminal vesicles, expanded uniformly by 5 mm. There were no margins added to organs at risk (OARs) to create PRVs. The ultrasound-based HDR-BT component of treatment was planned with Oncentra Prostate brachytherapy treatment planning system (Elekta Brachytherapy, Veenendaal, The Netherlands). Contouring of the prostate, urethra, and rectum was performed by a radiation oncologist on a three-dimensional (3D) ultrasound image that was acquired via sagittal rotation of the TRUS probe in the rectum, before ProGuide plastic catheters (Elekta Brachytherapy, Veenendaal, The Netherlands) were inserted under live ultrasound guidance. A second 3D ultrasound image dataset was obtained to finalize the contours, and the reconstruction of catheters was performed by radiation oncologist and medical physicist on this image set. The source dwell positions and times were determined using dose-volume histogram-based inverse optimization algorithm within the Oncentra Prostate treatment planning system [21], and followed by quality assurance checks. The brachytherapy component of treatment used a 0 mm margin from CTV-PTV. The treatment was delivered immediately. For additional information on the HDR-BT treatment and planning, this method was also described by Poder et al. [20]. All contours were generated by the same radiation oncologist, who also approved all patients’ treatment plans.
Prior to initiation of this study, rectum constraints were considered separately. The HDR-BT component was delivered first and was optimized to ensure that the maximum dose to the rectum was less than 80% of the prescription dose (15 Gy). The EBRT component was delivered after HDR-BT and was optimized to ensure that: V44 Gy < 15%, V41 Gy < 20%, V38 Gy < 25%, V35 Gy < 35%, and V29.5 Gy < 50%. In this study, for the combined modality treatment, the dose constraint of interest was a minimum dose to the most highly irradiated 2 cm3 of rectum tissue (D2cc). The aim was to achieve D2cc ≤ 75 GyEQD2 according to the GEC/ESTRO recommendations [6]. This investigation also considered D1cc for the rectum as a surrogate for maximum dose delivered to the rectum [13]. MIM software (v.6.9.2, MIM Software Inc., Cleveland, Ohio, USA) was used to investigate each of the described registration methods. The DIR algorithm had previously been commissioned using recommendations of the AAPM TG-132 [17], and performance tests demonstrated compliance within the recommended tolerances.

Radiobiology

HDR-BT doses were first converted to an equivalent dose in 2 Gy fractions (EQD2) on a voxel-by-voxel basis, using α/β function in MIM software, which is based on the linear quadratic model [7,12,22] and given by: , where n is the number of fractions, d is the dose per fraction, and α/β is the cell fractionation sensitivity [23,24]. The prostate α/β was taken to be 1.5 Gy [7,24,25], and 3 Gy for the rectum [7,12,13,15].

Partial parameter addition method (PPA)

The partial parameter addition method, assumes the EBRT dose distribution to be homogeneous throughout the target volume (prostate), and the dose to organs at risk (rectum) to be equal to the prescription dose, i.e. 46 Gy in this case. The HDR-BT D1cc and D2cc for the rectum were derived from the DVH planning data. These dose values were added to the prescription dose of 46 Gy from EBRT to obtain the accumulated dose.

Full parameter addition method (FPA)

This method does not assume a homogeneous EBRT dose to target and OARs, so the DVH parameters must be obtained from the EBRT dose distribution as well as the HDR-BT dose distribution. Using this method, areas of maximum dose from each of the components of the treatment are assumed to be in the same anatomical location. The corresponding dose parameters were added to obtain the dose accumulated in each relevant volume, i.e., EBRT D1cc + HDR-BT D1cc = total rectum D1cc, EBRT D2cc + HDR-BT D2cc = total rectum D2cc.

Rigid image registration (RIR)

For the RIR used in this study, the translation and rotation functions in the MIM software were used to manually match the prostate volumes in all views, while monitoring the alignment between the anterior rectal walls of CT and TRUS datasets (Figure 1).

Deformable image registration (DIR)

The MIM software contour-based DIR algorithm was used to deform voxels within the TRUS image (chosen moving image) onto the CT image (chosen stationary image). This algorithm works by minimizing surface differences between the TRUS prostate contour and CT prostate contour, with surrounding voxels following the deformation of contours. The optimization method used in this DIR algorithm was modified gradient descent [14,26]. Only the prostate contours were selected for contour matching, as the rectum contours in each modality were mismatched due to the TRUS rectum contour only extending anterior to and within the window of transducer, while the full rectum volume was contoured in the CT images. The original 3D ultrasound image set had voxel dimensions of 0.5 × 0.5 × 0.5 mm3, and the original CT image set had voxel dimensions of 1 × 1 × 2 mm3 (2 mm distance in the superior-inferior direction). The deformed 3D ultrasound was re-sampled onto the reference frame of CT, resulting in a deformed 3D ultrasound with voxel dimensions of 1 × 1 × 2 mm3. To determine rectal D1cc and D2cc accumulated doses, the TRUS-based HDR-BT EQD2 dose distribution and CT-based EBRT dose distribution were transferred to a dTRUS dataset, and these two distributions were summed on the dTRUS dataset.

Image registration quality metrics

Qualitative evaluation

To qualitatively assess the accuracy of image registration, visual assessment of anatomical landmarks was performed. This included consideration of the alignment between the TRUS and CT prostate posterior borders, anterior rectal walls in the mid-prostate region, and an evaluation of general anatomy in the region of deformed contours.

Quantitative evaluation

Optimization of prostate contour match between the TRUS and CT datasets was considered to be the best method to achieve anterior rectal wall overlap in this study due to the large mismatch of rectum contours in each modality (TRUS rectal contours only extending anterior to and within the window of transducer, and CT rectal contour comprising the full rectal volume). The quantitative metrics were based on the recommendations of the AAPM Task Group 132 [17], and included Hausdorff distance (HD), mean distance to agreement (MDA), Dice similarity coefficient (DSC), and Jacobian determinant (JD). The HD is useful for providing an indication of outlier points on the volume since it is defined as the maximum distance of “closest points” between the two volumes/contours [27]. Both the DSC and MDA assess how well contours from the stationary and moving image match. The DSC indicates the degree of overlap between two contour volumes, and ranges from 0 (volumes are not well matched) to 1 (perfect match) [17,27]. The MDA metric involves converting the stationary image contour into points, and finding the distance from each point to the closest point on the transformed moving image contour and vice versa for the moving image contour, and finally averaging these distances [17,27]. A JD < 0 indicates that the deformation has caused a non-physical or non-anatomically possible event, such as folding of anatomy [17,27]. Negative values can also indicate that the algorithm cannot adequately process the deformation [17].

Uncertainty in dose accumulation

The uncertainty in dose accumulation for each registration method was obtained by calculating the dose gradient for the accumulated dose across the posterior prostate and anterior rectal wall at the point in the rectum where the dose was the highest (usually around the mid-prostate according to visual inspection). All dose measurements were suitably distant from HDR-BT catheters to prevent abnormalities due to the dose spikes at the catheters, i.e. less than 5% of the highest accumulated dose. The dose gradient was then multiplied by the MDA of prostate contours to find the uncertainty in dose, see equation below: In an attempt to consider the uncertainty in dose due to all factors rather than dose gradient alone, we also included an estimate of uncertainty due to dose delivery and contouring. Table 1 provides details of values and associated references that were used in the estimate of uncertainty. The combined and total estimates of uncertainty represent the total of all relevant uncertainties added in quadrature. The PPA and FPA uncertainties were based on the total combined uncertainty from the EBRT and HDR-BT components of treatment [28,29,30], while the RIR and DIR uncertainties were based on the total combined uncertainty from the EBRT and HDR-BT treatment as well as the uncertainty in dose accumulation for the image registration methods (see above equation). A two-sided Wilcoxon signed-rank test (p < 0.05) was performed in R (v.3.5.3, R Foundation for Statistical Computing, Vienna, Austria) to compare the results of quantitative metrics between RIR and DIR as well as to compare the results from each method used to estimate the rectal doses (D2cc and D1cc).

Results

Image registration

Qualitative metrics

Visual assessment of the rigidly aligned TRUS and CT prostate contours was performed. Six of the ten patients’ posterior prostate borders and anterior rectal walls were aligned locally at the mid-prostate level (Figure 2A). For the remaining patients, there was overall alignment with prostate and surrounding anatomy; however, local alignment between the TRUS and CT anterior rectal walls was poor due to curvature of the CT rectal volumes (Figure 2B). The alignment between the patients’ anterior rectal walls on CT and TRUS were observed to become slightly offset after performing DIR, relative to RIR for cases, where the anterior rectal walls were compatible in RIR. This was a consequence of total image deformation following the prostate contour deformation. Eight of the ten patients’ anterior rectal walls were aligned locally at the level of mid-prostate, while the other two patients did not align locally due to a greater offset between the CT and TRUS rectal walls after DIR (misalignment ranging between 0 and 3 mm maximum).

Quantitative metrics

The DSC results are shown in Figure 3. The average DSC was found to be 0.78 ±0.06 (1 SD) and 0.93 ±0.01 for RIR and DIR, respectively. The MDA results are presented in Figure 4. The average MDA was found to be 2.50 ±0.70 mm (1 SD), and 0.69 ±0.06 mm for RIR and DIR, respectively. The HD results are demonstrated in Figure 5. The average HD was found to be 11.64 ±2.38 mm (1 SD) and 5.19 ±1.47 mm for RIR and DIR, respectively. The outlier values were in regions where the CT contour over-estimated the prostate volume; for example, for most patients, there were differences in prostate contours in the region of seminal vesicles. Also, since the anterior prostate is difficult to delineate in CT, this region was often overstated in the prostate contours compared to the TRUS prostate contours [31]. These outliers did not affect the anterior rectal wall alignment. For all quantitative metrics, there was a statistically significant difference between the RIR and DIR results (p = 0.002 for all metrics). For all patients, the JD of the deformed image was greater than 0, indicating that all deformable registrations were physically and anatomically viable.

Dose accumulation

The results for the accumulated D1cc and D2cc from each dose estimation approach are presented in the following section as well as their total combined uncertainty. The total combined uncertainties for the PPA and FPA methods were based on the total combined uncertainty from the EBRT and HDR-BT components of treatment, while the RIR and DIR uncertainties were based on the total combined uncertainty from the EBRT and HDR-BT treatment as well as the uncertainty in dose accumulation for the image registration methods, as described in the Material and methods section. The accumulated D2cc for the rectum for each registration method is shown in Figure 6. On average, the predicted D2cc was higher than that calculated by the DIR method by 6.59 GyEQD2 (range, –3.03 to 13.68 GyEQD2) for PPA, 4.88 GyEQD2 (range, –3.41 to 11.97 GyEQD2) for FPA, and 3.46 GyEQD2 (range, –5.50 to 7.90 GyEQD2) for RIR. The accumulated D1cc for the rectum for each registration method is presented in Figure 7. On average, the predicted D1cc was higher than that calculated by the DIR method by 5.77 GyEQD2 (range, –1.91 to 13.20 GyEQD2) for PPA, 4.72 GyEQD2 (range, –2.00 to 12.16 GyEQD2) for FPA, and 2.88 GyEQD2 (range, –4.60 to 5.82 GyEQD2) for RIR. Wilcoxon signed-rank tests indicated that when calculating the D2cc, RIR was significantly different to PPA but not FPA, while DIR was significantly different to all other methods. When calculating the D1cc, RIR was significantly different to both parameter addition methods, and DIR was significantly different to all other methods.

Discussion

In this study, RIR provided a closer match between the anterior rectal walls of TRUS and CT datasets than DIR in cases where the anterior rectal walls of the CT dataset did not have much curvature in the prostate region. However, it was found that DIR provided good anterior rectal wall registration overall. This difference between the two registration methods was due to the ability to manually optimize the RIR method in MIM software, while the DIR method relied on matching the prostate contours, with the rectum contours following the general movements of prostate contours. Only the prostate contours were selected as the rectum contours in each modality were mismatched due to the TRUS rectal contours only extending anterior to and within the window of the transducer, while the full rectal volume was contoured in the CT images. This reliance on prostate contours resulted in the curvature of the TRUS anterior rectal wall changing after deformation and no longer matching the CT rectum as well as it did in some of the RIRs. However, DIR improved anterior rectal wall alignment in other cases where anterior rectal wall alignment was impossible with RIR alone, for example, where the CT rectal anterior wall was too curved for RIR.
The AAPM Task Group 132 recommends DSC between co-registered contours to be greater than 0.8-0.9 [17]. It was found that 5/10 of RIRs had a DSC > 0.8, while for DIR, all registrations had a DSC > 0.9. For the MDA, the AAPM Task Group 132 report recommends a tolerance that is within the largest voxel dimension between the co-registered images, which in this case, was 2 mm [17]. All DIRs met this tolerance, while only 3/10 RIR MDAs were within 2 mm. The results for RIR were as expected, since the difference between the prostate volumes defined on CT and TRUS ranged from 1.91 cc to 22.59 cc, with a CT/TRUS volume ratio of 1.3. Such a large volume difference has been previously reported [32] and was due firstly to difficulty of delineating the anterior prostate in CT, which usually results in over-estimation of prostate volume in this region, and secondly, due to differences in prostate contours in the region of seminal vesicles [31]. Note that despite a CT/TRUS volume ratio of 1.3, the typical distance between the CT and TRUS contours when aligned by their center of mass was only on the order of 1-2 mm. The results for DIR MDAs were also as expected due to fact that feature-based registration was employed, where the prostate contours were the chosen feature. The RIR uncertainty ranged from 10.9 Gy to 39.3 Gy, while the DIR uncertainty ranged from 5.98 Gy to 10.1 Gy, and this disparity was due to larger MDAs in RIR. The statistical results showed a significant difference between DSCs and MDAs of rigid and deformable image registration indicating that deformable image registration performs significantly better than rigid. However in this study, when considering prostate contours, the assumption that prostate contour matching provided good anterior rectal wall matching was not always true.
On average, the accumulated dose found using the DIR method was lower than using the PPA and FPA methods in most cases. This is in line with previous studies from literature [13,15,33], and is due to the fact that the PPA and FPA methods assume that the high-dose regions from both treatment modalities overlap, hence, typically overestimating the dose. In contrast, RIR and DIR methods attempted to register the anatomy and volumes first before accumulating the D2cc/D1cc, taking into account the fact that the high-dose volumes for each modality may not necessarily overlap. The current investigation found that for the case of EBRT combined with TRUS HDR-BT for prostate, the dose estimation to the rectum using DIR was significantly different to parameter addition methods. This is a new result and different to the findings of previous studies, which investigated the use of DIR for dose estimation in other treatments and/or modalities, and found no significant differences [12,13,15,18,19,20]. This new result was likely due to the large extent of deformation in the rectum created by the TRUS probe, compared to the rectal filling and emptying that was investigated in previous studies.
The DIR quantitative and qualitative evaluations met the AAPM TG-132 recommended tolerances, suggesting that that the DIR method of accumulating dose was an acceptable method for accumulating dose to rectum between TRUS and CT for prostate HDR-BT and EBRT. Furthermore, DIR overcomes the issues of dose over-estimation prevalent in PPA and FPA methods, by considering the location of high doses within the patient. The total combined uncertainties in the RIR method were found to be high, making RIR unreliable for dose constraint estimation.
A limitation of this study was that only feature-based DIR was feasible in MIM software due to multi-modality images used (CT/TRUS). Further to being a feature-based method, only the prostate contour was examined due the large discrepancy in rectum volumes contoured in each modality. Feature-based DIR in MIM software was originally tested using both the prostate and rectum contours; however, this registration was anatomically and physically unacceptable due to the afore-mentioned rectal volume discrepancies. In addition, this study focused on the MIM proprietary software and in-built functions. Open-source software, such as Slicer [34], or exporting data for post-processing, may allow more flexibility with respect to which volumes or contours can be registered, and may also improve performance between CT-TRUS registrations. Another consideration in terms of algorithm design, potentially in post-processing, is that either during or after registration of the prostate, re-sizing of the relevant CT rectal volume and registration of the rectum contours should also be attempted to improve dose estimations to this organ.
Another limitation of this study was that the dose to the urethra was not estimated. This is because the urethra was delineated and contoured in the TRUS image but not on the CT image, due to difficulty of delineation of the urethra in CT images [31]. Similarly, the bladder was not contoured on the TRUS image but was present on the CT image; therefore, the dose to the bladder for this treatment was also not estimated.

Conclusions

The DIR method used in this investigation was found to provide dose estimations to rectum D2cc and D1cc, with a statistically significant difference to both RIR and to parameter addition methods. This indicates that DIR may be an improved method for estimating the accumulated rectal dose for combined EBRT and TRUS-based HDR-BT prostate treatment compared to alternative methods. Whilst this investigation suggests that DIR may provide improvements over parameter addition methods, further investigation is required to validate the accuracy of accumulated dose using DIR, with a necessary quality assurance before clinical implementation.

Disclosure

The authors report no conflict of interest.

References