3/2019
vol. 11
Original paper
Hybrid optimization based on noncoplanar needles for brachytherapy dose planning
Xiaodong Ma, Zhiyong Yang, Shan Jiang, Guobin Zhang, Bin Huo, Shude Chai
J Contemp Brachytherapy 2019; 11, 3: 267–279
Online publish date: 2019/06/23
Purpose
Malignant cancer is the leading cause of death in both men and women in the world [1]. By implanting
Iodine125 (125I) seeds into the target, the radioactive dose can be adapted to the shape of target volume. Since the dose intensity is inversely proportional to the square of the distance [2], the dose accumulated in normal tissues around the target is minimized. A number of studies on 125I seeds implantation for prostate cancer has achieved acceptable results [3,4,5,6]. Thus, 125I seeds implantation, with its satisfactory effectiveness, has been accepted as one of the methods of prostate cancer treatment.
With regards to 125I seeds implantation surgery, the preoperative dose planning is the key to successful operation. Dose planning in highdoserate brachytherapy is developing rapidly and at the moment, there are many wellknown inverse optimization procedures [7,8,9]. However, the inverse optimization method is rarely applied to 125I seed brachytherapy (lowdoserate brachytherapy), and dose planning is mostly completed manually by the physician.
In recent years, a number of researches on 125I seed brachytherapy based on the coplanar template have been published [10,11,12]. While, in clinical practice, a typical coplanar template is designed with parallel needle guiding holes, evenly arranged at 5 mm interval in a 13 by 13 matrix or a 17 by 17 matrix. Therefore, the needle puncture path is uniquely fixed. Considering dose distribution to the target, it is inevitable that occasionally the puncture path can be blocked by bones or other critical tissues, when the template is used in brachytherapy. In addition, the fixed path leaves a narrow space between the needle and other critical tissues or organs, which increases the risk of injury during the procedure. Some researches show that multiple independent templates can intensify flexibility in patient specific treatment planning and receive better treatment planning results [13]. Currently, there is no commercial treatment planning system available, which allows noncoplanar needle insertion and automatic treatment plans. Therefore, a novel method for an arrangement of needle with adequate dose distribution is urgently needed, combining flexible noncoplanar needle with effective dose planning.
In the last two decades, due to the development of computer technology, many optimization methods for dose planning in seed brachytherapy have been proposed [14,15,16]. However, the inverse planning technique is still investigated. Comparing with traditional optimization methods such as the Manchester [17] and Paris systems [18], the inverse treatment planning could achieve a balance between dose distribution and optimization time, with more accurate dose planning. The inverse optimization algorithm is different from the traditional dose planning. Firstly, it needs to set the goal to be achieved. Then, the objective function based on the constraint conditions should be determined. Finally, the iterative search method is used to find the optimal solution.
Dose volume histogram (DVH) curve has a weak connection with the geometries of organs. In numerous studies, the optimization is established by setting a series of constraints on DVH or matching the DVH to referenced ones [19]. Xing et al. established a dose optimization framework for incorporating metabolic information from functional imaging modalities into the intensitymodulated radiation therapy with inverse planning process [20]. Cotrutz’s group proposed a new inverse planning scheme with voxeldependent importance factors [21], in which the planning procedure can be divided into two steps. After the completion of a conventional trialanderror inverse planning procedure, the dose interval is identified to locate the voxels that receive the dose in the selected range. Mountris et al. proposed an inverse planning approach that can improve planning quality by combining DVH without compromising the overall execution times [22]. However, most DVHbased optimization algorithms are mostly used in external radiotherapy, and only a small number of specialists consider the evaluation of DVH with the automatic optimization of seeds distribution.
The aim of this article is to propose a new optimization method based on noncoplanar needle, which combines an inverse optimization algorithm with a DVH evaluation method called ‘the hybrid inverse optimization method’.
Material and methods
Planning procedure
A new tool was proposed to help the physician to acquire the treatment plan. The general workflow is presented in Figure 1.
Firstly, the geometric boundary of the tumor based on manual segmentation and the planning target volume (PTV) was obtained by the physician, and the skeleton based on computerized tomography (CT) value was automatically extracted. Secondly, the puncture needles were assigned and adjusted manually by the physician to cover the target. Since the noncoplanar needle was used in this process, there was no need to consider the direction of needle. Thirdly, the hybrid inverse optimization algorithm was applied to search for optimal seed distribution for each seed number (minimum distance between two seeds positions was 5 mm). Finally, the DVH evaluation method was used to determine the optimal seeds distribution.
Noncoplanar needle assignment method
The coplanar template is commonly used to install the puncture needle in ultrasonicguided brachytherapy of prostate cancer, and achieves good results, since there are no barriers for needle puncturing. The template ensures the exactly perpendicular location of puncture needle, which guarantees high precision and helps physician to perform the intraprocedure rapidly.
Dose distribution is often unable to meet the clinical demand, when using a coplanar template to treat a multifocal metastasis of lung cancer. Given that, the parallel needles in a coplanar template are difficult to cover all areas and dose distribution is not able to reach the prescription dose. To avoid bones and ensure dose coverage when the tumor is vertebra metastasis or in the shade of bones, we developed a dose planning method based on noncoplanar needle, with each needle to be adjusted in any direction.
The position adjustment of noncoplanar needle is illustrated in Figure 2. The initial needle was in the transverse plane (P1), parallel to Yaxis. Supposing, the length of the needle is l, its center point is C, and its tip point is T.
Firstly, the needle center was placed at C (x, y, z) point. Secondly, the initial needle was rotated around C by to
obtain an intermediate needle. Then, a new plane (defined as P2) perpendicular to P1 was created, containing the intermediate needle. Finally, the intermediate needle was rotated around C by to attain the final needle in P2.
Therefore, the needle center C and the needle tip T were chosen to control the position of the needle. The relative coordinates of needle tip are shown in Figure 3, which are respectively marked as oxyz, o1uvw, o2u′v′w′, and o3u″v″w″. The transformation process was described as follows. Firstly, the needle was translated along Xaxis with x, Yaxis with y, and Zaxis with z to acquire o1uvw. Secondly, the needle was rotated around the Waxis by to get o2u′v′w′. Finally, the needle was rotate around the U′axis by to reach o3u″v″w″. Consequently, the final position of the needle tip T3 can be obtained from the below formula [23]:
where T0 is the translation matrix of the puncture needle relative to the coordinate system oxyz. T1 and T2 are the rotation matrices of the puncture needle with respect to the Waxis and Uaxis, respectively.
Radiation dose calculation model
The 125I (6711) seed involved in this research was shaped in cylinder, with a length of 4.5 mm and diameter of 0.8 mm. The seed was wrapped with 0.05 mm titanium, and iodide ion adsorbed with radiation activity of 0.7 mCi as internal. The dose rate was calculated using the TG43 U1 protocol [24] published by AAPM, as shown in equation (2):
where Sk is the air kerma strength of the source, is the dose rate constant, G(r, θ) is geometry factor, gL (r) is radial dose function, and F(r, θ) is anisotropy function. r represents the distance from the center of the active source to the point of interest, r0 denotes the reference distance, which is specified to be 1 cm in the protocol, and θ signifies the polar angle specifying the pointofinterest relative to the source longitudinal axis. θ0 denotes the source transverse plane and is specified to be 90°.
Dose rate calculation model for spatial multiple seeds was based on scalar superposition, which means that the dose rate at a point in space is equal to the sum of the dose rate of all the seeds at that point. Suppose there was a point P(x, y, z) in the dose space, regardless of tissue heterogeneity and interseed attenuation, the total dose of n seeds can be calculated as equation (3):
The optimization framework
The purpose of dose planning for conventional brachytherapy is to obtain an ideal dose distribution. On the one hand, there is sufficient radiation dose in the target to terminate the tumor cells. On the other hand, the organs at risk (OARs) are prevented from being damaged by radiation dose. In this research, the hybrid inverse optimization algorithm was proposed to better achieve the purpose, in which an inverse optimization algorithm and a DVH evaluation method were involved. Firstly, the number (denoted by Ne) of seeds needed was calculated based on the Anderson nomogram [25], according to the geometry of the target area and the type of seed. Since the shape of the target was irregular in the actual situation, more seeds were needed to provide sufficient radiation dose. The upper limit of the number of seeds can be set to a large value, but this leads to increased computational time. Moreover, multiple simulations demonstrated that the maximum number of seeds after optimization was 1.3 Ne.
At that time, the hybrid inverse optimization method was carried out to find the optimal seed position. The number of seeds Ns, ranging from Ne to 1.3 Ne, was inputted to run the optimization program and then, each optimal solution was obtained. Finally, DVH evaluation method was included to determine the optimal result.
Objective function of the inverse optimization algorithm
Seed implantation aims to deliver the dose to cover PTV as much as possible. At the same time, organs at risk (OARs) should be irradiated minimally. Therefore, the objective function needs to consider all these factors by setting the weight coefficient of each organ’s objective function. For this purpose, a multiple objective function F with penalty value of each organ was given by equation (4):
where T refers to the target volume, p denotes to the prescription dose, and ω refers to the weighting factor of the target volume or OAR. n refers to the number of interest points in the target, nj refers to the number of interest points in the OAR, and nOAR refers to the number of OARs. In addition, Di refers to the dose of each interest point, and Dp denotes to the prescription dose of the target or the tolerated dose of OAR.
In terms of the contradictory goals that minimize the dose variation, the prescription on the surface of the tumor PTV will also minimize the dose to the OARs. All these factors were involved in equation (4), so that physicians can choose the best prescription dose for the target and the tolerated dose for OARs, according to their clinical experience. Then, the minimal F(M) can be determined through the inverse optimization algorithm. To determine the weighting factor and to avoid large number of trials, a multicriteria optimization (MCO) was involved [26].
Optimization by the inverse optimization algorithm
The inverse optimization algorithm used in this research is the simulated annealing [27], which is analogous to a crystal seeking its lowest energy, as it slowly decreases from a high initial value. As a kind of stochastic algorithm, it is suitable for an optimal dose allocation. This research improves the determination of the initial value in simulated annealing. In the simulated annealing algorithm, in order to obtain the global optimal solution, a large initial value needs to be established. However, the number of seeds can be determined based on the tumor volume (nomogram), thus it is easy to determine the initial value. In the iteration process of the seed distribution, each layout scheme has an evaluation function and a corresponding constraint condition. The inverse optimization algorithm determines whether to accept a new distribution according to the evaluation function, so that the evaluation function in the distribution model would be expressed as a quantitative index. The evaluation function was in fact the objective function we mentioned before. According to equation (4), F(M) > 0. Furthermore, it is clear that the smaller the value of F(M), the better the current state. According to the Metropolis criterion [28], F(M)i was the evaluation function of the current state of i. Then, a disturbance was made to the state i, and a new state j with its evaluation function F(M)j can be obtained. If the evaluation function F(M)j < F(M)i and j for the optimal condition, the radioactive seeds distribution will jump to state j from state i. If F(M)j ≥ F(M)i, with distribution state of j inferior to the state of i, the transition will be in accordance with a certain probability of Pi, which was called the acceptance probability:
where κB denotes the Boltzmann constant, n is the nth iteration, Tn is the parameter that controls the reduction of the initial value, and controls the reduction rate of the initial value. That is, the slower the rate, the more accurate the inverse optimization algorithm results, and vice versa.
In this algorithm, there are two ends conditions. The first one was that the initial value Tn decreased to zero through the equation (6). The initial value should be high enough to prevent the algorithm from falling into the trap of a local optimal solution. The second end condition was the equation (7). The algorithm terminated when the change in the result of two adjacent iterations was less than a certain value η. The smaller the value of η, the more accurate the optimization result, but the optimization time will be longer.
The DVH evaluation method
Clinically, DVH contains a number of important indicators that reflect the effect of treatment. DVH describes the cumulative dose received by a certain percentage of tumor volume or organ volume. In DVH, the Yaxis represents the volume (V) and the Xaxis represents the dose (D). Vx means the fraction of volume that received x% of the prescription dose, and Dx means the dose received by the fraction x of the volume. The relationship between tumor or organ and dose distribution can be observed from the curve. However, DVH cannot reflect the real geometric information of the tumor or organ. In addition, it cannot reflect the number and distribution of seeds. Hence, we can firstly build a reference DVH model and then, the result of the inverse optimization algorithm can be evaluated by comparing the practical DVH curve with the model.
In the seed implantation surgery for thoracoabdominal tumors such as lung and liver, there may be guidelines missing to evaluate the DVH. Dose irradiating a small volume (D0.1cc, D2cc), which may be correlated with the side effects, are generally reported. The dose homogeneity index (DHI) [29] is the target receiving a dose between 100% and 150% of the prescription dose, and that is DHI = (V100 – V150)/V100. To ensure the curative effect of tumor, the dose planning target requires that 90% of target volume receives at least 100% of the prescription dose (V100 ≥ 90%). In addition, V150 < 50% and V200 < 25% were set to avoid the delivery of excessive dose to the patient. In the DVH evaluation model, as shown in Figure 4, the red line represents the reference and the blue line denotes one of the current results of the inverse optimization algorithm. For PTV, V100, V150, and V200 were chosen as the turning point and four areas, V0V100 (S1), V100V150 (S2), V150V200(S3), and V200V300 (S4) were defined. For the target, the blue DVH was better than the red DVH, as shown in Figure 4A. S1 and S2 were above the reference line, showing that the results were better. The result’s quality is calculated from the size of S1 and S2 areas, where for increasing S1 and decreasing S2 areas, the quality increases. S3 and S4 were below the reference line, showing that the results were better. Similarly, the result’s quality is quantified from the size of S3 and S4 areas, where for increasing S3 and S4 areas, the quality increases. Therefore, we can define the quality of the PTV of current DVH:
where 1, 2, 3, and 4 are the weighting factors of each areas.
For OARs, in Figure 4B, S’1, S’2, and S’3 represent the areas where the DVH after optimization was better than the reference DVH, and the larger the area, the better the optimized DVH. The function for evaluating OAR was given as follows:
where 1, 2, and 3 are the weighting factors of each area.
Therefore, the evaluation equation, combining PTV and OARs, can be defined as follows:
where fi is the optimal result of number i seeds in the inverse optimization algorithm (the larger, the better). δ1 and δ2 are the weighting factors of each OAR, which are small numbers.
The reference DVH (based on noncoplanar needle) was obtained from one of the cases screened by physician that met the dose requirements. Comparing with the reference DVH, the physician managed to create a rank of selected DVHs manually, based on clinical experience. Then the best weighting factors of DVH evaluation function were selected to match the manual ranking results. After that, the weighting factors were incorporated into DVH evaluation function. Finally, all the results acquired from the inverse optimization algorithm were compared to evaluate a max fi as the ultimate optimization result. The final results were not influenced if the reference DVH was changed or if another physician was estimating the ranking. Firstly, the reference DVH should meet the requirements of V100 ≥ 90%, V150 < 50%, and V200 < 25%. Secondly, when compared with the reference DVH, the sizes of S1, S2, S3, and S4 should be considered. Therefore, if the reference DVH was changed and if another physician evaluated the ranking, the DVH of the optimal result will not be influenced by the change, but the sizes of S1, S2, S3, and S4 will be different.
Experimental platform construction
In order to implement our new planning method, a brachytherapy treatment planning system (BTPS) was designed using VTK, ITK [30,31], and the Microsoft Visual Studio 2010 software (Figure 5).
Case studies
All the images in the study were CT images (DICOM 3.0) of 5 mm spacing between slices, and the voxel resolution was 0.70 × 0.70 × 5.00 mm3. The first test was designed to examine the performance of the hybrid inverse optimization algorithm and to compare the difference between noncoplanar needles and coplanar needle planning. BTPS was run on a Dell computer with Intel Core i7 3.60 GHz CPU and 8 GB RAM. The second test was carried out to investigate the practical application of the new technique. In order to improve the computing speed of the method, the minimum cube of the wrapped target was extracted firstly and then extended 50 mm in every direction, as the dose optimization area. Additionally, the calculation range of single seed model was set to 50 × 50 × 50 mm3 by balancing the accuracy of dose and the speed of program. Moreover, the noncoplanar needle’s spacing was between 5 mm and 10 mm, and the needles with seeds were placed with 5 mm intervals. Furthermore, the optimization time was also recorded and included: finding the number of seeds based on nomogram (Ne), repeating the optimization by increasing the number of seeds (in the range from Ne to 1.3 Ne), calculating all DVHs and finding the best solution, based on the reference DVH.
Results
Lung adenocarcinoma case to test the method
Left lung adenocarcinoma case (Figure 6) with 21.22 cm3 of PTV was involved in this test. The spinal cord was chosen as OAR, and the type of seed was 6711 125I with strength of 0.7 mCi. The prescription dose for the tumor was 120 Gy. For OAR, the tolerate doses of the thoracic spine was 45 Gy [32]. In addition, it can be seen that the clinical data of lung at D2cc < 85 Gy had a lower risk in complication, so it was considered as the lowrisk limits [33]. The total number of the voxels in the series of DICOM images and the extracted region for this case were 499 × 499 × 80 and 99 × 99 × 40, respectively. 16 noncoplanar needles were initially assigned, and the optimization algorithm took 9.8 seconds, with 24 seeds and 13 needles. Particularly, the noncoplanar needles were assigned based on the experience of the physician, taking about 10 minutes. In order to accurately assign the initial needles, the minimum distance between two needles of 10 mm and the organs at risk avoidance should be considered. The initial needles should be assigned to cover the target. Therefore, the number of initial needles can be determined. Figure 7 shows the DVHs curves obtained from the optimization of the seeds’ distribution for different number of seeds. In addition, Table 1 and Table 2 shows the important parameters of the corresponding DVH curves. As can be seen from the graph, the curves with a small seed number were below the reference line. As the number of seeds increased, the DVH curve shifted upward. The curve (V100 < 90%, V150 < 50%, V200 < 25%) below the reference line demonstrated that some parts of the tumor did not receive enough dose. As for the top black line (V100 > 90%, V150 > 50%, V200 > 25%), most of the tumor received the prescription dose. The vast majority of the tumor received a sufficient dose; however, a large part was exposed to excessive amount of radiation. At this stage, increasing the number of seeds to allow more areas to receive sufficient doses would lead to more regional overdose, which should be avoided.
When the numbers of seeds are 22, 23, and 24, all the dose distribution generated by the inverse optimization algorithm can meet the clinical demand. In our study, through the DVH evaluation, 24 was chosen as the best seeds number. Clinically, on the premise of satisfying dose, the number of seeds should be as small as possible. Also, in order to reduce the injury of the puncture process and relieve the pain of the patient, the number of puncture needles should be the fewest possible as well, in view of meeting the dose requirements.
The coplanar puncture needles and noncoplanar puncture needles assigned to cover the target by physical therapists is presented in Figure 8. It was clear to see that the noncoplanar needle can easily avoid the bone to cover the target.
For coplanar template, 11 needles and 23 seeds were used to reach the best possible dose. However, the dose distribution did not reach the clinical requirement of terminating the tumor cells due to blocking rib. The dose parameters are described in Table 3.
Malignant tumor of spine case to test the technique
The tumor was 94.67 cm3 of PTV, and 63 seeds and 22 needles were used in the actual procedure. The simulation experiment was conducted under the same conditions as the surgery. The spinal cord was selected as organ at risk, and type 6711 125I seeds with strength of 0.5 mCi were used. The prescription dose for the tumor was 120 Gy. The total number of the voxels in series of DICOM images and the extracted region for this case were 499 × 499 × 170 and 99 × 99 × 90, respectively. According to nomogram, it was estimated that the number of seeds was 58, while the optimal result of the algorithm was actually 65 seeds. In addition, the procedure started with 30 needles but finally, 22 noncoplanar needles were used. The optimization algorithm took 20.5 seconds to obtain the ultimate optimal result.
As shown in Figure 9A, the tumor invaded the spinal cord, so the needle arrangement needed to avoid injuries to the spinal cord. The distribution of the noncoplanar needles is presented in Figure 9B, with spacing between 5 and 10 mm. Furthermore, it should be ensured that the needles cover all areas of the tumor and avoid interference. Figure 9C and Figure 9D are the dose coverage visualization of PTV by isodose curves and surfaces, from which we can directly screen the dose cold and hot points in the target area. As can be seen, the dose was able to meet the clinical requirements. In order to realize the clinical implementation of the new dose planning method, a threedimensional template was designed to guide and hold noncoplanar needles as depicted in Figure 9E and Figure 9F.
The dose parameters for the best optimal results are shown in Figure 10 and Table 4, in which the clinical dose requirements for PTV was met: V100 = 92.24 (V100 > 90%), V150 = 45.09% (V150 < 50%), and V200 = 14.63% (V200 < 25%). It should be noted that the tumor has greatly eroded the spine and the maximum dose to the spinal cord was 96 Gy. While in the case of high local dose, D90 = 1.3 Gy, V100 = 2.48%, and V200 = 2.06% were acceptable for clinical treatment.
In the optimization, the progress of the objective function value against the number of iterations is presented in Figure 11. The curve was smoothed to show the trend of F(M) more clearly. The result demonstrated that there was a local minima in the optimization process, and the algorithm can escape from the local minima and converge to the global minimum.
Comparison with other algorithms
To better evaluate the proposed method in this research, we compared the result of malignant tumor of spine case with existing algorithms [14,34,35,36,37]. Most of the current inverse optimization algorithms were developed for prostate cancer. D’Souza et al. used a mixedinteger linear programming and a branchandbound algorithm to generate treatment plans [14]. An inverse optimization planning process utilizing a biologicallybased objective was proposed [34]. The IPSA inverse planning algorithm was modified to include multiple dose matrices for the calculation of dose from different sources, and a selection algorithm was implemented to allow for the swapping of source type at any given source position [35]. The mixed integer programming was used to propose a volumebased objective function for dose optimization, which allowed for minimization of the number of under or overdosed voxels in selected structures [36]. McGeachy et al. [37] used a simple genetic algorithm for dose optimization. The comparison results are shown in Table 5, including V100, V150, and the optimization time of each optimization algorithm [37].
In addition, we made a comparison using different size of single seed model to test the calculation speed and accuracy. As such, single seed models of 20 × 20 × 20 mm3, 40 × 40 × 40 mm3, and 50 × 50 × 50 mm3 were involved. The results show that it takes 24 seconds for 20 × 20 × 20 mm3, 69 seconds for 40 × 40 × 40 mm3, and 2532 seconds for 50 × 50 × 50 mm3 in the optimization procedure. With reference to that, it also shows, the error is larger with smaller calculation range. The calculation error is more than 8% for 20 × 20 × 20 mm3, more than 5% for 40 × 40 × 40 mm3, and less than 5% for 50 × 50 × 50 mm3.
Discussion
In this paper, a hybrid inverse optimization method based on noncoplanar needles for CTguided brachytherapy was proposed. The planning tool developed is capable of avoiding bones, decreasing needles and seeds, and automatically obtaining an optimal seed distribution.
From the result of the first case, it can be found that the noncoplanar needle and the hybrid inverse optimization method performed appropriately on the dose distribution in thoracoabdominal tumor. The noncoplanar needle has great advantages over coplanar, for that it can puncture at any angle, which helps to cover the whole target. In addition, as shown in Figure 8, noncoplanar needles can decrease the injury of normal tissues and avoid bones. As shown in Figure 8A, the coplanar needles were blocked by ribs, which increased the difficulty to reach satisfactory dosage. In contrast, as shown in Figure 8B, the noncoplanar needles were able to avoid the ribs and organs at risk. Additionally, the noncoplanar needles could reach almost any places in the chest and obtain a perfect dose distribution, something that coplanar needle could not achieve.
The application of the algorithm was verified through the second case, in which the treatment used 63 seeds and 22 needles, and the optimized results were 65 seeds and 22 needles. The inverse optimization algorithm was able to handle large tumors in real time for only 20.5 seconds and achieved a clinically accepted plan with 94.67 cm3 PTV. In terms of running speed, the hybrid inverse optimization method was fast enough to meet realtime dose planning requirements.
The advantages of using noncoplanar needle in the treatment planning have been discussed; however, assigning the noncoplanar needles’ position remains to be a challenge for a physician. Given that, the seed candidate position is attached to noncoplanar needles; the greater the number of needles is, the more time is required to find the optimal seed position (but better dose distribution would be obtained). Clinically, the recommended needles’ spacing is between 5 mm and 10 mm, while an excess of needles will be excluded by the hybrid inverse optimization method. Also, the needles in the same plane should be kept parallel as much as possible to prevent interference. Hence, a physician could combine the characteristics of noncoplanar needles mentioned before with clinical target coverage to improve the speed of the noncoplanar needle arrangement.
The hybrid inverse optimization method performs well in calculating the optimal seed distribution. It is based on the location of candidate seeds, which are indeed determined by the needle. However, from the running time of the two experiments, we can see that different tumor volumes need different optimization time and the larger tumor volume, the more time consuming. Owing to that, a larger target requires a greater number of needles, which produce a bigger number of seed positions candidate; the hybrid inverse optimization method would be more time consuming to find each quantity of seeds’ optimal results. We found that the speed of the hybrid inverse optimization method can be visibly improved by imputing initial geometric distribution. Through many simulation experiments, we found that the optimal seeds distribution always satisfied a rule. The seeds were densely distributed around the tumor, while the seeds distribution inside the tumor was sparse, and the distance between the seeds was at least 10 mm. Therefore, for larger targets, the initiation of seeds distribution according to this rule would greatly reduce the optimization time. In addition, a reasonable arrangement of noncoplanar needles to avoid too much dense arrangement would also contribute to time saving.
In the prostate brachytherapy, the recommended values by AAPM are V100 > 95% and V150 < 50%. In lung cancer brachytherapy, our hospital recommends V100 > 90% and V150 < 50%. From Table 5, we can see that the algorithms 1 and 3 did not meet the requirement of V100 > 95%. The other three algorithms 2, 4, and 5 met the requirement of V100 > 95%, but their V150 was more than 50%. The algorithm proposed in the current research met all the requirements of V100 > 90% and V150 < 50%. In terms of optimization time, the current research took less time (20.5 seconds) than other algorithms. Certainly, the improvement of the optimization result in target coverage and the optimization time will be the subject of future research. In addition to that, we have already finished the work of implementing the noncoplanar needles in clinical treatment. We managed to design a virtual template based on the extraction of a part of skin surface and print it on a 3D printer (Stratasys Object30pro).
As it is mentioned before, the 50 × 50 × 50 mm3 single seed model was used for optimization. But in the final dose and DVH computations, the 40 × 40 × 40 mm3 single seed model was used. According to the results, the calculation speed is faster when the single seed model is smaller, but the error is larger. When the dose calculation range is larger than 50 × 50 × 50 mm3, the calculation error will be much smaller, but the calculation time takes 24 or more minutes. We also compared the results of the dose calculation with the commercial TPS, such as VariSeed™ (Varian), Panther TPS (Prowess). The results demonstrated that the dose calculation using the 40 × 40 × 40 mm3 single seed model was consistent with the commercial TPS. The dose calculation error using the 50 × 50 × 50 mm3 single seed model was less than 5%. Therefore, to improve the optimization accuracy, we used the 50 × 50 × 50 mm3 single seed model.
Unquestionably, there was a limitation in this algorithm. The dose calculation formula was defined based on TG43 line source (2D) formalism, in which all seeds were assumed to be a part of the needle with proper orientations. However, in clinical surgery, the initial seeds orientations are usually unknown, which would result in an error in the actual dose distribution and the optimized results. Nevertheless, currently, there is no reasonable method to solve this problem, and in a future research, this issue will be thoroughly investigated. In addition, this study is a preliminary attempt in the optimization of noncoplanar needles, and no Monte Carlo or other advanced dose calculation methods were applied. In future, we will incorporate advanced algorithms, such as Monte Carlo to our proposed method. Furthermore, we hope to establish a technique to evaluate whether the distribution of the initial needle is reasonable. Because the initial needle is assigned by the physician according to clinical experience, a better seed distribution may be required. We also plan to incorporate a homogeneity index into the optimization algorithm to effectively avoid hot and cold spots, and to improve the optimization speed.
Conclusions
A hybrid inverse optimization method based on noncoplanar needles for CT‑guided brachytherapy was proposed. The tool was able to arrange the noncoplanar needles to cover the tumor and avoiding organs at risk. In addition, the hybrid inverse optimization method has achieved a clinically acceptable dose distributions, which was effective for local tumor control. The hybrid inverse optimization last for 20.5 seconds in the tumor’s volume of 94.67 cm3, which is fast enough to be applied in practice.
We will continue our work on optimizing the hybrid inverse optimization method, which will be appropriate for multiple organs. We want to apply our method to intraoperation’s optimization by the detection of the real puncture needles and run the optimization based on the detected needles.
Acknowledgements
We gratefully acknowledge the research team at the Center for Advanced Mechanisms and Robotics, Tianjin University for their technical assistance. Department of Oncology, the Second Hospital of Tianjin Medical University provided us with great support of clinical experience.
Ethical approval
This article does not contain any studies with human participants or animals performed by any of the authors.
Funding
This study was funded by the National Natural Science Foundation of China (grant No. 51775368, grant No. 5171101938) and the Science and Technology Planning Project of Guangdong Province, China (Grant No. 2017B020210004).
Disclosure
Authors report no conflict of interest.
References
1. Torre LA, Bray F, Siegel RL et al. A 2015 Global cancer statistics. CA Cancer J Clin 2015; 65: 87108.
2. Nath R, Anderson LL, Luxton G et al. Dosimetry of interstitial brachytherapy sources: recommendations of the aapm radiation therapy committee task group no. 43. american association of physicists in medicine. Med Phys 1999; 26: 25142520.
3. Yu Y, Anderson LL, Li Z et al. Permanent prostate seed implant brachytherapy: report of the American Association of Physicists in Medicine Task Group No. 64. Med Phys 1999; 26: 20542076.
4. Li XA, Wang JZ, Stewart RD et al. Dose escalation in permanent brachytherapy for prostate cancer: dosimetric and biological considerations. Phys Med Biol 2003; 48: 27532765.
5. Wust P, von Borczyskowski DW, Henkel T et al. Clinical and physical determinants for toxicity of 125i seed prostate brachytherapy. Radiother Oncol 2004; 73: 3948.
6. Ahmad F, Cher M, Dalmia P et al. Comparison of I125 permanent prostate brachytherapy template guided implants versus interactive freehand technique. Cancer J 2001; 7: 550.
7. Siauw T, Cunha A, Berenson D et al. NPIP: A skew line needle configuration optimization system for HDR brachytherapy. Med Phys 2012; 39: 43394346.
8. Pouliot J, Foster W, Cunha A et al. SUFFT69: An inter‐fraction adaptive strategy for highdose rate prostate brachytherapy. Med Phys 2009; 36: 2535.
9. Cunha J, Siauw T, Garg A et al. WEABRB01: Robotic brachytherapy demonstration: implant of hdr brachytherapy needle configuration computeroptimized to avoid critical structures near the bulb of the penis. Med Phys 2012; 39: 3931.
10. Li J, Xie Q, Wang W et al. CTguided implantation of (125)I seeds (permanent brachytherapy) for metastatic tumors of the hepatic portal system: Effectiveness and safety in 13 patients. Brachytherapy 2016; 15: 224230.
11. Liu B, Zhou T, Geng J et al. Percutaneous computed tomographyguided iodine125 seeds implantation for unresectable pancreatic cancer. Indian J Cancer 2015; 52 Suppl 2: 6974.
12. Xiang Z, Li G, Liu Z et al. 125I Brachytherapy in locally advanced nonsmall cell lung cancer after progression of concurrent radiochemotherapy. Medicine 2015; 94: e2249.
13. Chaswal V, Thomadsen BR, Henderson DL. Template independent prostate brachytherapy simulation using adjoint sensitivity based treatment planning optimization technique. Brachytherapy 2010; 9: S41.
14. D’Souza W, Meyer R, Thomadsen B et al. An iterative sequential mixedinteger approach to automated prostate brachytherapy treatment plan optimization. Phys Med Biol 2001; 46: 297322.
15. Beaulieu L, Evans DA, Aubin S et al. Bypassing the learning curve in permanent seed implants using stateoftheart technology. Int J Radiat Oncol Biol Phys 2007; 67: 7177.
16. Tedgren ÅC, Carlsson GA. Specification of absorbed dose to water using modelbased dose calculation algorithms for treatment planning in brachytherapy. Phys Med Biol 2013; 58: 25612579.
17. Meredith WJ, Paterson R. Radium dosage: the Manchester system. Livingston Inc., Edinburgh 1967.
18. Pierquin B, Chassagne D, Perez R. Precis de curietherapie. Masson, Paris 1964.
19. Zarepisheh M, Shakourifar M, Trigila G et al. Amomentbased approach for DVHguided radiotherapy treatment plan optimization. Phys Med Biol 2013; 58: 18691887.
20. Xing L, Cotrutz C, Hunjan S et al. Inverse planning for functional imageguided intensitymodulated radiation therapy. Phys Med Biol 2002; 47: 35673578.
21. Cotrutz C, Xing L. Using voxeldependent importance factors for interactive DVHbased dose optimization. Phys Med Biol 2002; 47: 16591669.
22. Mountris KA, Visvikis D, Bert J. DVHbased inverse planning using Monte Carlo dosimetry for LDR prostate brachytherapy. Int J Radiat Oncol Biol Phys 2019; 103: 503510.
23. Farin G, Hansford D. Practical linear algebra: a geometry toolbox. 3rd ed. Taylor and Francis, London 2013; 199234.
24. Rivard MJ, Coursey BM, Dewerd LA et al. Update of AAPM Task Group No. 43 Report: A revised AAPM protocol for brachytherapy dose calculations. Med Phys 2004; 31: 633674.
25. Anderson LL. Spacing nomograph for interstitial implants of 125I seeds. Med Phys 1976; 3: 4851.
26. Milickovic N, Lahanas M, Papagiannopoulou M et al. Multiobjective anatomybased dose optimization for HDRbrachytherapy with constraint free deterministic algorithms. Phys Med Biol 2002; 47: 22632280.
27. Sloboda RS. Optimization of brachytherapy dose distribution by simulated annealing. Med Phys 1992; 19: 955964.
28. Hastings WK. Monte Carlo Sampling Methods Using Markov Chains and Their Applications. Biometrika 1970; 57: 97109.
29. Kataria T, Sharma K, Subramani V et al. Homogeneity Index: An objective tool for assessment of conformal radiation treatments. J Med Phys 2012; 37: 207213.
30. Hanwell M, Martin K, Chaudhary A et al. The Visualization Toolkit (VTK): Rewriting the rendering code for modern graphics cards. SoftwareX 2015; 12: 912.
31. Yoo TS, Ackerman MJ, Lorensen WE et al. Engineering and algorithm design for an image processing Api: a technical report on ITK – the Insight Toolkit. Stud Health Technol Inform 2002; 85: 586592.
32. Cefaro GA, Genovesi D, Perez CA. Delineating organs at risk in radiation therapy. Springer Milan 2013; 4973.
33. Murakami N, Kobayashi K, Nakamura S et al. A total eqd2greater than 85 gy for trachea and main bronchus d2ccbeing associated with severe late complications after definitive endobronchial brachytherapy. J Contemp Brachytherapy 2015; 7: 363368.
34. Haworth A, Mears C, Betts JM et al. A radiobiologybased inverse treatment planning method for optimisation of permanent l125 prostate implants in focal brachytherapy. Phys Med Biol 2016; 61: 430444.
35. Cunha JAM, Pickett B, Pouliot J. Inverse planning optimization for hybrid prostate permanentseed implant brachytherapy plans using two source strengths. J Appl Clin Med Phys 2010; 11: 3096.
36. Bednarz G, Michalski D, Anne PR et al. Inverse treatment planning using volumebased objective functions. Phys Med Biol 2004; 49: 25032514.
37. McGeachy P, Madamesila J, Beauchamp A et al. An opensource genetic algorithm for determining optimal seed distributions for lowdoserate prostate brachytherapy. Brachytherapy 2015; 14: 692702.
Copyright: © 2019 Termedia Sp. z o. o. This is an Open Access article distributed under the terms of the Creative Commons AttributionNonCommercialShareAlike 4.0 International (CC BYNCSA 4.0) License ( http://creativecommons.org/licenses/byncsa/4.0/), allowing third parties to copy and redistribute the material in any medium or format and to remix, transform, and build upon the material, provided the original work is properly cited and states its license.

