A. Simulation procedure

In order to derive DgN(E) values, we used the photon quanta per mm² per mR parameter, K(E), for a given energy E (in keV) of incident photon beam in a medium of air as derived by Johns et al. [25],

where, (μ_en/ρ)_air, in cm²/g, is the mass energy absorption coefficient in air.

The DgN(E) coefficient, in mGy/mGy, is given by [26, 27]:

where, E_dep is the energy deposited in the breast tissue, fg is the glandular fraction of the breast, ρ is the breast density, R₁=8 cm is breast radius, R₂= 7.6 cm is breast tissue radius excluding the skin, T_skin is the skin thickness, T is the breast thickness and G(E) is a factor describing the proportion of the absorbed dose of glandular tissue to the overall breast tissues and given by:

where, (μ_en/ρ)_g and (μ_en/ρ)_a are the mass energy absorption coefficients for glandular and adipose tissue, respectively, for the given energy. The unit of the constant 1.8352×10^-17 is (mGy.g)/keV.

Using the recent mass energy attenuation coefficients of photons through air, provided by NIST laboratory [14], a Monte Carlo-based simulation program was tailored for calculating absorbed energy in the breast model for a given configuration of breast tissue composition, thickness and for monoenergetic photon beam. Then, an additional procedure was added to that program, in order to compute the G(E) factor during the simulation run-time, which in turn outputs the DgN(E) coefficients for each studied case.

The Geant4-based simulation program was tailored to model monochromatic X-ray emission from a "point-like" source and to follow particle transportation through the breast. The physical processes used were the photoelectric effect, the Compton and the Rayleigh scattering for photons and the bremsstrahlung, the ionization and the "multiple scattering" for electrons. The production thresholds were set to 1 keV for photons and 10 keV for electrons. From the existing suite of physical packages, ready to use within the simulation program, the "PhysListEmStandard" was selected. Elemental compositions of adipose, glandular and skin tissue were based on the work of Ref. [27].

For each set of (monoenergetic) X-ray energy, breast size and breast composition, a simulation run was executed using d6 entrant photons providing sufficient statistical precision. In order to simplify analytical form of the fitting functions to be used and due to the small contribution to MGD parameters, we omitted the interval of energy of 1–11 keV during the fitting stage of the database. Notice that such contribution does not exceed 1% of the MGD for the studied X-ray spectra. DgN(E) coefficients for 11–120 keV X-rays were computed. At least a total number of 1582 (113×7×2) runs were carried out for the data-base generation purpose covering 113 energy values (11–120 keV), and seven thicknesses (2–8 cm) for glandularities of 100% and 50%. The absorbed dose by the breast tissue was recorded for direct derivation of the DgN(E) coefficients (in the SI units of mGy/mGy).

As we followed the same simulation setup used by Boone et al. [11], the same source-to-detector distance of 60 cm was used. Changing the distance to 80 cm or 40 cm can deviate the beam energy at the vicinity of the breast by less than 1%, which is logical due to the fact that the addition or subtraction of a 20 cm layer of air medium can affect the primary beam trajectory by 0.4% and 0.06% for 30 keV and 100 keV photon energy, respectively. Moreover, the assumption of modeling the geometry of the X-ray tubes focal spot as a point-like rather than a disc is valid to about 1%. The replacement of a point-like source with a disc of Φ40 μm or Φ300 μm (as found in clinical realistic tests), deviate the DgN(E) with 0.213% or 0.444% for a given couple of photon energy, glandularity and thickness of (50 keV, 50%, 2 cm) or (10 keV, 100%, 8 cm), respectively. Whereas the focal spot size affects contrast of the images [28].

We can consider the homogeneous distribution of the glandular tissues within the breast as a limitation of the current model, as mentioned by Dance et al. [29] that the absorbed energy by the glandular tissue depends on the glandular tissue position in the breast.

The simulation was carried out using the Geant4 Monte Carlo code version 9.4.p01 under a Red Hat Enterprise Linux 5 workstation on an Intel Core i7 running with a 1 GB RAM and a 3.40 GHz CPU. The statistical uncertainty (2σ) associated with all the Monte Carlo calculations presented in this work is less than 1%.

B. Parameterization procedure

The following parameterization procedure was followed to reduce the large amount of DgN(E) values while retaining the prediction precision to an acceptable level of accuracy and using simple analytical form of the fitting equation. For accuracy of the proposed method, the r² statistical coefficient value was calculated. We fitted the overall data set of DgN(T, E) using Eq. (4) for the energy interval of 11–120 keV, for glandularities of 100% and 50% separately

And the MGD was calculated using Eq. (5):

where, ϕ(E) is the photon fluence per exposure in units of mR per photons per mm² and E_min and E_max are the minimum and the maximum energy spectrum limits.

Physically spoken, it is clear from Eq. (4) that any DgN coefficient explicitly depends on the coming photon energy and the breast phantom thickness [30]. Thus, we prefer to fit those values with a bi-variant polynomial depending on E and T [13]. The quotient form can be explained by the fact that the DgN coefficient is a fraction of the absorbed dose by the breast model to the air medium, and each of them can be a polynomial energy dependent function. For a relative difference between computed and fitted values that does not exceed 3% on the MGD calculation for the majority of X-ray beam spectra. We were limited to adopt the second order polynomial fitting Eq. (4). More precision needs extension to higher order, with the penalty of complicating the fitting expression form.

To our knowledge, fitting methods using the parameter E was treated similarly by many authors [7, 26, 30], whereas the introduction of the parameter T as a second variable within its actual form is proposed for the first time herein.

A. X-ray Quanta per mm² per mR computation

We calculated the photon fluence per exposure conversion coefficient, using Eq. (1), for photon energy of 1–120 keV. Figure 1 shows the relative difference between K(E) calculated based on the most recent data of mass energy absorption coefficient of air, provided by NIST [14], and those tabulated by Boone et al. [26]. The major difference seen concerns photons with energy lower than 75 keV and can reach 12% for 27 keV, as seen in the inset. For photon energies over 75 keV, such difference does not exceed 1%.

Fig. 1.Full-screen View

Comparison of the photon quanta per mm² per mR parameter calculated in this work (denoted by NIST) against Boone et al. [26] data (denoted by Boone 97). The inset refers to their relative deviation as function of the energy.

B. K(E) variation effect on Boone’s look-up table

Table 1 illustrates a comparison of MGD calculations between literature (Boone et al. [5]) and present study using the updated K(E) values. The comparison concerns three different anode/filter combinations of Mo/Mo (30 μm), W/Rh (50 μm) and W/Pd (50 μm) for 23–35 kVp and HVL of 0.269–0.597 mm of Al. The breast thickness was 4 cm and the glandularity was 50%. We observed a maximum relative difference between literature and Geant4 data of about 1.9, 5.0 and 4.2%, and between calculated (parameterized) and Geant4 data of about 2.6, 3.8 and 4.1%, respectively, for the above anode/filter combinations.

Referring to the Mo and W experimental spectra in Ref. [31], we observed that the two characteristic peaks were located around 17 keV and 11 keV, respectively, as shown in Fig. 2. The spectra covered 23–35 kVp. Also, Fig. 2 illustrates the curves of the relative error (%) of K(E), DgN(E) and their sum. According to Eq. (5) and Fig. 3, we are intended to find the MGD relative difference corresponding to the tungsten anode twice or higher than those for molybdenum. This can be explained by coincidence of the tungsten peak and the maximum relative deviation. Furthermore, we can guess that the relative deviation corresponding to the rhodium anode spectra will be lower than those for the molybdenum, for the studied kVp interval.

Fig. 2.Full-screen View

Relative difference curves of DgN(E), K(E) and their sum vs. E, between those computed in this work to those derived using old cross section data-base and Monte Carlo simulation code found by Boone et al. [5, 22]. The two energy intervals (grey bands) indicate the localization of the main peak of the tungsten (W) and the molybdenum (Mo) experimental spectra.

Fig. 3.Full-screen View

(Color online) Simulated DgN (mGy/mGy) versus energy (keV) and thickness (cm) for 50% glandularity, generated using the Geant4-based simulation program.

C. Data set generation and parameterization

The overall normalized glandular dose data were generated using Eqs. (2) and (3) and the developed simulation code, for breast glandularities of 100% and 50%. The data set corresponding to f_g=50% is shown in Fig. 3. Plots corresponding to f_g=100% are similar to Fig. 3. Table 2 illustrates the parameters resulting from fitting the overall DgN values and using Eq. (4) for energy interval of 11–120 keV. There is no need to subdivide the energy and thickness intervals, for parameterization purposes. A statistical test p-value of 0.9999 was seen for both cases.

D. MGD prediction ability

In order to investigate the possibility of MGD prediction using the current data set of DgN following the fitting method, we used Eq. (5) and the photon fluence per exposure spectrum corresponding to the anode/filter combinations of: Mo/Mo (30 μm), W/Rh (50 μm) and W/Pd (50 μm) generated according to the MASMIP and the TASMIP Boone’s codes [26] for different kVp and HVL X-ray tube parameters. The heel effect associated to the spectra was not taken into account for this study, similar to Ref. [16]. Thus the beam is assumed to have a uniform intensity and uniform spectral quality within the field covering the entire breast. If we took it into account, the MGD would vary as it is proportional to the fluence spectra.

From the data at 35 kVp in Table 2, we can see that MGD values become close to each other. This is due to the small contribution of DgN for low photon energy and high kVp. The results are in acceptable agreement with those directly simulated.

E. Interpolation and extrapolation capability

The DgN data values in Table 3 demonstrate effectiveness of the interpolation and extrapolation procedures, in terms of glandularity. The direct calculation using the linear interpolation and extrapolation technique and the Monte Carlo simulation were conducted for different cases of glandularity and the data are tabulated. We find an acceptable agreement between the two sets of data allowing the straightforward interpolation and extrapolation of the MGD for any given glandularity from about 0% to 100%.