Twente Photoacoustic Mammoscope 2: system overview and three-dimensional vascular network images in healthy breasts

Abstract. We present the Twente Photoacoustic Mammoscope 2, a photoacoustic breast imaging system employing a tomographic configuration. It images one breast pendant inside an imaging tank filled with water while a woman lies prone on a bed. A dual-head laser (755 and 1064 nm) illuminates the breast with one beam directed at the nipple and nine beams directed at the sides. Ultrasound signals are detected using 12 arc-shaped arrays, each curving along the pendant breast. Each array comprises 32 piezocomposite elements each with a center frequency of 1 MHz. The imaging tank and the ultrasound arrays rotate around the breast in steps to obtain additional multiple projections. Three-dimensional images are reconstructed using a filtered backprojection algorithm. The system is described in detail, and measurements on a test object are presented. As part of a preliminary study to assess the system’s in vivo performance, the breasts of two healthy volunteers were imaged. These images show the breast contour, the nipple, and the vascular anatomy within the breast. In the nipple of one case, multiple high-intensity “hot spots” are observed, which we suspect are associated with the lactiferous ducts terminating in the nipple.


Introduction
Breast cancer is the most frequently diagnosed cancer as well as the leading cause of cancer death among women worldwide. 1 Imaging is an important aid in the detection, diagnosis, and treatment monitoring and follow-up of breast cancer. 2 Conventional imaging techniques currently used in the clinic are x-ray mammography (MMG), ultrasonography (US), and magnetic resonance imaging (MRI). 3 MMG is the only method utilized in organized screening programs in several countries around the globe since it is the only method proven to reduce mortality. 4 It does, however, require (often painful) breast compression, makes use of ionizing radiation, and the detection sensitivity in dense breasts is relatively low. 5 US is used only in the diagnostic trajectory and is often employed for targeted examinations of palpable lesions and lesions detected in screening. It is also the primary method used for guiding biopsies. It does not use ionizing radiation but has the drawback of being operator dependent. 6 MRI is used in the diagnostic trajectory, often as a problem-solving modality. It does not use ionizing radiation, has high resolutions, and possesses high sensitivity. However, it requires the administration of contrast agents, does not have a high specificity, is expensive, and is not widely disseminated. 3,7 Each technique has its strengths but also its weaknesses, leaving room for new imaging techniques that can potentially overcome the current limitations and drawbacks.
Photoacoustics (PA) is a technique where the absorption contrast of tissue constituents to light is combined with the resolution of ultrasound imaging. 8,9 Upon illuminating the tissue of interest with pulsed light, optical energy absorption at specific locations leads to a local pressure rise due to the effect of thermoelastic expansion. The pressure transient propagates as ultrasound, which can be measured outside the tissue with ultrasound transducers. Hemoglobin as an absorber provides functional contrast to allow visualization of (newly formed) vascularization associated with cancer progression. The production of new blood vessels is called angiogenesis; inducing angiogenesis is thought to be one of the biological capabilities of tumors and a hallmark of cancer. 10,11 With the first generation Twente Photoacoustic Mammoscope, the PAM 1 system, we have shown visualization of breast malignancies in vivo using PA. [12][13][14] In PAM 1, a woman lay prone on a bed with one breast pendant through an aperture. Below the aperture, the breast was slightly compressed between a glass plate on the cranial side of the breast and a flat detector array on the caudal side of the breast. The planar imaging configuration worked in forward mode, illuminating the breast through the glass plate with 1064 nm and detecting generated ultrasound with a two-dimensional (2-D) ultrasound detector array from the caudal side. In the latest clinical studies, 32 of 33 measured breast malignancies 13 were visualized with good contrast. One of the drawbacks of PAM 1 was the limited field of view (90 × 80 mm 2 detector surface area), precluding whole breast coverage in one measurement. Further, the PAM 1 configuration required slight compression of the breast for acoustic coupling between breast and detector and to obtain a uniform breast thickness. This compression is unfavorable with respect to patient comfort. Moreover, the flat detector array in the forward mode configuration made accurate three-dimensional (3-D) visualization of lesions challenging due to the limited views. To overcome these limitations, a tomographic system was designed and built, which we refer to as the second generation Twente Photoacoustic Mammoscope (PAM 2).
Several other research groups are working in the field of photoacoustic tomographic breast imaging. Matsumoto et al. 15 showed a (breast) imaging system (PAI-04) where a woman is prone on a bed with one breast inside a breast-holding cup. A hemispherical-shaped detector array (center frequency 4 MHz) is used to measure the ultrasound as the breast is illuminated by a large beam, directed at the ventral side of the breast. The system is scanned in a spiral pattern within a plane to provide full breast images. Lin et al. 16 have recently presented an instrument referred to as the single-breath-hold photoacoustic computed tomography breast imaging system, where again a woman lies prone on a bed. The breast to be imaged is slightly compressed upward against the chest wall. A full-ring detector array (2.25 MHz) is scanned elevationally to obtain a 3-D image. Also here, the breast is illuminated by one beam from beneath the breast, directed at the ventral side of the breast. Oraevsky et al. 17 co-registered anatomical US images with functional PA images in their LOUISA-3D system, which images a breast inside an imaging bowl while the woman lies in prone position on a bed. An arc-shaped detector array (50 kHz to 6 MHz) is rotated around the breast while multiple optical fibers positioned on an arc-shaped holder, which can rotate around the breast, are used for illumination.
We contribute to the field of photoacoustic breast imaging with our PAM 2 system, which is different from other systems recently presented [15][16][17] in terms of illumination and detector specifications. We show that adding more illumination beams and changing the ratio of delivered energy to these beams alters the observed results. Also we present images of absorbing features in the nipple that, to our knowledge, have not been reported earlier in photoacoustic breast imaging. The PAM 2 system contains an imaging tank filled with water, in which a breast is pendant while the woman lies in prone position on a bed. The system has low frequency (1 MHz), broadband detector arrays that are arc-shaped and placed around the breast inside the imaging tank, such that the arcs follow the breast contour. Illumination is done not only with a beam (which we call bottom beam) directed at the ventral side of the breast, as in Refs. 15 and 16, but also with nine smaller beams directed at the region of the breast close to the chest wall. The performance of the system is demonstrated on a test object and on the breasts of two healthy volunteers. Reconstructed images show clearly visible vascularization. Images obtained with two illuminating wavelengths (755 and 1064 nm) are qualitatively compared. We give an indication of the imaging depth. In the breasts of one healthy volunteer, we observe, for the first time, structures that might be linked to lactiferous ducts terminating at the nipple.

Light Delivery System
From in vivo measurements with PAM 1, operating at an excitation wavelength of 1064 nm, the absorption contrast could not be attributed to hemoglobin only. Water, lipids, and collagen were also expected to contribute to the contrast. 12 For PAM 2, a dual-wavelength approach was chosen. Excitation with 1064 nm was chosen for its deep penetration into tissue. 18 Water absorption is relatively high at this wavelength, 19 potentially a limiting factor in the penetration depth given the presence of water in the breast and the use of water as a coupling medium between breast and detectors. A second excitation wavelength is needed to distinguish contributions of two chromophores. Excitation with 755 nm was implemented for it has higher hemoglobin absorption and lower water absorption. 20 In PAM 2, the woman lies prone on a bed while one of the breasts is pendant through an aperture, hanging inside the imaging tank. The breast is illuminated with a large laser beam from the bottom (the ventral or nipple side of the breast). There is also illumination via nine optical fiber bundles, which are directed at parts of the breast close to the chest wall. A duallaser system (Qcombo, Quanta System, Milan, Italy) comprising an Nd:YAG and an Alexandrite Q-switched laser is used. The pulse durations are 5 ns at 1064 nm (Nd:YAG) and 60 ns at 755 nm (Alexandrite). The system can be operated in singlewavelength mode or in double-wavelength mode. In singlewavelength mode, one of the wavelengths is emitted with a repetition rate of 10 Hz. In double-wavelength mode, the system alternates between the two wavelengths at a 20-Hz repetition rate.
Each side illumination fiber bundle has a concave lens (focal length of 14.7 mm) to diverge the beam. The fiber bundles are directed slightly upward, under an angle of 5 deg, in order to illuminate the region of the breast close to the chest wall. The beam from the bottom is also diverged by a concave lens (focal length of 25 mm). Inside the laser system, the main laser beams are split into two beams-for bottom and side illumination-by means of a beam splitter (Quanta System, Milan, Italy). By exchanging the beam splitter for another, the energy distribution can be altered. In one setting, 67% of the laser light illuminates the breast through the large bottom beam; the remaining 33% is distributed equally over the nine fiber bundles (67/33 illumination ratio). Another energy division is possible, where 50% of the total energy illuminates from the bottom and the other 50% from the sides through the fiber bundles (50/50 illumination ratio). For practical reasons, only two configurations have been considered so far. The laser pulse energies (measured with StarBright power meter and PE50BF-DIFH-C sensor, Ophir Optronics Solutions Ltd., Jerusalem, Israel) of the bottom illumination beam at a 50/50 configuration are 160 and 130 mJ of the Nd:YAG and Alexandrite laser, respectively. At the 67/33 configuration, the bottom illumination energy was not measured at the same time point but can be expected to be ∼220 and 180 mJ for Nd:YAG and Alexandrite, respectively (with a tolerance of 8%). Laser fluence on the skin is below the safety limit of 20 mJ∕cm 2 at all times (NEN-EN-IEC 60825-1:2014). The smallest safe distances between skin and laser outputs inside the tank were calculated. The design of the imaging tank does not allow the breast to come that close to the laser outputs, ensuring a safe situation. Eye protection goggles are worn whenever the laser is emitting.

Ultrasound Detection
The single elements of the ultrasound detector system (Imasonic SAS, Voray sur l'Ognon, France) have a 3.5 × 3.5 mm 2 active surface area. The distance between two adjacent elements is 1.4 mm, resulting in a center to center distance of 4.9 mm. One array consists of 32 piezocomposite elements, placed on an arc with a radius of 120 mm (Fig. S1 in the Supplementary Materials). A total of 12 such arrays were installed in PAM 2.
A waterproof stainless steel case encloses an array and its integrated preamp (PA Imaging R&D B.V., Enschede, The Netherlands). The 32-channel low-noise ultrasound amplifiers with a 38-dB voltage gain employ a high-input impedance and filtering to prevent aliasing. The outer surface of the elements is covered with aluminum tape with a thickness of 30 μm, to prevent light from being absorbed by the detector surface and thus generating a photoacoustic signal.

Imaging Configuration
The fiber bundles and detector arrays are fixed inside the imaging tank. Detector arrays were placed in a circular arrangement inside the imaging tank, with an angle of 30 deg between two adjacent arrays. Fiber bundles and detector arrays are installed in an interleaved manner, but the flexible design of the tank allows other configurations also (e.g., transmission mode with all detectors on one side). The current configuration is depicted in Fig. 1. During a measurement, the tank with the fiber bundles and detectors rotates in steps around the breast to acquire more views (projections) around the breast for a tomographic scan. This increases the number of sampling locations and thereby improves the reconstructed image quality. At each position, an averaged signal is obtained from a number of laser pulses. This number is a parameter that can be chosen per measurement. After acquiring the average signal, the tank rotates a few degrees before stopping and obtaining another measurement. The choice for such a stop-and-go scanning method over a continuous scanning method 21 was made to limit the time needed for reconstruction (linearly increasing with the number of projections). The centers of the arcs of all detector arrays coincide, such that all single elements lie on the surface of a virtual hemisphere. The elements of six arrays are shifted 2.45 mm on the virtual hemisphere with respect to the elements of the other six arrays, to enhance the spatial sampling of the breast volume. The positions of all 384 detector elements were primarily deduced from the designed geometry of each arc-shaped array and the mounting positions inside the imaging tank (30 deg between two adjacent arrays). Subsequently, small corrections to the elements' positions were made based on photoacoustic calibration measurements with a small absorbing dot made from blank epoxy (G14250, Thorlabs, Inc., Newton, New Jersey) mixed with carbon black, suspended in the imaging tank. In every measurement, the exact angular rotation of the imaging tank between consecutive scanning steps is measured at each projection (MagLine system, SIKO GmbH, Buchenbach, Germany) since the actual motor rotation was found to deviate from the programmed rotation. In Fig. S2 in the Supplementary Materials, the positions (where tank rotation was measured with the MagLine system) of the centers of all 12 × 32 elements are plotted (Lambert azimuthal equal-area projection), for a tomographic scan obtained from 21 projection angles, covering a total angle of 60 deg.
Prior to a measurement, the subject climbs onto the bed via a step stool to lie prone on a bed with one of her breasts through the aperture with a diameter of 21 cm and pendant in the imaging tank. The imaging tank is filled with water to achieve acoustic coupling between the breast and the ultrasound detectors. The bed originally belonged to a prone breast stereotactic biopsy system (Hologic, Inc., Marlborough, Massachusetts). The imaging tank is double-walled with an inner diameter of 35.5 cm and an outer diameter of 40.5 cm. The gap between the two walls carries water heated externally in a heating circulator (F12-ED, Julabo GmbH, Seelbach, Germany). By heat transfer through the inner wall, the coupling water in the imaging tank is warmed up to maintain a temperature of 25°C to 30°C during the measurement. Figure 2 shows the PAM 2 system.

Data Acquisition
Photoacoustic data are acquired through a 12 times 32 channels analog to digital converter as part of a data acquisition (DAQ) system (PA Imaging R&D B.V., Enschede, The Netherlands). Sampling is done with a frequency of 25 MHz at a depth of 14 bits and for 167 μs following each laser pulse, to detect signals from over a distance of 25 cm [assuming a speed of sound (SOS) of 1500 m/s]. Average signals as well as raw signals following each laser pulse are acquired, transferred to and saved on a PC. Raw signals are not used for reconstructing images but are saved to enable more detailed analysis of the data, for example with respect to movement of the breast.

Signal Processing, Image Reconstruction, and Image Visualization
For the reconstruction, a filtered backprojection based on Algorithm 1 by Haltmeier et al. 22 was implemented in Deconvolution is implemented as a regularized division of the measured signal by the average frequency response of the transducers. Regularization refers to an attenuation of the effect of division for frequencies where the average frequency response is very low. The relative sensitivities of all transducer elements are used to scale the measurement data before being passed to the reconstruction algorithm. The SOS is assumed to be homogeneous inside the measurement volume (breast and surrounding water) and calculated from the measured temperature of the water surrounding the breast. 23 To keep reconstruction time within feasible limits, the core of the algorithm was implemented as a C++ extension module with support of OpenMP, a programming interface for parallel computation. 24 A singleimage reconstruction, i.e., for one breast or phantom and one laser wavelength (300 MB of averaged data), takes about 5 to 10 min on a dual-socket server platform with two octa-core Intel Xeon processors with 32 GB RAM. Without the C++ extension, calculation time would be over 5 h per reconstruction. Reconstructed 3-D volumes are by default visualized as maximum intensity projections (MIPs). Another way of visualizing a reconstructed volume is by means of local maximum intensity projections (LMIPs). 25 The first pixel above a set threshold is plotted, where the pixels are color-coded for distance from the observer (white is superficial and orange is further away). The depth corresponding to the orange color is different per image, depending on the size of the reconstructed volume along each axis.

Detector Relative Sensitivity
All detector elements were insonified with a 1-MHz V303 ultrasound transducer (Olympus NDT, Inc., Waltham, Massachusetts), driven at 1 MHz with a fixed pressure in water at room temperature. The V303 transducer was driven by an AFG3102 function generator (Tektronix, Inc., Beaverton, Oregon). The 32 detecting elements of each array were individually insonified by rotating the array around the V303 transducer, keeping the distance between detector and V303 fixed (120 mm). Alignment was optimized by translating the detector array to maximize the amplitude of the detected signal. The average was taken of 100 individual signals. From all measured signals, the amplitudes were extracted and normalized to the maximal value measured in any of the detectors to obtain relative sensitivities.

Detector Frequency Response
The frequency responses of the detecting elements were measured by acquiring signals induced by a custom-built laserinduced ultrasound (LIUS) setup. 26 This technique of LIUS is reviewed in Ref. 27. The frequency content of the LIUS signal was broad enough (frequencies up to 10 MHz were measured when the source was characterized) to assume the signal as an ideal acoustic calibration source. Hence, the measured PAM 2 detector response is considered to only represent the detector response and we did not need to correct for the source frequency content. The LIUS transmitter comprised polydimethylsiloxane mixed with carbon black and applied in a thin layer (29 μm) onto a 5-mm substrate of poly(methyl methacrylate). The LIUS transmitter was illuminated with a 10-mm diameter beam from a frequency-doubled Quanta-Ray Pro 250-10-Hz laser (Spectra-Physics, Santa Clara, California) operating at a wavelength of 532 nm. The 32 elements of each array were individually characterized by rotating the array around the LIUS transmitter, the distance between detector and LIUS transmitter being fixed at 120 mm. Detector arrays, LIUS transmitter, and illumination fiber were immersed in water at room temperature during these measurements. Each signal was averaged 100 times. The acquired LIUS signal was Fourier transformed to obtain the frequency response of the detector.

Resolution
The spatial resolution of the system was estimated by measuring the full-width at half-maximum (FWHM) of the point-spread function (PSF) of a subresolution strongly absorbing sphere. A batch of 26 drops made of two component blank epoxy (G14250, Thorlabs Inc., Newton, New Jersey) mixed with carbon black was produced. These drops were suspended on two perpendicular nonabsorbing nylon threads to create a test object (Fig. 3).

Healthy Volunteer Measurement
All human volunteer measurements were carried out as part of a study approved by an institutional review board (METC Twente, Medisch Spectrum Twente, Enschede) and registered in the Netherlands National Trial Register (TC 6508). To all subjects participating in our study, the procedure was explained and consent was obtained. To test the system and the measurement protocol, a measurement on a 29-year-old woman with a light skin tone was performed (healthy volunteer 1). The 67/33 illumination ratio was employed. A tomographic scan with 21 projection angles (total covered angle of 60 deg) and 100 laser pulses per average was done. Water temperature inside the imaging tank was 30°C. The scan took 4 min per breast. Images were reconstructed with a 0.4-mm grid size and are plotted as MIPs. Differences between images obtained with 755 nm versus 1064 nm were compared.

Optimizing the Illumination Scheme
In the PAM 2 system, illumination of the breast is spread over a total of 10 beams. The main beam illuminates the nipple side of the breast with a diverging beam from the bottom of the tank. The remaining energy is spread over nine fiber bundles that are installed at the top of the imaging tank, serving to illuminate breast tissue close to the chest wall. As mentioned, a beam splitter inside the laser splits the main beam and sets the energy division. In the original configuration, 67% of the energy was delivered to the nipple side of the breast. In order to find out whether a more homogeneous intensity distribution could be obtained, another configuration with more energy going to the side illumination was considered. Both configurations were compared by measuring on the breasts of a 25-year-old healthy volunteer with a light skin tone (healthy volunteer 2).
After the first measurement set with a 67/33 configuration, the beam splitter inside the laser system was replaced with a 50/50 beam splitter and a second measurement set was performed. The subject left the bed between these measurements and thus the breast was repositioned in between measurements. Leaving the bed was necessary since the system needed to be partially disassembled to reach the laser and replace the beam splitter, a procedure that took more than an hour. In all measurements, a tomographic scan over 21 projection angles (total covered angle of 60 deg) was done where signals were averaged over 100 laser pulses. The water temperature was 30°C. The acquisition time for one breast was 4 min. Reconstructed images (grid size = 0.4 mm) are plotted as MIPs. Another reconstruction was performed with a 0.25-mm grid size and plotted as LMIPs. Images were compared, where the focus was on the visibility of blood vessels near the chest wall and the homogeneity of signal intensity throughout the breast. Measurements on this volunteer were also used to estimate the imaging depth relative to the skin surface.

Detector Relative Sensitivity
The sensitivity of all 384 single elements was measured. All values were normalized and the interelement variability in terms of sensitivity is depicted in Table 1.

Detector Frequency Response
The mean impulse response of all LIUS signals, used for measuring the frequency response, is plotted (in red) together with the 384 (12 × 32) individual detector element impulse responses (in gray) in Fig. 4(a). The signals were shifted in time to start at t ¼ 0 s. The 384 corresponding Fourier transformed frequency responses (again in gray tones) can be seen in Fig. 4(b). The mean frequency response is plotted in red. The detector is most sensitive at a frequency of 1 MHz. The average −6-dB fractional frequency bandwidth of all detectors is 100%, ranging from 500 kHz up to 1.5 MHz. There is little variation in response among the individual detector elements (especially in the −6-dB frequency bandwidth). Therefore, the average frequency response was used in the deconvolution step of the reconstruction algorithm.

Resolution
Reconstructed images plotted as MIPs of the test object with epoxy beads can be seen in Figs  . With both wavelengths, the nipple is clearly visible (yellow arrowheads) and blood vessels around it. There is good anatomical correspondence between the images, due to the laser system alternating pulses of both wavelengths within one measurement. Some differences can, however, be observed. When illuminating with 755 nm, the portion of the breast illuminated by the bottom beam seems a bit larger than at 1064 nm [compare the structures indicated with blue arrowheads in Figs. 6(b) and 6(d)]. It is known that the bottom beam divergence of both wavelengths is not exactly the same, explaining this observation. In the 755-nm images, the skin surface around the nipple-areolar complex is better visible than at 1064 nm. This can be explained by the higher absorption coefficient of cutaneous melanosomes at 755 nm compared to 1064 nm. 20 Finally, illumination with 1064 nm shows structures close to the chest wall with more detail, indicated with green arrowheads in Figs. 6(b) and 6(d).

Optimization of the Illumination Scheme
In Fig. 7 Fig. 8(c)]. This could be partially explained by the fact that the x-position of the breast in the 50/50 configuration is more centrally located within the imaging volume. The y-position is, however, less favorable compared to in the 67/33 configuration. In the green box, the differences between the two configurations are small and both zoomed images are comparable. In the magenta box, the 67/33   (12) configuration shows more structures, probably due to the more favorable coordinates both along x as well as y.

Imaging Depth
In the reconstruction of the left breast of healthy volunteer 2, we assessed the imaging depth by tracing a vessel toward the center of the breast. Figure 9 shows a reconstruction done with a smaller grid size, namely 0.25 mm. Images were plotted as LMIPs. Markers (green crosses) were placed along a vessel in the transverse view [ Fig. 9(b)] up to the depth where the vessel was visible. The x-coordinates of those pixels were retrieved and subsequently the same vessel was marked in the sagittal view [ Fig. 9(a)]. From this, the imaging depth in these images was estimated to be 22 mm with respect to the skin surface.

Photoacoustic Features Inside the Nipple
The measurements on healthy volunteer 2, carried out to assess the illumination configuration, led to some interesting observations on the photoacoustic appearance of the nipple. In the images in Fig. 7 [particularly in Figs. 7(c)-7(d)], a number of small high-intensity features in the nipple show up. Figure 10 shows images of the left breast illuminated of the same volunteer (2) at 755 nm with a 50/50 illumination ratio.
Here these features are even more clear. The reconstruction was done with a smaller grid size, namely 0.25 mm, to further investigate the phenomenon. Images were plotted as LMIPs. Figure 10(a) shows the coronal LMIP. The area of interest is shown enlarged in Fig. 10(b). Sagittal and transverse LMIP images are plotted in Figs. 10(c) and 10(d), respectively.

Discussion and Conclusion
With PAM 2, we have noninvasively acquired 3-D images of vasculature in the human breast. Photoacoustic signals are detected in a tomographic scan by 384 single elements placed around the breast with a center frequency of 1 MHz and a −6-dB fractional bandwidth of 100%. The imaging capabilities were explored by measuring on a test object comprising of highly absorbing beads on a thread and by measuring on two healthy volunteers. A resolution of ∼1 mm along all axes was reached. PAM 2 was developed for noninvasively imaging vascularity in the female breast and thereby visualizing breast cancer. It is well-known from the literature that tumor growth induces new vessels to grow in the process of angiogenesis. 11 These blood vessel networks are disorganized; tumor vasculature contains vessels that are immature, tortuous, and hyperpermeable. The vessels are characterized by bulges and blind ends. 28,29 So far, we have shown images of healthy volunteers only. Since we are able to visualize these vessel networks in great detail, the expectation is to be able to identify regions with higher vessel densities and further find abnormal structures and thereby visualize breast cancer.  Additionally, photoacoustic imaging can potentially extract information about blood oxygen saturation using the knowledge that the absorption spectra of hemoglobin and oxyhemoglobin are markedly different at certain wavelengths. It is clinically relevant to assess oxygen saturation of blood vessels in or near a breast tumor since this is expected to be lower in tumor tissue than in healthy tissue. [30][31][32] Matsumoto et al. 15 adopted a method to extract the so-called S-factor from measurements at two  December 2019 • Vol. 24 (12) wavelengths (756 and 797 nm), in which a measure of oxygen saturation could be evaluated. Oraevsky et al. 17 have shown binary functional images with oxygenated and deoxygenated blood obtained by illumination with 757 and 797 nm. We have shown that images obtained with our PAM 2 system at two illumination wavelengths show a good anatomical overlap of visualized structures. The fact that our laser system alternates pulses of the two wavelengths leads to this correspondence between the images obtained at both wavelengths. The isosbestic point of oxy-and deoxyhemoglobin lies around 800 nm. Since our system employs illumination at 755 and 1064 nm, this opens possibilities for analysis of oxygen saturation. For oxygen saturation assessment, it is important to know (at least an estimate of) the light fluence distribution throughout the breast. Image intensities have not been related to absolute pressures nor have been corrected for fluence variations yet. We are planning to implement a method where a fluence distribution estimation for both wavelengths based on the diffusion approximation can be used to evaluate contributions of oxy-and deoxyhemoglobin to the absorption contrast throughout the reconstructed volume. Apart from the possibility for extracting functional information, the two wavelengths are complementary since they result in slightly different images of the same breast volume, for example, the different appearances of the skin surface in images obtained with the two wavelengths.
The breast is not only illuminated with a large bottom beam, but also with multiple side beams directed at the chest wall. Slightly more structures were visualized when delivering a higher amount of energy to the side illumination (at the cost of bottom illumination). It is, however, not easy to directly compare images with both illumination configurations [Figs. 7(a)-7(b) versus Figs. 7(c)-7(d)] since the breast was repositioned in between the measurements and thus its position was not the same. Moreover, the amount of movement of the volunteer in both measurements was difficult to control and complicates comparison. It would be of interest to study the added value of the side illumination and the optimal distribution of light energy and compare it to a situation with only a bottom illumination beam.
Since pendant breast diameters (at the chest wall) can be as high as 13.7 cm (mean for D/DD cup sizes), 33 reaching an imaging depth of around 6 to 7 cm would be required to image the entire breast. Imaging depth was estimated from reconstructed images obtained from healthy volunteer 2 ( Fig. 9) to be 22 mm, which is in the same range as recent results from Toi et al. 34 The estimation is based on the continuity of a vessel, so it does not necessarily reflect the deepest retrieved signal and is, therefore, meant as an indication. We believe that our estimated imaging depth is smaller than what our system would be capable of. Currently, the breast hangs freely inside the imaging tank and this may lead to (part of) the breast not being imaged optimally, when it is outside the region that is being illuminated and/or from where photoacoustic signals are detected. Moreover, the SOS utilized in the reconstruction algorithm is not physically correct. Currently, one SOS is assumed, based on the temperature-dependent SOS in the water surrounding the breast inside the imaging tank. We know, however, that the SOS in breast tissue is different from that of water and that within the breast, different tissue components all have their own SOS. The signals coming from relatively deep in the breast are most affected by a physically incorrect SOS; hence, the image quality deeper in the breast will be decreased most by assuming only one SOS, ultimately limiting our imaging depth. Finally, breathing and motion artifacts may be of influence given the measurement duration. In PAM 2, the breast is pendant in water. From investigating individual time-domain signals, it was observed that the breast drifts during the measurement duration of 4 min. Perhaps volunteers are tense at the beginning of a measurement and slowly relax toward the end of it. Small periodic movements due to breathing were found as well. Abrupt motion of the volunteer is also possible but is always registered by the PAM 2 operator and taken into account when analyzing reconstructed images. In order to enhance depth performance, individual time-domain signals might be further studied and employed to make corrections for motion in the future. Another option would be to immobilize the breast during a measurement. Oraevsky et al. 17 used a hemispherical, optically clear and acoustically thin plastic cup-stabilizer. Lin et al. 16 used an agar pillow to slightly compress the breast against the chest wall from below and measure within a single-breath hold to minimize breathing motion distortion. Matsumoto et al. 15 used a breast holding cup to immobilize the breast. When choosing any form of breast immobilization, a challenge would be to find a shape that supports every breast shape enough to reduce motion artifacts while at the same time having all the breast tissue anatomically available to be imaged. Especially imaging the areas of the breast close to the chest wall and axilla might be compromised when a suboptimal immobilization strategy is chosen. In designing such an immobilization strategy, it is important that a material is used that distorts the optical illumination and acoustic wave propagation as little as possible. Immobilizing the breast could also contribute to a larger portion of the breast being visualized by forcing the breast to always be in the center of the imaging tank. The position of the breast is currently not reproducible [compare breast positions in Figs. 7(a)-7(b) and 7(c)-7(d)] since the breast hangs freely inside the imaging tank. Furthermore, with an immobilization strategy, it will be easier to employ a simple two SOS model in the reconstruction algorithm once the watertissue boundary is known. This would improve our reconstruction algorithm and the visibility of deep features. Also having recognized the skin surface will simplify calculation of the depth of features with respect to the skin. In the long run, it will be beneficial to obtain a 3-D SOS map of the breast tissue simultaneously with the photoacoustic measurement, such that the map can be used to further improve the image reconstruction of that particular measurement. Efforts are made in this field, for example, by Matthews and Anastasio, 35 who proposed a method for a joint reconstruction of photoacoustic initial pressure and SOS distributions from photoacoustic and few-view ultrasound tomography measurements. PAM 2 depth performance will also be enhanced by employing an iterative reconstruction algorithm in addition to the currently used filtered backprojection algorithm. The relatively weak signals from deeper inside the breast will be retrieved from the background signal and enhanced. Both the improvement of the reconstruction algorithm as well as the immobilization of the breast are topics of current investigation. In volunteer 2, we observed intriguing features inside the nipples of both breasts. Figures 7(c)-7(d) show the features in the right breast, and Fig. 10 shows the same features with more detail in the left breast. We hypothesize that these features belong to the lactiferous ducts terminating from the nipple through orifices. The lactiferous ducts in the nipple are surrounded by smooth muscle cells, 36,37 which are vascularized. In general, ducts in the human breast are surrounded by a high number of microvessels. 38 We can observe a total of either seven or eight of these high-intensity features. From previous work, 39 the average breast was estimated to contain five to nine ductal orifices. The features are within an area with a diameter of 12 mm. The nipple (excluding the areola) of volunteer 2 was measured to be 13 mm. Given that these features are thus inside the nipple and the number matches the estimated number of ductal orifices, it is plausible that they indeed belong to the lactiferous duct orifices. This is worth exploring further in the near future. Moreover, it could be valuable to learn more about the anatomy of the nipple in a noninvasive manner, with PA. Knowledge about (ductal and vascular) nipple anatomy is of importance for ductal approaches to diagnosis, planning nipplesparing mastectomy 40 and for understanding normal breast development. 41 Finally, in the literature, there is still some discussion about and a large spread in the number of orifices in the nipple. 39,40,42,43 More measurements, if possible on the same volunteer but also on other volunteers, will be carried out to further investigate these findings.
With PAM 2, we have demonstrated 3-D images of vascular structures in human breasts. This vascular-morphological information combined with the potential of extracting functional information on oxygen saturation in the future could be of clinical use not only in breast cancer management but also in visualizing different arterial and venous disorders and abnormalities of the breast. 44 Although most of them uncommon, accurate diagnosis often demands multimodality imaging and pathologic correlation. Some benign tumors are vascularized (such as hemangiomas), which might complicate differentiation from malignant lesions. It is important to investigate different vascular disorders and abnormalities of the breast and whether PA could have an added value in this field.
PAM 2 system's parameters are unique in terms of the illumination scheme and detector characteristics, and it is worth exploring how it performs in measuring on breast tumors in the near future and how results compare with the other systems around the world. So far, promising results have been obtained in the form of 3-D images of vascular patterns in the breasts of healthy volunteers. The knowledge gained from exploring the system and measuring on healthy women will be used for imaging on breast cancer patients in a noninvasive and painless manner.

Disclosures
Authors Schoustra, Piras, Steenbergen, and Manohar have no relevant financial interests in this paper and no potential conflicts of interest to disclose. PA Imaging R&D B.V. is a 100% subsidiary of PA Imaging Holding B.V. PA Imaging Holding B.V. is a privately held company. Authors Huijink, op 't Root, Alink, and Muller Kobold (employees of PA Imaging R&D B.V.) do not have financial interests in the company.