Efficiency enhancements of a Monte Carlo beamlet based treatment planning process: implementation and parameter study

In recent work (Henzen et al 2014a, Mueller et al 2017, 2022, Guyer et al 2022), a fully MC based TPP as illustrated in figure 1(a) was developed to create treatment plans for IMRT, MERT, MBRT, VMAT, DTRT and dynamic mixed beam radiotherapy (DYMBER) (Mueller et al 2018). These treatment plans are all deliverable in the developer mode of a TrueBeam system (Varian Medical Systems, Palo Alto, CA) equipped with a Millennium 120 MLC (Varian Medical Systems, Palo Alto, CA). Source-axis-distance (SAD) of the TrueBeam system is 100 cm. DYMBER is a combination of MERT and DTRT (Mueller et al 2018). Plans created by this TPP for MERT, MBRT and DYMBER collimate the electron beams with the TrueBeam installed photon MLC. Thus, no cut-out nor applicator is used, but the patient is moved to a reduced source-to-surface distance (SSD) to reduce in-air scatter of electrons.

Zoom In Zoom Out Reset image size

This TPP starts by importing the CT data into a research version of the Eclipse treatment planning system (TPS) (Varian Medical Systems, Palo Alto, CA) and contouring the structures including the PTV, OARs and normal tissue. Afterwards, the fields are defined in Eclipse TPS. This step is treatment technique dependent. All photon and electron fields to be delivered with intensity modulated step-and-shoot apertures from static beam directions (called static fields) are defined manually as single fields in Eclipse TPS. Similarly, VMAT arcs are defined manually in Eclipse TPS. To setup photon dynamic trajectories in Eclipse TPS, the external path finding procedure described in Fix et al (2018) is used. In this work, conventional arcs and dynamic trajectories are called dynamic fields.

The next step is to calculate the beamlets for each beamlet collection (static field or control point of a dynamic field) within the framework of the Eclipse interfaced Swiss Monte Carlo Plan (SMCP) (Fix et al 2007, Manser et al 2018). For the purpose of this work, the part of the SMCP managing beamlet calculations was re-implemented in C++ to fulfill the following three functional requirements to steer the computational performance of beamlet calculation and plan optimization:

• The user can specify the desired statistical uncertainty of the beamlets.
• The user can specify a sparse dose threshold setting dose values below this threshold to zero.
• The user can specify distances from voxels to the PTV surface from which voxels are merged to larger voxels.

These parameters are described in more detail in the sections 2.1.1–2.1.3. The beamlet calculation is designed as a single multi-threaded process calculating all beamlets of all fields assigned to a treatment plan. This enables to perform common tasks for all beamlet collections and common tasks for all beamlets of a specific beamlet collection only once. Thus, computational overhead is reduced. The workflow is illustrated in figure 1(b). First, HU values of the CT data are transformed to physical density and material composition according to a dose scoring grid and CT calibration ramp. Next, as many photon and electron dose calculation engines are prepared as there are CPU cores available. Each engine is an instance of the Voxel Monte Carlo (VMC++) (Kawrakow and Fippel 2000) algorithm for photons or the electron Macro Monte Carlo (eMC) (Neuenschwander and Born 1992, Neuenschwander et al 1995, Fix et al 2013) algorithm for electrons. Afterwards, the contours of the structures are loaded and it is identified which voxels are merged together. It follows a first loop iterating over all beamlet collections of the fields assigned to the treatment plan starting with a preparation of a beamlet collection. This includes the definition of the rectangular beamlet grid in the middle plane of the MLC. The beamlet grid is conformal to the PTV in beams eye view (static field) or the smallest possible opening encompassing all control point specific conformal openings (dynamic field). An additional beam modality and SSD dependent margin is added to the beamlet grid size to account for the beam penumbra. The beamlet size is 0.5 × 0.5 cm² for the inner 0.5 cm wide leaves and 0.5 × 1.0 cm² for the outer 1.0 cm wide leaves, specified in the plane of the isocenter. Furthermore, the preparation of the beamlet collection also calculates the transformation of the beamlet collection from the beam coordinate system to the patient coordinate system. This is needed as the beamlet sources are patient-independent pre-calculated phase-spaces located at the treatment head exit plane and provided in the beam coordinate system, while the dose calculation engines simulate the particle transport in the patient coordinate system. Then, a second loop follows, iterating over all beamlets of the prepared beamlet collection in parallel using OpenMP with as many threads as CPU cores available. For each beamlet, first the number of particles to be simulated is estimated to reach the user-defined desired statistical uncertainty. Afterwards the phase-space is loaded and a third loop (not explicitly depicted as a loop in figure 1(b)) iterating over the particles is executed using a prepared dose calculation engine. The engine transports the particle and scores the dose values in the dose scoring grid (original voxel grid). Afterwards, the sparse dose threshold is applied, and then the dose is transformed to the prepared voxel merged format.

After beamlet calculation is completed, the treatment plan is optimized using a hybrid column generation (Romeijn et al 2005) and simulated annealing (Shepard et al 2002) direct aperture optimization (H-DAO) (Mueller et al 2022) resulting in a deliverable plan. The same H-DAO algorithm is applied for all treatment techniques considered in this work. The H-DAO considers mechanical constraints such as movement ranges of the MLC leaves and also transmission of photon beams through the MLC (Mueller et al 2022), while other impact of the MLC such as particle scattering is not considered. Thus, a final MC dose calculation is performed within the SMCP framework re-calculating the dose distribution of each static aperture and each control point of dynamic fields depending on the field type (static beam direction or dynamic field). For the final dose calculation of photon beams, the source is a pre-simulated phase-space located above the secondary collimator jaws. The particles of the photon beams are simulated through the secondary collimator jaws and MLC using VMC++. For the electron beams, the multiple source model ebm70 (Henzen et al 2014b) is used. The same algorithms VMC++ (photon beams) and eMC (electron beams) are used for the final dose calculation in the patient as for the beamlet calculation. Finally, the monitor units (MUs) of the static apertures and control points are re-optimized using a limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) algorithm with the same objectives as used for the H-DAO optimization. Purpose of this re-optimization is to mitigate the degradation from optimized to final dose distribution.

2.1.1. Statistical uncertainty

The user can specify the desired statistical uncertainty u_p of the beamlets for both photon and electron beamlets separately. p stands for the particle type x (photon beam) or e (electron beam) of the beam. For each specific beamlet, the number of particles n to be simulated to reach the desired statistical uncertainty u is estimated as a function of the actual beamlet size b, actual source-to-surface-distance SSD and the actual voxel volume V (only outlaid for cubic voxels). For this, pre-simulations of beamlets are performed to determine the number of particles ${n}_{{ref},b}$ to reach ${u}_{{ref}}\,$ = 1% statistical uncertainty in defined reference conditions (see table 1). b stands here either for small (0.5 × 0.5 cm²) or for large (0.5 × 1.0 cm²) beamlets. The reference conditions are at $SS{D}_{ref,p}$ with the beam perpendicular to a 30 × 30 × 30 cm³ water phantom discretized into voxels with a voxel volume of ${V}_{{\rm{ref}}}$= 0.25 × 0.25 × 0.25 cm³. The corresponding reference SSD values are $SS{D}_{ref,x}$ = 95 cm and $SS{D}_{ref,e}$ = 80 cm for photon and electron beams, respectively. The formula used to correct n for the specific beamlet is the following:

$$\begin{eqnarray}&&n\left(b,u,SSD,p,V\right)={n}_{{\rm{ref}},b}\cdot {\left(\displaystyle \frac{{u}_{ref}}{u}\right)}^{2}\cdot {\left(\displaystyle \frac{SSD}{SS{D}_{ref,p}}\right)}^{2}\cdot {\left(\displaystyle \frac{{V}_{{\rm{ref}}}}{V}\right)}^{2/3}.\end{eqnarray} \tag{ 1 }$$

Table 1. Number of particles needed to reach 1% statistical uncertainty with small 0.5 × 0.5 cm² (middle column) and large 0.5 × 1.0 cm² (right column) beamlets in the reference conditions. These conditions are an SSD of 95 cm for photon beams and 80 cm for electron beams with the beam perpendicular to a 30 × 30 × 30 cm³ water phantom with 0.25 × 0.25 × 0.25 cm³ voxel size.

Beam	Number of particles ${n}_{\mathrm{ref},\mathrm{small}}$	Number of particles ${n}_{\mathrm{ref},\mathrm{large}}$
6X	102 000	183 000
6E	4592 000	4455 000
9E	2840 000	2457 000
12E	2439 000	1894 000
15E	2535 000	1895 000
18E	3073 000	2288 000
22E	3349 000	2278 000

The quadratic correction for the statistical uncertainty is motivated by the history-by-history approach, in which we have a statistic over the number of independent histories (equal to n as all sampled particles here are independent). The quadratic correction for the SSD is motivated by the inverse square law and the correction for the voxel size is motivated by the change of the number of histories impinging on a cubic voxel (linear increase with the area of the voxel surface orthogonal to beam direction).

2.1.2. Sparse format

The total beamlet data is very memory-inefficient, because many voxels receive zero dose or very low dose. Thus, the beamlet data can be enormously reduced by neglecting dose values below a certain sparse dose threshold. In the implemented beamlet calculation framework, the user can specify a fractional sparse dose threshold t (e.g. 1%) to determine an absolute sparse dose threshold d_t in cGy/MU according to the following formula

$$\begin{eqnarray}&&{d}_{t}=t\cdot {d}_{{\rm{\max }}},\end{eqnarray} \tag{ 2 }$$

where d_max is the maximal dose value of beamlet j to any voxel in the whole body. Dose values in voxel i of beamlet j below this threshold (d_ij < d_t ) are then set to zero. The optimization does not consider any d_ij = 0.0 leading to reduced beamlet data and to reduced optimization time.

2.1.3. Voxel merging

A larger voxel size (in units of cm × cm × cm) leads to less beamlet data but also to less accurate beamlets. A promising compromise is to increase the voxel size further away from the PTV, because a small voxel size is probably more important close to the PTV and in the PTV. This concept of distance dependent voxel merging is implemented in the beamlet calculation framework. The user can specify two distances d_m and d_l . First, all original voxels are grouped in large voxels of 4 × 4 × 4 = 64 original voxels. Those large voxels with a distance, determined from their center to the PTV surface, larger than d_l are kept as large voxels. All others are re-grouped in medium voxels of 2 × 2 × 2 = 8 original voxels. Those medium voxels with a distance to the PTV surface larger than d_m are kept as medium voxels and all others as original voxels. However, there is one exception: The original voxel size is always used for structures with voxel volume smaller than 500 original voxels such that small structures like the chiasma or the lenses are represented by a reasonable large number of voxels. Voxels overlapping the PTV contours are always of original voxel size. The concept of voxel merging is illustrated in figure 2.

Zoom In Zoom Out Reset image size

The plan optimization accounts for the different voxel sizes and that a voxel of any size may only partly overlap with a structure. During optimization each voxel is weighted by its volumetric overlap with a specific structure for the calculation of the objective function (Mueller et al 2022). Furthermore, any dostrimetric evaluation during optimization (e.g. mean dose to a structure) also considers these volumetric overlaps of voxels and structures.

The compromise between final treatment plan quality and computational efficiency for the parameters listed in table 2 is studied for several parameter values. Each parameter is tested separately by running the whole TPP with the different parameter values listed in table 2. The parameter values of all other parameters stay with the most conservative value listed in table 2 (e.g. when d_m /d_l is tested for several values, u_x , u_e and t have always the values 2.5%, 5% and 0.1%, respectively). All beamlet calculations and optimizations are performed on a single Intel Broadwell CPU with 2 × 10 cores.

This parameter study is performed for five different treatment techniques, each planned for another clinical or academic case of another treatment site as listed in table 3. Original voxel size is always 0.25 × 0.25 × 0.25 cm³. The academic case consists of a 30 × 30 × 30 cm³ water phantom including contours of a wedge-shaped PTV and a single OAR surrounding the PTV. The PTV is located 0.5 cm from the body surface and the maximal depth from the closest body surface to the PTV is 5.5 cm.

The objective function, which is minimized by the H-DAO and also the MU re-optimization after final dose calculation, consists always of the same objectives per case. There are always a max-dose and min-dose objective for the PTV, a normal tissue objective penalizing dose values depending on their distance from the PTV and one generalized equivalent uniform dose (gEUD) (Niemierko 1999) objective for each OAR. The whole formalism of the objective function value is described in Mueller et al (2022). The value of the tissue-specific factor a_s used in the gEUD formula is OAR specifically selected based on the work of Luxton et al (2008). The weights and dose or gEUD values of all the objectives are identified by iterative adjustments of the human planner.

The dependency of the final treatment plan quality on the single parameter values is evaluated with the objective function value, dose homogeneity index HI in the PTV, the gEUD value averaged over all OARs ${\overline{gEUD}}_{OARs}$ using the tissue-specific factors a_s already used for the objectives and the mean dose ${\bar{D}}_{NT}$ in normal tissue. All these quantities are evaluated after MU re-optimization of the final dose distribution, which has always a statistical uncertainty of 1% or lower. Furthermore, there is never any sparse format nor voxel merging applied for the final dose distribution and the original voxel resolution is the same as for beamlet calculation.

HI is expressed by

$$\begin{eqnarray}&&\mathrm{HI}=\frac{{D}_{2 \% }-{D}_{98 \% }}{{D}_{p}},\end{eqnarray} \tag{ 3 }$$

where D_p is the prescribed dose and D_2% and D_98% the dose received by at least 2% and 98%, respectively, of the PTV. ${\overline{gEUD}}_{OARs}$ is expressed by

$$\begin{eqnarray}&&\begin{array}{l}{\overline{gEUD}}_{OARs}=\displaystyle \frac{1}{{N}_{str}}\displaystyle \sum _{s=1}^{Nstr}gEUD\left(s\right)\\ \,=\,\displaystyle \frac{1}{{N}_{str}}\displaystyle \sum _{s=1}^{{N}_{str}}{\left(\displaystyle \frac{1}{{V}_{s}}\cdot \displaystyle \sum _{k=1}^{{M}_{s}}{v}_{k,s}\cdot {\left({D}_{k}\right)}^{{a}_{s}}\right)}^{1/{a}_{s}},\end{array}\end{eqnarray} \tag{ 4 }$$

where gEUD(s) is the gEUD value of OAR s of in total N_str OARs. V_s is the volume and M_s is the number of voxels of OAR s. k is a voxel index with volume v_k,s overlapping with OAR s. D_k is the dose of the plan in voxel k.

The dependency of the computational efficiency on the single parameter values is evaluated by the photon and electron beamlet computation time and the optimization time and memory usageof the H-DAO.

Additionally, the electron contribution is evaluated for the mixed beam plans. It is calculated by the mean dose to the PTV delivered by electron beams divided by the mean dose to the PTV of the whole plan.

One more treatment plan (selection plan) is created with the same computer hardware for each case listed in table 3 with the following combined selection of parameter values: u_x = 5%, u_e = 15%, t = 0.1%, d_m = 1 cm, d_l = 2 cm. This plan is compared to the conservative plan with the most conservative values under consideration: u_x = 2.5%, u_e = 5%, t = 0.0%, no voxel merging. The parameter values of the selection plan aim to achieve a reasonable compromise between degradation in final plan quality and computational efficiency. They are motivated by the results of the parameter study described in section 2.2.

The dependency of the evaluation quantities on the values of the investigated parameters, namely statistical uncertainty of the photon beamlets and electron beamlets, fractional sparse dose threshold and distances for the voxel merging are shown in figures 3–6, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For each investigated case with photons, the statistical uncertainty of the photon beamlets has large impact on the final treatment plan quality as shown in figure 3. However, the degradation is low for 5% compared to 2.5% statistical uncertainty. The two cases which are studied with mixed beams are less affected. The increased electron contribution might explain the smaller impact, because electron beams can be used instead of photon beams to a certain degree. Beamlet computation time for the photon beamlets is substantially reduced by increasing statistical uncertainty from 2.5% up to 7.5%, but for higher statistical uncertainties the computation time converges quickly. More surprisingly, the optimization memory usageand optimization time is also substantially reduced with higher statistical uncertainty of the photon beamlets except for the breast case, where these evaluations are probably dominated by the electron beamlets.

The statistical uncertainty of the electron beamlets has smaller impact on the final treatment plan quality compared to the photon beamlets as shown in figure 4. However, there are still small discrepancies visible, especially in the normal tissue for the brain case, which is caused by the increased photon contribution. Like for the photon beamlets, beamlet calculation and optimization time and also optimization memory usageis substantially reduced by increasing the statistical uncertainty for the electron beamlets.

Any of the investigated fractional sparse dose thresholds leads to degradations in final treatment plan quality, especially sparing of OARs and normal tissue in general as shown in figure 5. As expected, the beamlet calculation time is not substantially influenced but the optimization time gets reduced by at least 30% by using any of the investigated fractional sparse dose thresholds. The memory consumption of the optimization is also substantially reduced by at least 70% except for the prostate case only about 20%. However, increasing the fractional sparse dose threshold to more than 0.1% does only lead to minor reduction of optimization time and also memory usage.

Surprisingly, any of the investigated voxel merging distances did not lead to substantial degradation in final plan quality (see figure 6). Beamlet calculation time is as expected also not substantially affected, but optimization time and memory usageare enormously reduced with every shortening of the voxel merging distances. These observations are consistent over all cases.

The evaluation quantities are compared between the conservative and selection plan in table 4. Furthermore, DVHs between the two plans are compared in figure 7. According to the DVH comparison, dose homogeneity is nearly identical between conservative and selection plan for all cases, while for the selection plan the sparing of the OARs is slightly worse for some cases (${\overline{gEUD}}_{OARs}\,$is maximally 0.7 Gy higher). ${\bar{D}}_{NT}\,$is similar between the conservative plan and selection plan for all cases, except for the brain case for which the selection plan has a 0.4 Gy higher value.

Table 4. Evaluated quantities of the conservative plan with the most conservative parameter values and the selection plan with selected parameter values to achieve a reasonable compromise between degradation in final plan quality and computational efficiency. The last column depicts the relative change from conservative to selection plan.

	Conservative plan	Selection plan	Relative change
Prostate, IMRT
Objective function value	0.068	0.079	+15.4%
HI	0.087	0.094	+8.2%
${\overline{gEUD}}_{OARs}$	27.4 Gy	27.6 Gy	+0.8%
${\bar{D}}_{NT}$	1.4 Gy	1.4 Gy	+1.0%
Photon beamlet calc. time	0.4 min	0.2 min	−52.1%
Optimization time	19.3 min	0.7 min	−96.4%
Optimization memory	14.8 GB	0.4 GB	−97.3%
Lung, VMAT
Objective function value	0.136	0.150	+9.8%
HI	0.088	0.082	−6.8%
${\overline{gEUD}}_{OARs}$	4.4 Gy	5.1 Gy	+17.4%
${\bar{D}}_{NT}$	3.6 Gy	3.7 Gy	+2.7%
Photon beamlet calc. time	20.5 min	8.6 min	−58.1%
Optimization time	345.1 min	10.8 min	−96.9%
Optimization memory	20.7 GB	0.6 GB	−97.7%
Academic, MERT
Objective function value	0.271	0.276	+1.7%
HI	0.164	0.167	+1.8%
${\overline{gEUD}}_{OARs}$	24.4 Gy	24.3 Gy	−0.4%
${\bar{D}}_{NT}$	1.2 Gy	1.2 Gy	+0.3%
Electron beamlet calc. time	36.2 min	4.5 min	−87.4%
Optimization time	9.6 min	0.3 min	−97.3%
Optimization memory	41.2 GB	0.8 GB	−98.0%
Breast, MBRT
Objective function value	0.252	0.260	+3.4%
HI	0.125	0.122	−2.4%
${\overline{gEUD}}_{OARs}$	10.6 Gy	10.7 Gy	+1.4%
${\bar{D}}_{NT}$	5.1 Gy	5.2 Gy	+2.0%
Photon beamlet calc. time	3.1 min	1.3 min	−56.8%
Electron beamlet calc. time	1120.1 min	131.3 min	−88.3%
Optimization time	47.8 min	3.8 min	−92.0%
Optimization memory	238.2 GB	10.3 GB	−95.7%
Electron contribution	18.8%	20.2%	+7.3%
Brain, DYMBER
Objective function value	0.069	0.074	+6.6%
HI	0.086	0.086	0.0%
${\overline{gEUD}}_{OARs}$	18.0 Gy	17.9 Gy	−0.7%
${\bar{D}}_{NT}$	7.0 Gy	7.4 Gy	+5.7%
Photon beamlet calc. time	8.4 min	3.3 min	−60.7%
Electron beamlet calc. time	151.6 min	17.4 min	−88.5%
Optimization time	98.6 min	11.2 min	−88.7%
Optimization memory	803.7 GB	38.2 GB	−95.2%
Electron contribution	39.7%	26.9%	−32.2%

Zoom In Zoom Out Reset image size

Performance improvements from the conservative to the selection plans are substantial and consistent over all investigated cases. The absolute numbers vary enormously between the different cases, because of the different number of fields. In more detail, computation time for photon beamlets is about 52% to 60% reduced and for electron beamlets about 88% reduced. For the selection plans, beamlet calculation time for the photons is between 0.2 and 8.6 min and for the electrons between 4.5 and 131.3 min. Optimization time and memory consumption are even more reduced for the selection plans by about 88% to 97% and by about 95% to 98%, respectively. In absolute units, optimization time and memory consumption are 0.3 to 11.8 min and 0.4 to 38.2 GB, respectively.