Virtual dosimetry study with three cone-beam breast computed tomography scanners using a fast GPU-based Monte Carlo code

Objective. To compare the dosimetric performance of three cone-beam breast computed tomography (BCT) scanners, using real-time Monte Carlo-based dose estimates obtained with the virtual clinical trials (VCT)-BREAST graphical processing unit (GPU)-accelerated platform dedicated to VCT in breast imaging. Approach. A GPU-based Monte Carlo (MC) code was developed for replicating in silico the geometric, x-ray spectra and detector setups adopted, respectively, in two research scanners and one commercial BCT scanner, adopting 80 kV, 60 kV and 49 kV tube voltage, respectively. Our cohort of virtual breasts included 16 anthropomorphic voxelized breast phantoms from a publicly available dataset. For each virtual patient, we simulated exams on the three scanners, up to a nominal simulated mean glandular dose of 5 mGy (primary photons launched, in the order of 1011–1012 per scan). Simulated 3D dose maps (recorded for skin, adipose and glandular tissues) were compared for the same phantom, on the three scanners. MC simulations were implemented on a single NVIDIA GeForce RTX 3090 graphics card. Main results. Using the spread of the dose distribution as a figure of merit, we showed that, in the investigated phantoms, the glandular dose is more uniform within less dense breasts, and it is more uniformly distributed for scans at 80 kV and 60 kV, than at 49 kV. A realistic virtual study of each breast phantom was completed in about 3.0 h with less than 1% statistical uncertainty, with 109 primary photons processed in 3.6 s computing time. Significance. We reported the first dosimetric study of the VCT-BREAST platform, a fast MC simulation tool for real-time virtual dosimetry and imaging trials in BCT, investigating the dose delivery performance of three clinical BCT scanners. This tool can be adopted to investigate also the effects on the 3D dose distribution produced by changes in the geometrical and spectrum characteristics of a cone-beam BCT scanner.


Introduction
Breast cancer is the main cancer in the female population (Siegel et al 2020), accounting for 30% of all female cancers.The widespread reference method for early diagnosis is digital mammography (DM), a projective radiography technique which employs a low dose of x-rays.To overcome limitations associated with this exam, where the overlap of healthy tissues can hide the presence of a massive lesion-in particular for dense breasts with a high glandular fraction-3D imaging techniques have been developed, including digital breast tomosynthesis (DBT) (Sechopoulos 2013) and computed tomography dedicated to the breast (BCT) (Boone et al 2001, Sarno et al 2015, Zhu et al 2022a, 2022b).BCT provides fully tomographic imaging of the complex breast anatomy, with respect to the limited quasi-tomographic imaging performance of DBT, at dose levels comparable to, but higher than, DBT and two-view DM exams (Shim et al 2023).
Though proposed and described as a breast-specific cone-beam computed tomography (CT) technique more than 20 years ago (Boone et al 2001), and first assessed clinically 15 years ago (Lindfors et al 2008), the adoption of BCT in the clinic has been a slow process.Despite its intrinsic and demonstrated high potential as a high-quality breast imaging modality for breast cancer diagnosis, at dose levels comparable to two-view DM (Boone and Lindfors 2006, O'Connell et al 2010, O'Connell and Kawakyu-O'Connor 2012), BCT is still considered an emerging modality, and the total number of patients scanned so far is relatively limited, estimated in the order of several thousand, worldwide.BCT is an imaging technology capable of providing impressive breast anatomical detail and fidelity, with isotropic resolution in the uncompressed pendant breast, without using the firm breast compression adopted in DM and DBT exams, so largely avoiding patient discomfort.Moreover, it can also be coupled to iodine contrast imaging (Zhu et al 2022b) and to hybrid imaging with PET (Bowen et al 2009) or SPECT scanners (Russo et al 2016, Shah et al 2017).
Various clinical BCT prototypes (based on energy-integrating detectors) have been developed in the last years at University California Davis (named Albion, Bodega, Doheny, all adopting a cone-beam irradiation geometry) (Gazi et al 2015).To our knowledge, two commercial scanners are currently available: a cone-beam system by Koning Health USA (https://www.koninghealth.com)employing an energy-integrating flat panel detector, and a recent fan-beam system by Advanced Breast-CT GmbH Germany (https://www.ab-ct.com),employing a photon-counting detector and spiral scanning.A third commercial platform inspirated from UC Davis research is under clinical test, developed by Izotropic Imaging Corp. Canada (https://izocorp.com).Under laboratory investigation there is a new technical approach of narrow-beam clinical BCT (Ghazi et al 2022), as well as a phase-contrast BCT scanner with microfocus x-ray tube and a high-resolution CMOS flat panel detector (Mettivier et al 2021(Mettivier et al , 2022)).Research efforts were directed also toward phase-contrast BCT imaging (Auweter et al 2014, Heck andHerzen 2020), using a microfocus x-ray tube (Sarno et al 2016), or using synchrotron radiation-based photon-counting phase-contrast BCT (Bliznakova et al 2015, Mettivier et al 2016, Delogu et al 2017, Taba et al 2018, Brombal et al 2019); or toward 4D contrast-enhanced BCT (Caballo et al 2019).
Research on BCT technology is ongoing in several labs worldwide, with significant technical developments under investigation over those initial setups.While the acquisition geometry is similar in all BCT scanners (under-table rotating gantry with the patient in prone position on a suitably shaped horizontal table, and the uncompressed breast hanging from a large opening in the bed), usage of different scanner parameters (such as x-ray spectrum, focal spot size, image magnification, source-to-isocenter distance, cone-beam or spiral scanning, detector pixel pitch, energy-integrating or photon-counting detectors, broad beam or narrow beam geometry) characterizes each setup.This technology variety produces different results in terms of imaging performance and dose efficiency.However, direct comparison of scanners' performance for the same patients, all scanned on each of the available units, is not an option, also given the imaging radiation burden (corresponding to a mean glandular dose, MGD, of about 5-6 mGy or higher in a single BCT scan for an average breast) and the related radiation risk (Boone et al 2001(Boone et al , 2005(Boone et al , 2006)).
Clinical trials are fundamental steps in the design, validation, and implementation of new medical imaging technologies.These trials represent a huge and critical effort, due to ethical limitations, expenses, and time requirements.A possible alternative are virtual clinical trials (VCTs).Several MC codes have been developed for breast imaging studies, for simulating clinical images and estimating organ doses (Delis et al 2005, Hunt et Massera and Tomal 2021).A particular attention was dedicated to the possibility of a fast simulation process, with techniques such as variance reduction or the use of graphical processing units (GPUs) (Jia et al 2011, 2012a, Bert et al 2013, Barufaldi et al 2018).
In this framework, in BCT research, VCTs could represent a valuable tool for technology development, assessment, and optimization.Virtual imaging of the compressed or uncompressed breast is implemented using computer simulation techniques that closely replicate clinical exams in terms of breast anatomy, system geometry, x-ray spectrum, imaging detection process and glandular absorbed dose (Erickson et al 2015, Badano et al 2018, Bakic et al 2018, Baneva et al 2017, Badal et al 2021, Barufaldi et al 2021).In silico imaging clinical trials in 3D breast imaging have been investigated recently, to replace lengthy, expensive and dose-risk associated in vivo clinical tests, via development of several platforms for VCTs carried out with computational numerical methods (Elangovan et al 2014, Barufaldi et al 2018, Boita and Mackenzie 2019, Sharma et al 2019, di Franco et al 2020, Marshall and Bosmans 2022, Mettivier et al 2022).Validation processes can also imply realization of physical replicas of the digital phantoms and subsequent scans with clinical DM, DBT and BCT imaging units (Bliznakova et al 2016, Hernandez and Boone 2017, He et al 2019, Varallo et al 2022) or synchrotron radiation.
This work aims to realize a virtual study for a dosimetric comparison of three different BCT scanners (two developed by UC Davis group, and the third one, commercially available) using the VCT-BREAST platform for VCT (di Franco et al 2020) in which we implemented a version of the code gCTD dedicated to breast imaging (Jia et al 2012a(Jia et al , 2012b)).gCTD is a fast GPU-accelerated Monte Carlo (MC) simulation code developed at U Texas Southwestern.The virtual trial compared the three setups for BCT of the uncompressed pendant breast, in terms of uniformity of glandular dose delivery performance, average dose to the breast and skin dose, in the presence of different scanner-specific protocols, different x-ray spectra and irradiation geometry.The comparison was based on a set of patient-derived anthropomorphic virtual breasts (digital voxelized breast phantoms), each taking a virtual BCT scan on the three scanners simulated, at fixed MGD.

Materials and methods
The three investigated BCT scan systems are based on the characteristics of: (1) the UC Davis Bodega scanner (Scanner 1), (2) UC Davis Doheny (Scanner 2), and (3) the Koning Health commercial scanner (Scanner 3), with their characteristics as described in the available literature.Their setup was replicated in MC-based computer simulations implemented on a desktop computer (Dell Alienware Aurora R12, Intel Core i9, 5.3 GHz, 256 GB RAM, 4 TB SSD data disk) equipped with a dedicated GPU card (NVIDIA GeForce RTX 3090) for GPUaccelerated computing with a dedicated software.

Simulated scanners
During a BCT exam in pendant geometry, the patient is prone on the exam table and the breast to be imaged is freely hanging through a shaped aperture in the table, without compression.The x-ray tube and the detector rotate 360°under table, around the breast.The irradiation geometry adopted is a half cone-beam geometry that avoids directly irradiating the tissues not subject to examination, while capturing in each view the shadowgram of the whole breast (for axial BCT scanners) or of a thick coronal slice (for spiral BCT scanners).In this work, we simulated the characteristics of three cone-beam BCT scanners for which the geometric orbit of the x-ray source is circular and the views are acquired equally sampled over a full rotation, as described in the following.We note that a virtual dosimetric study has been reported recently for the spiral BCT geometry on a large real patient cohort, for the nu:view commercial scanner (AB-CT-Advanced Breast-CT GmbH, Erlangen, Germany) (Shim et al 2023).
The scanners' specifications here adopted are briefly described in the following: Scanner 1: This is the second generation BCT scanner proposed by the UC Davis research group (named Bodega).The detector is a scintillator-based flat panel (Paxscan 4030B Varian) with 0.194 mm pixel size.The x-ray spectrum has an inherent filtration of 1.58 mm Al and typically, 350 views are recorded at 80 kV in a full circular orbit in a scan (Boone et al 2005).
Scanner 2: It is the latest UC Davis BCT scanner (named Doheny) featuring a higher spatial resolution (3.6 mm −1 at 10% MTF, at the center of the FOV), thanks to a finer-pitch CMOS flat-panel detector (Boone et al 2001, Gazi and Boone 2014, Gazi et al 2015).The detector pixel pitch is 0.075 mm, but here we simulated a pixel size of 0.150 mm, for a total of 1933 × 1533 detector pixels.The pulsed x-ray beam (60 kV) is filtered with 0.22 mm Gd and 300 views per scan are simulated.
Scanner 3: This BCT scanner, commercialized by Koning Health (www.koninghealth.com),adopts a pulsed x-ray tube operated at 49 kV and the same detector as Scanner 1, but with less magnification, for 300 views per scan.
Table 1 shows the main characteristics of the three BCT scanners.The x-ray spectra of these three scanners were reproduced using the TASMICS code (Hernandez et al 2017).The heel effect (with cathode-anode direction along the chest wall-to-nipple direction) was not considered in the simulations.

Anthropomorphic voxelized phantoms
The anthropomorphic patient-derived phantoms used for this study were the high-resolution, segmented, voxelized phantoms from the dataset available in the Zenodo public repository for research purposes (Sarno et al 2021).In this repository are reported 150 3D voxelized phantoms of the uncompressed breast, obtained via a segmentation process (di Franco et al 2020) applied to tomographic datasets of real patients' BCT exams performed at UC Davis Medical Center, in which the voxels assume a value code for either air, skin, glandular or adipose tissue.Figure 1 shows an example of the axial, sagittal and coronal views of a virtual breast of this study (Phantom 1).
From the 150 phantoms present in this repository, all showing no sign of mass lesions, we selected 16 breast phantoms (in Supplementary material, table S1) with different values of glandular mass fraction (average value 16%, minimum value 3%, maximum value 29%), to evaluate the simulation performance in the case of low-, medium-, and high-density breasts.The glandular mass fraction of the phantoms was derived from the clinical BCT of real patients scanned at UC Davis, from their volume glandular fraction.In terms of whole breast volume, the chosen subset contains 4 breasts of small volume (less than 250 cm 3 , termed V1 class), 5 breasts of medium volume (between 250 and 600 cm 3 , termed V2 class), and 7 breasts of large volume (greater than 600 cm 3 , termed V3 class), with a minimum breast volume of 129.07 cm 3 , a maximum volume of 996.88 cm 3 , and an average volume of 486.61 cm 3 (table S1 in Supplementary material).The densities ρ air , ρ s , ρ a , ρ g , adopted for air, skin, adipose, or glandular voxels, respectively, were (1.225 × 10 −3 g cm -3 , 1.09 g cm -3 , 0.93 g cm -3 , 1.04 g cm -3 ).For each phantom, the glandular fraction by mass was calculated from the total number of glandular (N g ) and adipose (N a ) voxels, as:

Simulation code
The code used is based on gCTD (Jia et al 2012a(Jia et al , 2012b)), a MC code running on a GPU architecture under CUDA (version 11.6) programming environment (CUDA 2023), developed by the team of the University of Texas Southwestern Medical Center.The code was run on the personal computer hosting a single GPU card type NVIDIA GEFORCE RTX 3090, featuring 10496 CUDA cores, clock speed of 1.70 GHz and 24 GB of dedicated GDDR6X memory.Several thousands of CUDA threads running in parallel for each run permit a high level of acceleration with respect to multi-core CPU based architecture (Jia et al 2010).The GPU-based MC code simulates the acquisition of BCT projections, from which a CT dataset can be reconstructed with a cone-beam filtered backprojection algorithm; for each projection, and for each phantom slice, the code also returns the dose  distribution map (i.e., the absorbed dose for each voxel) inside the irradiated phantom, as well as in the irradiated air volume between source and detector.By summing, over all views, the dose map derived in a single projection, we then derived slice-by-slice the 3D dose map in a virtual tomographic scan.An example is shown in figure 2, which illustrates the axial, coronal, sagittal views, and 3D rendering of the dose map of Phantom 4. The absorbed dose was calculated (in μGy, then expressed in mGy) on a per-voxel basis by scoring the energy deposited in the (air, skin, adipose, gland) voxels by each interacting photon, divided by the mass of the skin, adipose, glandular, or air voxel type.For glandular voxels, hence, the scored dose represents the glandular dose in the volume of the voxel, and the average value of the frequency distribution of glandular dose values represents the calculated MGD for the given phantom, delivered in the virtual BCT scan (equation (2) below).

Dose evaluation
We simulated virtual exams of all selected phantoms on the three BCT scanners.For all simulations the exact number of photon histories in the whole MC-simulated BCT scan (always exceeding 10 12 primary photons launched in each simulation) was set to determine a nominal MGD of 5 mGy in each scan, as derived from the calibration procedure reported in the Results section.We adopted this value of the fixed MGD in all virtual scans, because this is close to the standard MGD value adopted at UC Davis for BCT scans (5-6 mGy for the average breast) (Hernandez et al 2020), where 6 mGy is twice the maximum MGD value for the average breast in a single screening mammography scan in the US legislation) (i.e., corresponding to the glandular dose for two view screening mammography of the average breast).Moreover, a MGD of 5 mGy is two times the maximum MGD in a mammography view for the standard breast in EU (2.5 mGy) (European Commission 2006), considering a limit MGD in BCT as that of a screening mammography exam with two views per breast.
From the 3D absorbed dose map, by masking the air, skin, and adipose voxels, we derived the 3D map of the absorbed dose in the glandular tissue (D g ) in the breast phantom, as well as the distribution of skin dose (D s ) and absorbed dose in the whole breast, all tissues included (D b ).The mean value (MGD) and the standard deviation (σ g ) of the glandular dose distribution was derived as the average value of the histogram in the resulting 3D glandular dose map in the phantom.The procedure was applied also to derive the mean skin dose (MSD), the mean dose value at the skin, and its standard deviation, σ s , as well as the mean breast dose (MBD) in the whole breast (skin included), and its standard deviation, σ b .This definition of the MBD differs from that of Shim et al (2023) who exclude the skin layer from the calculation of their mean absorbed dose value.
With respect to DM and DBT exams, BCT scans are known to produce a markedly more uniform distribution of radiation dose in the breast volume, and of glandular dose as well, due to the full rotation of the x-ray source around the breast in a cone-beam BCT scan, and to the approximate rotational symmetry of the uncompressed shape of the breast (Russo et al 2010, Sechopoulos et al 2010, Sarno et al 2022).As a spread index of the glandular dose distribution, we computed the percentage coefficient of variation COV dose (%) of the frequency distribution of the absorbed dose in a glandular voxel, D g i j k , , (mGy), calculated for each glandular voxel with coordinates (i, j, k) in the phantom matrix.Then, the MGD for the virtual BCT exam was calculated from the average of D g i j k , , over all N g glandular voxels in the phantom: The relative glandular dose is given by D g i j k , , /MGD.The percentage COV dose was calculated as the ratio of the standard deviation of the D g i j k , , distribution in the phantom, σ g , to the average value of the glandular dose distribution (i.e., the MGD): the standard deviation, σ g , representingthe spread of the dose distribution: For a given breast, virtually scanned on one of the BCT scanners here investigated, a lower COV dose value indicates a more uniform distribution of the dose delivered to the glandular tissue during the exam.This value can be used also as a metric for analysing radiation-risk models of energy deposition in breast imaging exams, not based on the sole MGD metrics (Oliver and Thomson 2019).For a given breast (assumed positioned always centrally in the scanner FOV), the glandular dose distribution resulting from a BCT scan, at a fixed MGD, is dependent on the geometry of the scanner and on the x-ray spectral shape.Hence, the COV dose is adopted here as a figure of merit for comparing the effect on the uniformity of glandular dose delivery, due to the different protocols and hardware solutions of the scanners' implementation.

Air Kerma and D g N CT evaluation
In BCT dosimetry, the MGD for a given breast (size and glandular fraction) is proportional to the incident air kerma at isocenter, K air , via a MC-derived normalized glandular dose coefficient, D g N CT : Tables of D g N CT values are usually calculated with the homogeneous tissue approximation (i.e., where the phantom tissue is a homogeneous mixture of adipose and glandular tissues of given glandular fraction by mass), but they have been calculated also for anthropomorphic voxelized breast phantoms, where the voxels are either 100% adipose or 100% glandular, spatially distributed realistically in the breast phantom at a given value of glandular fraction (Hernandez et al 2017).For DM and DBT dosimetry, the former approach was also adopted in Massera et al (2021aMassera et al ( , 2021b)), Sarno et al (2021), Caballo et al (2022).Here, we adopt the latter method (as in Sarno et al 2022 and Mettivier et al 2022), using our set of patient-derived anthropomorphic breast phantoms.
By running a MC simulation with the given scanner setup but with only air (no phantom) in the scanner FOV, we calculated the voxelized 3D map of air kerma within the whole scanner FOV, for each setup.From these 3D maps we calculated K air by scoring the energy deposition in a 3 cm 3 volume simulating a 100 mm long CT ion chamber placed at the axis of rotation of the scanner.We used this value to evaluate the normalized glandular dose coefficient, D g N CT (mGy/mGy), according to equation (5).The knowledge of simulated 3D maps of air kerma in the scanner FOV-including the air kerma calculated for each air voxel of the phantom volume, K , for voxel coordinates (i, j, k) inside the phantom matrix-permitted also to derive, for each in silico BCT scan, 3D maps of the voxel-wise normalized glandular dose coefficient, D N , , , defined as: distribution of glandular dose in the specific breast phantom scanned: We then compared this value with that calculated with equation (5).
We note that the MGD is independent of the position of the breast in the scanner FOV-given the current MGD metric for BCT dosimetry which is based on equation (5) via spectrum-, breast size-, and glandularityspecific, pre-calculated D g N CT values for the given scanner geometry, for a given MGD to be delivered to the breast.Indeed, the air kerma K air is measured at isocenter during a BCT scan without the breast in the FOV.On the other hand, since the 3D map of K i j k air , , values is not spatially uniform over the whole scanner FOV, positioning the pendant breast in different lateral locations of the FOV might produce-for a given K air adopted in the (real) scan-a different distribution of the D g i j k , , values, and hence, a possible different value of the MGD.To investigate this issue, we performed a virtual scan of one phantom (no. 10 in table S2 in Supplementary material) shifted laterally by 2.5 cm in the coronal slices, with respect to the original position determined by the clinical CT scan.

N ph -to-MGD conversion factors
For each scanner, the number of primary photons (n ph ) per projection, necessary to obtain a given MGD value in the scan (here fixed to a nominal value of 5 mGy) in the simulated exams, depends on the number of projections (N views ), on the scanner's specification (geometry and x-ray spectrum) and on the characteristics of the breast, and specifically, on the 3D glandular distribution inside the breast.The total number of primary photons, N ph , launched in the simulation, is then equal to n ph × N views .We performed complete breast CT simulations with 10 8 , 10 9 and 10 10 photons/projection, for a virtual breast CT scan of each phantom with all scanners.These simulations were repeated for the phantoms no. 1, 2 and 3.

Imaging detector
While this work was focused on deriving 3D dose maps, rather than providing simulated BCT imaging datasets, within the VCT-BREAST platform the detector characteristics can be modeled effectively post-simulation.In each virtual scan the software provides in output both the map of absorbed dose in each voxel of the phantom, and the set of 300 (Scanners 2 and 3) or 350 (Scanner 1) raw projections.All three scanners adopt an energyintegrating, scintillator-based flat panel detector, with 0.075 mm (Scanner 2) or 0.197 mm pixel size (Scanners 1 and 3).However, a photon-counting detector can be simulated as well: as an output option, the software permits to record the energy of each single photon reaching any detector pixel, rather than recording the integrated energy of all photon incident on the pixel.For the present dosimetric study where no BCT virtual images are analysed, the imaging detector was simulated as an ideal energy integrating detector: the code scores, for each projection, the total energy deposited in each detector pixel.

N ph -to-MGD conversion factors
The MGD, σ dose , COV and MGD/N ph (N ph is the number of primary photons) values calculated from MC simulations, at increasing values of the photon histories per projection, for Phantoms 1, 2, 3 and 10, were calculated and reported in table S2.The number of primary histories determined a maximum uncertainty in the MC estimates (standard deviation of the mean) of about 10% (for 3 × 10 10 photons/scan), about 3% (for 3 × 10 11 photons/scan), and about 0.3% (for 3 × 10 12 photons/scan), based on the uncertainty estimation method of replicating 10 times the simulation with 1/10 of the photon histories and then calculating mean value and corresponding standard deviation of the mean of the repeated measurements (i.e., adopting the batch method, Walters et al 2002).From the data reported in table S2, in Supplementary material, we derived the N ph -to-MGD conversion factors as a linear fit coefficient through a linear fit of the data and, hence, the number of primary photons, N ph , to launch in the simulation for determining a fixed MGD (5 mGy nominal value) (Boone et al 2005, 2006, European Commission 2006) for each phantom and for each scanner (table 2).

Air kerma evaluation
As described in section 2.5, the air kerma at isocenter (K air ) was calculated in dedicated simulations, by scoring the energy deposition in a 3 cm 3 volume simulating a 100 mm long, 3 cm 3 CT ion chamber, placed at the axis of rotation of the scanner FOV.For this measurement we simulated only the three scanner setups but with just air (no phantom) in the FOV.The results are reported in table 3. The simulations were realized with the N ph necessary for a nominal MGD of 5 mGy (table 2).

Glandular dose distribution
As indicated in section 2.4, the simulation code provided the 3D dose distribution map inside the irradiated breast volume, keeping the voxel size of the corresponding voxelized phantom and CT image dataset.As an example, in figure 3(a) are reported the dose distributions for Phantom 1 in a central coronal, axial, and sagittal slice, for the three different scanners (nominal MGD fixed at 5 mGy), as well as the 2D dose profiles measured along the white line in the axial slices.As expected, based on the corresponding interaction coefficients, figure 3(a) shows that the highest dose is absorbed in the skin layer, and that the dose to the glandular tissue is higher than the dose to the surrounding adipose tissue.The cupping in the dose profile in coronal slices is evident, indicating the 'shielding' effect of outer layers (including skin) over tissues located at the center of the breast.In addition to showing that the dose decreases inward from the skin to the internal breast tissues, figures 3(a), (c) also show that it increases very significantly in the longitudinal direction from the chest wall to the nipple (up to 100% in the slice shown).This last trend is clearly shown in figure 4 for the example case of Phantom 4 and Scanner 2, as well as in axial and sagittal slices of figure 3 for Phantom 1 and Scanner 3.For all tissues, the dose in the longitudinal direction increases steadily from the chest wall towards the nipple, with the glandular tissues absorbing significantly higher dose than the MGD (e.g., up to the nipple, see also figure 2).This trend is quantitatively and qualitatively similar for all scanners and is likely related to the decrease of the tissue thickness in coronal planes, in the pendant breast.Moreover, figures 3(a), (b) show that Scanner 2 produces a dose distribution (in the coronal plane) very close (within few percent) to that determined by Scanner 1.On the other hand, the dose gradient for Scanner 3 is higher than for the two other scanners, with significantly higher dose to the skin (up to +30%), but with less dose to the adipose and to the glandular tissue located at the center of the breast and, correspondingly, higher dose for tissues located at the periphery of the breast, towards the skin layer.These global features of the 3D dose map are common to all virtual exams on the three BCT scanners, as shown in figure 5 for the central axial dose distributions for the Phantoms 1, 2, and 3, obtained with the three scanners.The comparison of 3D dose distributions for the glandular tissue only, for the three scanners, is shown in figure 6 in the case of Phantom 1.
A paired comparison of the glandular dose distribution for Scanner 2 (60 kV spectrum) vs. Scanner 1 (80 kV spectrum), and for Scanner 2 vs. Scanner 3 (49 kV spectrum), respectively, is reported in figure 6  to ±30%) for the case of Scanner 2/Scanner 3, and little spread (up to ±10%) for the case of Scanner 2/ Scanner 1.
As shown in figures 3 and 6, the dose delivered by the Scanners 1 and 2 is distributed more uniformly in the breast than for Scanner 3. A comparison of glandular dose histograms for breast Phantoms 1, 2 and 3 is reported in figures 7(a)-(c), showing that for Scanner 3 the maximum glandular dose can largely exceed the MGD, reaching values higher than 1.5 times the MGD.When scanning the three breasts on Scanner 1 (figure 7(d)), the maximum dose to the glandular voxels was 1.3 times the MGD.
The mean (MGD), standard deviation (σ dose ) and COV (100 × σ g /MGD) of the glandular dose distributions for each phantom in the dataset, are reported in table S3 in Supplementary material.Note that the average value of the MGD calculated over the cohort of breasts is constant: 5.005 ± 0.008 mGy for Scanner 1, 5.008 ± 0.007 mGy for Scanner 2, and 5.014 ± 0.012 mGy for Scanner 3. On the other hand, the distributions of dose delivered by Scanner 1 and 2 are peaked around this average value (Scanner 1: ̅ g s = 0.531 mGy, COV dose = 10.62%;Scanner 2: ̅ g s = 0.576 mGy, COV dose = 11.50%), while the distribution for Scanner 3 is wider ( ̅ g s = 0.958 mGy, COV dose = 19.10%),with a significant increase in COV .Data for MGD and COV contained in table S3 in Supplementary material, over the entire virtual cohort, are also reported in figure 8 as a function of the volume class (V1, V2 and V3).We note that with respect to a nominal MGD of 5 mGy for all scanners and all scans, the actual MGD in the set of virtual scans varied by only 0.01 mGy (1 std.dev.) (table S3).We note an increase ̅ COV dose values with an increase of breast volume for any given scanner, in the approximate range 8%-15%, 9%-17%, and 16%-25%, for Scanner 1, Scanner 2 and Scanner 3, respectively.

Whole breast and skin dose distributions
We analyzed the statistics of absorbed dose in the whole breast volume, and in the skin, adipose and glandular tissues.In table S4 and S5 and infigures 9 and 10 are reported the calculated mean (MBD and MSD), percentage difference (MBD-MGD)/MGD (Supplementary material, table S4) and percentage difference (MSD-MGD)/ MGD (Supplementary material, table S5) ratio of the 3D distribution of dose in the whole breast (MBD) and in the skin (MSD), respectively, obtained for all phantoms for the three scanners.In figure 9, the analysis is shown separated for the three volume classes V1, V2 and V3 for MBD and MSD, respectively.Figure 9 shows that the MBD is lower than the MGD by 10%-30%.On the other hand, the MSD is significantly higher than the MGD, with the Scanner 3 providing the highest values of the MSD/MGD ratio (figure 10).Specifically, for the V1, V2 and V3 volume classes, the MSD is 50%-80% higher than the MGD for Scanner 3, while it is only 10%-22% higher than the MGD for Scanner 1 and Scanner 2.
The same indications can be deduced, as an example, from figure 11, which shows the histograms of dose values for all the voxels and for the three scanners, in the case of whole breast Phantom 1 (figure 11 11(a)) shows three main peaks which, in the order of increasing dose, can be attributed, respectively, to the adipose voxels, to glandular voxels and to the skin voxels.In the case of Scanner 3 the dose distribution of glandular and adipose voxels overlap, to some extent, at variance with Scanner 1 and Scanner 2 (figures 11(b), (c)).

Dose distribution for the offset phantom
A virtual scan was realized for the Phantom 10 shifted laterally by 2.5 cm in the coronal slice with respect to the original position in the scanner FOV.The goal was to study the influence of the spatial non-uniformity of the 3D air kerma distribution in the scanner FOV, on the MGD.This might highlight any dosimetric effect of translational shifts in the breast position.Figure 12 shows the 3D dose distribution of Phantom 10 in the original position, and of the phantom 10 shifted laterally.In the figure the circle shows the center of rotation.The measured values of MGD, σ dose and COV for both distributions were reported in table 4. The simulations were  realized by simulating 10 9 photons per projection.The data in table 4 show that there is only a slight difference (less than 2.9% and 0.8%), in terms of MGD and COV values, respectively.

D g N evaluation
As described in section 2.5, D g N values for Phantoms 1, 2 and 3 with all three scanners were evaluated using two different methods.In the first method, we calculated K air by scoring the energy deposition in a 3 cm 3 volume simulating the sensitive volume of a 100 mm long CT ion chamber placed at the axis of rotation of the scanner.In the second method, we calculated the 3D maps of the voxel-wise normalized glandular dose coefficients.Table S6, in Supplementary material, reports the results obtained with the two methods.Figure 13(a) shows a comparison between the results obtained with the two different methods.We noted a difference of less than 4% between the two evaluation methods (slope of the linear fit equal to 1.038 ± 0.014), with the D g N estimates from 3D maps slightly higher (3.8%) than estimates derived from K air .A comparison with the D g N values reported in Hernandez et al (2019) is shown in figures 13(b) and (c), for estimates from K air and 3D dose maps, respectively.In these figures we observe a difference of about −25% for both methods (with values in this study less than the values in Hernandez et al 2019).Table 4. Calculated values of mean (MGD), standard deviation (σ dose ) and COV dose (100 × σ g /MGD) of the 3D distribution of glandular dose in Phantom 10, for the original scan and for the laterally shifted (by 2.5 cm) position in the scanner FOV.As many as 10 9 primary photons/projection were launched in the simulation.

Discussion
In diagnostic imaging dedicated to the breast, simulation techniques based on MC codes are increasingly adopted in VCT platforms (Elangovan et al 2014, Badano et al 2018, Bakic et al 2018).A noteworthy example is the recently completed VICTRE project (Virtual Imaging Clinical Trials for Regulatory Evaluation) (Badano et al 2018, Badal et al 2021), a large VCT study started by FDA in USA to virtually replicate the clinical trial of an FDAapproved clinical DBT system (FDA 2019).Our team is involved in the VCT-BREAST project, an ongoing MCbased VCT project for performing in silico x-ray examinations of the breast (Mettivier et al 2022).Simulations replicate the 2D or 3D exams and provide absorbed dose maps in the phantom (as well as image projections); they run on a GPU OpenCL based cross-platform adopting GPU accelerating engines for simulation of the imaging datasets and 3D maps of the glandular dose.Patient-derived, anatomically realistic, computational phantoms of the uncompressed breast were generated using a tissue classification algorithm on BCT images, as described elsewhere (Sarno et al 2021).In all virtual scans shown in this study, the dose decreases steeply at the skin-adipose tissue interface (figures 2, 3, 5).The larger gradient at the breast periphery occurs for the lower beam energy scan (Scanner 3, 49 kV), while the dose gradient decreases only slightly from the scan at 60 kV (Scanner 2) to that with Scanner 1 (80 kV) (figure 3(b)).The skin and the subcutaneous fat layer act as a shield for tissues at the core of the breast volume (see also figure 6(a)), by absorbing dose deposited by the incident radiation in the outer layers.The rotational geometry of the BCT scan distributes the dose radially around the center of rotation in the coronal slices (parallel to the chest wall), producing the characteristic dose cupping shown in all BCT dose maps (see figure 3(b)).This has an implication on the glandular dose map, since glandular tissue located closer to the axis of rotation receives a lower dose with respect to glandular tissue at the periphery.Indeed, for a given breast shape and fixed glandular fraction, and fixed air kerma at isocenter, the MGD in a clinical BCT scan depends on the actual spatial distribution of the glandular tissue in the patient's breast.
Another observation regards the dose distribution from the chest wall to the nipple: the related reduction of the breast diameter, due to the tapered shape of the pendant breast, produces an increase of the absorbed dose, for all three scanners (and beam energies), at a given value of incident air kerma at isocenter.Hence, the glandular tissue located in breast regions close to the nipple may receive a dose significantly higher than the MGD, in a cone-beam BCT scan.In cone-beam BCT, a way to attenuate this phenomenon could be to adopt a shaped beam compensator (bow-tie filter) with increased thickness of attenuation material from the side of the chest wall to the side of the nipple (i.e., in the cone direction, in addition to increasing the filter thickness in the fan direction as common for fan-beam CT scanners).Alternatively, one could adopt a shaped breast holder (as devised, e.g., by Sarno et al 2016), with a shape that could compensate for the 'missing' tissue towards the nipple.However, though bow-tie filters for each breast diameter have been tested for cone-beam BCT (Boone et al 2004), they are not commonly adopted for clinical scans, to our knowledge.In spiral BCT, on the other hand, where the same phenomenon occurs, a solution might consist in air kerma modulation, e.g., by reducing the tube current and/or the tube kilovoltage in the scanner rotation, when translating the scanner gantry vertically.A dedicated virtual study simulating the use of such scanning strategies in spiral BCT to compensate for this non-uniform glandular dose distribution, with possible reduction also of the MGD, is outside the reach of this study limited to cone-beam CT, but represents a foreseen addition of the VCT-BREAST GPU-accelerated platform.
Using COV dose as a figure of merit, this dosimetric virtual study showed that the glandular dose is distributed more uniformly with large breasts (V3 class, table S3) in the investigated dataset of phantoms (breast volume in the approximate range 100-1000 cm 3 and glandular fraction by mass in the range 2%-30%).Moreover, for a fixed MGD of 5 mGy for all virtual scans, the dose delivered is more uniformly distributed by scanning with Scanner 1 (80 kV) and 2 (60 kV) than using Scanner 3 (49 kV).Data shown in table S3 and figure 8 show that the mean value for σ dose and COV dose were, respectively, 0.532 ± 0.207 mGy and 10.62% ± 4.14% (Scanner 1), and 0.576 ± 0.210 mGy and 11.5% ± 4.19% (Scanner 2), while for Scanner 3, σ dose and COV dose were, respectively, 0.958 ± 0.296 mGy and 19.10% ± 5.90%.Furthermore, Scanner 3 provides a much higher dose value to the skin (≈ 8 mGy, for 5 mGy MGD) than that provided by other scanners (figures 10 and 11).From the values reported in table 4 and figure 12, we also noted that in the case of a laterally shifted phantom position (by 2.5 cm) in the FOV scanner, there is no practical difference in terms of MGD and COV values for Scanners 2 and 3, whereas a slight difference (2.9% in MGD) is present for Scanner 1.
As regards the comparison of the D g N values for BCT calculated in this study, and those reported e.g. by the UC Davis group, we provided a comparison with their recent D g N data in BCT obtained with heterogeneous breast phantoms (Hernandez et al 2019).When normalizing the MGD to the incident air kerma evaluated at isocenter, we observed a discrepancy in the order of 25%.In this respect, we note that the values for Scanner 3 simulated in Hernandez et al (2019) were obtained with an added filtration of 1.5 mm of Al, while a filtration of 1.76 mm of Al was adopted in this study, this representing a possible origin of the observed mismatch between the two methods.
These results can be promptly extended to a larger dataset of voxelized uncompressed breast phantoms, than the limited dataset (16 breast models) investigated here.This work was instrumental in testing and validating the whole procedure for 3D virtual dosimetry in BCT within the VCT-BREAST platform.A useful extension foreseen for VCT-BREAST is the digital reproduction of the photon-counting spiral BCT scanner developed by Advanced Breast-CT company, in line with recent reports (Germann et al 2021, Shim et al 2023).
The simulations-implemented on a single NVIDIA GeForce RTX 3090 graphics card embedded in a tower desktop personal computer, carried out by launching as many primary x-rays as to deliver an average glandular dose of about 5 mGy in the virtual BCT scans-allowed to obtain a complete imaging and dosimetric virtual study in about 3.0 h on a single high-performance GPU.This corresponds to about 3.6 s of total GPU run time per 10 9 primary photons launched in the simulation (including disk writing time for memorization of 3D dose maps as well as of projections for imaging); this value reduces by one half if projections are not saved to disk.This can be compared with Bert et al (2013), who report 650 ms/10 6 primary photons (i.e., 650 s for 10 9 photons) for cone-beam transmission tomography with their GPU code running on an NVIDIA GTX580 GPU card with 512 cores operated at a clock frequency of 1.23 GHz.This implies a reduction factor of 180 times for our gCTD code running on a NVIDIA RTX3090 GPU card with 10496 cores and a clock speed of 1.70 GHz, with respect to Bert et al (2013).The comparison of the hardware/software architecture of this work is also of interest with respect to Ghammraoui and Badal (2014), who report 2.9 × 10 8 x-rays/s (i.e., 3.4 s for 10 9 photons) processed by their MC-GPU code running in parallel on a cluster of 8 NVIDIA GTX580 cards.gCTD is also more computationally efficient (by a factor 10 2 ) with respect to the gMCDRR code of Qin et al (2022) written for cone-beam CT scatter correction, which run on the same hardware GPU platform.Using a single NVIDIA Quadro P5000 GPU card and the ImpactMC commercial MC dose simulation software, Shim et al (2023) report a computation speed equal to 1.1-2.6 × 10 7 primary photons s -1 for their simulations on patient-derived anthropomorphic breast phantoms, which is a performance two orders of magnitude worse than in our study.In their paper, Shim et al (2023) report a maximum statistical uncertainty of 1% when launching 2 × 10 7 photons per projection with 2000 projections per rotation (i.e., 4 × 10 10 photons in each virtual scan), against about 4 × 10 12 photons launched in each virtual scan of this study (table 2).Hence, by reducing by two orders of magnitude the number of photons launched in each simulation, the fast VCT-BREAST tool can produce 3D dose maps within less than 1 min of GPU time.Once the personalized breast phantom of the examined patient is ready (after the scan), obtained by application of a fast segmentation methods (e.g., within 12.8 min of average computational time per a spiral BCT scan when no silicone implant is present, as reported by Shim et al 2023), this performance will pave the way for real-time personalized BCT dosimetry after the scan.In this case, any mismatch between the prescan estimated MGD and the after-scan calculated MGD could be evidenced, and the actual personalized MGD estimate could be determined.Since computational methods exist for compressing the patient's digital 3D breast phantom (e.g., Sarno et al 2021), the above procedure might be adopted also for a personalized MGD estimate also in successive DM or DBT scans of the same patient.

Conclusions
A virtual dosimetric study was carried out in BCT, employing 16 virtual patients and the VCT-BREAST platform for in silico investigations, a GPU-based MC simulation tool for the realization of virtual imaging and dosimetric trials in the detection and diagnosis of breast cancer.The goal was to investigate the dosimetric performance of three scanners clinically used for BCT examinations.Computationally, the GPU-accelerated VCT-BREAST platform was highly efficient, returning a complete CT dataset and a 3D absorbed dose map within 3.0 h (with about 0.3% statistical uncertainty on voxel dose values, and 10 10 primary photons per projection), when delivering a virtual MGD of about 5 mGy, i.e., in the order of a typical BCT scan with clinical units.This timing included the production of the full dataset of image projections (not analysed in this dosimetric study), in addition to the 3D dose maps for skin, adipose and glandular tissues.Using COV dose as a figure of merit for homogeneity of glandular dose delivery, our virtual trial showed that, when the breast is imaged with the three different clinical scanners, similar uniform distributions of the glandular dose can be obtained with both Scanner 1 (operated at 80 kV) and Scanner 2 (operated at 60 kV).On the other hand, for Scanner 3 (49 kV) the dose gradient is larger and the dose is less uniformly distributed than with the other two scanners.The possibility of real-time personalized dosimetry in BCT, as well as in DM and DBT, is envisaged, when coupling this fast simulation performance with fast segmentation and digital compression methods already available, for producing a personalized computational breast phantom soon after the BCT exam.

Figure 1 .
Figure1.Central slices in axial, sagittal and coronal views of the uncompressed voxelized breast Phantom 1 (18.14%glandular fraction by mass).The voxels in the images correspond to air (black), skin (white), adipose (dark gray), and glandular (light gray) tissue, respectively.The phantom has 32 million voxels of size 0.3663 × 0.3663 × 0.2026 mm 3 , and a total volume of 666 cm 3 (see tableS1in Supplementary material).

Figure 2 .
Figure 2. The coronal, axial, sagittal, and volume rendered views of the 3D distribution of absorbed dose (mGy) in a central slice of the Phantom 4, after a virtual BCT scan with Scanner 3. Phantom 4 has 410 × 410 × 411 voxels of size 0.3490 × 0.3490 × 0.2026 mm 3 , with a glandular fraction by mass of 29.53%.The false-colour scale shows the decreasing levels of absorbed dose in going from the skin to the glandular and to the adipose tissue, respectively.

º
This metric considers the spatial variation of the air kerma in the scanned volume.Then, we calculated the 3D voxel map of D N g i j k CT , , for all glandular voxels with equation (6), and from its frequency distribution we derived the mean value, and the standard deviation, σ DgN ) of the D (a), by showing the ratio of the two corresponding maps of relative dose (i.e., dose per voxel in mGy divided by the MGD in mGy) (shown in figure 6(b)) in the selected central coronal slices of Phantom 1.Here, pixel values different from unity show differences in the dose distribution: the histograms shown in figure 6(b) indicate a large deviation (up

Figure 3 .
Figure 3. (a) Example of coronal (upper line), axial (middle line) and sagittal (lower line) slices of simulated 2D dose maps (mGy) showing the dose delivered to the breast tissues (skin, adipose and glandular) in Phantom 1, after a virtual BCT scan on the three scanners, at a fixed nominal MGD of 5 mGy.The highest dose values occur in the skin and in the glandular tissues.Dose decreases from the periphery to the center in coronal slices (dose cupping), and increases from the chest wall to the nipple, in axial and sagittal slices.(b) Dose profile (all tissues) along the line indicated in the three coronal slices.(c) Glandular dose profile in the longitudinal direction, along the white line shown in the axial profiles.In the three scans at different beam energies, the glandular tissues absorb more dose than the surrounding adipose tissue, which generates the 'absorbed dose contrast' seen in the volume dose maps. dose

Figure 4 .
Figure 4. Coronal slices (2.026 mm thick) of the 3D distribution of absorbed dose in the phantom tissues, at various axial depths in the direction from the chest wall to the nipple, for Phantom 4 imaged on Scanner 2. The calibration scale is in mGy per voxel.One can observe that the average glandular dose in the coronal slices increases from the chest wall to the nipple.The nominal MGD was 5 mGy.

Figure 5 .
Figure 5. Example of 2D dose distribution (mGy per voxel) in a central coronal slice of Phantom 1 (upper row), Phantom 2 (central row) and Phantom 3 (lower row), obtained by simulating a scan with a nominal MGD of 5 mGy with Scanner 1 (first column), Scanner 2 (second column) and Scanner 3 (third column).

Figure 6 .
Figure 6.(a) Relative glandular dose map for a central coronal slice of Phantom 1, scanned on the three BCT scanners.The relative glandular dose distribution ranges from 0.7 to 1.3, around the average value of 1.0.(b) Ratio of relative glandular dose maps for the two compared scanners: Scanner 2/Scanner 1, and Scanner 2/Scanner 3. The two dose ratio distributions have.For such a comparison, the histograms of the ratio of relative glandular dose values for all voxels in the phantom (lower right plot) show almost the same average value, but a largely different range of min-max values, with a three-times larger spread for the ratio Scanner 2/Scanner 3, with respect to the ratio Scanner 2/Scanner 1.

Figure 7 .
Figure 7. Distribution of dose values for glandular voxels, for a virtual BCT scan on Scanner 1 (thin black line), Scanner 2 (red line) or Scanner 3 (thick black line) at a fixed nominal MGD of 5 mGy, for Phantom 1 (a), Phantom 2 (b) and Phantom 3 (c).In (d) are shown the distributions for Phantom 1 (thin black line), Phantom 2 (red line) and Phantom 3 (thick black line), obtained in a virtual scan with Scanner 1.Note the large spread of dose values around the MGD value, for virtual scans on Scanner 3, with respect to the two other scanners.

Figure 8 .
Figure 8. Box plot showing mean, median line and values range for phantoms belonging to the V1 (small), V2 (medium) and V3 (large) classes of whole breast volume, of (a) mean glandular dose (MGD) and (b) COV (100 × σ g /MGD) of the 3D dose distribution, evaluated for each scanner.The ensemble mean MGD for all scans was 5.00 ± 0.01 mGy.

Figure 9 .
Figure 9. Box plot showing, for each scanner, the mean, median line and range in V1 (small), V2 (medium) and V3 (large) breast volume classes, of the mean average dose (MBD) of the 3D distribution in the whole breast; the corresponding percentage difference (MBD-MGD)/MGD is indicated on the right-hand scale.

Figure 10 .
Figure 10.Box plot showing, for each scanner, the mean, median line and range in V1 (small size breast), V2 (medium size breast) and V3 class (large size breast) of the mean skin dose (MSD) of the 3D distribution in the phantoms; the percentage difference (MSD-MGD)/MGD is indicated on the right-hand scale.

Figure 12 .
Figure 12.Example coronal slice of the 3D dose distribution (in mGy per voxel) simulated for Phantom 10 with Scanner 1, (a) in its original position in the scanner FOV, and (b) after a rigid lateral horizontal shift of 2.5 cm of the phantom within the scanner FOV.The white circle indicates the scanner isocenter; a square grid was overlaid for visual reference.

Figure 13 .
Figure 13.Comparison between (a) the D g N values calculated calculating K air in a 6 cm 3 air volume at isocenter (indicated as D g N from K air ), and the D g N values calculated simulating K air in each voxel (indicated as D g N from 3D); (b) the D g N from K air and the D g N values reported in Hernandez et al (2019); (c) D g N from 3D data versus D g N values reported in Hernandez et al (2019).

Table 1 .
table S1 in Supplementary material).System parameters of simulated cone-beam BCT scanners.

Table 2 .
N ph -to-MGD conversion factors (mGy/photon), number of primary photons n ph per projection, and total number of photons per scan N ph , launched in the simulations for determining a nominal MGD of 5 mGy in the virtual scan of sample Phantoms 1, 2 and 3.

Table 3 .
Air kerma (mGy) (mean ± std.dev.) calculated in a 3 cm 3 cylindrical volume simulating a 100 mm-long CT ion chamber, placed at the axis of rotation of the scanner FOV.The simulations were run with the number of photon histories/projection, n ph , necessary for a nominal (virtual) MGD of 5 mGy for each scanner, as indicated in table 2.