1.1.

Photoacoustic Tomography and Photoacoustic Tomography Phantoms

Breast cancer is the second leading cause of cancer-related deaths in American women.¹ X-ray mammography, the main screening method, requires harmful radiation and has low sensitivity for cancer detection in dense breasts.² Ultrasound can improve breast cancer detection in dense breasts³ and also distinguish benign from cancerous masses, but results in unnecessary biopsies due to low specificity.⁴ A possible remedy for this problem, the emerging imaging modality photoacoustic tomography (PAT), holds strong potential for mammography applications.⁵ PAT combines optical tissue excitation with acoustic detection to enable imaging of light-absorbing chromophores (melanin, hemoglobin, water, and lipids). Short-duration light pulses are delivered to tissue, where light absorption results in rapid thermoelastic expansion and subsequent emission of ultrasonic pressure waves. These waves can be detected at the tissue surface using acoustic transducers, such as clinical ultrasound probes.⁶ Beamforming algorithms are then used to backproject detected acoustic waves to their sources (light-absorbing tissue regions). Because sound is attenuated much more weakly than light in tissue, PAT can produce images of optical absorption in deep tissue (up to several centimeters).⁷ The primary endogenous source of PAT image contrast is hemoglobin found in blood, and several groups have proposed using multiwavelength PAT imaging to generate maps of oxygen saturation.^8,9 These images could be used to identify the abnormal vascular structure seen in cancers, resulting in higher specificity for cancer diagnosis than ultrasound alone. While there are currently no clinical PAT systems approved or cleared for use in the United States, several commercial systems, e.g., Seno Medical (San Antonio, Texas),¹⁰ iThera Acuity (München, Germany)¹¹, and Canon (Tokyo, Japan),¹² and research systems, e.g., University of Twente (Enschede, The Netherlands)¹³ and Washington University (St. Louis, Missouri),¹⁴ have been applied in patients. However, due to the nascent status of PAT, there are no well-established benchtop methods for predicting in vivo PAT system performance. The development of such approaches is essential for improving the understanding of complex tissue interactions and their impact on real-world clinical imaging system performance.

Consensus standards and best-practice documents for mature medical imaging modalities (e.g., mammography, magnetic resonance imaging, and ultrasound) identify tissue-simulating phantoms as essential tools for preclinical performance testing.^15,16 Standardized phantoms often implement an idealized geometry to enable uniform evaluation of the effect of one type of parameter on PAT image quality while keeping other parameters the same. A wide variety of tissue-mimicking materials (TMMs) have been studied for use as PAT phantoms, including traditional gelatin hydrogel,^17,18 poly(vinyl alcohol) cryogel,^19,20 and poly(vinyl chloride) plastisol (PVCP).^21–24 Most PAT phantoms presented in the literature use a homogeneous background medium to simulate healthy tissue. However, breast tissue is a highly heterogeneous mixture of fat and fibroglandular tissue. The fibroglandular tissue is surrounded by, admixed with, and covered by a layer of fatty tissue.²⁵ These two types of tissue have been reported to have substantially different optical and acoustic properties, and spatial variations of these properties in the breast may cause errors in image reconstruction or changes in image quality. In ultrasound imaging, the heterogeneity in acoustic properties between the fat and fibroglandular tissue results in significant degradation of ultrasound image quality and spatial resolution. Using acoustic detection, PAT also needs to evaluate this effect of acoustic heterogeneity on PAT image quality. Homogeneous phantoms with uniform optical and acoustic properties may be suitable for basic performance testing, but it is likely that heterogeneous phantoms are needed to simulate complex tissue interactions that affect in vivo imaging performance. To realize such phantoms, the spatial distribution of TMMs and their optical and acoustic properties need to mimic those reported for complex breast tissue.

1.2.

Optical and Acoustic Properties of Breast Tissue

Defining appropriate representative optical and acoustic properties of breast tissue is a challenging problem. The main concern is that both optical and acoustic properties vary significantly due to high variations in intrinsic breast tissue composition. Breast density, or the relative content of fatty tissue and fibroglandular tissue, is a major driver of this variation and may depend on many factors, including body mass index, age, and menopausal status.²⁶ In the field of mammography, breast density is used as a coarse estimate of the amount of fibroglandular tissue present in a breast. The Breast Imaging and Reporting Data System (BI-RADS)²⁷ is a well-accepted standardized format for reporting mammography cases. BI-RADS includes a categorization structure for mammographic density with four types: a (almost entirely fat), b (scattered fibroglandular densities), c (heterogeneously dense), and d (extremely dense). The wide variations in subject breast density have thus led to a broad variation in breast property measurements reported in the literature.^28–32

The most important optical properties governing PAT system performance are the optical absorption and scattering coefficients, which have been extensively studied in breasts for the past three decades. There is a significant volume of literature available on the optical absorption^{26,28,33–35} and scattering properties^26,34 of in vivo breast tissue measured using diffuse reflectance spectroscopy and diffuse optical tomography. One important limitation is that these measurements are averaged over large tissue volumes, and it is difficult to isolate the optical properties of pure fatty tissue from pure fibroglandular tissue. Therefore, the majority of optical properties reported are not well categorized according to BI-RADS breast tissue types. Some reported values for pure fatty or fibroglandular breast tissue are available from ex vivo measurements, but these may not be completely representative of in vivo tissue due to effects such as dehydration or exsanguination. Several studies have reported variations of in vivo breast optical properties with breast density (relative fat/fibroglandular tissue composition). However, these data show high variability due to other factors such as age and menopausal state, which may contribute to changes in breast density and relative water and lipid content.^26,28 Optical property ranges of BI-RADS breast tissue types at 750 nm are listed in Table 1. The values reported using the combination of the middle two types are labeled as b + c. Optical properties for additional wavelengths will be presented in Sec. 3.1.

Table 1

Reported optical and acoustic properties as functions of BI-RADS breast tissue type.

The critical acoustic properties affecting PAT system performance are speed of sound, frequency-dependent acoustic attenuation, and backscattering coefficient. Speed of sound describes the propagation velocity of acoustic pressure waves in the medium and is a key input parameter to image reconstruction algorithms, such as a delay-and-sum beamformer, where the received ultrasound signals are delayed to account for time-of-flight differences across the acoustic detector surface and summed at each reconstruction point. Both in vitro measurements^31,39 on excised breast tissues maintained at 37°C and in vivo measurements^32,36,40,41 reported sound speeds of less than $1460\mathrm{m}/\mathrm{s}$ in fat and sound speeds more than $1500\mathrm{m}/\mathrm{s}$ in fibroglandular tissue over the frequency range of 2 to 13 MHz. Speed of sound values of BI-RADS breast tissue types at 8 MHz are listed in Table 1 because this frequency is close to the central frequency, 7.5 MHz, of our imaging probe and there is one reference at this frequency. Note: the speeds of sound reported³⁶ at 8 MHz in Table 1 are lower compared to the corresponding values at the other frequencies in the figures in Sec. 3.2. However, the trend in speed of sound versus breast tissue types is similar. The spatial variation in breast speed of sound due to the heterogeneous distribution of fatty and fibroglandular tissue can degrade the spatial resolution and accuracy of spatial registration of ultrasound and photoacoustic images when a single uniform speed of sound is assumed in the image reconstruction.^42,43 Breast acoustic attenuation has been extensively characterized by several groups, particularly as a function of breast composition. Ex vivo^31,44–46 and in vivo measurements³⁸ have demonstrated lower acoustic attenuation in fatty tissue and higher attenuation in fibroglandular tissue over the frequency range of 1 to 10 MHz. Attenuation of BI-RADS a and b + c breast tissue types at 8 MHz is shown in Table 1. As with optical properties, both the speed of sound and acoustic attenuation of the whole breast also vary with menopausal status and age.³⁰

Because many PAT systems presented in the literature are bimodal,^47–49 allowing for overlay of coregistered PAT and B-mode ultrasound images, complex PAT tissue phantoms should be suitable for imaging in both domains. Thus, it is important to account for the acoustic backscattering coefficient, which describes the relative amount of energy reflected back toward an acoustic emitter. This mechanism represents the primary source of contrast in B-mode ultrasound images. It is well known that breast fatty and fibroglandular tissues have substantially different backscattering properties, as fibroglandular tissues produce high intensity in B-mode ultrasound images while fatty tissue is relatively hypoechoic. Ex vivo and in vivo studies have demonstrated that breast fibroglandular tissue may possess an order of magnitude higher backscattering coefficient compared to fatty tissue.^37,38,44 Backscattering of BI-RADS a and d breast tissue types at 8 MHz is shown in Table 1.

1.3.

Study Objectives

Stated succinctly, the overall objective of this research was to advance clinical translation of mammographic PAT through the development of a realistic tissue-simulating phantom and its implementation to elucidate practical issues impacting PAT image quality in the breast. The specific goals of this work included: (1) developing custom TMMs incorporating realistic optical and acoustic properties of breast fat and mixed fatty-fibroglandular tissue, (2) fabricating two-layer PVCP phantoms composed of these TMMs arranged in a tissue-relevant layered geometry with an undulating boundary, and (3) using this innovative breast phantom to assess the effect of tissue properties and reconstruction algorithm parameters on PAT image quality. By achieving these goals, we can take a major step toward filling critical gaps in the literature and facilitate standardization of this innovative technology.

2.1.

Phantom Fabrication

We have previously developed custom PVCP formulations for use in making PAT phantoms with tunable, biologically relevant optical and acoustic properties.²⁴ PVCP is a suspension of PVC resin in liquid plasticizers, which fuses into a gel after heating to high temperatures ($\sim 180°\mathrm{C}$). In this study, we utilize the same base materials, including a dispersion-grade PVC resin (Geon 121A, Mexichem Specialty Resins, Inc., Avon Lake, Ohio), benzyl butyl phthalate (BBP) (Sigma-Aldrich, St Louis, Missouri), and di(2-ethylhexyl) adipate (DEHA) (Sigma-Aldrich). For the first aim of the study, we were able to develop two TMM formulations simulating BI-RADS a and c type breast tissues, respectively. The acoustic properties of the two TMMs were kept within desired acoustic property ranges over frequencies from 4 to 9 MHz of the corresponding tissue types, while their optical properties were kept similarly within desired optical properties of BI-RADS b + c type breast tissue over a wavelength range of 650 to 850 nm. These two types of TMM were labeled as a-fat and c-fibrogland.+fat, respectively (Table 2). The speed of sound and acoustic attenuation were tuned by modifying the volume ratio between the two plasticizers (BBP and DEHA) and the mass/mass PVC concentration. Acoustic backscattering was adjusted by adding soda lime glass microspheres (Spheriglass A, Potter Industries LLC, Malvern, Pennsylvania). Microspheres with diameters of 38 to $60\mu \mathrm{m}$ and 63 to $75\mu \mathrm{m}$ were used in the formulations for a-fat and c-fibrogland.+fat tissue, respectively. Optical scattering was induced by adding anatase titanium dioxide (Sigma-Aldrich). Once the optical and acoustic properties fell within the desired range of the tissue type, the tuning process was stopped. There are two limitations to using the combination of BBP, DEHA, and PVC. One is that the current TMM is not able to achieve a speed of sound higher than $\sim 1510\mathrm{m}/\mathrm{s}$. However, the two TMM formulations utilized in this study are suitable because their acoustic properties are biologically relevant and present sufficient differences in acoustic properties to affect image quality. The other is that optical absorption of the current TMM cannot be further lowered to match that of fatty tissue. Since we previously showed this TMM to be optically tunable, we could have imparted distinct sets of optical properties to each layer, mimicking different breast tissue composition types. However, to specifically investigate the effect of acoustical heterogeneity on photoacoustic image quality, we tuned the two formulas to have similar, but still biologically relevant, optical properties.

Table 2

TMM formulations for stimulating breast fat (a-fat) and fatty-fibroglandular (c-fibrogland.+fat) tissue.

PVCP phantoms were fabricated as described previously.²⁴ Briefly, a large stock PVCP solution was prepared by mixing a binary mixture of plasticizers with 1% v/v calcium-zinc heat stabilizer (M-F Manufacturing Co., Fort Worth, Texas), adding PVC resin, then stirring for 30 min followed by degassing for 60 min. A 75-mL volume of PVCP solution was mixed with ${\mathrm{TiO}}_2$, sonicated at 40°C for 40 min, and then reintroduced to the original PVCP stock solution and stirred for 5 min. Glass microspheres were added to a 100-mL batch of PVCP solution and stirred for 5 min. This batch volume was then poured into a 200-mL round-bottom flask immersed in a magnetically stirred oil bath maintained at 190°C using a thermocouple. The flask, which contained a stir bar, was evacuated and stirred at $\sim 375\mathrm{rpm}$ for 15 min. The flask was then removed from the oil bath and stirred magnetically at a gradually reducing speed until the flask temperature was $\sim 110°\mathrm{C}$. After reaching this temperature, PVCP was poured into aluminum molds. This whole cooling protocol was used to increase PVCP viscosity at time of pouring, which would minimize settling of the glass microspheres. Samples for phantom material characterization were first produced, including 5-mm-thick, 38-mm-diameter disks for speed of sound, optical property measurements (Sec. 2.2), and acoustic attenuation (Sec. 2.3). A 20-mm-thick, 38-mm-diameter disk was produced for backscatter coefficient measurements (Sec. 2.3).

2.2.

Optical Characterization

As described in Ref. 24, PVCP disks were placed between 1-mm-thick glass slides (refractive index = 1.51), and diffuse transmittance and reflectance measurements were made over 600 to 1000 nm using an integrating sphere spectrophotometer (Lambda 1050, PerkinElmer, Waltham, Massachusetts). NIST-traceable Spectralon standards were used to normalize measurements. Optical absorption coefficients (${\mu }_{\mathrm{a}}$) and reduced scattering coefficients (${\mu }_{\mathrm{s}}^{\prime }$) were calculated using the inverse adding-doubling method,⁵⁰ which requires a priori knowledge of the anisotropy factor and refractive index of the sample. Using the Lorentz–Lorenz mixture rule for refractive index,⁵¹ a-fat and c-fibrogland.+fat TMMs are expected to have a refractive index of 1.485 and 1.528, respectively. The scattering anisotropy factor was assumed to equal 0.7 based on the Mie scattering theory of ${\mathrm{TiO}}_2$ (using open-source MATLAB^® software⁵²) for both TMMs.

2.3.

Acoustic Characterization

The sound speed of PVCP-based TMM was measured on 5-mm-thick, 38-mm-diameter disks^1,24,53 based on through-transmission⁵³ using a pair of identical broadband transducers (V320, Panametrics, Waltham, Massachusetts) with focal length of 3.81 cm, diameter of 1.27 cm, and center frequency of 7.5 MHz. The transducers faced each other in a water bath, with one transmitting, the other receiving, and the PVCP disk under test positioned at their focus. They were connected to a pulser/receiver (Model 5800PR, Panametrics), and the received ultrasound signals were digitized (8 bit, 50 MHz) using an oscilloscope (9310C, Teledyne LeCroy, Chestnut Ridge, New York). The speed of sound in a sample disk, $c_s$, was calculated as⁵³

Achieved concurrently with the sound speed measurement, the frequency-dependent acoustic attenuation coefficient was derived as¹

Backscatter coefficients were measured using the same setup, except in reflection mode involving only one transducer.^54–56 The central plane of the 20-mm-thick TMM disk under test was aligned with the focal plane of the transducer. Pulse-echo signals were recorded at 30 spatial positions across the face of the disk by laterally translating the TMM disk at 0.5-mm steps. At each position, 20 acquisitions were averaged. The signals were time gated to isolate echoes of 1.5 mm in depth centered at the focal plane, and the mean, $\overline{P_s(f)}$, of the corresponding power spectra was calculated. A reference spectrum, $\overline{P_r(f)}$, was obtained using the same approach but replacing the TMM disk with a planar reflector made of low-density polyethylene at the transducer focus. This type of material was chosen to avoid signal saturation during the measurements. The normalized power spectrum, $P_o(f)$, was calculated as⁵⁶

2.4.

Homogeneous and Heterogeneous Phantoms

A phantom ($86\times 60\times 35\mathrm{mm}$, 180.6 mL) with homogeneous background (no undulating boundary) was constructed by pouring PVCP (c-fibrogland.+fat formulation) into a mold containing six 0.5-mm-diameter steel wires (A-M Systems, Sequim, Washington) spaced 5 mm apart vertically and horizontally, as shown in the left side of Fig. 1(a). Metallic wires were chosen because they are visible in both ultrasound and PAT images and easily aligned due to their high tensile strength. The diameter of the wire is similar to the size of a brachytherapy seed (metallic implant used for lesion identification and treatment) and in the diameter range of a blood vessel.^57,58 Since the phantom was larger than our 100-mL batch volume, the phantom was produced by sequential layer pouring. Prior to pouring the second layer, a heat gun was applied to the surface of the first layer to produce a thin layer of remelted PVCP. As the second layer was poured, this remelting prevented formation of air gaps and poor PVCP bonding, which would cause photoacoustic artifacts at the boundary. Both layers used the c-fibrogland.+fat formulation to mimic dense breast. A heterogeneous phantom with the same layout was created using the same mold, but pouring the a-fat PVCP formulation for the first layer and the c-fibrogland.+fat formulation for the second layer. This phantom also had an undulating boundary between its two layers to mimic the boundary between superficial fatty tissue and deeper fibroglandular tissue in breast. An in vivo ultrasound B-mode image is shown in Fig. 1(b) as an illustration of this tissue structure. To create this feature, an aluminum block with machined concave surfaces [Fig. 1(d)] was inserted into the mold [Fig. 1(c)]. PVCP solution for the first (fatty) layer was poured onto the plate, and after cooling, the plate was removed, resulting in convexities at the surface of the layer onto which the deep (fibroglandular) layer of TMM was poured. The boundary layer consisted of four segments whose radii of curvature were designed to reflect curvatures found by analysis of the fat-glandular tissue boundary in the ultrasound B-mode images of 10 patients acquired for a different research study (an IRB-approved ultrasound elasticity study conducted at George Washington University). All the images included in this analysis present a clear boundary that can be approximately fitted with an arc. Based on the fitting results, we then modeled the boundary between superficial and deep breast layers as a cascade of convexities covering a reasonable range of radii of curvature. In Fig. 1(a), the boundary convexities (from left to right) have radii of curvature of 30, 24, 16, and 16 mm and boundary angles of 58 deg, 66 deg, 81 deg, and 41 deg.

Fig. 1

The cross-sections of two heterogeneous phantoms: phantom 1 with wires on the left side of (a) and phantom 2 with tubes on the right side of (a) with boundary angles theta of 58 deg, 66 deg, 81 deg, and 41 deg (from left to right), (b) an example of ultrasound breast image with an undulating boundary between fatty and fibroglandular tissue containing a cancer (circled), and an example of (c) an aluminum mold with (d) an undulating plate.

Imaging was performed using a previously developed custom PAT system composed of a cart-based tunable near-infrared pulsed laser (Phocus Mobile, Opotek, Inc., Carlsbad, California) and a research-grade ultrasound system (Vantage 128, Verasonics, Inc., Kirkland, Washington). Phantoms were imaged at 750 nm and a radiant exposure of $20\mathrm{mJ}/{\mathrm{cm}}^2$. Acoustic sensing was performed using a 128-channel ultrasound linear array transducer with a 7.5-MHz center frequency, 7.0-MHz bandwidth, and 38.1-mm length (L11-4v, Verasonics). A significant clutter artifact encountered using this approach was caused by photoacoustic generation at the transducer face resulting from high fluence at the phantom surface and optical absorption of transducer surface material. To reduce this effect, the transducer surface was covered with aluminum foil and coupled to the transducer surface with a thin layer of acoustic coupling gel. While the foil can still cause clutter artifacts, the overall image quality is significantly improved. The transducer was then brought into contact with the phantom top surface using a thin layer of water as the couplant. One-way delay-sum beamforming was applied to the received radio-frequency ultrasound signals to create photoacoustic images with dynamic receiving focusing. Data were apodized using Hanning windows, and the $F$ number was set at 0.9. The reconstruction speed of sound was varied between 1437 and $1500\mathrm{m}/\mathrm{s}$ in the heterogeneous phantom and between 1440 and $1560\mathrm{m}/\mathrm{s}$ in the homogeneous phantom. The spatial size of the metal wires with the same diameter at different depths across the phantom was estimated and compared in the reconstructed PAT imaging domain. The axial and lateral spatial sizes of these imaged targets were defined as the full-width half-maximum (FWHM). All PAT images generated in this study were normalized using the same maximum intensity value. B-mode ultrasound images were also acquired using the same PAT system, applying two-way delay-sum beamforming to received data from plane wave transmission. Gray levels for all ultrasound images were normalized to a single maximum intensity value.

Target signal-to-noise ratio (SNR) and contrast-to-noise ratio (CNR) were also computed using uncompressed image intensity data to study the effects of the choice of reconstruction speed of sound on these two imaging quality parameters. To calculate these parameters, a $2\times 2\mathrm{mm}$ rectangular region of interest (ROI) was first centered at each target. A masked subset of this ROI using an FWHM ($-6\mathrm{dB}$) was selected as a target ROI. The background ROI was the same size and immediately to the right of the target ROI. Target SNR values were calculated as the mean of the target ROI divided by the standard deviation of the background ROI. Target CNR was calculated as the difference between the mean of the target ROI and background ROI over the standard deviation of the background ROI.

2.5.

Undulation Effects in Heterogeneous Phantoms

To evaluate the effects of the boundary undulation between superficial (fatty) and deeper (fibroglandular) layers on PAT image quality, a second complex phantom was constructed using the same process as the heterogeneous phantom described in Sec. 2.4. Instead of wire inclusions at different depths, a horizontal array of six polytetrafluoroethylene (PTFE) tubes (Component Supply Co., Fort Meade, Florida) was positioned at a depth of 20 mm, as shown in the right side of Fig. 1. These tubes have an inner diameter of 0.559 mm and a wall thickness of 0.150 mm and were filled with an India ink solution. We used a 1% aqueous solution of India ink (Speedball, Statesville, North Carolina) with an optical absorption of $4.6{\mathrm{cm}}^{-1}$ at 750 nm. This value was chosen to simulate venous blood with a hemoglobin concentration of $15\mathrm{g}/\mathrm{dL}$ and oxygen saturation of 70%.⁵⁹ Similarly, as with the target size measurements in Sec. 2.4, the spatial extent of these photoacoustic signals is also estimated using their FWHM in the corresponding direction. There are signals close to the top and bottom of the tubes. Those close to the top wall were used for spatial size and intensity estimation. Axial size was calculated by doubling the half width at half maximum of the rising part of the bell-shaped signal closest to the surface of the phantom since the falling portion of the bell-shaped signal was affected by the signal close to the bottom wall of the tubing.

3.1.

Optical Properties

Optical properties from the literature have been classified according to four BI-RADS categories. Optical properties of b and d breast tissue types are included to delineate the property range of the other two breast tissue types because the change of optical properties follows the change of breast density. The same rule is also applied to acoustic properties reported from the literature. Sandell and Zhu³⁵ summarized the optical properties of fatty breast tissue measured with four in vivo experimental methods. The resulting 95% confidence intervals are shown in Fig. 2. Jacques³⁴ focuses on the optical properties of mixed breast tissue, which can be represented by the fitted curves. The mean and standard deviation of these curves are also shown in Fig. 2 as shaded areas with the mean at the center. The range of values in Jacques’s review paper overlapped with the lower values in Sandell and Zhu’s review paper. The large variation might be due to the variation in a range of biological factors, such as breast density across subjects, and the measurement variation among experimental methods. Optical properties for four tissue types measured in vivo on human subjects are reported by Taroni et al.²⁶ The values shown in Fig. 2 are calculated from the supplementary material of Taroni et al.’s⁶⁰ paper provided online. The standard deviations were shown only for tissue type c for clarity. The optical properties of our two phantom formulations were designed to be similar to each other, within the optical property range of breast tissue, and also capable of elucidating the effect of heterogeneous acoustic properties on image quality. As shown in Fig. 2(a), the optical absorptions of two formulas are similar to each other except for the wavelength range 900 to 940 nm. These values are slightly higher than the results reported by Taroni et al. over 670- to 920-nm wavelengths, but still in the range of values reported by Sandell and Zhu. The reduced scattering coefficients of both phantom materials are similar and consistent with the values of BI-RADS c type breast tissue reported by Taroni et al. over the wavelength range of 600 to 1000 nm.

Fig. 2

(a) Optical absorption coefficient and (b) reduced scattering coefficient.

3.2.

Acoustic Properties

Speed-of-sound values reported in the literature^{31,32,36,39–41} are categorized into four BI-RADS tissue types (a, b, c, and d) according to the description of tissue types in the literature (fat, fatty-fibroglandular, less fatty-fibroglandular, and fibroglandular tissue in order), as shown in Fig. 3. These values are reported as either $\text{mean}\pm \text{standard deviation}$ or as a range of values. There are large variations in speed of sound for each tissue type. The general trend of increasing speed of sound is reported with increasing amounts of fibroglandular tissue. Our formulations are designed to achieve a speed of sound of $1437\mathrm{m}/\mathrm{s}$ for a-fat and $1500\mathrm{m}/\mathrm{s}$ for c-fibrogland.+fat. These values fall in the sound speed range of tissue types of a and c and are consistent with the references.

Fig. 3

Speed of sound measured for four types of breast tissue from references compared with those of the two TMMs formulated. s1: Kossoff et al.,³² s2: Carson et al.,⁴⁰ s3: Glover,⁴¹ s4: Greenleaf and Bahn,³⁶ s5: Bamber,³¹ s6: Foster and Hunt.⁴⁵

Acoustic attenuation values reported in the literature^31,38,44,45 are shown in Fig. 4. Shaded areas represent either the range of values or standard deviation around the plotted mean values. The trend of increasing sound speed with increasing glandular tissue content is also observed in acoustic attenuation values, with lower values in fat and higher ones in fatty-fibroglandular tissue. The acoustic attenuation of the phantom formulations is consistent with those in the references and is representative of tissue types a and c.

Fig. 4

Acoustic attenuation coefficient measured for four types of breast tissue from references compared with those of the two TMMs formulated. r1: Bamber,³¹ r2: Foster and Hunt,⁴⁵ r3: D’Astous and Foster,⁴⁴ r4: D’Astous and Foster,⁴⁴ r5: Nasief et al.³⁸ Acoustic attenuation of TMMs a-fat and c-fibrogland.+fat are $(0.082\pm 0.005)\times f^{2.14\pm 0.03}\mathrm{dB}/\mathrm{cm}$ and $(0.252\pm 0.008)\times f^{1.90\pm 0.01}\mathrm{dB}/\mathrm{cm}$, respectively, where $f$ is the frequency in MHz.

Acoustic backscatter coefficients reported in the literature^37,38,44 are shown in Fig. 5 for tissue types a and d. Tissue type d has significantly higher backscatter coefficients than type a. The a-fat formulation gives a backscatter coefficient consistent with that of type a breast tissue reported. The c-fibrogland.+fat formulation gives backscatter coefficients slightly lower than the reported values for d, which is consistent with the dependence of acoustic properties on the proportion of fibroglandular tissue present.

Fig. 5

Backscatter coefficients for two types of breast tissue from references compared with those of the two TMMs formulated. g1: D’Astous and Foster,⁴⁴ g2: Anderson et al.³⁷ using 7.5 MHz transducer; g3: Anderson et al.³⁷ using 10 MHz transducer; g4: Nasief et al.³⁸

3.3.

Homogeneous and Heterogeneous Phantoms

The TMM formulations we have developed enable the construction of homogenous and heterogeneous phantoms. Figures 6(a) and 6(b) are the reconstructed ultrasound and PAT images of a heterogeneous phantom with six wires using the speed of sound of the top layer, $1437\mathrm{m}/\mathrm{s}$. Figures 6(c) and 6(d) are the reconstructed ultrasound and PAT images of a homogeneous phantom with six wires using a speed of sound of $1500\mathrm{m}/\mathrm{s}$. The heterogeneous phantom can be used to evaluate the image quality of the dual-modality photoacoustic/ultrasound system. The ultrasound images shown in Fig. 6(a) provide valuable anatomical information mimicking the clinical imaging environment. This is due to a realistic difference in backscatter coefficient for the two layers. While we simulated one specific breast heterogeneity case with a specific set of layer properties, this TMM and phantom method could be readily tuned/adapted to simulate different tissue composition ratios and geometries.

Fig. 6

(a) Ultrasound and (b) PAT images of heterogeneous and (c) ultrasound and (d) PAT images of homogeneous phantoms containing six 0.5-mm-diameter wires.

3.4.

Effects of Reconstruction Speed of Sound on Photoacoustic Tomography Image Quality

The heterogeneous phantom described above provides a realistic imaging environment to evaluate the effects of boundary undulation and reconstruction speed of sound on PAT target size. Our conventional delay-and-sum reconstruction algorithm assumed a homogeneous speed of sound, and the mismatch between this assumed speed of sound value and the tissue’s true sound speed will lead to overestimation in PAT target size due to incorrect signal alignment. Using a homogeneous phantom, Fig. 7(a) demonstrates that axial FWHM of targets with the same diameter at different depths is not sensitive to the reconstruction speed of sound. The axial FWHM of the first three targets is smaller than the wire diameter; because of the high optical absorption coefficient of the metal target, target light penetration is shallow and thus photoacoustic signals are generated only near the target surface.⁵⁷ Figure 7(b) shows that lateral FWHM is sensitive to reconstruction sound speed. This is because incorrect speed of sound leads to residual phase modulation across presummation signals even after beamforming delays are applied, and such a phase term increases lateral FWHM. Reconstruction sound speeds higher than the true values result in greater degradation than do sound speeds lower than the true value. For example, underestimating sound speed by $40\mathrm{m}/\mathrm{s}$ ($1460\mathrm{m}/\mathrm{s}$, green curve) results in approximately one-third of the lateral resolution degradation compared with overestimating sound speed by $40\mathrm{m}/\mathrm{s}$ ($1540\mathrm{m}/\mathrm{s}$, gold curve) for depths between 10 and 25 mm. At all speed-of-sound values, both axial and lateral FWHMs increase as the depth increases. This is mainly due to the TMM’s nonlinear acoustic attenuation frequency spectrum,^61,62 which causes a decrease in the spectral bandwidth as shown in Fig. 7(c) and a reduction of central frequency as shown in Fig. 7(d).

Fig. 7

(a) Axial target signal FWHM and (b) lateral FWHM estimated using reconstruction sound speed from 1460 to $1540\mathrm{m}/\mathrm{s}$ in a homogeneous phantom. (c) Spectrum FWHM and (d) central frequency estimated using true sound speed of $1500\mathrm{m}/\mathrm{s}$ in the delay-and-sum PAT reconstruction.

The heterogeneous phantom contains a top layer with a sound speed of $1437\mathrm{m}/\mathrm{s}$ and a bottom layer with a sound speed of $1500\mathrm{m}/\mathrm{s}$. The choice of reconstruction speed of sound affects the estimation of spatial size of the imaging targets in the two layers differently. As in the homogeneous phantom, axial resolution is similar across all reconstruction speeds in the heterogeneous phantom [Fig. 8(a)]. Lateral resolution is more severely degraded along the depth for the first three targets (located in the top layer) than for the fifth and sixth targets (in the bottom layer) [Fig. 8(b)]. Two factors contribute to this phenomenon. First, applying overestimated reconstruction sound speed in the top layer can cause more degradation in lateral FWHM calculation of imaging targets than applying underestimated reconstruction sound speed in the bottom layer. Second, the underestimated sound path length in the deeper layer can be partially compensated by the overestimated path length in the superficial layer when the reconstruction speed falls between the true speeds of the two layers. The lateral FWHM of the fourth target located in the deeper layer shows deviation from an expected monotonic trend and is typically larger than the deeper fifth and sixth targets; this may be due to undulating boundary effects, which will be discussed in Sec. 3.5.

Fig. 8

(a) Axial FWHM and (b) lateral FWHM estimated using reconstruction sound speed from 1437 to $1500\mathrm{m}/\mathrm{s}$ in a heterogeneous phantom.

Target SNR values are shown in Fig. 9(a) for a homogeneous phantom over a reconstruction speed of sound from 1450 to $1540\mathrm{m}/\mathrm{s}$. As expected, target SNR decreases with depth. The highest SNR is achieved using the true speed of sound at $1500\mathrm{m}/\mathrm{s}$. For heterogeneous phantom images reconstructed over a reconstruction speed of sound from 1437 to $1500\mathrm{m}/\mathrm{s}$, shown in Fig. 9(b), target SNR decreases with the increase of the depth, except for the first target, which may have been due to strong ringing effects for the first target. Setting the reconstruction speed of sound close to that of the top layer provides higher target SNR, except for the first target. Target CNRs have similar plots and are not shown here. As an example, target CNRs at $1500\mathrm{m}/\mathrm{s}$ in the homogeneous phantom are provided here. The values are 46.0, 45.3, 43.5, 36.8, 28.8 in dB and decrease as the depth of the targets increases.

Fig. 9

Target SNR estimated in (a) a homogeneous phantom over a reconstruction speed of sound from 1450 to $1540\mathrm{m}/\mathrm{s}$, where true speed of sound is $1500\mathrm{m}/\mathrm{s}$ and (b) in heterogeneous phantom over a reconstruction speed of sound from 1437 to $1500\mathrm{m}/\mathrm{s}$ (where the bounds represent the true values in the top and bottom layers, respectively).

Mismatch between reconstruction sound speed and true speed of sound also causes incorrect estimation of target depth. Figure 10(a) shows that reconstruction sound speeds higher than the true sound speed cause overestimation of target depth as expected (because reconstruction depth is the product of reconstruction sound speed and time delay), whereas values lower than the true sound speed cause underestimation of target depth in a homogenous phantom. The error in reconstruction depth is directly proportional to the depth. In this analysis, we considered target depths reconstructed using the true speed-of-sound value as ground truth, and the depth of the first target was considered as a reference point. Similar depth overestimation can also be seen for the first three targets in the top layer of the heterogeneous phantom, as shown in Fig. 10(b). Underestimated depths for the second target may be due to slight misalignment during the construction process. The estimation error might be partially offset by the accumulative error from the top layer. The estimation errors for the deepest targets (fourth and fifth) are due to a combination of three factors: reconstruction speed mismatch, compensation from the overestimated path length in the superficial layer, and undulating boundary effects.

Fig. 10

Error in axial location estimation for the first five wires in (a) the homogeneous phantom and (b) the heterogeneous phantom using the location of the first target as reference and location estimated using true sound speed of $1500\mathrm{m}/\mathrm{s}$ in the homogenous phantom as ground truth.

3.5.

Boundary Undulation Effects on Photoacoustic Tomography Image Quality

Figure 11 shows ultrasound (a) and PAT (b) images of a heterogeneous phantom containing six PTFE tubes at the same depth as the fourth wire in the above heterogeneous phantom. Assuming that a two-layer phantom with a flat boundary will result in equal size and intensity of each tube, a phantom containing an undulating boundary can illustrate the changes of the dimensions and intensity of the imaged targets as a function of lateral position relative to the undulation pattern. Figure 12(a) indicates that, as expected, axial FWHM values are not significantly different from each other for all reconstruction sound speed values. Figure 12(b) shows that lateral FWHM values of the targets at the same depth depend on the reconstruction speed of sound. When reconstruction speed of sound is between 1437 and $1460\mathrm{m}/\mathrm{s}$, the lateral FWHM values of the middle four targets are similar across the transducer surface. When the reconstruction speed of sound is between 1460 and $1500\mathrm{m}/\mathrm{s}$, the lateral FWHM values are related to the locations of the targets relative to the boundary undulation pattern. The first and sixth targets from left to right have the worst lateral FWHM. This is mainly because (1) reconstruction is performed with dynamic subaperture, (2) only half of the normal subaperture can be used for the targets close to the edge, and (3) lateral FWHM near the focal plane is inversely proportional to subaperture width, which is 23 mm for a focal ratio ($F$ number) of 0.9 at the imaging depth of about 21 mm. Across the different speeds of sound, the second and fifth targets show similar dependences of lateral FWHM values on reconstruction sound speed. This might be due to the similar positions of the two targets relative to the boundary undulations. The fourth target FWHM is more sensitive to the reconstruction sound speed. This might be due to this target being directly under one of the undulation valleys, which has a focusing effect and caused the worst signal alignment in image reconstruction due to the refraction of sound across this valley. According to Snell’s law, when the incident angle changes from 0 deg to 60 deg in the deeper layer with a speed of sound of $1500\mathrm{m}/\mathrm{s}$, the corresponding refraction angle in the top layer with a speed of sound of $1437\mathrm{m}/\mathrm{s}$ changes from 0 deg to 56 deg. In other words, a beam can be bent inward by 0 deg to 4 deg for these incident angles. Such bending was ignored in reconstruction of the PAT image, which was based on delay-and-sum and assumed constant speed of sound and wave propagation along a straight line, causing inaccurate delay compensation before beam summation across the probe surface and a subsequent increase of the lateral FWHM of a target. To further illustrate the effects of sound speed mismatch, we simulated a region ($-7.6\mathrm{mm}<x<12.4\mathrm{mm}$, $0\mathrm{mm}<y<20\mathrm{mm}$) in Fig. 11(a) centered around the fourth target and calculated the beamformer delay error profiles across 64 ultrasound probe elements for this target. For reconstruction sound speed between 1437 to $1460\mathrm{m}/\mathrm{s}$, the low delay error shown in Fig. 12(c) explains the small lateral FWHMs shown in Fig. 12(b), while the largest delay error from using $1500\mathrm{m}/\mathrm{s}$ correctly predicts the largest lateral FWHM of this target. In addition to sound speed variation, the delay error also depends on geometry. Figure 12(d) shows beamformer delay error profiles for the same simulated region with the convex boundary replaced by a planar boundary positioned at 12 mm below the probe surface. The profiles are different from those in Fig. 12(c) as well as those (not shown) generated from the same convex boundary at the same location as the planar boundary. In summary, the speed difference across a boundary, the shape of the boundary, the location of the boundary relative to the ultrasound probe, and the light absorber all play a role in signal misalignments in image reconstruction. This finding also suggests that the undulation also contributes to the large variation in lateral FWHM estimation for the fourth wire at the same depth in the first heterogeneous phantom, as discussed in Sec. 3.4.

Fig. 11

(a) Ultrasound and (b) PAT images of a heterogeneous phantom containing six tubes. Target 1 is the signal farthest to the left, and target 6 is on the far right.

Fig. 12

(a) Axial FWHM, (b) lateral FWHM, beamformer delay error profiles across 64 ultrasound elements for (c) the fourth target for a convex boundary, (d) for a planar boundary positioned at 12 mm below the ultrasound probe, (e) peak intensity, and (f) target SNR in dB of the photoacoustic signal close to the tube top wall reconstructed using sound speed from 1437 to $1500\mathrm{m}/\mathrm{s}$ in a heterogeneous phantom containing six tubes at the same depth. Note: the values of the curves in (a) to (d) increase with the increase of the speed of sound, while the values of the curves in (e) and (d) decrease with the increase of the speed of sound.

Figure 12(e) shows that the peak intensity of the six target signals at the same depth is also sensitive to both the choice of reconstruction speed of sound and acoustic attenuation distribution across the undulating boundary. For the former, for example, the peak intensities of the fourth target decreased with the increase of reconstruction speed of sound. This trend is the opposite of that of the lateral FWHM values versus reconstruction speed of sound. This opposite relationship may signify that the decrease in peak intensity is due to averaging effects in PAT imaging resulting from the corresponding increased target’s lateral size due to mismatched reconstruction sound speed. The second dependence can be observed crossing the middle four targets at one fixed reconstruction sound speed, such as $1460\mathrm{m}/\mathrm{s}$. Even though the lateral and axial FWHMs of the fourth and fifth targets were similar to each other at this speed, the photoacoustic signal intensity of the fifth target reduced to about 83% of that of the fourth target. This decrease may be explained by the relatively larger portion of highly attenuating fibroglandular TMM in the signal propagation path for the fifth target. For example, averaging 1-mm extra fibroglandular TMM thickness can cause about a 12% intensity decrease using the acoustic attenuation difference $5.5\mathrm{dB}/\mathrm{cm}$ between two PVCP formulations at the central frequency of 7.5 MHz. The actual variation of PAT intensity due to acoustic attenuation can be more complex due to the acoustic heterogeneity caused by the undulating boundary. This complexity might also explain the different PAT intensities between the second and fifth targets even though they have similar geometries. As expected, the target SNR shown in Fig. 12(f) has a similar trend as the peak intensity for the middle four targets.