Journal of Cancer Research and Therapeutics Medknow Publications on behalf of the Association of Radiation Oncologists of India (AROI) ISSN: 0973-1482 EISSN: 1998-4138 Vol. 1, Num. 1, 2005, pp. 12-20

Journal of Cancer Research and Therapeutics, Vol. 1, No. 1, January-March, 2005, pp. 12-20

Review Article

Advantages of multiple algorithm support in treatment planning system for external beam dose calculations

Animesh

Kirloskar Theratronics Pvt. Limited, Mumbai

Correspondence Address:Kirloskar Theratronics Pvt. Limited, Mumbai - 400057
animesh@kirloskarmedical.com

Code Number: cr05005

Abstract

The complexity of interactions and the nature of the approximations made in the formulation of the algorithm require that the user be familiar with the limitations of various models. As computer power keeps growing, calculation algorithms are tending more towards physically based models. The nature and quantity of the data required varies according to the model which may be either measurement based models or physical based models. Multiple dose calculation algorithm support found in XiO Treatment Planning System can be used to advantage when choice is to be made between speed and accuracy. Thus XiO allows end users generate plans accurately and quickly to optimize the delivery of radiation therapy.

Keywords: Treatment Planning System, Dose Calculation Algorithm, Convolution and Superposition

The functionality and quality of any Treatment Planning System (TPS) is dependent on the type of algorithms used in the different steps of the planning process. An algorithm is defined as sequence of instructions that operates on a set of input data, transforming that information into a set of output results that are of interest to the user. A number of algorithms are used in the treatment planning process. The most well known algorithm is the dose calculation algorithm that generates the dose at any point within the patient while taking into account the patient and beam (or source) characteristics.

Many different types of dose calculation algorithm are used in modern TPSs. Early TPS calculation models were based on a simple tabular representation of the dose distribution that was obtained directly from beam measurements. Over the years, as calculation models have become more sophisticated and computer power has grown, TPS calculation algorithms have progressively evolved towards more physically based models. The most advanced current algorithms now are based on the Monte Carlo approach, in which the histories of many millions of photons are traced as they interact with matter using basic physics interactions. There is a full range of possibilities between table based models and Monte Carlo models. For every algorithm, the quality of the dose representation is strongly dependent on the data or parameters used by the algorithm.

The nature and quantity of the data required varies according to the model. Typically, for measurement based models a lot of tables are required, whereas for physical based models only some parameters are necessary.

Knowledge at some level of the various algorithms used within the TPS can help the user understand the capabilities and limitations of the specific algorithms; also, can help the user diagnose TPS problems and in developing a QA protocol.

It is of paramount importance to understand the general principles of the model and its implementation details. The model parameters and input data have a significant impact on the accuracy of the calculation results. Each model suffers from limitations; for example, a model may not be able to deal with inhomogeneity corrections, and if the issue is considered to be important a manual correction for the inhomogeneities must be used. On the other hand, an inhomogeneity correction might be part of the model and give acceptable results in a number of circumstances but may deviate significantly in other situations. Finally, even if the model is able to account for a given physical effect, the actual implementation in the treatment planning software is often simplified, leading to inaccurate or unexpected results for certain situations.

XiO TPS from CMS Inc., USA, is a combination of teletherapy with IMRT capability and brachytherapy treatment planning system.

In teletherapy module dose calculation algorithms that are used to calculate the dose throughout the calculation area, volume, or set of interest points are:

· Modified Clarkson Sector Integration
· Fast Fourier Transform (FFT) Convolution
· Superposition
· Fast Superposition
· Electron 3D Pencil Beam
· Proton Broad Beam and Pencil Beam

The XiO Fast-Fourier Transform (FFT) Convolution algorithm and the Superposition (Wiesmeyer and Miften)[1] algorithms are similar in that they both compute the dose by convolving the total energy released in the patient with Monte Carlo-generated energy deposition kernels computed by Mackie et al.[2] The major differences are that FFT Convolution does not calculate dose as accurately as Superposition in the presence of tissue in-homogeneities. Proton beam algorithm is not being discussed as it has little relevance with respect to the radiotherapy being practiced in this part of the world.

EXTERNAL BEAM DOSE CALCULATION ALGORITHMS IN XIO SYSTEM

Clarkson Sector Integration Algorithm

The most basic Modified Clarkson Sector Integration Algorithm uses patient data, treatment machine data and setup information to simulate dose distributions inside the patient. The patient data consists of relative electron density information that represents a section of the patient. These will have previously been generated by assigning density values to areas enclosed by traced contours or by applying the CT to relative electron density conversion to patient image data. Generally speaking, for a 3-D plan, only transverse patient data is used in dose calculation.

The treatment machine data required for the Clarkson algorithm consists of a set of Tissue Ratios (TXRs) which are TAR or modified TPR values, a set of diagonal Off-Center Ratios (OCDs), various machine (energy) specific constants and information about the various beam modifiers.

If the beam has a customized port or multileaf collimator, the outlines and transmissions for the blocked areas need to be entered. The algorithm takes into account primary dose correction for inhomogeneities in the patient and bolus and the transmission by the wedge, compensators and blocks. Additionally, the Clarkson algorithm takes into account scatter modifications due to field shape. It does not take into account scatter modifications due to differences in field intensity (wedges), patient density or surface curvature. The Clarkson scatter calculation does not accurately model dose for bifurcated structures because it does not account for the scatter reduction from the air space between the separated contours. Each beam is calculated separately and then the beams are summed together. The total dose from all beams represents the dose to the patient.

Physics and Assumptions of the Clarkson Algorithm

The scatter at a point from one irregular-shaped field is calculated by dividing the field into sectors and integrating the scatter reaching the calculation point from each of the sectors.

These following assumptions are made in this scatter calculation.

i) The Clarkson algorithm assumes a homogeneous patient for the purpose of scatter calculation.
ii) The primary component of the radiation beam is uniform across the beam. Since the scatter contribution only depends on the distance between the dose point and the field edge, and not on the relative intensity of the beam at the source of scatter, the scatter contributions may be inaccurate for cases where the beam profile is not flat.
iii) None of the block contours overlaps the other blocks. If overlapping regions do exist, the scatter in these regions will be considered twice.
iv) There is a homogeneous scattering medium over the region defined by the field. Therefore, if the dose point lies at the edge of a vertical boundary of the patient, the scatter will be overestimated at this point because of missing tissue.

FFT Convolution Algorithm

The energy deposition kernels of Mackie et al.[2] must be interpolated from spherical to Cartesian coordinates on a common grid with the TERMA (Total Energy Released per unit MAss) to perform FFT convolution. Sampling and interpolation of kernels from spherical to Cartesian coordinates is complicated by steep kernel gradients. Adaptive quadrature techniques ensure that the correct energy at and near the interaction point is represented in the Cartesian coordinates.

Comparison of calculation results indicates that incorrect doses are obtained if the effect of scatter from neighbors is omitted over a large enough region. It is important that patient data be represented over a 3-D volume because the scatter calculated at a point is based on a 3-D volume of scattering medium.

The required volume over which scatter or kernel contributions must be included, and the maximum volume used on the XiO system, is about 30 cm in the forward direction, 5 cm in the backscatter direction, and twice the field size dimension laterally (essentially the contributions from all interaction points must be accumulated). Sharpe and Battista[3] report using these ranges, and Mackie et al4],5] report the same required lateral range.

Including dose contributions over such a large area requires significant computation time. This computation time can be reduced by performing two separate calculations:

· One with the primary kernel for which the calculation is performed at high-resolution but over a small region
· The other with the scatter kernel, where the calculation is performed at a lower resolution but over a larger area, as proposed by Mackie et al.4],5]

This approach is possible since the primary kernels have extremely large gradients close to the interaction point, but they make no contributions beyond a few centimeters from the interaction point, whereas the scatter kernels have smaller gradients but contribute dose over a much larger range.

The XiO system performs a separate high and low-resolution FFT calculation for the primary and scatter kernels, achieving a time savings of about 65 percent over performing a single calculation at high-resolution.

Superposition Algorithm

The XiO superposition dose deposition method is an adaptation of "collapsed cone" dose calculation method.[6] As with FFT convolution, all superposition calculations are done in beam coordinates, and dose in beam coordinates is interpolated to the user specified calculation volume.

Spherical kernel representation

It is possible for superposition algorithms to directly emulate the kernel calculation process; that is, to calculate deposited energy by spreading energy released (TERMA) at interaction points to points in the volume of interest according to the distribution implied by the kernel. This method is known as the "interaction point of view."

However, most do not do this for reasons of efficiency. To calculate dose for these few points using this direct method, a complete distribution of doses for every TERMA point would have to be calculated. In the process, dose at all points in the volume would be calculated, not just the few of interest.

To calculate dose only at the points in which we were interested, the kernel probability distribution function must be inverted . By inversion it is meant that the forward scatter part of the kernel points towards the source, rather than away from the source. Using the inverted kernel, for any dose point, the probability distribution function embodied by the kernel is used to collect TERMA from all appropriate points in the volume-this is known as the "dose point of view." Thus, dose may only be calculated at the specified minimum number of points.

Kernel and inverted kernel.

As in FFT convolution, superposition dose deposition is in two parts: calculation of primary electron dose using "primary" kernels, and calculation of scattered photon dose using "scatter" kernels.

Spherical kernels, or "dose spread arrays" are cylindrically symmetric and defined in terms of rays traced along zenith and azimuth angles. Solid lines represent rays and dashed lines represent the boundaries of solid angles (or cones) that are implied.

Theoretically, any number of evenly-spaced azimuth angles greater than four may be specified. The combination of zenith and azimuth angles determines the direction of ray traces that in turn determine the implied collapsed cones upon which calculations are based. The lower bound of four azimuth angles is required to ensure that there is coverage of dose for a full 360 degrees around the axis of symmetry of the kernel.

Accuracy improvements were made possible by optimizing the distribution of ray traces. Since most higher-energy kernel distributions are forward peaked, we have increased ray tracing in the forward direction and decrease it in the backward direction without significantly affecting accuracy in the backward direction.

Calculation of deposited dose requires that spherical kernels be recast in cumulative form; that is, energy for each spherical "cone" is integrated from the interaction point outwards. Energy summed in such a way for these kernels quickly reaches a plateau, after which no more additional energy is added.

Fast Superposition Algorithm

The spherical kernel computation has been augmented with the ability to combine (select and sum) adjacent zenith rays in the kernel. Thus, it is possible to limit the number and direction of zenith rays for the purpose of optimizing speed/accuracy tradeoffs: the more rays, the slower and more accurate the calculation; the fewer rays, the faster and less accurate the calculation. Control of both the direction and number of zenith rays and azimuth rays is possible, although azimuth angles must be evenly-spaced.

The fast mode provides a fast Superposition dose calculation with a speedup factor of ~2.5 at the cost of small loss in accuracy, compared to the "standard" superposition calculation.

Specifically, for certain clinical situations, the fast Superposition dose is less accurate than the Superposition dose by 1-2%. Also, the fast Superposition monitor unit (MU) is less accurate than Superposition MU by 1-2%.

Acceleration/Speed up of Superposition Algorithm

Even using methods such as kernel inversion, the dose deposition stage of superposition is still the most time-intensive aspect of typical volumetric dose calculations, often requiring thousands of point calculations. As a result XiO has implemented in the superposition algorithm optimization techniques that are able to speed the algorithm while preserving its accuracy.

Required points are determined by back projecting the points in the patient-coordinate calculation volume to corresponding points in the beam-coordinate volume in which the superposition calculations are done. A margin is added so that when the beam coordinate volume is interpolated to the calculation volume there will be no "edge effects." A calculation mask defined according to the beam-coordinate calculation volume (used also in Multigrid) is created that determines which points in beam coordinates must be calculated.

The calculation mask is a point-by-point representation of all dose points that could be calculated for a beam, patient, and calculation volume. It contains information about which points must be calculated using superposition and which azimuth angles will be ray-traced for each point.

Multigrid Dose Calculations

The purpose of multigrid dose calculation is to limit the number of points calculated with superposition, because superposition calculations are expensive. The idea is simple: only dose is calculatd using superposition at points where it is necessary, and at all other points interpolation is used to get a reasonable estimate of dose.

Dose points that are in beam edges or tissue density gradients must be calculated at a higher resolution using more superposition calculations.

Doses will change greatly at beam edges due to entry dose build up and exit dose build down and penumbra. Dose also changes greatly at tissue density gradients due to "secondary" build up (less dense to more dense tissue) or build down (more dense to less dense tissue).

Regions of high resolution calculations are represented as binary "edges" in the calculation mask. This is the same calculation mask used to limit the actual points calculated, and to control which azimuth rays are traced for each point. Multigrid edge detection occurs after back projection from calculation volume to beam coordinates and before azimuth angle masking.

Thus, the number of points calculated depends not only on the beam geometry, but on the nature of the tissue for which the calculations are done.

Implementation of the FFT Convolution and Superposition

Accurate three-dimensional dose computations for radiation therapy using convolution/superposition algorithms have been investigated by many researchers.4],5],7],8] Unlike measurement-based methods, the convolution/superposition method can account for

· the transport of primary and secondary radiation inside the patient.
· the variation of beam intensity across the patient surface.
· the effects of tissue inhomogeneities on the dose.
· irregular blocked or multileaf (MLC) shaped fields.

The Monte Carlo kernels upon which FFT convolution and Superposition algorithms are based - it is calculated by simulating the effects of photons interacting with water in a spherical 60 cm radius "water tank".[2] From the interaction point, the paths of primary electrons and scattered photons are traced until the energy they possess is deposited in the tank or the photons leave the volume. For millions of these "histories," the energy deposition statistics are tabulated, producing probability distributions for energy released at the interaction point and distributed in the water sphere.

The Fast-Fourier Transform (FFT) Convolution algorithm and the Multigrid Superposition[1] algorithm are similar in that they both compute the dose by convolving the total energy released in the patient with Monte Carlo-generated energy deposition kernels computed by Mackie et al.[2] The major differences are that FFT Convolution does not calculate dose as accurately as Multigrid Superposition in the presence of tissue inhomogeneities.

TERMA CALCULATION - CENTRAL-AXIS AND OFF-AXIS ACCELERATOR SPECTRA

The accelerator spectrum is the primary physical data required in convolution and superposition calculations. The machine spectrum cannot be directly measured, so accurate determination of the spectrum requires Monte Carlo simulation and extremely detailed and accurate specification of the physical properties of the accelerator treatment head, as has been described by several authors.9],10],4],6]

The mean energy of the spectrum has been shown by Zhu and van Dyk[11] to be a good predictor of the overall effect of spectrum changes. The system displays the mean energy of the spectrum to aid in the commissioning process and spectrum evaluation. Heuristics for spectrum modification have been formulated on the system based on published results of spectrum behavior.6],9],10] Electron contamination and corrections for kernel-hardening and the parallel-kernel approximation have been implemented, which reduce the complexity and degree of spectrum modification required on the TPS. Implementation of these features also leads to a reduced magnitude in the apparent spectrum changes with field-size required to obtain agreement with measured data.[12]

The XiO system is shipped with default spectra for a variety of machine energies and accelerator types. Off-axis spectra have been created for all the default spectra using published changes in average energy as a function of off-axis distance[9]

Fanned Beam Grid

The FFT Convolution and Superposition algorithms create a fanned beam grid up to a distance of 10 cm outside the field in the lateral direction and over the entire patient extent in the depth direction. It will be roughly equal to the smallest dimension of the calculation volume grid at the isocenter plane. The radiological path for each point in the beam grid is found by performing a ray-trace through the patient electron density matrix to the fan line/depth line intersection points.

Focal-Radiation (Finite Source size) Penumbra

The fluence that has passed through the collimator jaws is set equal to the collimator transmission defined for the FFT convolution and Superposition algorithms. This fluence is then convolved in both the x and y directions with a Gaussian source distribution function to model the penumbra due to the finite source size. The sigma (full-width at half maximum) of the Gaussian distribution function, defined by the user, is scaled to the isocenter plane. For open fields, the fluence is convolved with the user-defined collimator sigma. For block or MLC shaped fields, the fluence is first attenuated by the block or MLC transmission, respectively. The resultant fluence is then convolved with the user-defined block or MLC sigma. Since the transmission through these devices is very small, spectral hardening through them is usually not performed.

Extra-Focal Radiation (Head Scatter)

Several authors10],13-17] have shown that diffuse scatter generated from the flattening filter can contribute up to 12 percent of the total fluence. Whereas the focal radiation or primary fluence is independent of the field size, the extra-focal radiation (head scatter) fluence contribution is strongly dependent on field size, and can be described as an aperture effect.

The main effect of the head scatter contribution on dose distributions is the contribution to dose outside the field edge. The head scatter contribution is very small and can be described as the convolution of the primary fluence with a very broad Gaussian. Measurements by Jaffray et al[13] and Chaney and Cullip[10] indicate that the source distribution is approximately Gaussian and produces a distribution with full-width half-maximum of 2.0 to 2.4 cm at a plane positioned near the flattening filter.

The head scatter fluence contribution is found by convolving the primary fluence with a Gaussian source distribution, the sigma of which is specified by the user. The magnitude of head scatter fluence relative to primary fluence is approximated to be 12 percent for a 35 x 35 cm2 field and 0 percent for a 4 x 4 cm2 field, with a linear relation between these fields.

Although measurements of dose output in air (S c ) provide information about the ratio of head-scattered to primary fluence as a function of field size, these values cannot be used to derive the head scatter contribution since they may contain backscatter into the monitor chamber in linear accelerator.

The head scatter fluence distribution is scaled relative to the primary fluence so that the relative fluence (unattenuated by the jaws) equals the above factor. The head scatter fluence is then added to the primary fluence outside the jaws only.

Modification of the Fluence through Treatment Aids

The energy fluence of each energy component is independently attenuated through wedges and blocks. Therefore a large fraction of the wedge factor will be predicted by the FFT convolution and superposition algorithms without the need for measured data. However, the convolution and superposition algorithms do not model the effects of scatter from treatment aids. Therefore, a wedge factor correction factor is applied to doses of wedged fields equal to the ratio of the measured wedge factor to the calculated wedge factor. A correction factor is generated by the software for each field size for which a wedge factor has been defined. The wedge factor correction factor used for a setup containing a block is the factor interpolated for the blocked equivalent square. Fluence through blocks is modified by the block thickness only.

Dose Deposition Calculation - Generation of Polyenergetic Kernels

Polyenergetic kernels are formed by TERMA weighting monoenergetic kernels using the surface spectrum for open fields and the surface hardened spectrum along the field center fan under the wedge for wedged fields. This weighting method correctly accounts for the variation in attenuation coefficient with energy but ignores kernel hardening with depth, as described by Hoban et al.[18] Papanikolaou et al,[19] Metcalfe et al,[20] and Miften et al.[21] Hoban states that TERMA-weighting is required to account for the variation in attenuation coefficient with energy, and suggests that the effect of not using TERMA-weighting is similar to the effect of using a spectrum which is too hard or too soft. Comparisons of calculated and Monte Carlo-generated depth doses for a 10MV beam indicate that absolute dose errors are largest at the depth farthest from the depth at which the spectrum is composed, but differences in relative doses using a 0 cm, 20 cm or 40 cm spectrum were found to be negligible, as is consistent with the findings of Metcalfe. A kernel-hardening correction is applied on the XiO system to correct for changes of the spectrum with depth, as described below.

Kernel-Hardening Corrections

A physics based kernel-hardening correction factor is applied to both the primary and scatter TERMA, according to the method proposed by Hoban et al.[18] Hoban shows that errors in depth doses of a few percent can be introduced if corrections are not made for the hardening of the kernels, even if the TERMA is hardened accurately as emphasized by Metcalfe. The correction factor is determined based on the collision KERMA-to-TERMA variation with depth. This correction accounts for the increase in the linear energy absorption coefficient to linear attenuation coefficient ratio as the beam hardens.

Parallel-Kernel Approximation Corrections

Convolution algorithms require an invariant kernel, meaning that the kernel orientation for all interaction points must be assumed to be parallel to the beam′s central axis. In reality, kernels should be oriented along the direction of the beam divergence rather than in a parallel orientation. This is a well-known approximation that has been mentioned in many publications.[[3],[4],[5],19]

Kernel transformations can be applied to properly model kernel tilting. However, this rigorous, computationally expensive model can only be used in superposition calculations, since it creates a variant kernel.

Sharpe and Battista[3] and Miften et al,[21] have shown that the parallel-kernel approximation can lead to large dose errors (greater than 3 percent) on the central axis, off-axis, and outside the field edge.

Discrepancies are greater for larger fields and smaller SSDs (that is, in cases where larger beam divergence is present). The nature of the discrepancy is to cause the calculated doses to be more penetrating than the measured doses along the central axis, with diminishing error with off-axis distance. Dose in the penumbra and outside the field edge is likewise underestimated since the kernels are oriented parallel to the beam′s central axis rather than along the divergent beam edge.

A correction for the discrepancies introduced by this approximation along the central axis has been proposed by Papanikolaou et al.[19] The correction involves applying the inverse square correction to the dose after convolution/superposition rather than applying the inverse square correction to the fluence, as is theoretically correct. This correction leads to greatly improved agreement on the central axis, since it causes the calculated doses to be less penetrating, resulting in a scaling to approximate the correct dose on the central-axis. A TPS may use the Papanikolaou correction factor in all dose calculations to correct for errors due to the parallel-kernel approximation. Testing has revealed that the effects of the correction are well-behaved and predictable. However, the corrections are not theoretically rigorous and some discrepancies are still likely to be present in some setups.

Electron Contamination

Since convolution/superposition calculations describe the transport of electrons and photons set in motion by incident photons only, electron contamination should be modeled independently. Electron contamination has been shown by many authors[[23],[24],[25],[26],[27]] to contribute significant dose at the surface. Therefore incorporation of an electron contamination component is extremely important. The electron contamination component is approximately exponential, strongly energy-dependent, field-size dependent, and accelerator-dependent (Sjogren shows a 20 percent discrepancy at 2 cm for a 25 x 25 cm, 20MV microtron). The XiO system automatically determines the electron contamination contribution as a function of field size for your accelerator using entered measured depth doses.

The software computes electron contamination based on the following steps:
· Photon dose is calculated without electron contamination in the same setup for which depth doses have been defined.
· Calculated photon and measured doses are normalized to 10 cm (beyond any electron contamination).
· Calculated dose is subtracted from measured dose and the difference in dose is attributed to electron contamination.
· The resultant electron contamination dose is "unnormalized", and the inverse square dependence is removed, to obtain a contribution which does not contain inverse square dependence, and which is independent of the normalization point.

During validation, the software computes the electron contamination component at field sizes for which measured depth doses have been defined.

The electron contamination dose contribution is computed as a function of depth. An electron contamination dose array is created in fan coordinates by assuming that the depth dependence applies to each fan, starting at the fan entry point. The doses are interpolated from the fanned grid onto a Cartesian grid in the beam orientation, and then added to the dose computed by the FFT or Superposition. Finally, this total dose is interpolated onto the user-requested rectilinear grid.

The behavior of the electron contamination component in the presence of treatment aids is extremely complex and dependent on numerous parameters. A comprehensive electron contamination model accounting for behavior with treatment aids has not yet been implemented. Currently, the electron component is added to all portions of the beam, independent of the presence of treatment aids, on the basis that electrons absorbed by the material are replaced by electron generated in the material. This is approximately the behavior exhibited in Clarkson dose calculations, since measured depth doses or TARs/TPRs include electron contamination and no additional component is added or removed in the presence of treatment aids.

For physical (fixed or motorized) wedges, the amount of electron contamination can be adjusted on per-wedge basis.

Tissue inhomogeneities and Superposition

Tissue inhomogeneities will cause the greatest shape change of the energy deposition kernels and of the resultant dose distribution compared with those predicted in an all water patient. Therefore, an important challenge is to account for tissue inhomogeneities in dose calculation algorithms.

In the XiO algorithms, tissue which lay along the primary radiation ray lines directly influence the TERMA at each point in the patient. All XiO dose calculation algorithms account correctly for this variation. The primary beam penetration is calculated by ray tracing through voxel densities in the 3D volume along divergent beam ray path.

Reliance only on this simple inhomogeneity correction of the primary radiation is inadequate in some clinical applications of high energy photon therapy.

Unlike the FFT convolution algorithm, the superposition algorithm energy deposition kernels can be modified to account for variations in electron density. The density scaling method based on O′Connor′s theorem[28] is used to distort the kernels by finding the average density along the straight-line path between the interaction site and dose deposition site. Density scaling is a good approximation for scattered photons because the photons travel in straight lines and the mass attenuation coefficient scales linearly with material density (assuming that the atomic number remains unchanged).

Electron 3D Pencil Beam Algorithm

The electron pencil beam algorithm as developed by Hogstrom in 1981[29] is used to predict electron dose distributions in the presence of inhomogeneous media for use in radiation treatment planning. An electron beam is modeled as a collection of forward-directed "pencils" at the final plane of collimation at the Source-to-Cone Distance (SCD). The electron pencil beams at subsequent planes are redistributed in a Gaussian distribution due to both scatter occurring above the SCD (air scatter) and scatter in the medium due to Multiple Coulomb Scattering (MCS).

The lateral spread of electrons occurring below the SCD due to scattering occurring above this plane (air scatter) is determined based on the change in the lateral spread of the penumbra of air profiles measured at various source-to-chamber distances. The lateral spread of electrons as they traverse inhomogeneous media is determined according to the Fermi-Eyges theory of thick-target multiple Coulomb scattering (MCS). The air and MCS Gaussian scatter distributions serve as kernels with which the pencils are convolved.

The first step in the calculation is to convolve the initial pencil beam intensity distribution as projected to a given depth with the air Gaussian at that depth. The next step is to determine the central-axis term using measured percent depth dose (PDD) data.

A photon component (approximated as the PDD value at 5mm beyond the practical range for depths above this and the PDD for depths below this) is removed from the measured PDDs to separate the electron and photon components of the beam. The central-axis term is then determined as the dose deposited at depth by an unsmeared pencil beam. This quantity is equal to the electron component of the PDD data divided by the amount of dose scattered away from the central axis which does not get scattered back by surrounding pencils - conceptually, this is effectively de-convolving the measured PDD data.

Each point in the air-convolved distribution is multiplied by the central-axis term for the point′s corresponding effective depth. Each point is also multiplied by an inverse square factor accounting for intensity variations between the distance between the source and the calculation point and the distance of the central axis SSD plus the calculation point′s effective depth.

The next step is to convolve the points with the position-dependent MCS Gaussian. The lateral distribution of photons is determined by multiplying the photon component of the PDD data by a penumbra term dependent on the lateral spread of the air Gaussian. Finally, the electron and photon distributions are added back together in order to obtain the final dose distribution.

DISCUSSION

It is much useful to have different dose calculation algorithm so that it can be invoked as per planning needs.

Significant computation time savings can be gained by performing FFT convolution calculations rather than superposition calculations. However, FFT convolution calculations require the approximation that the energy deposition kernels are invariant over the patient. There are several cases where this assumption is not strictly valid, although in most cases correction factors have been applied to overcome some of these limitations. It is to be noted that many of these approximations such as the polyenergetic and parallel kernel approximations are used in the multigrid superposition algorithm to save computation time.

The 3-D algorithms such as FFT Convolution/Superposition calculate dose to points in the calculation volume based on the scatter contributions from adjacent points. When a 3-D dose calculation is desired, but there are too few patient slices, it is common practice to duplicate cross-sections on a 2-D studyset until sufficient patient data is available. The resulting studyset is sometimes not extended beyond the field size being calculated, doses on the outermost patient crosssections are too low (due to lack of contribution scatter). In order to avoid this problem, the patient studyset ought to be simply extended beyond the field size being calculated (or preferably to the extent of the dose calculation fan lines) for accurate 3-D dose calculations.

Most treatment-planning systems employ pencil-beam algorithms. Dose distributions of photon beams with finite field sizes are calculated by a 2-dimensional convolution of these kernels. It is difficult to modify pencil beams for patient shape and heterogeneity.

For expedience, empirical scaling methods are employed. Because the energy transport of secondary electrons, which may travel up to a few centimeters in water, is not taken into account by scaling methods, dose calculations with pencil-beam algorithms lead to errors near surfaces and inhomogeneities. In addition, Z -variations of different kinds of tissues or materials are not taken into account.

Many of the approximations in the present formulation of the pencil beam algorithm for electron beams are based on analysis of the trade-offs between greater accuracy and such issues as increases in the computation time, greater amounts of input data which must be measured and maintained by users, and greater software complexity.

The most well-known shortcoming of the pencil beam algorithm is its inability to accurately model the effects of inhomogeneities which vary in the lateral direction. The pencil beam formulation is based on the assumption that any inhomogeneity is infinite in lateral extent. Hogstrom (1991)[29] made the following statement regarding the accuracy limitations of the pencil beam model due to this assumption: "... Pencil-beam algorithms have been proved useful for broadbeam electron dose calculations in radiotherapy. However, all pencil-beam models encounter significant limitations (i.e., lateral discontinuities as well as those in depth) when applied to nonslab geometries. Recently Lax[30],31] studied the limitations of the semi-infinite slab approximation for inhomogeneity corrections in electron-beam dose planning. The approximation was applied to both the Gaussian model and a generalized Gaussian model30],31] of the pencil beam. Results show that substantial errors (10-20%) are made in that portion of the dose distribution located on the central axis behind the air cavity when a long and narrow air cavity is placed in the phantom. This error results from the basic assumption that the different materials along the central axis of each pencil beam are infinite in their lateral extent for the purpose of determining the dose to the patient from that pencil beam. ..."

Brahme (1985)[32],33] warns, "... due to the complexity of the interactions of high energy electrons with matter, quite severe errors in the calculated dose distributions may result, particularly when the often considerable limitations of the programs are not fully realized by the user."

Hogstrom (1985)[34] has also stated some of the problems in specific geometries. Bone heterogeneities are frequently found near the skin surface, and have a significant influence on the resulting dose distribution, particularly in treatment of head and neck tumors. The influence of air cavities is also significant, particularly the nasal passages and sinuses in the head, and the trachea in neck irradiations. With respect to the former, Hogstrom et.al. (1984)[35] have made measurements in tissue substitute phantoms which were designed from patient CT scans. Each phantom was made with a constant cross section (i.e., 2-D phantom) because it was technically feasible and that matched the anatomy which the computer code had assumed. These measurements showed that under ideal circumstances the algorithm was accurate to within ±4% in high dose regions of the beam and that isodoses above 10% agree to within ±4 mm in the region of falloff. When air cavities were long and parallel to the beam, errors as great as 13% were observed. These results were consistent with what was expected, that the existing pencil beam algorithm would become less valid for heterogeneities deep to the surface and/or having edges long and parallel to the beam.

By modeling the effects of range straggling and large angle single scattering, Brahme demonstrated significant improvement in the calculated dose distributions.

Various approaches have been taken to include additional interactions ignored in the Fermi-Eyges theory in order to more accurately model the lateral distribution. Shiu and Hogstrom, for example, modeled the variable energy spectrum of electrons over depth in their redefinition of algorithm of 1990, which yields kernels which deviate from the Gaussian distribution predicted by the Fermi-Eyges theory, and produces significant improvements in accuracy.

As for the external photon beam dose calculation model such as FFT convolution and superposition algorithms, they use fundamental physics describing photon and electron interactions and transport. Therefore, the algorithms provide more accurate dose calculation with less parametrization than measurement- based algorithms.

References

Copyright 2005 - Journal of Cancer Research and Therapeutics