© EDP Sciences 2018

1 Introduction

The reliability of Monte Carlo (MC) simulations of clinical photon beams in radiation oncology depends on the accuracy of the linac head model including parameters of the primary electron beam. For the Varian Clinac accelerator, the manufacturer recommends to approximate the primary electron beam as monoenergetic with a two-dimensional Gaussian spatial distribution. Mean energy and FWHM of the initial electron beam are commonly estimated by comparing the simulation results with measured dosimetric data. Percentage depth dose curves (PDDs) and dose profiles are commonly used for the comparison. The estimates of the primary beam parameters can be found in numerous studies (Sheikh-Bagheri and Rogers, 2001a, 2001b; Ding, 2002; Keall et al., 2003; Pena et al., 2007; Hedin et al., 2010; Chibani, 2011; Tsiamas, 2011) and their summary is presented in Table 1 here. Conclusions of the available studies are strongly dependent on the number of measured and simulated fields and on the scoring function used. No attempt has been made in the published studies to determine the confidence region of the optimized parameter values.

The aim of our study was threefold:

• to design a quantitative, robust, and simple method for optimization of parameters of the primary electron beam, i.e. of the electron energy, E, and of the beam geometric full width at half maximum, FWHM, using experimentally measured dose data;

• to design a quantitative, robust, and simple method for uncertainty estimation of the optimized parameter values;

• to apply both the proposed methods to a set of measured clinical accelerator data.

2 Materials and methods

2.1 Linear accelerator head

We modelled the Varian Clinac 2100C/D accelerator head and its beams for the photon beam energy of 6 MV. We used the BEAMnrc, v.BEAM2013 (Rogers et al., 2005) system which is based on EGSnrc, v.v4-r2-4-0 (Kawrakow and Rogers, 2003) system for Monte Carlo modelling. The model included the target, primary collimator, exiting window, flattening filter, monitoring chamber, mirror, secondary collimator jaws, multi-leaf collimator (MLC), and cross-hair plastic. Shading was not considered in the model. Material and geometry data for the linac head structures were kindly provided by the manufacturer.

2.2 BEAMnrc and EGSnrc calculation parameters

The phase space plane was placed at the distance of 100 cm from the source. We simulated PDDs and dose profiles for the fields 5 × 5, 10 × 10 and 40 × 40 cm² in water. Dose profiles were modelled in depths 1.5, 5, 10, 20 and 30 cm. For water phantom simulation, we used the system DOSXYZnrc (Walters et al., 2007). The water phantom was simulated with dimensions 40 × 40 × 40 cm³ for the fields smaller than 40 × 40 cm² and 60 × 60 × 40 cm³ (l × w × d) for the field 40 × 40 cm². For the PDDs, we used voxels 1 × 1 × 0.2 cm³ from the surface to 30 mm depth followed by 1 × 1 × 0.5 cm³ voxels for the rest of the central axis. For dose profiles, we chose voxel size 0.5 cm in the region out of field, 0.2 cm in the penumbra region and 0.5 cm in the homogenous part of the profile, always with 0.5 cm depth.

The number of initial electrons was about 10⁷. The final number of particles in the phase space depended on the chosen energy and on the FWHM of the entrance electron beam. In the DOSXYZnrc for dose profile, we simulated 10⁸ particles for the fields 5 × 5 and 10 × 10 cm² and 10⁹ particles for the field 40 × 40 cm². For PDDs, we modelled 10⁷ particles for the fields 5 × 5 and 10 × 10 cm² and 10⁸ particles for the field 40 × 40 cm².

Cutoff energy for electrons (ECUT) was chosen as 0.7 MeV and cutoff for photons (PCUT) as 0.01 MeV, in agreement with the EGSnrc code authors (Sheikh-Bagheri and Rogers, 2001a, 2001b; Kawrakow et al., 2004; Rogers et al., 2005; Kawrakow, 2005; Kawrakow and Walters, 2006). For the bremsstrahlung splitting, we used the choice directional bremsstrahlung splitting (DBS). Number of particles originated from 1 photon was set to 1000, according to the recommendation in the BEAMnrc manual (Rogers et al., 2005). The radius of the DBS splitting field was chosen 6 cm for the clinical field 5 × 5 cm², 10 cm for the field 10 × 10 cm² and 35 cm for the field 40 × 40 cm² in the distance of 100 cm from the target. The electron splitting was turned on with redistributing charged particles in a radially-symmetric manner around the Z axis. The splitting plane was placed within the flattening filter at the distance of 12.765 cm from the source. The Russian Roulette plane was chosen several millimeters above the electron splitting plane, 12.3 cm from the source. The rejection plane was placed at a distance of 90 cm from source.

First, we had made several simulations with different parameters of the entrance electron beam (6.3 MeV and 1, 1.5, 2 mm FWHM; 6 MeV and 1, 1.5 mm FWHM; 5.7 MeV and 3.5, 2.3, 2 mm FWHM; 5 MeV and 2.5 mm FWHM) in order to get acquainted with the model behavior. Based on these initial simulations and with respect to the previously published studies cited in Table 1, four energy values were selected, namely 5 MeV, 5.5 MeV, 6 MeV, 7 MeV, along with four FWHM values, 2 mm, 3 mm, 4 mm and 5 mm. The resulting matrix of sixteen points in the (E; FWHM) plane spans over the expected region of optimum parameter values. All subsequent simulations were performed at these sixteen simulation points.

2.3 Measured data

Percentage depth-dose curves and dose profiles for 6 MV beam of the linear accelerator Varian Clinac 2100C/D were measured in the Faculty Hospital Motol for three field sizes, namely 5 × 5, 10 × 10 and 40 × 40 cm². The profiles were measured in 5 depths (1.5, 5, 10, 20 and 30 cm). Water phantom PTW MP3 with dimensions 63.6 × 63.4 × 52.3 cm³ and ionization chamber PTW semiflex 31010 with the volume of 0.125 cm were used for the measurements. PTW Mephysto software and its standard measuring protocol (medium resolution, stepwise measurement, dwell time 0.2 s) were used for the phantom control and for data acquisition.

2.4 Cost function

Gamma index (2%/2 mm criterion) evaluated according to literature (Ju et al., 2008) was used to quantify the agreement between measured and simulated data. PDDs were normalized to the dose at maximum for each field size. The simulated profile data were symmetrized and normalized to the dose at central axis in each depth. The measured data were interpolated to the points of MC simulations (voxel centres). The resulting score function, F, has the form: $$F=\frac{1}{15}{\displaystyle \sum _{i=1}^3}{\displaystyle \sum _{j=1}^5\sum _{\mathit{k}}}\frac{g_{i,j,k}}{N_{i,j}},$$ where g_i,j,k is the gama index for kth point in jth depth for ith field size, and N_i,j is the number of evaluated dose points for a dose profile in jth depth for ith field size.

Thus, we obtained one score value, F, for each simulation point (E; FWHM), i.e. sixteen score values in total. The score values are listed in the line “big simulation” of Table 2.

2.5 Parameter optimization

Our model has two input parameters, the energy of electrons, E, and the electron beam geometrical width, FWHM. The cost function, F, is a non-linear function of these two parameters, E and FWHM, with continuous first and second order partial derivatives. There are standard software tools available which can be used for function minimization and thus for optimum parameter values estimation with multiparametric differentiable non-linear models, e.g. (James and Roos, 1975; Fournier et al., 2012; Matlab, 2017).

All these software packages perform high numbers of repeated evaluations (so called “calls”) of the cost function for varying values of the input parameters. The number of calls might be several hundred for a complex two-parameter model and increases markedly with any additional free parameter. In our case, each call of the cost function, F, requires a complete simulation of the beam and of its dose profiles in a water phantom. The parameter optimization process is therefore quite time demanding.

While the above mentioned obstacle can be overcome by increasing the CPU time and/or extending the computational power used, there is a second principle reason why the commonly used parameter optimization software packages cannot be used in our case. The optimization algorithms assume in general that two subsequent enumerations of the cost function which are performed for exactly the same combination of parameter values would result in exactly the same value of the cost function. In other words, they expect the cost function to be constant in time. If this assumption is not fulfilled, the derivative-based methods for function minimization would fail and also other methods, e.g. variants of the simplex method, might not converge. Some of the minimization software packages even perform a reproducibility test at the very beginning of their parameter optimization. The cost function, F, is enumerated several times with equal values of the input parameters and the program execution is terminated if all the cost function values are not equal. However, for any model with a built-in stochastic element, all Monte Carlo models included, two subsequent model calls with exactly the same values of the input parameters would result in slightly differing values of the cost function. Stochastic models are not constant in time.

Instead of exploring which nonlinear function minimization software package can cope with this situation, we used a simple parametric approach. A two-dimensional second-order surface over the (E; FWHM) plane was fitted to the set of our sixteen score values. The values of E and FWHM at the minimum point of the quadratic surface can be considered as estimates of the optimum values of the entrance electron beam parameters.

2.6 Common 95% confidence region of parameter values

The estimated optimum parameter values calculated in the way described above have unknown uncertainties. However, the estimated uncertainty should be stated with each measured value, regardless of the complexity of the model of measurement. From the measurement uncertainty point of view, our phantom data analysis can be treated as a measurement of the E and FWHM values, even if the model of measurement is a complex one.

Most of the software packages for non-linear function minimization provide dedicated tools for parameter uncertainty estimation. These tools are usually based on the cost function second-order partial derivatives and would thus fail in case of stochastic models. In order to estimate the two-dimensional 95% confidence region in the (E; FWHM) plane, we used a simple statistical approach described in following 3 steps.

2.6.1 1600 “small” simulations

We defined a matrix of 16 points (E_i; FWHM_j; i,j = 1 to 4) in the (E; FWHM) plane around the expected optimum point (E_opt; FWHM_opt). We repeated the enumeration of the cost function 100 times in each of the 16 (E_i; FWHM_j) matrix points, the total number of “small” simulations thus being 1600. The number of simulated particles in each run was chosen high enough so that the resulting statistical uncertainty of the MC calculated dose values was < 2% (in our case the number of simulated particles is 100 times smaller than in the first set of simulations). The seeds in the random number generator were changed for every “small” MC simulation. Thus, we obtained a distribution of 100 score values (average value of the gamma index) for each of the 16 (E_i; FWHM_j) matrix points. We call these score values pre-calculated, because they are calculated first and used for further statistical analysis subsequently.

2.6.2 Distribution of 1000 score surfaces and of their minimum points

We generated a distribution of 1000 score surfaces and of their minimum points in the following way:

• we selected in random one out of the pre-calculated 100 score values at each of the 16 (E_i; FWHM_j) matrix points;

• the resulting mesh of the 16 randomly selected score values was fitted by a second order 2D polynomial function;

• the minimum point of the 2D polynomial function was found.

This procedure was repeated 1000 times. In this way, we obtained a random distribution of a 1000 two-dimensional score surfaces generated from the sixteen one-dimensional score distributions which had been pre-calculated at the 16 (E_i; FWHM_j) matrix points. The 2D distribution of the minimum points of these 1000 surfaces was used for the confidence region estimation in the next step.

2.6.3 Confidence region estimation

The resulting 2D distribution of the 1000 minimum points (E_min,k; FWHM_min,k; k = 1 to 1000) in the (E; FWHM) plane was analyzed using standard statistical tools for confidence region estimation of 2D distributions.

3 Results

3.1 Dose profiles

The resulting average values of gamma index for the dose profiles are shown in Table 2. As mentioned above, these values were fitted with a second order polynomial surface (Fig. 1) in the Matlab Curve Fitting Tool (Matlab, 2014). For the original big simulation, the minimum of the polynomial fit was found with the function fminunc at the coordinates E = 5.32 MeV and FWHM = 3.55 mm. Example of a dose profile for parameters near optimum (E = 5.5 MeV; FWHM = 3 mm) is shown in Figure 2. The most influencing data for the parameter optimization are the dose profiles of the field 40 × 40 cm², in agreement with literature (Aljarrah et al., 2006). The relative standard uncertainty of the simulated dose values was about 0.5% of the maximum dose.

Fig. 1 Two-dimensional quadratic fit of average gamma index for the big simulation set (3–D view and 2-D contours).

Fig. 2 Dose profile for 6 MV beam of Varian Clinac 2100C/D, E = 5.5 MeV, 3 mm FWHM, field 40 × 40 cm2 in depth 5 cm and 10 cm – big simulation (MC) and measured data.

3.2 PDDs

Our calculations confirmed that the PDDs are considerably less sensitive than the profiles to changes in the parameters of the primary electron beam, in accordance with previous reports (Tzedakis et al., 2004; Aljarrah et al., 2006). Therefore, we neither simulated PDDs for all sixteen (E; FWHM) combinations nor performed any small PDD simulations. We only simulated PDDs at FWHM = 3 mm for four values of E, namely 5.0 MeV, 5.5 MeV, 6.0 MeV and 7.0 MeV in order to avoid any unexpected PDD behavior. Example of these PDDs is shown in Figure 3. The relative standard uncertainty of the simulated doses beyond the maximum of the PDD curves was less than 0.5%.

Fig. 3 PDDs for 6 MV beam of Varian Clinac 2100C/D, field 10 × 10 cm2 – big simulation (MC) and measurement.

3.3 Repeated smaller simulations

The relative statistical uncertainty of the dose values in the profiles was less then 2% of the maximum dose. An example of the dose profile for a small simulation is shown in Figure 4, together with the corresponding error bars. The distribution of positions of minima of the one thousand best-fitting second-order surfaces is shown in Figure 5. The ellipse was calculated using the R-cran package “ellipse” (Murdoch and Chow, 2015) and represents a statistical estimate of the 95% confidence region. The centre of the ellipse is located near the point E = 5.3 MeV and FWHM = 3.6 mm.

Fig. 4 Comparison between big and small simulation of dose profile for a field 5 × 5 cm2.

Fig. 5 Minima distribution for small simulations and their 95% confidence region (ellipse).

4 Discussion

Using measured water phantom data and MC simulated data, we optimized the parameter values (energy; FWHM) of the primary electron beam for a clinical accelerator Varian Clinac 2100C/D and for its 6 MV photon beam. We found that the dose profiles for large fields are highly sensitive to changes in the E and FWHM values. By contrast, the percentage depth dose curves are quite insensitive. This finding is in accordance with previous studies (Aljarrah et al., 2006). It means that dose profiles for large photon fields, ideally 40 × 40 cm², are essential for efficient optimization of the primary electron beam parameters.

To find the optimum E and FWHM values, we enumerated the score function, F, in a semi-regular (E; FWHM) grid. Then we fitted the obtained score values with a two-dimensional second-order polynomial fit. In order to obtain a robust fit, it was necessary to extend the calculation grid sufficiently, so that the score values at the grid borders were multiples of the minimum score value. With our measured dose data, the energy interval 5 MeV to 7 MeV and the FWHM interval 2 mm to 5 mm were satisfactory.

We applied a simple probabilistic approach in order to estimate a 95% confidence region in the (E; FWHM) plane. Determination of the confidence region was based on a high number of precalculated simulations, each one with a reduced number of simulated beam particles. We found that using 1/100 particles of the original big simulation resulted in a statistical dose uncertainty of 2% for the small simulations which was sufficient for the optimization purposes.

The resulting 95% confidence region is represented by a tilted ellipse (see Fig. 5). The ellipse centre represents optimal parameter values which are E = 5.3 MeV and FWHM = 3.6 mm. The tilt of the ellipse reveals a marked correlation between the optimum values of E and FWHM.

Our results suggest a lower electron energy and a higher FMHM value compared to the results of some previously published studies. We have to keep in mind, however, that most of the published studies utilized measured beam data for small and medium-sized fields only, which might be a source of significant bias. It is also impossible to quantify the agreement/disagreement between our results and the previously published values since no confidence intervals are given for the previously published data by their authors.

5 Conclusions

The proposed approach which is based on a set of precalculated “small” simulations and on subsequent analysis of their random combinations can be used quite generally for estimation of the primary electron beam parameters and of their joint confidence region. The method described above is mathematically approved, objective, easy to implement, and sufficiently robust. The approach respects stochastic nature of the Monte Carlo models and works well with them. To our knowledge, no comparable method of beam data analysis has been published so far.

Once the primary electron beam parameters and their confidence regions are derived from the measured photon beam data for a particular treatment machine, the results can be used in subsequent MC simulations of more complex clinically relevant irradiation techniques.

Acknowledgements

The authors wish to acknowledge Jana Meluzínová and Lukáš Kotík from the National Radiation Protection Institute in Prague for their statistical advice. The authors are also thankful to the team of medical physicists from the Faculty Hospital Motol in Prague for their practical arrangements and their assistance with the water phantom dose measurements.

References

All Tables

All Figures

	Fig. 1 Two-dimensional quadratic fit of average gamma index for the big simulation set (3–D view and 2-D contours).
In the text

	Fig. 2 Dose profile for 6 MV beam of Varian Clinac 2100C/D, E = 5.5 MeV, 3 mm FWHM, field 40 × 40 cm2 in depth 5 cm and 10 cm – big simulation (MC) and measured data.
In the text

	Fig. 3 PDDs for 6 MV beam of Varian Clinac 2100C/D, field 10 × 10 cm2 – big simulation (MC) and measurement.
In the text

	Fig. 4 Comparison between big and small simulation of dose profile for a field 5 × 5 cm2.
In the text

	Fig. 5 Minima distribution for small simulations and their 95% confidence region (ellipse).
In the text