A point kernel algorithm for microbeam radiation therapy

Microbeam radiation therapy (MRT) is a treatment approach in radiation therapy where the treatment field is spatially fractionated into arrays of a few tens of micrometre wide planar beams of unusually high peak doses separated by low dose regions of several hundred micrometre width. In preclinical studies, this treatment approach has proven to spare normal tissue more effectively than conventional radiation therapy, while being equally efficient in tumour control. So far dose calculations in MRT, a prerequisite for future clinical applications are based on Monte Carlo simulations. However, they are computationally expensive, since scoring volumes have to be small. In this article a kernel based dose calculation algorithm is presented that splits the calculation into photon and electron mediated energy transport, and performs the calculation of peak and valley doses in typical MRT treatment fields within a few minutes. Kernels are analytically calculated depending on the energy spectrum and material composition. In various homogeneous materials peak, valley doses and microbeam profiles are calculated and compared to Monte Carlo simulations. For a microbeam exposure of an anthropomorphic head phantom calculated dose values are compared to measurements and Monte Carlo calculations. Except for regions close to material interfaces calculated peak dose values match Monte Carlo results within 4% and valley dose values within 8% deviation. No significant differences are observed between profiles calculated by the kernel algorithm and Monte Carlo simulations. Measurements in the head phantom agree within 4% in the peak and within 10% in the valley region. The presented algorithm is attached to the treatment planning platform VIRTUOS. It was and is used for dose calculations in preclinical and pet-clinical trials at the biomedical beamline ID17 of the European synchrotron radiation facility in Grenoble, France.

by the kernel algorithm and Monte Carlo simulations. Measurements in the head phantom agree within 4% in the peak and within 10% in the valley region.
The presented algorithm is attached to the treatment planning platform VIRTUOS. It was and is used for dose calculations in preclinical and petclinical trials at the biomedical beamline ID17 of the European synchrotron radiation facility in Grenoble, France.
Keywords: microbeam radiation therapy, dose calculation, synchrotron radiation, point kernel algorithms (Some figures may appear in colour only in the online journal)

Introduction
Microbeam radiation therapy (MRT) (Slatkin et al 1992) is a preclinical radiotherapy treatment approach that modulates beam intensities on a micrometre scale creating extremely high peak doses that are separated by low dose regions with doses below the tissue tolerance level. Abundant preclinical data show that these irradiation patterns control and even ablate fast growing tumours on the one hand (Dilmanian et al 2002, Bouchet et al 2010, Miura et al 2006 while sparing normal tissue on the other hand (Laissue et al 2001, Serduc et al 2006, 2008, Bouchet et al 2010. This differential effect gives great hope for the treatment of brain tumours in young children where operative removal is impossible or likely to affect the brain development (Laissue et al 2007). Another promising target of microbeam treatment may be lung cancers as recent, yet unpublished studies show. Until now MRT is bound to large 3rd generation synchrotrons being the only sources of almost parallel photon beams with sufficiently high dose rates in the relevant energy range. The European synchrotron radiation facility (ESRF) in Grenoble (France) is one of the leading institutions in MRT research. At its biomedical beamline ID 17 a multislit collimator (MSC) situated 41.7 m downstream from the wiggler source produces peak dose regions of around 50 µm width and a few cm height with 400 µm centre-to-centre peak distance (ctc) (Brauer-Krisch et al 2009). The photon spectrum extends from around 40 to 200 keV with a mean energy of approximately 100 keV (Slatkin et al 1992, Siegbahn et al 2006. The first MRT clinical trials in pet patients were carried out at the ESRF in November 2013.
Accurate dose calculation is an important prerequisite for a future clinical application of MRT as treatment modality. Dose calculations in MRT imply additional challenges as compared to dose calculations in conventional radiation therapy: dose needs to be determined on a micrometre grid and over a dynamic range of several orders of magnitude. Moreover, photon energies in the spectrum of the synchrotron wiggler are far below the typical photon energies in conventional radiotherapy treatment which renders existing fast kernel based dose calculation algorithms inapplicable.
Until now dose calculations in MRT are based on Monte Carlo simulations (Siegbahn et al 2006, Nettelbeck et al 2009, Martínez-Rovira et al 2010. Slatkin et al (1992) calculated MRT doses for monochromatic x-ray beams of 50, 100 and 150 keV in a cylindrical water phantom of 16 cm diameter and 16 cm length representing a simple model of a human head with the EGS4 Monte Carlo code. They compared cylindrical and planar microbeams of different ctc and beam width. Beams were assumed to have no beam divergence. Siegbahn et al (2006) studied microbeams with simulations in Penelope. They analyzed spectral differences between the radiation in peak and valley regions and studied the influence of ctc, beam width and photon energies on the peak to valley dose ratio (PVDR). In their study parallel microbeams were assumed and the dose distribution of a single simulated microbeam was superimposed in order to obtain the dose distribution of microbeam fields. Persisting differences between measured and simulated doses led to the investigation of influences from beam divergence and beamline components on the dose distribution (Nettelbeck et al 2009). Martínez-Rovira et al (2011) succeeded finally with the simulation of the whole beam line and analyzed the phase space of the source in detail. Bartzsch et al (2014) showed that the phase space in the plane of the MSC can be described by a simplified model of parallel beams that accounts for leakage radiation and corrects variations in the beam intensities.
Monte Carlo simulations require a large amount of particle histories for accurate dose calculations on a micrometre grid and therefore the calculation times are too long for treatment planning (Siegbahn et al 2006) and optimization of treatment parameters such as beam direction, energy spectrum (filtering), beam distances etc. Future clinical applications of MRT necessitate faster dose calculation methods. Apart from MRT other applications of keV-range x-rays such as small animal irradiators (Mackie et al 2003, Xing et al 2006, Prajapati 2015 and diagnostic x-ray dose estimates (Kalender 2014, Pawlowski andDing 2014) would benefit from fast dose calculation algorithms.
In conventional radiotherapy with linear accelerators fast dose calculation algorithms were developed on the basis of point scatter kernels (Mackie et al 1985). The superposition or convolution of dose kernels over the energy transfer in primary interactions can be used to deduce dose distribution estimates in patients during therapy planning. Even faster are pencil beam algorithms (Mohan et al 1986) that use radiological depth scaling of pencil beam kernels reducing dose calculation to 2D convolutions (Bortfeld et al 1993). However, these algorithms are based on O'Connors theorem (O'Connor 1957) of dose scaling with electron density. This theorem is valid for the Compton effect only and hence for photon energies where the Compton effect is the prevailing interaction (Alaei et al 1999). In the for MRT interesting energy range photoelectric absorption becomes more important. There have been attempts to adapt kernel based dose calculation algorithms to keV photons. Alaei et al (2000Alaei et al ( , 2009 used an empirically motivated density scaling but got up to 16% deviations and close to material interfaces they observed even higher deviations. Another possibility to obtain fast dose calculation is accelerating existing Monte Carlo based dose calculation algorithms by parallelization on GPUs (Badal andBadano 2009, Jia et al 2010). Jia et al (2012) achieved dose calculation on a realistic geometry given by CT-data in just 17 s. However, they calculated on voxels of 4 mm side lengths, which is far from micrometer sized resolution required in MRT.
Here a convolution based dose calculation algorithm is presented that uses analytic dose kernels described in previous work (Bartzsch and Oelfke 2013). The code is implemented in C++ and requires just a few minutes for a complete 3D dose calculation in MRT. The program is linked to VIRTUOS, a platform for 3D radiotherapy that was developed at the German Cancer Research Centre (DKFZ) in Heidelberg (Bendl et al 1994) and was used for the treatment planning of the first MRT pet patient at the ESRF.

Methods
The two dominant effects at 40 to 200 keV photon energy in photon matter interaction are Compton effect and photoelectric absorption. In these interactions secondary electrons are produced, which get around 10% of the photon energy when produced in Compton scattering or all of the primary photon energy when produced in photoelectric absorption. Secondary electrons lose their energy very rapidly in inelastic scattering and their range in water is less than 10 µm for Compton electrons and less than 300 µm for photoelectrons, respectively. In a microbeam field peak doses are almost exclusively built up by electrons produced in interactions of the primary unscattered photon beam, whereas the valley dose has two constituents, energy absorbed from scattered photons and energy transported to the valley by electrons created in interactions of the primary photon beam in the peak. Energy transport of electrons created in interactions of scattered photons can be neglected since electron ranges are much shorter than scattered photon ranges. Hence in order to calculate the therapeutically important peak and valley doses it is sufficient to consider scattering of electrons produced in primary photon interactions and photon scattering.
Therefore the authors present a dose calculation approach that completely separates electron energy transport on a micrometre scale from photon energy transport on a millimetre scale and the dose is calculated in two steps. Firstly, the primary and scattered photon dose, i.e. the energy per mass element transferred to electrons from primary and scattered photons, are calculated for each millimetre sized voxel (e.g. grid of a CT) in a photon point kernel algorithm. Secondly an electron point kernel algorithm is used to compute the transport of electrons produced by the primary photon beam in the microbeam peaks on a micrometre grid within each voxel. The results of both steps are then combined to calculate dose profiles, peak doses and valley doses.

Treatment of photons
In the photon point kernel algorithm the energy absorption is separated into different scattering orders, which allows an analytic calculation of dose kernels using the differential scattering cross sections of the Compton effect. Bartzsch and Oelfke (2013) present a comprehensive derivation of photon point and pencil beam kernels of photons with a kinetic energy between 40 and 200 keV and describe methods to transform kernels between energies and materials by simple geometrical scaling.
The photon dose is calculated on the millimetre sized voxel grid. If E 0 denotes the energy of the primary photon, q(E 0 , M) the probability that the first interaction in material M is a Compton interaction and p(E 0 ) the average energy transferred to the Compton electron in a Compton scattering event, then the average energy transferred to electrons in the primary photon interaction is and is absorbed within a voxel length from the point where the interaction takes place. The scattered photon dose is calculated in a kernel superposition algorithm. Scatter point kernels K E,M ( r) are pre-calculated for various materials M and energies E on the millimetre grid. Currently implemented are the materials presented in table 1 at the photon energies E 1 = 45, E 2 = 55, E 3 = 65, E 4 = 80, E 5 = 105 and E 6 = 135 keV. This choice of energies proved to yield reliable results for the ESRF ID17 beamline spectrum.
The spectral energy of the primary photon beam is assigned to weighted contributions of the six energies E 1 , . . . , E 6 and for each energy E the interaction strength S E is determined 4 , defined by where the cartesian coordinate system with x, y, and z-axis is defined such that the z-axis points in beam direction. B is the set of points that are within the microbeam treatment field in the x-y-plane. The absorption coefficient μ is obtained from the voxel material (for details see section 2.3), and the integral in the exponent is approximated as a discrete sum on the millimeter sized voxel grid. Each voxel is assigned to the closest matching material M in M = 1, 2, . . . , 10 (table 1). The kernel superposition is performed separated by energy and material and the individual dose contributions are calculated as a convolution of the interaction strength and the photon scatter kernels K E,M , Here D E,M is the dose contribution of material M and energy E and r = (x, y, z) the position vector. S E,M is the interaction strength for the material M, which we define as S E for all voxels that were assigned to M and 0 otherwise. After the calculation of all D E,M (currently 6 · 10 = 60) the total scattered dose is computed as a weighted sum where the weighting takes into account the local absorption coefficient and transferred energy fraction at position r in analogy to equation (1). The primary photon dose, defined as the dose absorbed by the primary unscattered photon beam is obtained from the interaction strength S E and the primary energy absorption fraction in equation (1), V Voxel is the volume of a single voxel and ρ the mass density at r . For broad beams the total dose D tot is the sum of scattered and primary photon dose D tot = D scattered + D primary

Treatment of electrons
Secondary electrons lose their energy by radiative loss (bremsstrahlung) and scattering. The fraction of radiative energy loss is approximately given by (see for example Krieger (2007)), where S is the stopping power, Z the atomic number and E the electron energy. For oxygen atoms ( Z = 8) and E = 150 keV radiative energy loss is responsible for only 0.2% of the total electron energy dissipation and can hence be neglected. Depending on the distance of an electron trajectory from a scatterer atom (impact parameter) and the kinetic energy of the electron elastic, inelastic collision, Coulomb scattering and δelectron production lead to electron deviation and deceleration. For energies below 1 MeV inelastic collisions are the dominant interactions. The collisional stopping power as a function of energy was described in the non-relativistic case by Bethe (1930) Here e is the positron charge, ρ the density, M A the atomic mass, 0 the dielectric constant, I the atomic ionisation energy, and N A the Avogadro constant. From S(E) the electron range can be calculated in the continuous slow down approximation (Podgoršak 2010) In order to derive analytic electron scatter kernels it is useful to simplify the expression for the stopping power in equation (7). Identifying ζ with and η with we can rearrange equation (7) to and write the derivative of ζ as The ionisation energy is in the order of 10 eV and hence η will be between 9.5 and 11.5 for kinetic energies between 40 and 200 keV. When the electron slows down |dζ/dη| will decrease. However, when dζ/dη starts to change significantly the electron has already travelled most of its path. An electron of 150 keV for example has a CSDA range of 282 µm in water. When |dζ/dη| has decreased by 1% the kinetic energy of the electron will be approximately 55 keV and its range has decreased to 51 µm (Berger et al 2005). Hence it has already travelled 80% of its total distance. Therefore it is reasonable to assume that dζ/dη is constant over the distance travelled and to approximate ζ linearly with constants m and n by with constants K and λ. K and λ can be derived from ionisation energies. We have, however, pursued another approach and fitted experimental data from Berger et al (2005) to determine K and λ. It was found that λ is approximately 0.71 whereas the value of K depends on the material.
In order to derive analytic electron scatter kernels we propose and use the following assumptions.
1. Primary interactions of monoenergetic photons create monoenergetic electrons. The energy of the electron will be the photon energy in the case of photo electric absorption, i.e. the ionisation energy of the atom is neglected. In the case of Compton scattering a mean energy transfer is assumed for all electrons, irrespective of the scattering angle. 2. The electrons are expected to travel on straight lines and the statistical process of slowing down is neglected. 3. Since electrons do not travel very far the material is assumed to be locally homogeneous.
I.e. material variations are on a scale larger than 300 µm. 4. Production of bremsstrahlung is ignored. 5. The stopping power is approximated as S = K · E −λ . 6. Electrons scatter isotropically in all directions.
Electrons that travel a distance greater than 5 µm have kinetic energies which are much larger than the ionisation energy of smaller atoms. The kinetic energies of Compton electrons depends on the scattering angle. However, the energy and range of Compton electrons is much lower than that of photoelectrons. Hence variations in the kinetic energy of Compton electrons will only affect the beam penumbras very close to the microbeams, which justifies the first assumption. Electrons change permanently directions due to elastic scattering and hence the assumption of isotropic scattering (assumption 6) can be justified. This may seem to contradict the assumption of straight electron paths (assumption 2). However, an average over many electrons is performed and in order to account for the meandering electron path the projected electron range is used in the kernel calculation (see below).
Under these assumptions the electron fluence Φ in a distance r from the primary interaction can be described as where N 0 is the number of primary electrons, e r the radial unit vector and σ the projected CSDA range of the electrons. Dose is related to the divergence of the energy fluence Ψ = E( r) · Φ( r) by D = − 1 ρ ∇ Ψ. In spherical coordinates this yields 5 Integration of equation (14), yields This allows to write the dose as a function of distance r from the electron production event by Equation (19) is the electron point scatter kernel and gives the average dose per electron in distance r from its creation. More complicated electron sources can be treated as a set of point sources and the total dose can be obtained from where ν( r) describes the distribution of the electron source and is defined such that ν( r)d r 3 is the volume V of the considered voxel. The material is considered to be homogeneous within the range of the electrons (i.e. approximately 300 μm), which is a reasonable assumption since structural information on smaller spatial scales are usually not available. If we further neglect the beam attenuation of the primary photon beam within the range of the electrons, ν becomes a function of one dimension perpendicular to the microbeam planes only, that may be denoted with y. It is therefore possible to rewrite equation (21) in cylindrical coordinates (s cos φ, y, s sin φ) as where we write r for y 2 + s 2 . Substitution of x = 1 − 1 σ s 2 + y 2 allows to simplify The integrand of the integral over x does not depend on the material or photon energy any more and it is therefore convenient to introduce the function I α (q) = q 0 x α /(1 − x)dx 6 . The factor N0E0 ρV can be identified as the primary photon dose D Primary and the electron scatter kernel for an infinitesimal plane of electron emitters obtains the shape The dose distribution of microbeams is retrieved by convolving this kernel over the primary fluence of photons and adding the scattered photon dose as background, This equation illustrated that the valley dose comprises of electrons scattered into the valley and scattered photons. In the peak the primary dose is usually outraging the scattered dose by several orders of magnitude. Due to the short range of electrons the calculation of peak doses, 6 Although I α (q) only exists in its integral form, it is possible to show that valley doses or penumbras widths can be performed on a micrometre scale within the voxel of a CT-image, where scattered and primary photon dose are considered to be constant.

The dose calculation and therapy planning system
Photon absorption and scattering of the photon beam in the traversed material is determined by the Compton scattering and photoelectric absorption coefficients µ C and µ Ph . They depend on the atomic number Z, mass number A, density ρ and energy κ = E/(m e c 2 ) (m e is the electron mass and c is the speed of light). In the given energy range the absorption coefficients are approximately given by Bartzsch and Oelfke (2013) µ Ph = α Ph κ 3.27 F 1 and (26) Here α Ph = 1.65 · 10 −9 m 2 kg −1 , α C = 3.549 · 10 −2 m 2 kg −1 and β C = 3.083 · 10 −2 m 2 kg −1 are constants, obtained by fitting absorption and scattering coefficients from Berger et al (2010) in the for MRT relevant energy range. Material dependencies are described by the parameters for photoelectric absorption and Compton scattering, respectively. They depend on the mass density ρ and the weight fractions ω i of the elemental compounds i. Z i and A i are the atomic number and mass number of component i. Material composition and density are computed from the Hounsfield units (HU) of a diagnostic CT image using the method described by Schneider et al (2000). For the dose calculation the CT cube is rotated into the beams-eye-view such that the beam propagation axis is parallel to the z-axis. This is possible since the radiation field can be assumed to consist of parallel beams only (Bartzsch et al 2014). The CT-voxel size is adapted to the voxel size of the photon scatter kernels, which are calculated in advance. Interaction strength and primary photon dose are calculated, and afterwards the photon scatter dose is computed by multiple three dimensional convolutions in Fourier space, sequentially for different energies and materials according to equation (3). The convolution results are added to the total scatter dose (equation (4)). All Fourier transformations use the fftw libraries (Frigo and Johnson 2005).
For all CT voxels within the treatment field microbeam peak and valley dose values are calculated by convolution of the electron scatter kernel with the primary photon fluence in real space. Electron kernels comprise of a Compton and photo electron part and are calculated in each voxel separately depending on the voxel material and on the energy of the primary photons.
The described point kernel algorithm is attached to a treatment planning system (TPS) in order to facilitate MRT therapy planning at the ESRF. The authors implemented the algorithm in C++ and attached it to the DKFZ treatment planning routine VIRTUOS, enabling CT image visualization and organ segmentation as well as the creation of treatment plans and the evaluation of dose distributions. The treatment plan defines relevant parameters of the beam set-up, such as irradiation angles, beam isocentre, field size, centre-to-centre distance and width of the microbeams. VIRTUOS passes the relevant data, e.g. the CT of the patient and treatment plan, to the dose calculation algorithm, which performs the dose calculation. The resulting peak, valley dose and PVDR distributions are imported as binary data into VIRTUOS and can be displayed and evaluated.

Monte Carlo simulations
For comparison Monte Carlo (MC) simulations are performed in the Geant4 tool kit version 9.5 using the livermore low energy libraries for photons and electrons. Analogue to the analytic dose calculation it is based on VIRTUOS plan files and CT cubes. The conversion of Hounsfield units to material parameters is identical to that used in the convolution based dose calculation program. However, due to memory restrictions only a smaller number of materials is simulated. The number of simulated materials is set to 100, i.e. 100 different, equidistant Hounsfield values between the highest and lowest value found in the CT are distinguished in the simulation.
The MC simulation is performed in the reference system of the CT-image. The radiation source is defined in a beam reference frame S spanned by the three pairwise orthogonal unit vectors e xs , e ys and e zs . This reference frame is at the same time the laboratory frame. In this reference system the beam propagates along e zs and the microbeams are parallel to e xs . The origin of the beam reference system is the centre of the microbeam field at the MSC. If we denote the width of the microbeam field by W Field and the width of the microbeams by w then the n th microbeam (n ∈ N ) is simulated with a rectangular profile from y s = −W Field /2 + (n − 1) · ctc to y s = −W Field /2 + w + (n − 1) · ctc. The number of particles per beam are computed according to the source model in Bartzsch et al (2014). The energy distribution of simulated photons correspond to the spectrum of the ID17 beam line at the ESRF (Martínez-Rovira et al 2010) and photons are linearly polarized in the y S direction. Coordinates in the beam reference frame are converted into the reference frame of the CT-image r = (x, y, z) by r = x s · e xs + y s · e ys + z s · e zs (29) Peak and valley doses are scored on the voxel grid defined by the CT image. Since the microbeam width is usually smaller than the voxel size, dose accumulation has to be separated into a peak and a valley part. To achieve this the interaction point P is transformed into the beam reference frame by projection onto the vectors e xs , e ys and e zs , x s = < e xs , P >, y s = < e ys , P >, z s = < e zs , P > .
Absorbed energy is assigned to the peak if the interaction occurs in a microbeam path while it is assigned to the valley if the interaction occurs in the central 80% of the inter-beam region, avoiding the beam penumbra regions. The calculation of radiation doses requires mass and hence volume of peak and valley scoring region, as illustrated in figure 1(a). Since all voxels are identically oriented the fraction of the voxel covered by the peak and the valley scoring region depends only on the position of the voxel parallel to y s . This fraction is calculated in a numerical integration prior to the Monte Carlo simulation. Figure 1(b) shows the changing volume fraction of the peak region depending on the position of the microbeams for two beam orientations.

Example beam geometries
Dose calculation and treatment planning is demonstrated for the following selected examples and the performance of the kernel based dose calculation is compared to Monte Carlo simulations. The width and centre to centre distance of the microbeams is always set to 50 and 400 μm, respectively.
1. Adipose tissue, white matter and cranium are three materials of very different radiological properties. For a monoenergetic microbeam field of 20 × 20 mm 2 field size the calculation of peak and valley doses at 55 keV in adipose tissue, 135 keV in white matter and 105 keV in cranium is performed. 2. A microbeam profile is calculated in white matter for polychromatic photons with the convolution based algorithm and is compared to Monte Carlo results. The beam spectrum follows the ESRF ID17 spectrum as it is used for most of the preclinical trials. 3. An anthropomorphic head phantom (Radiosurgery head phantom, CIRS, Norfolk, USA) is used to demonstrate the dose calculation in a realistic geometry and a calibrated CT image of the phantom is used to compute material properties. Dose calculation results are compared to Monte Carlo simulations and film dosimetry.

Film dosimetry
Film dosimetry was carried out at the biomedical beamline ID17 of the ESRF using the protocols described in previous work (Bartzsch et al 2015). Radiochromic HD-V2 films were cut into 3 × 3 cm 2 sized pieces. Always 8 film pieces can be obtained in each row from a single film sheet. Four of these pieces were used for calibration and were irradiated with a homogeneous dose of 90, 120, 150 and 180 Gy, the other four films were used to separately measure peak and valley doses at different positions in the human head phantom. The peak entrance dose was adjusted such that either peak or valley doses were in the range of the calibration. As an example figure B shows the volume fraction for microbeams that are perpendicular to the z-axis and microbeams that are perpendicular to direction (1, 1, 1) as illustrated in A. The calculation was done for a ctc of 400 µm and a beam width of 50 µm.
Films were kept in darkness for at least 48 h before they were scanned with a Zeiss Axio Vert.A1 microscope with a 5x/0.16 EC Plan-Neofluar objective lens. The microscope settings were kept constant between calibration and measurement films. The calibration data was fitted with a rational fit (Lewis 2011) and used to convert scanned pixel values into doses. Only doses within the calibration range were used for analysis. Figure 2 shows a comparisons of Monte Carlo and kernel based dose calculation in adipose tissue, white matter and cranium for a 20 × 20 mm 2 microbeam field of 50 µm wide and 400 µm spaced beams. Highest deviations between Monte Carlo and convolution based dose calculation are observed in the beam entrance region, in particular at high energies. For depths larger than 2 cm differences are below 3% and 5% of the maximum in the peak and valley, respectively. In the first 2 cm from the beam entrance peak dose differences are up to 3% in cranium and up to 8% in white matter and valley doses calculated in the kernel based dose calculation are between 10 and 15% higher than Monte Carlo simulated doses. The reason for these differences is that the kernel based dose calculation overestimates the photon backscattering. This does mainly affect the valley dose which consists to a large part of absorbed energy of scattered photons.

Homogeneous materials
Dose calculations are repeated for a cube of white matter and the photon spectrum of the ESRF, ID17 beamline. Beam profiles are calculated in an additional Monte Carlo simulation with a finer grid with voxels measuring 5 µm perpendicular to the microbeams. In the centre of the field, averaged over a depth range between 5 and 7 cm figure 3 compares beam profiles of Monte Carlo and kernel based dose calculations. The kernel algorithm overestimates the valley dose by approximately 5%. Otherwise the dose profiles are in good agreement within the 95% confidence interval of the Monte Carlo simulation. The beam intensity drops from 95% to 5% of the peak dose within 59.3 µm and 60.5 µm in the kernel and Monte Carlo algorithm. Figure 3(B) shows the electron kernel used for the convolution. Contributions at different energies and therefore different electron ranges of the kernel in equation (24) are added up to obtain the polychromatic electron kernel in this figure. Both, photoelectrons and Compton electrons contribute to the final kernel shape. Figure 4 shows a CT slice of the human head phantom. The phantom has inlets for radiochromic films. Film dosimetry is carried out at four positions: 50, 60, 70 and 80 mm from the back of the phantom. The treatment field lies in the sagittal plain and points horizontally in rostral direction. The beam path is drawn into the CT slices in figure 4. Results of the dose calculation in the human head phantom are shown in figure 5. Deviations between peak doses calculated by Monte Carlo and the developed convolution algorithm are smaller than 3% in soft tissue and smaller than 8% in bone. Deviations in the valley are up to 8% in soft tissue and 20 % in the skull. The differences between measurement and convolution algorithm are below 4% in the peak and below 15% in the valley. A comparison between Monte Carlo and film dosimetry gives similarly up to 4% differences in the peak and up to 10% differences in the valley.

Human head phantom
The convolution based dose calculation algorithm is less accurate for the valley dose in small bone structures, such as the skull in the human head phantom. The shape of the scatter   kernels does only depend on the material at the primary interaction point and hence changes in the photon absorption and backscattering are neglected. As a consequence scattering contributions in the skull of the human head phantom from adjacent less dense materials are overestimated and the valley dose in the skull is higher than predicted by Monte Carlo.

Computation times
Calculation times strongly dependent on the performance of the computer system and therefore the data presented in the following apply to the specified hardware configuration, only. The Monte Carlo simulations in this study are based on 10 8 particle histories. On a 3.4 GHz processor with 8 GByte RAM the simulations take around 107 processor hours. The kernel based dose calculation of the primary and scattered photon dose take around 3.2 processor minutes. Another 0.5 min are required to rotate and interpolate the CT and dose cube before and after the convolution. The calculation of peak and valley dose values takes 10.2 µ s per voxel. The dose calculation is performed on a 2 mm side length grid, a phantom depth of 200 mm and for a field size of 20 × 20 mm 2 and hence peak and valley doses are calculated for 10 4 voxels, which takes 102 ms.

Discussion
A dose calculation algorithm for microbeam radiation therapy is presented that is based on analytically calculated electron and photon scatter kernels and is capable of calculating microbeam peak and valley doses in less than 5 min on any state-of-the-art desktop PC. In the calculation photon scattering is separated from electron scattering and the calculations are carried out successively. While photon scattering contributes to the valley dose, electron scattering is responsible for the microbeam shapes and penumbras. In contrast, Monte Carlo simulations require several hours for a complete dose calculation. Parallelisation of the code, which is not yet implemented, could further reduce calculation times of the point kernel algorithm.
Calculated doses agree with Monte Carlo simulations within around 5% in the peak and around 8% in the valley. Close to bone soft tissue interfaces and in the beam entrance region valley dose differences can locally be much higher and reach up to 20%. Microbeam profiles computed by Monte Carlo calculations and by the kernel based approach resemble each other closely. Beam penumbra widths calculated in Monte Carlo simulations and the convolution algorithm match within 2%. For future clinical trials the uncertainty in the peak and valley dose delivered to a target regions should be less than 10%, a challenge also with respect to dosimetric validation.
The presented convolution algorithm accurately reproduces microbeam profiles calculated in Monte Carlo simulations as demonstrated in figure 3 for photon energies that are typically used in MRT. What are the limits of the developed model in terms of photon energy?
The assumption that radiative energy loss can be ignored (assumption 4) is valid up to very high energies. Even at 1 MeV the contribution of radiative energy loss to total energy loss is smaller than 1% (Berger et al 2005). From around 5 MeV radiative energy loss cannot be ignored anymore. The assumption of constant stopping power (assumption 5) is also valid for a wide range of photon energies. However, at energies above 500 keV range straggling will lead to differences between model prediction and observation close to the end of the electron range. The assumption of homogeneous material (assumption 3) can only be applied as long as the electron range is short compared to the size of structures in the phantom. In fact even for photon energies at around 100 keV the energy absorption in organs with strong microscopic inhomogeneities, such as lung, may not be described correctly and requires further investigation. To account for the microscopic dose absorption the determination of dose to water instead of dose to material in MRT treatment planning may be advisable, especially in view of the low kinetic energies of secondary electrons (Enger et al 2012).
The assumption of isotropically, mono-energetically emitted electrons (assumption 6 and 1) that propagate on straight paths (assumption 2) can only be made in a limited range of photon energies. The approximation of straight paths becomes worse at lower electron energies, but there the range of these electrons is very short. At higher photon energies Compton electrons are emitted increasingly anisotropic with a broad spectrum of kinetic energies. Therefore the electron convolution algorithm becomes inaccurate above photon energies of around 300 to 400 keV. By including spectrum and direction of the Compton electrons into the model, an application of the algorithm at higher energies is conceivable, though.
Synchrotron radiation is linearly polarized, which influences the differential scattering cross section of the Compton effect. Currently, the presented algorithm does not take polarization effects into account. However, in-field dose and PVDR differences between microbeams produced by polarized and unpolarized photon beams are expected to be less than 3% (Bartzsch et al 2014). Larger differences are only expected in the out-of-field scatter dose. In future work the unpolarized photon and electron scatter kernels will be replaced by polarized scatter kernels.

Conclusions
The developed algorithm forms part of the treatment planning system at the biomedical beamline ID17 of the ESRF. It is applied to calculate peak and valley doses for preclinical and pet-clinical trials. In a comparison with film dosimetry measurements in an anthropomorphic head phantom at the ESRF, differences between calculated and measured peak doses are smaller than 4%, however, up to 15% in the valley. Also Monte Carlo simulations show around 10% lower valley doses when compared to measurements. Film dosimetry of microbeams is very delicate and sensitive to various readout settings and conditions (Bartzsch et al 2015). Repeated measurements are required in the future. This will reveal whether systematic differences in the order of 10% persist between dosimetry and dose calculation.
The main weakness of the presented dose calculation is the description of scattering at material interfaces. In order to overcome these inaccuracies a local deformation of the scattering kernels would be necessary. A superposition algorithm with point kernels that account for tissue inhomogeneities and related changes in the photon scattering in an efficient way could further improve the accuracy of the dose calculation and will be part of future developments.
Currently the treatment planning system allows only the definition of rectangular microbeam fields and dose calculations for planar microbeams. Methodically there is, however, no limitation to apply the developed methods for irregularly shaped fields, cross firing beam geometries or even pencilbeam (Schültke et al 2013) irradiations. A translation of the methods into the more common treatment planning system Eclipse (Varian Medical Systems) is envisaged and will allow treatment planning for more complex beam geometries in the future.