A deep learning-based dose calculation method for volumetric modulated arc therapy

Abstract

Background

Volumetric modulated arc therapy (VMAT) planning optimization involves iterative adjustment of numerous parameters, and hence requires repeatedly dose recalculation. In this study, we used the deep learning method to develop a fast and accurate dose calculation method for VMAT.

Methods

The classical 3D UNet was adopted and trained to learn the physics principle of dose calculation. The inputs included the projected fluence map (FM), computed tomography (CT) images, the radiological depth and the source-to-voxel distance (SVD). The projected FM was generated by projecting the accumulated FM between two consecutive control points (CPs) onto the patient’s anatomy. The accumulated FM was calculated by simulating the movement of the multi-leaf collimator (MLC) from one CP to the next. The dose, calculated by the treatment planning system (TPS), was used as ground truth. 51 head and neck VMAT plans were used, with 43, 1 and 7 cases as training, validation, and testing datasets, respectively. Correspondingly, 7182, 180 and 1260 CP samples were included in the training, validation, and testing datasets.

Results

This presented method was evaluated by comparing the derived dose distribution to the TPS calculated dose distribution. The dose profiles coincided for both the single CP and the entire plan (summation of all CPs). But the network derived dose was smoother than the TPS calculated dose. Gamma analysis was performed between the network derived dose and the TPS calculated dose. The average gamma pass rate was 96.56%, 98.75%, 98.03% and 99.30% under the criteria of 2% (tolerance) -2 mm (distance to agreement, DTA). 2%-3 mm, 3%-2 mm and 3%-3 mm. No significant difference was observed on the critical indices including the max, mean dose, and the relative volume covered by the 2000 cGy, 4000 cGy and the prescription dose. For one CP, the average computational time of the network and TPS was 0.09s and 0.53s. And for one patient, the average time was 16.51s and 95.60s.

Conclusion

The dose distribution derived by the network showed good agreement with the TPS calculated dose distribution. The computational time was reduced to approximate one-sixth of its original duration. Therefore the presented deep learning-based dose calculation method has the potential to be used for planning optimization.

Introduction

Volumetric modulated arc therapy (VMAT), as one of the most promising radiotherapy technologies emerged in recent decades, has the potential to improve the treatment quality and the delivery efficiency [1]. During the VMAT delivery, the aperture shape (formed by multi-leaf collimator, MLC) and the beam intensity were modulated with the precise synchronization of gantry rotation. Correspondingly, the VMAT optimization is challenging because of the large optimization size and complex constraints. In order to get a deliverable VMAT plan, the leaf movements between adjacent gantry angles cannot exceed the product of maximal leaf speed and gantry rotation time. During planning optimization, numerous parameters are iteratively adjusted and the dose is recalculated [2]. Therefore a fast and accurate dose calculation method is the prerequisite of VMAT optimization.

The early dose calculation methods are referred to as correction-based methods [3]. The dose is calculated by interpolating or extrapolating based on the basic measurements in water. The empirical calculation is fast, but the calculation is inaccurate especially in the heterogeneous region [4]. Later developed model-based methods [5,6,7] convolute the primary photon energy fluence with a kernel that describes the contribution from scattering photons and electrons. The accuracy is significantly improved but at the cost of increased computational time. The principle-based method, commonly known as Monte Carlo method [8, 9], following the fundamental physics principle, simulates the radiation transport and energy deposition of individual particles. With enough simulated events, it can be used as the benchmark to verify the accuracy of other dose calculation algorithms. Correspondingly, the computational speed of Monte-Carlo (MC) is slow and the process is very time-consuming.

In order to balance the computational time and accuracy, modern commercial treatment planning systems (TPS) use the model-based dose calculation methods in the clinic, like the analytical anisotropic algorithm (AAA) (Varian Medical System) and the collapse cone convolution (CCC) algorithm (Pinnacle System). But still, the computational time is rather long to calculate the intermediate results during the optimization process due to the repeatedly adjusted planning parameters. So the simplified dose calculation method is used, like the singular vector decompose method of Pinnacle, which may cause the discrepancy between the intermediate dose and the finally calculated dose.

The recently emerged deep learning technique has provided an alternative “data-driven” approach for dose calculation. The initial application was to “denoise” the gross dose distribution calculated using the model-based method, so as to approach the accurate dose distribution calculated using the MC simulation [10]. Recently, Fan et al. [11] and Kontaxis et al. [12] developed the deep-learning based dose calculation methods for the static IMRT, which used the CT and fluence map (FM) as input and direct yielded the dose distribution. The network itself could be interpreted as a dose engine. These studies proved the feasibility of using deep learning for dose calculation. The network-based dose engine saved the trouble of modeling the complicate physics process, and also significantly increased the computational speed. However, these methods were developed for the static IMRT. For static IMRT, the FM was simple and straight, it was directly formed by the segment. But for VMAT, the MLC was modulated in the dynamic mode. This issue has not been addressed. Another issue was the divergence effect when projecting the FM onto the 3D volume, which may cause the omission of certain voxels within the 3D volume.

In this study, we proposed to expand the deep learning-based dose calculation method from IMRT to VMAT. The FM between the sequential control points (CPs) was accumulated to handle the dynamic MLC mode. A new projection method was developed to address the issue of divergence effect. The objective was to achieve rapid computation of the intermediate dose distribution during the process of VMAT planning optimization.

Methods and materials

Patient dataset

51 head and neck (H&N) VMAT plan data were used, with 43, 1 and 7 cases as training, validation, and testing datasets, respectively. The VMAT plans were optimized using the direct machine parameter optimization (DMPO) method of Pinnalce system. The MLC coordinates and MU at each CP were exported. The FM was derived by accumulating the fluence of every two sequential CPs. The corresponding dose was also exported. Each plan contained two full arcs. With 3°or 4°CP spacing, the number of CPs was 241 or 181, and the number of the FM and corresponding dose was 240 or 180, correspondingly. In total 7182, 180 and 1260 samples were included in the training, validation, and testing datasets.

Data preparation

FM accumulation

Figure 1 illustrated the FM accumulation between two sequential CPs with a simple example. The example supposed the leading MLC was static, and the tracking MLC moved from the position plotted in solid line to the position in dashed line. With this example, it was easy to deduct more complicate cases. Since the magnitude of the MLC movement was marginal, we supposed the MLC moved in a uniform speed for the sake of simplicity. The intensity at any position (x) was calculated as:

$$I(x)=\left\{ {\begin{array}{*{20}{c}} 1&{{x_2} \leqslant x<{x_3}} \\ {\frac{{x - x{}_{1}}}{{x{}_{2} - {x_1}}}}&{{x_1} \leqslant x<{x_2}} \\ 0&{x<{x_1},x \geqslant {x_3}} \end{array}} \right.$$

(1)

Fig. 1

Illustration of FM accumulation

Figure 2 showed the accumulated FM and its corresponding dose of the CPs at 181°, 185°and 189°gantry angle. The MLC position was exported from the dicomRT plan file. The leaf gap and transmission were also considered, which were set to 0.25 mm and 1% according to the commissioning data. The corresponding dose were shown in coronal view. Due to the leaf transmission, low dose within the body was observed at the region where the MLC was closed.

Fig. 2

Examples of the accumulated FM and corresponding dose. The corresponding dose were shown in coronal view

FM projection

As shown in Fig. 2, the accumulated FM was in 2D format. But the desired dose distribution was in 3D format. It was difficult for the DL network to use the 2D FM as input and to derive the 3D dose distribution. The 2D FM needed to be projected onto 3D volume representing the patient body for further proceeding. That was to assign corresponding FM intensity to the 3D volume. The commonly used 3D-DDA [13] or Bresenham [14] traversal methods first constructed the vector connecting the source point and the FM pixel, and then calculated the intersection voxels within the 3D volume. These methods were sensitive to the resolution of FM due to the divergence effect. The divergence effect was significant especially for the FM pixels far away from the center. If the resolution of FM was coarse, the divergence effect may cause the omission of certain voxels within the 3D volume. In order to handle this issue, we constructed the vectors by connecting the source point and the voxels. The intersection points on the FM plane were calculated. And the voxels were assigned to the value of corresponding intersection point. The comparison of the 3D-DDA and Bresenham methods against our method was shown in Fig. 3.

Fig. 3

FM Projection. 3-a and 3-b ploted only one row of voxels and one row of FM pixels for clarity. 3-c showed our method on a patient case. One coronal layer of the 3D volume was shown. The example voxel was plotted with blue dot, and the intersection point and source point in orange and red dots

Model input

The input of the network included the projected 3D FM, CT, radiological depth and the source to voxel distance (SVD). All inputs were in 3D format. Figure 4 showed the inputs at one axial layer. The CT value was converted to Hounsfield unit (HU) value. The radiological depth was calculated using the ray-tracer method proposed by Siddon [15]. The SVD was calculated as the magnitude of the vector connecting the source point and the voxel with 3D volume. The range of 3D FM and HU was 0 to 1 and 0 to around 2, respectively. In order to keep the data consistence, the radiological depth and SVD was normalized by 10 and 100, respectively.

Fig. 4

Network input

Model architecture

We used the classical 3D UNet [16] in this study. The architecture was shown in Fig. 5. UNet was currently one of the most widely used networks in the field of medical image processing. The network included encoder and decoder modules. The encoder module sequentially performed convolution, rectified linear unit (Relu), and down-sampling operations, and the decoding module sequentially performs up-sampling, convolution, and Relu operations. Additionally, UNet useed skip (copy and crop) connections to fuse the features at different scales to improve the network’s performance.

Fig. 5

Network architecture. The number of channels was denoted at the bottom

Training and evaluation

All plans were delivered on the Varian Novalis linac. The width of MLC was 5 mm at the center, and 10 mm at the top and bottom. The resolution of accumulated FM was set to 2.5 mm. The dose was calculated using the adaptive convolution (AC) method of Pinnalce system. The resolution was set to 4 mm×4 mm×4 mm.The dimension was cropped to 128 × 80 × 80, which was wide enough to cover the regions of interest (ROIs) of all patients. The network inputs were also sampled and cropped to align with the dose volume. The FM accumulation, projection, radiological depth and SVD calculation were implemented using Matlab software. The network training was implemented with pytorch [17] on a desktop computer with Intel i9-11900 K processor and NVIDA GeForce RTX 3090 GPU. The batch size was set to 10. The learning rate was set to 1.0e-04 with 1.0e-6 weight decay.

Results

Qualitative evaluation

Figure 6 showed the comparison of the TPS calculated and network predicted dose at one single CP. The positions of the comparison profile were selected to intersect the region with steep dose gradient. In general, the dose profiles coincided with each other. The calculated dose profile showed a sharper peak, and the predicted dose was smoother. This was because the major operation of the network was convolution operation, which tended to produce a smoothed result.

Fig. 6

Comparison of TPS calculated and predicted dose at one single CP. The unit of the color bar in figure a-b was cGy. The locations of the horizontal and vertical profiles were annotated in red and green dash lines

Figure 7 showed the comparison of the total dose of one testing patient. The positions of the comparison profiles were also selected based on the dose gradient. Figure 7-c showed the isodose line of the absolute dose difference. For the major regions, the difference was lower than 100 cGy. The discrepancy greater than 200 cGy was observed at only two small areas at the edge region. The profile comparison showed the same trend as the single CP: two profiles coincided with each other well, expect the predicted profile was smoother. This also explained the target curves (PTV1, PTV2 and GTVnd) of the calculated dose were slightly steeper than the predicted curves. Expect this, the rest DVH curves displayed good coincidence.

Fig. 7

Comparison of TPS calculated and predicted total dose. The unit of the color bar in figure a, b and c was cGy. The locations of the horizontal and vertical profiles were annotated in red and green dash lines. Figure f plotted the curves of calculated and predicted dose in solid and dash lines, respectively

Quantitative evaluation

Figure 8-a showed the gamma analysis of the calculated and predicted dose under different criteria. The calculated and predicted dose were normalized to the max dose of the calculated dose (Dcal max). The average gamma pass rate was 96.56%, 98.75%, 98.03% and 99.30% under the criteria of 2% (tolerance) -2 mm (distance to agreement, DTA). 2%-3 mm, 3%-2 mm and 3%-3 mm. Figure 8-b showed the relative error of the calculated and predicted dose. In order to keep consistent with the gamma analysis, the relative error was also normalized to Dcal max.

Fig. 8

Relative error and gamma analysis

Table 1 listed the max (D_max), mean dose (D_mean), the relative volume irradiated by the dose greater than 2000 cGy (V₄₀) and 4000cGy(V₄₀) of the critical ROIs. For the target, the relative volume irradiated by the prescription dose (V_DP) and the mean dose were listed. The average dose calculation time for the testing patients of the network and TPS was 16.51s and 95.60s. And the average time for one single CP was 0.092s and 0.530s.

Table 1 Comparison of critical indices. The unit of D_max and D_mean is in cGy, and the unit of V_DP V_{20 and} V₄₀ is %. The 95% confidence interval is also listed

Discussion

The motivation of this study was to develop a deep-learning based dose calculation method to calculate the intermediate dose distribution during planning optimization. The AC method, which was used to calculate the final dose distribution in the clinic, was used as the ground truth. This could avoid the dose discrepancy between the intermediate and final dose. However, as a model-based method, sharp peaks were observed on the AC dose profiles. On the other hand, the deep learning network tended to produce an averaged or smoothed results due to the basic operation of convolution. This contributed the main difference between the predicted and calculated dose. Using the MC simulation as ground truth, or a finer dose grid may reduce the discrepancy. With improved accuracy, the dose calculation method has the potential to be used for final dose calculation.

Another limitation of this work was the dataset lacks of diversity: all plans were delivered on the same linac. The network training feed the network with the compulsory data for dose calculation, and then minimized the difference with the TPS calculated dose. The physical principle for dose calculation was acquired by the network in the training process. In other words, the training process could be interpreted as the commissioning process of the model-based method. Compared to the commissioning process, the training process is more generic. Because the basic physical principle is universal, it is easy to apply the well-trained network to different linac using the transfer learning strategy. In addition, all training data were head and neck plans. The trained model was not validated on other sites. The accuracy may be acceptable for the abdomen plans, but not for the thorax plans due to the heterogeneous tissues. But we would like point out that the main contribution of this work was to expand the deep learning-based dose calculation method from IMRT to VMAT. This study proved the feasibility for head and neck VMAT plan. In the future, we plan to train and validate the model using more data from various sites.

The difficulty of the deep learning-based VMAT dose calculation method is how to get the paired sample of accumulated FM and corresponding dose, i.e. the input and labeled ground truth for the network. In this study, we simplified the FM accumulation process by supposing the MLC moves linearly between adjacent CPs. In addition, the gantry rotation between adjacent CPs was also ignored. In practice, we found that the CP dose was only stored in memory. It would be lost once log out. So the dose was recalculated before exporting the CP dose. The results demonstrated that the simplification did not affect the accuracy. But for greater gantry gap, this should be taken into consideration. On the other hand, one VMAT plan contained around 180 paired samples. So limited patient cases were sufficient for the network training.

Conclusion

The proposed learning-based dose calculation method reduced the computational time to one-sixth of the TPS without substantial loss of accuracy. It has the potential to be used for VMAT planning optimization.