MU-Tomo: Independent dose validation software for helical tomotherapy

A software program, MU-Tomo, has been developed to perform an independent point dose calculation and compare it to the dose calculated from the TomoTherapy (TomoTherapy, Inc., Madison, WI) treatment planning system (TPS). Input parameters required for this software include: archived tomotherapy patient files, QA plan image coordinates, tomotherapy-calculated point dose and machine-specific dosimetric parameters such as the off-axis ratios (OARx and OARy), tissue phantom ratios (TPR) and output functions (Scp). The software was validated on four phantom models and fifty tomotherapy patient plans representing various anatomical sites. Our results indicate that MU-Tomo can perform in a few seconds an independent dose calculation accurately and provide a secondary check for a point dose validation of helical tomotherapy plans.


Introduction
The TomoTherapy ® Hi-Art ® unit (TomoTherapy, Inc., Madison, WI) is able to deliver intensity modulated radiation therapy (IMRT) via a helical trajectory about the patient using a fan beam irradiation geometry (Mackie et al., 1993;Mackie, 2006). Helical radiation delivery is simultaneous with couch translation and gantry rotation and has certain advantages in delivering conformal dose distributions (Zacarias et al., 2006;Mackie et al., 2003;Beavis, 2004;Yang et al., 1997). The American Association of Physicists in Medicine (AAPM) recommends that doses to cancer patients need to be carefully designed and verified (Kutcher et al., 1994). Therefore, especially for IMRT, it is common to have in addition to a computer generated treatment planned a patient specific quality assurance (QA) measurement and an independent monitor unit or dose calculation to verify the treatment plan.
For independent monitor unit or dose calculation, several studies have been performed for different modalities in the years past. Ayyangar et al. (2001a, b; reported on their independent dose calculations based on multi-leaf collimator (MLC) for Corvus ® , Peacock ® , and ADAC ® systems from 2001 to 2003. Beck et al. (2004) developed an algorithm for an independent verification of Gamma Knife ® treatment plans. Chen et al. (2002) provided an independent monitor unit calculation using the MIMiC™ multileaf collimator. The MLC log information was also studied as a means to perform second checks and has been reported by several research groups (Chen et al., 2005;Lorenz et al., 2007Lorenz et al., , 2008; Luo et al., 2006;Yang et al., 2003;Zhu et al., 2003). Georg et al. (2007a, b) developed a fluence-based dose calculation software MUV for independent dose verification in IMRT planning. Jansen et al. (2003) imported the planning parameters from the CadPlan ® treatment planning system into the Pinnacle 3 treatment planning system as a means to perform an independent second check.
Another approach for an independent dose calculation can be established using a Monte Carlo simulation method. Edimo  In this study, we report on the development of in-house software called MU-Tomo which performs an independent point dose validation of helical TomoTherapy ® treatment plans. Our method is an extension of that proposed by Gibbons, et al. (2009). Our implementation is an analytical, correction-based second check that was developed on the Matlab platform. The input data needed for the calculation include the archived patient plan and dosimetric functions, such as off-axis ratio (OAR x and OAR y ), tissue phantom ratio (TPR), and output function (S cp ). Innovations in our method include the averaging of OAR x and OAR y over the positive and negative side of the profile axis, the optimization of OAR y to minimize fluctuations, the application of Mayneord factor to correct for SSD changes in the percent depth dose, automatic depth calculation to the point of measurement (obtained from the patient specific QA calculation in the Cheese Phantom), and automated retrieval of gantry starting angle, delivery duration, initial position of IEC-Y of tumor setting, field size and pitch from the archived documents. Furthermore, in our implementation the point of measurement can be located anywhere inside the tumor as mapped to the Cheese Phantom and is not limited to be at the geometric center of the phantom.

MU-Tomo software framework
The purpose of this in-house developed independent dose calculation software is to calculate quickly a point dose and compare the calculated dose with the planned patient dose provided by the TomoTherapy ® TPS. The software has three input components: the archived patient documents, the initial set up coordinates and dose to the QA point, and the machine-specific dosimetric functions. The archived patient files provide core information for the patient's plan such as the initial gantry start angle, treatment time, pitch, initial IEC-Y position of the target, field width size and end of planning (EOP) sinogram. These parameters are retrieved automatically from the archived documents through the Extensible Markup Language (XML) file associated with each patient. The initial gantry start angle, initial IEC-Y position of the target and sinogram are used to localize the first projection. The TomoTherapy ® planning system determines the location of the first projection in a complex manner. Initially, the first possible projection starts with the gantry at zero degrees (up right) and couch at an IEC-Y position such that the beam is centered on the first slice of the planning CT. Subsequently, the planning station calculates an initial sinogram that covers the entire planning CT volume. During the optimization, the user's selection of target contours determines which of the possible projections are actually used. Because of the variation in selection and position of the target volume, the first projection does not necessarily start at zero degrees and this will solely depend upon the gantry angle where the edge of the field width first encounters the initial IEC-Y position of the target.
The Hi-Art ® unit has 32 pairs of multileaf collimator (MLC) leaves. The sinogram file, originally in hexadecimal format, is converted into a decimal sinogram map using an in-house code. The sinogram file shows the relative open time of each individual leaf as a function of leaf number and projection. By knowing the total plan delivery time, the absolute open time for every leaf can be calculated-an important parameter for the dose calculation during the segmentation of the sinogram.
Additional input parameters for the second check software include the initial setup coordinates and the tomotherapy-calculated point dose. For the current version of the TomoTherapy ® software (Version 3.1), these parameters must be acquired from the treatment planning system. MU-Tomo utilizes these parameters to calculate the point dose independently and compares the value to the tomotherapycalculated point dose. Additionally, dosimetric functions, such as the off-axis ratio along IEC-X and IEC-Y (OAR x and OAR y ) directions, tissue phantom ratio (TPR), and output function (S cp ), are embedded into the MU-Tomo software. Details regarding these functions are discussed in the next section. (Here Figure 2. A sample sinogram from a patient treatment plan).
Based on the input parameters discussed previously, the software starts a projection along the sinogram map. The Hi-Art ® system calculates treatment plans in helical mode by computing 51 static projections for each gantry rotation. Figure 2 shows an example of a sinogram map for a sample patient. The relative time intensities for all leaves in each projection provide the foundation for the segmentation of the leaf opening in the sinogram (Gibbons et al., 2009). Finally, the second check software executes the dose calculation with a series of mathematical combinations of dose rate, dosimetric functions, and time-dependent segmentations.

Dosimetric functions
Dosimetric functions were embedded into the software from the TomoTherapy ® unit commissioned data. Currently, commissioned data from two jaw settings, nominal 2.5 cm and 5.0 cm, were utilized in the software development. The lateral and longitudinal beam profile data (OAR x and OAR y ) was used to calculate the off-axis dose dependence. The scatter function, S cp , represents the output factors for different equivalent square field sizes using the MLC. For the tissue phantom ratio (TPR), values were calculated from the output factor, S cp and the percent depth dose (PDD) profiles. The PDD data was also acquired from the commissioned data tables. (Here Figure 3. Normalized, lateral beam profile, OAR x ).
A normalized, lateral beam profile, OAR x , for the 2.5 cm and 5.0 cm field width settings is shown in Figure 3 (A). In our software development, OAR x profiles at a depth of 10.0 cm were utilized for calculations. Moreover, the OAR x ratios were calculated using the   With the depth of the measurement point known, the OAR y value at the specific depth was identified by interpolation or extrapolation using the discrete measurement profiles. Ultimately, the OAR y was determined by both the off-axis distance in the longitudinal direction as well as the depth. This dual dependency is an important characteristic of OAR y and will be reflected in the final dose calculation formula. From Figure 3 (B) and (C), the sensitivity of the OAR y to field width size is shown. This sensitivity provides different OAR y values upon the definition of segmentation along the sinogram map. In the current software, the full leaf open field was used to perform an OAR y optimization and thereby resulting in a better dose differences in over 50 patient plans. (Here Figure 4. Output scatter factor function vs equivalent square field size).
The output factor, S cp , was determined by measuring the data at the depth of 10.0 cm for the two different field width settings investigated. S cp is field size dependent as shown in Figure 5, particularly for small field sizes. In the current software, the equivalent square field size was determined during segmentation, and the corresponding S cp was applied for the calculation. (Here Figure 5. Percent depth dose profiles).
The tissue phantom ratio (TPR) was not measured directly. The specific TPR values were calculated using the known PDD and S cp values ( Figure 5). The approximate relationship among them was obtained from [28]: where d denotes the depth of the point of measurement, r and r d denote the equivalent square field sizes at the surface and at 10 cm depth, respectively. and PDD (d,r,SSD=75) represents the profile with a SSD= 75 cm. Since our PDD measurements were acquired using SSD= 85 cm, the Mayneord factor was used to convert PDD profiles:   where SSD i represents the source to surface distance for each projection and d m is the depth with maximum PDD. The TPR depends on both the depth of measurement and equivalent square field size.

Segmentation of sinogram projection
The sinogram was converted to provide the open times for each MLC leaf. From the sinogram file, the time intensity distribution for the 64 leaves in any projection was found and shown to be mostly non-symmetrical. The goal of the sinogram segmentation was to bin multiple patterns with constant intensity in each pattern within any projection so that the superposition of these patterns produced the composite intensity modulation. The segmentation was accomplished by averaging the time distribution symmetrically around the central leaf position for each projection. The central leaf position was defined as the MLC leaf directly normal to the point of measurement. Subsequently, the outer-most leaf, which was identified as the leaf with the maximum distance from the central leaf, was determined. Starting from the outer-most leaf, the first evenly distributed segmentation pattern was determined. All leaves in the first pattern were assigned the same open time as the outer-most leaf. Gradually, the segmentation patterns moved from the first outer-most leaf to the central leaf. The open time for any given segmentation pattern was obtained by subtracting the previous segmented time from the current leaf open time. Once the iteration reaches the central leaf, the segmentation is complete.
Within each segmentation pattern, the equivalent square field size was determined. The S cp for each segmented pattern was obtained from the relationship between the equivalent square field size and S cp as discussed above. The software calculated the dose from all these segmentation patterns within one projection and accumulated the total dose from all the projections.

Dose calculation algorithm
The dose calculation formula used in our second check software was first proposed by Gibbons et al. (2009) and is shown in Eq. 3. The objective of this calculation is to utilize the dosimetric functions, archived patient documents and coordinates of the point of measurement along with the segmented sinogram to accumulate the dose to a point from all the segmented patterns of all projections.
Equation (3) shows how the dose at a point P is calculated based on the information entered into MU-Tomo. D 0 is the dose rate under normalization conditions and is measured at a depth of 10.0 cm in water with a field size of 40×5 cm 2 and SAD= 85.0 cm. SPD is the source to point of measurement distance. For each projection, the SPD was calculated as follows: where P x , P y , and P z are the relative coordinates from the point of measurement to the gantry axis; X i , Y i , and Z i are the relative coordinates from the source to the gantry axis. P x , P y , and P z are acquired from the initial input data, and X i , Y i , and Z i are computed as follows: Where  i is the angle between the gantry at the i th projection and zero degree. d m is the distance of movement along the gantry projection. In the static mode, d m is always zero because no movement exists. In helical mode, d m is computed as follows: Using the pitch and field width from the patient archive, the distance of movement per projection can then be calculated.
Both OAR x and OAR y values were normalized at the maximum value and interpolation and extrapolation procedures were applied. The software can derive the OAR x value for a given off-axis distance along the lateral direction. OAR x values were extrapolated to ±30.0 cm. The OAR y values were extrapolated to±50.0 cm in the longitudinal direction. To identify the OAR y value, both the off-axis longitudinal distance and depth of the point of measurement are required. Using Equation (1), the TPR values can then be calculated. For each segmentation pattern, the S cp,j value is the output factor for the j th segmentation pattern and t ij is the delivery time in the j th segmentation pattern within the i th projection.
From Equation (3), it can be seen that both TPR and OAR y values depend on the depth between the surface and point of measurement. In our software, the second check calculation is based on the patient plan as calculated on a cylindrical water equivalent phantom (Standard Imaging, Middleton WI). The depth is obtained automatically from the plan and is unique for each patient and for the point that was selected for the second check calculation. (Here Figure  6. The geometrical setup using a cylindrical phantom).
The cylindrical phantom, also commonly known as the "cheese" phantom, is a quality assurance device that is included with the Hi-Art™ unit QA package. The radius of the phantom is 15.0 cm. Figure  6 shows an illustration of the geometry used in the software to calculate the depth for any given projection. In Figure 6, the line BP is the projection of the depth in the IEC-X and IEC-Z plane. The green lasers represent the isocenter and gantry rotation axis of the Hi-Art™ unit, and the red lasers represent the localization of point P. Point O is the geometrical center of the phantom. Point S represents the source position. Angles  1 and  2 are calculated as: With these parameters known, the length of BP can be derived from with the following relationship: The final depth from the surface to the point of measurement, P, is then deduced from BP.

OAR y optimization
OAR y is a parameter dependent on three factors: the off-axis distance in the longitudinal direction, the depth of the point of calculation in the patient or phantom and the segmented field size. Interpolation and extrapolation procedures were used to derive the value of OAR y for randomly selected points of measurement. This can be a very time consuming procedure when performed on the fly.
To optimize this process and to gain substantial time savings (tenfold) in the calculation, we implemented an optimization method for where  is the Gaussian variation factor between the two segmentation patterns, LO 1 and LO 2 are the leaf openings, and y is the off-axis distance. In the MU-Tomo software, f 2 denotes the profile from the fully open field size and f 1 denotes the specific segmentation field. Fluctuations of the tail region were separated into two subregions: one between the 10% to 5% OAR y , and the other from 5% to the end of the OAR y profile. In the tail region, the two sub-regions were described with equation (10) of two different lambda values, denoted as exponential factors  and  respectively. Optimization was executed by varying  and  at 50 different channels from 0 to 49, separately, and on 50 different patients. Thus, a total of 50×50×50 dose differences were calculated to complete the optimization. For the 50×50 exponential factor channel pairs, the absolute values of dose differences from all patients were summed together to provide a volume map. For example, if d1%, d2%, …., d50% represent the dose differences from the 50 patients with one specific  - pair, the volume for this specific channel is equal to the sum of |d1|, |d2|,… and |d50|. The minimum volume is used in order to obtain the optimized parameters.

MU-Tomo benchmarks
Two methods were used to benchmark the software. First, treatment plans with static beam deliveries were calculated with the MU-Tomo software, and the dose was measured using the Hi-Art™ unit. For these plans, the measured and calculated doses were compared for benchmarking purposes. For the second part of testing, fifty previously treated TomoTherapy ® treatment plans were randomly selected for a point dose comparison using the planned versus MU-Tomo calculated point dose. Differences between the TomoTherapy ® TPS and MU-Tomo calculated doses were evaluated.

Software validation
The second check software, MU-Tomo, has been validated on the following four different phantom studies: (1) Dose calculation using a phantom composed of a 20.0 cm thick, rectangular slab of Virtual Water™, (2) Step valley dose calculation in the same Virtual Water™ phantom, (3) Dose calculation of a monthly output check procedure in the cheese phantom, and (4) Dose calculation of seven patient dose validation (DQA) treatment plans in the cheese phantom.
In the first validation phantom study, the phantom was composed of multiple rectangular slabs of Virtual Water™ and arranged such that the point of ion chamber measurement was centered in the phantom  The geometrical setup using a cylindrical phantom-also referred to as the cheese phantom-for depth calculation to the point of measurement, P, as a the cheese phantom-for depth calculation to the point of measurement, P, as a function of projection angle. S represents of the source location, A represents function of projection angle. S represents of the source location, A represents the isocenter of the HiArt™ unit, and B represents the surface point on the the isocenter of the HiArt™ unit, and B represents the surface point on the phantom. phantom.  The third phantom study validated the helical radiation delivery on the cheese phantom. The treatment plan delivered is used routinely as IMRT dose output verification and is part of our monthly QA checks. The plan is a helical IMRT delivery with treatment parameters of: FW= 2.5 cm, pitch= 0.287 and time= 262.4 seconds. The ion chamber point of measurement was 0.5 cm beneath the geometric center of the cheese phantom. The dose difference between the MU-Tomo calculated point dose and ion chamber measured point dose was -0.31%.
The fourth phantom study was for validation of patient IMRT DQA plans on the cheese phantom. Six IMRT DQA treatment plans were created. The DQA treatment plans were created using the optimized patient-dependent treatment fluence and mapped onto the cheese phantom. These same treatment plans were imported into MU-Tomo and a point dose was calculated. Figure 7 (B) shows the result for this validation. For all six plans, dose differences between DQA planned point dose and MU-Tomo calculated point dose were within 2.0% with a mean of 0.08% and standard deviation of 1.3%.

Algorithm optimization and second check on fifty patient plans
The OAR y optimization was quantified with a Gaussian variation factor. Based on the commissioning data of the OAR y optimization, we can determine the best values for the two exponential factors, which directly decide the variation amplitudes at the profile's tail region: one value models the 10% to 5% region of the relative profile and the other value models the region from the 5% value to the end of the profile. (Here Figure 8. The dose difference volume map).
Each channel at the  - plane in Figure 8 represents an iteration of the absolute dose difference summed over all 50 patients' second check results. Volumes were generated from absolute values of percentage dose differences and multiplied by 100. The minimum volume was used to obtain the optimized parameters by relating the selected channel to the corresponding exponential factors. Treatment plans were sorted by target location. For the optimization study, 50 patients were selected from a broad range of sites: 15 prostate, 8 lung, 14 head & neck and chest, 7 abdomen/pelvis, 2 liver, and 4 brain cancers. They were all planned using a 2.5 cm field width. (Here Figure 9. Second check result summary). Figure 9 shows percentage calculated point dose differences between TomoTherapy ® TPS and MU-Tomo for the 50 treatment plans when evaluated before and after the algorithm optimization. Before the optimization, the percentage dose difference was up to ±7.0% with a mean value of -0.99% and standard deviation 2.9%. After optimization, the results improved and all 50 cases were within a 5.0% dose difference-49 of them within 3.3%. The mean dose difference was 0.22% with a standard deviation of 1.77%.

Discussion
MU-Tomo utilizes a correction-based analytical dose calculation method requiring dosimetric functions, archived QA plans, and the point dose measurement. The aim of the software development was to perform an independent dose calculation that can be completed within one minute. Most of the second check procedures performed for this study finished within twenty seconds using a computer with an Intel dual-core CPU@3 GHz processor and 3 GB RAM. Because of this rapid calculation time, MU-Tomo can serve as a quick and accurate secondary check method as compared to a full Monte Carlobased calculation.
Dosimetric functions embedded in the software are obtained from the Hi-Art™ commissioned data. It is important to note that these data were obtained during the commissioning of our Hi-Art™ unit, and the software may need to be re-commissioned for use with other Hi-Art™ machines. The software assumes a source to axis   distance of 85 cm for the tomotherapy unit. During the software development, the off-axis ratio profiles were averaged around the centroid symmetrically, so that it is more reliable for both positive and negative coordinates.
The point dose calculation that we perform is subject to any errors that are inherent to the dosimetric functions that we extracted from the planning system. The OAR and PDD were measured per manufacturer's recommendation and were checked by tomotherapy at the time of commissioning. However these data are subject to approximately 1% measurement error. In addition, every time we change the target or any other beam control components in our unit, we do have about 1% discrepancy between our commissioned and current clinical data. Those discrepancies are not significant but could introduce as much as 3% uncertainty in the dose prediction from our software. If the user updates the software parameters with the corresponding dosimetric functions every time there is a change in the tomotherapy unit, the agreement should improve. We recommend that each clinic using the software imports their site specific data and performs an acceptance of this software.
Optimization of the OAR y using a Gaussian variation method benefits MU-Tomo in two aspects: (1) the calculation accuracy improves and (2) only the OAR y from the fully open field size is required, improving the calculation efficiency by a magnitude of ten. The reason for defining the tail region from 10% of the maximum value and below is because most fluctuations occurred in that region. The reasons to divide the tail region into two sub-portions in this study are that fluctuations within those two portions of the tail show different amplitudes; the first portion from 10% to 5% affects the dose calculation more than the second portion from 5% to the end. We did not pursue further segmentation of the profile tail region, as that would increase computation time with no significant improvement in the calculation accuracy. The optimization applied in this study may not be the best overall method, but it yielded good results with fast calculation times.
An additional feature of this second check software is the usage of archived patient documents. Initial gantry starting angle, initial position of IEC-Y coordinate of the target, and the sinogram are read in, enabling the software to determine the initial radiation projection position. Plan delivery duration, pitch of each specific plan, field width and the sinogram enable the software to accumulate the dose to the calculation point. To perform a secondary point dose check, a cylindrical phantom was embedded into the software. The symmetrical shape of the phantom simplifies the calculation of depth for any given projection angle. Because of the phantom geometry, the dose calculation point for the independent dose calculation can be placed anywhere inside the cylindrical phantom and is not limited only to the geometrical center. If a different phantom were used, recommissioning of the depth calculation method would be necessary.
With regards to the study performed by Gibbons et al., the MU-Tomo software has several novel developments that increase calculation accuracy and efficiency. The averaging of OAR x and OAR y profiles over the positive and negative axis sides removes fluence fluctuations, the optimization of OAR y values improves the dose calculation accuracy and efficiency, the application of the Mayneord F factor corrects the percent depth dose values for varying SSD, and a cylindrical phantom geometry allows fast and automatic depth calculation to the user defined dose calculation point. Furthermore, the dose calculation point for an independent dose calculation using MU-Tomo can be located at any point inside the cylindrical phantom, which makes our point dose verification more general. Although in this work we present a point based calculation comparison, it is possible to extend our method to include several points in three dimensions. This is currently under investigation.

Conclusion
An independent dose calculation software, MU-Tomo, has been successfully developed and benchmarked. Fifty treatment plans from different treatment sites were evaluated with our software and dose differences between measurement and calculation were found to be less than 5.0% (49 plans were within 3.3%). For all the cases, the mean dose difference was 0.22% with a standard deviation of 1.77%. Results show that the MU-Tomo software is able to perform independent dose calculations accurately and quickly and may be used to satisfy the clinical need for a secondary dose calculation of TomoTherapy ® treatment plans.

Disclosure
This project was supported in part from a research grant by Oncology Data Systems, Inc., Oklahoma City, OK, USA