The effect of statistical uncertainty on inverse treatment planning ...

Home

Search

Collections

Journals

About

Contact us

My IOPscience

The effect of statistical uncertainty on inverse treatment planning based on Monte Carlo dose calculation

This article has been downloaded from IOPscience. Please scroll down to see the full text article. 2000 Phys. Med. Biol. 45 3601 (http://iopscience.iop.org/0031-9155/45/12/307) View the table of contents for this issue, or go to the journal homepage for more

Download details: IP Address: 129.78.32.23 The article was downloaded on 10/05/2012 at 16:17

Please note that terms and conditions apply.

Phys. Med. Biol. 45 (2000) 3601–3613. Printed in the UK

PII: S0031-9155(00)14652-4

The effect of statistical uncertainty on inverse treatment planning based on Monte Carlo dose calculation Robert Jeraj†‡! and Paul Keall§ † Joˇzef Stefan Institute, Ljubljana, Slovenia ‡ Department of Medical Physics, University of Wisconsin–Madison, Madison, WI, USA § Department of Radiation Oncology, Medical College of Virginia Hospitals, Virginia Commonwealth University, Richmond, VA, USA E-mail: [email protected] Received 12 June 2000 Abstract. The effect of the statistical uncertainty, or noise, in inverse treatment planning for intensity modulated radiotherapy (IMRT) based on Monte Carlo dose calculation was studied. Sets of Monte Carlo beamlets were calculated to give uncertainties at Dmax ranging from 0.2% to 4% for a lung tumour plan. The weights of these beamlets were optimized using a previously described procedure based on a simulated annealing optimization algorithm. Several different objective functions were used. It was determined that the use of Monte Carlo dose calculation in inverse treatment planning introduces two errors in the calculated plan. In addition to the statistical error due to the statistical uncertainty of the Monte Carlo calculation, a noise convergence error also appears. For the statistical error it was determined that apparently successfully optimized plans with a noisy dose calculation (3% 1σ at Dmax ), which satisfied the required uniformity of the dose within the tumour, showed as much as 7% underdose when recalculated with a noise-free dose calculation. The statistical error is larger towards the tumour and is only weakly dependent on the choice of objective function. The noise convergence error appears because the optimum weights are determined using a noisy calculation, which is different from the optimum weights determined for a noise-free calculation. Unlike the statistical error, the noise convergence error is generally larger outside the tumour, is case dependent and strongly depends on the required objectives. (Some figures in this article are in colour only in the electronic version; see www.iop.org)

1. Introduction When performing Monte Carlo calculations, the results are always burdened with some statistical uncertainty. From the central limit theorem it is known that the true value of the calculated quantity is distributed with the Gaussian probability around the estimated mean value with the standard deviation (statistical error) σ . It should be noted that there is also some systematic error, but this is neglected in our consideration. From the statistical analysis, it is also known that the statistical error of Monte Carlo calculations√is proportional to the inverse square root of the number of simulated particle histories, 1/ N . For radiotherapy planning, the calculated results are given in terms of absorbed doses in different parts of the patient geometry. It has recently been theorized (Sempau and Bielajew 2000, Deasy 2000, Jeraj 1999) and shown (Keall et al 2000) that the absolute statistical error of Monte Carlo dose ! Address for correspondence: Department of Medical Physics, University of Wisconsin–Madison, 1530 Medical Sciences Center, 1300 University Ave, Madison, WI 53706, USA. 0031-9155/00/123601+13$30.00

© 2000 IOP Publishing Ltd

3601

3602

R Jeraj and P Keall

calculation, √ σd , is approximately proportional to the square root of the dose in the scoring region, d. Combining both dependences, the following relations can be obtained ! d σd = C (1) N where d is the dose in the scoring region and N the number of simulated histories in the problem. The constant C was estimated by Sempau and Bielajew (2000) to be C = (Sr /ρ)η/$V 2/3 , where Sr /ρ is the restricted mass stopping power, η a numerical parameter dependent on the mean path length of the electrons in a scoring region (approximately 1) and $V is the volume of the scoring region. It should be stressed that this simple relation should be used with caution in Monte Carlo based intensity modulated radiation therapy (IMRT), where more beams with different weights and number of simulated histories are combined, because it is valid for one beam only. When more beams of different weight, which may also be calculated to different precision, are combined, N should be seen as a weighted sum of the number of particles simulated in each beam contributing to the dose in the scoring region, as will be shown later. When Monte Carlo transport is used to calculate dose distributions in conventional treatment planning, the statistical error arising from the statistical nature of the Monte Carlo simulation is present (Keall et al 2000, Buffa et al 1999). But when using Monte Carlo transport to calculate dose distributions for inverse treatment planning, one has to deal with an additional error (Jeraj and Keall 1999b). Along with the statistical error, which is the same as that in conventional treatment planning, a so-called noise convergence error also occurs. The noise convergence error should arise because optimization converges to the ‘optimal’ solution for the noisy dose calculations, which is different from the optimal solution for the noise-free dose calculations. It is worth mentioning that the convergence error should be present in any inverse treatment plan which uses either inaccurate dose calculation (e.g. pencil beam) or imprecise dose calculation (e.g. Monte Carlo). The main aim of this paper is to demonstrate the existence of statistical and noise convergence errors in inverse treatment planning based on Monte Carlo dose calculation and to investigate characteristics of these errors. 2. Statistical and noise convergence errors Before we continue, some terminology needs to be introduced. A solution of inverse treatment planning (intensity distribution) may be denoted as wN . This is a vector consisting of weights (or importances) of each beamlet. A beamlet is the dose distribution resulting from one beam element, or pencil beam, normalized arbitrarily. The index N should remind us that a particular solution of the inverse treatment planning is dependent on the number of simulated histories in the (noisy) beamlet B N . wN is calculated by minimizing a given objective function, OF: wN = min(OF(B N )). w∈P

(2)

P represents the space of possible solutions. The noise-free solution, which is obtained with the noise-free beamlets B ∞ , shall be denoted as w∞ and similarly in the solution for noisy beamlets it is calculated by minimizing a given objective function: w∞ = min(OF(B ∞ )). w∈P

(3)

One should be aware that there are no real noise-free solutions in any Monte Carlo calculation. They can be obtained only when N → ∞. However, in practice, when N is very large and statistical uncertainty almost negligible, the solutions can be considered as almost noise-free and so were in our study too.

3603

Errors in Monte Carlo based inverse treatment planning

The dose distribution vector is a product of the weight distribution and beamlet dose distributions. The ‘optimal’ dose distribution for noisy beamlets, dNN , or noisy dose distribution is thus obtained by summing up the noisy beamlets according to the calculated solution of inverse treatment planning, so that it can be written dNN = wN B N .

(4)

∞ d∞ = w∞ B ∞ .

(5)

dN∞ = wN B ∞ .

(6)

$stat = dNN − dN∞ = wN (B N − B ∞ ).

(7)

Throughout this paper the superscript refers to the number of particles in each beamlet calculation, and the subscript refers to the number of particles in each beamlet used to obtain ∞ the solution of the plan. The noise-free dose distribution, d∞ , obtained with the noise-free beamlets, can be similarly expressed as a product of the noise-free solution and noise-free beamlets: When calculating the errors of Monte Carlo based inverse treatment planning another dose distribution will be needed. This is the so-called recalculated dose distribution. The recalculated dose distribution is obtained when the noisy solution is recalculated with the noise-free beamlets, which can be written as The statistical error in the inverse treatment plan appears because of the statistical nature of the Monte Carlo calculation used to calculate dose distributions of individual beamlets. For the study of the statistical error of the inverse treatment plan one needs to know the difference between the noisy and recalculated dose distribution. Hence the statistical error of the plan, $stat , can be determined by subtracting the recalculated dose from the noisy dose distribution: It should be stressed that the statistical error of the inverse treatment plan is very similar to the statistical error of the conventional treatment plan and appears because of the statistical uncertainty (or noise) of the Monte Carlo dose calculation. However, because the inverse treatment plan typically consists of many fields some care should be taken in evaluating the statistical error of the whole plan. If the number of particles contributing to the dose score is large enough, so that the central limit theorem is valid, it can be assumed that the noisy beamlets B N can be obtained by adding to the individual voxels of the noise-free beamlets B ∞ some value sampled from the Gaussian distribution G(N ): B N = B ∞ + G(N ).

(8)

√

According to equation (1) the standard deviation of the Gaussian is Cd/N, where d is the dose in a voxel and N the number of histories simulated in a particular beamlet. Equation (7) can thus be rewritten as (9)

$stat = wN G(N ).

This shows that the statistical error is reflected in the inverse treatment plan as a random statistical noise, dependent on the dose level. More than the differences themselves, the standard deviations of the distributions $stat are interesting. Because beamlets are not correlated, the standard deviations of the differences, σstat , can be calculated via 2 σstat

=

'$2stat (

2

= '(wG(N )) ( =

M " j =1

wN2 j σj2

=C

M w2 " Nj j =1

Nj

B Nj .

(10)

M is the number of beamlets in the problem and Nj the number of simulated histories in the j th beamlet. In general, σstat is case dependent because it depends on a particular solution of the

3604

R Jeraj and P Keall

inverse treatment plan w and the statistical error does not necessarily decrease with the inverse square root of the number of simulated histories as one might expect. In our study, all beamlets were calculated to the same precision, so that Nj = Ntot /M, where Ntot is the total number # 2 of histories simulated in the problem. Since the sum M j =1 wj dj is approximately constant (deviations are due to a combined effect of statistical and noise convergence error, which is forced to be zero in the tumour if the objective √ function requires uniform dose distribution), the statistical error is dependent only on 1/ Ntot , which is an expected result. However, if the number of simulated histories is not chosen appropriately, for example when low weight beamlets are calculated to high precision √ and vice versa, it may happen that the statistical error of the plan does not decrease as 1/ Ntot , as expected for a one-field Monte Carlo dose calculation. The noise convergence error appears as an indirect consequence of the noise, because the inverse treatment planning algorithm converges to the ‘optimal’ solution for the noisy beamlets, which is different from the optimal solution for the noise-free beamlets. To estimate the noise convergence error of the plan one needs to know how much is the difference between the recalculated dose distribution and the optimal dose distribution. Hence the noise convergence error for the plan, $conv , was determined by subtracting the optimal dose distribution, obtained by using the noise-free beamlets, from the recalculated dose distributions: ∞ $conv = dN∞ − d∞ = (wN − w∞ )B ∞ .

(11)

In this way, the statistical error was eliminated. Because the solution of the inverse problem is obtained through optimization, it is not possible to predict the magnitude of the noise convergence error in general. However, as long as the statistical error is small, the shape of the solution space (defined with the objective function) is only slightly perturbed. Therefore one might expect that noise convergence error approximately follows the statistical error. On the other hand, a large statistical error may significantly change the shape of the solution space, and hence it may be expected that the noise convergence error starts to diverge. Two methods have recently been proposed to reduce the effects of statistical uncertainty on dose volume histograms (Sempau and Bielajew 2000, Jiang et al 2000). Though these methods may be advantageous in some situations, they will not reduce the noise convergence error. Similarly, any noise removal procedure (for example as proposed by Deasy (2000)) would avoid the problem of the noise convergence error only if the statistical uncertainty was removed completely and no systematic error was introduced by the denoising procedure. It should be mentioned that the noise convergence error is meaningful only if the convergence close to the global minumum is achieved, in other words, if the systematic error of the optimization algorithm is sufficiently small. 3. Method and materials The effect of Monte Carlo statistical uncertainty was examined on a lung tumour example. MCI, a Monte Carlo based inverse treatment planning algorithm (Jeraj and Keall 1999a) for photon beams, was used in the study. The geometry was specified on a 0.5 × 0.5 × 10 cm3 grid (MCI assumes invariant geometry in the superior/inferior dimension). The algorithm was directed to optimize seven equispaced 6 MV beam positions with 21 allowed beamlet (narrow beam) orientations at each position. The beamlet size was chosen to be 1 × 1 cm2 at the isocentre. EGS4 (Nelson et al 1985) with usercodes BEAM (Rogers et al 1995) and DOSXYZ (Ma et al 1995) were used to calculate the dose distribution for each beamlet. MCI uses fast simulated annealing (Szu and Hartley 1987) as the optimization algorithm.


3605

Various objective functions were chosen for the comparison, in either of the following forms, to estimate their effect on both the statistical and noise convergence errors. One objective function is in the square form, with quadratic penalties for all tissues " CT " OF = (di − DT )2 (12) N T T i∈T while in the second form, the quadratic penalties in the tumour are replaced with exponential penalties on the basis of the linear-quadratic model (Fowler 1989) " CT * " CT " |di −DT | (e − 1) + (di − DT* )2 . (13) OF = * NT i∈T N T * * T i∈T CT is the weight of a tissue type T and di is the dose in a voxel i of type T . DT is the desired dose in a tissue type T (required mean dose in tumour and zero for other tissue types) and NT the number of voxels of type T . In the exponential/quadratic form, T represents tumour tissue and T * all other tissues. Throughout this work the normal and critical structures are penalized using square objective functions. To differentiate the objective functions used in this work, those in the form of equation (12) are referred to as quadratic(tumour), while those in the form of equation (13) are referred to as exponential(tumour). The weighting factors, CT , were generally kept to 1, except when the critical tissue was spared where CT * was set to 10. Cases with different types of objective functions were examined because the ‘optimal dose distribution’ depends very much on the chosen objective function and convergence is sensitive to the shape of the objective function in the vicinity of the global minimum. In addition to the above so-called soft constraints, hard constraints were also imposed. The dose in the tumour was required to be within the ICRU limits of 95–107% (ICRU 1993) and in the case of sparing, the upper dose limit was set for the critical structures. Three different objective functions were examined in our study. In the first case the exponential(tumour) objective function was used, and no critical structure sparing was requested. In the second case (exponential(tumour) with sparing), the same objective function was used as in the first case, except that the right lung was spared by limiting the maximum dose in the lung so as not to exceed 60% of the maximum dose. In the third example, the objective function for the target was changed into the square form, quadratic(tumour). No critical structure sparing was required. It should be mentioned that in our study only two types of objective functions with predefined importance factors and physical (dose) constraints have been considered. Alternatively, the importance factors for different structures could be optimized with a higher level optimization as in Xing et al (1999) or biological rather than physical objectives used, as in Niemierko (1997) or Söderström and Brahme (1993). Therefore, our study should be seen as an example revealing general errors that appear when using Monte Carlo dose calculations in inverse treatment planning. To establish correlation between the statistical uncertainty of the Monte Carlo calculations and optimization results, beamlets B N were calculated to different levels of precision, giving approximate absolute statistical errors from 0.5 to 4% at Dmax . The results obtained with noisy beamlets were compared with the results calculated with (almost) noise-free beamlets B ∞ , which were calculated to very high precision with the absolute statistical error σstat less than 0.2% at Dmax . The test for determining convergence to the global minimum of the objective function for the optimization algorithm used (simulated annealing) was made. It was performed by optimizing dose distributions for the exponential(tumour) without a sparing objective using two different initial weight vectors and using noise-free beamlets. The first case initially used

3606

R Jeraj and P Keall

Figure 1. Optimal dose distribution for the exponential(tumour) objective function without sparing with indicated geometry structures.

uniform weights for all beamlets, and the second initially used the weight vector calculated to be optimal for noisy beamlets. Optimization was then performed to determine the final weights, and the two final distributions were compared. Similar intensity maps and thus dose distributions show convergence to the same minimum of the given objective function. Though this test is by no means exhaustive or complete, nevertheless similar distributions indicate that the optimization method (simulated annealing) converged to near the global minimum. The standard deviation of the differences gives also an estimate of the systematic error of the optimization algorithm used in MCI. 4. Results As an example of our inverse treatment plan the dose distribution for the case with ∞ exponential(tumour) objective function without sparing using d∞ is shown in figure 1. 4.1. Statistical error The effect of statistical uncertainty of the Monte Carlo calculated dose distributions on inverse treatment planning was examined by using beamlets calculated to different precision in the MCI algorithm. To determine the statistical errors, all noisy dose distributions, dNN , were recalculated N with noise-free beamlets. It is interesting to note that the recalculated dose distributions, d∞ , revealed several pitfalls, the most severe of which was underdose in the tumour in some cases. An example of the original and recalculated distributions for the exponential(tumour) case and the corresponding dose volume histograms is presented in figure 2. Note that the largest differences are seen within the tumour volume, while in the normal tissue and critical structure (lung) the differences are almost negligible. Normal tissue is regarded as all tissue other than the tumour and lung. Changes in meeting the posed requirements were additionally examined by monitoring the minimum and maximum dose to the tumour of the originally optimized and recalculated plans. The results for different levels of statistical uncertainty are presented in table 1. Note that while being successfully optimized to satisfy the required uniformity between 95 and 107%, the recalculated dose distributions reveal under/overdose. The increased precision of calculation (more Monte Carlo histories simulated) reduces the observed under/overdose.


3607

(a)

(b)

(c) Figure 2. Original dNN (a) and recalculated dN∞ (b) dose for the exponential(tumour) without sparing objective function with approximately 3% statistical error at Dmax . Note smoothing of dose distributions and underdose in some cases (the light region in the tumour, which indicates a region below 95%). In the corresponding dose volume histograms (c) of the original (full curve) and recalculated (dotted curve) dose distributions note how, with noisy beamlets apparently successfully optimized, the recalculated results are indeed unacceptable because the dose coverage of the tumour is insufficient.

3608

R Jeraj and P Keall Table 1. Minimum and maximum dose to the tumour volume. The numbers represent min/max dose in the tumour when the plans were recalculated with noise-free beamlets, dN∞ . Numbers in parenthesis represent min/max dose in the tumour in the original inverse treatment plans dNN . Note that while being successfully optimized to satisfy a required uniformity between 95 and 107%, recalculated dose distributions reveal under/overdose. The quoted errors are approximate statistical errors at Dmax . σstat at Dmax

Minimum tumour dose

Maximum tumour dose

Noise-free 0.5% 1% 2% 3% 4%

95.0 93.8 (95.0) 92.4 (95.0) 91.7 (95.0) 87.0 (94.1) 90.0 (92.5)

106.9 107.2 (106.9) 107.1 (106.8) 109.0 (106.9) 108.4 (107.0) 109.0 (107.0)

These differences between the original and recalculated dose distributions are presented in figure 3. These differences represent statistical uncertainty in the dose calculation. It can be observed that the differences increase as the statistical uncertainty of the dose calculation (beamlets) increases, and that the errors are larger towards the tumour volume as theorized in the introduction. This increase occurs because the dose delivered to the tumour is on average higher than the dose delivered elsewhere. Examination of the cases where different objective functions were employed showed that the statistical error of a plan is rather insensitive to the objective function used because of the overall similarity of the dose volume histograms of different plans. 4.2. Convergence error The convergence to the global minimum was tested before evaluation of the noise convergence error of the inverse treatment plan. The treatment plan calculated with noisy beamlets was reoptimized with noise-free beamlets. The resultant plan was then compared with the originally obtained plan with noise-free beamlets. The histogram of the differences between both plans is presented in figure 4. It is clear that the convergence error of the reoptimized plan is much lower, especially in the tumour volume, where in more than 90% of the tissue the differences are smaller than 1%. Outside the tumour volume, more than 90% of the tissue has differences smaller than 2%. The established convergence error (standard deviation) of approximately 0.5% can be ascribed to the systematic error of MCI for the optimization parameters as used in the study (number of iterations in simulated annealing optimization, cooling schedule). The differences for the beamlets calculated to different precisions indicating the noise convergence error are presented in figure 5. Differences indicated in figure 5 should be understood as differences in the optimized intensity profiles given in terms of the dose distribution. From the results it can be noted that the differences increase as the statistical uncertainty of the dose calculation (beamlets) increases, and that the errors are generally larger outside the tumour volume. These trends are because the dose in the tumour is forced to be uniform while there is no such strict requirement outside. The fact that the noise convergence error strongly depends on the form of the objective function, which will be shown in the following section, is also very important. A similar convergence error would also be present if inaccurate calculations (those burdened with systematic error) were used (only imprecise calculations (those burdened with statistical error) are studied here), as the weights determined for the inaccurately calculated beamlets would be different from the weights determined for the accurately calculated beamlets.


(a)

(b)

(c)

(d)

3609

Figure 3. Differences between the original and recalculated dose distributions $stat indicating the statistical error of the plans. Approximate statistical errors of the beamlets are (a) 0.5%, (b) 1%, (c) 2% and (d) 4% at Dmax . Note that the differences are larger in the tumour volume than outside because of the higher dose delivered in the tumour. The figure shows the absolute differences normalized to the maximum dose in the plan.

Figure 4. Test of the convergence to the global minimum. The noise convergence error of the reoptimized plan (with noise-free beamlets) is approximately 0.5% (right), compared with almost 5% of the original plan (left). This indicates successful convergence near the optimum for the given objective function.

4.3. Dependence of errors on the Monte Carlo statistical uncertainty level As has been seen, the statistical uncertainty introduces two errors, statistical and convergence, which depend on the number of simulated histories. However, they are also dependent on the specified objective function. Both errors were examined for three different objective functions as functions of the number of simulated particle histories in each beamlet. The cases differed

3610

R Jeraj and P Keall

(a)

(b)

(c)

(d)

Figure 5. Differences between the recalculated and optimal dose distributions $conv indicating the noise convergence error of the plans. Approximate statistical errors of the beamlets are (a) 0.5%, (b) 1%, (c) 2% and (d) 4% at Dmax . Note that the errors are larger outside the tumour volume. The figure shows the absolute differences normalized to the maximum dose in the plan.

in the penalty function for the tumour (exponential or quadratic) and whether the critical tissue (right lung) was spared or not. The results are presented in figure 6. Examination √ of the results shows that the statistical error of the plan has approximately the expected 1/ N behaviour, where N is the number of simulated particle histories. Also notable is that the statistical error is not dependent on the form of the objective function, while the noise convergence error is very sensitive to the specified objectives. It should be mentioned that even though the dose volume histograms (DVHs) for exponential(tumour) and quadratic(tumour) objective functions are very similar, the behaviour of the noise convergence error is very different. It is also interesting to note that when the statistical uncertainty in the dose calculation becomes significant (of the order of 2%), the optimization solution starts to diverge considerably from the optimal solution, whereas when the statistical uncertainty in the dose calculation is moderate (below approximately 2%), the noise convergence error approximately follows the statistical error. The numbers given are only approximate and depend on the objective function used. The monotonic behaviour of the noise convergence error in all cases in figure 6 infers that optimization close to the global minimum was achieved for each case. If the solution for some situations was trapped in a local minimum, we would expect larger variations in the behaviour of the noise convergence error as a function of the number of particles per beamlet. The convergence of both errors towards zero implies the expected conclusion that when the statistical uncertainty of Monte Carlo dose calculation is small enough, the solutions of the inverse treatment planning can be considered to be the same as those which would be obtained with idealized calculations without statistical error.


3611

Figure 6. Statistical and noise convergence errors of the plans as functions of the number of simulated particle histories. The errors shown here are standard deviations of the differences in the plans as shown in figures 3 and 5. Top: exponential(tumour) objective function. Middle: exponential(tumour) objective function with right lung sparing. Bottom: quadratic(tumour) objective function.

5. Conclusions When using a Monte Carlo-driven dose calculation, one should always consider statistical uncertainty, which arises because of the finite number of simulated particle histories. For this reason, in Monte Carlo based conventional treatment planning, all the results are burdened

3612

R Jeraj and P Keall

with some statistical error. However, when using Monte Carlo transport to calculate optimal dose distributions in inverse treatment planning, another error is introduced. This is a so-called noise convergence error. In the current work the statistical and noise convergence errors for Monte Carlo based inverse treatment planning were studied by analysing results from plans optimized with both different objective functions and beamlet (narrow beam) statistical uncertainty levels. It was determined that the choice of the objective function significantly affects the noise convergence error, but has little impact on the statistical error. The effect of statistical uncertainty on inverse treatment planning was first examined by comparing optimized plans with noisy beamlets and recalculated plans with noise-free beamlets. The recalculated plans revealed several pitfalls, the most severe of which was underdose in the tumour in some cases. Apparently successfully optimized plans with noisy beamlets, which satisfied the required uniformity of the dose within the tumour, showed an underdose of up to 7% when recalculated. The statistical error in inverse treatment plans appears because of the statistical nature of the Monte Carlo calculation used to calculate the dose distributions of individual beamlets. It was observed that the differences due to statistical error are larger within the tumour volume because of the higher dose delivered there than outside. Examinations of cases where different objective functions were employed showed that the statistical error is rather insensitive to the √ objective function used and displays the expected 1/ N behaviour, where N is the number of simulated histories. The convergence error is a specific feature of any inverse treatment planning method which uses either inaccurate (e.g. pencil beam) dose calculation or imprecise (e.g. Monte Carlo) dose calculation. When using Monte Carlo based inverse treatment planning, the noise convergence error occurs because the inverse treatment planning algorithm converges to the global minimum solution for the noisy beamlets, which is different from the global minimum for the noise-free beamlets. From our study it was concluded that the absolute differences due to the noise convergence error are generally larger outside the tumour volume, because the dose in the tumour is forced to be uniform, while it is not so strictly prescribed outside the tumour. It should be stressed that the noise convergence error strongly depends on the objective function used and is case dependent. For any clinical implementation of Monte Carlo based inverse planning, a similar study to the one performed here should be conducted spanning several different tumours and treatment sites and using the objective functions that the institution uses for IMRT. However, the question remains as to how well the objective functions describe the clinical objectives. It may be that the noise convergence error is lower than the error due to the use of a given objective function. The difficulty in determining the objective function for inverse treatment planning is the variability in the response of the population to radiation, as well as inaccuracy and lack of reliable clinical data. Also important in such discussion is the magnitude of the systematic error of the dose calculation algorithm, and its relation to the statistical and noise convergence error as well as the systematic error of the optimization procedure employed in the inverse treatment planning.

Acknowledgments Robert Jeraj acknowledges financial support from Ministry of Science and Technology, Slovenia, and Paul Keall acknowledges the financial support of NIH Grant CA 74158-02.


3613

References Buffa F, Nahum A and Mubata C 1999 Influence of statistical fluctuations in Monte Carlo dose calculations on radiobiological modelling and dose volume histograms (abstract) Med. Phys. 26 1120 Deasy J 2000 Is Monte Carlo noise removal feasible? The Use of Computers in Radiation Therapy ed W Schlegel and T Bortfeld (Heidelberg: Springer) pp 403–5 Fowler J F 1989 The linear-quadratic formula and progress in fractionated radiotherapy Br. J. Radiol. 62 679–94 ICRU (International Commission on Radiation Units and Measurement) 1993 Prescribing, recording and reporting photon beam therapy ICRU Report 50 (Washington, DC: ICRU) Jeraj R 1999 Monte Carlo based inverse treatment planning PhD Thesis University of Ljubljana Jeraj R and Keall P J 1999a Monte Carlo based inverse treatment planning Phys. Med. Biol. 44 1885–96 ——1999 Effect of noise on Monte Carlo based inverse treatment planning (abstract) Med. Phys. 26 1099 Jiang S B, Pawlicki T and Ma C-M 2000 Removing the effect of statistical uncertainty on dose-volume histograms from Monte Carlo dose calculations Phys. Med. Biol. 45 2151–61 Keall P J, Siebers J V, Jeraj R and Mohan R 2000 The effect of dose calculation uncertainty on the evaluation of radiotherapy plans Med. Phys. 27 478–84 Ma C-M, Reckwerdt P, Holmes M, Rogers D W O, Geiser B and Walters B 1995 DOSXYZ users manual Report PIRS-0509(b) (Ottawa: NRCC) Munzenrider J E, Brown A P, Chu J C H et al 1991 Numerical scoring of treatment plans Int. J. Radiat. Oncol. Biol. Phys. 21 147–63 Nelson W R, Hirayama H and Rogers D W O 1985 The EGS4 code system SLAC Report SLAC-265 (Stanford, CA: SLAC) Niemierko A 1997 Reporting and analyzing dose distributions: a concept of equivalent uniform dose Med. Phys. 24 103–10 Rogers D W O, Faddegon B A, Ding G X, Ma C-M, We J and Mackie T R 1995 BEAM: a Monte Carlo code to simulate radiotherapy units Med. Phys. 22 503–24 Sempau J and Bielajew A 2000 Towards the elimination of Monte Carlo statistical fluctuation from dose volume histograms for radiotherapy treatment planning Phys. Med. Biol. 45 131–57 Söderström S and Brahme A 1993 Optimization of the dose delivery in a few field techniques using radiobiological objective functions Med. Phys. 20 1201–10 Szu H and Hartley R 1987 Fast simulated annealing Phys. Lett. A 122 157–62 Xing L, Li J G, Donaldson S, Le Q T and Boyer A L 1999 Optimization of importance factors in inverse planning Phys. Med. Biol. 44 2525–36