Automated calculation of point A coordinates for CT-based high-dose-rate brachytherapy of cervical cancer

Purpose

Locally advanced cervical cancer is routinely treated with concurrent chemotherapy and external beam radiation treatment (EBRT) followed by high-dose-rate (HDR) brachytherapy using a tandem and either ovoid pair or ring applicators [1]. For disease involving lymph nodes or parametrial extension, an additional EBRT boost is prescribed [2]. The American Brachytherapy Society (ABS) recommends reporting point A dose for treatment planning for HDR cervical brachytherapy [3]. When 3D imaging is used for image-guided brachytherapy (IGBT), magnetic resonance imaging (MRI) can be used to delineate the high-risk clinical tumor volume for volumetric planning [4,5,6]. When MRI scans of implanted HDR applicators are not readily available for target delineation, computed tomography (CT) is typically used to define the implant geometry, dose prescription points (point A), and organs-at-risk [3,7,8,9]. For dose calculation, the ABS recommends use of point A in their reports [1,3], located 2 cm perpendicularly to the tandem plane, and 2 cm plus the radius of the ovoids superiorly along the tandem from a line connecting mid-dwell positions of the ovoid pair. Even for volume-based planning, ABS [3], the Groupe Européen de Curiethérapie and the European Society for Radiotherapy and Oncology (GEC-ESTRO) [7] require reporting of dose to point A: “The dose delivered to point A should be reported for all brachytherapy applications regardless of treatment-planning technique” [3]. Defining a consistent dose point relative to the applicator would facilitate correlation of dose with clinical outcomes as clinics transition from 2D planning to volumetric 3D planning, according to GEC-ESTRO’s recommendation “to strictly follow the rules of a certain system with an enduring clinical tradition (for dose prescription)” [7]. With more advanced imaging technologies and brachytherapy modalities, volumetric prescription to the clinical target volumes is becoming prevalent for IGBT [9]. However, as of 2015, nearly 50% of centers continue to prescribe to point A despite the predominate use of CT-based planning [8]. Automated point coordinate calculations reduces dependence on either the user or the treatment planning system (TPS), thereby alleviating potential human errors in dose calculation point placement. Such errors were identified as one of the high-risk failure modes in HDR planning procedures [10,11,12,13]. Automated calculation of point A could also increase efficiency in HDR brachytherapy treatments of cervical cancer [11,14]. We have developed an automated point coordinate calculation algorithm for Fletcher-Suit-Delclos (FSD) tandem and ovoid (T&O) implants (Varian Medical Systems, Palo Alto, CA, USA) in order to expedite the process while simultaneously improving the precision of point localization.

Material and methods

High-dose-rate brachytherapy treatment planning data

This study utilized 34 CT scans with 2 mm slice thickness from 7 cervical cancer patients implanted with titanium FSD T&O for HDR treatment at our institution in 2014. High-dose-rate treatment plans were performed in BrachyVision v.11 (Varian Medical Systems, Palo Alto, CA, USA). To reduce applicator positioning uncertainties due to imaging artifacts, we utilize the solid applicators with pre-defined geometries provided by the manufacturer. The solid applicators are inserted into the CT scans, and the user translates and rotates them to match the implanted FSD T&O. All plans are generated with the 0.5 cm length 192Ir source (VariSource, Varian Medical Systems, Palo Alto, CA, USA) using a contiguous 0.5 cm step size and a 0.3 cm retraction from the distal tip to account for wire over-travel. The total source dwell length in the tandem is determined depending upon the cervical sleeve length inserted into the patient and the location of the cervical os. A standard source dwell length of 2 cm (i.e., 4 active dwell positions) is utilized in each ovoid. After insertion of the solid applicators, the user manually places points A by translating and rotating to multiple 3D views intersecting the applicators.

Automated calculation of point A coordinates

A stand-alone Matlab (Mathworks, USA) executable application calculates point A locations described in the previous section. The GUI-based application reads in a text file specifying the geometry of the implanted applicators, including the dwell positions in the applicators and the diameter of the ovoids, which is generated from the plan. The application first calculates the “tandem plane” in the sagittal plane, in which all tandem source dwell positions reside. The “tandem plane” is defined from the tip, middle, and most inferior source dwell positions in the tandem. The point A coordinate calculation begins by defining the ovoid central plane (Figure 1 – A), in which both ovoid dwell lines, assumed to be coplanar, are contained. The ovoid central plane is determined from any combination of three points in the ovoid dwell positions, one from one ovoid and two from the other ovoid. For our calculation, the 2nd dwell position in the left ovoid (LO2), and the 1st (RO1) and 3rd (RO3) dwell positions in the right ovoid are selected. The vector normal to the ovoid central plane (nˆ, Figure 1 – C) is calculated as in Eq. 1. (Eq. 1) where MIDOV1 and MIDOV2 are the vectors connecting two dwell points between the ovoids. The “plane at ovoids’ top surface” (Figure 1 – B) is determined from the three points, P0, P1, and P2 that are calculated from Eq. 2, where Rovoid is the radius of the ovoids. P0 = Rovoidnˆ + LO2, P1 = Rovoidnˆ + RO1, P2 = Rovoidnˆ + RO3 (Eq. 2) The “point at the ovoids’ top surface”, Ovtop (marked as × in Figure 1) is the point on the tandem intersecting with the plane at ovoids’ top surface. The point A is located 2 cm superiorly to Ovtop and 2 cm laterally in the normal direction to the tandem plane. To define Ovtop, the algorithm first finds the tandem segment vector, which crosses the plane at the ovoids’ top surface. The tandem segment vectors are defined between consecutive tandem dwell positions, vi = ti+1 – ti. Here ti and ti+1 are the ith and (i + 1)th dwell positions from the most inferior source position in the tandem. Each tandem segment line is determined by the general parametric equation for a line, ti + vi × s, where vi × s is the difference vector between ti and ti + 1, and represents the distance from ti to a point on the line along vi. Ovtop is the point where a tandem segment line intersects with the plane at the ovoids’ top surface. This surface is defined using the general equation for a plane, P0 + (P1 – P0)u + (P2 – P0)w, where P0, P1, and P2 are the three points as described above, and u, w  R, are parameters of a point in the plane formed by P1 – P0 and P2 – P0. To search for this intersecting point Eq. 3 is solved to determine s for each tandem segment line. ti + (ti+1 – ti)s = P0 + (P1 – P0)u + (P2 – P0)w (Eq. 3) The tandem segment line, which crosses the plane is identified by s  [0, 1]. The resulting crossing point, calculated from ti + vi × s, is assigned as Ovtop. Once Ovtop is determined, the algorithm locates the point on the tandem 2 cm superior from Ovtop. The “vector along the tandem”, Talong, defined between Ovtop, and a tandem dwell position typically 2-2.5 cm above Ovtop, traces the superior-inferior direction of the tandem as shown in Figure 1 – D as a dotted line vector. The algorithm finds the point 2 cm away from Ovtop along Talong, and projects it back onto the actual tandem. To find this point, first the tandem segment vectors vi are each projected onto Talong. The projection of vi on Talong is vicosi, where cosi is the cosine of the angle between vi and Talong. The projection is computed and accumulated one segment vector at a time until the final (fth) segment is reached. The final segment and its projection are computed as follows. The projections of individual segment vectors are accumulated as until the next vi+1 projection to be added would cause the sum dsum, to exceed 2 cm. For the final (f+1)th segment projection, the remaining distance dremain= 2 cm – dsum on Talong is projected back onto the actual tandem to locate the corresponding point by solving . Finally, the points A are located 2 cm away from this point on the tandem described above in the direction normal to the tandem plane.

Data comparison and statistical analysis

To evaluate the accuracy and precision of the algorithm, a single experienced user (HK) retrospectively and manually placed points A for each plan to compare with the automatically calculated point coordinates. Because retrospective manual point placement by a single user avoids introducing inter-user variability and alleviates time constraints/pressure, points were considered the “gold standard”. To evaluate intra-user variability, the same user (HK) repeated point A placement in a second but independent session. To characterize inter-user variability within a time-pressured clinical setting, points A defined by various users for these clinical plans are also compared. Doses were calculated to check the dosimetric effects of point coordinate differences. The time span required to place points A manually were tabulated when the “gold standard” was placed, and used to estimate the increased efficiency that could potentially be gained by automatic point coordinate calculation. Differences from the gold standard coordinates were calculated for automatically placed points, the repeated placements by a single user, the clinically used coordinates placed by multiple users, and used to assess the algorithm accuracy, intra-user variability, and inter-user variability, respectively. Coordinate differences in each dimension, in 3D magnitude, and associated dose ratios were compared for the three groups using a non-parametric, paired Friedman test, which does not require data be normally distributed. If differences in a group were found, a paired Wilcoxon test was used to identify the pair of groups that differed. To compare the variances of these coordinate differences among the same three groups, the nonparametric counterpart to the F-test, the Ansari-Bradley test, was used. All statistical tests were assessed at a significance level of 0.05.

Results

Table 1 summarizes the signed coordinate differences (Δx, Δy, Δz) and their 3D magnitudes (Δ3D) as well as dose ratios averaged over left and right points A from 34 HDR CT scans. All differences are referenced to the gold standard and thus summarize algorithm accuracy (first row), intra-user variability (second row), and inter- user variability (third row). For 60 of 68 points, automatically calculated coordinates agree within 0.5 mm in each direction with the gold standard. Differences larger than 0.5 mm occur only in the S-I direction. The 3D magnitude differences between the automatically-calculated and gold standard points agree to 0.35 mm on average, with a maximum difference less than 0.7 mm. In the majority of cases (82%), the dose to point A determined using the automated calculation agrees with the gold standard to within 1%, with a maximum dose difference in the cohort of 2.3%. The intra-user variability in 3D magnitude is 0.38 mm on average with a maximum dose difference of 3.6%. However, the inter-user variability in 3D magnitude increases to 0.71 mm on average, with a maximum dose difference of 7%. Statistical comparisons between each row of Table 1 indicated no differences among mean values of Δx, Δy, or Δz. However, the variances in each dimension were statistically larger (p < 0.05) for inter-user variability compared to either intra-user variability, or overall algorithm accuracy. When comparing Δ3D, both the mean and variance for inter-user variability were significantly different from either intrauser variability or overall algorithm accuracy. No significant differences in means or variances were found when comparing algorithm accuracy to intra-user variability. Dose ratios were 1-2% on average and were not significantly different among the three groups, although the maximum dose differences did increase from 2.3%, to 3.6%, to 7% for the studies of the algorithm’s accuracy, intra-user variability, and inter-user variability, respectively. To further analyze the relationship between point A localization accuracy and calculated dose, a simulation was used to shift point A in L-R and S-I directions. Larger dose differences, up to 5%, were found for L-R shifts than for S-I shifts, which only comprised differences to within 1%. These dose differences were within the range of those shown in Table 1.

Discussion

We developed a stand-alone application to automate point A coordinate calculation based solely on T&O implant geometry for HDR brachytherapy. This tool expedites planning and reduces user-dependence from the dose point localization process. The point placement time was reduced to less than 3-5 minutes, including generating the text input file and typing the point coordinates back into the TPS. The calculation time of the algorithm is less than one second. Because the initial planner and the second physicist checking the plan must both localize points A independently, the mean time added to the overall planning process is on the order of 10 minutes. An advantage of automating point A coordinate calculations is to minimize the potential for human error and reduce the variability among various planners. Automation can improve efficiency and result in consistent point A localization even for planners who are not familiar with a TPS. In this study, the automatically calculated point A coordinates show submillimeter localization agreement to 0.35 mm on average, and dosimetric agreement to 1% on average compared to retrospective point insertion by a single experienced user (i.e., the gold standard). While a single user localizing point A in a retrospective and non-rushed setting could mimic the high accuracy and precision of the algorithm, this was not true when multiple planners inserted point A into clinical plans in a time-pressured setting. Our study demonstrated that for multiple users, the 3D magnitude differences were significantly larger as was the variance of the coordinates in x, y, z, and 3D. The localization accuracy achieved by a single user in a retrospective setting could have been achieved at the expense of time, potentially leading to an overestimation of the time efficiency afforded by the point placement algorithm. Although the coordinate localization accuracy was decreased among multiple planners, the doses were not significantly impacted to within 1-2% on average. While not significant, differences in dose were up to 7% as shown in Table 1. In addition, we simulated the effects of 1 mm-shifts in L-R and S-I on dose and found that dose differences were larger, up to 5%, when point A was shifted laterally to the tandem compared to when it was shifted in S-I, for which differences were within 1%. The accuracy of both the coordinates and doses at the automatically placed points A demonstrates that the algorithm is able to perform robustly, thus reducing potential human errors arising from clinical time pressure. Studies on failure mode in high-dose-rate dosimetry identifies incorrect dose point placement as a mode with high severity [10]. Other studies have shown that failure modes associated with the highest risks originate most frequently from the treatment planning process [12] and from human failures [13]. Automated point A calculation assumes that all ovoid central dwell positions are contained on the ovoid central plane (Figure 1 – A), and thus, any combination of three dwell positions should describe the same plane. This assumption could be violated in a clinical implant, particularly when the cervical anatomy of patient is atypical leading to asymmetric ovoid geometry. Non-coplanar ovoids can result in different central ovoid planes depending on which three dwell positions are selected to define the plane. To estimate the effects of non-coplanar ovoids and errors in localizing the ovoid tip on our methodology, point A was re-calculated with the dwell positions shifted by 3 mm along the ovoid lines, and also by swapping the combination of dwell positions, two dwells from the left and one from the right ovoids, instead of one from the left and two from the right. The uncertainties of the point A coordinates from these potential variations in ovoid localization was typically less than 0.1 mm. Because the algorithm compares the dwell positions within a curved tandem to the straight line drawn from the top of the ovoids plane, and a dwell position 2-2.5 cm superiorly in the tandem (i.e., dashed line shown in Figure 1), the algorithm did fail in 3-4 instances when the curved tandem approximated the straight dashed line (e.g., 15o). The calculation algorithm could be easily applied to other TPS as long as a text file specifying the applicator geometry is available and can easily account for either more or fewer dwells in either the tandem or ovoids. The algorithm could also be extended to tandem and ring applicators due to the similarities in the definition of point A [3]. We recommend that TPS manufacturers incorporate such algorithms into their software system automating both calculation of point A coordinates and insertion into the 3D plan. While this would ensure that plans use correctly defined dose points without potential human errors in the data transfer between TPS and standalone point placement applications, verification by the end user would always be necessary. Such implementation would be the first step to standardization of dose point definition for brachytherapy.

Conclusions

By utilizing the automated dose point calculation algo­rithm, points are defined consistently between treatment plans without inter-user variability that might affect dosimetric evaluation and reporting. Not only does this provide compliance with the ABS recommendation of consistent point A placement and dose reporting [1,3] but it also addresses one of the failure modes identified in HDR treatment planning [10,11,12]. Automated calculation has the potential to reduce human errors in radiotherapy planning, particularly in clinics with a single planner for HDR brachytherapy [11].

Disclosure

Authors report no conflict of interest.

References

1. Nag S, Erickson B, Thomadsen B et al. The American Bra­chy­therapy Society recommendations for high-dose-rate brachytherapy for carcinoma of the cervix. Int J Radiat Oncol Biol Phys 2000; 48: 201-211.
2. Yin G, Wang P, Lang J et al. Dosimetric study for cervix carcinoma treatment using intensity modulated radiation therapy (IMRT) compensation based on 3D intracavitary brachytherapy technique. J Contemp Brachytherapy 2016; 8: 221-232.
3. Viswanathan AN, Thomadsen B, American Brachytherapy Society Cervical Cancer Recommendations Committee, American Brachytherapy Society. American Brachytherapy Society consensus guidelines for locally advanced carcinoma of the cervix. Part I: general principles. Brachytherapy 2012; 11: 33-46.
4. Pötter R, Fidarova E, Kirisits C et al. Image-guided adaptive brachytherapy for cervix carcinoma. Clin Oncol (R Coll Radiol) 2008; 20: 426-432.
5. Viswanathan AN, Erickson BA. Three-dimensional imaging in gynecologic brachytherapy: a survey of the American Brachytherapy Society. Int J Radiat Oncol Biol Phys 2010; 76: 104-109.
6. Pavamani S, D’Souza DP, Portelance L et al. Image-guided brachytherapy for cervical cancer: a Canadian Brachytherapy Group survey. Brachytherapy 2011; 10: 345-351.
7. Pötter R, Haie-Meder C, Van Limbergen E et al. Recommendations from gynaecological (GYN) GEC ESTRO working group (II): concepts and terms in 3D image-based treatment planning in cervix cancer brachytherapy-3D dose volume parameters and aspects of 3D image-based anatomy, radiation physics, radiobiology. Radiother Oncol J Eur Soc Ther Radiol Oncol 2006; 78: 67-77.
8. Grover S, Harkenrider MM, Cho LP et al. Image Guided Cervical Brachytherapy: 2014 Survey of the American Bra­chy­therapy Society. Int J Radiat Oncol Biol Phys 2016; 94: 598-604.
9. Otter S, Franklin A, Ajaz M et al. Improving the efficiency of image guided brachytherapy in cervical cancer. J Contemp Brachytherapy 2016; 8: 557-565.
10. Wilkinson DA, Kolar MD. Failure modes and effects analysis applied to high-dose-rate brachytherapy treatment planning. Brachytherapy 2013; 12: 382-386.
11. Damato AL, Devlin PM, Bhagwat MS et al. Independent brachytherapy plan verification software: improving efficacy and efficiency. Radiother Oncol 2014; 113: 420-424.
12. Mayadev J, Dieterich S, Harse R et al. A failure modes and effects analysis study for gynecologic high-dose-rate bra­chy­therapy. Brachytherapy 2015; 14: 866-875.
13. Wadi-Ramahi S, Alnajjar W, Mahmood R et al. Failure modes and effects analysis in image-guided high-dose-rate bra­chy­therapy: Quality control optimization to reduce errors in treatment volume. Brachytherapy 2016; 15: 669-678.
14. Damato AL, Lee LJ, Bhagwat MS et al. Redesign of process map to increase efficiency: Reducing procedure time in cervical cancer brachytherapy. Brachytherapy 2015; 14: 471-480.