An automated dose verification software for brachytherapy

Purpose

Brachytherapy is characterized by a high-dose near the source and a steep dose falling away from the source. It can protect normal tissue while increasing the target dose compared with external beam radiotherapy (EBRT). Brachytherapy is commonly used in the treatment of cervical, prostate, and head and neck cancers [1,2,3]. Because of the small number of fractions and the large dose per fraction, it is difficult to assess and compensate for the consequences if there is a deviation between the planned dose and delivered dose [4]. Various regulatory agencies recommend to independently verify the dose calculated by the treatment planning system (TPS) prior to treatment [5,6,7].
Dose verification software packages have been reported in papers. Cohen et al. [8] described an independent dose-to-point calculation software for the verification of brachytherapy treatment planning. The software takes the data file used by the high-dose-rate (HDR) remote afterloader as primary input, and an additional user input is required for template-based treatment plans. Lachaine et al. [9] described an in-house software used for verifying HDR dosimetry; the software adds negligible additional waiting time for the patient, regardless of the number of applicators, paths of the applicators, or complexity of the dwell times and positions. Carmona et al. [10] and Safian et al. [11] presented spreadsheet-based software packages to verify the dose of an individual brachytherapy plan. However, only dose points calculations are performed in these software packages.
There are also some commercial software packages that allow the user to perform a fully independent calculation check of their brachytherapy treatment plan, such as BrachyCheck (ROS, California, USA) and DIAMOND for Brachytherapy (PTW, Freiburg, Germany). The plan information of BrachyCheck can be imported via a DICOM RT plan file or entered manually. BrachyCheck also allows the creation of additional pseudo-structures for inclusion in the DVH computation; however, this software does not provide information on the point dose. DIAMOND quickly performs point dose calculations and generates three-dimensional (3D) dose that could be stored in a DICOM RT dose file, but the RT dose file have to be exported to another software in case of further 3D analysis.
Currently, with modern medical imaging, there is a huge interest of 3D treatment planning for brachytherapy. The 3D brachytherapy technique conforms the dose to the anatomy of each patient, and takes into account both the tumor shape and the position of organs at risk (OARs). It has been shown that 3D brachytherapy significantly improves the dose distribution and overall survival [12]. However, most previously mentioned dose verification software packages are specific to point-dose verification or inconvenient to use for 3D brachytherapy dose analysis. In this article, we described an in-house software that performs 3D independent dose verifications for HDR treatment plans. The dose verification software (DVS) imported the dwell positions and dwell times from DICOM RT files and used them to calculate the dose. Dosimetric parameters and  analysis were used to verify the dose calculated by the DVS.

Material and methods

Data reading and writing

The DVS was written in C++ format. In order to increase the DVS applicability to be used in various treatment planning systems, the parameters required by the DVS, such as the coordinates of the region of interest, dwell position, and dwell time, were read from DICOM RT files. The DVS adopted a modular structure design to reduce the complexity of the software, and to simplify the design, debug, and maintenance operations. Each functional module had a corresponding class (Table 1).

Dose calculation

In our department, we use Ir-192 mHDR-v2 source (Elekta microSelectron units), Elekta microSelectron treatment unit, and Oncentra Brachy version 4.3 TPS. The TPS of this version is equipped only with the American Association of Physicists in Medicine, Task Group 43 (AAPM TG-43)-based dose calculation algorithm. The dose calculation algorithm of the DVS was the same as the TPS, and both used the general 2D formalism
(1),

where D· is the dose rate, r is the distance between the dose point P and the center of the source, r0 = 1 cm, θ is the angle subtended by the central axis of the source and the line connecting the center of the source and the dose point, , Sk is the air kerma strength, Λ is the dose rate constant, G is the geometry function, gL is the radial dose function, and FL is the anisotropy function.
The dosimetric data of the source in the TPS used in our clinic is based on the dataset recommended in reference [13], Λ is 1.108 cGyh-1U-1, Table 1 and Table 2 in reference [13] are gL and FL data of the TPS, respectively. In order to ensure the independence of the DVS, the dosimetric data of the source in the DVS were not obtained from the TPS but from the database recommended by the AAPM and the European Society for Radiotherapy and Oncology (ESTRO) [14]. The Λ is 1.109 cGyh-1U-1, Table 5 and Table 8 in AAPM report No. 229 are gL and FL data of the DVS, respectively. The air kerma strength of the source at a given time was manually entered by the user at the source exchange date and the DVS internally calculated the source decay. The half-life time of Ir-192 was consistent with that of the TPS; both were 73.83 days.
The mHDR-v2 is a line source. The line source needs to consider the source orientation when calculating the dose. Figure 1 is a diagram of the source coordinate system and the patient coordinate system. In this study, the orientation of the source in the patient coordinate system at the i-th dwell position was defined by the vector of this dwell position Si and the next dwell position Si+1. The orientation of the last dwell position was the same as the previous one.
In order to increase the dose calculation speed under the premise of dose calculation accuracy, we first calculated a 2D dose rate distribution table T(m, n) that was aligned with the longitudinal direction (V direction) and perpendicular to it (U direction) before the beginning of the dose calculation. This table was stored in the computer memory, and the resolution took 0.1 cm in the U and V directions. Considering the symmetry of the dose around the source, the table range was 10 cm in the U direction and from –10 cm to 10 cm in the V direction.
NA is the total number of source channels, NS is the total number of source positions in a single channel, dj,i is the dose contribution at the i-th dwell position to P in the j-th source channel, and tj,i is the dwell time of the i-th dwell position in the j-th channel. The dose at position P in the human coordinate system is calculated in the following way:
Step 1: To calculate the distance r using formula (2), where r is the distance between the dose point P and the center of the source Si;
Step 2: To calculate the angle between SiSi+1 and SiP using formulas (3), (4), and (5);
Step 3: To calculate the distance  using formula (6), where  is the distance between the dose point P and the center of the source Si in the V direction;
Step 4: To calculate the distance u using formula (7), where u is the distance between the dose point P and the center of the source Si in the U direction;
Step 5: To look up the dose distribution table T(m, n). The corresponding dose rate dj,i(u, ) is calculated using formula (8) if m * 0.1 cm ≤ u < (m + 1) * 0.1 cm and n * 0.1 cm ≤  < (n + 1) * 0.1 cm. In formula (8), w1, w2, w3, and w4 are the weights corresponding to T(m, n), T(m, n + 1), T(m + 1, n), and T(m + 1, n + 1), respectively, in the dose distribution table.
Step 6: To calculate the dose D(x, y, z) at P(x, y, z) using formula (9):
r = (x – xi)2 + (y – yi)2 + (z – zi)2 (2)
SiSi+1 = (xi+1 – xi, yi+1 – yi, zi+1 – zi) (3)
SiP = (x – xi, y – yi, z – zi) (4)
(5)
 = r cos θ (6)
u = r2 – ʋ2 (7)
dj,i(u,  ) = w1T (m, n) + w2T (m, n + 1) + w3T (m + 1, n) + w4T (m + 1, n + 1) (8)
(9).
This method has the advantage that the dose rate contribution dj,i(u,  ) is obtained by a simple bilinear interpolation of the data extracted from the computer memory without the need to repeat the calculation using formula (1).

Patient data

To illustrate the results obtained by the DVS, we randomly selected 15 patients with radical cervical cancer who had completed treatment from March 2018 to April 2018. The patient images and plan information are as follows: each patient underwent a CT scan after the insertion of a Fletcher applicator (Elekta part no. 189.730). The CT resolution was 0.1 cm × 0.1 cm × 0.3 cm. The ESTRO states that CT slice thickness ≤ 3 mm gives the best visualization [15]. The CTV and the OARs were generated using the Oncentra Brachy version 4.3 TPS. The OARs included the bladder, rectum, and sigmoid. The prescription dose was 6 Gy for the CTV. All the plans were optimized using the IPSA + graphic approach, optimized with IPSA first and then adjusted using the graphic method. The source step size was 0.25 cm.

Dose evaluation

The dose grid resolution was 0.1 cm × 0.1 cm × 0.1 cm for the TPS and the DVS. The dose calculated by the DVS was compared with the dose calculated by the TPS. We defined Dx% as the dose expressed in Gy received by x% of the total volume and Dycc as the dose expressed in Gy received by y cm3. We reported D100% and D90% for the CTV and D0.1cc, D1cc, and D2cc for the OARs [16]. The dose deviation (Dev) was calculated using formula (10):
De = 100% × DTPS – DDVS ,
DTPS (10) where DTPS and DDVS were the dosimetric parameters of the TPS and the DVS, respectively.
For EBRT,  analysis is widely used in quality assurance [17]. The DVS used 3D  analysis to verify the dose calculated by the TPS. The distance criteria were 0.1 cm, and the dose criteria were 5%. In order to speed up the value calculation, we ignored the point where the dose was less than 1% of the prescribed dose, that means if the dose at a certain point was less than 0.06 Gy, the  value was not calculated, nor it was involved in the calculation of the  pass rate.

Workflow

From the perspective of a physicist, the DVS should be as automated as possible, preferably without much manual intervention. The workflow of the DVS is as follows: 1. After the patient treatment plan is approved, the physicist saves the patient’s DICOM files into a fixed folder on the computer; 2. The DVS automatically scans this folder and reads the information from DICOM files for dose calculation; 3. The dose is calculated automatically; 4.
The dose from the TPS is automatically compared with the dose from the DVS; 5. The comparison result is loaded to a user-named txt file for future analysis; 6. The DVS automatically deletes the DICOM files in the fixed folder. This workflow has greatly improved treatment planning efficiency.

Results

For the CTV, the average difference of D100% of the 15 patients is 0.07% (minimum –0.63%, maximum 0.57%). The maximum and minimum difference of D90% is 0.53% and –0.39%, respectively, and the average value is 0.07%. For the bladder and rectum, the maximum deviation of D0.1cc, D1cc, and D2cc is less than 0.6%, and the sigmoid is 0.8%.
Figure 2 shows the  pass rate for each patient. The average  pass rate for all of the patients is 98.97%, with a minimum value of 98.24%. Although the  pass rate is high, it is still necessary to pay attention to the points where the  values are greater than 1. Taking one patient as an example, Figure 3 shows the positions that failed the  verification ( value > 1). Points with  values larger than 1 are mainly concentrated in two regions: 1. The dwell positions at the top end and the connection end of the applicator; 2. Close to the wall of the applicator.

Discussion

This article reports on the implementation method and results of the DVS developed in our center. As an independent dose verification tool, the DVS can automatically calculate the 3D dose after loading the dwell position and dwell time information exported from the TPS, and then use dosimetric parameters and  analysis to verify the dose calculated by the TPS. For the 15 randomly selected patients in this study, the difference in the DVH parameters was less than 0.8% and the  pass rate was greater than 98.24%.
The parameters of the  calculation of EBRT generally adopt 0.3 cm and 3% but the parameters of this study adopted 0.1 cm and 5%. This is because there is a large dose gradient in brachytherapy. The dose calculation grid resolution for EBRT is usually 0.3 × 0.3 × 0.3 cm3, whereas it is 0.1 × 0.1 × 0.1 cm3 for brachytherapy, and dose accuracy within a 5% level is acceptable. We also collected the data of the  analysis when the criteria were 0.3 cm and 3%. The average  pass rate for all of the patients in this study was 99.42%, with a minimum value of 99.04%. This shows that, in brachytherapy, 0.1 cm and 5% is a more conservative parameter than 0.3 cm and 3%.
Although the DVS and the TPS use the same dose calculation algorithm, there are still individual points where the  value is larger than 1. Two main reasons contribute to this: 1. The dosimetric data of the source has a tiny difference between the DVS and the TPS, also the data resolution is different. The DVS used a linear interpolation for gL and a bilinear interpolation for FL in the T(m, n) calculation, and the dose rate contribution dj,i(u, ) is obtained by a bilinear interpolation of the T(m, n). The results of the DVS and the TPS after interpolation at different resolutions could vary. This reason has a particularly significant effect on the dose near the applicator since the high-dose gradient of brachytherapy; 2. The method of source orientation determination used in the DVS differs from the one used in the TPS. That is the reason why Figure 3 shows that the locations of the  values greater than 1 are at the tip end and connection end of the applicator.
This DVS can be used for various treatment planning systems. DICOM is an international standard for medical images and related information. The DVS obtains calculation parameters from DICOM files. The AAPM and the ESTRO have recommended dosimetric parameters of commonly used brachytherapy sources. For sources without dosimetric parameters, the dosimetric parameters can be simulated using the method reported in papers [18,19]. The method described in this article can be used to conduct dose verification for the brachytherapy plan.
Since the dosimetric parameters of the DVS are independent from the TPS, errors that are usually not easily detected can be found such as: accidental changes or damage of the radiation source parameter database, calibration date, and activity input error after the radioactive source is replaced. However, it should be noted that there are several uncertainties and errors in brachytherapy [20,21]. Dose verification is only one part of a complete quality assurance system. Given that, the dose verification software and the TPS use the same dose calculation algorithm. The coordinate system, dwell point position, and dwell weight also have limitations, such as failure to detect the wrong applicator reconstruction and incorrect prescriptions, among others. These errors require a physicist to carefully examine the treatment plan.

Conclusions

Assessing the accuracy of dose calculation is critical to ensure the accuracy of brachytherapy. The DVS presented in this study provides a simple and rapid tool that can perform dose verification before treatment, and is expected to reduce the chance of an error that may result in a significant detriment to a patient. The DVS can be easily modified for other treatment units and TPS. The entire verification process of this DVS does not require much human-machine interaction; thus, it improves staff work efficiency in our department.

Acknowledgements

This research has been financed by the Key Research and Development Project of Sichuan, No. 2017GZ0322 and No. 2018GZ0197.

Disclosure

Authors report no conflict of interest.

References

1. Haie-Meder C, Pötter R, Van Limbergen E et al. Recommendations from Gynaecological (GYN) GEC-ESTRO Working Group (I): concepts and terms in 3D image based 3D treatment planning in cervix cancer brachytherapy with emphasis on MRI assessment of GTV and CTV. Radiother Oncol 2005; 74: 235-245.
2. Mazeron JJ, Ardiet JM, Haie-Méder C et al. GEC-ESTRO recommendations for brachytherapy for head and neck squamous cell carcinomas. Radiother Oncol 2009; 91: 150-156.
3. Martinez AA, Demanes JC, Vargas C et al. High-dose-rate prostate brachytherapy: an excellent accelerated-hypofractionated treatment for favorable prostate cancer. Am J Clin Oncol 2010; 33: 481-488.
4. Tanderup K, Menard C, Polgar C et al. Advancements in brachytherapy. Adv Drug Deliver Rev 2017; 109: 15-25.
5. Kubo HD, Glasgow GP, Pethel TD et al. High dose-rate brachytherapy treatment delivery: report of the AAPM Radiation Therapy Committee Task Group No. 59. Med Phys 1998; 25: 375-403.
6. Erickson BA, Demanes DJ, Ibbott GS et al. American Society for Radiation Oncology (ASTRO) and American College of Radiology (ACR) practice guideline for the performance of high-dose-rate brachytherapy. Int J Radiat Oncol Biol Phys 2011; 79: 641-649.
7. Nath R, Anderson LL, Meli JA et al. Code of practice for brachytherapy physics: Report of the AAPM Radiation Therapy Committee Task Group No. 56. Med Phys 1997; 24: 1557-1598.
8. Cohen GN, Amols HI, Zaider M. An independent dose-to-point calculation program for the verification of high-dose-rate brachytherapy treatment planning. Int J Radiat Oncol Biol Phys 2000; 48: 1251-1258.
9. Lachaine ME, Gorman JC, Palisca MG. A fast, independent dose check of HDR plans. J Appl Clin Med Phys 2003; 4: 149-155.
10. Carmona V, Perez-Calatayud J, Lliso F et al. A program for the independent verification of brachytherapy planning system calculations. J Contemp Brachytherapy 2010; 2: 129-133.
11. Safian NAM, Abdullah NH, Abdullah R et al. Verification of Oncentra brachytherapy planning using independent calculation. J Phys Conf Ser 2016; 694: 012003.
12. Derks K, Steenhuijsen JLG, van den Berg HA et al. Impact of brachytherapy technique (2D versus 3D) on outcome following radiotherapy of cervical cancer. J Contemp Brachytherapy 2018; 10: 17-25.
13. Daskalov GM, Loffler E, Williamson JF. Monte Carlo-aided dosimetry of a new high dose-rate brachytherapy source. Med Phys 1998; 25: 2200-2208.
14. Perez-Calatayud J, Ballester F, Das RK et al. Dose calculation for photon-emitting brachytherapy sources with averageenergy higher than 50 keV: Report of the AAPM and ESTRO. Med Phys 2012; 39: 2904-2929.
15. Hellebust TP, Kirisits C, Berger D et al. Recommendations from Gynaecological (GYN) GEC-ESTRO Working Group: Considerations and pitfalls in commissioning and applicator reconstruction in 3D image-based treatment planning of cervix cancer brachytherapy. Radiother Oncol 2010; 96: 153-160.
16. 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 2006; 78: 67-77.
17. Low DA, Harms WB, Mutic S et al. A technique for the quantitative evaluation of dose distributions. Med Phys 1998; 25: 656-661.
18. Selvam TP, Bhola S. Technical note: EGSnrc-based dosimetric study of the BEBIG Co-60 HDR brachytherapy sources. Med Phys 2010; 37: 1365-1370.
19. Taylor RE, Rogers DW. An EGSnrc Monte Carlo-calculated database of TG-43 parameters. Med Phys 2008; 35: 4228-4241.
20. Kirisits C, Rivard M, Baltas D. Review of clinical brachytherapy uncertainties: analysis guidelines of GEC-ESTRO and the AAPM. Radiother Oncol 2014; 110: 199-212.
21. Palmer A, Bradley D, Nisbet A. Physics aspects of dose accuracy in high dose rate (HDR) brachytherapy: source dosimetry, treatment planning, equipment performance and in vivo verification techniques. J Contemp Brachytherapy 2012; 4: 81-91.