A comparison of anatomical and dosimetric variations in the first 15 fractions, and between fractions 16 and 25, of intensity-modulated radiotherapy for nasopharyngeal carcinoma.

The purpose of this study was to compare anatomical and dosimetric variations in first 15 fractions, and between fractions 16 and 25, during intensity-modulated radiotherapy (IMRT) for nasopharyngeal carcinoma (NPC). Twenty-three NPC patients who received IMRT in 33 fractions were enrolled. Each patient had two repeat computed tomography (CT) scans before the 16th and 25th fraction. Hybrid IMRT plans were generated to evaluate the dosimetric changes. There was a significant decrease of the transverse diameter of nasopharyngeal and neck as well as gross tumor volume (GTV) in the primary nasopharyngeal carcinoma (GTVnx) and involved lymph nodes (GTVnd) during the first 15 fractions, and between fraction 16 and 25 (p < 0.05). Consequently, there was a significant reduction of the percentage of the volume receiving the prescribed dose (V100) of CTV1 and GTVnd, which was more prominent after the first 15 fractions treatment compared to that between fraction 16 and 25 (p < 0.05). Additionally, there was a significant increase in the mean dose (Dmean) and percentage of volume receiving ≥ 30 Gy (V30) to the bilateral parotid in the first 15 fractions (p < 0.05), but not between fraction 16 and 25. While the maximum dose to the spinal cord was significantly increased both in the first 15 fractions, and between fraction 16 and 25 (p < 0.05), the increase of the percent of spinal cord volume receiving ≥ 40 Gy (V40) was significantly higher in the first 15 fractions compared to that between fraction 16 and 25 (p < 0.05). Based on the dose constraint criterion in the RTOG0225 protocol, a total 39.1% (9/23) of phantom plan 1 (generated by applying the beam configurations of the original IMRT treatment plan to the anatomy of the second CT scan) and 17.4% (4/23) of phantom 2 (generated by applying the beam configurations of the replan 1 to the anatomy of the third CT scan) were out of limit for the dose to the normal critical structures. In conclusion, our data indicated that anatomic changes resulted in more predominant dosimetric effects in the first 15 fractions, and between fractions 16 and 25, of IMRT.

I. INTRODUCTION

The fundamental physical quantity of interest for relating radiation treatment to its outcome is the absorbed dose. Furthermore, a significant component of the absorbed dose at a point is due to the scattering of the primary beam; therefore, it is essential to include the correct amount of scattering in any dose calculation algorithm used in treatment planning system (TPS). These algorithms are correction‐based or model‐ (or convolution–superposition‐) based. Correction‐based algorithms use parameters of dose measured in water phantom and correct the data to apply to the patient's specific situation. This requires the percentage depth dose (PDD) for a number of square fields, a set of profiles for a number of square fields measured at some standard depths, and phantom scatter factor curve at reference depth as a function of the field size of the square fields.( , Model‐based algorithms directly compute the dose to the patient by modeling the beam and interactions in the patient, which require some measurements to set parameters and verify the model. Monte Carlo methods, implemented to mimic the basic processes in a straightforward way, have served many purposes in medical physics. However, they have not yet become suitable for routine treatment planning of photon beams due to their huge requirement for CPU time. Consequently, analytical methods which would be reliable within an acceptable limit of error and for reasonable range of parameters to calculate PDD and profile at any depth are increasingly desirable. Increased requirements on standards for safety and quality assurance during treatment have also emphasized the important role of simple dose calculation methods for independent checks of the output from treatment planning systems.

For a given energy spectrum incident on a homogeneous medium and assuming lateral electronic equilibrium, the primary component of central axis depth dose for any field shape will be the same, and only differences in the scatter component will affect the final shape of the central axis depth dose; therefore, the concept of equivalent field size based on the separation of the primary and scatter radiation was proposed by Day and Aird. For regular fields, tabulated data was presented by Day and Aird, or some empirical formulas by fitting were carried out, such as the equal area‐to‐perimeter ratio (A/P).

Various methods have been used for the prediction of off‐axis ratio for a symmetric open field. Fermi‐Dirac distribution function suggested by Kornelson and Young and Wood‐Saxon term applied by Pal et al. represent the off‐axis ratio (OAR) in the SAD and SSD techniques, respectively. Usually these methods need data fitting at several depths.

In the 1990s, some studies have addressed asymmetric and wedged asymmetric fields using symmetric field data. , , , The method based on work by Thomas and Thomas generates asymmetric field profiles by computing the off‐center ratio (OCR) of the asymmetric field while using output factors and PDDs of the equivalent symmetric field. A second method, proposed by Kwa et al., applies to the situation where only one of the independent jaws is closed down to form an asymmetric field of smaller width or length than the original symmetric field. This method uses the original symmetric field profile corrected point‐by‐point by a correction factor.

Accordingly, we put forward a correction‐based dose calculation algorithm based on equivalent field concept for the fixed source‐to‐surface distance (SSD) formalism and develop proper analytical expression which could equip the computer with this formula for direct calculation of the profiles at any depth for some fields (e.g., symmetric or asymmetric field with or without wedge and blocks).

This paper introduces an analytical method to calculate the equivalent field in regard to the central axis, first. It will then follow by generalizing this method to a more general case in which the equivalent field is calculated concerning off‐axis points. Then, some profiles will be measured to set the imperial correction factors required to the calculations. Finally, the calculated and measured profiles will be compared for some practical fields.

A. Theoretical background

A.1 Central axis percentage depth dose (PDD)

In this study, to calculate PDD for a custom field, the equivalent field method is used. A set of PDD is tabulated for a number of square fields measured along the central beam axis for both energies 6 and 18 MV and for open fields and every wedge angle separately. For any other square field size calculated, one can interpolate from the table of PDDs or utilize semiempirical equations for PDD depending on field size and depth.( ,

For any fixed point at fixed depth on the central axis in the medium, the primary component of the dose will be the same for all fields. It therefore follows that equivalency between standard and nonstandard fields is determined by the equal dose contribution along the central axis from scattered photons for the two fields. The method used to determine the scatter contribution the equivalent field in this study is further explained.

Consider a reference plane normal to the central axis, placed at the fixed source‐to‐surface distance (e.g., ). The field bound generated by the collimators can be projected over this plane as X1, X2, Y1, and Y2. The origin (e.g., ) is on the central axis. Suppose that there is a parallel beam striking the surface of the reference plane. Hence, each surface element on that plane acts as a source of scattered radiation. The amount of scatter radiation reaching the central axis is inversely proportional to the square of the distance between surface elements and the origin; therefore, an asymmetric irregular field (e.g., an asymmetric field with some shielded parts by cerrobend blocks or multileaf collimators (MLC)) to circular field equivalence at the central axis is calculated using the following equation: where , which shows the transmission of the blocks or the leaves of MLC in which is the block attenuation coefficient and is its thickness. shows the fraction of primary photons absorbed by the block when they are passing through it. Therefore, in the absence of the block, will be 1 and the second term on the left side will be removed. It should be noted the attenuation of the water phantom using the factor of was considered on both sides of Eq. (1), where was the water phantom attenuation coefficient for scattered photons. However, because of its negligible effect, the attenuation of the water has been eliminated.

Let us use a polar system and eliminate the small area at the origin of coordinate with radius e (refers to the ion chamber radius) to overcome singularity; therefore, by dividing the field into four parts regarding the origin and integrating both sides of Eq. (1) with some mathematic operations, it can be rewritten as: where the exponential parts are computable by numerical calculation using MATLAB software (The MathWorks, Natick, MA). As expected, e will be eliminated because it is present on both sides of Eq. (2). As a special case, the circular field equivalent to a square field is obtained by substituting S (side of square) in Eq. (2): where the integral is equal to using numerical integration. After calculating the size of the field side, the percentage depth dose at any depth can be interpolated from data tabulated for the square fields.

A.2 Computation of profile

Dose distributions along the beam central axis give only part of the information required for an accurate dose description inside the patient. Dose distributions in 2D and 3D are determined by central axis data in conjunction with off‐axis dose profiles. Megavoltage X‐ray beam profiles consist of three distinct regions: central, penumbra, and umbra. The central area represents the central portion of the profile extending from the beam's central axis to within 1–1.5 cm from the geometric field edges of the beam (e g., the 50% dose level points on the beam profile). In the penumbral region of the dose profile, the dose changes rapidly and depends also on the field defining collimators, the finite size of the focal spot (source size), and the lateral electronic disequilibrium. Umbra is the area outside the radiation field, far removed from the field edges. The dose in this district is generally low and results from radiation transmitted through the collimator and head shielding. The flattening filter also produces a differential hardening across the transverse direction of the beam which results in off‐axis peaking at a depth shallower than 10 cm. However, this model does not include the variations in off‐axis beam quality affected by a flattening filter.

Consider a symmetric open field, regardless of the differential hardening effect of the flattening filter in the central region, the farther from the central axis of the radiation field, the less amount of phantom scatter; therefore, to calculate the amount of the phantom scatter reaching the point , the equivalent square regarding is used (Fig. 1). Thus the total amount of the scattered radiation from the irradiated field reached to is considered to be equal to the scattered radiation from an equivalent square field to the central axis and the side of this square is calculated. In this way, by dividing the field to four parts as regards the point , Eq. (2) can be modified as follows: where illustrates the calculated equivalent field regarding the point . It can also be shown that by substituting in Eq. (4), it will be transformed to Eq. (2). Now, this method is evaluated for the crossline profile at depth 10 cm for 6 MV energy (Fig. 2). At first, is calculated for some points on the crossline using Eq. (4). Subsequently, the percentage depth dose at depth 10 cm for 6 MV energy on the central axis can be interpolated for from data tabulated for the square fields (see Fig. 2, , dash line). As can be seen, it is in good agreement with the results of measurements up to 1.5 cm from the beam edge; therefore, for points near the edge of the beam, should be corrected.

Figure 1

An asymmetric field with shielded parts by cerrobend blocks or MLC is equivalent to the circular field with radius , both projected on the phantom surface. ε is the ion chamber radius at . The points (0, 0) or illustrate the central axis of the beam and any point in the irradiated field, respectively.

Figure 2

Profiles measured (solid) and calculated (circle) at depth 10 cm and 6 MV energy. (dash) shows PDDs on the central axis extracted from the data stored in the PDD table with regard to the equivalent square calculated from Eq. (5) for the point .

As previously mentioned, in the penumbra region, the dose changes rapidly around the geometric beam edge and creates a sigmoidal shape. Furthermore, the photon source has a Gaussian distribution (normal distribution). Accordingly, the 3D photon influences the distribution of the primary source in air, which can be calculated by analytical integration of the Gaussian functions leading to an error function. Therefore, to calculate the profile (Fig. 2 circle) at depth, d, should be multiplied by the correction factor ofjaws as follows: where is the correction factor of the X1 jaw at which consists of three regions: infield, near outfield (penumbra), and far outfield (umbra). In the near outfield, both the transmission of the jaw and the correction of the electronic disequilibrium play a role. However, in the far outfield, the transmission of the jaw just presents (see Eq. (6)). where is the transmission through the collimator jaws, is an empirical correction factor for electronic disequilibrium, and . are the variance of inside and outside of the field, respectively, which are determined empirically for any energy. Due to beam divergence, in order to draw in‐depth (d) profile, magnification ( should be considered. The for X2, Y1, and Y2 are defined as Eq. (6) similarly. Obviously, the empirical constant mentioned above is independent from the depth, and field size due to includes the effects of both of them.

In the wedged field case, due to beam hardening, should be extracted from the table of measured PDD for several wedged square fields. Also, the primary beam is attenuated due to the wedge thickness variation. Hence, the wedge correction factor for any point off‐axis is described as the ratio of primary beam attenuation at that point to the central axis, as follows: where is the wedge attenuation coefficient, and are the wedge thickness corresponding to the points and (0,0), respectively. Finally, by multiplying Eq. (5) in Eq. (7), the calculated wedge profile is obtained.

It should be noted that equals to one for the open field.

For irregular fields, the blocks act like the collimators with variable thickness. So Eq. (6) can be rewritten for the block correction factor: where represents the coordinate of the point on the edge of block projected on the phantom surface and the crossline profile passes through it. shows the block transmission where is the block attenuation coefficient and is its thickness. Eq. (9) can be used for MLC, where represents the coordinate of the leaves edge of MLC projected on the phantom surface and shows the leaves transmission.

II. MATERIALS AND METHODS

Varian 2100 C/D and Siemens Primus Plus with 6 and 18 MV photon beam were used for measurements. The treatment units were equipped with independent jaws assigned to as Y1 and Y2 for the upper and X1 and X2 for the lower jaws. A Scanditronix blue phantom (Wellhofer, Germany) with two 0.13 CC ionization chambers (IBA, Germany) were employed for the measurement. Omni‐AcceptPro 6.5 software (Wellhofer, Germany) was connected to the interface and utilized for collecting and recording data from two chambers.

At first, PDDs were measured for a large number of square fields for both energies and for 45° wedge which were employed for profiles in any depth in this model. Then, to access the accuracy of the model for profile prediction, several profiles were measured for some special treatment fields such as symmetric, asymmetric, wedged asymmetric, and irregular fields (shown in 3, 4, 5, 6, 7).

Figure 3

Measured and calculated crossline profiles and γ‐index (a) in depths 5, 10, and 15 cm for a symmetric open field 10 by 10 cm, Varian linac, and 6 MV energy; (b) measured and calculated crossline profiles and γ‐index in depths 5, 10, and 15 cm for a symmetric open field 10 by 10 cm, Varian linac, and 18 MV energy.

Figure 4

Measured and calculated crossline profiles and γ‐index in depth 10 cm for a symmetric wedged field 10 by 10 cm, wedge angle of 45°, Varian linac, and two energies 6 and 18 MV. The wedge attenuation coefficient was 0.5 and 0. 45 1/cm for energy of 6 and 18 MV, respectively.

Figure 5

Measured and calculated the asymmetric profiles and γ‐index in depth 10 cm for asymmetric open field 10 by 10 cm with 3 cm offset at 100 cm SSD, Siemens linac, and 18 MV energy.

Figure 6

Measured and calculated asymmetric wedged profiles and γ‐index in depth 10 cm for field 10 by 10 cm with 3 cm offset at 100 cm SSD, wedge angle of 45°, Siemens linac, and two energies 6 and 18 MV. The wedge attenuation coefficient was 0.3 and 0.24 1/cm for energy of 6 and 18 MV, respectively.

Figure 7

Measured and calculated symmetric shielded profiles and γ‐index (a) for three blocks dimensions of 7 by 7cm with 2.5, 4, and 7.9 cm thickness in the right corner of 10 by 20 cm field in depth 10 cm, Varian linac, and 6 MV energy, with the attenuation coefficient of block as 0.44 (1/cm); (b) measured and calculated symmetric shielded profiles and γ‐index for three blocks dimensions of 7 by 7 cm with 2.5, 4, and 7.9 cm thickness in the right corner of 10 by 20 cm field in depth 10 cm, Varian linac, and 18 MV energy, with the attenuation coefficient of block as 0.42 (1/cm).

Several dose distribution comparison methods have been developed based on various combinations of doses and spatial acceptance tolerances, including the simple dose difference (DD) test and the distance‐to‐agreement (DTA) test.( , The gamma index calculation and modified dose difference (MDdiff) evaluation are dose comparison methods which produce a quantitative measure based on both dose and spatial criteria. , , In this paper, the μ–index evaluation was utilized which was accorded a 3% dose criterion and a 3 mm dose‐to‐agreement (DTA) acceptance criterion.

Finally, a homemade computer program was developed in MATLAB 7.14 on Windows platform to process the data quickly, to plot the profiles, and to compare the calculations with measurements (see Appendix).

III. RESULTS

At first, to set the empirical correction factors, profiles were calculated for a series of estimated correction factors for a standard field 10 cm by 10 cm in depth 10 cm and both energies 6 and 18 MV. Then, the calculated profiles were compared with the measured profiles using γ‐index evaluation. This procedure was repeated for different correction factors until the best ones were selected (e.g., γ ≤ 1). The final correction factors were obtained as follows:

As shown in Figs. 3(a) and (b), these correction factors could be applied for depths 5 and 15 cm with acceptable γ‐index due to interpolating PDDs at 5 and 15 cm. In other words, the calculated correction factors for a given depth can be generalized to other depths. As can be seen in Figs. 3(a) and (b), the γ is less than 1 for most points. However, in the penumbra, γ‐index values represent larger numbers. The maximum γ values were 1.00, 0.93, and 1.00 for 6 MV. and 1.10, 0.75, and 0.73 for 18 MV at depth 5, 10, and 15 cm, respectively. It seems that the misalignment of the chamber axis plays role in the relatively high gamma values visible in the penumbral region. To clarify this, one can slightly shift the measured profile. Consequently, the γ values will be reduced in the penumbral region.

In the second place, the calculated correction factors were used to the wedged field. Figure 4 shows the measured and calculated profiles in depth 10 cm for a symmetric wedged field 10 by 10 cm, wedge angle of 45°, Varian linac, and two energies 6 and 18 MV. As expected, the calculated profiles show good agreement with measurements by the γ less than 1 for most points, like the symmetric open field. Using the PDDs table regarding the wedged square fields will obviate the beam hardening which can affect the shape of profile at different depths. The maximum γ value was 0.96 and 1.05 for 6 MV and 18 MV, respectively.

For asymmetric fields, Fig. 5 represents the profiles measurement and calculation in depth 10 cm for an asymmetric open field 10 by 10 cm, with 3 cm offset at 100 cm SSD, Siemens linac, and 18 MV energy. As can be seen in Fig. 5, the asymmetric field has larger γ values in the central region, unlike the symmetric field. The maximum γ value was 1.09 in the central area. As shown in Fig. 6, the measured and calculated profiles were plotted in depth 10 cm for an asymmetric wedged field 10 by 10 cm with 3 cm offset, wedge angle of 45°, Siemens linac, and two energies 6 and 18 MV. As was mentioned for asymmetric field, here there are some points at which γ‐indexes are greater than 1 in the center district. This issue seems to be due to the lack of the flattening filter effect on the primary beam in this model, and this effect was larger in the profiles of Siemens linac.

The profiles for a rectangular field () with three blocks () with 2.5, 4, and 7.9 cm thickness in the right corner were evaluated in depth 10 cm for two energies 6 and 18 MV (Figs. 7(a) and 7(b)). Note that, to avoid clutter, these figures do not show the γ calculations. For more details, the mean γ‐index values and 1 standard deviation are given in Table 1 for all fields.

Table 1

Mean γ‐index and standard deviations of calculated and measured dose profiles for the fields used in this study

(Fig. 3(a), 3(b)) 10×10−open
	γ‐index–6MV				γ‐index – 18 MV
Depth	Infield		Outfield		Infield		Outfield
(cm)	AVR	STD	AVR	STD	AVR	STD	AVR	STD
5	0.20	0.2	0.30	0.29	0.21	0.21	0.32	0.23
10	0.16	0.17	0.56	0.25	0.20	0.19	0.32	0.19
15	0.08	0.09	0.61	0.33	0.17	0.16	0.45	0.18
(Fig. 4) 10×10−wedge45°
10	0.36	0.23	0.51	0.26	0.25	0.20	0.38	0.22
(Fig. 5) 10×10−open(offset3cm)
10					0.25	0.22	0.47	0.31
(Fig. 6) 10×10−wedge45°(offset3cm)
10	0.51	0.24	0.44	0.27	0.45	0.27	0.37	0.31
(Fig. 7(a)),7(b)) 10×20shielded
10(small)	0.43	0.73	0.51	0.26	0.44	0.44	0.46	0.53
10(med)	0.41	0.64	0.50	0.29	0.64	0.77	0.36	0.16
10(large)	0.49	0.62	0.56	0.24	0.47	0.67	0.36	0.17

IV. DISCUSSION & CONCLUSION

Today, the use of radiotherapy treatment planning systems (TPSs) is inevitable. Beam profile and PDD are the parameters used to verify the dose calculation algorithms of TPS; therefore, a patient‐independent model calculating beam profile and PDD can be used to minimize the number of measurements for verification processes.

The algorithm used in this study calculates the beam profile in the water phantom by separation of the primary and the scatter beams, for irregularly blocked fields of open or wedged photon beam. The requirements are only the PDDs for several squares for each energy and wedge angle and some profiles to set the empirical correction factors. As can be seen in Table 1 and 3, 4, 5, 6, 7, the comparison with measurements using γ‐index shows that the accuracy of the calculated dose distributions fits well in a 3% error in low‐dose gradient region, except for asymmetric fields. In general, the primary dose rate at shallow depths in the phantom may actually increase at distances away from the central axis (called horns) as a result of flattening filter effect on the radiation beam. , , A flattening filter correction that depends on depth in a phantom and radial distance from the central axis is required to model the increase in dose rate away from the central beam axis that is not included in this paper. The comparison also shows an approximate agreement in a 3 mm isodose shift in the penumbra region. The dose near the edge of the beam changes rapidly and depends also on the field defining collimators, the source size, and the lateral electronic disequilibrium. Since the photon source has a Gaussian distribution (normal distribution), the dose falloff around the geometric beam edge is sigmoid in shape, for which an error function has been employed. However, considering the fact that the Gaussian distribution alone, like the models proposed by Pal et al. The Kornelsen and Batho model suffer from lack of electronic disequilibrium; therefore, a correction factor, , has been considered here.

The correction‐based algorithm( , and the previous methods mentioned( , need beam profile data for a large number of depths to predict the off‐axis ratio. In addition, a general problem with empirical scatter scaling techniques is that they are developed for open, not modulated beams, and there is a great need to improve these models to include effects from modulations. Nevertheless, the findings of the current study show the empirical correction factors are independent on depth. The depth independence is due to which includes the effects of the depth. In addition, the correction factor of jaw, , acts similarly not only for symmetric fields but also for asymmetric fields. The depth dose differences between symmetric and asymmetric field at any point are indicated by unlike the previous studies. , , , ,

This method is valuable due to the calculated profile in the presence of the blocks with variable thickness. The variable thickness shields can protect organ at risk inside the field. This means that organ at risk and normal tissue receive desirable tolerance dose, whereas the lethal dose is delivered to the tumor; further research should be done to investigate this in compensator‐based IMRT.

Consequently, a general algorithm was proposed to calculate the profile at any depth for symmetric and asymmetric, wedged or open photon fields. The advantages of this algorithm are depth independence, minimum measurement data requirements, quick‐run, low cost, easy‐to‐use, and fairly accurate. Moreover, this algorithm can be carried out to plot 2D or 3D isodose for multiple field treatment planning. A number of important caveats need to be considered about the calculation algorithm. First, it is limited to rectangular and triangular block shape. Second, it does not include the curvature of body surface and inhomogeneity. Finally, the flattening filter effect should be considered as another correction factor. It is recommended that further research be undertaken to solve these problems.

ACKNOWLEDGMENTS

The authors thank the reviewers for their comments that help improve the manuscript. The authors would also like to thank Dr. Rahmani for his unsparing support. This study was funded by the office of Vice‐Chancellor for Research of Joundishapour University of Medical Science Ahvaz, and the Radiotherapy and Oncology Center of Ahvaz Golestan Hospital.

APPENDICES

Appendix A: The Procedure to Transition Eq. (2) from Eq. (1)

Using a polar system and eliminate the small area at the origin of coordinate with radius e and dividing the field into eight parts regarding the origin and integrating both sides of Eq. (1), it can be written as:

Integration with respect to r leads to natural logarithm as follows:

In the next step, integration respect to θ and by more simplification, it will be rearranged as Eq. (2).

The algorithm to calculate profile was as follows:

PDDs were measured for a number of square fields along the central beam axis, for open fields and 45° wedge and for both 6 and 18 MV energies separately. The PDDs were tabulated for desired interval (e.g., 0.5 cm) and stored in the MATLAB program.

A set of required parameters such as field size, wedge angle, and blocks properties (e.g., thickness, location and size), was fed to the calculation algorithm.

For open symmetric field, the equivalent square, was calculated for any point placed on crossline using Eq. (4).

To plot the profile at depth 10 cm, was interpolated from stored data for both energies at depth 10 cm.

According 5, 6, profiles were calculated for a series of estimated correction factors for a standard field 10 cm by 10 cm in depth 10 cm and both energies 6 and 18 MV.

The γ‐index was calculated for every point; and were changed until .

Finally, and were calculated to plot profiles for asymmetric, wedged, and irregular fields at any depth.

Supplementary Material

Click here for additional data file.

Supplementary Material

Click here for additional data file.

Supplementary Material

Click here for additional data file.

Supplementary Material

Click here for additional data file.