A fast and accurate treatment planning method for boron neutron capture therapy

Radiotherapy and Oncology 46 (1998) 321–332

A fast and accurate treatment planning method for boron neutron capture therapy Cornelis P.J. Raaijmakers*, Iaı¨n A.D. Bruinvis, Eric L. Nottelman, Ben J. Mijnheer Radiotherapy Department, The Netherlands Cancer Institute/Antoni van Leeuwenhoek Huis, Plesmanlaan 121, 1066 CX Amsterdam, The Netherlands Received 26 February 1997; revised version received 21 July 1997; accepted 10 September 1997

Abstract Background and purpose: The aim of this study was to test the applicability of conventional semi-empirical algorithms for the treatment planning of boron neutron capture therapy (BNCT). Materials and methods: Beam data of a clinical epithermal BNCT beam obtained in a large cuboid water phantom were introduced into a commercial treatment planning system (TPS). For the calculation of thermal neutron fluence distributions, the Gaussian pencil beam model of the electron beam treatment planning algorithm was used. A simple photon beam algorithm was used for the calculation of the gammaray and fast neutron dose distribution. The calculated dose and fluence distributions in the central plane of an anthropomorphic head phantom were compared with measurements for various field sizes. The calculation time was less than 1 min. Results: At the normalization point in the head phantom, the absolute dose and fluence values agreed within the measurement uncertainty of approximately 2–3% (1 SD) with those at the same depth in a cuboid phantom of approximately the same size. Excellent agreement of within 2–3% (1 SD) was obtained between measured and calculated relative fluence and dose values on the central beam axis and at most off-axis positions in the head phantom. At positions near the phantom boundaries, generally in low dose regions, local differences of approximately 30% were observed. Conclusions: A fast and accurate treatment planning method has been developed for BNCT. This is the first treatment planning method that may allow the same interactive optimization procedures for BNCT as applied clinically for conventional radiotherapy.  1998 Elsevier Science Ireland Ltd. Keywords: Boron neutron capture therapy; Treatment planning; Dosimetry; Epithermal neutrons; Pencil beam algorithm

1. Introduction Boron neutron capture therapy (BNCT) is a form of radiation therapy which is based on the high cross-section of the 10 B nucleus for capturing thermal neutrons compared with other elements present in tissue. After selectively loading the malignant cells with 10B an irradiation with thermal neutrons may result in tumour cell destruction by the heavy charged particles released in the 10B(n,a)7Li reaction, while healthy tissue without 10B is spared. The principles, history and recent developments of BNCT have been extensively described in various review articles [1,6,7]. Three epithermal neutron beams are currently available for clinical use of BNCT [15,23,26]. The energy of the

* Corresponding author.

neutrons in such beams ranges from approximately 0.5 eV to about 10 MeV. When travelling through tissue, epithermal neutrons (0.5 eV to 100 keV) loose their energy, mainly by the proton recoil process and reach thermal energies. As a result of this process, the maximum thermal neutron fluence in tissue is obtained at a depth of 2–3 cm. Besides the boron neutron capture dose several other undesired but nonnegligible dose components are present in boron-containing tissue irradiated with an epithermal neutron beam. Thermal neutrons interact not only with 10B but also with nitrogen and hydrogen present in all tissues. The particles resulting from the 10B(n,a)7Li and the 14N(n,p)14C reaction have ranges in tissue of less than approximately 10 mm. Consequently, the macroscopic dose from these reactions can be calculated from the local thermal neutron fluence and the local macroscopic concentrations of boron and nitrogen. Although the thermal neutron fluence is measured in this

0167-8140/98/$19.00  1998 Elsevier Science Ireland Ltd. All rights reserved PII S0167-8140 (97 )0 0183-7

322

C.P.J. Raaijmakers et al. / Radiotherapy and Oncology 46 (1998) 321–332

work, the expression boron neutron capture (absorbed) dose will also be used occasionally to refer to this component. Epithermal neutrons deliver dose through the proton recoil process. The main contribution to the dose produced by the proton recoil process, however, originates from fast neutrons (100 keV to 10 MeV) which are unavoidably present in the beam [21]. This dose component will therefore be referred to as the fast neutron dose. Gamma rays are unavoidably present in the beam but are also generated in tissue by neutron capture of hydrogen through the 1H(n,g)2H reaction. The contribution of the generated 2.2 MeV capture gamma rays to the total gamma-ray dose in a phantom is approximately 60–80% [21,26]. For the treatment planning of BNCT using an epithermal neutron beam, the distribution of three distinct components has to be calculated: (1) the thermal neutron fluence distribution, (2) the gamma-ray dose distribution and (3) the fast neutron dose distribution. The resulting dose components (boron neutron capture dose, nitrogen neutron capture dose, gamma-ray dose and fast neutron dose) have largely different relative biological effectiveness (RBE) values. In the case of the boron neutron capture dose this value depends on the compound used as boron carrier [18]. The various distinct dose components have different radiation transport characteristics. This fact and the assumption that explicit 3D calculations with a complete treatment of particle scattering are required have led to the opinion by various research groups that the conventional semiempirical approach for treatment planning is of limited use for BNCT [6,19,28]. The Monte-Carlo technique has become the major tool for treatment planning of BNCT. This technique is, in principle, capable of solving the complex radiation transport problems of BNCT but rests on the assumption that all parameters influencing the dose distributions including reactor core, beam tube, beam filters and the irradiated object can be modelled in sufficient detail. Three treatment planning systems (TPSs) for BNCT have been described, all based on the Monte-Carlo technique [25,27,28]. A drawback of these systems is that the results of detailed calculations are subject to relatively large statistical uncertainties unless an excessive amount of calculation time is used [17,27,28]. Calculation times are reduced at the cost of accuracy by decreasing the spatial resolution of the calculations and by ignoring electron transport [5,28]. With the development of faster computers, calculation times are becoming less problematic. Nevertheless, the calculation times currently needed to obtain a reasonable accuracy are in the order of hours and do not permit interactive treatment planning, i.e. comparing various beam configurations to achieve an optimum treatment plan [27,28]. Experimental validations of TPSs for BNCT are scarce and restricted to comparisons between measurements and calculations along the central beam axis. Experimentally determined correction factors are necessary to obtain agreement between measurements and calculations of the absolute dose and fluence values [28].

TPSs for conventional external beam radiotherapy are generally based on empirical beam data obtained under reference conditions, usually measured in large cuboid water phantoms. Calculation times for conventional TPSs are generally of the order of minutes or even seconds. A disadvantage of this approach is that the influence of phantom and beam geometry as well as phantom inhomogeneities on the dose distribution are taken into account using semi-empirical correction procedures which have their limitations. Calculations of dose distributions resulting from external photon beams are, in their most simple forms, based on interpolations of basic beam data. Semi-empirical algorithms correct for the differences between the reference conditions and the actual geometry. Generally, these corrections are performed along lines connecting the point of calculation to a virtual source. In these essentially 1D corrections, scattering properties are not fully taken into account. Consequently, such calculations are not well suited for the calculation of thermal neutron distributions, which are mainly determined by scattering events. More sophisticated algorithms for conventional treatment planning are based on pencil beam approximations. Although these algorithms are currently also used in photon beam TPSs, they were originally developed for the calculation of electron beam dose distributions [2]. The fluence distributions resulting from electron beams are highly influenced by multiple scattering and have some similarity with the thermal neutron fluence distributions produced by epithermal neutron beams [23]. Pencil beam algorithms are based on the principle that the dose resulting from an external irradiation field is equivalent to the dose distribution of an infinite number of infinitesimally small beams (pencil beams) determining the same field. Once the dose distribution of a single pencil beam has been theoretically or experimentally determined, the dose at a point in an arbitrary irradiation geometry can be calculated by a summation of the contribution of the individual pencil beams. The concept of pencil beams applied for electron beam dose calculation has been extensively described elsewhere [2,11]. The fluence distribution of an electron pencil beam is obtained from an analytical solution of the Boltzman equation under two different sets of simplifying assumptions. First, it is assumed that the electrons are scattered only over small angles, whereas the second assumption is based on quasi-isotropic diffusion. Both approaches, however, lead in first approximation to Gaussian shaped fluence distributions [2]. Using the second set of assumptions, the Boltzman equation is reduced to the age-diffusion equation which was developed to describe the diffusion of neutrons in a medium [13]. Therefore, the Gaussian pencil beam approach might also be a valuable tool for describing the thermal neutron distribution in a medium irradiated with an epithermal neutron beam. It is the purpose of this work to investigate the applicability to BNCT of the semi-empirical algorithms used for

C.P.J. Raaijmakers et al. / Radiotherapy and Oncology 46 (1998) 321–332

conventional treatment planning in order to develop a treatment planning method for BNCT that allows the application of interactive treatment planning procedures as used clinically for conventional treatment planning.

2. Materials and methods 2.1. Basic beam data Basic beam data, similar to those required for semiempirical photon and electron beam algorithms, were obtained for the clinical epithermal neutron beam (HB11) of the High Flux Reactor in Petten, The Netherlands. Thermal neutron fluence, gamma-ray dose and fast neutron dose distributions were measured in phantoms at various phantom to beam exit distances (PBDs). The reference data were acquired in a 30 × 30 × 30 cm3 water phantom positioned at a PBD of 30 cm for the three available field sizes. The field size can be varied by changing a 19 cm long cylindrical aperture. The apertures are conically shaped with exit diameters of 8.0, 12.0 and 15.0 cm. The field size will be referred to by the exit diameter of the aperture. The distance from the beam exit to the exit side of the aperture is approximately 10 cm. The influence of the PBD on the absolute value of the various dose components on the central beam axis of a phantom can be modelled using a virtual source located at a distance of 3.0 m from the beam exit. Details of the magnitude and shape of the fluence and dose distribu-

323

tions under reference conditions have been reported elsewhere [23]. 2.2. Phantoms In order to compare the treatment planning calculations with measurements in a clinically relevant situation, an anthropomorphic head phantom allowing insertion of various detectors at predefined positions was constructed. For this phantom, the CT contours of a head phantom often used in radiotherapy were used (Alderson radiation therapy phantom, Radiology Support Devices). The CT data (5 mm slice thickness) from the phantom were, if necessary, manually edited and the holes for detector positioning were indicated. The resulting data were translated by a conversion program into a format readable by a computer controlled milling machine and each slice was automatically milled out of a 5 mm thick polymethyl-methacrylate (PMMA) slab. Holes for detector positioning were milled and the slices were stacked on small fixation rods to form a head phantom (Fig. 1). A rod through the centre of the phantom is used to clamp the slices together. All parts of the phantom consist of PMMA. A plane of irradiation (central plane) was chosen 8 cm below the top of the skull, approximately through the middle of the brain. The most superficial measurement position in this plane was located at a depth of 2.7 cm. This position, where the maximum thermal neutron fluence is expected, was used as the normalization point in this work. On the

Fig. 1. Human head phantom based on CT contours of a standard anthropomorphic head phantom. The PMMA phantom allows the insertion of detectors at various positions and depths.

324

C.P.J. Raaijmakers et al. / Radiotherapy and Oncology 46 (1998) 321–332

central beam axis, five measurement positions were used, 2.5 cm distant from each other. Six positions at 5 cm and one position at 8.5 cm distance from the central beam axis were available. Since the treatment planning calculations were based on measurements in a water phantom, while the head phantom consisted of PMMA, the influence of phantom material on the dose and fluence distributions was investigated. This was done by comparing the various dose components measured in a 15 × 15 × 15 cm3 phantom of distilled water [23] with those in a solid PMMA phantom of the same dimensions. In the solid phantom, various positions along the central beam axis and two positions 5 cm off-axis at 2 cm depth are available.

HB11 beam. Thus, fluctuations in time of the beam output parameters were taken into account. A long-term reproducibility of approximately 1% (1 SD) is obtained in the thermal neutron fluence per monitor unit under output reference conditions. For the gamma-ray dose and the fast neutron dose this value amounts to approximately 3% (1 SD) [24]. When the difference between a measured and a calculated value is larger than the measurement uncertainty, it will be expressed as a percentage of the measured value at that point. In other words, when a relative value of 50% is calculated at a position where 40% is measured this will be expressed as a difference of 25%.

2.3. Phantom measurements

In the Gaussian pencil beam approximation the beam profiles of a pencil beam are Gaussian shaped and have a width that depends on the depth in the phantom. An empirical function is introduced to describe the depth dependence of the fluence along the central pencil beam axis. The fluence distribution, f(r,v,z), of a pencil beam perpendicularly incident on a flat surface of an infinite, homogeneous medium at a point with co-ordinates (r,v,z) = (0,0,0) can then be written as [4] f (z) f(r, v, z) = ∞ exp[ − r2 =r2 (z)Š (1) pr2 (z)

The phantoms, positioned at a reference PBD of 30 cm, were irradiated with a single horizontal field using the HB11 epithermal neutron beam [23]. The phantoms were irradiated with the central beam axis coinciding with the central measurement axis. The head phantom was positioned horizontally on a therapy table, on a support normally used for head positioning in the clinic. This irradiation position will probably also be used for the treatment of brain tumour patients with the BNCT facility at Petten in the future. Delineation of the phantoms was performed using four line lasers. Since the phantom size has a considerable influence on the dose and fluence values in the phantom [23], the absolute dose and fluence values at the normalization point in the head phantom can not be deduced directly from measurements in large cuboid phantoms, as is generally done in conventional radiotherapy. Therefore, these absolute values were compared with those at the same depth in the cuboid PMMA phantom that is approximately the same size. For an estimation of the dose and fluence values at the normalization point at 2.7 cm depth in the head phantom, the measured data in the cuboid PMMA phantom at 2 and 3 cm depth were interpolated using the shape of the depth dose curves determined in the water phantom. The thermal neutron fluence rate in the phantoms was determined using AuAl (5% wt. Au) and MnNi (88% wt. Mn) activation foils and a home-made silicon alpha-particle detector with a 6Li converter plate. The gamma-ray dose rate and fast neutron dose rate were determined using a magnesium ionization chamber flushed with argon gas and a tissue-equivalent ionization chamber flushed with tissue-equivalent gas. Details of the measurement procedures can be found elsewhere [21–23]. The uncertainty in the experimentally determined relative thermal neutron fluence and gamma-ray dose is estimated to be 2% (1 SD). The estimated uncertainty in the relative fast neutron dose amounts to 3% (1 SD) [21]. The beam output parameters were monitored with a monitor system installed in the collimator assembly of the

2.4. Gaussian pencil beams

where z is the depth in the phantom, r is the distance to the central beam axis, f∞ (z)=pr2 (z) is the central axis depth fluence curve of the pencil beam normalized in a convenient way as will be elucidated later and r2 (z) is the mean square radial spread of the particles at depth z. The square root of this parameter will be referred to as the pencil beam width. In age-diffusion theory the mean square radial spread is replaced by 4t, where the age-parameter t corresponds to the distance travelled by the neutrons to acquire a certain energy. For a collimated, uniform irradiation field of area A, the fluence distribution in a phantom is obtained by integrating the pencil beam distribution between the external field boundaries as follows … f∞ (z) exp‰ − (r − r′)2 =r2 (z)Šr′dr′dv (2) fA (r, v, z) = A pr2 (z) A purely empirical approach was chosen to determine the pencil beam width as a function of depth and the depth fluence function, f∞(z). In this method, previously applied to electron beams by Bruinvis et al. [3], only central beam axis data for two fields were used. For a circular field with radius R, the fluence on the central axis of the field, fR (z), can be analytically obtained from Eq. (2) yielding fR (z) = f∞ (z)(1 − exp‰ − R 2 =r2 (z)Š)

(3)

If R → ∞ we have fR (z) = f∞(z) and thus f∞(z) equals the central axis depth fluence curve for an infinitely large field. The ratio of the fluence values on the central beam axis for

C.P.J. Raaijmakers et al. / Radiotherapy and Oncology 46 (1998) 321–332

two circular fields with radii R 1 and R 2 is given by fR 1 (z) 1 − exp[ − R 21 =r2 (z)] = fR 2 (z) 1 − exp[ − R 22 =r2 (z)]

(4)

From measured fluence ratios, the pencil beam width was derived by numerically solving Eq. (4) for various depths. For the determination of the pencil beam width, the central axis depth fluence curves under reference conditions for the 8 and 15 cm beam apertures were chosen. The radii of the fields at the position of the phantom were determined, assuming that the neutrons originate from a point source. The distance from the phantom surface to the virtual source in combination with the exit diameter of the aperture then geometrically determines the field size. In this approximation, the diameter of the field at the surface of the phantom for a PBD of 30 cm amounts to 9.1, 13.6 and 17.1 cm for the 8, 12 and 15 cm diameter beam apertures, respectively. Once the pencil beam widths have been determined, the central axis depth fluence curve for a circular field of arbitrary size, including an infinitely large field, can be calculated using Eq. (4) and a measured depth fluence curve for a circular field. 2.5. Treatment planning system For the calculation of the relative distributions of the various dose components, a commercially available TPS was used (Plato version RTS 1.1.5 and 1.4, Nucletron, The Netherlands), which is designed for the calculation of external photon and electron beam dose distributions. The TPS runs on a Silicon Graphics Indigo machine with a 33 MHz R3000 processor. With the photon beam model an isodose pattern in the plane of calculation with a dose grid of 4 mm is generated within 5 s. With the electron beam model using a pencil beam grid of 2 mm an isodose pattern is generated in approximately 30 s. 2.5.1. Basic beam data For the thermal neutron component, the Gaussian pencil beam algorithm of the electron beam model was used. Parameters that describe the pencil beam width as a function of depth in the phantom and the depth fluence curve for an infinitely large field were entered manually. This set of data was completed by introducing the virtual source position. For the calculation of the gamma-ray dose and fast neutron dose, the photon beam model of the TPS (version 1.1.5) was used. The basic beam data for these components consisted of percentage depth dose (PDD) curves along the central beam axis, beam profiles at 2, 5 and 10 cm depth for every field size and the virtual source position. The basic beam data are independent of the actual calculation geometry and have to be introduced only once for a specific beam. 2.5.2. Calculations under reference conditions The thermal neutron fluence distributions under reference

325

conditions were calculated according to Eq. (2). The integration of the pencil beam contributions was approximated by a discrete summation over the area of the field. The field shape at the entrance of the phantom is introduced in the TPS by digitizing circular field shapes with the diameters mentioned in Section 2.4. Also the maximum fluence for a particular field size relative to that of an infinitely large field is obtained from the pencil beam summation. Calculations of the gamma-ray and fast neutron dose distribution under reference conditions consist merely of interpolations of the basic beam data. In the photon beam model, no algorithm to correct for the field shape was used since the dose distributions were explicitly measured for all available field shapes. 2.5.3. Calculations in phantoms of arbitrary shape Patient or phantom data can be introduced in the TPS using CT scans or manually digitized contours in multiple parallel planes. In both calculation models, the dose distribution is calculated per plane using a 2D approximation of the patient anatomy. In this approximation it is assumed that the patient contours in all planes are identical to that in the plane of calculation. In the pencil beam model, the phantom contour is taken into account by adjusting the point of incidence of each individual pencil beam. The intensity of a pencil beam is corrected as a function of the distance from the virtual source to the point of incidence according to the inverse square law. In the photon beam model, the TPS uses a geometrical method to account for the difference in contour shape between the patient and the reference phantom, known as the effective source-surface distance (SSD) method [12]. In this method, the relative dose at a point in a homogeneous phantom is assumed to be equal to the relative dose in a large cuboid phantom at the same off-axis distance and depth, when this phantom is positioned at the effective SSD. The effective SSD is defined as the distance between the phantom surface and a virtual source along the line from the point to the virtual source. The depth is defined as the distance between the point and the phantom surface along this line. The inverse square law is used to correct for the difference between the effective SSD and the reference SSD.

3. Results 3.1. Phantom composition The thermal neutron fluence rate at 2 cm depth on the central beam axis was 3% lower in the cuboid PMMA phantom than in the water phantom. No significant difference was observed between the shape of the central axis depth fluence curve for both phantoms (Fig. 2). The gamma-ray dose rate was 12% lower in the PMMA phantom compared with the water phantom while a somewhat higher relative

326

C.P.J. Raaijmakers et al. / Radiotherapy and Oncology 46 (1998) 321–332 Table 1 Measured fluence rate and dose rate values at a depth of 2.7 cm in a PMMA head phantom and a cuboid phantom of approximately the same size, irradiated with an epithermal BNCT beam Field size (cm) 8 12 15

Cuboid Head Cuboid Head Cuboid Head

fthermal (108 cm − 2 s − 1)

Dgamma (Gy h − 1)

Dfast (Gy h − 1)

5.02 4.89 6.98 7.01 7.83 7.70

1.74 1.64 2.31 2.25 2.51 2.50

0.53 0.53 0.58 0.55 0.58 0.54

(0.97) (1.00) (0.98)

(0.94) (0.97) (1.00)

(1.00) (0.95) (0.93)

The numbers in brackets are the fluence rate and dose rate values in the head phantom relative to those in the cuboid phantom.

Fig. 2. Boron neutron capture dose (X, solid line), gamma-ray dose (♦, broken line) and fast neutron dose (O, dotted line) on the central beam axis in a 15 × 15 × 15 cm3 water phantom (lines) and a PMMA phantom (symbols) of the same dimensions. The values have been normalized to 100% at the maximum value obtained in the water phantom.

gamma-ray dose at depth was observed in PMMA. The fast neutron dose rate at 2 cm depth in PMMA was 3% higher than in water. No significant difference between the fast neutron PDDs in both phantoms could be detected. The ratio of the dose and fluence values at 5 cm distance from the central beam axis and those on the central beam axis at a depth of 2 cm was the same in both phantoms for all components. The influence of the phantom composition on the dose and fluence distributions was independent of field size. 3.2. Pencil beam width The pencil beam width as a function of depth in the phantom could be described by a linear function and

Fig. 3. Measured central beam axis thermal neutron fluence curves in a 30 × 30 × 30 cm3 water phantom irradiated with an epithermal neutron beam for a 15 cm field (solid line), a 12 cm field (broken line) and an 8 cm field (dotted line). Also shown are the calculated fluence curves obtained using a Gaussian pencil beam model for a 12 cm field (O) and an infinitely large field (X). No normalization was performed.

increased from 4.5 cm at 2 cm depth to 7.0 cm at 10 cm depth. As a first test of the Gaussian pencil beam approximation, the depth fluence curve for the 12 cm field was calculated using Eq. (4). The shape of the calculated curve as well as the absolute fluence values were in excellent agreement with the measured values (Fig. 3). The depth fluence curve for an infinitely large field with complete lateral scattering equilibrium needed by the TPS to calculate the pencil beam depth fluence distribution differed slightly in shape from that of the 15 cm field. The maximum thermal neutron fluence was approximately 4% higher for the infinitely large field compared with the 15 cm field (Fig. 3). 3.3. Fluence and dose value measurements In general, the absolute dose and fluence values at the normalization point in the head phantom agreed with those in the cuboid phantom within the measurement uncertainty. The values in the head phantom are in general slightly lower. Only for the gamma-ray dose using the 8

Fig. 4. Comparison of the measured (symbols) and calculated (lines) thermal neutron profiles at 2 cm (X), 5 cm (O) and 10 cm (♦) depth in a large cuboid water phantom irradiated with the 12 cm diameter epithermal BNCT beam. The calculations have been performed by a conventional treatment planning system using empirical data obtained on the central beam axis for 8 and 15 cm fields.

C.P.J. Raaijmakers et al. / Radiotherapy and Oncology 46 (1998) 321–332

cm field does the difference (6%) appear to be significant (Table 1). The 7% difference between the fast neutron dose in the head phantom and the cuboid phantom for the 15 cm field might be explained from the relatively low accuracy of the measurements for this dose component. 3.4. Treatment planning calculations The relative thermal neutron fluence distributions under

327

reference conditions, calculated by the TPS, were in good agreement with the measured distributions at all phantom positions for all field sizes. The field size dependence of the absolute thermal neutron fluence was exactly predicted by the TPS. It should be noted that for the calculation of beam profiles only central beam axis data for the 8 and 15 cm fields were used (Fig. 4). In the head phantom, the agreement between measurements and calculations was within the measurement uncer-

Fig. 5. Thermal neutron isofluence lines in a human head phantom calculated using the Gaussian pencil beam algorithm compared with measurements in the phantom for an 8 cm (a), 12 cm (b) and 15 cm (c) diameter epithermal BNCT beam. The calculated position of the maximum fluence is indicated (M). The dotted lines indicate the field size and the beam direction.

328

C.P.J. Raaijmakers et al. / Radiotherapy and Oncology 46 (1998) 321–332

tainty of approximately 2% (1 SD) along the central beam axis for all field sizes (Fig. 5). The deeper penetration of thermal neutrons in the head phantom compared with the cuboid phantom, caused by the curved surface of the head phantom, was exactly predicted by the TPS. For the 15 cm field, where this effect is most prominent, a relative thermal neutron fluence of 83% was measured at 5.2 cm depth while a value of 84% was calculated (Fig. 5c). In the large cuboid water phantom a value of 73% was obtained at this depth (Fig. 3). Good agreement between measurements and calculations was obtained at most off-axis positions. Near the lateral phantom boundary, however, an overestimation of the thermal neutron fluence by the TPS of approximately 35% was observed. At the most distant position from the central beam axis, at the 10–20% level, the TPS also overestimated the fluence by approximately 35%. The difference between the measured and the calculated relative gamma-ray dose distribution in the head phantom was at most positions smaller than 5% for the various field sizes. A small underestimation of the gamma-ray dose component on the central beam axis was observed (Fig. 6). Larger differences between measurements and calculations were observed at positions close to the lateral phantom boundary. At these positions the TPS overestimates the gamma-ray dose by approximately 10–20%. At the most

Fig. 7. Fast neutron isodose lines in a human head phantom irradiated with a 12 cm diameter epithermal BNCT beam calculated using a conventional photon beam model. Measured values at specific points are also given. The dotted line indicates the field size and the beam direction.

distant position from the central beam axis, an overestimation of the gamma-ray dose by the TPS of approximately 30% was observed. No significant difference was observed between the measured and the calculated fast neutron dose distribution (Fig. 7). An occasional difference up to approximately 20% can be attributed to the relatively low reproducibility of the measurements of this dose component especially at depths larger than 5 cm.

4. Discussion

Fig. 6. Gamma-ray isodose lines in a human head phantom irradiated with a 12 cm diameter epithermal BNCT beam calculated using a conventional photon beam model. The calculated position of the maximum dose is indicated (M). Measured values at specific points are also given. The dotted line indicates the field size and the beam direction.

A treatment planning method for BNCT has been developed that uses a commercially available TPS. The system is capable of calculating the distribution of the various dose components in a phantom or a patient irradiated with an epithermal neutron beam within 1 min. The short calculation times and the various viewing options of such a commercial system allow the application of the same treatment planning procedures for BNCT as used clinically for conventional treatment planning. The agreement between measurements and calculations in the central plane of a homogeneous anthropomorphic head phantom was within the measurement uncertainty at most positions. Significant differences between measurements and calculations were

C.P.J. Raaijmakers et al. / Radiotherapy and Oncology 46 (1998) 321–332

only observed near the phantom boundary, generally in low dose regions. Consequently, treatment planning methods as generally applied in conventional radiotherapy are suitable for BNCT. The type of algorithms are wide-spread and thus treatment planning for BNCT can probably be performed in most radiotherapy centres, provided a set of beam data under reference conditions is available. At positions close to the lateral phantom boundary or at relatively large distances from the central beam axis, the TPS overestimates the thermal neutron fluence and gamma-ray dose by approximately 30%. In many cases this might be acceptable under the condition that this limitation is recognized, since the overestimation occurs mainly in low dose regions. However, the presence of organs at risk at such positions, such as the eyes, optic chiasm or brain stem, might require a higher accuracy. The various dose components at these positions could be accurately determined by a simulation of the actual irradiation using a head phantom of approximately the same size as that of the patient with measurement positions at the organs at risk. A limited set of phantoms will probably be sufficient to cover most head sizes. The Gaussian solution of the age-diffusion equation does not directly describe the overall distribution of thermal neutrons but describes the distribution of neutrons that just have acquired thermal energy. After slowing down to thermal energies, the neutrons become distributed according to diffusion laws. Theoretically, the diffusion kernel function should, therefore, be applied over the Gaussian pencil beam distribution to represent the thermal neutron distribution. Only minor differences were, however, observed between measurements and calculations under reference conditions (Fig. 4), indicating that the Gaussian pencil beam approximation is valid for the calculation of thermal neutron distributions. The overestimation of the thermal neutron fluence near the phantom boundary can qualitatively be explained from ‘missing tissue’ effects. Since neutrons leak out of the phantom near the edges while no neutrons are scattered back into the phantom, the thermal neutron fluence at these positions is considerably lower than the fluence at positions at the same depth and off-axis distance in an infinitely large phantom. This effect is not taken into account by the pencil beam model which has included the assumption that the shape of the pencil beam fluence distribution is only dependent on the material traversed by the central axis of the pencil beam. Furthermore, the TPS assumes that the phantom contours are identical in all planes. Because of the large neutron pencil beam widths, this is a crude approximation in a spheroid object like the head. Since several commercially available TPSs use algorithms similar to those used in this work, it can be expected that the calculations can be further extended to TPSs which take the phantom contour into account in three dimensions. A considerable missing tissue effect was, however, already observed for thermal neutron distributions measured in cuboid phantoms of various size

329

[23]. This indicates that the observed overestimation of thermal neutron fluence near the phantom boundary results probably mainly from the limitations of the applied algorithm, rather than from the approximation of the phantom shape. Another algorithm will be needed to improve the performance of the TPS at these positions. Pencil beam models as applied to electron beam treatment planning have known limitations especially in dose calculations near inhomogeneities at depths where the angular distribution of the electrons is large [14]. These limitations are essentially due to the same phenomenon as the missing tissue effect, i.e. the approximation that the pencil beam fluence distribution is only influenced by the material traversed by the central axis of the pencil beam. The width of the neutron pencil beams is considerably larger than the width of electron pencil beams which ranges from approximately 0.2 to 2 cm, depending on the depth in the phantom and the energy of the electrons. Consequently, the volume of surrounding tissue influencing the fluence at a point is much larger for thermal neutrons than for electrons and the limitations of the pencil beam approach are observed more clearly when these algorithms are applied to epithermal neutron beams. Improved calculation models under development for electron beam treatment planning [10] might be useful for BNCT treatment planning calculations. It has been proposed to apply the principles of BNCT to fast neutron therapy [20]. In tissue irradiated with fast neutron beams, a non-negligible thermal neutron fluence is created which, in combination with a high boron concentration in the tumour, could result in a substantial BNCT boost dose to the tumour treated with fast neutron therapy. As demonstrated by Matzen et al. [16], the Gaussian pencil beam model is also well suited for the calculation of thermal neutron distributions for this form of radiotherapy. For the calculation of the gamma-ray dose distribution, use has been made of a simple geometrical calculation model. Nevertheless, the differences between measurements and calculations were within 5% at most phantom positions. The overestimation of the gamma-ray dose near the phantom boundaries can be explained by missing tissue effects. The effective SSD method is physically not completely correct since the induced gamma rays are generated in the phantom and do not originate from a distant virtual point source. However, since the reference conditions are rather similar to the actual irradiation geometry, the corrections are relatively small and a good estimation of the gamma-ray dose distribution is obtained using this method. This example illustrates the inherent advantages of semiempirical treatment planning methods over other calculation techniques. A logical improvement of the calculation model for the gamma-ray dose component would be the application of a photon pencil beam model, as is implemented in the latest version of the Nucletron TPS and which is currently being tested. A physically more correct approach to calculate the

330

C.P.J. Raaijmakers et al. / Radiotherapy and Oncology 46 (1998) 321–332

gamma-ray dose distribution uses the thermal neutron distribution as a source for capture gamma-ray production. By applying an appropriate kernel function over this distribution, the capture gamma-ray dose distribution can be calculated as reported by Gupta et al. [8]. In such an approach the incident gamma-ray dose must be calculated separately. Empirical data on the distribution of the separate incident gamma-ray dose in a hydrogen-rich phantom are, however, difficult to obtain due to the presence of the capture gammarays. For the fast neutron component, scattering effects are of little importance. The simple contour correction method as applied for external photon beams can, therefore, well be applied and further improvement of the treatment planning calculations at the moment seems unnecessary. It should be noted that the fast neutron dose component reaches its maximum value at a very superficial depth where no measurement position was available in the head phantom. At 2.7 cm depth the fast neutron dose amounts to only approximately 50% of its value at 1 cm depth (Fig. 2). Measurements of the fast neutron dose at superficial depths in a head phantom might be needed to validate the treatment planning calculations in the high dose region. A major difference between the Petten epithermal neutron beam and other existing and proposed beams is the large distance between the beam filter and the beam exit resulting in a high degree of collimation of the Petten beam. For poorly collimated beams, the virtual source approximation might be less valid which might necessitate adjustment of the algorithms when applied to other epithermal neutron beams. It can, however, be expected that also for such beams semi-empirical algorithms can be applied for treatment planning purposes by adjusting the position and extension of the virtual source(s). Obviously, the basic beam data and the pencil beam parameters have to be determined experimentally for each specific beam. Only minor differences were observed between the relative distributions of the various components in the cuboid PMMA and water phantoms. Therefore, relative measurements in the PMMA head phantom were directly compared with the calculations based on empirical data obtained in a water phantom. At this moment, the influence of the elemental composition and density of the phantom material and of tissue inhomogeneities on the distribution of the various dose components can not be taken into account in the calculations. The photon and electron beam model of the TPS have algorithms implemented for this purpose, which may be adapted for epithermal neutron beams. The use of these algorithms implies empirical knowledge on the effect of phantom composition on the dose distributions. These effects are currently under investigation. Theoretical calculations in a canine head have indicated that effects of macroscopic tissue heterogeneities are generally small but nonnegligible, in particular not near void regions [17]. Since the human brain has a reasonably uniform composition, it can be expected that semi-empirical correction algorithms can

account for tissue composition and inhomogeneities effects to an acceptable level of accuracy. A conventional TPS is an excellent tool for optimizing the irradiation geometry in order to arrive at the best treatment plan for an individual patient. Once the optimal geometry has been determined, time consuming Monte-Carlo calculations could be performed for a final calculation. In these TPSs, inhomogeneities can be taken into account in a full 3D calculation [19,28]. Currently, however, no experimental results have been reported on the accuracy nor on the necessity of such calculations. For an estimation of the expected healthy tissue and tumour response, the various dose components have to be weighted and then added. This can be performed by the TPS using its conventional beam weighting options. Obviously, the weighting factors at the normalization point have to be calculated separately using the absolute dose values and the appropriate RBE values. Dose prescription for BNCT is a complicated issue. The relative contribution of the various dose components to the total biologically weighted dose varies largely and depends, among others, on the position of the prescription point, the irradiation geometry, the boron concentration in the various tissues, the boron compound, the chosen RBE values and the fractionation scheme. Gupta et al. [9] reported an investigation of the effect of the different treatment variables on the biologically weighted dose values using fitted dose curves on the central beam axis in a cylindrical phantom. Due to the short calculation times and high accuracy near the central beam axis, the TPS presented in this work is an excellent tool for the investigation of such parameters in clinically relevant situations and will be useful in visualizing some of these difficult aspects of dose prescription for BNCT. Zamenhof et al. [28] reported a comparison between calculations using a Monte-Carlo based TPS and measurements on the central beam axis in a reference phantom irradiated with the epithermal neutron beam of the Massachusetts Institute of Technology for one irradiation geometry. For 6 million particle histories, requiring 60 h of calculation time, the statistical uncertainty (1 SD) was estimated to be 4, 30 and 12% for the thermal neutron fluence, the gamma-ray dose and the fast neutron dose, respectively. The calculation time could be reduced significantly if faster computers were used. Differences between measurements and calculations of the relative thermal neutron fluence of approximately 10% at the 80% level were explained from the statistical calculation uncertainty and a relatively large measurement uncertainty of 7% (1 SD). Similar differences were reported at the 40% level for the relative fast neutron dose. Part of the incident gamma-ray dose was estimated. No other experimental validations have been reported for TPSs for BNCT. In this work, the differences between relative measurements and calculations on the central beam axis were within the measurement uncertainty of approximately 2–3% (1 SD), which is smaller than the statistical uncer-

C.P.J. Raaijmakers et al. / Radiotherapy and Oncology 46 (1998) 321–332

tainty reported for the Monte-Carlo calculations. It can, therefore, be concluded that besides the obvious advantages of semi-empirical treatment planning calculations, regarding calculation times and availability, also the accuracy on the central beam axis in a homogeneous phantom is similar to the accuracy currently obtained using a Monte-Carlo based TPS [28] for BNCT.

5. Conclusions A fast and accurate treatment planning method has been developed for BNCT using an epithermal neutron beam. The method uses conventional dose calculation algorithms implemented in a commercially available TPS. The agreement between measurements and calculations in the central plane of a homogeneous anthropomorphic head phantom was excellent both on the central beam axis as well as at most off-axis positions. Only at positions near the phantom boundaries or relatively far away from the central beam axis were local differences of approximately 30% observed. The system is the first treatment planning system for BNCT that allows interactive treatment planning and will be a useful tool for the optimization of treatment plans for patients treated with BNCT.

[7] [8]

[9]

[10]

[11]

[12] [13] [14]

[15]

[16]

[17]

Acknowledgements The authors would like to thank T. Minderhoud for assistance with the treatment planning system, L. de Mooij and E. Kos for constructing the head phantom, W. Voorbraak and A. Paardekoper for preparing and analyzing the activation foils, K. Ravensberg and F. Stecher-Rasmussen for technical assistance and for operating the epithermal neutron beam and R. Moss for providing beam time. This work was financially supported by the Dutch Cancer Society (NKB Grant NKI 94-778).

[18]

[19]

[20]

[21]

References [22] [1] Barth, R.F., Soloway, A.H. and Brugger, R.M. Boron neutron capture therapy of brain tumors: past history, current status, and future potential. Cancer Invest. 146: 534–550, 1996. [2] Brahme, A. Current algorithms for computed electron beam dose planning. Radiother. Oncol. 3: 347–362, 1985. [3] Bruinvis, I.A.D., Van Amstel, A., Elevelt, A.J. and Van der Laarse, R.V. Calculation of electron beam dose distributions for arbitrarily shaped fields. Phys. Med. Biol. 28: 667–683, 1983. [4] Bruinvis, I.A.D., van Amstel, A., Elevelt, A.J. and van der Laarse, R.V. Dose calculation for arbitrarily shaped electron beams. Acta Radiol. 364 (Suppl.): 11–19, 1983. [5] Coderre, J.A., Bergland, R. and Chadha, M. et al. Boron neutron capture therapy of glioblastoma multiforme using the p-boronophenylalanine-fructose complex and epithermal neutrons. In: Cancer Neutron Capture Therapy, pp. 553–562. Editor: Y. Mishima. Plenum Press, New York, 1996. [6] Dorn III, R.V. Boron neutron capture therapy (BNCT): a radiation

[23]

[24]

[25]

[26]

331

oncology perspective. Int. J. Radiat. Oncol. Biol. Phys. 28: 1189– 1201, 1994. Gabel, D. Present status and perspectives of boron neutron capture therapy. Radiother. Oncol. 30: 199–205, 1994. Gupta, N., Niemkiewicz, J., Blue, T.E., Gahbauer, R. and Qu, T.X. Effect of head phantom size on 10B and 1H(n,g)2H dose distributions for a broad field accelerator epithermal neutron source for BNCT. Med. Phys. 20: 395–404, 1993. Gupta, N., Gahbauer, R.A., Blue, T.E. and Wambersie, A. Dose prescription in boron neutron capture therapy. Int. J. Radiat. Oncol. Biol. Phys. 28: 1157–1166, 1994. Huizenga, H. and Storchi, P.R.M. Numerical calculation of energy deposition by broad high energy electron beams. Phys. Med. Biol. 34: 1371–1396, 1989. Jette, D. Electron beam dose calculations. In: Medical RadiologyRadiation Therapy Physics, pp. 95–120. Editor: A.R. Smith. Springer-Verlag, Berlin, 1995. Khan, F.M. The Physics of Radiation Therapy, 2nd edn. Williams and Wilkins, Baltimore, MD, 1994. Lamarch, J.R. Introduction to Nuclear Reactor Theory. AddisonWesley, Reading, MA, 1966. Lax, I. Inhomogeneity corrections in electron-beam dose planning. Limitations with the semi-infinite slab approximation. Phys. Med. Biol. 31: 879–892, 1986. Liu, H.B., Brugger, R.M., Greenberg, D.D., Rorer, D.C., Hu, J. and Hauptman, H.M. Enhancement of the epithermal neutron beam used for boron neutron capture therapy. Int. J. Radiat. Oncol. Biol. Phys. 28: 1149–1156, 1994. Matzen, T., Lu¨demann, L., Schmidt, R. and Scobel, W. Calculation of the thermal neutron flux in a 14.1 MeV neutron beam facility for application in BNCT. Z. Med. Phys. 6: 170–175, 1996. Moran, J.M., Nigg, D.W., Wheeler, F.J. and Bauer, W.F. Macroscopic geometric heterogeneity effects in radiation dose distribution analysis for boron neutron capture therapy. Med. Phys. 19: 723–732, 1992. Morris, G.M., Coderre, J.A., Hopewell, J.W., Micca, P.L. and Rezvani, M. Response of rat skin to boron neutron capture therapy with p-boronophenylalanine or borocaptate sodium. Radiother. Oncol. 32: 144–153, 1994. Nigg, D.W. Methods for radiation dose distribution analysis and treatment planning in boron neutron capture therapy. Int. J. Radiat. Oncol. Biol. Phys. 28: 1121–1134, 1994. Poller, F., Sauerwein, W. and Rassow, J. Monte Carlo calculation of dose enhancement by neutron capture of 10B in fast neutron therapy. Phys. Med. Biol. 38: 397–410, 1993. Raaijmakers, C.P.J., Watkins, P.R.D., Nottelman, E.L., Verhagen, H.W., Jansen, J.T.M., Zoetelief, J. and Mijnheer, B.J. The neutron sensitivity of dosimeters applied to boron neutron capture therapy. Med. Phys. 23: 1581–1589, 1996. Raaijmakers, C.P.J., Konijnenberg, M.W., Verhagen, H. and Mijnheer, B.J. Determination of dose components in phantoms irradiated with an epithermal neutron beam for boron neutron capture therapy. Med. Phys. 22: 321–329, 1995. Raaijmakers, C.P.J., Konijnenberg, M.W. and Mijnheer, B.J. Clinical dosimetry of an epithermal neutron beam for neutron capture therapy; dose distributions under reference conditions. Int. J. Radiat. Oncol. Biol. Phys. 37: 941–951, 1997. Raaijmakers, C.P.J., Nottelman, E.L., Konijnenberg, M.W. and Mijnheer, B.J. Dose monitoring for boron neutron capture therapy using a reactor based epithermal neutron beam. Phys. Med. Biol. 41: 2789–2797, 1996. Reinstein, L.E., Ramsay, E.B., Gajewski, J., Ramamoorthy, S. and Meek, A.G. SBNCT: a 3-dimensional treatment planning system for boron neutron capture therapy. In: Advances in Neutron Capture Therapy, pp. 171–175. Editors: A. Soloway, R.F. Barth and D.E. Carpenter. Plenum Press, New York, 1993. Rogus, R.D., Harling, O.K. and Yanch, Y.C. Mixed field dosimetry

332

C.P.J. Raaijmakers et al. / Radiotherapy and Oncology 46 (1998) 321–332

of epithermal neutron beams for boron neutron capture therapy at the MITR-II research reactor. Med. Phys. 21: 1611–1625, 1994. [27] Wheeler, F.J., Wessol, D.E., Babcock, R., Nigg, D.W., Atkinson, C.A. and Evans, J. Improvements in patient treatment planning systems. In: Cancer Neutron Capture Therapy, pp. 289–294. Editor: Y. Mishima. Plenum Press, New York, 1996.

[28] Zamenhof, R., Redmond, E., Solares, R. et al. Monte Carlo based treatment planning for boron neutron capture therapy using custom designed models automatically generated from CT data. Int. J. Radiat. Oncol. Biol. Phys. 35: 383–397, 1996.