We propose the following expression to model the 3D distribution of PPD (in mGy/MU):

where the coordinates x, y, and z  (cm) indicate the position of each calculation point in a coordinate system with the origin at the treatment isocenter (whichever it is). x, y, and z go in the anterior–posterior, left–right, and caudal–cranial directions, respectively (see Figure 1 ). The model is intended to be used for dose estimation outside the 5% isodose surface from where TPSs are not accurate enough (3, 7). In agreement with other published works (14, 21, 25), our model assumes that for distances to the isocenter larger than 40 cm, the main contribution to PPD comes from the leakage (L _u ), which is considered constant for the purposes of this work.

Equation 1 has some similarities with the model previously proposed by Sánchez-Nieto et al. (13). As in that work, the following correction factors considered here are

● F(f): It corrects the field size when it is different from the one used in the reference treatment (see Reference treatment on the ICRP 110 male phantom) .  This correction is essential as the scattered radiation is field size dependent (25–27). In this work F(f)=F _U (f)/F _R , where F_U and F_R are the areas representing the field sizes used in the user and reference treatment plan, respectively. For field size calculation, we propose to take the average of the areas inside the 50% isodoses at the coronal and sagittal planes of the 3D dose distribution at the isocenter level. The estimated value for F_R was 149.2 cm².

● ϵ(MU): It corrects the number of monitor units (MUs) when they differ from the reference treatment plan (see Reference treatment on the ICRP 110 male phantom). This correction accounts that the PPD depends on the number of MUs corresponding to each treatment. In this work, ϵ ( M U ) = E U ( M U ) E R , where E_U represents the efficiency of the user treatment (in terms of the prescribed dose at the isocenter per MU) and, E_R is the treatment efficiency of the reference treatment plan. E_R was calculated, for the calibration conditions of the linac for which the reference treatment was created (1cGy/MU at Source to Surface Distance (SSD), at d_max ), as the MU that delivers 2 Gy to the isocenter of the ICRP 110 phantom as if it was made of water ( E R = 2 · 10 3 G y 278   M U = 7.2 m G y M U )

● L_u : It corrects the leakage value whenever it is different from the one used in the reference treatment (L _r ). This quantity should be measured (in mGy/MU) for every accelerator, but if this parameter is not available, we recommend using the value in this work as an approximation (see Results).

The values of the fitted coefficients  A ₁(mGy cm ² UM ⁻¹) ,  A ₂(mGy cm UM ⁻¹), and  A ₃(cm ⁻¹)  were obtained by fitting the model to the 3D PPD distribution simulated with MC for the reference treatment plan (see Reference treatment on the ICRP 110 male phantom for more details).

In summary, to use this model, the user requires for each calculation point (coordinates in cm ), E _U (MU) in m G y M U , F_U (f) in cm² , and L_u (when available) in mGy/MU. If the absolute absorbed dose is needed, the total MU will be additionally required.

The reference treatment was an equally spaced eight-field isocentric plan centered at the mid-abdomen of the adult reference computational phantom ICRP 110 (28), with 10×10cm² open fields. The whole-body dose distribution was generated by an MC simulation (BEAMnrc code) of an Elekta Axesse with the Agility collimation system, up to 40 cm from the isocenter. The technical details of the MC simulation can be found in Sánchez-Nieto et al. (7). The MC simulation of the ICRP 110 considers the electronic density of each voxel. The uncertainty of the MC dose to points is given directly by the BEAMnrc code within the dose output file (“*.3ddose”) as a relative error value array in row 6 of the file.

Parameters A ₁, A ₂, and A ₃ were obtained by fitting Equation 1 to the 3D MC dose distribution corresponding to the reference treatment (ϵ =1 and F =1) using the fminsearch function in MATLAB^® (version R2021a). No information about the electronic density is considered by Equation 1 but the spatial position of the voxels. Only phantom voxels outside the 5% isodose surface were considered for the fitting, as the TPSs accurately estimate the dose distribution inside (7). Voxels representing the body contour were also excluded for the parameterization due to possible electron contamination, which is not considered by this model. As the geometry of the MC simulation did not include the gantry’s shielding, MC data were only used up to 40 cm from the isocenter, and, farther than this point, our measurement of leakage was used instead ( L r = L u = 0.001 m G y M U ). The fitting process gave the values of the constant coefficients. A ₁(mGy cm ² UM ⁻¹), A ₂(mGy cm UM ⁻¹), and A ₃(cm ⁻¹).

The model in Equation 1 depends on the three Cartesian coordinates; therefore, when the calculation is made on a whole-body CT, the model generates a 3D out-of-field dose cube from which the dose–volume histogram (DVH) of the contoured organs can also be extracted.

Representative isodoses of the reference treatment plan on the ICRP 110 phantom are shown in Figure 2 .

The fitted constant coefficients are (Model calibration)   A 1 = 37.890   ±   1.415   ( m G y   c m 2 M U ) , A 2 = 0.679 ± 0.074   ( m G y   c m M U ) , A ₃ = 0.007 ± 0.004 (cm ⁻¹).

The final (calibrated) version of the model can therefore be written (for points outside the 5% isodose) as in Eq (2).

The uncertainty of the model was calculated considering the absolute percentage differences between the doses given by the model and the ones given by MC, relative to MC. In total, 95% of all points (x, y, and z) presented an absolute percentage difference<23.2% (the average percentual difference was 7.84%).

Hereafter, the model in Equation 2 will be named Periphocal 3D. Figure 3 depicts the peripheral dose to points relative to the isocentric dose, estimated by Periphocal 3D and MC used for calibration. Note that this model does not use any electronic density information.

Although Periphocal 3D represents an improvement compared to its previous version, there are some limitations to take into account. Regarding geometry and X-ray attenuation, it is worth noting that Periphocal 3D was calibrated with a nearly symmetrical eight-field treatment (see Figure 2 ) and Equation 1 has spherical symmetry. Thus, the model should become less accurate when non-symmetrical isodoses, usually associated with plans with fewer beams, as in Figure 5 (lung case), are generated. Luckily, VMAT and IMRT treatments usually consider several beam entries around the isocenter. Another possible issue is the attenuation of X-rays inside the patients when significant tissue heterogeneities are present. Since Periphocal 3D was calibrated in the abdomen, which is mainly composed of soft tissue, its application to the treatments of tumors in tissues with different X-ray attenuation characteristics will necessarily imply worse accuracy. This can be observed in Figure 6 , where calculations within the soft tissue after the lung shows an overestimation of the dose because it does not take into consideration the smaller backscatter contribution from the lung tissue. Additionally, in bone, as Periphocal 3D assumes that everything is water, the calculated absorbed dose underestimates the MC calculation.

It is worth noting that the large discrepancy between MC and Periphocal 3D calculations in the lung plan is mainly associated with the use of small fields (5 × 5 cm²) for which a sharp drop at 21 or 22 cm approximately from the isocenter occurs in correspondence with the edges of the primary collimator that provides additional shielding. The same phenomenon has been described elsewhere (30, 31). This drop is much softer for larger beams, as seen in Figure 3 for the 10 × 10 cm², in agreement with Kaderka et al. and Jagetic and Newhauser (30, 31). That is, the full MC simulation recreates a profile shape that the Periphocal 3D does not (see Comparison with another analytical model to verify how a more complex physics-based model also fails in fitting the decrease due to the additional shielding of the primary collimator).

The A ₁, A ₂, and A ₃ coefficients were established by fitting the model to the dose distribution generated by the MC simulation of an Elekta Linac. Stoval et al. (18) showed in Appendix A how peripheral dose profiles depend on the design and construction of the machine head and collimators. However, our measurements—figure 1.a in Sánchez-Nieto et al. (13)—of peripheral doses for the same treatment delivered on the same phantom using different combinations of energy and linacs allowed us to conclude that the observed variability was within the model´s uncertainty. Those who can generate a 3D dose distribution of their specific linac from a full MC model or thorough experimental measurements can obtain their specific A ₁, A ₂, and A ₃ coefficients.

Periphocal 3D does not work for skin dose calculations. A peripheral skin dose can be separately estimated based on other previously published works (32).

The availability of 3D dose distribution allows for the calculations of DVHs for peripheral organs (input for some models of secondary cancer risk), dose profiles along any axis or 2D dose distributions on any plane. The 3D dose distribution may also be helpful, for example, when the patient has a pacemaker, an intern defibrillator at positions where TPSs are not accurate, or for any of the harmful effects listed by Mazonakis and Damilakis (3).

Figure 7 shows the dose to organs for a lung treatment in comparison to PERIPHOCAL and MC. This is a case for which Periphocal 3D presents some limitations, as discussed in the previous section. Even so, the model offers an improvement compared to PERIPHOCAL for all organs but the thyroid. Both PERIPHOCAL and Periphocal 3D overestimate the dose, which can be explained by the highly non-symmetrical geometry of the dose distribution in this area. For the prostate, despite Periphocal 3D performing much better than PERIPHOCAL, a significant difference compared with MC is still present. As we already mentioned, the geometry of our MC did not include the gantry’s shielding and, therefore, might not correctly account for leakage. Thus, the dose for the prostate (farther away than 40 cm from the isocenter) given by MC may be underestimated. A detailed study of the effect of linac’s shielding on the leakage is being conducted. For the urinary bladder, close to the prostate but closer to the isocenter than 40 cm, Periphocal 3D and MC agree.

This model can be used retrospectively and prospectively (for example, using the virtual whole-body CT generated by our home-made software (24) to calculate in a systematic way dose to peripheral organs and, together with clinical follow-ups, detecting possible secondary cancers, creating a database for a more accurate parameterization of secondary cancer models (3–6).

Sánchez et al. (13) showed that PPD does not significantly change with energy (differences within the model´s uncertainty), and thus, Periphocal 3D can be used even above 10 MV. However, neutron contamination might become relevant above this threshold. For those dealing with energies >10 MV, the total peripheral dose to organs should include the peripheral neutron dose to organs. The model published by Irazola et al. (33) can be used to estimate the neutron dose to out-of-field organs.

Schneider et al. (22) developed a nice physics-based analytical model of the total absorbed dose for the primary, scattered, and leakage radiation of square fields of 6 MV at any arbitrary point in a phantom. That work is one of the latest models published but tested only for square fields. They mention the validation of the model for arbitrary MLC aperture to determine the model applicability to IMRT treatments as a future step. Thus, to our knowledge, no other analytical model has been developed and tested for intensity-modulated treatments with which we can compare.

Figure 8 depicts the comparison between the dose as calculated by our model and Schneider´s model for irradiation with a 10×10cm ² field using a 6 MV beam. The parameter field size F_U (f) at the isocenter was estimated from the width of the profile at 50% of the isocenter dose as F _U (10)=103.63 cm ² . It was also assumed that ϵ _U = ϵ _R

Both models agree within our model´s uncertainties. Curiously enough, both models fail to reproduce the additional shielding of the primary collimator as predicted by Schneider´s MC model. Note that this effect is sharper for smaller fields, as in the case of Figure 6 , corresponding to a treatment plan of 5×5 cm².

In this work, the ICRP 110 phantom and the CT of the ATOM phantom were used for the 3D calculations. However, for personalized clinical applications, a whole-body CT of each patient is needed. We have solved this problem without actually irradiating the patient (which would be unacceptable due to the unjustified additional dose), following a methodology presented in the companion article. This methodology uses the always-available planning CT to generate an approximate patient-specific whole-body CT based on a rigid 3D image registration algorithm. The input for calculating the 3D dose distribution is the whole-body CT ignoring the differences in electronic density between the voxels but with all organs’ contours considered in the ICRP110.

A GUI was created in MATLAB^® (version R2021a) to ease the use of Periphocal 3D ( Figure 9 ). The output of our whole-body CT software is a MATLAB array. Periphocal 3D’s GUI takes less than 10 s to calculate the whole-body dose distribution after loading those MATLAB arrays. However, Periphocal 3D’s GUI also accepts CT DICOM files. Additionally, for all organs segmented in the input CT, a DVH is created. Figure 8 depicts an example of the mentioned GUI.