A novel greedy heuristic‐based approach to intraoperative planning for permanent prostate brachytherapy

This paper presents a greedy heuristic‐based double iteration and rectification (DIR) approach to intraoperative planning for permanent prostate brachytherapy. The DIR approach adopts a greedy seed selection (GSS) procedure to obtain a preliminary plan. In this process, the potential seeds are evaluated according to their ability to irradiate target while spare organs at risk (OARs), and their impact on dosimetric homogeneity within target volume. A flexible termination condition is developed for the GSS procedure, which guarantees sufficient dose within target volume while avoids overdosing the OARs. The preliminary treatment plan generated by the GSS procedure is further refined by the double iteration (DI) and rectification procedure. The DI procedure removes the needles containing only one seed (single seed) and implements the GSS procedure again to get a temporary plan. The DI procedure terminates until the needles number of the temporary plan does not decrease. This process is guided by constantly removing undesired part rather than imposing extra constrains. Following the DI procedure, the rectification procedure attempts to replace the remaining single seeds with the acceptable ones within the existing needles. The change of dosimetric distribution (DD) after the replacement is evaluated to determine whether to accept or to withdraw the replacement. Experimental results demonstrate that the treatment plans obtained by the DIR approach caters to all clinical considerations. Compared with currently available methods, DIR approach is faster, more reliable, and more suitable for intraoperative treatment planning in the operation room. PACS number: 87

This paper presents a greedy heuristic-based double iteration and rectification (DIR) approach to intraoperative planning for permanent prostate brachytherapy. The DIR approach adopts a greedy seed selection (GSS) procedure to obtain a preliminary plan. In this process, the potential seeds are evaluated according to their ability to irradiate target while spare organs at risk (OARs), and their impact on dosimetric homogeneity within target volume. A flexible termination condition is developed for the GSS procedure, which guarantees sufficient dose within target volume while avoids overdosing the OARs. The preliminary treatment plan generated by the GSS procedure is further refined by the double iteration (DI) and rectification procedure. The DI procedure removes the needles containing only one seed (single seed) and implements the GSS procedure again to get a temporary plan. The DI procedure terminates until the needles number of the temporary plan does not decrease. This process is guided by constantly removing undesired part rather than imposing extra constrains. Following the DI procedure, the rectification procedure attempts to replace the remaining single seeds with the acceptable ones within the existing needles. The change of dosimetric distribution (DD) after the replacement is evaluated to determine whether to accept or to withdraw the replacement. Experimental results demonstrate that the treatment plans obtained by the DIR approach caters to all clinical considerations. Compared with currently available methods, DIR approach is faster, more reliable, and more suitable for intraoperative treatment planning in the operation room.

I. INTRODUCTION
Transrectal ultrasound (TRUS)-guided permanent prostate brachytherapy (1) has become a standard cure for localized prostate cancer. (2) As the energy emitted by commonly used seeds attenuates rapidly along with the distance, the dosimetric distribution (DD) within regions of interest (ROIs) greatly depends on the positions of implanted seeds. Hence treatment plan, which determines the positions of implanted seeds, directly affects operation outcome.
The treatment plan can be designed several days prior to (preplanning) or just before immediate execution of the plan (intraoperative planning). (3) Compared with preplanning, intraoperative planning avoids the need for two separate TRUS procedures and a reproducible patient positioning. It is also reported that intraoperative planning-based implants achieve better biological therapeutic effect and lower possibility of morbidity. (4)(5)(6) Currently available optimization methods can be classified into stochastic, deterministic, and heuristic approaches. (7) Stochastic approaches are proposed by Pouliot et al., (8) Yu and Schell, (9) Yu et al., (10) and Yang et al. (11) Based on either simulated annealing (SA) or genetic algorithm (GA), these approaches can generate a feasible treatment plan quickly, and can be used for intraoperative planning.
Deterministic approaches were reported by Lee et al., (12) Lee and Zaider, (13) and D'Souza et al. (14) Treatment planning is modeled as a mixed-integer programming (MIP) problem. And the problem is solved by branch-and-bound (BB) approach. In practice, the goal is often set to search a feasible solution that satisfies all constraints, not necessarily the optimal solution. Otherwise, the optimization process would take several days to be completed. (15) Heuristic approaches were reported by Yoo et al. (15,16) and Chaswal et al. (17,18) For each step, the potential seeds are evaluated according to their ability to irradiate target volume while spare organs at risk (OARs), and the optimal seed is selected until sufficient dose is delivered to target volume. In this process, in order to prevent selected seed from congregating, an isodose surface-based constraint is used to exclude the potential seeds which are close to the selected ones. In order to limit the number of used needles, Yoo and colleagues confine the search space within existing needles when needles number reaches the predetermined threshold. Chaswal and colleagues define a penalty function for the seed requiring adding a new needle.
The approach presented in this paper is also based on greedy heuristic. In our previous work, we developed an improved seed evaluation criterion. (19) We continued the work and developed an adaptive termination condition. The termination condition guarantees sufficient target coverage while avoids overdosing OARs. Another innovation of this approach is the DI and rectification procedure, which reduces the puncture needles without degrading the quality of treatment plan. Experiments show the algorithm can generate a satisfactory treatment plan in about 30 sec. This approach provides an alternative option for intraoperative treatment planning.

II. MATERIALS AND METHODS
The flowchart of this approach is shown in Fig. 1. First, the dose value delivered to every voxel of prostate, margin, urethra, and rectum volume (PMUR) by each potential seed is calculated. Then the greedy seed selection (GSS) procedure is implemented, which yields a preliminary plan (plan P). Then plan P is refined by the double iteration (DI) procedure. The refined plan (plan R) is further modified by the rectification procedure.

A. Dose calculation
The set of potential seeds is denoted as S p , and the set of selected seeds is denoted as S s . At the beginning stage, S p contains all potential seeds, and S s is empty. The dose value (D ij ) delivered to every voxel of PMUR by each potential seed is calculated according to Eq. (1): (1) where i is the potential seed in S p and j is the voxel of PMUR; r is the distance between the positions of i and j. Detailed description of other parameters can be found in Nath et al. (20) The mean dose delivered to the organs of PMUR ( , , and ) by each potential seed is: (2) where N p , N u , N r , and N m are the voxel number of prostate, urethra, margin, and rectum volume, respectively. Given a set of selected seeds S s , the dose value delivered to PMUR is: where N s is the number of the selected seeds.

B. The GSS procedure
The flowchart of the GSS procedure is shown in Fig. 2. The GSS procedure selects (transferred from S p to S s ) one seed at a time until termination condition is satisfied. For each time, the potential seeds are assessed by the evaluation criterion (C i ). The seed with minimum value is considered as the currently optimal seed (o, o S p ). If the termination condition is not reached, o is selected, and D j ′ (j ε PMUR) is updated by adding the dose value delivered to PMUR by o (D oj , j ε PMUR); see Eq. (4)). Otherwise, the GSS procedure terminates and the seeds in S s constructs plan P. (4)

B.1 Seed evaluation criterion
It has been demonstrated that the mean dose delivered to PMUR by the potential seed reflects its ability to irradiate these organs. (15) As the objective is to deliver sufficient and uniform dose to target volume while spare OARs, C i is defined as: where Dev i is the deviation of dose value within prostate volume after the seed i is selected (6) where is the mean value of the current dose value within target volume (7)

B.2 Termination condition
The termination condition of the GSS procedure is defined based on the target coverage and OAR damage. The target coverage (P 100 ) is the percentage of target volume within which the dose value is equal to or greater than the prescription dose (D p ). The OAR damage is assessed by the percentage of urethra (U 120 ) and rectum (R 80 ) volume within which the dose value is equal to or greater than 120% and 80% D p , respectively. P 100 , U 120 , and R 80 are calculated as follows: where Θ(x) is the Heaviside function. Through experiments we found that when P 100 reaches nearly 90%, the target region receiving lower dose is closely adjacent to urethra and rectum. In this case, the gain of P 100 is at the expense of rapid increase of U 120 and R 80 . (14) In order to guarantee sufficient target coverage while avoiding overdosing urethra and rectum volume, the termination condition is defined as follows: Stage 1: when P 100 is below 95% (this threshold value is determined according to the criterion recommended by AAPM (3,21) ), the currently optimal seed (o, o ε S p ) is selected. Stage 2: when P 100 is greater than 95%, o is selected and the increase of P 100 (ΔP 100 ), U 120 (ΔU 120 ), and R 80 (ΔR 80 ) is calculated: Whether to keep or withdraw o is determined by the ratio (R) of the sum of ΔU 120 and ΔR 80 to ΔP 100 : (14) A greater value of R indicates that o will cause relatively severe OAR damage but marginal target coverage increase. In this case, o is removed and the GSS procedure terminates. Otherwise, o is kept and the GSS procedure continues.
In our work, the threshold value of R (T R ) is fixed to 10, which means one percent increase of P 100 should cause no greater than ten percentage of the total increase of U 120 and R 80 . The threshold value is suitable for the ten tested patient cases used in this paper. It could be adjusted to cater to different priorities.

C. Double iteration procedure
When the GSS procedure terminates, the selected seeds (S s ) construct plan P. Plan P contains considerable needles carrying only one seed (single seed). The DI procedure aims to refine plan P by reducing these single seeds.
The flowchart is shown in Fig. 3. Plan P is the input of the DI procedure. First, all the single seeds are removed, and D j ′ (j ε PMUR) is updated by subtracting the corresponding dose value. Then the GSS procedure is implemented, which yields a temporary plan (plan T). The needles number of plan T (Needle t ) is compared with plan P (Needle p ). If plan T uses fewer needles than plan P (Needle t < Needle p ), it is used as the input of the next iteration procedure. Otherwise, the DI procedure terminates and plan T is the output.
For each step of the GSS procedure, the optimal seed for the current DD is selected. During DI procedure, the removed single seeds may not be optimal for the updated DD, and these seeds may not be selected. If the selected seed is in the existing needles -in other words, does not require adding a new needle -it is kept and remain valid.
The efficiency of the DI procedure is testified on ten patient data. These data are also used to test the entire algorithm. As shown in Fig. 4, for the ten cases, after 3 to 5 iterations single seeds are reduced significantly, so is the number of used needles.   4. The needles numbers of the treatment plans before and after the DI procedure. The gray bar is the number of needles carrying one seed. The blue bar is the number of the needles carrying at least two seeds.

D. Rectification procedure
The rectification procedure aims to further refine plan R by replacing the remaining single seeds with the acceptable seeds in the existing needles. The rectification procedure attempts to replace one single seed at a time; it terminates when all single seeds are processed. For each time, the replacement is evaluated by the change of P 100 , U 120 , R 80 and dose nonuniformity ratio (DNR) within target volume. DNR is defined as the ratio of P 150 to P 100 . P 150 is the percentage of target volume within which the dose value is equal to or greater than 150% D p : (22) The flowchart is shown in Fig. 5. One single seed is removed, then D j ′ (j ε PMUR) is updated by subtracting the corresponding dose value. The potential seeds in the existing needles are evaluated by C i , and the optimal seed is selected. The total increase of U 120 , R 80 , and DNR (ΔURD) is calculated: where P 100 0 , U 120 0 , R 80 0 , DNR 0 , P 100 ′, U 120 ′, R 80 ′, and DNR′ are the values of P 100 , U 120 , R 80 , and DNR before and after the replacement, respectively. The criteria used to determine whether to keep or withdraw the replacement is defined as: 1. P 100 ′ should be greater than 95%. 2. ΔURD should be below the threshold (T urd ).
The value of T urd is determined by pondering between the concerns about OAR damage, DNR, and the number of puncture needles. In our work, the value of T urd is fixed to 0.1, which is suitable for the tested patient data. It could be adjusted to sustain different practical variations.

III. RESULTS
Ten patient data were used to test the efficiency of the DIR approach. Among currently available approaches, the heuristic approach is significantly faster and the deterministic approach can generate more satisfactory treatment plan. So we implemented the heuristic approach reported by Yoo et al. (YH) (15,16) and the improved deterministic (BB) approach (14) using the same patient data for comparison.
The detailed description of the patient data is listed in Table 1. The margin region is the 3 mm layer outside the PTV of prostate. The voxel grid is set to 1 × 1 × 5 mm 3 . And the voxel numbers of the PMUR are listed in Table 1. The positions of seeds and needles are the determined by the 5 mm grid template. For the three approaches, only the seeds and needles within target volume are considered as the potential seeds and needles. The numbers of potential seeds (N ps ) and needles (N pn ) are listed in Table 1. The same type of radioactive seed, 0.4 mCi 125 I, is used for the three approaches, and D p is all set to 145 Gy. All the three approaches require one to calculate D ij at the beginning stage, and the dose calculation time is also listed in Table 1.
Both DIR and YH approach are implemented with C++ language. And the BB approach is implemented with the CPLEX MIP solver of the General Algebraic Modeling System (GAMS). (23) All the three algorithms are executed on the computer with AMD Athlon (tm) II ×4 (3.00G Hz) processor. Note: we limit the computational time of BB approach to two hrs; otherwise, it may take a rather long time to complete. Table 1. Detail description of patient data: the volume (cc) and voxel numbers of prostate, margin, urethra, and rectum; the numbers of the potential needles and seeds; and the dose calculation time (sec).

A. The YH approach
The YH approach is also based on greedy heuristic. The YH approach adopts an isodose-based strategy to constrain the search space by eliminating the potential seeds which are closed to the selected ones. And the search space is further confined within the existing needles when the number of used needles reaches the threshold (determined by the size of target volume). If there are no available seeds while the target coverage is not sufficient, all the selected seeds are removed, the constraint parameters are adjusted, and the procedure starts from scratch. This process is later referred as one iteration. The seed evaluation criterion is defined as: (18) where α, β, and γ are the weighting factors, and other parameters are the same as Eq. (5). It is reported that the weighting factors reflect the concern about the protection of OARs. (16) For example, if one intends to reduce the dose value delivered to urethra volume, the value of α should be increased. At first, we set the value of α, β, and γ to 1 and found that the dose value delivered to rectum volume is acceptable, but the dose value delivered to urethra is much higher than acceptable level. So we increased the value of α in order to achieve better protection of urethra volume. The critical parameters of the generated treatment plans for Patients 1 and 2 are shown in Table 2.
Although the value of α is increased, the value of U 120 does not change as expected. The reason may be the strategy used to constrain the search space. In the process of seed selection, the search space is determined by the selected ones. If the value of α is adjusted, the selected seed of each step will be different, which results in different search space for the next step. For each step, the YH approach always selects the optimal seed from the current search space without considering the dose value within OARs, and it terminates when sufficient target coverage is achieved. Therefore, the generated treatment plans may be uncontrollable.
For the ten tested patient data, the value of α is set to 1-6 to generate six treatment plans. And the most satisfactory treatment plan is chosen for comparison.

B. The BB approach
The BB approach models the potential seed as a binary variable (B i ), and the dose value delivered to prostate, urethra, and rectum volume calculated as: (19) where N seed is the number of potential seeds. We adopt the model reported by D'Souza et al. (14) The objective function and the constraints of this model are: In Eq. (20), α, β, γ, and δ are the weighting factors, N needle is the number of used needles. In D'Souza et al., (14) the values of these weighting factors are determined through experiments. In our work, the values of α, β, γ, and δ are determined by the principle reported by Lee et al. (12)

C. Experiment results
The primary criteria for treatment planning recommended by AAPM (3) are listed in Table 3. D 90 , D 10 and D 2cc are defined as the minimum dose in the "hottest" certain percentage (90% and 10%) or certain size (2cc) of the volume. For target volume, if P 100 is greater than 95%, D 90 will be surely greater than D p . (3) The two criteria are equivalent to each other.
Apart from these parameters, we also compared DNR, U 120 , and R 80 of the generated plans. These critical parameters, along with the numbers of needles, seeds, iterations, and computational time, are all listed in Table 4. The computational time does not include the time consumed by dose calculation. The value of α of the YH approach is also listed in Table 4.
The treatment plans generated by the three approaches all cater to the planning criteria, which indicates all the three algorithms could generate feasible treatment plans. Table 3. The criteria for the DD within target volume and OARs.

Prostate Urethra Rectum
Criterion D 90 ≥D p =145 Gy (P 100 ≥95%) Because of the flexible termination condition, the DIR approach achieves better OAR protection at the expense of one or two percentage lower target coverage. However, the lowest target coverage of the DIR plans is 96.4%, and the D 90 is greater than D p which is still sufficient and feasible. Figure 8 shows DNR values of these plans. The DNR values of the DIR plans are significantly lower than that of the BB and YH plans. This indicates the DD within target volume of DIR plans is more uniform.    Figure 9 shows the computational time of the DIR and YH approaches. The computational time of the DIR approach is proportional to the size of target volume. For the largest prostate (47.1 cc), the computational time is 46.5 sec. The YH approach involves considerable iterations, which is to remove all selected seeds and start from scratch; thus it takes longer time to obtain the final solution.  Figure 10 shows the needles numbers of these plans. The DIR approach and the YH approach have similar performance on this category. In general, the BB plans use slightly fewer needles than the DIR and YH plans. However, the treatment plan for patient 8 uses far more needles than the plans obtained with the other two approaches. This indicates that the set of weighting factors are not suitable for this patient. Figure 11 shows the seed numbers of these plans. The DIR plans use slightly fewer seeds than the BB and YH plans. This coincides with the fact that the DIR plans achieve better homogeneity within target volume.   Figure 12 shows the isodose line and the selected seed of the DIR and BB plans on each slice. The 100% D p isodose line demonstrates that both plans achieve satisfactory target coverage. The 120% D p isodose line of the DIR plan has better conformity than the BB plan. And the 150% D p isodose line demonstrates that the higher dose region of the DIR plan is significantly lower than that of the BB plan.

IV. DISCUSSION & CONCLUSIONS
Compared with the YH approach, the DIR approach yields better treatment plan in a shorter time. The improvements result from the improved seed evaluation criterion, the flexible termination condition, and the DI and rectification procedure.
Compared with the BB approach, the DIR plans achieve much better homogeneity. The BB approach protects OARs more effectively. However, the most obvious obstacle of using the BB approach for intraoperative planning is the computational time. In D'Souza et al., (14) the approach is implemented using the 2D data. The computational time is 20-45 min, and it would be much longer if 3D data are used. In Yoo et al., (16) the computational time is limited to 2 hrs. The BB approach fails to generate treatment plans for two out of ten patients. Although, the hardware and the algorithm have been greatly improved since then, we assert it is very difficult to get a feasible treatment plan using the BB approach within several minutes.
Another feature of the BB approach is that the quality of generated treatment plan is greatly relied on the values of the weighting factors. In practice, the values need to be adjusted to get a satisfactory plan. As the case of Patient 8, the value of δ needs to be increased to reduce the number of used needles. Meanwhile, other parameters (such as P 100 , U 120 , and R 80 ) also need to be considered in this process. Given the computational time, this process will be rather time-consuming.
In summary, the DIR approach presented in this paper is able to generate a feasible treatment plan for prostate brachytherapy quickly. The generated treatment plan achieves satisfactory dosimetric distribution. And the DIR approach is faster and more robust than the currently available approaches. Thus, the DIR approach is potential for intraoperative planning.