Feasibility assessment of the interactive use of a Monte Carlo algorithm in treatment planning for intraoperative electron radiation therapy

The physical properties of the electron beam for each energy and applicator diameter are characterized in the form of a PHSP model defined at a plane normal to the geometrical axis of the applicator. The PHSP describes the properties of the beam as a combination of multiple annular electron sources with nominal energy E_i and radius R_i, radiating along a line defined by the PHSP plane at the angles (θ_i, φ_i) (Herranz 2013). A weight factor w_i established the contribution of each elemental source to the beam. The model for the treatment unit, i.e. the computation of the individual weights w_i, was generated using a robust iterative maximum likelihood expectation maximization method (Herranz et al 2013), which provided a plausible set of weights for the existing experimental measurements. The iterative procedure consisted of two main steps: forward and backward projection. The forward projection provides, by the principle of linear superposition of doses, the dose D^model produced by the current PHSP as a weighted linear combination of the individual dose footprints d_i generated by each elemental source. The measured dose D^exp was then compared with D^model, and correction factors were obtained for each voxel of the volume being considered. The correction factor for each bin is computed based on the ratio of measured to modelled dose.

To efficiently generate particles out of the PHSP during simulation, elemental sources were sorted and stratified based on the value of their individual weight factors. The n₁ microsources with the highest weight factor values w_i belong to the first cluster, and similarly for the other microsources. The number of microsources n_j in the cluster C_j was selected so that it had a predefined a priori probability Pj.

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{C}_{j}}=\left\{{{w}_{{{s}_{j}}}}\ldots {{w}_{{{s}_{j}}+{{n}_{j}}}}\right\},\forall \left(i < k\right):{{w}_{i}}\geqslant {{w}_{k}} \\ {{P}_{j}}=\sum\limits_{i={{s}_{j}}}^{{{s}_{j+1}}-1} {{w}_{i}} \\ \end{array}\end{eqnarray} \tag{ 1 }$$

The actual source k was selected in two steps. Firstly, a random uniform number was generated to select the cluster j, and then a random element k was taken from the subset C_j using the rejection function r(k).

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} r\left(k\right)=\frac{{{w}_{k}}}{\text{max}\left\{{{w}_{k}}\right\}}> \eta ,{{w}_{k}}\in {{C}_{j}} \\ \eta ::U\left(0,1\right) \\ \end{array}\end{eqnarray} \tag{ 2 }$$

Particle variables (E, θ, φ) were generated randomly based on the nominal value of the randomly selected microsource k and the discretization binning (ΔE, Δθ, Δφ). Particle position (x, y) in the PHSP plane was sampled directly from a ring with nominal radius R_k using the following expressions:

$$\begin{eqnarray}&&\begin{array}{*{20}{l}} {\xi ::U(0,1),\varsigma ::U(0,1)}\\ {{r_2} = {R_k} + \frac{{\Delta R}}{2},{r_1} = {R_k} - \frac{{\Delta R}}{2}}\\ {\left. {\begin{array}{*{20}{l}} {\theta = 2\pi \varsigma }\\ {\rho = \sqrt {{r_1} \cdot {r_1} \cdot (1 - \xi ) + {r_2} \cdot {r_2} \cdot \xi } } \end{array}} \right\}\begin{array}{*{20}{c}} {x = \rho \cdot {\rm{cos}}(\theta )}\\ {y = \rho \cdot {\rm{sin}}(\theta )} \end{array}} \end{array}\end{eqnarray} \tag{ 3 }$$

The PHSP was located above the exit plane of the applicator in a plane normal to its geometrical axis. The applicator was inserted in the volume as a methacrylate hollow cylinder, as shown in figure 2, overriding any material information coming from the CT image. With this approach, the effect of the bevel on the dose became part of the simulation, reducing the number of PHSP models to a single plane per nominal energy level and applicator diameter.

Zoom In Zoom Out Reset image size

Each voxel of the simulation volume was automatically assigned a tissue type and a density using a stoichiometric method (Schneider et al 1996, Vanderstraeten et al 2007) that generated device-specific conversion tables with the aid of a calibration phantom. In our implementation, model fitting was achieved by a constrained nonlinear multivariable optimization that searched for the set of cross-section values (K_ph, K_coh, K_KN) that minimized the maximum error between the model and the observed Hounsfield values in the calibration phantom. Constraints were imposed to avoid the convergence of cross-section fitting constants into negative values observed in the original algorithm (Kock and Schreuder 1996). Non-constrained results may provide a good fit but have no physical meaning.

A calibration phantom was acquired with an Mx8000 CT Scanner to validate the model-fitting algorithm developed for the scanner calibration. For this scanner, the algorithm provided the fitting constants: K_ph/K_KN = 1.001  ×  10⁻⁸ and K_coh/K_KN = 1.212  ×  10⁻⁵. These values provided the relationship between simulation tissues and Hounsfield values for this particular machine.

ROs need to define the expected surgical frame on the TPS by removing the tissues that will be resected or displaced during surgery from the simulation volume, as well as including new elements not found in the CT volume but present during dose delivery. Beam modifiers, such as bolus, are placed at the end of the applicator or on the patient surface, and shielding is used to protect critical structures (Beddar et al 2006). Twelve different protection materials, including lead, brass, aluminium and polymethyl methacrylate, among others, are currently available to define protection shields.

Protection systems and modifiers are 3D structures for which the user needs to specify dimension, orientation, shape and coordinates in the volume. During the planning procedure, these structures are inserted in the simulation volume, overriding the local tissue type and density given by the CT for each voxel and replacing it with the corresponding material type.

Dose to the medium was computed using a custom implementation of the dose planning method (DPM) (Sempau et al 2000), a 'mixed'class algorithm that runs in a fraction of the time of the original computer program. DPM uses an analogue simulation of photon interactions and a class II condensed history method (Kawrakow and Bielajew 1998) for electron simulation, treating large energy-transfer collisions in an analogue sense and using the continuous slowing-down approximation (CSDA) to model small energy-loss collisions. Cross-sections and mean free paths for each energy and material type were computed and tabulated offline. During simulation, these physical parameters were estimated by interpolation.

In our implementation, computational acceleration is achieved by optimizing the algorithm flow and efficiently using the hardware resources available in a modern multicore computer. The graphic processor unit, which can also be used to speed up the computation of dose (Jia et al 2011, Hissoiny and Ozell 2011), is devoted to image rendering in the TPS.

2.4.1. Parallelization.

2.4.2. Algorithm.

2.4.3. Compiler.

2.4.4. Hardware.

Particles are generated and simulated in blocks to increase data and code locality, and to improve throughput. The memory layout of interpolation spline coefficients is arranged to have them stored in consecutive memory positions, so that they are fetched from external memory in a single access and stored in a single cache line, thus reducing the number of memory accesses and cache misses.

The statistical uncertainty of the simulated dose is determined using the history-by-history method (Walters et al 2002), which provides the estimated uncertainty σ_j² for each voxel as:

$$\begin{eqnarray}&&\sigma _{j}^{2}=\frac{1}{N-1}\left(\frac{{\sum}_{i=1}^{N}D_{i}^{2}}{N}-{{\left(\frac{{\sum}_{i=1}^{N}{{D}_{i}}}{N}\right)}^{2}}\right),\end{eqnarray} \tag{ 4 }$$

where N is the number of simulated histories and D_i is the deposited dose per voxel. Periodically during the simulation, the maximum uncertainty is calculated for those voxels whose dose is higher than 50% of the maximum dose (D_max). Simulation execution ends when the maximum uncertainty is lower than a value specified by the user, which in our simulations was usually set to 2%.

At the end of the simulation, the worst-case uncertainty σ_max² is calculated for those voxels where the dose was higher than 50% of D_max. Simulation efficiency ε, is defined as

$$\begin{eqnarray}&&\varepsilon =\frac{1}{\sigma _{\text{max}}^{2}T},\end{eqnarray} \tag{ 5 }$$

where T is the CPU time required to run the simulation, is reported after every simulation.

The accuracy of the underlying algorithm to reproduce dose distributions, obtained either by using other MC codes (Sempau et al 2000) or measured experimentally (Chetty et al 2007), has already been shown by other authors. However, the proposed implementation introduces multiple optimizations and modifications that had to be verified.

The verification of dose distributions was performed in three steps. First, DOSRZnrc (Kawrakow and Rogers 2000), which is widely used by the medical physics community for the calculation of dose deposition, was used as a reference to assess the accuracy of the MC simulation. Second, experimental measurements in a water tank were compared with simulated dose profiles using a PHSP model of the accelerator obtained with our iterative method. Finally, monitor units (MU) required to deliver a given prescribed dose with different bevel angles and field sizes were computed following the standard protocol used in the clinic, and these values were compared with the dose profiles obtained using the MC method.

Dose distributions for a monochromatic parallel source, with beam energy between 4 and 20 MeV, in water and slab phantoms were simulated to compare our MC simulation with DOSRZnrc. Two alternative slab phantoms, consisting of water, aluminium and lung, were used to stress the code with a high-energy beam and a sharp inhomogeneity (Rogers and Mohan 2000).

PHSP models for the Varian 21EX system, with nominal energy in the 6–20 MeV range and an applicator diameter of 6 cm, were generated using the iterative method described in section 2.1. For each energy, a PHSP model was estimated using water measurements with a right-angled applicator. The 6 cm applicator was used as an example to show the capability of reproducing dose distributions at different energies and bevel angles. The computation of dose footprints d_i was carried out offline in a computer cluster, with a computational effort equivalent to 10 months of a single-core processor. Once the dose footprints were computed, the iterative algorithm found a set of weights fitting the estimated dose to the measured dose in water for this specific accelerator. A full iteration of the fitting algorithm was completed when the multiplicative factor had been computed for all bins, and all bins in the PHSP were updated using these factors. The entire fitting procedure typically required around 200 full iterations before convergence, which took a few minutes to fit each energy and cone diameter using a single-core processor.

MC simulations were then used to obtain the number of MUs required to deliver a given prescribed dose with different bevel angles and energies, and these values were compared with those obtained through measurement and hand computation (Cygler et al 1997).

Following the recommendations of the Spanish Medical Physics Society, absorbed dose to water protocols were used for calibration (IAEA2000). Relative measurements were taken in water and air with an MP3 radiation field analyser (PTW, Freiburg, Germany), equipped with a p-type electron diode detector (Scanditronix Medical AB, Uppsala, Sweden). Another diode detector was placed in the periphery of the radiation field during scanning to acquire reference conditions. The TANDEM electrometer (PTW) was used for all scanning measurements. The uncertainty of such a measurement protocol is less than 1% in dose and 0.5 mm in distance. P-type electron diodes were also used to measure cone output factors in water (Björk et al 2004), with an uncertainty lower than 0.5%.

MUs were calibrated for each energy to provide 1 cGy at the maximum of the percentage depth dose (PDD) in water, per MU, with the standard 10 cm × 10 cm applicator and a source-to-surface distance (SSD) of 100 cm. Output factors (OF) were computed (Almods 1999) at the maximum of the clinical axis with three different SSDs for each bevel-diameter combination and intermediate distances were estimated by interpolation, taking a 10 cm × 10 cm applicator as a reference. To estimate MU with the MC method, a gain factor was computed for each PHSP so that the simulation output was in the same units as the water measurements. In this way, MUs were estimated for a 127 cm applicator using the OFs computed previously.

Dose–depth profiles in a water phantom (figure 3) and slab phantoms (figures 4 and 5) were simulated using a monoenergetic beam at different energies. Dose profiles in water of figure 3 confirmed agreement between the reference (DOSRZnrc) and computed profiles, with a difference between curves lower than 1.5% of the maximum dose at max.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The slab phantom used in figure 4 demonstrate the accuracy of the MC model to reproduce dose profiles in the energy range of interest for procedures with protection around the lung, with a relative error below 1.5%. Additionally, the slab phantom described in figure 5 was used to stress the code with a high-energy beam and a sharp inhomogeneity (Rogers and Mohan 2000). The results show good agreement between the profiles, with a slight mismatch along the lung. In fact, this mismatch was smaller than that found using the original computer program (Sempau et al 2000), as shown in the upper-right corner of the figure 5. The differences were attributed to the use of a finer grid for the precomputed cross-sections, as well as to the algorithmic changes introduced.

These results show that the MC method implemented accurately predicted dose profiles in the presence of sharp inhomogeneities, which is a common scenario in IOERT where high atomic number (Z) materials are used to shield organs at risk. In addition, the computations required a fraction of the time needed by the original computer program, as discussed in section 3.3.

Figure 6 compares the experimental and simulated PDD curves for a 6 cm right-angled applicator for nominal energies between 6–20 MeV, while figure 7 compares the experimental and simulated lateral profiles at 16 MeV for a 6 cm bevelled applicator, with bevel angles between 15° and 45°. Figure 7 shows the ability to model the effect of the bevel and to predict the shape of the lateral dose profiles. The fact that measurements with bevelled applicators were not employed to fit the PHSP must be taken into account, and thus they are a good indicator of the validity of the PHSP + MC algorithm for this configuration. However, the ability to reproduce measurements depends greatly on the PHSP model, and, therefore, a complete validation of the iterative methods used to model the phase space is required; such a study falls beyond the scope of this work and has been addressed by (Herranz et al 2013).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The gamma function (Low et al 1998) was used to compare the experimental profiles with the simulation output for all energies and bevel angles, with the result that the gamma function was smaller than that corresponding to a dose difference of 3% and to estimating distance to an agreement of 3 mm, which is well within the needs of IOERT.

As a further test of the suitability of the MC + PHSP for this case, MUs were measured and computed for dose distributions in a water phantom with an SSD of 127 cm and bevels of 0°, 15°, 30° and 45° to verify that the resulting values could be trusted. Table 1 summarizes the required number of MUs estimated by the MC tool for each energy and bevel, and enabled comparison of this value with manually computed MUs using the measurement data.

The ratio between values calculated by the MC tool and the measurements confirmed their consistency, with a mean value of 1.004 and a standard deviation of 0.013 for the 6 cm applicator. By averaging these values for all measured applicator diameters (3, 4, 5, 6, 7, 8 and 9 cm), the mean value of the ratio was 0.999 and the standard deviation was 0.016, values that are comparable with those reported by other authors (Cygler et al 2004) for electron beams with a different MC-calculation engine.

The algorithm was profiled to identify and progressively optimize bottlenecks. Particular attention was paid to the generation of random numbers and the sampling of different probability density functions, portions of code used intensively in the MC.

The actual load distribution depends on the object being modelled, as well as the energy and applicator size. With the aim of reducing the computation time in worst-case scenarios, only configurations with wide applicators and high energy were considered for code profiling. The distribution of workload in a water phantom with an 8 cm applicator at 20 MeV, summarized in table 2, is very similar to the distribution in the reference implementation of the algorithm (Sempau et al 2000). Simulation effort was divided into four main modules: particle generation, simulation management, particle transport and energy recording. The simulation of electron physics took most of the computational effort (71.7%). The simulation of electrons was split into the implementation of the CSDA and the simulation of discrete events, such as the generation of bremsstrahlung photons. In both cases, the computational effort was divided into the actual transport of particles within the volume (geometry) and statistical sampling of the different physical processes to compute displacement and rotation angles.

The sampling of the PHSP was costly without stratification. Four layers were enough to reduce the effort of particle generation from 30% to 3.9% of the total simulation time. Moreover, random number generation, which was used at multiple points, accounting for 12% of the computing time. This fact justifies the use of an efficient re-entrant combination of linear generators.

To assess the impact of the different modifications introduced to the original implementation, the efficiency gain provided by the MC tool in every run was evaluated for a fixed number of particles after algorithmic improvements (CODE OPT), memory layout and compiler considerations (MEM OPT) and simulation parallelization into multiple threads (THREAD). These values were compared with the efficiency of original code. The results, summarized in figure 8, show that the efficiency doubled after algorithmic and memory optimizations, and an additional factor of four was achieved by execution on multicore architecture. Because the underlying physics were the same, the efficiency gains were due mostly to a reduction in the computation time.

Zoom In Zoom Out Reset image size

The motivation behind these optimizations was the intention to use the TPS and the MC to perform treatment preplanning and intraplanning. The intraoperative simulation helps to take into consideration the actual setting in the operating room and to validate hypotheses assumed during planning before the actual delivery. In clinical practice, where the patient is anaesthetized and undergoing surgery, the aim is that the MC computation should not delay the procedure, and it is necessary that the simulation be completed in less than 3 min. In dose simulations with a water phantom with 1  ×  1 × 1 mm³ voxels, it was observed (figure 9) that such a restrictive time constraint was not met when wide applicators and high energies were used, even for a relatively lax uncertainty requirement of 4% of the dose at the maximum, which required the simulation of more than 30 million particles.

Zoom In Zoom Out Reset image size

In figure 9 it is observed that, if we restrict ourselves to the combinations used most frequently in clinics (where energies above 12 MeV are rarely used) and consider a larger voxel size (8 mm³ instead of 1 mm³), accurate dose distributions, with uncertainties lower than 2%, could be completed in less than 3 min in most scenarios.

Two representative scenarios of treatment planning in breast cancer are presented to illustrate the utility and limitations of the workflow that was developed and implemented. In case 1, the preoperative image would typically be used to analyse possible treatment approaches. In case 2, the intraoperative CT image would be used to confirm whether the actual set-up in the operating room, including the applicator and position of the protection, meets the prescribed requirements before dose delivery.

CT studies were acquired at Hospital Gregorio Marañón in Madrid on a Toshiba Aquilion LB scanner (Toshiba Medical Systems Corporation, Tochigi, Japan) with the following parameters: voltage 120 KVp, current 40 mA, slice thickness 3 mm and pixel spacing 1.8  ×  1.8 mm² for case 1, and slice thickness 5 mm and pixel spacing 1  ×  1 mm² for case 2. The intraoperative CT was obtained as part of a research protocol in which a specifically designed couch was used to transport the patient from the OR to the CT scanner under anaesthesia, with the applicator treatment position fixed using an articulated arm. This protocol was designed to evaluate the difficulties and contributions of intraoperative CT imaging in IOERT scenarios. The optimal clinical protocol definition is still under investigation.

Case 1 emulates a preplanning workflow. In this scenario, the RO would use the preoperative image to segment the planning target volume (PTV), organs at risk and fluids on the tumour bed; select the treatment configuration (applicator position, diameter and bevel angle, as well as the nominal energy); and place the beam modifiers depending on several factors (size of the treatment target, anatomical structures, candidate areas to be protected). In this simulated scenario, shown in figure 10, the RO would prescribe 21 Gy to irradiate the PTV around the tumour bed delineated as a truncated cone shape, with a base of 40 mm and a height of 15 mm. The user, in this particular case an RO with IOERT experience following the standard clinical protocol that he follows in the OR, positioned a 5 cm diameter applicator and inserted an 8 cm diameter/5 mm thick cylindrical shield (3 mm aluminium/2 mm brass) at the location where he considered that the ribs and lung would be best protected. With these settings, the RO and the radiation physicist expected to protect the lungs and ribs completely, while the delineated target volume would receive at least 90% of the prescribed dose. Multiple simulations, with slight variations in the applicator, protection positions and angles, were run to identify the optimal applicator location in terms of PTV coverage. This case used a two-layered disc configuration, which is frequently utilized in practice; the top disc with the lower Z slows down beam particles and absorbs backscattered particles (Catalano et al 2007), and the bottom disc with the higher Z stops the residual electrons.

Zoom In Zoom Out Reset image size

Dose was estimated by considering a 6 MeV electron beam generated by a Varian 21EX linear accelerator. To obtain a statistical dose uncertainty σ_max² below 2%, approximately 11.2 million particles, with an execution time of less than 35 s, were simulated in each MC run. This execution time would allow for a certain degree of interactivity in the planning process.

After several iterations, the RO determined that in the best-case scenario (figure 11, top left) the PTV was not completely enclosed by the 90%-isodose curve obtained with the 5 cm applicator, and it was estimated that, at best, 88% of the tumour volume was inside the 90%-isodose curve. The same setting was simulated with a 6 cm applicator (figure 11, top right), and in this case the delineated PTV was completely enclosed by the 90%-isodose curve, allowing for a margin of safety in the OR for position of the applicator and/or protection misalignments with respect to the optimum location. These simulation results are in accordance with the new standard treatment proposed by the European Institute of Oncology where the recommended applicator size is now 6 cm, with an occasional use of 5 cm for very small lesions (Dall'Oglio 2013).

Zoom In Zoom Out Reset image size

In case 2, the applicator was placed in the TPS to reproduce the real applicator position, as displayed on the CT image (figure 12). The CT scan was obtained with the protection in place, although other clinical protocols for breast cancer treatment avoid the use of any protection (Fastner et al 2013). The reconstructed CT image showed severe artefacts due to the high density of the shielding materials. To simulate dose deposition, the applicator was positioned so that it overlapped the end of the applicator visible in the image, protection was placed on top of the protection with artefact visible on the image, and the tissue between the applicator and the protection was segmented as soft tissue. This workaround overcomes the limitations caused by the artefacts by replacing pixel values with the expected materials. However, the correspondence between the simulated dose and the actual delivered dose is uncertain. The simulation required 6.2 million histories and took 18 s to meet a 2% level of uncertainty. This computation time is small enough to not interfere with the surgical procedure.

Zoom In Zoom Out Reset image size