Short title. Monte Carlo perturbation methods in radiotherapy

A Monte Carlo based perturbation method for electron radiotherapy planning

R. P. Hugtenburg
Department of Medical Physics, Queen Elizabeth Medical Centre, University Hospital Birmingham NHS Trust, Birmingham B15 2TH, UK

Abstract. A perturbation method for radiotherapy treatment planning (RTP) is described that determines departures from measured dose distributions caused by variations in treatment set-up and patient anatomy. The calculation is a Markov chain Monte Carlo method, the efficiency of which depends on the degree of a priori knowledge of the radiation source and the appropriateness of the measured data. Monte Carlo generated dose distributions are calculated to emulate in phantom measurements. Concurrently an in vivo dose distribution is calculated and correlated sampling is used if the local geometry of the particles is comparable. The method is tested for an obliquely incident beam with good agreement to measurements made with a diode detector. The method is demonstrated for a test object constructed from connected facets. The calculations show that convergence to the in phantom measurement ensures an adequate enough representation of the beam to compute an in vivo dose. Measured data, such as would be acquired during the commissioning of a therapeutic linear accelerator, is shown to adequately characterise the radiation source with basic prior assumptions.

Introduction

The aim of a perturbation method is to use a measured dose distribution to provide a first approximation to the dose distribution in the patient and to infer the source characteristics of the beam so as to enable higher order corrections to be made. This is made possible using a Markov chain method to determine a set of most likely phase space contributions to the beam. Electron radiotherapy planning is complex due to the strong interactions of electrons with solid materials leading to complicated scatter effects in the accelerator componentry and patient. The Monte Carlo method is an increasingly practical tool for electron radiotherapy problems.

In Monte Carlo based radiotherapy planning [1,2,3], an in vivo dose distribution is computed without reference to the commissioned data and so requires detailed knowledge of the source configuration and careful benchmarking. Concurrent computing of an in phantom dose distribution is feasible in that there will frequently be correlations with the in vivo dose distribution over proportions of the volume. The correlated sampling technique used in this context [4] will lead to a speed up in the calculation depending on the degree of correlation between the measured and patient scenarios and even immediate convergence if the patient set up is equivalent to the measurement conditions.

Correlated sampling

Source particles are transported through both the regular phantom used for data commissioning (typically a semi-infinite geometry) and the patient anatomy. Firstly and importantly, this enables the evaluation of a perturbation to the absolute dose.

The computation time will depend on the degree of geometric correlation of the patient set-up with that of the measurement data set. Electron planning prescriptions are often limited to a small set of variations of the measured data set and a department may be inclined to measure output factors for any new field shape. The algorithm is most efficient where the treatment is a small variation on the measured data.

Correlated sampling for radiotherapy treatment planning as introduced by Holmes et al takes the form,

$$D_v = D_0 \frac{\sum_n D_v^n}{\sum_n D_0^n},$$ (1)

where D₀ is a measured distribution. D_vⁿ, D₀ⁿ are the nth histories for Monte Carlo computed distributions in vivo (v) and in phantom (0). In the limit, $n\rightarrow\infty, \sum_n D_0^n/n \rightarrow D_0$, so a variant on this expression,

$$D_v = D_0 + \left(\sum_n D_v^n - \sum_n D_0^n\right)/n,$$ (2)

has also been explored and is advantageous where regions of low dose and high variance in the denominator of equation 1 will compromise the convergence. The correlated sampling method enables the transport of radiation in vivo to be treated as a perturbation of transport in phantom. Where possible, the local environments of the incident particles are arranged to be equivalent. For example, incident particles are translated to the surface of the patient where the air/water equivalent interface occurs. Locally this is the same as the surface of the phantom. Dose deposition occurs at the translated point and will not necessarily correspond spatially with the phantom but will still contribute to the overall coherency of the distribution.

Source characterisation.

A Markov chain Monte Carlo method [5] has been used to obtain the distribution of incident particles based on the measured data. The measured data should therefore be beam characterising and relevant to the range of patient set-ups. A typical commissioning data set usually fulfils these two criterion.

Each consecutive result forms a Markov chain where the transition to the next state is dependent wholly on the current state. The randomness in the summation of dose serves to drive the convergence of the accumulated dose towards that of the measured distribution by rejecting those contributions which are the less likely to contribution. The likelihood of a particle contributing to the dose distribution, determined to within a confidence interval, is dictated by the sum of squares of differences between the two sets compared to the variance.

The finite number of histories used per iterative step introduces a bias towards outlying results that achieve statistical significance and is a function of the random error (variance) of the step. It is desirable therefore to have steps that contain a large batch of histories. The greater the size of the batch the better the sensitivity of the sampling to the actual beam fluence.

The acquisition of a multidimensional phase space set is served well by a Markov chain-Monte Carlo method. An alternative is to build a discrete array of Monte Carlo generated dose kernels d_ij for the expected range of energy and angle of incidence. A set of incident particles can be computed by solving an equation of the form,

$$D_i = \sum_j d_{ij}{w_j},$$ (3)

where a set of weightings, w_j, are then computed for a set of jparticles that best achieve a distribution, D_i, the set of dose measurements. An inversion method such as the maximum likelihood method [6,7], simulated annealing or simplex [8] will enable a solution set of w_j's to be determined for this equation where w_j>0 for all j. Such techniques are collectively known as the non-iterative inverse Monte Carlo method [9,10]

With a Markov chain method it is not necessary to discretise the fluence distribution nor store a dose kernels. This is particularly desirable if the incident particles have a large number of degrees of freedom (in which case j would need to be large). In the case of electron radiotherapy, the distribution of particle trajectories and energies is a complicated function of position. None-the-less, several workers have used these techniques to analyse spectral content in a broad-parallel electron beam [11,12,13].

A full in phantom dose distribution such as would be acquired during a typical commissioning of a therapeutic linear accelerator has been used for the purposes of beam characterisation. An electron beam has a considerable scatter component which is evident if measurements are performed over a range of source-surface distances (SSD) and leads to an effect that gives a reduced effective source distance (ESD) [14]. In fact electrons scattered in the applicator cones can acquire trajectories that are highly oblique and even convergent.

To resolve electron trajectories more precisely, in phantom data has been acquired over a range of SSD. The influence of this extra data can be introduced into the calculation very simply and does not require the transport of particles through the phantom at all of the measured SSD. If the electrons are assumed to maintain a reasonably straight trajectory through the intervening air volume, a small lateral translation can be used to determine the dose deposition at the extended SSD.

The translation to the point of intersection, $\delta{\bf x^{l}}$ is given by

$$\delta{\bf x^{l}}=(l-l_0)\frac{{\bf u}}{u_z},$$ (4)

where l is length of the extended SSD, l₀ is the standard SSD and ${\bf u}=(u_x,u_y,u_z)$ is the direction of the incident particle as a cosine to axes x, y and z.

The particles then begin their transport at z=0 for all the SSD. The translation in the x and y positions need only be considered for the sake of energy deposition as the phantom is semi-infinite and the transport independent of their values.

Method

Commissioning of measured data

A set of data has been acquired for 10 MeV electrons on the Philips SL15 with a 10$\times$10 cm² applicator cone using a diode detector (Scanditronix). While practical electron dosimetry requires some information about the energy of the electrons as a function of position (depth) and is often ill-defined, here, the energy deposited can be immediately converted into an ionisation defined by the material constituency of the detector.

  

Figure: Measured percentage depth dose at 96, 100 and 105 cm SSD, for 10 MeV electrons with a 10$\times$10 cm² applicator cone. The peak dose (22 mm) as a function of SSD are 97.5%, 88.4% and 78.7%, respectively, of the peak dose for zero stand-off (95 cm SSD) (an ESD of 76.3 cm).

The measured data set consists of percentage depth dose (PDD) and off-axis ratios (OAR) at three phantom positions 95 cm, 100 cm and 105 cm SSD. The PDD are measured in 2 mm graduations to a depth of 100 mm (figure 1). OAR are measured in a single plane in 2 mm graduations to a range of 10 cm from the central axis at depths of 10, 20, 30, 40 and 50 mm (figure 2). A total of 918 measurement points seems excessive for a single applicator cone and energy but is comparable to a typical data set required for photon planning. The acquisition of this data is easily automated with a field analyser (RFA-300, Scanditronix).

  

Figure: Measured (peak) off-axis ratios at 96, 100 and 105 cm SSD, for 10 MeV electrons with a 10$\times$10 cm² applicator cone.

Implementation of the Monte Carlo method.

The Monte Carlo transport has been implemented on EGS4 ( ESTEPE = 0.01, PCUT = 0.100 MeV, ECUT = 0.611 MeV) with the FIXTMS option. Transport within the phantom is straightforward with a single surface boundary and cut-offs beyond 20 cm depth and radius. The patient anatomy has been implemented using sets of connected facets to form a closed surface. Anatomical structures are also represented by enclosed volumes and variations in density are catered for.

The transport within the volumes is uninhibited by the dose scoring grids. Instead deposited energy ( EDEP) is binned according to the location. A potential problem occurs with the condensed history method used for electrons. In EGS4 the dose is deposited at the beginning of the transport step. This will lead to biases in the deposition particularly if there is a correlation between the surface position and the dose grids. The problem is resolved by depositing dose at a random position along the electron step-length as determined by TVSTEP, a variable that can be interrogated in the AUSGAB energy scoring routine by including the common block, EPCONT.

Initialisation

The energy and direction of the incident particle is sampled from a prior domain. The energy is simply sampled from a uniform distribution up to a maximum energy of 12 MeV. The trajectory of the incident particle is sampled in the range $-\pi/4$ to $\pi/4$ from the central axis. A more complicated prior can be used to limit the sampled domain to a detailed description of the beam such as would be provided by a beam model [2]. Bremsstrahlung contamination is observed in practice and so photons are also sampled. More highly variant particles such as photons should be sampled more frequently to avoid longer than necessary 1convergence times.

An intersection point is determined for the particle incident upon the patient. This is a useful speed-up because if the air gap is modelled the electrons will interact several times in an intervening volume. The interaction serves to destroy the correlation but will influence the trajectory of the electron only slightly. The particle therefore begins its transport at the interface of the surface which is in the first approximation a planer geometry akin to the phantom.

The output factor of a blocked field which isn't already a part of the measured data set can be calculated using this method. In this case only particles from the in-phantom data set are transported and are effectively subtracted from the measured data set.

Energy deposition

The measured dose distribution is worked into a set of bins, $D_{0,l}(\bf x_i)$, as a total absorbed energy. The deposited energy from the convergent Markov chain Monte Carlo computation of the dose distribution, $\sum_nD^n_{0,l}(\bf x_i)$, is binned accordingly for each SSD. It has been useful to work the measured data into cumulative sums that correspond to large volumes of relatively uniform dose of approximately equal orders of magnitude.

An octet symmetry of uniform square fields has been assumed for the purposes of simplification. The bins consist of enclosed rectangles of width 20 mm, 36 mm, 48 mm, 56 mm, 60 mm, 64 mm ... out to 200 mm, a total of 40 regions of reasonably equivalent volume.

A ``burn-in'' period will be necessary where data is rejected until such a point that equilibrium in the Markov chain sequence has been achieved. The energy deposition is then recorded in two new grids for an in phantom dose distribution at an appropriate SSD and the in vivo dose distribution.

Particle splitting and the in vivo geometry

When a particle reaches a boundary that is not shared by the two geometries, the particle is split into two, each will now follow separate paths. The particles are given separate identities via a flag, LATCH. In the simplest circumstance, any interaction with a surface will mean that one particle is leaving the geometry. The particle is recorded as a discard in the geometry it is leaving but continues to be transported in the other.

Two test in vivo geometries have been examined. The first is a 45⁰ oblique surface. This is a typical problem associated with electron radiotherapy set-ups and a dose enhancement effect has been well reported [15,16,17]. Calculations will be compared with a set of measurements. The boundary crossing procedure in this circumstance is trivial.

  

Figure: Test object - a 20 cm³ cube with one corner missing. Half of a 10$\times$10 cm² incident beam projects onto a normal surface the other half projects onto the oblique section.

The second example is an enclosed volume, which in a general form can be used to describe the patient. A cube with one corner removed (figure 3) has been considered with the beam incident upon the object as shown. Half of the beam intersects with a perpendicular face while the other half is incident on a highly oblique (45⁰) face.

Boundary crossings are determined in an enclosed volume from the intersection of the particle trajectory with a set of three pointed facets defining the surface of the volume. A patient geometry could be formed from a set of separations or taken from a series of anatomical CT slices. An algorithm for the facetting of sequential CT slices has been composed in-house. With larger numbers of facets, an heuristic based on direction should be used to eliminate unlikely intersections.

Results

  

Figure: Percentage depth dose for a 10 MeV electron beam with a 10$\times$10 applicator cone, incident at 45 degrees, calculated by the perturbation Monte Carlo method.

Figure 4 shows the depth dose distributions calculated for a 45⁰ oblique incidence compared to normal incidence. This and the previously acquired data are compared to measurements at several depths made with a diode detector in a water phantom. The normalisation for the Monte Carlo calculation and measurement are relative to the measurements at normal incidence and so offer a realistic impression of the dose enhancement effects described previously. The full calculation takes 13 hours on a Sun Spark Ultra II however only 12.3% of incident particles sampled are used which amounts 1.6 hours of computation time if the phase space data set is pre-computed.

  

Figure: Dose distributions on a 0.5 cm resolution grid calculated for 10 MeV electrons with 10$\times$10 applicator cone blocked to an 8$\times$8 field impinging on the test object. 18 million histories are generated in 13 hours and 40 minutes on a Sun Spark Ultra II. Statistical uncertainty at peak is 2% (2sd).

Figure 5 presents the dose distributions computed for the object shown in figure 3 using direct Monte Carlo computed from the acquired phase space distribution and via equations 1 and 2. 180000 incident particles were acquired using the Markov chain method with 100 histories per iterative step. The calculation time of 13 hours 40 minutes doesn't include particles rejected during the Markov chain process (about 90%). The statistical error at peak is comparable for each of the three cases although with correlated sampling (equations 1 and 2) the regions that are common to both the phantom and in vivo geometries, i.e. the upper right hand side, have negligible statistical variance. The limited size of the phase space data set leads to spikes that are most pronounced in the direct Monte Carlo calculation but are cancelled to some extent by correlated sampling. Equation 1 leads to large statistical variances in dose distribution beyond the practical range of the electrons in phantom.

Several techniques have been used to improve the degree of correlation in the two geometries. In the case of dose distribution computed for an oblique plane, the perturbation is treated as a rotation about the isocentre ensuring that many particles along the central axis will deposit energy within the same cell. For a general anatomy, the translation of particles to the surface of the patient enables transport to be performed in tandem and allows spatial correlations to occur where the translation is small. Non-homogeneous anatomies could be considered by scaling individual track-segments to maintain correlations. This has yet to be investigated but shares similarities with the superposition Monte Carlo algorithm [18] that scales track lengths and angles of deflection in the kernel according to relative stopping powers of the local tissue type.

Conclusion

The use of a correlated sampling method enables a measured set of data to contribute to the convergence of a Monte Carlo calculation. Computation times of several hours for an electron dose distribution are feasible using either the direct Monte Carlo method or perturbation Monte Carlo method with the pre-computation of a phase space data set. The perturbation method will achieve faster convergences depending on the degree of correlation between the measured data and the treatment set-up. Calculations of obliquity, applicator and stand-off factors can be performed even more rapidly and in a time frame that will be useful in a clinic.

The use of a Markov chain Monte Carlo method in combination with a set of measured data, such as would be acquired during the commissioning of a therapeutic linear accelerator, offers an alternative to detailed modelling of the accelerator head. Faster convergence and acquisition rates are achieved with a better knowledge of the source characteristics and beam modelling is therefore complimentary.

Bibliography

1

D. W. O. Rogers, B. A. Faddegon, G. X. Ding, C. M. Ma, J. Wei, and T. R. Mackie, ``BEAM: A Monte Carlo code to simulate radiotherapy treatment units,'' Med. Phys., vol. 22, no. 5, pp. 503-524, 1995.

2

C. M. Ma, B. A. Faddegon, D. W. O. Rogers, and T. R. Mackie, ``Accurate characterization of Monte Carlo calculated electron beams for radiotherapy,'' Med. Phys., vol. 24, no. 3, pp. 401-416, 1997.

3

A. Kapur, C.-M. Ma, E. C. Mok, D. O. Findley, and A. L. Boyer, ``Monte Carlo calculations of electron beam output factors for a medical linear accelerator,'' Phys. Med. Biol., vol. 43, pp. 3479-3494, 1998.

4

M. A. Holmes, T. R. Mackie, W. Sohn, P. J. Reckwerdt, T. J. Kinsella, A. F. Bielajew, and D. W. O. Rogers, ``The application of correlated sampling to the computation of electron beam dose distributions in heterogeneous phantoms using the Monte Carlo method,'' Phys. Med. Biol., vol. 38, pp. 675-688, 1993.

5

W. K. Hastings, ``Monte Carlo sampling methods using Markov chains and their applications,'' Biometrika, vol. 57, no. 1, pp. 97-109, 1970.

6

J. Llacer, ``Inverse radiation planning using the dynamically penalized likelihood method,'' Med. Phys., vol. 24, no. 11, pp. 1751-1764, 1997.

7

G. H. Olivera, D. M. Shepard, P. J. Reckwerdt, K. Ruchala, J. Zachman, E. E. Fitchard, and T. R. Mackie, ``Maximum likelihood as a common computational framework in tomotherapy,'' Phys. Med. Biol., vol. 43, pp. 3277-3294, 1998.

8

W. H. Press, Numerical recipes: The art of scientific computing.
New York: Cambridge University Press, 1986.

9

W. L. Dunn, ``Inverse Monte Carlo analysis,'' J. Comp. Phys., vol. 41, pp. 154-166, 1981.

10

W. L. Dunn, V. C. Boffi, and F. O'Foghludha, ``Applications of the inverse Monte Carlo method in photon beam physics,'' Nucl. Instrum. Meth. A, vol. 255, pp. 147-151, 1987.

11

M. D. Altschuler, P. Bloch, E. L. Buhle, and S. Ayyalasomayajula, ``3D dose calculations for electron and photon beams,'' Phys. Med. Biol., vol. 37, no. 2, pp. 391-411, 1992.

12

R. Martinelli, G. Ghiso, P. Boccacci, F. Foppiano, and L. Andreucci, ``Determination of the energy of the electrons incident upon a water phantom to build composite energy beam kernels for 3D dose calculation,'' in Proc. Int. Conf. Comp. Radiat. Ther. XI, (Manchester), pp. 124-125, 1994.

13

R. P. Hugtenburg, J. R. Turner, and P. H. Butler, ``Numerical reconstruction of electron therapy beams,'' Med. Phys., vol. 22, no. 6, p. 948, 1995.

14

S. P. Baggerley, D. M. Mannering, R. P. Hugtenburg, and J. R. Turner, ``Field size dependence of effective source distance and field factors of a 9 MeV electron beam,'' Australian Physical and Engineering Sciences in Medicine, vol. 17, no. 2, pp. 84-87, 1994.

15

P. J. Biggs, ``The effect of beam angulation on central axis percent depth dose for 4-29 MeV electrons.,'' Phys. Med. Biol., vol. 61, no. 4, pp. 465-472, 1984.

16

ICRU, Report 35. Radiation dosimetry: Electron beams with energies between 1 and 50 MeV.
International Commission on Radiation Units and Measurement, 1984.

17

P. M. Ostwald and T. Kron, ``Surface does measurements for highly oblique electron beams,'' Med. Phys., vol. 23, no. 8, pp. 1413-1420, 1996.

18

P. J. Keall and P. W. Hoban, ``Super-Monte Carlo: A 3-D electron beam dose calculation algoritm,'' Med. Phys., vol. 23, no. 12, pp. 2023-2034, 1996.

About this document ...

A Monte Carlo based perturbation method for electron radiotherapy planning

This document was generated using the LaTeX2HTML translator Version 98.1p1 release (March 2nd, 1998)

Copyright © 1993, 1994, 1995, 1996, 1997, Nikos Drakos, Computer Based Learning Unit, University of Leeds.

The command line arguments were:
latex2html -split 0 -no_reuse pertmc.

The translation was initiated by Richard Peter Hugtenburg - RTP on 2001-04-12