The physical 3DCM has two basic components: a non-uniform SOBP modulator customized to the various thickness of the target volume along the beam axis, and a range-modulating compensator customized to the target contour at the distal end (Figure 1).
To demonstrate the principle and feasibility, a theoretical mono-energetic proton beam with 184 MeV (no energy spread, Bragg peak at water equivalent depth of 22.5 cm) was used with Lucite being the material for the modulator, and water media for dose computation.
To generate the optimal geometry of 3DCM and the proton fluence map, mini pencil beamlets with a 0.01 cm radius, rather than actual clinical pencil beams are used. The dose kernel of the mini pencil beamlet was generated by the Monte-Carlo (MC) simulation program (FLUKA version 4-0.0) (6-9) using the PRECISIO mode, which defaults particle transport threshold at 100 keV. The resultant dose kernel was stored in a matrix with 1 × 1 × 1 mm3 resolution.
2.1 Non-uniform SOBP modulator
The non-uniform SOBP modulator is an assembly of pyramid-like ridge filters. Each pyramid is formed by stacking cuboids with decreasing side length. A single-energy broad beam (to be formed by array of pencil beamlets), with the beam axis perpendicular to the pyramid base, would penetrate through different thicknesses of the material. The dose distribution along the beam axis is the integration of all sub-beamlets passing through various heights of the pyramids, forming a SOBP.
The geometry of the pyramid is derived through multiple iterations. Each pyramid with its unique shape is expected to produce a homogenous SOBP region that would not only cover the thickness of the target along the beam axis but also merge smoothly with the neighboring SOBPs produced by adjacent pyramids. The critical elements in the geometry of the pyramid array include the shape of the pyramids (base, height, lateral graduation) and the spacing between them.
2.1.1 Critical size of the pyramids and their spacing
To obtain the critical size and spacing of the pyramids, the mini pencil beamlets were used in hexagon formation, in which the mini pencil beamlets were located at the center and at the vertices of each hexagon, where the central axis of the pyramids are expected to situate. The critical spacing is defined and obtained by gradually changing the size of the hexagon until the composite lateral dose distribution at the depth of the Bragg peak achieves desired uniformity, with individual pencil beamlets dose pattern vanished. The obtained critical spacing was then used as a guide to determine the distance between the centers of each pyramid, which would also be the maximum side length of the pyramid base. Note that no pyramid was yet used in this calculation. The purpose of this procedure is to assure all beamlets going through various parts of a pyramid and through the neighboring pyramids will merge smoothly, when the side length of the pyramid base is equal or less than the critical spacing.
2.1.2 Geometry of the SOBP pyramid
The shape of a pyramid is formed by stacking cuboids of various thicknesses (levels), which leads to a spread out energy distribution of outgoing beams. Each level’s net area, a circumference of beam cross-section corresponds to the weighting factor of a specific out-going energy associated with the thickness of that level. The design criterion of pyramid levels is to produce a uniform plateau of SOBP. Accordingly, we generated another set of 40 dose kernels using the FLUKA software, each corresponding to a 184 MeV mini pencil beamlet (0.01 cm radius) passing through Lucite of thickness 0.2 cm to 8 cm with 0.2 cm increment. These dose kernels would be integrated with a certain derived weighting assignments to obtain the final depth dose curve in the direction of beam axis. The weightings of all levels (cuboids) of a pyramid were inversely calculated using a modified version of the basic gradient descent algorithm (10):
where θ is the array of weighting factors corresponding to each beamlet, that is forced to be greater than or equal to zero to ensure physical validity; α is a user-defined step size to ensure fast convergence to a local minimum of the cost function; x is an array that contains all the depth dose curves; y is the optimization goal, a set of the anchor points of an ideal SOBP curve. In this study, the anchor points were assigned to be 1 for the depths from 15 cm to 20.2 cm and 0 for the depths from 21.5 cm onward, and the initial value of θ were assigned to be zero. These assignments ensure the calculated weighting factors would yield a smooth plateau between depth 15 cm and 20 cm with a sharp distal dose fall-off. The setting would keep the entrance dose as low as possible. J is the cost function that tracks the difference between optimization results and the preset goal.
2.1.3 SOBP modulator and its design validation
To demonstrate the validity of the described method, we constructed a uniform SOBP modulator with a cross-sectional size of 5 × 5 cm2, formed by arranging individual pyramids abutting each other side by side with each row offset by half of a pyramid. The side length of the pyramid base was set to the critical spacing described in 2.1.1. We then conducted MC simulation of the dose distribution in water produced by a 184 MeV broad proton beam with a field size of 4 × 4 cm2 penetrating through the modulator. The MC simulation was performed in the PRECISIO mode, with resolution of 1 × 1× 1 mm3. The resultant percentage depth dose along the central axis was compared with the intended shape of SOBP to validate the inverse weighting factors calculated by the modified version of the basic gradient descent algorithm. The lateral dose homogeneity of the SOBP plateau region was analyzed to validate the determination of the critical spacing and the side length of pyramid base.
2.2 Generation of 3DCM with optimized 3D conformal dose distribution
With the inverse algorithm for the optimal structure of SOBP modulator (pyramid-based) at hand, we proceed to combine energy and range modulations to produce 3D conformal modulator using a 184 MeV mono-energy proton beam.
The test target volume was a sphere of 4 cm in diameter. We first designed a virtual “mesh” perpendicular to the beam axis. The size of each cell of the mesh is 2 mm × 2 mm (less than the aforementioned critical pyramid spacing), and would serve as a placeholder for 2 mm x 2 mm pyramids. The rows of the mesh offset with each other by 1 mm. We created a new sphere by expanding the target sphere radius by one cell width. The entire mesh area was defined by the maximum cross-section of the expanded sphere.
Using the inverse optimization method described by Equation (1), the set of 40 beamlet dose kernels passing through Lucite of thickness from 0.2 cm to 8 cm with 0.2 cm increment were placed to the center of each cell. The dose kernels with their Bragg-Peaks located outside the expanded target sphere were removed. The remaining dose kernels constitute the parameter x. The array of the weighting factors θ, corresponding to each of the 40 dose kernels, would be used to define the graduating shape of each individual pyramid. The optimization goal y would be the anchor points representing the desired dose at specific locations, and would be adaptively adjusted based upon the optimization staging. The dose calculation uses a 2 mm × 2 mm × 2 mm grid spacing.
We developed an adaptive ring optimization technique that consists of three stages. The first stage focuses on the target coverage. Target coverage in this stage was defined as the percentage of the target volume covered by 90% isodose surface. The initial inputs for θ would be zero. The first optimization goal y would be set to one within the spherical target volume, without any constraint outside the target volume. During the optimization iteration, the weighting factor array θ resulting in the highest target coverage would be cached, and the iteration would stop once the target coverage converged.
The second stage aims to improve the conformity index by introducing adaptive ring constraints:
where Vtarget denotes the target volume. Vtreated denotes the treated volume, which is the volume enveloped by the prescribed isodose surface (90% in this study). The conformity index nCI is the product of two ratios, one being between the target volume and the intersectional volume of target and treated volume, and the other ratio being between the treated volume and the intersectional volume of target and treated volume. The conformity index defined by Equation (2) takes into consideration the intersectional volume of target and treated volume. When the target volume and the treated volume are completely overlapped, nCI = 1 is the optimal index value. In order to regulate the optimization, we created three sets of “rings” perpendicular to the beam axis at various depths surrounding the target sphere. Each set of rings contains ten rings with different diameters at ten depths. The first set of rings has diameter 10 mm larger than the target sphere’s cross-sectional diameter at the corresponding depths; the second and third set of rings have diameters 16 and 24 mm larger than the target cross-section, respectively. In addition to the original optimization goal of setting 100% target volume coverage, constraints on rings would be implemented into y. In the first iteration, the weighting factor θ computed from stage one was used to calculate the dose distribution. Both nCI and target coverage were assessed. The relative dose on each ring was evaluated. The dose-goal for the first set of optimization rings was set to be 10% less than the prescribed dose, the second set to be 6% less, and the third set to be 3% less. The gradient descent algorithm would then be used to optimize θ, the array of the weighting factors. During the first iteration, the θ corresponding to the highest nCI with the coverage greater than 95% was saved, and the iteration would stop once the cost function stabilized. After the first iteration, the optimization goals of all rings would then be adaptively adjusted based on the difference between the current dose distribution from the updated θ and that of the previous iteration. The program would subsequently proceed to optimize θ using gradient descent algorithm and would cache the optimal θ that yield the highest nCI while keeping the target coverage greater than 95%. The adaptive ring optimization would continue until the target coverage dropped below 95%. The cached θ resulting highest nCI would be the final optimal weighting factors.
The optimal θ values contain the weightings of all dose kernels that would result in the desired width of SOBP covering the target along that beam axis. The weightings of those dose kernels passing through the same cell on the virtual mesh would be used to reconstruct the pyramid located on that cell. Essentially the shape (height and graduation) of individual pyramid was determined by the values of θ on its corresponding cell. The sum of the weightings belonging to each cell was assigned to a new array, which represents proton fluence weight going through the pyramid located on that cell.
Figure 2 presents a flow chart of the aforementioned optimization process.
The 3DCM geometry and the proton fluence map would be imported into FLUKA to compute the dose distribution in water using MC simulation (in the PRECISIO mode, with resolution of 1 × 1 × 1 mm3). The nCI and the target coverage would then be compared with those calculated by the inverse optimization process.