1 Introduction
One of the main advantages of proton therapy is the depth dose profile, which has been named Bragg peak. The few-mm-wide pristine Bragg peak allows to irradiate a well localized region within the body, keeping low the dose release in the proximal (entrance plateau) and distal (tail) regions. If the tumor cells are bigger than the width of the Bragg peak, superposition of several Bragg peaks from the different beam energies is introduced to make a uniform radiation dose over the tumor.
For pencil beam scanning, the target is irradiated layer by layer with the different beam energy. However, switching time between energy layers is a major contribution to the overall time. Few energy layers may result in large ripple produced on the flat-top depth-dose curve, which becomes serious for low energy. The ripple filter, which is of great interest in recent years[1‒5], is used to address this problem. Fig. 1 shows the vertical cross-section of the commonly used ripple filter. After inserting the ripple filter perpendicular to the beam, the Bragg peak of the proton beam will become smooth[1,3].
Although Weber and Kraft[1] have presented a least square (LS) objective function to seek the width of the ripple bar, they have not mentioned the method to solve it. In Ref.[2], iterative-least-square technique (ILST) is introduced to optimize the weights. In this paper, a new quadratic objective function is proposed for designing the width of the ripple filter bar. Besides, the ILST is also used to optimize the LS function for the comparison.
2 Optimization methods
The purpose of the treatment planning is to provide a desired biological dose which is calculated via multiplying the physical dose by the relative biological effectiveness (RBE). As recommended by the ICRU[6], RBE is set to 1.1 in this paper.
The principle of the analytical computation for the ripple filter design used in Ref.[1] is adopted in this study. When the proton beam passes through the ripple filter with the thickness tr, the dose distribution will be shifted to a certain water-equivalent value that corresponds to tr along the beam direction. Therefore, when the proton beam passes through the ripple filter perpendicularly to the beam direction, the pristine dose Bragg peak will be transformed by the ripple-bars into a superposition of displaced Bragg curves. Gaussian Bragg peaks were chosen to generate homogeneous depth dose profiles with a low ripple for the smoothing effect after superposition[1,3,5]. To form the dose profile close to Gaussian shape in the peak region, the weight of proton component has to be calculated. Since the weight of the proton component passing through the ripple-bar step is proportional to the width of the ripple-bar step, the weight determines the shape of the ripple bar.
Additionally, the difference of proton fluence caused by the nuclear interaction in the ripple filter and water is taken into account for the computation. After proton beam passing through the ripple filter, the biological dose is modulated as follows:
where dmod(zi) represents the modulated biological dose at a depth zi; ∆z is the step size of the dose points; λ is the period of the groove structure; jmin and jmax are defined as the value of the minimum and maximum water-equivalent thickness of ripple bar divided by ∆z, respectively; dphy(zi+j) is the pristine physical dose at a depth zi+j; wj is the weight of the jth segment of proton beam; fj , the fluence difference introduced by the jth step of ripple-bar, as described in Eq.(2)[7]:
where tj=t(zj) = 1.19tr(zj) is the water equivalent thickness of ripple filter in z direction; λrf = 69.54 cm is the nuclear interaction length of plexiglass; mp = 0.012 cm-1; R is the proton range in water.
To simplify the analysis, (1) can be further written as the following matrix form:
where Dmod is modulated biological dose; aij = 2Δz·dphy(zi+j) fj /λ; m is the dimension of Dmod; n=jmax – jmin+1, which is the dimension of weight W; A is a m×n matrix; W is a n×1 vector and wj≥0. Since it is often the case that m>n, then (3) is an over-determined equation thus there might be no unique solution. Besides, the coefficient matrix A is usually ill-conditioned so that the solution of (3) is derived by solving an optimization problem. In this paper, two optimization methods with different objective function are used and compared for the ripple filter bar design.
2.1 Iterative least-squares technique
ILST is a commonly used method for solving the least square objective function[8‒10] with relatively fast converge rate and good performance of averaging out the noisy data. In hadron therapy range modulation field, Schaffner et al.[2] utilized the ILST to produce the large, biologically uniform spread-out Bragg peak (SOBP) depth dose from the small SOBPs. The weight updating in two consecutive iterations (kth iteration and (k−1)th iteration) can be calculated by
where g is set to 1 inside of the SOBP and 0 outside; P(zi) is the prescribed biological dose at a depth zi.
2.2 Quadratic programming method
A large condition number of the matrix A, i.e., larger than 100, may result in the different order of magnitudes of the derived weights value wj. This will lead to the manufacture difficulties of the ripple filter because the normalized value of wj is proportional to the width of the jth ripple bar. For example, if wj+1 is 100 times smaller than wj and the width of the jth segment is 10 μm, then the required width of the (j+1)th segment is 0.1 μm which might beyond the manufacture precision (1 μm). To prevent the weights from growing too large, a new objective function is proposed in this study, it consists of a least square item and an extra weight decay item[11].
where ‖·‖ represents Frobenius norm operation; P (m×1) is the prescribed biological dose; γ (γ ≥ 0) is a parameter that determines how strongly large weights are penalized. The larger γ is, the more uniform weights are. Particularly, when γ = 0, Eq.(5) becomes
The problem (6) is a special case of (5). Problem (6) is a least square optimization problem with nonnegative constrains. Problem (5) can be written in the standard quadratic form as follows:
where I is a n×n identity matrix. Problem (7) can be solved by a quadratic programming method (QPM) [12]. Adding the weight decay term in ripple bar design can improve generalization[11].
2.3 Performance evaluation
The optimized results with ILST[2] and QPM are compared in this paper. Moreover, the robustness (anti-noise ability) of two methods which specified as the dose difference was also examined.
where dn and d0 represent the dose derived from the noise data and original data, respectively; r is the relative dose variation, which represents the robustness of the filter. The smaller r is, the stronger robustness is.
3 Results and discussion
In this study, a typical beam delivery system for proton therapy is simulated with Geant4. Bragg peaks obtained from the simulation are used to evaluate the methods. Since the Bragg peaks of proton beams becomes narrow with the decrease of energy, the relative low energy of 70 MeV is chosen in the simulation.
3.1 Simulated model
As shown in Fig. 2, the schematic of the simulated model consists of a kapton window located at the end of the beam transport line, a nozzle, and a water tank phantom. The components of the nozzle are scanning magnets, vacuum chamber, beam monitors, a ripple filter and etc. The scanning magnets are used to steer the beam. And the vacuum chamber is designed to reduce the beam scattering. The beam monitors, are ionization chambers with an equivalent water thickness of 1.1 mm. The ripple filter is used to smooth the flat-top region of the SOBP. The distance between the isocenter and the ripple filter is 50 cm. The momentum spread for protons is set as Δp/p=0.05%.
3.2 SOBP dose uniformity with ripple filter designed by QPM
Fig.3 shows the dose uniformity with and without ripple filters for proton beams of 70 MeV. As illustrated in Fig.3(a), for the same scanning step (3.2 mm), the dose uniformity is as high as 24.7% without ripple filter and the dose uniformity sharply falls down to less than 2.5% when a 3 mm ripple filter was added to the system. In Fig.3(b), the dose uniformity is more than 20.7% without ripple filter, while with ripple filter the maximum dose scanning error decreases to less than 2.5% with the same step (3 mm).
As described in Ref.[3], the widths of distal fall-off (80%–20% distal falloff) are both moderate (less than 3 mm).
3.3 Comparison between ILST and QPM
Fig.4 illustrates the simulation results with ILST and QPM when pencil beam scanning step is 2.14 mm (a random value). As shown in Fig.4(a), the depth dose distribution with QPM (solid line) is more uniform and closer to the ideal value (1.00 in the figure) than that with ILST (dashed line). Particularly, the maximum dose uniformity with QPM is 0.5%, which is less than that with ILST 0.7%. In addition, the shape of ripple filter with QPM is closer to triangle than that with ILST as shown in Fig.4(b). The reason is that the extra weight decay item prevents the value of the individual weight which corresponds to the width of ripple filter bar from growing too large. Therefore, the width of ripple filter bar tends to have the same order magnitude. Therefore, its corresponding ripple filter bar is close to the triangle shape and so that is convenient for manufacture.
In the real application, the dose is measured in the experiment instead of deriving by the simulation. Then the unavoidable measurement error of dose should be considered. The noise with Gaussian distribution was added to the simulation dose in this study. Fig.5 shows the relative dose variation and the ripple filter shape curve difference while the 5% noise is added in the simulation. The dose difference value with QPM and ILST is shown in Fig.5(a). The calculated standard deviation of the relative dose uniformity with QPM (0.0039) is smaller than that with ILST (0.0042). Fig.5(b) further compares the filter shape difference, the horizontal axes represents the shape difference with and without noise data. The smaller values derived by QPM indicates its stronger robustness. In other words, the proposed QPM is capable of obtaining the similar ripple filter even using the data with noise, which is preferred in the real application. The good robustness is due to good generalization ability of the QPM. The results of scanning step of 3.2 mm with the use of 3 mm ripple filter are similar.
Fig.6 compares the individual proton beam Bragg curves and the corresponding lateral dose distribution at different positions with and without ripple filter. The relative dose of the modulated curves are higher than that of the pristine Bragg curve at the depth of modulated Bragg peak as illustrated in Fig.6(a). This phenomenon is also reflected from Fig.6(b), in which the lateral dose distribution of proton beam without ripple filters are both lower than that with ripple filters. Besides, the curves for A and B is close to each other in Fig.6(b), which indicates there is slight difference of lateral dose distribution at Bragg peak between two different ripple filters. This is due to the proton beam is scattered by different ripple filters mainly according to Moliere theory[13,14].
4 Conclusion
A new design approach for ripple filter bar was proposed in this paper. The new approach uses the new objective function with QPM. It has been successfully applied to calculate the weights, which were further converted to the ripple filter dimension. Compared to the ILST, the proposed method not only achieves better performance, i.e., better depth dose uniformity and ripple filter’s shape for manufacture, but also has better robustness characteristics. This ripple filter can be manufactured more accurately and achieve a good flatness of the SOBP at the same time.
