Search tips
Search criteria 


Logo of brjradiolSubmitSubscribeAboutBJR
Br J Radiol. January 2017; 90(1069): 20160426.
PMCID: PMC5605017

A kernel-based dose calculation algorithm for kV photon beams with explicit handling of energy and material dependencies



Mimicking state-of-the-art patient radiotherapy with high-precision irradiators for small animals is expected to advance the understanding of dose–effect relationships and radiobiology in general. We work on the implementation of intensity-modulated radiotherapy-like irradiation schemes for small animals. As a first step, we present a fast analytical dose calculation algorithm for keV photon beams.


We follow a superposition–convolution approach adapted to kV X-rays, based on previous work for microbeam therapy. We assume local energy deposition at the photon interaction point due to the short electron ranges in tissue. This allows us to separate the dose calculation into locally absorbed primary dose and the scatter contribution, calculated in a point kernel approach. We validate our dose model against Geant4 Monte Carlo (MC) simulations and compare the results to Muriplan (XStrahl Ltd, Camberley, UK).


For field sizes of (1 mm)2 to (1 cm)2 in water, the depth dose curves show a mean disagreement of 1.7% to MC simulations, with the largest deviations in the entrance region (4%) and at large depths (5% at 7 cm). Larger discrepancies are observed at water-to-bone boundaries, in bone and at the beam edges in slab phantoms and a mouse brain. Calculation times are in the order of 5 s for a single beam.


The algorithm shows good agreement with MC simulations in an initial validation. It has the potential to become an alternative to full MC dose calculation.

Advances in knowledge:

The presented algorithm demonstrates the potential of kernel-based dose calculation for kV photon beams. It will be valuable in intensity-modulated radiotherapy and inverse treatment planning for high precision small-animal radiotherapy.


Mimicking the spatial and temporal dose distributions delivered in state-of-the-art patient treatment in a pre-clinical setting is expected to advance dose–effect studies and radiobiological investigations.1 To this effect, a number of research groups have developed precision radiotherapy systems for small animals.27 In contrast to conventional pre-clinical irradiations, these systems allow positioning and irradiation of tissues of interest with submillimetre accuracy. Common features of high-precision irradiation devices consist of a 225-kVp X-ray tube mounted on a rotating gantry, a precise collimation system down to 0.5-mm beam size, a robotic couch and on-board cone beam CT (CBCT) and planar imaging.

From a hardware point of view, these irradiators are capable of delivering complex dose distributions. However, sophisticated dose calculation and treatment planning systems (TPSs) are still required to fully exploit their potential.8 Clinically available tools designed for patient treatment are not suitable for small-animal radiotherapy due to the different energy regime (kV instead of MV) and much smaller field sizes.9 In addition, the workflow in pre-clinical irradiations differs from patient treatment. The anaesthetised animal is imaged, planned and treated subsequently in one anaesthesia setting of approximately 15 min, which requires especially short treatment planning times.

The introduction of pre-clinical precision irradiators has sparked the development of dedicated dose calculation methods for treatment planning tools.9 Different approaches were pursued, from a measurement-based method10 and a superposition-convolution algorithm1113 to Monte Carlo (MC)-based dose calculation algorithms.1416 The available planning tools for pre-clinical radiotherapy currently support forward-planned dose delivery techniques with open treatment fields from multiple beam directions or arcs, with the recent addition of beam-on time optimization.17

We intend to further advance pre-clinical radiotherapy by developing intensity-modulated dose delivery techniques based on inverse planning strategies. The small-animal-specific workflow combined with inverse optimization requires an accurate and fast dose calculation. Although MC-based dose calculation algorithms are the gold standard in terms of accuracy, they generally suffer from long calculation times. Kernel-based methods permit fast dose calculations but introduce inherent uncertainties at tissue boundaries, and care has to be taken to correctly model the kV photon interactions with matter. In this work, we present a superposition–convolution dose engine, with explicit handling of energy and material dependencies.


Dose calculation

When comparing photon beam interactions of kV and MV energy spectra, three main differences have to be accounted for. First, the importance of the photoelectric effect for kV energies is increased. The cross-section of the photoelectric effect scales with the atomic number as Z3−4, introducing a strong material dependence of the absorption coefficients. Therefore, knowledge of the electron density alone is no longer sufficient, but information on the tissue composition is also required. Additionally, the cross-section depends strongly on the photon energy. Second, the angular distribution of the differential cross-section for Compton scattering changes. Instead of being strongly forward-peaked, the angular distribution of scattering angles is more isotropic. Finally, the range of secondary electrons in water is much shorter for kV energies, in the order of a few micrometres to a few hundred micrometres. Thus, energy is transported through matter by scattered photons, rather than secondary electrons, and the build-up effect is small.

These peculiarities have to be accounted for and modelled in any dose calculation algorithm for kV photons. Our approach to do so is based on a superposition–convolution algorithm. Superposition–convolution is a two-step process. First, the local photon fluence at any point x is determined via the total energy released per unit mass (TERMA)


with the initial photon fluence ϕ0, the energy E, the attenuation coefficient μ, the mass density ρmass and the path of the photons γ. The dose is then calculated as a convolution of the TERMA with a dose deposition kernel.

Dose model

Our dose model is based on previous work on microbeam radiation therapy.18,19 In microbeam radiation therapy, targets are irradiated with parallel kV photon microbeams with a width of a few micrometres. In this approach,18 interactions of photons and electrons with matter are treated separately due to their different ranges in tissue. The dose deposition kernel is split into a photon kernel and an electron kernel. For small animal radiotherapy, this can be further simplified. Electrons with energy of 80 keV, the mean energy of the commercially available systems, have a range of approximately 100 μm in water. This is of the same order of magnitude as the voxel size. We therefore assume local energy deposition at the primary photon interaction site.

For simplicity, we briefly describe the algorithm assuming monoenergetic beams in homogeneous media. The dose deposition kernel can be separated into energy deposited locally by the primary interaction and a scatter kernel. At the primary interaction point, the average fraction


of the photon energy is transferred to the medium.19 qP, qC and qR denote the probabilities for the photoelectric effect, Compton and Rayleigh scattering, respectively, defined as the ratio of absorption coefficients μ, qinteraction=μinteraction/μtotal. pC is the energy fraction transferred to Compton electrons.19 The primary dose can thus be expressed as


All further energy deposition by scattered particles is treated in a point kernel approach. The scatter dose contribution is calculated as a convolution of TERMA with a pre-calculated dose deposition kernel (refer the Kernels section)


The total dose deposited is the sum of primary and scatter contributions Dtotal = Dprimary + Dscatter.

To account for the energy and material dependencies of the absorbed dose, the calculations outlined above are performed for a discrete set of monoenergetic photons and homogeneous materials. The absorbed dose in inhomogeneous media irradiated by polychromatic beams is calculated as an energy-weighted sum of these contributions:


where wE is the weight of the energy bin in the spectrum with wE=1, and NE and NM are the numbers of energy and material components.

Handling of energy and material dependencies

As discussed in the Dose model section, the dose is calculated as a weighted sum of discrete energy and material contributions. The number of chosen calculation points balances accuracy and runtime. We currently employ 6 photon energies and 10 different materials.

The energy sampling points are spaced equally in the cumulative energy-weighted spectrum. The weights wE are given by the spectral integral of the energy bin and satisfy wE=1. The energy values and the respective weights for the small animal radiation research platform (SARRP; XStrahl Ltd., Camberley, UK) are given in Table 1.

Table 1.
Sampling points of the energy spectrum and their respective weighting factors wE

The dose calculation requires knowledge of the absorption coefficients of the irradiated materials, which depend on the tissue composition in addition to the density for kV dose calculation. The tissue parameters are extracted from the Hounsfield units (HUs) of CT according to the method described by Schneider et al.20 It has to be noted that the tissue parameters are based on human tissues, since very limited data on small-animal tissue compositions is available. These continuous tissue parameters are directly used in the calculation of the TERMA and the primary dose through analytical approximations of the absorption coefficients.18,19

The scatter dose, on the other hand, is calculated for a set of 10 materials. For each voxel, we determine the base material closest to the voxel material and use the appropriate kernel in the convolution. The 10 base materials are given in Table 2. They do not correspond to common tissues but are chosen heuristically based on the absorption coefficients, as follows: The material-dependent part of the absorption coefficient for the photoelectric effect has an approximately exponential relationship to that of the Compton scattering. In order to limit the variation in both absorption coefficients within each material segment and thus limit the dose calculation error, we use an exponential segmentation


where α = 2.5 cm3g−1. Each material is assigned to a material segment M{1,2,...,10} based on the material-dependent part of μC:


with the Avogadro constant NA, the elemental weights wi, the atomic number and mass Zi and Ai, and the mass density ρ.

Table 2.
List of materials (M) for which the scatter dose kernels are pre-calculated


We generated the scatter dose deposition kernels with the MC engine Geant421,22 version 10.01p01, using the PENELOPE electromagnetic physics list. For all 60 energy and material combinations, 5 × 108 to 109 particle histories were simulated in homogeneous cubes with monoenergetic photon beams. A biasing region of one-eighth of the voxel size was defined at the centre of the cube, and a mono-energetic, unidirectional photon point source was placed on the border of this region. All primary photons were forced to interact within this biasing region. The dose was scored into a cube with 2513 voxels with a size of (0.275 mm)3. Energy deposited by secondary electrons originating from the first interaction was not scored, as this is accounted for in the primary dose.

The geometry of the pre-clinical irradiators, combined with the small treatment fields used in small animal radiotherapy, leads to small beam divergences. For the SARRP system, the maximal beam divergence of a (1 cm)2 field is 1.15°. In addition, the angular dependence of the Compton scattering cross-section leads to almost isotropic scatter kernels. We therefore neglect the beam divergence within a single beam. An exploratory analysis showed that even with kernels tilted by 45° against the beam direction, the errors in the dose distribution are small. For example, for a (5 mm)2 irradiation field in water, a rotation of the kernels by 45° led to local deviations <0.5% in both the depth dose curve and the lateral profile. Therefore, we do not rotate the kernels to the beam direction. Instead, we determine the main axis closest to the beam direction and flip the kernel by 0°, 90°, 180° or 270° accordingly. The TERMA for all beams within a quadrant is summed, and the scatter dose calculated once for each quadrant. This reduces the number of convolutions from Nbeam × 60 to a maximum of 4 × 60. This enables the calculation of complex dose deliveries such as many-beam treatments or deliveries with a variable collimator, and we do not foresee any issues with the implementation of conformal arcs in the future.


In order to evaluate the performance of the algorithm, we developed a basic stand-alone TPS, which calculates the dose for a given CT and treatment plan as illustrated in Figure 1.

Figure 1.
Workflow of the developed dose calculation algorithm.

In a first step, the CT Hounsfield units (HUs) are converted to material properties (refer the Handling of energy and material dependencies section). The TERMA is pre-calculated for all beams and energy contributions using a Siddon ray tracer,23 which calculates the photon attenuation in depth by ray tracing from a point source to a grid of target points behind the CT volume.24 The actual dose calculation is then split into primary and scattered dose as described above.

Currently, the dose engine is adapted to the geometry of the SARRP. In our current implementation, we only account for the primary radiation emerging from the focal spot of the radiation source but neglect scatter or leakage radiation from the collimator, similar to the approach described by Granton et al.25 The 3 mm focal spot is approximated by a two-dimensional Gaussian distribution.

The cross-platform compatible implementation of the dose calculation method runs on a multicore central processing unit (CPU) environment. All modules of the software have been implemented for modern parallel processors. The convolutions of TERMA and kernel for all energy and material combinations are calculated in Fourier space and implemented using the fast Fourier transformations of the Intel MKL library.26 To prevent cross-correlations due to circular convolution of non-periodic signals in Fourier space, the input is typically zero padded by the size of the kernel. However, this increases the size of the convolution problem and hence the calculation time significantly. The sharp fall-off of the kV dose kernels limits the aliasing errors on the dose distribution. We therefore reduced the zero padding to half the size of the kernels.

Monte Carlo reference dose

MC simulations are generally accepted as the gold standard in terms of accurate dose calculation. We therefore compared our dose engine with MC simulations to validate the calculation accuracy. As for the kernel calculation (refer the Kernels section), we used Geant4 toolkit version 10.01p01,21,22 with the PENELOPE electromagnetic physics description. The same simplified geometry of the SARRP system used in the superposition–convolution dose engine is implemented. The photon energy is sampled from the spectrum, without the discretization needed in the presented dose engine. The same beam parameters are used as input for the simulation and our kernel-based dose calculation, and the conversion method of HUs to material properties are identical.


In addition to MC simulations, we compared our dose engine to Muriplan 2.0.1 (XStrahl Ltd, Camberley, UK), the TPS supplied with the SARRP. Muriplan was developed at Johns Hopkins University, Baltimore, MD. It uses a superposition–convolution dose engine implemented on a graphics processing unit (GPU), based on their MV dose calculation algorithm.1113 The entire small animal treatment is controlled from its 3DSlicer27-based user interface.28 The acquired CBCT is transformed into material properties by defining five discrete windows for the materials air, lung, fat, tissue and bone. For phantom experiments, presets for air, cork, water, graphite and aluminium are available. The dose engine uses look-up tables for the absorption coefficients of these materials, according to the National Institute of Standards and Technology database.29 The tissue density and tissue composition are taken from the International Commission on Radiation Units and Measurements (ICRU) Report 4430 (T Kanesalingam, XStrahl, 08/2016, personal communication). The energy dependence is integrated into the TERMA calculations, using 21 energy bins.11 During the installation of the SARRP system on site, Muriplan is commissioned using measurement data. To the best of our knowledge, no detailed information on the performance of Muriplan or the dose engine has been published.


We validated our developed dose calculation software by comparison with Geant4 MC simulations (refer the Monte Carlo reference dose section) and compared it with the commercial TPS Muriplan (refer the Muriplan section). We evaluated the dose distributions in a virtual homogeneous water phantom, virtual slab phantoms of soft tissue with bone and lung inserts and those calculated on a mouse CBCT. The voxel size was (0.275 mm)3; the standard voxel size at the SARRP. The dose is reported as dose to medium in all cases.

First, we demonstrate the performance of the new dose engine in a (7.04 cm)3 homogeneous water phantom for three different square treatment fields of (1 mm)2, (5 mm)2 and (10 mm)2. To evaluate the performance in inhomogeneous materials, we inserted a 2 mm bone slab or a 5 mm lung slab into a soft-tissue phantom, positioning the top of the slab at a depth of 5 mm. In both cases, we used a (5 mm)2 square treatment beam.

For our dose comparison, we created virtual CT cubes with 2563 voxels. The HU–tissue association as described by Schneider et al20 resulted in HUwater = 0, HUsofttissue = 43, HUlung = −741 or HUbone = 1542 for the ICRU-44 materials used in Muriplan. The same CTs were imported into Muriplan, and the materials assigned accordingly. In the MC simulations, we created the virtual phantom directly in Geant4 from the material compositions.

Finally, we calculated the dose distributions on a mouse CBCT acquired at the SARRP. It has to be noted that the output of the on-board CBCT of the SARRP are not Hounsfield units but arbitrary units, as the system is not calibrated for absolute attenuation coefficients. Therefore, the uncalibrated CT numbers of mouse CBCT were converted to HUs prior to the dose calculation using piece-wise linear interpolations between the five materials used in Muriplan. The converted CT was used in the new dose engine as well as in the MC simulations. In Muriplan, the standard 5 material workflow was used. We defined an artificial target in the mouse brain and planned a 5-field irradiation with equidistant beams of (3 mm)2 at 0°, 72°, 144°, 216° and 288°.

The runtime of our dose calculation method was measured on a shared memory system configured with two Intel Xeon E5-2697v3 CPUs.


Dose distributions


The dose distributions in water calculated with our kernel-based algorithm, Muriplan, and MC simulations are compared in Figure 2. For three different treatment field sizes ((1 mm)2, (5 mm)2 and (10 mm)2), the depth dose curves and the lateral profiles at a depth of 6 mm are shown. Both our kernel-based algorithm and the MC simulations have not yet been calibrated in terms of machine output and therefore no direct correlation between delivery times and absolute dose is possible. Therefore, the depth dose curves are normalized to the point at 2 cm depth, as this is the standard reference point in kV dosimetry protocols. The lateral profiles are normalized to the integral of the depth dose curve.

Figure 2.
Comparison of the analytical dose engine with Muriplan and Monte Carlo (MC) simulations in water. Depth dose curves (a,c,e) and lateral profiles at a depth of 6 mm (b,d,f) are shown for three different field sizes [(1 mm)2, a,b; (5 mm) ...

For our dose engine (shown in red), the mean agreement of the calculated depth dose curves with the simulations is 1.7% for all collimator sizes in water. Larger deviations of up to 4% are observed in the first voxel, i.e. in a 0.275 mm deep entrance region. In addition, the deviation between the new analytical algorithm and MC simulations increases with depth but does not exceed 5%. The lateral beam profiles were evaluated at different depths (a depth of 6 mm was chosen for demonstration purposes in Figure 2). Within the beam, the simulations and our superposition–convolution algorithm show a mean agreement of 3% for all beams, apart from the beam edge. The out-of-field doses, which are defined by values <10% of the maximum dose, agree within 12% (mean) for the (5 mm)2 and (10 mm)2 beams, up to a distance of 34 mm, which corresponds to the size of the kernel. For the smallest treatment field, the mean agreement is 22%.

The depth dose curves calculated with Muriplan (shown in green) agree with the simulations within 3% (mean) for the (5 mm)2 and (10 mm)2 beams, and 12% (mean) for the smallest field. The dose fall-off is consistently steeper in Muriplan than in the MC simulations, and this effect is most pronounced for the smallest beam size. In the lateral profiles, the mean agreement is 3% within the beam and 20% (mean) in the out-of-field regions for all but the smallest fields. For the (1 mm)2 field, the in-beam region agrees within 11%, and the out-of-field doses to 43% (mean). It has to be noted, however, that a different head model is used in the MC simulations and Muriplan, which probably accounts for the large differences in the smallest field size.

Figure 3 depicts the depth dose curves in slab phantoms of water with either 2 mm of bone or 5 mm of lung inserted at a depth of 5 mm. For the new dose engine, the mean agreement in the bone case is 1.5%, with the largest deviations up to 7% at the boundary of the bone slab.

Figure 3.
Comparison of the depth dose curves of (5 mm)2 fields, normalized to the value at 2 cm depth, calculated with the analytical dose engine, Muriplan and Monte Carlo (MC) simulations in slab phantoms of soft tissue with bone (a) and lung ...

In Muriplan, the dose in the bone insert is overestimated by up to 43% compared with the MC simulations. For the lung phantom, simulations and the Muriplan results agree within 3% (mean), whereas our analytical calculations underestimate the dose in lung by 9%. However, if the same CT and material conversions as in the analytical dose calculation are used in the MC simulation (Schneider conversion, shown in pink), our new dose engine and simulations agree within 2% (mean). The same tendencies regarding the build-up region, larger depths and the slope of the depth dose curves as in the water phantoms are observed in the slab phantoms.

Mouse cone beam CT

The dose distributions from a five-field irradiation of a mouse brain, obtained with our kernel-based dose engine, MC simulations and Muriplan, are shown in Figure 4a,b,c, respectively. They are normalized to the dose at the centre of the target. The MC simulations were calculated to a statistical uncertainty of 0.5% in the target region. In the bottom row of Figure 4, the dose difference maps between the simulated dose distribution and our dose engine [Figure 4d] and Muriplan [Figure 4e] are shown. In both cases, the largest differences can be observed in bone and at the beam edges. Muriplan overestimates the dose in bone by 53% (mean), whereas our dose engine shows an underdosage of 14% (mean).

Figure 4.
Dose distributions from a five-field (3 mm)2 irradiation calculated on a mouse cone beam CT using the new dose engine (a), Monte Carlo (MC) simulations (b) and Muriplan (c). The absolute dose differences are shown in (d) and (e).

Calculation times

The calculation times depend on the size of CT and the number of treatment beams. In our implementation, the calculation time does not change significantly from four beams onwards, since we perform the scatter dose calculation for four quadrants instead of each beam. As described in the Implementation section, we reduced the zero padding of the scatter dose kernels to reduce the runtime. We found that this influences the results by <10−4%, whereas saving roughly 40% of the calculation time. The calculation times of the single-beam plans on the cubic water CT with 2563 voxels and the five-beam plan on mouse CT are given in Table 3.

Table 3.
Calculation times for different CT and beam configurations for our dose engine (using loss-less full padding and the reduced padding approximation) and Muriplan


In this study, we presented a superposition–convolution-based dose engine for small animal radiotherapy. Owing to the fact that pre-clinical irradiators operate in the kV energy regime instead of MV, the increased importance of the photoelectric effect and the resulting increased material and energy dependencies have to be taken into account.31 We split the calculation into energy deposited by the primary interaction and all subsequent scattering events. The primary dose is calculated analytically based on statistical considerations. The scatter dose is calculated as a convolution of TERMA with the appropriate MC-generated kernel from a set of 60 energy and material combinations.

We implemented the algorithm into a basic stand-alone TPS. The cross-platform compatible software is optimized for modern multicore CPU systems. It has been adapted to the geometry and energy spectrum of the SARRP system but can be adapted to any pre-clinical irradiator.

In comparisons with MC simulations, the depth dose curves in water show a mean disagreement of 1.7% (Figure 2). A larger deviation of up to 4% is observed in the first voxel. This behaviour is expected, since the kernel-based approach does not correctly model the build-up effect and overestimates the dose in the entrance region. For the energies used in pre-clinical research, the electron ranges and hence the build-up region are of the same order of magnitude as the voxel size.

One factor causing the increasing deviation in depth is our energy-sampling approach. As discussed in the Handling of energy and material dependencies section, the TERMA and hence the attenuation in depth is calculated separately for each of the six energy contributions. Beam hardening is therefore implicitly accounted for. However, we only regard six discrete energy values, and a careful analysis of the chosen sampling points could further improve the calculations. In the lateral profiles, large deviations at the field edge can be observed. This is a discretization effect of ray tracing through a voxelized geometry. The out-of-field dose values agree within 12% for all but the smallest beam. Improving the dose calculation accuracy at the beam edges is expected to also increase accuracy distant from the beam. This issue will be addressed in future work.

The dose distributions in water calculated with Muriplan show systematic differences to our MC simulations. It has to be noted that in contrast to the MC simulations and our dose engine, Muriplan has been commissioned to measurements. The differences in slope of the depth dose curves are most pronounced for the 1 mm field. We assume this is due to differences in the head model. For the smallest beam sizes, the focal spot distribution becomes important.32 This highlights the importance of a thorough validation of the treatment dose calculation by comparison to measurements, especially for the smallest beam sizes. Future work includes the commissioning of our dose engine and an adaptation of the head model where necessary.

In the slab phantoms, larger deviations between our dose calculation and MC simulations are observed at tissue interfaces (Figure 3). This is an inherent issue of kernel-based approaches. Our approach of choosing the appropriate kernel for each voxel based on the voxel material assumes local homogeneity. Kernel stretching is not straightforwardly adapted from MV dose calculations to the kV energy regime due to the increased importance of the photoelectric effect. This results in a smearing of the scatter dose at sharp tissue boundaries. The primary dose is not affected since there is no material approximation in its calculation.

In addition, discrepancies are observed within the tissue slabs. In bone, Muriplan overestimates the dose by up to 43%, whereas our dose engine underestimates dose to the lung by up to 9%. In these cases, we used the ICRU material definitions of the Muriplan base materials in the MC simulations. Our dose engine, however, transforms HUs to material parameters according to Schneider et al.20 If the same CT and material conversion as in the analytical dose calculation are used in the MC simulation, the mean disagreement in lung reduces to 2%. The overestimation in bone observed with Muriplan, however, cannot be explained by the tissue segmentation. In Muriplan, the kernels are calculated in water and the dose to water is converted to dose to medium using the ratio of mass energy absorption coefficients in water and tissue (T Kanesalingam, XStrahl, 10/2016, personal communication). Although this approach is suitable in MV dose calculations, it introduces errors in the kV energy regime which will be resolved in the future using correction factors for dose to medium (T Kanesalingam, XStrahl, 10/2016, personal communication).

The same tendencies as in water and slab phantoms can be observed in the 5-field irradiation of a mouse brain (Figure 4). Discrepancies between MC and analytically calculated dose distributions are largest within the bone, at the beam edges and at tissue–bone interfaces.

The deviations introduced by the material segmentation are not surprising. Bazalova et al31 suggest that 29 materials are needed to calculate the dose to within 2% for 225-kV photon beams on noisy CT images. A four-tissue segmentation similar to the one used in Muriplan lead to dose differences of 57% in the ribs,31 which, in combination with the discrepancies in the bone slab described above, corresponds well to the overestimation of dose in bone that we observe in the skull for Muriplan. Using eight tissues in the segmentation leads to errors in the order of 7%.31 We currently use 10 base materials for the scatter dose kernels as a compromise between calculation time and accuracy. However, this segmentation does not affect the primary dose calculation where we calculate the absorption coefficients for each voxel based on the local material parameters (refer the Handling of energy and material dependencies section). This reduces the influence of the material segmentation on the dose distributions but requires accurate material definitions. The discrepancies between the doses in lung calculated with different material definitions suggest that improvements in the HUs to material conversion are required. This will be addressed in future work.

We follow the current practice by assigning human tissue parameters to mouse tissues, which will lead to further inaccuracies in small animal treatments.31 Only very limited data on mouse tissue composition is currently available. Once more data on mouse tissues are available, further analysis of the number and choice of materials for the pre-calculated kernels and the HU to material conversion will be necessary.

Our dose engine is implemented on a multicore CPU environment and has been optimized for calculation speed. Calculation times for a (5 mm)2 field single-beam plan in the water phantom are 5.2 s and are 18 s for a five-beam irradiation on mouse CBCT. The dose calculation in Muriplan for the same treatment plans required 5.4s and 18 s, respectively. Halving the size of the zero padding has significantly reduced the calculation times without reducing calculation accuracy. Based on these results, we will implement a smart padding, which chooses the smallest possible zero padding without compromising the accuracy based on the CT and beam plan, to reduce the calculation times further.

Kernel-based approaches generally come with a loss of accuracy compared with full MC simulations but offer faster computation times. The commercial TPS SmARTplan (SmART scientific solutions B.V., Maastricht, Netherlands) employs a MC dose engine. Dose calculation times of 243 s have been reported for an irradiation of the left kidney of a rat with four (1 cm)2 beams.17 In addition, MC dose calculations are generally performed for a specific statistical uncertainty in the target region. Although this avoids an increase in computation times with the number of treatment beams, it results in a loss of statistics in low dose regions. Our dose engine is a fast alternative to general full MC simulations with an acceptable accuracy, and it is promising for the use in an inverse treatment planning environment.


We have presented a dose calculation engine for small-animal radiotherapy, which is based on a superposition–convolution approach adapted to kV photon beams. Good agreement of the calculated dose distributions with MC simulations could be shown in water, in slab phantoms and for a mouse CBCT. Runtimes are in the order of 5 s for a single beam. Future work will focus on the energy and material sampling points, as well as improvements to the head model and the adaptation of input parameters of our dose algorithm to the treatment device. We believe the presented dose calculation algorithm will be valuable in inverse treatment plan optimization, which we identify as one of the main steps to further advance high-precision small-animal radiotherapy.


Research at the ICR is also supported by Cancer Research UK under Programme C33589/A19727. MFF is supported by Cancer Research UK under Programme C33589/A19908.


We acknowledge support from an Engineering and Physical Sciences Research Council strategic equipment grant (EP/N015266/1).


1 . Butterworth KT, , Prise KM, , Verhaegen F.. Small animal image-guided radiotherapy: status, considerations and potential for translational impact. Br J Radiol 2015 ; 88: 20140634. doi: [PMC free article] [PubMed]
2 . Wong J, , Armour E, , Kazanzides P, , Iordachita I, , Tryggestad E, , Deng H., et al. . High-resolution, small animal radiation research platform with x-ray tomographic guidance capabilities. Int J Radiat Oncol Biol Phys 2008 ; 71: 1591–9. doi: [PMC free article] [PubMed]
3 . Zhou H, , Rodriguez M, , van den Haak F, , Nelson G, , Jogani R, , Xu J., et al. . Development of a micro-computed tomography-based image-guided conformal radiotherapy system for small animals. Int J Radiat Oncol Biol Phys 2010 ; 78: 297–305. doi: [PMC free article] [PubMed]
4 . Clarkson R, , Lindsay PE, , Ansell S, , Wilson G, , Jelveh S, , Hill RP., et al. . Characterization of image quality and image-guidance performance of a preclinical microirradiator. Med Phys 2011 ; 38: 845–56. doi: [PubMed]
5 . Stojadinovic S, , Low DA, , Hope AJ, , Vicic M, , Deasy JO, , Cui J., et al. . MicroRT-small animal conformal irradiator. Med Phys 2007 ; 34: 4706–16. doi: [PubMed]
6 . Tillner F, , Thute P, , Löck S, , Dietrich A, , Fursov A, , Haase R., et al. . Precise image-guided irradiation of small animals: a flexible non-profit platform. Phys Med Biol 2016 ; 61: 3084–108. doi: [PubMed]
7 . Verhaegen F, , Granton P, , Tryggestad E.. Small animal radiotherapy research platforms. Phys Med Biol 2011 ; 56: R55–83. doi: [PubMed]
8 . Bazalova M, , Nelson G, , Noll JM, , Graves EE.. Modality comparison for small animal radiotherapy: a simulation study. Med Phys 2014 ; 41: 011710. doi: [PubMed]
9 . Verhaegen F, , van Hoof S, , Granton PV, , Trani D.. A review of treatment planning for precision image-guided photon beam pre-clinical animal radiation studies. Z Med Phys 2014 ; 24: 323–34. doi: [PubMed]
10 . Marco-Rius I, , Wack L, , Tsiamas P, , Tryggestad E, , Berbeco R, , Hesser J., et al. . A fast analytic dose calculation method for arc treatments for kilovoltage small animal irradiators. Phys Med 2013 ; 29: 426–35. doi: [PubMed]
11 . Jacques R, , Smith D, , Tryggestad E, , Wong J.. GPU accelerated real time kV/MV dose computation. Proc. ICCR 2010 ; 2010: 2–5.
12 . Jacques R, , Taylor R, , Wong J, , McNutt T.. Towards real-time radiation therapy: GPU accelerated super-position/convolution. Comput Meth Programs Biomed 2010 ; 98: 285–92. [PubMed]
13 . Jacques R, , Wong J, , Taylor R, , McNutt T.. Real-time dose computation: GPU-accelerated source modeling and superposition/convolution. Med Phys 2011 ; 38: 294–305. doi: [PubMed]
14 . Chow JC, , Leung MK.. Treatment planning for a small animal using Monte Carlo simulation. Med Phys 2007 ; 34: 4810–7. [PubMed]
15 . Graves EE, , Zhou H, , Chatterjee R, , Keall PJ, , Gambhir SS, , Contag CH., et al. . Design and evaluation of a variable aperture collimator for conformal radiotherapy of small animals using a microCT scanner. Med Phys 2007 ; 34: 4359–67. [PubMed]
16 . van Hoof SJ, , Granton PV, , Verhaegen F.. Development and validation of a treatment planning system for small animal radiotherapy: SmART-Plan. Radiother Oncol 2013 ; 109: 361–6. doi: [PubMed]
17 . Balvert M, , van Hoof SJ, , Granton PV, , Trani D, , den Hertog D, , Hoffmann AL., et al. . A framework for inverse planning of beam-on times for 3D small animal radiotherapy using interactive multi-objective optimisation. Phys Med Biol 2015 ; 60: 5681–98. doi: [PubMed]
18 . Bartzsch S.. Microbeam radiation therapy physical and biological aspects of a new cancer therapy and development of a treatment planning system. PhD Thesis. Germany: University of Heidelberg; 2014 .
19 . Bartzsch S, , Oelfke U.. A new concept of pencil beam dose calculation for 40–200 keV photons using analytical dose kernels. Med Phys 2013 ; 40: 111714. doi: [PubMed]
20 . Schneider W, , Bortfeld T, , Schlegel W.. Correlation between CT numbers and tissue parameters needed for Monte Carlo simulations of clinical dose distributions. Phys Med Biol 2000 ; 45: 459–78. doi: [PubMed]
21 . Agostinelli S, , Allison J, , Amako K, , Apostolakis J, , Araujo H, , Arce P., et al. . Geant4 - a simulation toolkit. Nucl. Instrum. Methods Phys. Res., Sect. A. 2003 ; 506: 250–303.
22 . Allison J, , Amako K, , Apostolakis J, , Araujo H, , Dubois P, , Asai M., et al. . Geant4 developments and applications. IEEE Trans Nucl Sci 2006 ; 53: 270–8. doi:
23 . Siddon RL.. Fast calculation of the exact radiological path for a three-dimensional CT array. Med Phys 1985 ; 12: 252–5. doi: [PubMed]
24 . Siggel M, , Ziegenhein P, , Nill S, , Oelfke U.. Boosting runtime-performance of photon pencil beam algorithms for radiotherapy treatment planning. Phys Med 2012 ; 28: 273–80. doi: [PubMed]
25 . Granton PV, , Verhaegen F.. On the use of an analytic source model for dose calculations in precision image-guided small animal radiotherapy. Phys Med Biol 2013 ; 58: 3377–95. doi: [PubMed]
26 . [Cited November 7 2016]. Available from:
27 . Fedorov A, , Beichel R, , Kalpathy-Cramer J, , Finet J, , Fillion-Robin JC, , Pujol S., et al. . 3D Slicer as an image computing platform for the quantitative imaging network. Magn Reson Imaging 2012 ; 30: 1323–41. doi: [PMC free article] [PubMed]
28 . Cho NB, , Kazanzides P.. A Treatment Planning System for the Small Animal Radiation Research Platform (SARRP) based on 3D Slicer. MIDAS J 2012 ; 1–8. [Cited November 7 2016]. Available from:
29 . [Cited November 7 2016]. Available from:
30 . ICRU. ICRU Report 44. Tissue Substitutes in Radiation Dosimetry and Measurement; 1989.
31 . Bazalova M, , Graves EE.. The importance of tissue segmentation for dose calculations for kilovoltage radiation therapy. Med Phys 2011 ; 38: 3039–49. doi: [PubMed]
32 . Tryggestad E, , Armour M, , Iordachita I, , Verhaegen F, , Wong JW.. A comprehensive system for dosimetric commissioning and Monte Carlo validation for the small animal radiation research platform. Phys Med Biol 2009 ; 54: 5341–57. doi: [PMC free article] [PubMed]

Articles from The British Journal of Radiology are provided here courtesy of British Institute of Radiology