Int. J. Radiation
Biol.
Phys., Vol. 36. No. 3, pp. 71 l720, 1996 Copyright 0 1996 Elsevier Science Inc. Printed in the USA. All rights resewed 03603016/96 $15.00 + .OO
PII: SO3603016(96)003653
ELSEVIER
l
Oncology
Technical Innovations and Notes OPTIMIZATION
JEAN POULIOT,
OF PERMANENT 12% PROSTATE SIMULATED ANNEALING PH.D.,*
DANIEL TREMBLAY, PH.D.,* SANTO FILICE, M.D.*
IMPLANTS
USING FAST
JEAN ROY, M.D.*
AND
*Centre Hospitalier Universitaire de QuCbec, Department of Radiation Oncology, Pavillon L’H&elDieu de QuCbec, QuCbec, Canada Purpose: Treatment planning of ultrasoundguided transperineal lx1 permanent prostatic implants is a timeang task, due to the large number of seeds used and the very large number of possible source arrangements within the target volume. The goal of this work is to develop an algorithm based on fast simulated annealing allowing consistent and automatic dose distribution optimization in permanent “‘1 prostatic implants. Methods and Materials: Fast simulated annealing is used to optimize the dose distribution by limlmg the best seed distribution through the minimization of a cost function. The cost function includes constraints on the dose at the periphery of the planned target volume and on the dose uniformity within this volume. Adjustment between peripheral dose and the dose uniformity can be achieved by varying the weight factor in the cost function. Results: Fast simulated annealing algorithm finds very good seed distributions within 20,000 iterations. The computer time needed for the optimization of a typical permanent implant involving 60 seeds and 14 needles is approximately 15 min. An addition& 5 min are necessary for isodose distribution computations and miscellaneous outputs. Conclusion: The use of fast simulated annealing allows for an efficient and rapid optimization of dose distribution. This algorithm is now routinely used at our institution in the clinical planning of “‘1 permanent transperineal prostate implants for early stage prostatic carcinoma. Copyright 0 1996 Elsevier Science Inc. Brachytherapy,
Permanent prostate implant, Treatment planning, Optimization, Fast simulated annealing.
tatic carcinoma while maintaining very low morbidity, high patient acceptance, and low cost (24, 12, 21, 31). Ultrasoundguided transperineal “‘1 permanent im
INTRODUCTION Retropubic ‘251 permanent seed implantation for localized prostatic cancer has been marked by inconsistent results due to inherent technical difficulties in source implant disposition and inadequate patient selection (3,7, 14, 16,29). In 1983, Holm et aZ.(9) described a technique by which permanent “‘1 source placement within the prostate was achieved using a transperineal percutaneous approach assisted by transrectal ultrasound. This technique was later refined and popularized by Blasko et aZ.(5). More recently, the group at Memorial SloanKettering Cancer Center described a modified method of transperineal implantation using tomodensitometryplanning (18, 24, 30). Postimplant source distribution analysis demonstrates that the transperineal approach yields better source distribution within the planned target volume as compared to the previous retropubic approach (3, 9, 10, 19, 24, 25, 27, 29). This improved
source configuration
plants for earlystage
prostatic
carcinoma
has recently
been implemented at our center. The methodology of implantation used at our institution has been thoroughly described by Blasko and Radge et aZ.(5, 22). Patients are prescribed a 160 Gy minimum peripheral dose (mPD) to the planned target volume. The custom treatment planning of permanent 12’1seed prostate implants based on 3D imaging involves the positionning of a large number of radioactive seeds within the target volume to achieve a proper dose distribution while preserving surrounding normal tissues. Although there is no conceptual problem linked to the dosimetry, the clinical implementation is complicated by the very large number of possibilities of seed positions to form an adequate distribution. In fact, the number of possibilities is factorially large, such that they cannot be explored exhaustively. For instance, the
is related to a better
distribution
of 60 seeds on 160 positions
(16 needles X
control of the technique and to planning dosimetry and has resulted in excellent short term clinical and biochemical local control in selected patients with early stage pros
10 available positions per needle) gives more than lo4 different seed arrangements. This is comparable to the to
Reprint requests to: Jean Pouliot, Centre Hospitalier Universitaire de QuBbec, Department of Radiation Oncology, Pavillon
L’H&elDieu de QuCbec, QuCbec, Canada GlR 2J6. Accepted for publication 12 July 1996. 711
712
I. J. Radiation
Oncology
0 Biology
0 Physics
tal number of gas molecules contained in the earth atmosphere and even with the most efficient of computers, it would require more than the age of the universe to test every seed distribution for the best one. To obtain a seed arrangement that can be used for an actual treatment, one may choose an intuitive approach based on trial and error. Tools to assist the treatment planning have been developed to compute the dosimetry with homedesigned codes or to convert the data to a format suitable to commercial radiation treatment planning systems (8). Although good distributions can be obtained, the repetition of the procedure can require several hours to complete and is highly dependent on the skill of the individual planner. Reduction in treatment planning time and maintenance of procedure consistency can only be achieved if the planning process is automated. It has been reported that the use of a leastsquare optimization procedure (rninimization of the squared difference between the desired and calculated values) for permanent prostate implants can reduce planning time by a factor of 10 compared to the manual approach (24). This represents a substantial improvement and clearly indicates the direction to be followed. However, optimization by means of a local search algorithm cannot guarantee that the identified source configuration will be the best configuration. To illustrate this, we have performed a cost function minimization of an implant. The costfunction is a quantity used to define the quality of the implant and will be explained in detail later in the text. Figure 1 presents the values of these cost functions at each iteration during the minimization process. The two curves represent calculations that differ only by their initial seed distribution. In both cases, the algorithm optimized the costfunction by attempting small changes to the seed arrangement. In at least one case, the calculation stopped at a suboptimal solution because no distribution with a lower cost function could be found nearby the current one. Clearly, the algorithm was trapped in a local minimum. The average harmonic doses on the prostatic contour (l/average [Z doseP’], where the sum is taken over all dose points) (34), were found to differ by 8 Gy. This is clinically significant and indicates that one solution is better than the other. In cases where a cost function has multiple minima, an optimization technique allowing for wider sampling and hill climbing for escaping from local minima is required. Simulated annealing (SA), first introduced by Kirkpatrick et al. (13), and other related algorithms (Fast SA, Very Fast SA, or others) are optimization techniques that can process cost functions with arbitrary boundary conditions with the statistical guarantee of finding an optimal solution. Equally important, SA class algorithms are easily coded compared to other nonlinear optimization algorithms. A multitude of problems arising from very different fields have been solved using SA algorithms (11). In radiation therapy, this technique can be used to define the parameters of an external beam treatment (15, 33). Recently,
Volume
36, Number
3, 1996
a.7
7.2
l
i
a7. 0
10000 20000 30000 40000 50000 Iteration Number
Fig. 1. Result of costfunction least square optimization: the two curves differ only by their initial random seed distribution. One of the curves becomes trapped in a local minimum. Units are arbitrary.
this approach has been successfully
adapted to brachy
therapy to optimize the dose distribution of up to ten seeds in a cylindrical vaginal vault applicator (20 possible positions) (26).
In this article, we present the optimization of transrectal ultrasound guided permanent ‘251 prostatic implants with the aid of a fast simulated annealing algorithm. After a brief description of the fast SA algorithm and the procedure, the link between clinical requirements and the cost function is described. An example is used througout the article to illustrate implementation of the method.
METHODS
AND
MATERIALS
Fast simulated annealing algorithm
The simulated annealing algorithm (SA) has been described on several occasions. Only the aspects essential to the comprehension of the procedure are presented here. The algorithm is rooted in an analogy between minimization problems and the cooling of a liquid. When a liquid is cooled slowly,
its internal energy is lowered
and the
liquid eventually freezes with its atoms placed in a highly structured lattice in a state of minimum energy. Even with a huge number of atoms, nature finds a way to do it, provided cooling is sufficiently slow (annealing). However,
Optimization of permanent 9 prostate implants 0 J. PouLIoTet al.
if the cooling is too fast (quenching), the system will end up in an unstable amorphous or polycrystalline state with higher internal energy. In the present case, the energy is replaced by a cost function describing the adequacy of the implant. Following the same analogy, the liquid temperature is replaced by a parameter that plays the same role. The algorithm slowly cools the system to find a minimal cost function solution. A detailed description can be found in Aarts and Korst (1). According to Press et aZ.(20), the four major components of an SA algorithm are: (a) a concise description of the problem (how to generate any possible seed arrangement); (b) a transition mechanism (given a seed arrangement, how to change it to obtain a new one); (c) an objective function to be minimized [a cost function E(k) playing the role of internal energy]; and (d) a cooling schedule for the temperature parameter T(k) at annealing step k. The algorithm runs as follows: beginning with a seed arrangement at step k (the current one), a transition is generated (step k+l). The new seed arrangement (K+l) is obtained by moving at least one seed at random from an occupied position to an unoccupied one. Then, the algorithm must decide if the new seed distribution is accepted to become the current one and, thus, replace the former one. The criterion used is based on the costfunction values for the two arrangements: E(k + 1) and E(k). The probabity P(AE) of accepting the transition depends on the difference AE = E (k + 1)  E(k). For a transition resulting in a lower value of the cost function (AE < 0), the new distribution is always accepted as the current one. A change that results in an increase of the cost function, i.e., AE > 0, is accepted with a probability given by P(AE) = exp[AE/T(k)]. Thus, for a given positive AE, transitions are accepted more easily at higher temperatures. The temperature parameter T(k) is reduced at each step K until the system no longer evolves. The difference between fast simulated annealing and simulated annealing is discussed below in the cooling schedule section. In the case where T would be set to zero for all k, the probability of accepting a transition with a positive AE would be zero and the optimization would behave as a regular local search algorithm. Our algorithm uses the following steps: 1. Choice of a set of needle locations within the planned target volume (PTV). 2. Choice of an initial number of lz51 seeds of given activity. 3. Generation of an initial random seed distribution; k set to 1. 4. Initialization of the initial temperature parameter TO. 5. Transition of the current seed distribution by moving one (or more) seeds; k incremented to k+ 1 6. Comparison of the cost functions between the seed distributions k and k+ 1. (a) If E(k+ 1) < E(k), then the distribution k+l is accepted and becomes the current
713
seed distribution. (b) If E(k+l) 2 E(k), then the distribution k+l is accepted with probability P(AE). Upon rejection, distribution k remains the current distribution. 7. Reduction of T(k) according to cooling schedule; T(k) = TJ(k+l). 8. Return to step 5 until optimization is satisfactory. Data input Optimization of permanent prostate implants using SA algorithm requires 3D information regarding prostatic shape, extent of local disease and relative position of all doselimiting normal tissue structures. A set of transverse images are acquired using a transrectal ultrasound with the patient in lithotomy position. A computer interfaced to a frame grabber is connected directly to one of the video outputs of the echographic imaging system. The first image (plane # 0) is taken at the base of the prostate and sequential images are taken at 5 mm intervals, until the prostatic apex is reached. This usually yields between seven to nine images. Prostatic contouring is obtained from each image and digitized to produce the xyz coordinates (Fig. 2a). This represents the planning target volume (PTV). The figure uses continuous lines for illustration purpose only. Prostate contours define the target volume to be treated and allows for prostatic volume calculation. A maximum number of seeds of given activity is specified (see Eq. 5) and subsequently distributed along the selected needles. The number of needles selected and needle distribution within the PTV is chosen such that an adequate prostatic volume coverage is achieved. Clinical experience often dictates the needle number and configuration within the PTV. Typically 12 to 18 needles are required. Needle positions must coincide with predefined positions on a universal template used to guide needles during the implantation. Holes on the template allow for insertion of the needles at 5 mm increments horizontally and vertically. However, a 1 cm gap between needle positions is generally preferred. The program optimizes the number of seeds and their positions, but to limit the number of variable parameters, the needle locations are specified. Thus, the prostate contours, the needle locations, and the maximum number of seeds represents the basic set of data needed by the program. Seed conjiguration space The configuration space is defined as the set of all permitted seed positions. These potential seed positions are distributed along each needle at 5 mm intervals from the base (plane # 0) to the apex. Seed positions coincide with contour planes (0, 1, 2, . . .) and are located either inside, or outside of the prostate contours. In this latter case, a margin must be specified as a parameter. Although the margin can take any value, we generally use 5 mm or less.
714
I. J. Radiation Oncology l Biology l Physics
Volume 36, Number 3, 1996
Permitted seed positions are represented by the open circles. In the configuration space, some positions are occupied by a seed, while others remain unoccupied. The set of positions occupied by seeds defines a seed arrangement (or distribution). Figure 2c illustrates a seed distribution. Occupied seed positions are marked with black dots. Given an actual distribution, one can define a neighboring distribution by moving a seed from an occupied position to an unoccupied one. All arrangements obtained in this manner define the immediate neighborhood of that distribution. One can also expand the neighborhood horizon by moving more than one seed. Cost function
The cost function E(k) is used to compare a given seed distribution k relative to the next k+ 1. This cost function must be constrained by the physician’s prescription to reflect the clinical requirements. The clinical criteria used for the optimization are: (a) the delivery of a minimum prescribed peripheral dose (mPD) on prostate contour points; (b) a uniform dose distribution within the planned target volume, i.e., no hot spot or cold spot; and (c) a limited dose to the urethra and the rectum. The degree of fulfillment of the first clinical criterion is measured by the use of the dosecontour potential (DCP).
where Di is the dose contribution from all seeds at a given dosecontour point i (DCpoints) and D, is an adjustable value to take into account the prescribed dose. The exponents p and a are chosen empirically; p must take an
(A)
(B)
DCP Fig. 2. (a) 3Dobliqueview of the prostatecontours.(b) Needles areaddedto (a) with the permittedseedpositions(opencircles) and the uniformity potential (UP) points (crosssymbols).(c) Optimizedseeddistributionfor needlelocationsgiven in (b) and parametersof Table 1. Black circles indicate the presenceand positionof a seed.Prostateapex is on the lefthand side.
UP 30
25
25
20
‘: 3 9
15
m
10 5
This implies that seeds can be placed outside the prostate only in the regions of the base and the prostatic apex for we do not allow needles to be inserted outside the middle contours (in that case, no contour would be intercepted by a needle). A typical seed configuration space provides between 150 to 200 potential seed positions. Figure 2b illustrates the needle positions relative to prostate contours.
0 I
L 200
Dose
(Gy)
Dose
400 (Gy)
Fig. 3. (a) Dose contour potential fH (DCP). (b) Uniformity potential gB (UP). The parametersusedfor the potentialsare presentedin Table 1.
Optimization of permanent ‘*‘I prostate implants 0 J. POULIOT et al.
even value. The greater the p value, the steeper the potential wall. The exponent a yields an assymetric shape to the well: doses below D, are more penalized than the doses higher. Equation 1 is evaluated on each contour point. An example of a typical potential is shown on Fig. 3a using p and a values of Table 1. The second clinical criterion requires a uniform dose distribution throughout the target volume. Cold spots and hot spots must be avoided; the latter could be created if two or more seeds occupy consecutive positions on a needle. This criterion can be satisfied by requiring a rather uniform dose on points j located at middistance between any two seed positions on a needle. These points are shown in Fig. 2b (x symbol). We define a uniformity potential (UP) calculated on points j:
DU () Bj=‘
P
R
Q
/D2\a
(W. 2)
Points j are called uniformity check points. With a proper choice of parameters D,, p, and a, Eq. 2 gives a high Bj value at uniformity check pointj either if Dj is too low (cold spot) or too high (tight needle loading). The parameters a and p play the same role as in the DCP, although they can take different values in DCP and UP. The value of Bj is evaluated at each midpoint between permitted seed positions. Examples of typical potential are shown in Fig. 3 using the parameter values of Table 1. The cost function E(k) of the seed distribution k is given by the sum of the two terms representing the criteria (a) and (b), respectively. The third clinical criterion (c) is discussed in the Discussion section. The costfunction value may be expressed as: (Eq. 3) where N and M represent the number of dose points on the prostate contours (DCpoints) and the number of UPpoints, respectively. In the case where one were to change the p and a values in DCP and/or UP, we define f and g normalization factors to make the relative normalization between the two terms independent of the shape of the potentials (a and p values chosen). The normalization factors are:
f
= (0,
6P
2
+
+
4”
(S)P (DC

;
715
Table 1. Cost function parameters DoseContour Potential (DCP) 210 Gy DC a 4 4 P s 40 Gy Uniformity Potential (UP) 170 Gy D” a 4 2 F3 40 Gy Weighting factor W 0.02
With the use off and g, the average value of the product f X H,,i evaluated at doses DCS and Dc+S is the same as the average of g X Bj evaluated at D,6 and D,+S. The importance of a clinical criterion over the other is adjusted with the weighting factor W. Increasing W favors dose uniformity within the prostate gland over the dose adequacy on its contour. Variable seednumber
The search for the best dose distribution requires that both the number and seed positions be allowed to change during the optimization process. Although it is possible to vary both parameters during the annealing process, this requires the introduction of a weight factor in the decision to either change the number of seeds or to modify their position (26). The approach followed in this work is somewhat different. The total number of seeds is fixed at an initial value. However, during the definition of the configuration space, eight new possible positions are added. These positions are located at an infinite distance from the implant such that the presence of a seed in one of these positions has no effect on the dose distribution (in practice, a finite distance of 10 m is sufficient). No uniformity point is defined between these distant positions. During the annealing process, if one of these far positions is occupied by a seed, this effectively reduces the number of seeds inside the implant. The program can add or remove seeds from these far positions with the same probability than any other possible positions. For instance, if the optimization process begins with 70 seeds, the optimized distributions could have between 62 and 70 seeds. The initial total number of seeds (Ni”ir) is estimated from the prostatic volume V (in cubic centimeters) and the activity per seed A (in mCi), using the empirical relation: Ninit = 4 + [(4.674 X V”.562)/A]
(Eq. 3
sy
This equation represents a fit over the data accumulated from a set of patients, with an additional four seeds added that can be used if needed. If one chooses to increase D, (Table l), for instance, to change the margin around the prostatic contour, a greater total activity will be needed.
716
I. J. RadiationOncology 0 Biology l Physics
80.0
r
Volume 36, Number 3, 1996 The same value is also used for the probability of three seedsbeing displaced as compared to two. The cooling sequenceis stopped when a large number of iterations have been performed without any change in the costfunction value.
RESULTS
5000 10000 15000 Iteration Number
20000
Fig. 4. Reductionof costfunctionvalue with fast simulatedannealing.The cost function evolvesfrom a high value to its minimumvalue. No changein cost function is observedafter 12000 iterations.Units are arbitrary.
An example of implant optimization by minimizing the costfunction value with iteration number is shown on Fig. 4. Note that asthe iterations proceed, someconfigurations involving an increaseof the cost function (decreaseof the implant quality) are accepted, especially at the beginning (high “temperature” region): this is the essenceof simulated annealing. The final or optimized seeddistribution is shown on Fig. 2c. The initial temperature TO was set to 100 and reduced according to the FSA algorithm (Eq. 6). This value has been found adequate in comparison to the costfunction values obtained with random initial seeddistributions and ranging between 50 and 150. Typically, 20,000 iterations, a number sufficient to minimize the cost function, can be computed with the code written in C language in approximately 15 min on a work station.’
Implant Cooling schedule In the annealing schedule typical of standard SA, the temperature is decreasedonly after a large number of distributions have been tested, assuringthat thermal equilibrium has been reached, at least statistically. This approach gives consistent results in the search of the global minimum. However, as stated in the literature, standard SA tends to be too slow (11). The fast simulated annealing (FSA) used here differs from SA in two respectsand can be used as an alternative to standard SA: (a) The temperature parameter T(k) used to calculate P(AE) is reduced at every iteration according to the relation
T(k) = T,l(k + 1)
(Eq. 6)
where k is the iteration number and T,, the initial temperature; (b) transition to a new configuration can be obtained by moving more than one seed before Eq. 3 is used to evaluate the cost function for the new seed distribution. In the present case, this implies the possibility to probe the neighborhood of a present distribution which differs by more than one seedposition. The cooling schedule of Eq. 6 is admissible only if a small number of long jumps are allowed (28). Therefore, the relative probability of moving two seedsas opposedto one has been set at 25%. ’ Spare5SUN workstation,Mountain View, CA.
evaluation
In order to evaluate the quality of the preimplant dosimetry, several outputs are generated by the program. Figure 5 illustrates the needle loading obtained from the optimum seed distribution. The isodose distributions are then calculated and displayed simultaneously with the prostate contours and the seedpositions for each transverse plane. In Fig. 6, it can be seen that the potential parameters of the cost function presented in Table 1 would typically result in an average harmonic dose of 210 Gy with a minimum peripheral dose close to 160 Gy. Dosimetric data underlying dose calculations are from TG45 (23). In addition to the isodosedistributions, Chang et aL(6) suggestedthat the graphical displays of dose area could be an efficient tool to assistin the evaluation of the implant. We have, therefore, included area diagrams. Figure 7 illustrates the area covered on each plane by the prostate (black squares), and the areas covered by the 160 Gy (open squares)and 80 Gy (gray diamonds) doses,respectively. This graph provides a quick way to verify that proper margins were set around each contour. Dosevolume histograms (17), both cumulative and differential, are also made available by the program. Matched peripheral dose, can also be computed. To further characterize the preplan dosimetry and help evaluate the implant quality, the program also reports minimum dose, average dose, and harmonic dose at the prostate surface for each plane.
Optimization of permanent “‘1 prostate implants 0 J. POULIOT et al.
#:
1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13:
x
c d c E b c d e b c D E F
Y
3.0 3.0 2.5 2.5 2.0 2.0 2.0 2.0 1.0 1.0 1.0 1.0 1.0
:I0
9
8
7
6
5
4
717
3
2
10
. . . . .. . . . . .. . . .. . . . . .. .. :
Fig. 5. Needle loading for the optimum seed distribution obtained from the source space configuration X and Y coordinates correspond to template hole positions.
DISCUSSION To demonstrate the relative contributions of the two potentials DCP and UP to the cost function, we have performed an optimization with the weighting factor W set to zero. This has the effect of turning off the UP potential and the optimization takes only into account the doses on the prostate contours. Figure 8 shows the needle loading obtained with the same space configuration as before (Fig. 2b) but with W = 0. In this case, the algorithm ended with 49 seeds to produce the optimum solution, compared to 47 when both DCP and UP were activated in the cost
AaRbCcDdEeFfG
of Fig. 2b.
function (W = 0.02). More importantly, one can observe on the needle loading presented in Fig. 8 that without the UP potential, the seeds tend to agglomerate themselves at the center of the implant. Although the minimum peripheral dose condition is still fulfilled (Fig. 9), this situation creates very high doses within the implant and is clinically unacceptable. In contrast, it would be easy to forbid the
*Orno m 15.0
5.0
0.0 Fig. 6. Isodose distribution of Plane #3 obtained from the optimum seed distribution. The reference isodose (160 Gy) is in black and others are in gray. Diamonds illustrate the intersection of a seed with the plane. Results from other planes are similar.
I 0
H Prostate Ml6OGY e*60Gy
\
2 4 6 810 Transverse Plane Number Fig. 7. Dose area histograms.
718
I. J. Radiation Oncology l Biology l Physics
#: I: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13:
x
Y 3.0 3.0 2.5 2.5 2.0 2.0 2.0 2.0 1.0 1.0 1.0 1.0 1.0
Fi C E b ii e b C D E F
:I0 .. ..
I
Volume 36, Number 3, 1996
9
8
7
6
5
4
3
2
10
I
I
I
I
I
I
I
I
I
I
. . .. . . . . .. .. : : .. . . :
Fig. 8. Needle loading for an optimum seed distribution obtained with the weight factor W set to zero. This implies that the uniformity potential is not taken into account in the cost function.
occupation of two consecutive seed positions by assigning an arbitrary large value to IV. The average prostatic volume implanted was 30 cc (range: 2045 cc) with an average total implanted activity of 31.3 mCi (range: 2439 mCi). The seed activity is chosen in the range of 0.49 to 0.67 mCi per seed. This results in an average of 58 sources implanted (range: 4266). This number of seeds and activity per seed result in adequate dose distributions despite seed migration, while minimizing total number of seeds and needles implanted. We believe this provides a low costtobenefit ratio. Within this activity range, parameters shown in Table 1 yield adequate results, and are the
AaBbCcDdEeFfG
Fig. 9. Isodose distribution of Plane #3 for an optimum seed distribution obtained with the weight factor W set to zero (see caption of Fig. 8).
parameters used at our institution. However, the cost function is independent of the seed activity and the algorithm can run with any activity. With the present version of the program, the dose to the rectum is limited by assuring that the lowest row of needles implanted is not located too close to the rectal wall (27). Also, the central needles are placed so as to avoid the region around the urethra. In fact, to minimize the risk of longterm morbidity, efforts are made to maintain the rectal surface dose below 100 Gy and the central urethral dose below 400 Gy (32). The latter limit is also prevented by the presence of the uniformity potential that penalizes excessive doses along needles. With the constraints imposed on needle locations, the third clinical criterion is satisfied. However, the cost function could easily be modified to include specific constraints to the urethral and rectal doses. On certain occasions, the treatment may require an additional optimization step. Initially, a set of needle locations is chosen, taking into account the shape of the prostate contour, and the position of the urethra and rectum. Then, the maximum number of seeds and the seed activity are chosen. Seed positions are optimized using the FSA algorithm. Once completed, the needle loading is evaluated. If a needle is too closely packed with seeds, one may consider adding an extra needle nearby and restarting the optimization process. In contrast, if a needle is too sparsely loaded (for instance, only one seed on the needle), the needle may be removed and optimization restarted. Needles with only one seed are rarely used thereby minimizing the total number of needles implanted and the consequent trauma of implantation. The existence of local minima with significantly higher costfunction values depends on the exact expression of the
Optimization of permanent lz51 prostate implants 0 J. POULIOT et al.
cost function used (Eq. 3) and the way it is defined (contour and uniformity check points, number of distant seed positions, and so on). In our study, the FSA algorithm was found to result in a more optimal solution than the downhill search algorithm. When proper parameters are selected, the algorithm yields a very good solution within a reasonable time frame. Although application of the algorithm, as detailed in this article, does not allow one to certify that the best solution is always obtained, it does, however, efficiently avoid suboptimal solutions with higher costfunction values (see Fig. 1). In fact, near the neighborhood of the absolute minimal costfunction solution, there exists a few tens of local minimal solutions with nearly identical costfunction values. They differ typically by only a few source positions near the center of the implant and are clinically equivalent. To get the very best one with good confidence, one would have to begin the search with a value of To considerably larger than the one used here and iterate hundreds of thousands of times. This would take hours to compute, with little or no clinical gain. CONCLUSION Fast simulated annealing has been easily adapted to a broad range of problems. In this study, the development of
719
an algorithm based on FSA has proven to be an efficient tool to generate the preimplant dosimetry by optimizing the seed distribution for permanent lz51prostatic implants. Not only will this approach find the best distribution, but the planning is obtained much faster, sometimes in a few minutes compared to other semiautomatic or manual planning evaluations which may require hours to days to complete. The parameters of the cost function can be modified to take into account clinical criteria specific to an institution. For our set of parameters used in the cost function (Table l), the algorithm gives good dose distributions in several minutes. One of the major advantages of FSA is its great adaptability. The cost function can be modified to account for different constraint specifications. For instance, we are currently working on a modified version of the code starting the computation with a large number of initial needle positions to optimize the implant with a cost function that gives a penalty according to the total number of needles and to their loading scheme. In this way, the optimum needle arrangement could be achieved automatically. Another example would be to introduce in the cost function the contribution from rectum and/or urethra dose points thereby minimizing overdosage to these critical normal structures.
REFERENCES 1. Aarts, E.; Korst, J. Simulated annealing and Boltzmann 2.
3. 4.
5.
6. 7.
8. 9.
10.
machines. Chichester: Wiley; 1989. Arterbery, V. E.; Wallner, K.; Roy, J.; Fuks, Z. Shortterm morbidity from CTplanned transperineal 1125 prostate implants. Int. J. Radiat. Oncol. Biol. Phys. 25:661667; 1993. Biasko, J. C.; Grimm, P. D.; Radge, H. Brachytherapy and organ preservation in the management of carcinoma of the prostate. Semin. Radiat. Oncol. 3:240249; 1993. Blasko, J. C.; Radge, H.; Grikk, P. D. Transperineal ultrasoundguided implantation of the prostate: Morbidity and complications. Stand. J. Urol. Nephrol. Suppl. 137: 113l 18; 1991. Blasko, J. C.; Radge, H.; Schumacher, D. Transperineal percutaneous iodine125 implantation for prostatic carcinoma using transrectal ultrasound and template guidance. Endocuriether. Hypertherm. Oncol. 3:131139; 1987. Chang, F.; Chang, P.; Winton, L.; Share, F. A dosearea study of ultrasoundguided iodine125 prostate seed implants. Endocurie. Hypertherm. Oncol. 8:119125; 1992. DeBlasio, D.; Hilaris, B.; Nori, D.; Fuks, Z.; Whitmore, W.; Fair, W.; Anderson, L. Permanent interstitial implantation of prostatic cancer in the 1980s. Endocuriether. Hypertherm. Oncol. 4:193201; 1988. Feygelman, V.; Noriega, B. K.; Sanders, R. M.; Friedland, J. L. A spreadsheet technique for dosimetry of transperineal prostate implants. Med. Phys. 22:97100; 1995. Holm, H. H.; Pedersen, J. F.; Hansen, H.; Stroyer, I. Transperineal I125 iodine seed implantation in prostatic cancer guided by transrectal ultrasonoghraphy. J. Urol. 130:283286; 1983. Holm, H. H.; TorpPederson, S.; Myschetsky, P. Transperineal seedimplantation guided by biplanar transrectal ultrasound. Urology 36:248252; 1990.
11. Ingber, L. Simulated annealing: Practice vs. theory. J. Math. Comput. Modelling 18(11):2957; 1993. 12. Kaye, K. W.; Olson, D. J.; Lightner, D. J.; Payne, J. T. Improved technique for prostate seed implantation: Combined ultrasound and fluoroscopic guidance. J. Endourol. 6:6166; 1992. 13. Kirkpatrick, S.; Gelatt, C. D., Jr; Vecchi, M. P. Optimization by simulated annealing. Science 220(4598):67 l680; 1983. 14. Kuban, D. A.; ElMahdi, A. M.; Schellhammer, P. F. I125 interstitial implantation for prostate cancer: What have we learned 10 years later? Cancer 63:24152420; 1989. 15. Morril, S. M.; Lam, K. S.; Lane, R. G.; Langer, M.; Rosen, I. I. Very fast simulated reannealing in radiation therapy treatment plan optimization. Int. J. Radiat. Oncol. Biol. Phys. 31(1):179188; 1995. 16. Morton, J. D.; Peschel, R. E. Iodine125 implants vs. external beam therapy for stages A2, B, and C prostate cancer. Int. J. Radiat. Oncol. Biol. Phys. 14: 1153l 157; 1988. 17. Niemierko, A.; Go&en, M. Dosevolume distributions: A new approach to dosevolume histograms in threedimensional treatment planning. Med. Phys. 21(1):3l 1; 1994. 18. Osian, A. D.; Anderson, L. L.; Linares, L.; Nori, D.; Hiltis, B. Treatment planning for permanent and temporary percutaneous implants with custom made templates. Int. J. Radiat. Oncol. Biol. Phys. 16:219223; 1989. 19. Osian, A. D.; Nori, D. Conformal brachytherapy for carcinoma of the prostate. Endocurie. Hypertherm. Oncol. 10: 1524; 1994. 20. Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; FIannery, B. P.Numerical recipies, 2nd ed. Cambridge, England: Cambridge Univ. Press; 1992 21. Priestly, J. B., Jr; Beyer, D. C. Guided brachytherapy for treatment of confined prostate cancer. Urology 40:2732; 1992.
720
I. J. Radiation Oncology l Biology l Physics
22. Ragde, H.; Blasko, J. C.; Schumacher, D. Use of transrectal ultrasound in transperineal Iodine125 seeding for prostate cancer: Methodology. J. Endourol. 3:209218; 1989. 23. Ravinder, N. ; et al. Dosimetry of interstitial brachytherapy sources: Recommendations of the AAPM Radiation Therapy Committee Task Group No. 43. Med. Phys. 22(2):209233; 1995. 24. Roy, J. N.; Wallner, K.; ChiuTsao, S.; Anderson, L. L.; Ling, C. C. Ctbased optimized planning for transperineal prostate implant with customized template. Int. J. Radiat. Oncol. Biol. Phys. 21:483489; 1991. 25. Roy, J. N.; Wallner, K.; Harrington, P. J.; Ling, C. C.; Anderson, L. L. A CTbased evaluation method for permanent implants, application to prostate. Int. J. Radiat. Oncol. Biol. Phys. 26:163169; 1993. 26. Sloboda, R. S. Optimization of brachytherapy dose distributions by simulated annealing. Med. Phys. 19(4):955964; 1992. 27. Stock, R. G.; Stone, N. N.; Wesson, M. F.; DeWyngaert, J. K. A modified technique allowing interactive ultrasoundguided threedimensionnal transperineal prostate implantation. Int. J. Radiat. Oncol. Biol. Phys. 32(1):219225; 1995.
Volume 36, Number 3, 1996 28. Szu, H.; Hartley, R. Fast simulated annealing. Phys. Lett. A122(3,4):157162; 1987. 29. Wallner, K. Brachtherapy for prostate cancer. ASTRO, San Francisco, California, October 4, 1994. 30. Wallner, K. Iodine 125 brachytherapy for early stage prostate cancer: New techniques may achieve better results. Oncology 5:115122; 1991. 3 1. Wallner, K.; ChiuTsao, S. T.; Roy, J. An improved method for CTplanned transperineal iodine125 prostate implants. J. Urol. 146:9095; 1991. 32. Wallner, K.; Roy, J.; Zelefsky, M.; Fuks, Z.; Harrison, L. Dosimetry guidelines to minimize urethral and rectal morbidity following transperineal 125I prostate brachytherapy. Int. J. Radiat. Oncol. Biol. Phys. 32:465471; 1995. 33. Webb, S. Optimization by simulated annealing of threedimensional conformal treatment planning for radiation fields defined by a multileafcollimator. Phys. Med. Biol. 36: 12011226; 1991. 34. Yan Yu, F. M.; Waterman, F. M.; Suntharalingam, N. Quantitative assessment and dose specification of iodine 125 permanent implants. Med. Phys. 21(6):959; 1992.