Development of a digital breast phantom for photoacoustic computed tomography

: Photoacoustic (PA) imaging provides morphological and functional information about angiogenesis and thus is potentially suitable for breast cancer diagnosis. However, the development of PA breast imaging has been hindered by inadequate patients and a lack of ground truth images. Here, we report a digital breast phantom with realistic acoustic and optical properties, with which a digital PA-ultrasound imaging pipeline is developed to create a diverse pool of virtual patients with three types of masses: ductal carcinoma in situ, invasive breast cancer, and fibroadenoma. The experimental results demonstrate that our model is realistic, flexible, and can be potentially useful for accelerating the development of PA breast imaging technology.


Introduction
Worldwide, breast cancer is the leading type of cancer in women, contributing 25.4% of the total number of new cases diagnosed in 2018 [1]. It has been proved that the key to breast cancer survival is early diagnosis [2][3][4]. Mammography is currently the gold standard for breast cancer screening which, however, has lower sensitivity in women with dense breasts [5,6]. Ultrasound (US) is seen as an adjunct to mammography, but suffers from speckle artifacts and low specificity [7,8]. Magnetic resonance imaging (MRI) poses a large financial burden and sometimes requires the use of intravenous contrast agents that can potentially cause damage to the human body [9][10][11]. It is also hampered by often excluding patients due to claustrophobia, pacemakers, etc. In recent years, photoacoustic computed tomography (PACT) is growing rapidly as a non-invasive imaging modality for breast cancer screening and diagnosis, which has overcome many of these limitations [12][13][14][15][16][17].
In PACT, the biological tissue being imaged is illuminated by a pulsed laser. Scattered photons are selectively absorbed by chromophores in the tissue, giving rise to fast thermal expansion, which generates ultrasonic waves [18]. The photoacoustically generated sound waves are then received by ultrasonic transducers and digitally processed to reconstruct images of light absorption contrast. Since scattered photons are employed for contrast generation while ultrasound waves are used for signal localization, the imaging depth of PACT can reach several centimeters, way beyond the optical diffusion limit of 1 mm in typical biological tissues [19,20]. By scanning the wavelength of the illumination light, the absorption spectrum of each image voxel can be measured, which contains important functional information, such as saturation of oxygen (sO 2 ) [13,16]. To date, several clinical PACT systems for breast imaging have been built with vastly different configurations and image characteristics [12][13][14][15][16][17].
For the following reasons, in silico studies are important for continued development and improvement of PACT imaging devices: (1) In PACT, the properties of the signals are coupled with image features such as blood vessels' orientation and spatial density. A computer-generated phantom with various anatomically and physiologically relevant properties can help engineers optimize system designs. Moreover, building an imaging system and doing clinical experiments are costly and time consuming. It can be very inefficient to iteratively optimize the system via real experiments. In contrast, system design based on numerical phantoms is fast and cost-effective. (2) The digital phantom provides a gold standard for both anatomical and functional imaging, thus making system evaluation more quantitative. (3) A fair comparison among different imaging systems requires that the same object be imaged and evaluated, thus all systems need to be present on site to image the same patient. However, different systems can be compared more easily by in silico experiments using data generated by the same numerical phantom. (4) As computing and storage capacities increase, the difficulties of performing accurate simulations reduce.
For the above reasons, researchers have conducted several numerical studies [21][22][23][24][25][26][27][28] to simulate the PA imaging process of human breasts, yet these models have limitations that remain unaddressed. Current PA breast phantoms are based either on segmentation of MRI/mammography dataset [21][22][23][24][25] or on the combination of a few simple objects [26][27][28]. Details of these studies are provided in Table 1. Compared with the existing PA breast phantoms, we aim to cover their advantages and bring several improved features: (1) Realistic phantoms for human breast and easy-to-use package tailored to simulate realistic PACT systems, which enables users to freely modify parameters. To meet these goals, we applied a newly reported software, known as the VICTRE breast phantom [29]. This software generates random voxelized digital breast phantoms step by step that contain skin, nipple, lactiferous duct, terminal duct lobular unit (TDLU), interlobular gland tissue, fat, suspensory ligament, muscle, artery, and vein. The general characteristics of the breast (volume, gland fraction, etc.) are selected by the user or randomly sampled from a population-based distribution [30]. A model for breast masses has been introduced as a supplement of the digital breasts [31], this mass model, however, was oversimplified for PACT. Following the same modeling strategy, we make further improvements to the mass model based upon anatomical and physiological properties of the breast. Breast cancers are classified as non-invasive breast cancer and invasive breast cancer (IBC). Non-Invasive breast cancers are those which do not spread, derived mostly from the mammary ductal epithelium [32]. Thus our breast mass models are classified into three major types: (1) benign fibroadenoma, (2) ductal carcinoma in situ (DCIS), and (3) IBC. Utilizing the new phantom, we simulated PA wave propagation and detection in two dimensions (2D) and three dimensions (3D). Then we compared PA images of different types of healthy breasts. In addition, we performed studies using our breast tumor models, and presented overlaid PA and ultrasound (US) images containing both anatomical and functional information. PA/US dual-modality imaging provides complementary information from different contrast mechanisms, which enhances each individual contrast for improved clinical performance. Figure 1(a) presents an overview of the VICTRE breast model. To account for different patient positions (i.e., standing, supine, or prone), we adjusted the simulation parameters to change the profile of the breast, as shown in Fig. 1(c) where the left profile corresponds to the standing position while the right one mimics the supine position. The model also allows for the control of breast density, which is known to affect lesion detection in digital mammography (DM) [5,6]. The breasts of the female population fall into four density categories: extremely dense (0.548 glandular volume fraction (GVF)), heterogeneously dense (0.339 GVF), scattered fibroglandular densities (0.143 GVF), and almost entirely fat (0.071 GVF). According to this, we created four types of breast models corresponding to different breast densities (shown in Fig. 1(d)) [33]. Once the model parameters are fixed, breast tissue textures are randomly generated with pre-determined statistical values. In this way, multiple breast models of each type can be made (with different microscopic tissue realizations) by simply changing a random number seed [30].

Methods
To set up a breast model, a series of deformations are performed to create a breast with a particular volume and shape, starting from a base quadratic surface. Then the Voronoi technique [34] is applied to generate the initial glandular compartments and grow a ductal tree-like structure with terminal duct lobular units from the nipple into each compartment. The next step involves using the Perlin noise function [35] to create a realistic gland/fat tissue interface and the Cooper's ligament network. Then, to simulate a realistic vascular network, random trees (one for each major vein and artery in human breast) originating from the backing muscle layer of the phantom are generated. A cost function is used to yield a preferential growth direction of the main branch in the skin surface/chest wall direction. These main branches spawn smaller secondary branches and the process proceeds recursively to fill the breast volume. In general, the generated vascular network tends to have bigger vessels near the skin surface and the chest wall, while smaller vessels are in between ( Fig. 1(b)). After finalizing the breast model, the stochastic Gaussian random sphere model [36] is used to generate a tumor mass in the center, and the iterative fractal branch algorithm is performed to add complex needle-like structures simulating microvascular clusters around the breast mass [31]. Studies have identified distinct patterns in PA images, corresponding to various vasculatures associated with different types of tumors [37][38][39]. For example, a study [39] reported that the centripetal blood vessel structure was observed in the PA images of breasts diagnosed with IBC, whereas intratumoral spotty signals were observed in the DCIS cases. In another study [37], the PA images showed absence of hemoglobin within or around fibroadenoma. Moreover, using diffuse optical spectroscopy and tomography, researchers have found that malignancies are characterized by a substantially higher concentration of blood, with significantly lower oxygenation than normal or benign tissue [40][41][42]. We modified the mass model taking into account the above observations, and produced three types of mass models: fibroadenoma, DCIS, and IBC (Fig. 2). In our phantom, fibroadenomas are modeled by a regular oval-shaped mass without vessels inside or microvasculature around. Meanwhile, IBC's are irregularly shaped and have a centripetal blood vessel structure as well as surrounding microvasculatures. For DCIS, we add a few thin vessels inside an irregular mass, to simulate the observed intratumoral spotty patterns. We also decrease the oxygenation level within malignant masses. The digital mass is embedded into the breast phantom in the end. TDLUs are the primary site for cancer formation. As a result, we choose the locations for the insertion of masses where a TDLU has been generated, including the chest wall. In our study, MCmatlab was used as the simulation platform. The outer boundary for the simulation zone is a rectangular cuboid whose internal volume is uniformly divided into cubic voxels. The user assigns to each voxel a medium or tissue type, described by its absorption coefficient µ a , its scattering coefficient µ s , and its Henyey-Greenstein scattering anisotropy factor g at the applied optical wavelength. Then an input light beam is simulated by launching photon packets and calculating their energy deposition in the simulated volume using a Monte-Carlo model (MC model) [43].
Near-infrared (NIR) light is preferred for PACT due to the relatively low optical scattering and absorption [44]. Hence we conducted our simulation at 700 nm and 900 nm using the optical parameters listed in Table 2 [45,46]. The optical properties of lactiferous duct, TDLU, and suspensory ligament are approximately equal, we define them all as the fiber tissue [47]. As it is computationally challenging to assign individual blood vessel elements with different optical properties [44], we assume average blood oxygenation (sO 2 ) of 75% for vein and 95% for artery. In this study, we simulated PA imaging using linear transducer arrays. Linear arrays have been extensively used in PACT for breast imaging due to their commercial availability and compatibility with traditional ultrasound scanners [44]. Figure 3 illustrates the schematic of the imaging setup in our simulations. A 128-element linear array (2.5 MHz central frequency, -6 dB bandwidth of 70%) is placed on top of the breast of a patient in the supine position. A linear illumination pattern is chosen because it can provide relatively homogenous optical fluence at the imaging plane [44]. Light sources are placed outside the model and provide a uniform fluence at the breast surface with an area of 12 × 12π mm 2 . The light incident position and the center line of the transducer imaging plane are in coincidence. The background medium is non-scattering and non-absorbing, and has no refractive index mismatch with the breast model. We employed the k-Wave toolbox [48] for all the simulations regarding PA wave propagation and detection. The initial PA pressure was simulated based on the optical fluence distribution from the 3D MC model. The initial pressure gave rise to PA waves that propagated in tissues which were acoustically heterogeneous (Table 3) [49][50][51] with negligible acoustic attenuation. The PA signals were detected by transducer elements that were curved in one dimension to create elevational focusing. This allowed us to simulate the acoustic sectioning of a real handheld probe. The detected signals were then band-pass filtered to account for the finite bandwidth imposed by the ultrasound detector and its electronics. To form an image, we used the simplest delay-and-sum (DAS) reconstruction algorithm. Since the optical absorption spectra of oxygenated hemoglobin (HbO 2 ) and deoxygenated hemoglobin (Hb) are different, sO 2 can be determined from PA images acquired at multiple wavelengths. In the simplest case when the spectral coloring effect is negligible, and Hb and HbO 2 are the dominant absorbers, total concentration of hemoglobin (tHb) and sO 2 can be estimated by acquiring images at two wavelengths, followed by linear spectral unmixing [52]. To maximize the estimation accuracy, the two selected wavelengths should be on the opposite side of the isosbestic point [37]. Wavelength-dependent light scattering has a substantial influence on the optical fluence maps. For the breast phantom shown in Fig. 4(a), the light fluence distribution at 700 nm and 900 nm are shown in Fig. 4(b), the distribution of the optical fluence has become quite different after propagating just a few millimeters (Fig. 4(c)). The artifacts in the PA images are mainly caused by the limited detection angular coverage, the heterogeneous acoustic properties, and insufficient spatial sampling frequency. The existence of the artifacts increases the difficulty of sO 2 estimation, and thus we employ the region growth method for image segmentation and accordingly apply a mask to selectively calculate sO 2 in the area of interest.
To display the reconstructed PA images on top of gray-scale ultrasound images (as many researchers do in their experiments), we also simulate pseudo-B-mode ultrasound images using the same numerical phantom. These images are termed 'pseudo' because they are not created by rigorous simulations; rather, they are directly rendered from the echogenicity feature maps of the numerical phantoms. The image is formed by first adding motion blur to the input acoustic impedance map and then multiplying it with a speckled noise pattern. Such speckle patterns are generated from our in vivo US images by the exemplar-based texture synthesis method [53].

Results
All our MC simulations of light fluence distribution were performed in 3D, while for acoustic propagation and detection we did both 2D and 3D simulations. 2D simulations are much faster, while 3D simulations are in principle more accurate. PA images corresponding to the three types of breast mass models (shown in Fig. 2) are given in Fig. 5, from which we conclude that the 2D-and 3D-models are quite similar in terms of revealing anatomical structures, and both suffer from artifacts due to the limited-view detection [54,55]. For the 2D model, the reconstructed PA image can better recover the actual anatomical structures, because out-of-plane artifacts are absent. In the following simulations, we choose to use the 2D model for higher efficiency. An extra simulation was carried out to investigate if various breast density affects the PA image features, especially imaging depth. Four breast models corresponding to different breast densities as displayed in Fig. 1(d) were used in the simulation. Considering the optical fluence distributions could be remarkably affected by the vascular distribution, we applied the same vascular structure. MC simulation was performed at a wavelength of 900 nm and light intensity on the breast surface was about 150 mw cm −2 . After that we employed a 128-element linear ultrasonic transducer array (0.75 mm pitch, 0.25 mm inter-element spacing) for signal detection. The received PA signals were reconstructed by DAS algorithm assuming a uniform speed of sound of 1500 m/s. Figure 6 shows that breast density does not have a significant influence on the optical fluence distribution. The penetration depths were all around 2 cm for the four types of breasts with dramatically different density, according to the MC simulation. As expected, in the right column of Fig. 6, PA features reconstructed at the same depth were similar for all four cases. This observation agrees with experiments, that is, compared to mammography, PA image quality is less sensitive to radiographic density of the breast [12].
Using our simulation platform, in silico imaging can be easily performed and compared with experimental ones. As an example, we simulated PA/US dual-modality imaging of a dense breast with an embedded lesion whose center was located 1 cm below the skin. To measure sO 2 , we acquired PA images at both 700 nm and 900 nm and then employed linear spectral unmixing for demodulation. However, linear unmixing only works when the excitation spectrum stays unchanged at different depths [56]. As we have demonstrated in Fig. 4(c), the excitation spectrum is affected by wavelength-dependent light scattering and absorption. This could induce substantial error in the sO 2 measurement. Thus for deep-seated targets optical fluence compensation is essential for accurate sO 2 imaging [56,57]. Here we simply assumed that perfect compensation has been made to create the gold standard sO 2 maps. By applying the mask segmented from the reconstructed PA images before calculating, we can eliminate most artifacts without loss of the features. Figure 7 illustrates three images corresponding to three types of tumors with their respective PA features, the anatomical structure of the digital phantom has been shown in Fig. 5(a). The embedded mass had a diameter of 8 mm, representative of an early-stage lesion, which can be clearly seen from the US images as hypoechoic masses. These images exhibit the common features of PA/US dual-mode imaging where US provides the grayscale structural background while superimposed PA images indicate functional information, such as tHb and sO 2 [37]. The figures show the possible advantage of PA/US imaging over US imaging alone. For example, in Fig. 7(c), the US image feature may suggest a benign lesion such as a fibroadenoma. However, the centripetal vessel with relatively low oxygenation, as well as increased blood supply around the mass, indicate possible malignancy. Figure 7(b) provides another example of how PA imaging may be useful for differentiating among subtypes of malignant mass.

Discussion
Malignant tumors are often associated with angiogenesis and reduced sO 2 due to their aggressive growth [37]. Therefore, PA imaging is a potential candidate for breast cancer diagnosis and screening [12] thanks to its high sensitivity to morphological and functional changes of the vascular network. The development and clinical translation of the PACT technology can be accelerated by well-designed digital phantoms which have realistic physiological, anatomical and physical properties. Here, we introduce a breast phantom created from the computational methods. The phantom is not only realistic, but also flexible enough to allow adjustments of its density, shape, and detailed structural realization. We performed 3D MC simulation and used k-Wave to simulate acoustic propagation and detection. In our demonstrations, we chose to apply acoustic simulation in 2D for better efficiency without losing much accuracy. Based on the digital phantom, we also generated grayscale US images upon which pseudo-color PA images were superimposed.
An advantage of this computational model is that it is designed specifically for photoacoustic imaging. The simulated tissue types, anatomical structures, and optical/acoustic properties are all chosen to be realistic. Based on the breast phantom, the simulated PA/US images have shown features that are commonly found in experimental images. This model also offers the capability of creating a large population of virtual patients having benign or malignant tumors, and with dramatically diverse tissue contents (density) and breast shapes. Used in conjunction with the optical and acoustic simulation software reported above, this phantom allows researchers to investigate a number of technical problems conveniently (since no real experiments are required) and quantitatively (since the gold standard is available). For example, different light illumination and detection schemes can be quantitatively compared and optimized, if the simulated results could emulate the real images. Similarly, the reflection artifacts and spectral coloring effect in PA imaging can also be analyzed. By comparing PA images generated from breasts that are acoustically heterogeneous (according to Table 3) and homogenous (speed of sound 1500 m/s, density 1000 kg/m 3 ), the generation of reflection artifacts was clearly observed in Fig. 8(a-b). Likewise, we calculated the sO 2 map around the mass (IBC) by linear unmixing, but without fluence compensation. In contrast with Fig. 7(c), the inaccurate estimation of sO 2 (Fig. 8(c)) clearly reveals the spectral coloring effect [37,56]. Comparing our simulation results with the in vivo ones (obtained using a custom-made hand-held PA imaging probe, central frequency 2.5 MHz, bandwidth 55%, 128 channels), the synthetic images show highly similar features (Fig. 9). The digital phantom can be further improved iteratively using clinical data as feedback. In practice, due to limitations in bandwidth, spatial sampling density, and viewing angle, important vascular features may be distorted or even absent in the reconstructed image; due to acoustic reflections and reverberations, artifacts can be generated and superimposed on true structures; due to the spectral coloring effect, using two colors for sO 2 demodulation can be highly unreliable. Optimization of the imaging hardware can be systematically pursued in silico using the digital phantoms, taking into account all the above-mentioned complexities. Now, open source dataset for PA breast imaging is insufficient, while a large amount of data can be generated digitally to train neural networks for various tasks such as image reconstruction, vessel segmentation, lesion detection, artifacts removal, quantitative imaging, etc.. Another advantage of using digital phantoms for algorithm development is that gold standards are readily available.
The mass models used in this study were constructed based on a few cases of photoacoustic imaging results. These models are still quite limited to accurately represent the structural and functional diversity of real tumors and tumor microenvironments. The optical and acoustic properties and the exact 3D shapes of the solid tumors, as well as the morphology and physiological state of the vascular networks inside and surrounding the tumor due to angiogenesis, all need to be more accurately modelled in the future.

Conclusion
In this paper, we present a realistic, 3D digital breast phantom which exhibits the fine structures observed in in-vivo PA/US images. We find that although 3D acoustic model is more precise, 2D model is more efficient with sufficient accuracy. The numerical phantom is useful for researchers to develop PA systems for breast cancer imaging, and to investigate problems such as reflection artifacts and the spectral coloring effect. On the one hand, using a phantom with the real complexity of the human breast, it is possible to study the spectral coloring effect in detail. On the other hand, 3D modelling of photon transportation in digital breasts could provide information for more sophisticated optical fluence compensation. Using our phantom in conjunction with prior structural information (e.g., acquired using US imaging), and by properly adjusting tissue optical properties, we could customize a personalized compensation scheme for each patient. This is another potential application where our phantom can contribute to better functional PA imaging. Owing to insufficient number of cases reported of PA imaging of the breast, the mass model for different types of breast cancer is still rudimentary [38]. Our future research will focus on improving the mass model to better reflect the diversity and complexity of real tumors. Disclosures. The authors declare no conflicts of interest.