INSTITUTE OF PHYSICS PUBLISHING Phys. Med. Biol. 51 (2006) 5785–5807 PHYSICS IN MEDICINE AND BIOLOGY doi:10.1088/0031-9155/51/22/005 Comparison of dose calculation algorithms for treatment planning in external photon beam therapy for clinical situations Tommy Knöös1, Elinore Wieslander1, Luca Cozzi2, Carsten Brink3, Antonella Fogliata2, Dirk Albers4, Håkan Nyström5 and Søren Lassen5 1 Radiation Physics, Lund University Hospital, S-221 85 Lund, Sweden Medical Physics, Oncology Institute of Southern Switzerland, Bellinzona, Switzerland 3 Laboratory of radiation Physics, Odense University Hospital, DK-5000 Odense C, Denmark 4 Department of Radiotherapy and Radio-Oncology, University Medical Center Hamburg-Eppendorf, Germany 5 Department of Radiation Physics, 3994, University Hospital of Copenhagen, Rigshospitalet, DK-2100, Copenhagen, Denmark 2 E-mail: tommy.knoos@med.lu.se Received 1 May 2006, in final form 22 September 2006 Published 24 October 2006 Online at stacks.iop.org/PMB/51/5785 Abstract A study of the performance of five commercial radiotherapy treatment planning systems (TPSs) for common treatment sites regarding their ability to model heterogeneities and scattered photons has been performed. The comparison was based on CT information for prostate, head and neck, breast and lung cancer cases. The TPSs were installed locally at different institutions and commissioned for clinical use based on local procedures. For the evaluation, beam qualities as identical as possible were used: low energy (6 MV) and high energy (15 or 18 MV) x-rays. All relevant anatomical structures were outlined and simple treatment plans were set up. Images, structures and plans were exported, anonymized and distributed to the participating institutions using the DICOM protocol. The plans were then re-calculated locally and exported back for evaluation. The TPSs cover dose calculation techniques from correction-based equivalent path length algorithms to modelbased algorithms. These were divided into two groups based on how changes in electron transport are accounted for ((a) not considered and (b) considered). Increasing the complexity from the relatively homogeneous pelvic region to the very inhomogeneous lung region resulted in less accurate dose distributions. Improvements in the calculated dose have been shown when models consider volume scatter and changes in electron transport, especially when the extension of the irradiated volume was limited and when low densities were present in or adjacent to the fields. A Monte Carlo calculated algorithm input data set and a benchmark set for a virtual linear accelerator have been produced which have facilitated the analysis and interpretation of the results. The more sophisticated 0031-9155/06/225785+23$30.00 © 2006 IOP Publishing Ltd Printed in the UK 5785 5786 T Knöös et al models in the type b group exhibit changes in both absorbed dose and its distribution which are congruent with the simulations performed by Monte Carlo-based virtual accelerator. S This article has associated online supplementary data files (Some figures in this article are in colour only in the electronic version) 1. Introduction The evaluation of dose calculation algorithms for radiotherapy treatment planning systems (TPSs) in external photon beam therapy is usually based on a number of measurements performed in simple geometries such as slabs of various density and/or composition. According to the technical report by IAEA, these actions of verification can be divided into benchmark data, generic beam data and the user’s beam data (IAEA 2005). From a logistic point of view the validation procedure starts at the vendor and continues with an increasing transfer of responsibility to the local user until the final decision of using the system clinically. A large number of papers have been published in this field (see, e.g., AAPM Report 85 (2004) for a list of references) with only limited studies with several TPSs applied in clinical situations. For example, Cheng et al (1997) presented a study with 11 different systems calculating a tangential treatment for conservative breast cancer. In the present study, the full implementation of a TPS, for some of the most commonly used and commercially available systems, was compared using identical patient data sets based on CT. For all systems, exactly the same treatment plans were applied. The main goal of this study is focused on the handling of electron transport especially in non-water conditions and how the systems manage to consider volume (phantom) scatter integration. During the last 40 years a number of papers have been published on characterizing and modelling in volume scatter in the irradiated medium starting with the paper by Cunningham (1972). Cunningham introduced the concept of scatter–air ratios in combination with the Clarkson sector integration (Clarkson 1944). The widely known program code for scatter integration, IRREG, was based on this concept (Cunningham et al 1972). Other examples include the equivalent tissue–air ratio method (ETAR) which is based on the O’Connor rectilinear scaling theorem (O’Connor 1957). This correction algorithm has been the most sophisticated approach to model phantom scatter for many years (Sontag and Cunningham 1978a, 1978b). The latest developments using pencil and point kernels for scatter integration have been reviewed together with the earlier methods by Ahnesjö and Aspradakis (1999) and also in the AAPM Report 85 (2004). The limitations of applying the O’Connor theorem for primary electron transport in photon beams have been shown by, e.g., Woo and Cunningham (1990). Several methods have been proposed to explicitly account for changes in electron transport between media, but to date none of these have been implemented in a commercial TPS; mainly due to the prohibitive increases in computation time (Keall and Hoban 1995, Yu et al 1995, Keall and Hoban 1996). All current commercial dose calculation algorithms (other than Monte Carlo approaches) use rectilinear density scaling to approximate dose contributions originating from changes in primary electron transport. The most advanced models are based on collapsed-cone convolution techniques with a separation of primary and scatter dose point kernels. These kernels are described approximately with a number (∼100) of lines along which the energy/dose is transported. A Dose calculation algorithms for treatment planning in external photon beam therapy 5787 rectilinear scaling following the O’Connor theorem is performed separately for the two point kernels along these lines with different ‘attenuation’ (Miften et al 2000, Ahnesjö et al 2005). The purpose of this work is to study how well these approximations for electron transport as well as scatter models for photons perform using clinical treatment planning systems of today. Four patient cases have been selected, representing different degrees of complexity regarding heterogeneity. The patient cases were evaluated with five commercial systems representing nine different implementations and/or calculation algorithms. As a golden standard or benchmark for the evaluation, the same patient/irradiation geometry was also set up in a Monte Carlo simulation environment. 2. Material and methods The basic requirement for including a commercial TPS in this study has been either the availability of the DICOM protocol with its extensions for radiotherapy, including RTPLAN, RT-STRUCTURES, RT-IMAGES and RT-DOSE6 , or the RTOG7 protocol. With these protocols it was possible to export/import all necessary data to perform the planning and calculations as well as to export the resulting data, i.e., dose matrices for analysis. Rather simple beam arrangements have been applied to minimize problems porting the plans to the included systems including the Monte Carlo environment. Thus, only open rectangular beams without wedges (physical or software based) or shaped beams (blocks or MLC) have been used. This restriction to simple open beams limits the influence from different scattering models present in the TPSs handling wedge filters, blocks, collimator jaws and MLCs. 2.1. Patient data All patient data were acquired at the same CT scanner8 and the information was transferred via DICOM to a TPS. The patient cases were selected from the clinical database at one of the participating centres. All studies were performed with abutted 0.5 cm thick slices and all identification tags were removed to avoid future reference to the original patient. A sufficient number of structures including the patient outline were defined for all cases. These structures were similar to the ones used for clinical planning based on the treatment policy at one of the institutions. The included cases represent a large proportion of the most common treatments at a radiotherapy department. The following cases were chosen: (a) Prostate: a conventional four-field box technique with equal weight on all four beams (field size 8.6 × 7.6 cm2). Gantry directions: 0◦ , 90◦ , 180◦ and 270◦ . (All geometric parameters are according to IEC (1996)). (b) Head and neck: two lateral and equally weighted fields covering the primary tumour and the bi-lateral nodes (field size 16 × 18 cm2). Gantry directions: 90◦ and 270◦ . (c) Breast: two equally weighted, almost opposed beams with the medial field border aligned to avoid unnecessary irradiation of the underlying lung tissue (field size 9 × 17 cm2). Gantry directions: 309◦ and 134◦ . (d) Lung: a five-field technique with different weights on all beams to treat a boost volume in the left lung (field size 11 × 10 cm2). Gantry directions: 0◦ , 72◦ , 144◦ , 216◦ and 288◦ . 6 7 8 See http://medical.nema.org/ for information regarding DICOM. See http://rtog3dqa.wustl.edu/ for information regarding the RTOG file format. General Electric HiSpeed NX/i Pro at Lund University Hospital, Sweden. 5788 T Knöös et al All beams are set up isocentric and the beam weight is equal to the relative number of monitor units (m.u.), i.e. the same beam weight, e.g., for a four-field box technique means that the relative number of m.u. is equal. 2.2. Treatment planning systems The included TPSs in this study are from five of the leading vendors in the radiotherapy community, i.e., Nucletron: Oncentra MasterPlan (OMP), Philips: Pinnacle (Pinnacle), two systems from Varian: Eclipse (Eclipse 1 and 2), Elekta: PrecisePlan (PPLAN) and CMS: XiO (XiO).9 These systems are based on a broad spectrum of algorithms for volume (phantom) scatter integration and heterogeneity modelling. The OMP system has two different models: both are descendents from the older Helax-TMS system—pencil beam convolution/superposition (OMP/PB) and collapsed-cone convolution (OMP/CC). The models are based on energy fluence and also include head scatter modelling. The first model is based on a two-dimensional pencil beam convolution for volume integration. Inhomogeneities are handled by an equivalent path length correction (EPL) for the primary dose contribution and a one-dimensional convolution along fan lines with an exponential for scattered radiation (Ahnesjö and Trepp 1991, Ahnesjö et al 1992, 2005). The pencil beam is parameterized at every 0.25 cm along the direction of propagation down to 50.25 cm. The second model in OMP is a collapsed-cone convolution approach where a ray-trace procedure through the irradiated object is utilized to get the TERMA (Total Energy Released per unit MAss, Ahnesjö 1989) at all points in the dose calculation matrix. The TERMA is separated into a primary part (collision kerma) and a scatter part (i.e., the difference between TERMA and collision kerma which is named scerma) which each are transported separately along 106 lines from the interaction point. The energy from each voxel intersected by a fan line in the irradiated medium is collected and deposited according to the elemental composition of the medium and density variations along the fan line (Ahnesjö 1989, Ahnesjö et al 2005). In this way, the collapsed-cone convolution as implemented in OMP actually calculates dose to tissue which was shown in an earlier paper (Wieslander and Knöös 2003). Both the pencil beam and the collapsed cone are based on the same point dose kernels produced by the EGS4 Monte Carlo code (Mackie et al 1988). These two dose calculation models have been extensively benchmarked by several scientific groups (e.g. Hurkmans et al 1995, Knöös et al 1995, van’t Velt 1997, Aspradakis et al 2003, Nisbet et al 2004, Haedinger et al 2005, Paelinck et al 2005), thus the performance and limitations are fairly well known. The system from Philips (Pinnacle/CC) is also based on a fluence model using collapsedcone convolution but this model is developed by the Madison group in Wisconsin. In fact, the development of this model and the one in the OMP occurred in parallel, thus many similarities exist, e.g. the point dose kernels (Mackie et al 1988), but different approaches in the solution of the problem have been used. For example, the kernel is not separated in a primary and a scatter part during convolution as in the other two collapsed-cone models studied in this work. Details for the collapse-cone convolution algorithm in Pinnacle are given in, e.g., Papanikolaou et al (1993), McNutt et al (1996) and Starkschall et al (2000). To facilitate computation speed, faster implementations of this model are also available in Pinnacle. These are the ‘fast convolve’ and ‘adaptive convolution’ methods. In the present study, only the collapsed-cone convolution was included in this study since it is the only one used for curative patients at the participating institution. 9 Nucletron Oncentra MasterPlan (OMP) V1.4.1.2, Philips Pinnacle3 V7.4f, Varian Eclipse system 1—V7.3.10 with algorithm V7.1.63 and system 2—V7.3.10 with algorithm V7.5.14.3, Elekta PrecisePlan Release 2.02, Computerized Medical Systems XiO Release 4.2.0. Dose calculation algorithms for treatment planning in external photon beam therapy 5789 For the systems from Varian (Eclipse 1 and 2) three different models or algorithms are studied. Two of the algorithms have been developed as an extension to a previously developed algorithm for rectangular fields which is based on the Milan–Bentley method. A pencil beam convolution/superposition is added to consider irregular shaped beams. The same pencil beam integration method and data are used for the first two models. This pencil is parameterized at five depths and accordingly the convolution with the incoming fluence is only performed at these depths; for depths in between, interpolation is performed. The pencil beams are extracted from measurements in a similar way to in the scatter formalisms, e.g., TPR/SPR approaches. Details are given in Storchi and Woudstra (1995, 1996) and Storchi et al (1999). This calculation model is named single pencil beam (SPB). For the SPB-based integration model two heterogeneity models have been analysed—the first one is the modified Batho (Webb and Fox 1980). This method only considers density variations along the fan line, i.e. it is an EPL-based method. A more sophisticated method which accounts for the density distribution in three dimensions is the second method, i.e., the ETAR method (Sontag and Cunningham 1978b). The ETAR implementation follows, unfortunately, the original papers where the full volume is approximated by a single scattering slice/plane to increase calculation efficiency instead of performing the calculation in the full three-dimensional volume (Sontag and Cunningham 1978a, 1978b). None of these models calculate the dose based on variations in density at positions lateral to the calculation point. These models will be abbreviated Eclipse/ModBatho and Eclipse/ETAR. The third model, named anisotropic analytical algorithm (AAA), present in the second Eclipse system, is also based on a pencil beam convolution technique. The pencil beam is compiled from Monte Carlo calculations and adjusted to fit measurements. In this case, the longitudinal distribution of the pencil beam is scaled according to EPL. The lateral extension of the pencil beam is, unlike the other two algorithms implemented in Eclipse, scaled with the density relative to water in directions normal to the pencil beam. The lateral description of the pencil beam is a sum of three terms where the first represents short-range dose depositions which approximate the primary electrons set in motion along the pencil beam and the other two scattered photons. When scaling, the densities in the previous layer are slightly considered when the lateral distribution is adjusted according to the radiological distance to the calculation point. In this way, the changes in lateral transport of energy are modelled when the density varies in the irradiated object (Ulmer and Kaissl 2003, Ulmer et al 2005). The next system is the Elekta PrecisePlan (PPLAN) which was earlier known as RenderPlan 3D and was developed primarily by W D Renner in the late 1980s. The system was later taken over by Elekta and re-written to a more modern GUI-like environment. Basically, however, the algorithms are probably the same as in the original code. The dose contribution to a point is divided into a scatter component and a primary component. The primary component is computed directly taking into account the path to the point of calculation, i.e. EPL correction. The scatter component is computed summing the contribution from pixels covering the areas of the beam, taking into account the fluence that reaches each pixel and the distance to the skin surface. The division of the beam into a scatter component and primary component is somewhat artificial. Since only the primary component is corrected for inhomogeneity, this division can affect the inhomogeneity correction. The required data are extracted from and/or fitted to measurements (Renner 2006). The XiO system from CMS has two dose calculation models implemented, an FFTbased convolution algorithm (XiO/Conv) and a more exact multigrid fast superposition method (XiO/Super). The implementation of these algorithms has been extensively described (Miften et al 2000, 2001). In short, the FFT convolution uses two invariant point kernels representing the dose transport by primary and scatter which is convolved with the ray-traced 5790 T Knöös et al two-component TERMA. The convolution is performed in the frequency domain to significantly decrease the calculation time. The point dose kernels come from the same data set that OMP and Pinnacle use, i.e., from Mackie et al (1988). For the multigrid superposition method, the same point dose kernels and TERMA are used but during integration the kernels which are now represented in spherical coordinates are modified according to density changes in the irradiated medium. The model is an adaptation of the collapsed-cone convolution discussed above for OMP/CC and Pinnacle/CC and uses 128 directions during the ray tracing. The term multigrid comes from the use of a varying resolution for the calculation grid depending on potential dose gradients, e.g., at high density gradients and at beam edges a finer grid is used. The multigrid superposition model is similar to the adaptive convolution present in the Pinnacle system (not included in this study). The algorithms in this study have been divided into two groups: (a) Models primarily based on EPL scaling for inhomogeneity corrections. Changes in lateral transport of electrons are not modelled. The algorithms in this group are Eclipse/ModBatho and Eclipse/ETAR, OMP/PB, PPLAN and XiO/Conv. (b) Models that in an approximate way consider changes in lateral electron transport. The models in this group are Pinnacle/CC, Eclipse/AAA, OMP/CC and XiO/Super. Further, these will be referred to as type a and type b models, respectively. For the type b models, electron modelling is not explicitly performed; however, since the kernel(s) represents the energy from both primary electrons and scattered photons which are scaled rectilinearly with electron density according to the theorem by O’Connor, an approximate modelling is included. The patient data have in all systems been translated to density (and media when appropriate) from Hounsfield numbers according to the local lookup tables present in the planning systems. All algorithm input data for a study like this should of course be identical; however, due to the complexity of characterization and implementation of beam data for a certain linear accelerator in all the used systems this was not feasible. Instead we decided to use what was available at each clinic; a dual energy x-ray machine with 6 and 15 or 18 MV was chosen. To estimate the uncertainty introduced by not using the same data in all systems, the TPR20,10 was calculated using the relation in IAEA TRS-398 which is based on the dose at 20 and 10 cm depth for a reference situation. In this study, we have chosen an SSD of 100 cm and field size equal to 10 × 10 cm2. Additionally, the lateral or off-axis variation among the implemented beam data was checked at 10 cm depth for a 20 × 20 cm2 beam. 2.3. Benchmark data created by Monte Carlo simulations Ideally all the accelerators used for the planning should be modelled by Monte Carlo to be able to construct a comprehensive benchmark set. This is, however, not easily accomplished, especially since we have studied plans produced for linear accelerators from three different vendors. Instead, an idealized linear accelerator has been constructed, a so-called virtual linear accelerator which has been presented earlier (Wieslander and Knöös 2000, 2003). For this unit, a complete algorithm input data set has been simulated consisting of depth doses, profiles, fluence distribution, output factors, etc, required by the OMP system. These data have been processed and implemented by the vendor according to their standard routine. After validation that algorithm input data and output from the TPS agreed within tolerances, these MC-based virtual accelerators were ‘commissioned’ for producing plans for the included patients in this study. Both the pencil beam and the collapsed-cone convolution algorithms were used in OMP. These models implemented in the OMP system based on MC simulations of input data Dose calculation algorithms for treatment planning in external photon beam therapy 5791 were named OMP/PB-VL and OMP/CC-VL, respectively. A benchmark set for the studied patients was also created directly by MC simulations (MC-VL). The Monte Carlo simulations both for the algorithm input data as well as the benchmark set have been performed with EGSnrc as the engine for BEAMnrc and dosxyznrc from the NRCC group in Ottawa, Canada (Rogers et al 1995, Kawrakow 2000a, 2000b, Kawrakow et al 2004). The virtual linac has a source emitting photons sampled from spectral distributions which produces depth doses with the required penetration. These photons are transported through an ideal collimator with no transmission or scatter production. BEAMnrc was used to create the virtual linac and phase spaces, i.e. the distribution of particles (photons, electrons and positrons), were sampled at 70 cm from the source. The number of particles in the phase space files was from 66 to 95 millions for the used fields. During the patient simulations with dosxyznrc, approximately 250 × 106 of particles were used per field to yield a statistical uncertainty of around 0.5% in the centre of the PTV, i.e., the isocentre. To accomplish this, particles in the phase space files had to be reused several times. Electrons and photons were traced until their kinetic energies fell below the cut-off energies which were 189 and 10 keV, respectively. All other transport parameters were according to the default settings of the code system. The CT data for each case were converted from DICOM to a suitable format for dosxyznrc with the DICOM-RT toolbox (Spezi et al 2002). This software is available within the software package CERR from Washington University, St Louis (Deasy et al 2003). The dose grid resolution during MC was as in most of the TPSs, i.e. the slice spacing in the cranial/caudal direction was 5 millimetres and 5 by 5 millimetres in the transversal plane. CT numbers were converted to density and materials which were available in the BEAMnrc package (used materials were air, lung, soft and bone tissue). One has to note that MC as it is set up in this work calculates the absorbed dose to the tissue contrary to most TPS algorithms that calculate the dose to water of various densities. The OMP/CC algorithm as described above also gives dose distributions as dose to tissue. 2.4. Data comparison and analysis All plans were set up to give an intended absorbed dose per fraction of 2.0 Gy at the specification point, i.e. the isocentre for all cases. The relative number of monitor units (m.u.) for the individual fields was constant for all the TPSs. The number of m.u. for each plan and field was scored. The aim of this study was to study the effects from the modelling in the patient, i.e., the performance of the algorithm. Thus to be able to correct for differences in calibration procedures, a normalization of the data was needed. Each institution was asked to make a ‘reference measurement’ of the dose (Drefernce) and number of m.u. (Mreference) for each field used, SSD = 100 cm at a depth of 10 cm. Based on these numbers, a relationship between absorbed dose and monitor units could be established. This number was then used to normalize the dose per monitor units for all treatment plans: 
N Dreference i=1 Disocentre,i . (1) dplan,normalized =
N Mreference i=1 Mi The summation is performed over all fields i = 1 to N. In this way, differences in m.u. calculation, head scatter modelling and fluence modelling, etc, were removed or at least the influence from these contributions was minimized. But most important is that differences in, for example, loss of scatter for a specific algorithm which is normally compensated by an increased number of m.u. will show up as a lower absorbed dose in the plan compared with 5792 T Knöös et al Table 1. Evaluation of TPR20,10 for the studied TPSs. The source of algorithm input data is given as the type of linac and the nominal x-ray energies. All TPR values are estimated from calculated doses at 10 and 20 cm depth in a reference geometry (SSD = 100 cm, 10 × 10 cm2) by each studied TPS/model. Off-axis ratios are at 10 cm depth which are normalized at the central axis. Treatment planning system/algorithm Type of linac—energy OMP/PB Philips/Elekta SL25 —6/18 MV Philips/Elekta SL25 —6/18 MV 1-Clinac 23C/D —6/18 MV 2-Clinac 2100C/D —6/15 MV 1-Clinac 23C/D —6/18 MV 2-Clinac 2100C/D —6/15 MV 2-Clinac 2100C/D —6/15 MV Primus—6/18 MV Clinac 2100C —6/18 MV Siemens Primus —6/15 MV Siemens Primus —6/15 MV OMP/CC Eclipse/ModBatho Eclipse/ETAR Eclipse/AAA Pinnacle/CC PPLAN XiO/Conv XiO/Super Average SD (%) TPR20,10 Off-axis ratio at 2.5, 5.0 and 7.5 cm 6 MV 15/18 MV 6 MV 15/18 MV 0.665 0.750 1.009; 1.000; 0.980 1.008; 1.017; 1.015 0.655 0.743 1.010; 1.001; 0.982 1.009; 1.016; 1.015 0.647 0.757 1.004; 1.004; 0.990 1.010; 1.018; 1.005 0.647 0.725 1.005; 1.006; 0.994 1.014; 1.021; 1.008 0.647 0.757 1.004; 1.004; 0.990 1.010; 1.018; 1.005 0.647 0.725 1.004; 1.005; 0.993 1.014; 1.021; 1.008 0.647 0.734 1.006; 1.009; 0.989 1.019; 1.027; 1.009 0.650 0.654 0.744 0.756 0.995; 0.991; 0.969 1.002; 0.998; 0.975 1.004; 1.017; 0.997 1.004; 1.005; 0.997 0.654 0.740 1.002; 1.009; 0.989 1.009; 1.034; 1.019 0.652 0.738 1.003; 1.013; 0.993 1.010; 1.036; 1.022 0.653 1.1 0.743 1.4 1.004; 1.004; 0.986 0.4; 0.6; 0.8 1.010; 1.021; 1.009 0.4; 0.9; 0.8 an algorithm that does not model loss of scatter to the same accuracy. This methodology is in line with what is suggested by ESTRO (Mijnheer et al 2005) which was termed m.u. or output normalization. This methodology results in plans that are only comparable within the patient case and the used beam quality. To improve the comparability all plans are normalized against the average dose within the PTV for all type a models for each case and beam quality. In this way, we have kept the variation in dose per m.u. but rescaled the values to the same level, i.e., they vary around 100% which facilitates the interpretation. The benchmark sets from the virtual linac calculations are normalized to the Monte Carlo data, i.e., the average dose in the PTV for the MC-VL is set to 100%. All data have been exported from the various systems either by the DICOM or the RTOG protocol and imported into the CERR software. Evaluation is based on descriptive dose statistics for the PTV and a high priority OAR for each lesion when present. Absorbed doses also reported during the analysis are Di which is the dose received by i% of the structure. 3. Results and discussion A total of 104 treatment plans for the clinical patient cases have been produced which have all been imported into the CERR software for analysis. The evaluated radiotherapy plans include different treatment machines. In table 1, the TPS, accelerator used and the beam quality parameters are presented. The choice of accelerator has mainly been based on the available x-ray energy to be able to perform comparisons with Dose calculation algorithms for treatment planning in external photon beam therapy 5793 the same energies (6 and 15/18 MV x-rays). All planning systems have used algorithm input data which were processed according to the manufactures directives or standards. For one of the systems, the vendor itself has processed and implemented the beam data. As can be seen in table 1, differences exist in beam quality: the TPR20,10 is on average 0.653 and 0.743 for the low (6 MV) and the high energy (15 or 18 MV), respectively. The variation among the systems is close to 1% for both energies in terms of a standard deviation (SD). Interesting to note is that for the SPB models in the Eclipse systems it does not matter which algorithm is chosen, the beam quality index is exactly the same. This supports the idea that this system is basically a corrective-based algorithm, i.e. the dose is modelled from measured data and corrected for scatter and inhomogeneities in certain ways depending on model selected. Thus, a calculation in a water tank should give exactly the same result as a measured depth dose curve. For the AAA model, a different data processing is used and consequently the fitting process leads to small differences. This can also be seen for the other systems where we have more than one algorithm available, i.e. OMP and XiO differences are present since these are model-based algorithms, i.e. no measured depth doses are used during the dose calculation only parameters extracted during the preprocessing are in use. Off-axis ratios at 10 cm depth for a 20 × 20 cm2 field setup with SSD equal to 100 cm in water are also given in table 1. The SDs for the three points and both energies are less or equal to 1%. Of most importance is the point 2.5 cm off-axis which has an SD of 0.5% or less since all used beams include this point. The cases that are influenced by the larger variation for the 5.0 and the 7.5 cm points are the head and neck case and the tangential breast case where the beam sizes are 16 × 18 cm2 and 9 × 17 cm2, respectively. One of the systems (Eclipse) has been used at two different centres and from these results one can perform a rough estimate of the accuracy and feasibility of the approach for comparison used in this study. Choosing the 6 MV x-ray plans (the beam qualities are close to identical) and looking at the ratio of average dose to the PTV between the two systems gives 0.996 (for ModBatho) and 1.001 (for ETAR) for prostate, 1.004 and 1.008 for head and neck, 1.005 and 1.012 for breast and 1.002 and 1.009 for lung. Remembering that the algorithm input data sets are collected by two different groups for different accelerators and equipment which are entered into the system and processed at the different sites, the results are excellent support for the notion that the approach used will have the potential to reveal differences between the systems. For the virtual linac we have the following TPR20,10 ratios: 0.661 for OMP/PB-VL, 0.653 for OMP/CC-VL and 0.669 for MC-VL for the lower energy. For the high energy the ratios are 0.749, 0.741 and 0.744. These values agree rather well with the average value for the included linear accelerators at the participating institutes. A comment regarding the benchmark set based on the virtual accelerator is that the profiles are flatter than for clinical beams, which can have an impact on the average dose in certain structures especially when elongated rectangular beams were used, as for the head and neck and the tangential breast cancer cases. Thus, these data should be used to get the trend between a type a and a type b and how they compare to MC calculation. All treatment plans for both energies and all systems including the implementation of the virtual accelerator and MC calculations are presented in the electronical supplement available at stacks.iop.org/PMB/51/5785. Figures in the supplement will be referred to as Figcmp. 3.1. The prostate case The selected treatment for the prostate gland was a standard four-field box technique with equal weight expressed as number of m.u. Since this site consists of a deep-situated target 5794 T Knöös et al volume (PTV) surrounded by fairly homogeneous tissues, i.e., soft tissues rather similar to water and some smaller regions with denser bone, no dramatic differences were expected when analysing the results from the different TPS/algorithms. All models are expected to perform fairly well in this case. The bone tissue attenuates the beams more than water which should be handled rather similarly in all systems; however, some influence may be present due to the higher atomic number of bone especially if dose to bone tissue is scored. Treatment plans for all models are presented in Figcmp 1 and 2 for both x-ray qualities. Qualitatively, the plans are rather similar independent of the model used for calculation. The average PTV dose given as relative dose per m.u. (cf table 2) for all systems ranges from 98.4 to 101.9% for 6 MV with a mean dose 0.5% higher for the type b compared to the type a calculation models10 . D95 for the PTV range from 96.2 to 100.3% for all studied systems with an average value of 98.2%. For D5, the interval is from 100.3 to 103.4% for the PTV. There was no significant difference between the two types of calculation. Looking at 15/18 MV, similar figures were found with slightly less spread among the systems. The absorbed dose in the PTV was not influenced by the type of algorithm used, nor were any of other the parameters evaluated for the prostate treatment. Having these figures in mind, the conclusion is that for a simple geometry like a four-field box technique in the pelvic area the choice of model/system is not critical for the final dose or dose distribution. Detailed study of each TPS reveals that a move to the more sophisticated model within a system does not have any significant change. This is in agreement with e.g. Aspradakis et al (2003), who showed only minor differences in the pelvic region between a pencil beam and a collapse-cone convolution model implemented in the same TPS. The largest change was found for the OMP system with about 1.5% lower average PTV dose when collapsed-cone convolution was used instead of pencil beam integration. This difference might partly be explained by the fact that OMP/CC contrary to the PB model calculates the dose to the media instead of dose to water since the dose to soft tissue is about 1–2% lower than in water (see further discussion in section 3.6). 3.2. The head and neck case Most treatment planning systems use algorithm input data measured in large water phantoms which differ a lot from the narrow and much smaller neck region. The beam setup chosen for this study represents the first phase in a curative treatment where both sides of the cervical nodes together with a GTV are included in the PTV. The treatment plan consists of two opposed laterals treating the PTV and is intended to be followed by a boost technique to further increase the dose to the GTV (this latter plan is not included in this study). The results for all models are presented in table 3 and in Figcmp 4 and 5. The average dose per m.u. to the PTV decreases by 1% for the low energy if the more sophisticated models of type b were used. This difference was not present for the higher energy which may be due to less scatter in the high energy beam. Looking at the low and the high dose volumes, i.e. D95 and D5, shows no change when moving to the type b models neither for 6 nor 15/18 MV x-rays. The type a models are dominated by pencil beam integration where missing tissue beyond the patient, i.e. the lack of back scatter, is not modelled. This deficiency has been demonstrated fairly well by, e.g., Hurkmans et al (1995). In this geometry, a spread in dose exists among the type a models (cf table 3) which can be seen easily for the Eclipse system with two models of this type included in the study. The ETAR method results in a dose closer to the type b models due to the improved scatter integration which considers the 3D extension of the volume more accurate. The difference 10 The average dose for all the type a models is by definition 100%. Eclipse/ Eclipse/ Mod Eclipse/ Mod Eclipse/ OMP/ XiO/ Batho 1 ETAR 1 Batho 2 ETAR 2 PB PPLAN Conv Average Average values for Eclipse/ OPM/ Pinnacle/ XiO/ values for OMP/ OMP/ type a AAA 2 CC CC Super type b PB-VL CC-VL MC-VL 6 MV 99.5 PTVmean PTV D95 97.2 PTV D5 101.4 PTV D5 − D95 4.2 99.8 98.2 101.4 3.1 99.9 98.2 101.4 3.1 99.7 97.2 101.4 4.2 101.9 100.3 103.4 3.1 98.4 96.2 100.3 4.2 100.7 100 99.3 98.1 102.4 101.7 3.1 3.6 100.6 99.3 102.4 3.1 100.5 99.3 102.4 3.1 100.0 98.2 101.4 3.1 101.1 99.3 102.4 3.1 100.5 99.0 102.1 3.1 100.6 98.6 101.8 3.1 99.1 97.6 100.7 3.1 100.0 97.6 101.8 4.2 15/18 MV 100.9 PTVmean 99.5 PTV D95 102.4 PTV D5 2.9 PTV D5 − D95 100.7 99.5 102.4 2.9 99.7 97.6 101.4 3.8 99.5 97.6 101.4 3.8 100.6 99.5 101.4 1.9 98.8 96.6 100.5 3.8 99.7 100 98.5 98.4 101.4 101.6 2.9 3.2 99.8 98.5 101.4 2.9 98.8 97.6 99.5 1.9 99.4 98.5 100.5 1.9 99.8 98.5 101.4 2.9 99.4 98.3 100.7 2.4 101.0 99.3 102.3 2.9 99.4 97.4 100.3 2.9 100.0 97.4 101.3 3.9 Dose calculation algorithms for treatment planning in external photon beam therapy Table 2. Prostate: data for the PTV for all treatment planning systems and the MC benchmark. The parameters shown are the mean dose to the PTV (PTVmean) and the minimum dose to 95% and 5% of the PTV (PTV D95 and PTV D5). The variation of the dose in the PTV is expressed as PTV D5 − D95. The average value for the dose in PTV calculated by the type a models is by definition 100%. 5795 5796 Table 3. Head and neck: data for the PTV for all treatment planning systems and the MC benchmark. The parameters shown are the mean dose to the PTV (PTVmean) and the minimum dose to 95% and 5% of the PTV (PTV D95 and PTV D5). The variation of the dose in the PTV is expressed as PTV D5 − D95. The average value for the dose in PTV calculated by the type a models is by definition 100%. Eclipse/ Eclipse/ Mod Eclipse/ Mod Eclipse/ OMP/ XiO/ Batho 1 ETAR 1 Batho 2 ETAR 2 PB PPLAN Conv Average Average values for Eclipse/ OPM/ Pinnacle/ XiO/ values for OMP/ OMP/ type a AAA 2 CC CC Super type b PB-VL CC-VL MC-VL 6 MV 101.5 PTVmean PTV D95 73.2 PTV D5 110.4 PTV D5 − D95 37.2 99.7 70.1 109.6 39.5 101.1 72.4 109.6 37.2 98.9 69.3 108.9 39.5 100.0 72.4 109.6 37.2 98.4 68.6 109.6 41.1 100.5 100 70.1 70.9 111.2 109.8 41.1 39.0 99.6 67.8 111.2 43.4 98.6 67.8 108.9 41.1 96.5 65.5 108.9 43.4 101.1 70.1 111.9 41.8 99.0 67.8 110.2 42.4 103.3 75.4 113.6 38.3 101.8 70.5 113.6 43.2 100.0 71.3 111.2 39.9 15/18 MV 98.7 PTVmean 81.2 PTV D95 106.0 PTV D5 PTV D5 − D95 24.8 98.5 79.6 107.6 28.0 101.5 81.2 108.4 27.2 100.5 79.6 109.2 29.6 100.1 82.0 107.6 25.6 99.0 80.4 106.8 26.4 101.6 100 80.4 80.6 110.0 107.9 29.6 27.3 100.6 78.0 109.2 31.2 98.7 78.0 106.0 28.0 97.8 76.4 106.8 30.4 102.1 79.6 110.0 30.4 99.8 78.0 108.0 30.0 103.2 82.4 112.6 30.2 101.6 79.8 111.7 31.9 100.0 80.7 110.0 29.3 T Knöös et al Dose calculation algorithms for treatment planning in external photon beam therapy 5797 between ETAR and ModBatho is less for the higher energy where the amount of scatter is less than for 6 MV. For the XiO system one would expect a change in the same direction going from the FFT-based convolution to the superposition; however, there is about 0.5% increase in dose when moving to the sophisticated model. 3.3. The breast case Studying the tangential treatment of breast after conservative surgery adds another dimension of complexity, i.e. only part of the irradiating beam is actually impinging on the patient plus the fairly limited thickness of the irradiated volume, i.e. the lack of scattering media exists both beyond as well as lateral to the treated volume. A small influence of the low density lung in proximity to the target volume may also add to the complexity. Data for all models are presented in table 4 and plans are shown in Figcmp 6 and 7. The average PTV dose decreases with 0.7% and 1.6% for 6 and 15/18 MV, respectively, when replacing type a with type b models. Studying the dose to the adjacent left lung, especially the high dose area, one can note that D5 decreases significantly for both energies, for 6 MV with 9.5% and for 15/18 MV with 12.2%. The larger change for 15/18 MV may be attributed to the longer range of electrons especially in lung tissue of low density which transport energy further into the lung. Some effect of the same phenomenon may also be present superficially along the patient outline. Altogether these effects give a larger change in average dose to the PTV for the high energy when comparing type a and type b models. Dose–volume histograms (DVHs) for 6 MV treatment of the breast are presented in figure 1. The distribution of average dose among the systems is slightly larger for the type a models probably due to different scatter volume integration techniques. For the type b the variation in average dose is less. One can also note that the DVHs for the PTV are steeper for type a indicating a slightly higher homogeneity which, compared with the more accurate models, is only an illusion. For the estimation of lung dose, the choice of model has no impact on the DVH. Within the type a models a sub group exists based on three-dimensional scatter corrections (i.e. Eclipse/ETAR) and fortunately an EPL-based algorithm (ModBatho) is also implemented in the same systems. By comparing these two algorithms, approximately 3% lower average dose per m.u. to the breast is found for the improved scatter correction. Measurements within anthropomorphic phantoms have been published verifying this discrepancy for EPL-based algorithms (Knöös et al 1986, van Bree et al 1991). The width of the penumbra increases when passing the low density lung region which is nicely modelled by the type b but ignored by the type a implementations, see Figcmp 6 and 7 where plans for the high energy are presented. The distance from the thoracic wall to the 50% isodose is, however, identical for all plans. This is also supported by the volume of the lung receiving 50% or higher dose (cf the DVH in figure 3). Notable also for type b plans is the concave curvature of the isodose(s) just below the medial part of the breast which are present already for the low energy but very pronounced for the high energy. When changes in electron transport are considered a slightly more conformal distribution is obtained. The change for the systems where both types of models are present clearly shows that a more sophisticated model gives up to about 3% lower doses per m.u. to the PTV. As for the head and neck case, the change to an ETAR-based model leads to results similar to the type b models for both low and high energies. The dose to the lung is clearly transported further away with the sophisticated models as D5 indicates. The dose to this volume, for both energies, decreases by 10% or more when type b is used. 5798 Table 4. Breast: data for the PTV and adjacent lung for all treatment planning systems and the MC benchmark. The parameters shown are the mean dose to the PTV (PTVmean) and the minimum dose to 95% and 5% of the PTV (PTV D95 and PTV D5). The variation of the dose in the PTV is expressed as PTV D5 − D95. For the organ at risk, i.e. the left lung, the minimum dose to 5% and 50% (Pulm sin D5 and Pulm sin D50) of the structure is shown. The average value for the dose in PTV calculated by the type a models is by definition 100%. Eclipse/ Eclipse/ Mod Eclipse/ Mod Eclipse/ OMP/ XiO/ Batho 1 ETAR 1 Batho 2 ETAR 2 PB PPLAN Conv Average Average values for Eclipse/ OPM/ Pinnacle/ XiO/ values for OMP/ OMP/ type a AAA 2 CC CC Super type b PB-VL CC-VL MC-VL 6 MV 102.0 PTVmean 93.0 PTV D95 110.7 PTV D5 PTV D5 − D95 17.7 Pulm sin D5 93.0 Pulm sin D50 2.8 99.2 89.8 108.3 18.5 93.8 2.8 101.5 93.0 109.1 16.1 91.4 2.8 98.0 88.9 106.7 17.7 91.4 2.8 100.3 91.4 109.1 17.7 93.8 3.6 97.4 88.1 106.7 18.5 89.8 3.6 101.7 100 93.0 91.0 111.5 108.8 18.5 17.8 95.4 92.6 4.4 3.3 100.3 91.4 110.7 19.3 83.3 4.4 98.9 89.8 108.3 18.5 82.5 3.6 97.6 88.9 106.7 17.7 80.9 3.6 100.3 91.4 109.9 18.5 85.7 4.4 99.3 90.4 108.9 18.5 83.1 4.0 100.9 91.3 110.3 19.0 91.3 2.9 99.6 89.6 109.5 19.8 79.7 2.9 100.0 91.3 108.7 17.4 86.3 2.1 15/18 MV 99.7 PTVmean 92.1 PTV D95 105.1 PTV D5 PTV D5 − D95 13.0 93.0 Pulm sin D5 Pulm sin D50 1.2 99.6 88.9 105.9 17.0 96.2 1.2 101.3 95.4 106.8 11.4 92.1 2.0 100.1 93.0 106.8 13.8 94.6 2.0 100.2 93.8 106.8 13.0 96.2 2.8 97.5 92.1 101.9 9.7 91.3 5.3 101.5 100 95.4 93.0 108.4 105.9 13.0 13.0 97.0 94.4 3.7 2.6 99.9 93.8 106.8 13.0 84.8 2.8 96.3 90.5 101.9 11.4 80.8 3.7 97.5 90.5 103.5 13.0 79.2 2.8 100.2 94.6 107.6 13.0 84.0 3.7 98.4 92.3 104.9 12.6 82.2 3.2 102.3 93.4 109.6 16.2 95.1 2.1 98.8 90.0 105.4 15.4 80.6 3.0 100.0 90.0 107.1 17.1 85.8 3.0 T Knöös et al Dose calculation algorithms for treatment planning in external photon beam therapy 5799 1.0 Eclipse/Mod Batho Eclipse/ETAR OMP/PB PP XiO/Conv Eclipse/Mod Batho Eclipse/ETAR OMP/PB PP XiO/Conv 0.9 Fraction of volume 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 0 10 20 30 40 50 60 70 80 90 100 110 120 80 90 100 110 120 Absorbed dose (%) 1.0 0.9 Eclipse/AAA OMP/CC Pinnacle/CC XiO/Super Eclipse/AAA OMP/CC Pinnacle/CC XiO/Super Fraction of volume 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 0 10 20 30 40 50 60 70 Absorbed dose (%) Figure 1. DVH for the PTV and the adjacent lung for the 6 MV tangential treatment of the breast calculated with type a (upper panel) and type b (lower panel) models. Solid lines are the PTV and the dashed lines are for the lung. 3.4. The pulmonary case The large boost volume in the thorax is treated with five non-equal weighted beams (different number of monitor units for the individual beams). In Figcmp 9 and 10, the plans for 6 MV x-rays and in Figcmp 11 and 12 for the higher energy are shown. The average dose per m.u. to the PTV when 6 MV is used decreases with 2.5% when the models of type b are used. Changing the energy to 15/18 MV increases this difference to 3.7% (cf table 5). This figure is also reflected in the systems where we have both types of models available; Eclipse has about 3% difference between AAA and the average of the modified Batho (2.8%) and the ETAR (3.3%), OMP 4.6% from pencil beam to collapsed cone and XiO has a difference of 5800 Table 5. Lung: data for the PTV, the adjacent (sin) and the contralateral (dx) lung for all treatment planning systems and the MC benchmark. The adjacent lung volume does not include the PTV. The parameters shown are the mean dose to the PTV (PTVmean) and the minimum dose to 95% and 5% of the PTV (PTV D95 and PTV D5). The variation of the dose in the PTV is expressed as PTV D5 − D95. For the organs at risk, i.e. the lungs, the minimum dose to 1%, 5% and 50% (Pulm sin/dx D1, Pulm sin/dx D5 and Pulm sin/dx D50) of the structures is shown. The average value for the dose in PTV calculated by the type a models is by definition 100%. Average Average values for Eclipse/ OPM/ Pinnacle/ XiO/ values for OMP/ OMP/ type a AAA 2 CC CC Super type b PB-VL CC-VL MC-VL 100.1 93.2 105.5 12.3 7.4 31.1 44.2 12.2 103.0 106.7 99.8 92.3 106.4 14.0 7.4 30.2 44.2 12.3 104.4 108.7 99.9 93.2 105.5 12.3 7.4 31.1 44.2 15.7 102.0 105.9 98.9 91.5 105.5 14.0 7.4 30.2 43.3 15.7 102.6 106.8 100.3 93.2 106.4 13.1 8.3 31.1 44.2 10.8 105.1 108.4 97.0 89.7 102.8 13.1 6.6 29.3 41.6 14.4 100.4 104.2 103.9 95.8 111.6 15.8 8.3 33.7 45.1 18.5 108.6 114.4 100 92.7 106.2 13.5 7.6 30.9 43.8 14.2 103.7 107.9 97.8 91.5 102.8 11.4 7.4 30.2 41.6 20.1 96.2 100.0 96.7 90.6 102.0 11.4 7.4 29.3 41.6 19.4 94.6 98.3 96.5 89.7 102.0 12.3 7.4 29.3 40.7 17.8 95.2 98.9 99.0 93.2 104.6 11.4 7.4 30.2 43.3 21.5 99.1 102.7 97.5 91.3 102.8 11.6 7.4 29.8 41.8 19.7 96.3 100.0 101.6 94.7 108.3 13.6 7.7 30.3 43.9 13.6 107.5 111.9 98.2 91.9 103.7 11.8 7.7 28.5 41.2 19.5 96.8 101.2 100.0 93.8 105.5 11.8 7.7 29.4 42.1 17.2 102.0 105.9 15/18 MV PTVmean 99.9 96.2 PTV D95 PTV D5 103.0 PTV D5 − D95 6.8 Pulm dx D50 7.2 29.2 Pulm dx D5 40.3 Pulm dx D1 Pulm sin D50 14.0 99.7 Pulm sin D5 Pulm sin D1 102.2 101.5 96.2 106.4 10.2 8.1 30.1 42.0 14.4 103.3 106.5 99.9 95.4 103.8 8.5 7.2 30.1 42.0 15.8 100.0 103.0 100.4 95.4 105.5 10.2 8.1 30.1 42.0 16.1 102.1 105.4 99.2 94.5 103.8 9.3 8.1 30.1 41.1 13.7 101.9 104.4 96.6 92.0 100.4 8.5 8.1 29.2 41.1 17.8 97.2 99.6 102.4 97.1 108.1 11.0 8.1 31.8 42.8 18.4 104.9 109.1 100 95.2 104.4 9.2 7.8 30.1 41.6 15.7 101.3 104.3 97.1 92.8 100.4 7.6 7.2 29.2 40.3 20.3 93.9 97.3 94.5 89.4 97.9 8.5 7.2 27.5 38.6 21.6 89.6 93.4 95.7 90.3 99.6 9.3 7.2 27.5 38.6 18.9 90.2 94.7 97.8 93.7 101.3 7.6 8.1 29.2 40.3 21.7 93.8 98.1 96.3 91.5 99.8 8.3 7.4 28.4 39.4 20.6 91.9 95.9 103.7 98.3 108.2 9.9 42.7 31.0 42.7 15.3 106.8 110.1 98.7 93.0 102.8 9.9 39.1 27.4 39.1 22.3 93.8 98.3 100.0 94.7 103.7 9.0 40.0 28.3 40.0 21.6 97.2 101.4 6 MV PTVmean PTV D95 PTV D5 PTV D5 − D95 Pulm dx D50 Pulm dx D5 Pulm dx D1 Pulm sin D50 Pulm sin D5 Pulm sin D1 T Knöös et al Eclipse/ Eclipse/ Mod Eclipse/ Mod Eclipse/ OMP/ XiO/ Batho 1 ETAR 1 Batho 2 ETAR 2 PB PPLAN Conv Dose calculation algorithms for treatment planning in external photon beam therapy 5801 4.6% from FFT convolution of an invariant kernel to superposition of a scalable point kernel. The high dose volume within the PTV decreases with 3.4 and 4.6% moving to type b models for the two studied energies. The doses reported for the adjacent lung exclude the volume occupied by the PTV. Looking closer at the isodoses in proximity to the target reveals that they have a larger separation when more sophisticated models are used, especially for the high energy (cf Figcmp 11 and 12). This is also supported by the D1 and the D5 parameters for the left lung which decreases by approximately 8% for both energies. The results for the pulmonary case support the very common practice to not treat patients with tumours in the lung with high energy photon beams due to the inaccurate dose calculation models present in most TPSs available up to now. When the more accurate algorithms of type b are getting more common, the higher energies may be more attractive due to their higher penetration. The widening of the penumbra, etc, may, however, still be a limitation which can result in larger irradiated volumes. 3.5. Monte Carlo benchmark with the virtual accelerator Plans for all cases with the virtual accelerator are presented in the electronic supplement for the implementation in OMP (OMP/PB-VL and OMP/CC-VL) as well as for the Monte Carlo calculations (MC-VL). 3.5.1. Prostate. The calculations with the virtual accelerator produced plans which are within ±1% of the MC-VL for the OMP/PB-VL and OMP/CC-VL algorithms, cf table 2. Translating this information to the evaluation of the clinical systems gives that further improvements in calculations in this region of the body are not expected. Thus, it may be concluded that the choice of algorithm is not a question for treatments in the pelvic region regarding neither the relative distribution nor the absolute dose per m.u. 3.5.2. Head and neck. The benchmark test gives that OMP/PB-VL is 3.3% and 3.2% higher than the MC-VL for both energies. The difference for the OMP/CC-VL algorithm is less (on average 1.7%). The difference is an indication that the scatter integration for both TPS models probably is inferior. This is more pronounced for the PB model which is preintegrated approximately to 50 cm in the longitudinal direction; thus it does not consider the lack of backscatter from the missing tissue on the exit side. The point-kernel-based CC model performs better. Interesting is that the difference between type a and type b models is less but looking in particular at the clinical OMP system (based on real accelerator data) reveals that there is a difference between OMP/PB and OMP/CC of about 1.4% for both energies. This is in close agreement with the benchmark calculations. 3.5.3. Breast. For 6 MV x-rays, the benchmark set (cf figure 2, Figcmp 8 and table 4) gives that the OMP/PB-VL and the OMP/CC-VL result in average PTV doses within 1% (+0.9% for PB and −0.4% for CC) of the MC-VL. The results for the higher energy differ slightly more (+2.3% and −1.2%). The result for the higher energy and the PB model is probably due to a slight disequilibrium in electron transport along the lung and soft tissue interface which is not considered by the algorithm. The clinical system, i.e. OMP/PB, differs with approximately the same figures, i.e. 1.4% for 6 MV and 3.9% for 15/18 MV. 3.5.4. Lung. The MC benchmark set based on the virtual linac for the pulmonary case yields that the difference between OMP/PB-VL and OMP/CC-VL for 6 MV is 3.4% compared 5802 T Knöös et al Figure 2. The results for the MC benchmark for 6 MV (upper row) and 15/18 MV (lower row) for the breast cancer treatment. OMP/PB-VL: left panels, OMP/CC-VL: middle panels and MC-VL: right panels. The PTV is outlined in light blue. with 3.6% for the clinical OMP/PB and OMP/CC models. The figures for 15/18 MV are 5.0% and 4.6%, respectively. The PB plans are about 2% and 4% higher than the MC plan and CC is approximately 1.5% lower. For the breast case above, a lower dose was also found for the CC model combined with high energy compared to the MC simulations. This difference may be a consequence of the approximate handling of electron transport that exists even in a sophisticated method like the collapse-cone convolution. It has been shown that the CC model in OMP (or in those days Helax-TMS) for high energy x-rays combined with low density has a tendency to spread out the energy too much resulting in a wider penumbra compared to measurements (Aspradakis et al 2003). Nisbet et al (2004) have showed that CC may overestimate the dose in low density regions up to 6% outside the field when 15 MV x-rays are used. Francescon et al (2000) pointed out that the collapsedcone convolution approximation could ‘potentially corrupt the accuracy of the convolution– superposition principle, especially in lateral electronic non-equilibrium conditions’. The main reason is probably the simplification of using rectilinear scaling along the collapsed cones of the kernel representing the dose contributions from electrons. 3.6. Dose to tissue versus dose to water All dose calculation algorithms usually calculate the dose to a cavity of water fulfilling the Bragg–Gray criteria. In this case, we have MC calculations and these are calculating the dose to the actual tissue. A straight comparison between TPS and MC can consequently not be performed without considering this conceptual difference. To complicate this further, the collapsed-cone convolution in the OMP system also calculates dose to tissue as shown in, e.g., Wieslander and Knöös (2003). To evaluate this, the prostate plan with 15/18 MV x-rays has been re-calculated in our MC system where all materials were replaced with water. The electron density given by the CT number was kept, thus the MC simulation will be performed in a geometry similar to what is used in most TPSs. The difference in absorbed dose between tissue and water is given by the mass-stopping power for the electrons set in motion by the x-ray photon interaction with the irradiated media. Soft tissue is often approximated by water; however, a small difference in elemental composition yields a small difference Dose calculation algorithms for treatment planning in external photon beam therapy 5803 Figure 3. The four-field box technique for the prostate case with MC calculations for all tissues mapped as tissues (upper left) and mapped as water (upper right). The difference is shown in the lower panel. The isodoses shown are 115, 112.5, 110, 107.5, 105, 102.5, 100, 97.5, 95, 90, 70, 50, 30, 10, 5 and 2%. The dose in cortical bone at the blue dot is −4% and the average dose in the PTV is −0.6% compared to the dose to water plan. in mass-stopping power of about 1% for all energies applied in external radiotherapy (cf NIST website11 for mass-stopping powers). Looking at the prostate simulations we found a difference in dose of about 0.6%, see figure 3. Taking a point anterior to the prostate in a bony area where cortical bone dominates gives a difference of about 6% and in a less cortical bone area 4% which is also supported by the difference in mass-stopping power. According to data from NIST, the mass-stopping power ratios between dense bone, soft bone and soft tissue to water are 0.90, 0.97 and 0.99, respectively. For head and neck treatments large air cavities are present. The mass-stopping powers for air and water are significantly different mainly due to a difference of several magnitudes in density. This should show up as a lower dose in regions where air is present. In figure 4, a cross section at a level where large air cavities are present is shown. Two point doses are shown within the ‘low density’ volume with doses of 100% (point (a)) and 87% (point (b)). The difference between these two points which are in a low density region is the elemental composition of the voxels. The voxel at (b) is mapped as air where the other at (a) is soft tissue. The importance of the mapping and rebinning of the CT information in the case of MC calculations and/or an algorithm calculating dose to tissue must not be underestimated (see also Verhaegen and Devic (2005)). It is of importance doing comparisons like this that the difference in absorbed dose is considered in the evaluation. For soft tissue, where we are within the uncertainty of the calculations, this is not a big problem, but in regions such as dense bone where the elemental composition differs significantly from water with the high Z components this is of high importance. Other regions of importance are when cavities of air are included in the PTV 11 http://physics.nist.gov/PhysRefData/Star/Text/contents.html. 5804 T Knöös et al (a) (b) Figure 4. Dose distribution for the MC benchmark set for 6 MV in the head and neck region. The cross section is positioned at 6 cm cranial of the isocentre. The absorbed dose at point (a) is 100% and at (b) 87%. The displayed isodose lines are 90, 95, 100 and 105%. since air has a rather different elemental composition to soft tissue. This problem has been discussed extensively by, e.g., Siebers et al (2000). 4. Conclusions The differences between the models are in many situations small and, if users are aware of limitations in their TPS, fairly accurate dose distributions will be predicted also with the type a-based systems. For planning of prostate and probably other lesions in the pelvic area, in principle any system will do, since the differences between type a and type b models are insignificant (see also Aspradakis et al (2003)). Increasing the complexity, as in the case for head and neck lesions where a fairly homogeneous volume exists but is limited in size, will favour moving to a more sophisticated model. If three-dimensional scatter integration is part of a type a model (e.g. ETAR), the results are close to those by the type b group. Looking at the tangential breast treatment, a 1–3% difference exists between the two groups of different calculational approaches. For the most complicated situation, the pulmonary boost technique, about 4% difference exists if 15/18 MV x-rays are used, but this is decreased to 2.5% if a lower energy (i.e. 6 MV) is used instead. The MC calculations in the benchmark suggest that models based on convolution with rectilinear scaling in the lateral directions which we have in the type b group are more accurate than those based on corrections along the fan line only (type a). One may also conclude that in certain situations the approximate electron modelling used in these type b methods results in Dose calculation algorithms for treatment planning in external photon beam therapy 5805 non-accurate plans. For example, the prolonged range in low density media is overestimated compared to MC calculations. One should have in mind that the results found in this study can be influenced by the patient cases selected. The design of the study was performed to avoid this, but drawing general conclusions of the performance in, e.g. the head and neck region or any other site studied should still be performed carefully. For example, the amount of air-filled volumes and their locations may have an unexpected influence not revealed in this study. Another example is the thoracic treatment where the results may be case specific. To be able to keep the uncertainty in the final dose to the patient as low as discussed in several publications (see e.g. Mijnheer et al 2005 or the IAEA TRS 430), which is in the order of 3%, all contributing sources have to be minimized. The dose calculation is one of the major contributors and moving to type b models will in general reduce the uncertainty in the delivered dose to the patient. When evaluating treatment planning systems one must remember that they are very complex systems which may approach the problem in different ways. In this study, we have tried to isolate only the heterogeneity and the scatter volume handling. Other aspects are head scatter modelling, penumbra, handling of wedges (hard and soft) and other beam modifiers which are not studied in this work. Today one also has to consider the models/algorithms for optimizing intensity-modulated beams both for dynamic and static MLC delivery. Especially, IMRT fields comprising of small segments may bring up new problems for the dose calculation models. Acknowledgments We hereby greatly appreciate the contribution from the Lund University Hospital Research funds. All colleagues at the participating centres are also acknowledged for their contributions and continuous support. The reviewers are also acknowledged for their valuable contributions and proposed changes of the manuscript. References AAPM Report 85 2004 Tissue inhomogeneity corrections for MV photon beams Report of Task Group No. 65 of the Radiation Therapy Committee of the American Association of Physicists in Medicine (Madison, WI: Medical Physics Publishing) Ahnesjö A 1989 Collapsed cone convolution of radiant energy for photon dose calculation in heterogeneous media Med. Phys. 16 577–92 Ahnesjö A and Aspradakis M M 1999 Dose calculations for external photon beams in radiotherapy Phys. Med. Biol. 44 R99–155 Ahnesjö A, Saxner M and Trepp A 1992 A pencil beam model for photon dose calculation Med. Phys. 19 263–73 Ahnesjö A and Trepp A 1991 Acquisition of the effective lateral energy fluence distribution for photon beam dose calculations by convolution models Phys. Med. Biol. 36 973–85 Ahnesjö A, Weber L, Murman A, Saxner M, Thorslund I and Traneus E 2005 Beam modeling and verification of a photon beam multisource model Med. Phys. 32 1722–37 Aspradakis M M, Morrison R H, Richmond N D and Steele A 2003 Experimental verification of convolution/ superposition photon dose calculations for radiotherapy treatment planning Phys. Med. Biol. 48 2873–93 Cheng C W, Das I J, Tang W, Chang S, Tsai J S, Ceberg C, De Gaspie B, Singh R, Fein D A and Fowble B 1997 Dosimetric comparison of treatment planning systems in irradiation of breast with tangential fields Int. J. Radiat. Oncol. Biol. Phys. 38 835–42 Clarkson J R 1944 A note on depth doses in fields of irregular shape Br. J. Radiol. 14 255 Cunningham J R 1972 Scatter–air ratios Phys. Med. Biol. 17 42–51 5806 T Knöös et al Cunningham J R, Shrivastava P N and Wilkinson J M 1972 Program IRREG-calculation of dose from irregularly shaped radiation beams Comput. Programs Biomed. 2 192–9 Deasy J O, Blanco A I and Clark V H 2003 CERR: a computational environment for radiotherapy research Med. Phys. 30 979–85 Francescon P, Cavedon C, Reccanello S and Cora S 2000 Photon dose calculation of a three-dimensional treatment planning system compared to the Monte Carlo code BEAM Med. Phys. 27 1579 Haedinger U, Krieger T, Flentje M and Wulf J 2005 Influence of calculation model on dose distribution in stereotactic radiotherapy for pulmonary targets Int. J. Radiat. Oncol. Biol. Phys. 61 239–49 Hurkmans C, Knöös T, Nilsson P, Svahn-Tapper G and Danielsson H 1995 Limitations of a pencil beam approach to photon dose calculations in the head and neck region Radiother. Oncol. 37 74–80 IAEA 2005 Commissioning and quality assurance of computerized planning systems for radiation treatment of cancer TRS-430 (Vienna: IAEA) IEC 1996 Radiotherapy Equipment—Coordinates, Movements and Scales International Electrotechnical Commission CEI/IEC 1217 Kawrakow I, Rogers DW and Walters BR 2004 Large efficiency improvements in BEAMnrc using directional bremsstrahlung splitting Med. Phys. 31 2883–98 Kawrakow I 2000a Accurate condensed history Monte Carlo simulation of electron transport: I. EGSnrc, the new EGS4 version Med. Phys. 27 485–98 Kawrakow I 2000b Accurate condensed history Monte Carlo simulation of electron transport: II. Application to ion chamber response simulations Med. Phys. 27 499–513 Keall P and Hoban P 1995 Accounting for primary electron scatter in x-ray beam convolution calculations Med. Phys. 22 1413–8 Keall P J and Hoban P W 1996 Superposition dose calculation incorporating Monte Carlo generated electron track kernels Med. Phys. 23 479–85 Knöös T, Ahlgren L and Nilsson M 1986 Comparison of measured and calculated absorbed doses from tangential irradiation of the breast Radiother. Oncol. 7 81–8 Knöös T, Ahnesjö A, Nilsson P and Weber L 1995 Limitations of a pencil beam approach to photon dose calculations in lung tissue Phys. Med. Biol. 40 1411–20 Knöös T, Ceberg C, Weber L and Nilsson P 1994 Dosimetric verification of a pencil beam based treatment planning system Phys. Med. Biol. 39 1609–28 Mackie T R, Bielajew A F, Rogers D W O and Battista J J 1988 Generation of photon energy deposition kernels using the EGS Monte Carlo code Phys. Med. Biol. 33 1–20 Mackie T R, Scrimger J W and Battista J J 1985 A convolution method of calculating dose for 15-MV x rays Med. Phys. 12 188–96 McNutt T R, Mackie T R, Reckwerdt P, Papanikolaou N and Paliwal B R 1996 Calculation of portal dose using the convolution/superposition method Med. Phys. 23 527–35 Miften M, Wiesmeyer M, Kapur A and Ma C M 2001 Comparison of RTP dose distributions in heterogeneous phantoms with the BEAM Monte Carlo simulation system J. Appl. Clin. Med. Phys. 2 21–31 Miften M, Wiesmeyer M, Monthofer S and Krippner K 2000 Implementation of FFT convolution and multigrid superposition models in the FOCUS RTP system Phys. Med. Biol. 45 817–33 Mijnheer B, Olszewska A, Fiorino C, Hartmann G, Knöös T, Rosenwald J C and Welleweerd H 2005 ESTRO Quality assurance of treatment planning systems. Practical examples for non IMRT photon beams ESRTO Booklet no 7 (Brussels: ESTRO) Nisbet A, Beange I, Vollmar H S, Irvine C, Morgan A and Thwaites D I 2004 Dosimetric verification of a commercial collapsed cone algorithm in simulated clinical situations Radiother. Oncol. 73 79–88 O’Connor J 1957 The variation for scattered x-rays with density in an irradiated body Phys. Med. Biol. 1 352–69 Paelinck L, Reynaert N, Thierens H, De Neve W and De Wagter C 2005 Experimental verification of lung dose with radiochromic film: comparison with Monte Carlo simulations and commercially available treatment planning systems Phys. Med. Biol. 50 2055–69 Papanikolaou N, Mackie T R, Meger-Wells C, Gehring M and Reckwerdt P 1993 Investigation of the convolution method for polyenergetic spectra Med. Phys. 20 1327–36 Renner W D 2006 Personal communication Rogers D W, Faddegon B A, Ding G X, Ma C M, We J and Mackie T R 1995 BEAM: a Monte Carlo code to simulate radiotherapy treatment units Med. Phys. 22 503–24 Siebers J V, Keall P J, Nahum A E and Mohan R 2000 Converting absorbed dose to medium to absorbed dose to water for Monte Carlo based photon beam dose calculations Phys. Med. Biol. 45 983–95 Sontag M R and Cunningham J R 1978a Clinical application of a CT based treatment planning system Comput. Tomogr. 2 117–30 Dose calculation algorithms for treatment planning in external photon beam therapy 5807 Sontag M R and Cunningham J R 1978b The equivalent tissue–air ratio method for making absorbed dose calculations in a heterogeneous medium Radiology 129 787–94 Spezi E, Lewis D G and Smith C W 2002 A DICOM-RT-based toolbox for the evaluation and verification of radiotherapy plans Phys. Med. Biol. 47 4223–32 Starkschall G, Steadham R E Jr, Popple R A, Ahmad S and Rosen I I 2000 Beam-commissioning methodology for a three-dimensional convolution/superposition photon dose algorithm J. Appl. Clin. Med. Phys. 1 8–27 Storchi P, van Battum L J and Woudstra E 1999 Calculation of a pencil beam kernel from measured photon beam data Phys. Med. Biol. 44 2917–28 Storchi P and Woudstra E 1995 Calculation models for determining the absorbed dose in water phantoms in off-axis planes of rectangular fields of open and wedged photon beams Phys. Med. Biol. 40 511–27 Storchi P and Woudstra E 1996 Calculation of the absorbed dose distribution due to irregularly shaped photon beams using pencil beam kernels derived form basic beam data Phys. Med. Biol. 41 637–56 Ulmer W and Kaissl W 2003 The inverse problem of a Gaussian convolution and its application to the finite size of the measurement chambers/detectors in photon and proton dosimetry Phys. Med. Biol. 48 707–27 Ulmer W, Pyyry J and Kaissl W 2005 A 3D photon superposition/convolution algorithm and its foundation on results of Monte Carlo calculations Phys. Med. Biol. 50 1767–90 van Bree N A, van Battum L J, Huizenga H and Mijnheer B J 1991 Three-dimensional dose distribution of tangential breast treatment: a national dosimetry intercomparison Radiother. Oncol. 22 252–60 van’t Velt A 1997 Analysis of accuracy in dose and position in calculations of a treatment planning system for blocked photon fields Radiother. Oncol. 45 245–51 Verhaegen F and Devic S 2005 Sensitivity study for CT image use in Monte Carlo treatment planning Phys. Med. Biol. 50 937–46 Webb S and Fox R A 1980 Verification by Monte Carlo methods of a power law tissue–air ratio algorithm for inhomogeneity corrections in photon beam dose calculations Phys. Med. Biol. 25 225–40 Wieslander E and Knöös T 2000 A virtual linear accelerator for verification of treatment planning systems Phys. Med. Biol. 45 2887–96 Wieslander E and Knöös T 2003 Dose perturbation in the presence of metallic implants: treatment planning system versus Monte Carlo simulations Phys. Med. Biol. 48 3295–305 Woo M and Cunningham J 1990 The validity of the density scaling method in primary electron transport for photon and electron beams Med. Phys. 17 187–94 Yu C X, Mackie T R and Wong J W 1995 Photon dose calculation incorporating explicit electron transport Med. Phys. 22 1157–65