A tutorial in photoacoustic microscopy and tomography signal processing methods

Section II will provide a broad overview of photoacoustic (PA) imaging, briefly discussing its physical principles, imaging implementations (categorized by the spatial resolution), rudimentary signal processing, technical requirements for building microscopy, and tomography systems and a highlight on biological applications. Section III will focus on signal analysis methods that can be utilized in both optical and acoustic resolution PA imaging. We discuss the role that a PA signal plays in the reconstruction of images, the structural and functional information encoded within it, and highlight several applications in microscopy and tomography. Finally, we provide our perspective on how one can perform PA tissue characterization across multiple biological length scales.

The origin of PA imaging (also referred to as optoacoustic imaging) can be traced back to Bell's accidental 1880 observation that chopped sunlight incident upon a solar cell resulted in the production of audible sound.¹ Acoustic waves are produced through the PA effect via three energy conversions,^2–4 as illustrated in Fig. 1. The incident electromagnetic energy is converted to thermal energy through rapid absorption by chromophores. This thermal energy is subsequently converted to mechanical energy through the thermoelastic effect, creating a pressure wave that propagates away from the chromophores. For the efficient conversion of the incident optical energy to an acoustic wave, the laser pulse length must be shorter than two important time constants in order to satisfy thermal and stress confinement (also known as heat and acoustic relaxation times).^5,6 These time constants quantify the time required for the heat and stress to propagate across the absorbing structure, respectively. Outside of these time scales, heat and stress start propagating before the pulse can be absorbed, resulting in a PA wave with a lower amplitude than what is optimally possible.^7,8

Since the late 1980s, advances in light delivery technologies have contributed to PA imaging systems that are more mobile, tunable, compact, and affordable.⁹ The scalability of PA imaging illumination and acoustic detection schemes has given rise to multiple configurations that broadly can be divided into optical or acoustic resolution.¹⁰ The focusing of the laser source, the choice of the acoustic transducer, and the physical orientation of these components relative to one another dictate the biological length scale that can be probed with a given configuration.¹¹ To that end, we will divide our description of PA imaging into following two broad categories:

1. Photoacoustic microscopy (PAM), where a focused laser beam is used to achieve high spatial resolution imaging at shallow depths; and

2. Photoacoustic tomography (PAT) where broader illumination schemes are employed for deeper penetration and the ultrasound transducer dictates the spatial resolution.

Both PAM and PAT can be performed with either an optical or an acoustic resolution configuration, where the spatial resolution is dictated by the optical focal spot or the transducer/array imaging characteristics, respectively. Figure 2 depicts typical examples of PAM and PAT setups and representative images acquired using each configuration. In optical resolution photoacoustic microscopy (OR-PAM) [Fig. 2(a)], the optical excitation volume is smaller than the acoustic resolution volume, and the lateral spatial resolution is defined by the laser spot size at the objective focus.¹² This approach can achieve sub-micrometer lateral spatial resolution, but it is limited to a few hundred micrometers in depth in order to maintain a small optical focal spot from which the PA signal is generated.^13,14 Such exquisite lateral resolution comes across in the image of mouse ear vasculature shown in Fig. 2(b), where vessels as small as 4.2 μm can be resolved through optically focusing the laser beam.¹⁵ This representative image is one of many examples of OR-PAM that can resolve individual vessels,¹⁶ with some variants reporting resolutions down to 88 nm.¹⁷ In these cases, the optical resolution volume can be thought to be defined by cylinder whose base is the optical focal spot and whose height is the acoustic axial resolution [Fig. 2(c)]. When the optical excitation field is greater than the acoustic resolution volume, the spatial resolution of the system is dictated by the acoustic transducer or the acoustic resolution volume in the imaging approach.¹⁸ This variant is referred to as acoustic resolution photoacoustic microscopy (AR-PAM) and will be described in this Tutorial (Sec. III C). Several PAM configurations exist, where the laser and the transducer are parallel [e.g., as in Fig. 2(a)], opposite each other, or the laser can be focused through the transducer itself.

In PAT, deeper volumetric imaging is achieved through broad illumination schemes [Fig. 2(d)]. Numerous configurations of optical illumination and acoustic detection for PAT exist.²¹ They can be broadly categorized by their use of circular/hemispherical transducer arrays^22,23 or linear/hand-held, multi-element planar transducers where broad optical illumination can achieve sub-millimeter resolution at depths of several centimeters through beamforming.^24,25 The hemispherical transducer arrays overcome the limited view artifacts that arise from linear arrays and can minimize the distortions of the target morphology.²⁶  Figure 2(e) shows a representative image from the vasculature of a mouse bearing a breast cancer subcutaneous tumor.²⁰ Unlike OR-PAM, where individual vessels can be resolved, the broad illumination and acoustic detection schemes employed in this setup cannot typically resolve the smaller vessels of the vasculature; however, the lower acoustic frequency and broader illumination greatly increase the depths from which PA signals can be acquired. As we will demonstrate in Sec. III C, these non-resolvable structures can still be evaluated by extracting encoded information in the PA signals using signal analysis techniques adapted from quantitative ultrasound (QUS).

PAM has made significant contributions to fundamental life sciences research at biological length scales ranging from single blood vessels to single cells and their subcellular constituents [Fig. 2(a)].^27,28 OR-PAM can map the spatial variations of hemoglobin oxygen saturation (sO₂) in the subcutaneous microvasculature of rats in vivo.²⁹ High-speed, high frame rate OR-PAM can measure the release of oxygen from single red blood cells in the mouse brain.³⁰ It can also explore the intrinsic absorption of DNA and RNA to perform label-free, histology-like assessment of tumor specimen margins.³¹ Functional PAM can be used to map the resting-state neuronal functional connectivity in the mouse brain.³² Genetically encoded, photo-switchable molecular probes with near infrared absorption can be used as contrast agents in PAM.^33,34 Furthermore, PAM can be used quantify the intratumoral metabolic heterogeneity at a single cell level³⁵ and study early tumor vascular malformations.³⁶ The left side of Fig. 3 provides an overview of the PAM applications that will be discussed in this Tutorial. Readers are directed to comprehensive PAM reviews for additional applications.^37–39 

PAT has made significant strides toward clinical translation, achieving imaging depths far beyond those of OR-PAM.^9,40,41 Breast cancer diagnosis through PA imaging of hypoxia and angiogenesis has the potential to replace mammography.^41–44 Portable imaging systems have deployed the superficial detection of oral cancers,⁴⁵ skin cancer/eczema/dermatitis⁴⁶ and determining the excision margin for improving the prognosis of melanomas.⁴⁷ Acoustic resolution PAT can be used to examine thyroid tumor vasculature to a depth of 2 cm.⁴⁸ Multispectral PAT has shown clinical utility in the assessment of Crohn's disease by evaluating the hemoglobin content present in intestinal wall inflammation.⁴⁹ More recently, a hand-held PA imaging probe allowed for a 2 cm³ field of view imaging of the human carotid artery, opening doors to early assessment of stroke risk through plaque imaging.²³ Lastly, the use of clinically approved optical contrast agents such as indocyanine green for PA imaging aids in the identification of sentinel lymph node metastases from melanoma patients.⁵⁰ The right side of Fig. 3 summarizes the PAT applications in this Tutorial.

The waves emitted from an absorber due to the PA effect are broadband in the time domain.⁴ In conventional PA systems, these waves are detected by an ultrasonic transducer operating in receive mode and are converted into an analog electrical signal via the piezoelectric effect. To use the analog signals from the transducer for practical signal processing techniques, it is necessary to digitize them using a data acquisition system such as an analog–digital converter (ADC) or digitizer. While specific recommendations for the characteristics of these devices for different types of PA systems are described in the subsequent section, here we briefly define three of the most important characteristics for PA systems: the sampling rate, bit depth, and dynamic range.

The sampling rate is the number of times that the analog signal is sampled per unit time. This is typically measured in millions of samples per second (MS/s) or billion samples per second (GS/s). Using a digitizer with an appropriate sampling rate is crucial, as if it is too low, unwanted effects such as aliasing can occur. The bit depth of the digitizer is the number of quantization levels that are available to characterize the signal amplitude. A one-bit digitizer would transform an analog signal into a digital signal that only attained two values (e.g., 0 or 1). Each additional bit increases the number of quantization levels for the signal amplitude by a factor of two. Thus, more bits provide a more accurate digital representation of the original analog signal. The dynamic range sets the maximum and minimum bounds for the digitized signal amplitude. If the analog signal exceeds either of these extremes, the digitizer will saturate, and the amplitude of the signal will be clipped. Using a larger dynamic range will reduce the chances of clipping; however, for the same fixed bit depth, it will also result in coarser amplitude quantization and may increase the amplitude of the noise in the digitized signal.

The digitized PA signal can be thought of as an N × 1 vector of values representing the amplitude of the PA wave reaching the transducer surface at a given point in time after the excitation laser pulse. In this form, the signals can be referred to as being “raw”—no additional processing has yet been performed. A common first step in signal processing is to apply techniques that can be used to reduce the prominence of noise in the signal, thereby increasing the signal-to-noise ratio (SNR). This is usually achieved by acquiring several digitized PA signals originating from the same spatial location and then averaging them. This operation results in “averaging out” of the random fluctuations in the signal and thereby emphasizing the real signal components. Artifacts in the RF lines that are constant in time, such as those caused by reflections within the transducer⁵¹ or from the laser impinging directly on the piezoelectric element⁵² can be characterized and subtracted from the signal when performing the analysis. Conventional ultrasound signal processing techniques, such as thresholding, gating, or enveloping, can be applied to further improve PA signals,⁵³ and software bandpass filters can be applied to reduce noise from vibrations or from non-zero initial or terminal values of the signal (which could be perceived as a discontinuity and add high-frequency signal artifacts). This latter approach enhances the appearance of the digitized PA signal at the expense of discarding some genuine signal energy. While these techniques are commonly used in PA signal analysis approaches, they can also be performed prior to using the signals to generate an image (e.g., beamforming, MAP) or in more specialized applications such as spectral unmixing.^54,55

In some cases, it is of interest to examine the features of the PA wave in the frequency domain. As shown in Sec. III, while the infinite bandwidth representation of signals in the time domain can look drastically different for objects of different size/shape, when convolved with the system response of the transducer, the signals are often indistinguishable in the time domain. However, this is not the case in the frequency domain, and many features that can be used to characterize the absorber are readily apparent. To obtain this frequency data, a discrete Fourier transform of the digitized time domain signal is calculated, most commonly using the fast Fourier transform (FFT) algorithm. The FFT algorithm will produce a vector containing the amplitudes of sine or cosine waves of a given frequency that can be added to reconstruct the original digitized time domain signal. The description of the FFT algorithm is beyond the scope of this work; however, comprehensive references on this topic in the context of signal processing exist,⁵⁶ and most commercial computational software has native implementations of this algorithm built in.

The spectral resolution of the FFT (i.e., the width of each frequency bin) is dependent on the length, N, of the digitized signal. If the sampling frequency (in samples/s) is denoted as f_s, applying the FFT to a signal of length N will produce a spectrum with a resolution of f_s/N that ranges from −f_s/2 to f_s/2. Thus, while larger sampling frequencies give larger ranges of usable signals, they also provide coarser spectral resolution for the same N. Increasing the length of the signal (typically to a power of 2) by appending zeros to the end of the digitized time domain vector is known as zero padding. Since the maximum frequency that can be obtained is f_s/2, increasing N by zero padding has the same effect as interpolating the frequency spectrum by splitting its spectrum into smaller frequency bins. However, while this gives the spectrum a smoother appearance, it does not increase the spectral resolution of the FFT. The ability of the FFT to display two unique closely spaced frequency components is dependent on the original number of samples in the signal (i.e., N) even in the case of zero padding. If two spectral features cannot be discerned from one another in the original spectrum due to insufficient spectral resolution, no amount of zero padding will result in both features being recovered.

This section will provide an overview of the technical requirements for a PA imaging system that uses piezoelectric transducers for acoustic detection. All-optical detection systems such as Fabry-Pérot sensors⁵⁷ or non-interferometric PA remote sensing⁵⁸ have been successful in a wide range of PA imaging applications.^59,60 The reader is directed toward relevant reviews on the subject.^61–64 

Transducers: Transducers that are used for PAT applications typically operate between 1 and 70 MHz to reconstruct a PA backprojection image of the field of view.⁶⁵ Technical specifications in PAT are determined by the type (e.g., linear array) and geometry (e.g., spherical, cylindrical, and planar) of the acoustic transducer in relation to the sample. The spherical and cylindrical detection geometries involve the use of transducer arrays, arranged in a circular or hemispherical fashion using 128, 256, or 512 individual detection elements.^11,66–69 Often made of piezoelectric zirconate titanate, these transducers operate at 2–5 MHz center frequencies and require 20–40 MS/s sampling rates (i.e., the emitted photoacoustic waveform must be temporally sampled 20–40 million times per second or at a rate of 20–40 MHz). The sampling rate required for a given system is dictated by the maximum frequency of the transducer being used, and at a minimum must be higher than the Nyquist sampling rate. For this reason, using the same system with a higher frequency transducer will usually require a higher sampling also. For example, PA systems operating in the 1–15 MHz range require digitization at 40–60 MS/s, while higher frequencies in the 15–70 MHz range require sampling rates exceeding 300 MS/s. Their applications are limited to breast⁴² or small animal imaging.⁷⁰ The most commonly used inversion approach is the universal back-projection algorithm.⁷¹ Planar transducer detection geometries offer limited view reconstructions but are far more versatile and provide access to a greater range of anatomical targets. These PA probes resemble conventional clinical ultrasound (US) transducers in terms of both dimension (128 or 256 linear array elements) and frequency range (1–70 MHz).^41,49,72,73 QUS techniques typically depend on high-speed analog to digital sampling of the acoustic signals, typically on the order of ten times the central frequency of the transducer. This ensures that there is adequate sampling for the analysis of the power spectra of the ultrasound backscatter signals. It is possible to combine a laser with commercially available transducers, thus achieving acoustic coupling with the target and real-time reconstruction while scanning the probe.^24,74,75 This is a relatively inexpensive means of implementing PA imaging since it exploits hardware beamformers and acquisition electronics advances that are commonly used in diagnostic US imaging.

Most PAM systems use US transducers in the 1–100 MHz range with tightly focused lasers to achieve good imaging resolution and sufficient penetration depth [Fig. 4(a)]. The acoustic attenuation increases as a function of transducer frequency and the penetration depth decreases. When using analysis techniques that require PA frequencies greater than 100 MHz, such as those applied when analyzing micrometer-sized samples, wide transducer bandwidths are required, typically in the hundreds of MHz. These wide bandwidths are important for two reasons. First, axial resolution is inversely proportional to the transducer bandwidth, enabling axial resolution on the order of 1–10 μm at frequencies up to 1000 MHz.⁷⁶ Second, the bandwidth also determines the range of usable frequency content in the PA signal (i.e., the minimum and maximum frequencies); a wide bandwidth improves the frequency-domain signal analysis techniques, such as those described in subsequent sections. Several companies produce ultra-high frequency (UHF) acoustic microscopy systems, but most cannot be used for PA imaging. The equipment required to achieve UHF-PAM can be a barrier to entry, as the UHF transducers, digitizers, and other hardware can be expensive. Transducers with frequencies over 100 MHz are available typically as a single element and can be made from a variety of piezoelectric materials.^77,78 Up to about 100 MHz, PVDF is used as the piezoelectric material, and commercial transducers are relatively inexpensive (<$1000 each) and readily available. These transducers also lend themselves to in-house fabrication, as the thin PVDF film is flexible and easily conforms to shape.^79–81 Over 100 MHz, transducers use LiNO₃^82–84 or ZnO^85–88 as piezoelectric materials and typically have precisely designed single-crystal buffer rods and glass matching layers. Due to high acoustic attenuation that occurs at these frequencies, UHF transducers must also be designed to have extremely high sensitivity. The added complexity means that these UHF transducers are difficult to manufacture and can cost up to $20 000 at central frequencies up to 1000 MHz.

A key requirement for UHF signal acquisition is a high sampling rate. The Nyquist theorem states that the sampling rate must be at least twice the maximum resolvable frequency within the signal.^89,90 However, for typical signal analysis techniques, sampling at 8–10× the maximum frequency is needed to accurately capture the PA time-domain signal with high fidelity required for signal analysis techniques. For example, if using a 200 MHz transducer with 100 MHz bandwidth (i.e., 150–250 MHz), the minimum system sampling rate required to satisfy the Nyquist theorem is 500 MS/s; however, in practice, 1.5–2.0 GS/s is desirable to faithfully recreate the time domain signal. Recent studies utilizing compressed sensing are currently being developed to help ease these requirements.⁹¹ The bit depth is also an important parameter with analog-to-digital converters (ADCs). Higher bit depth provides better dynamic range and SNR, particularly when high amplitude signals are acquired. ADCs with 10- or 12-bit rates over a ±5 V range are ideal. However, ADCs with both GS/s sampling rates and high bit depth can be very expensive, costing up to $40 000. Additionally, a computer capable of fast data transfer is required, as the transfer rates between the ADC and computer can approach hundreds of MB/s.

Lasers: Another crucial component of any PA system is the laser used for generation of the PA signal. In PAT, diffuse illumination schemes are typically used [Fig. 2(d)] at energies of ∼10–20 mJ/pulse, with 10–20 Hz pulse repetition rates and <10 ns pulse width.⁷² In UHF PAM, samples are typically microscopic and thus require nanosecond or picosecond pulse widths to satisfy the stress and thermal confinement conditions. Here, samples are typically scanned in a raster pattern, which requires moving the sample or the laser/transducer to create the image. High repetition rates (>1 kHz, ideally >10 kHz) are required to ensure reasonable scan times. The laser energy depends on the type of imaging performed. When the laser is focused to a diffraction limited spot for OR-PAM, low energies of 1–100 nJ/pulse are generally adequate, assuming sufficient chromophore absorption and concentration. When larger areas require illumination, then greater energies are required, typically 0.1–1 μJ/pulse or more.

Several types of lasers are available for PA imaging. Solid state lasers such as Nd:YAG typically have excellent pulse energy (microjoules to millijoules per pulse), short pulse widths (<10 ns), and are capable of achieving high repetition rates (1 kHz) but are limited to 1064 nm and harmonics (e.g., 532 and 355 nm).^92–94 Lasers that satisfy these requirements and can cost between $10 000 and $50 000. Diode lasers are much less expensive and can be found through a range of wavelengths throughout the visible and near-IR spectrum but typically have longer pulse widths (over 100 ns), low energy (<100 nJ/pulse), and require complex circuitry to enable rapid nanosecond-switching times. Recent advances in diode lasers have made them attractive for PA imaging.^95–99 Dye lasers are also available but tend to have very high energies (millijoules per pulse) and low repetition rates (<100 Hz), are expensive, and are less ideal compared to other available lasers. For PAT applications, the laser pulse repetition frequency often limits the imaging frame rate. Most optical parametric oscillators employed for multiwavelength PAT applications in the near infrared range are limited to about 10–20 Hz laser repetition rates.²⁴ The need for signal averaging to increase SNR often results in low Hz frame rates. Another important consideration is the choice of laser wavelength. This becomes particularly important for biomedical applications of both PAM and PAT where a balance between penetration depth and chromophore absorption must be carefully considered.^27, Figure 4(b) shows how the laser penetration in skin as a function of wavelength dictates the type of application that PA imaging can target. For instance, in the UV ranges (100–400 nm), only PAM performed at submicrometer depths is possible and some applications target the nucleic acids as the structure that generates the PA contrast.³¹ Within the visible range (400–800 nm), most PAT applications take the advantage of lower water absorption. At the near infrared range (>800 nm), the optical absorption of lipids and collagens becomes the dominant source of the PA signal and drives the respective applications.¹⁰⁰ 

Dual US-PA systems: It is often advantageous to construct a dual US and PA system that utilizes the same transducer. US imaging can provide structural/anatomical information about the target, complementing the multiwavelength PA imaging that typically captures functional information (e.g., red blood cell oxygenation). It is possible to implement dual US-PA systems in both PAM and PAT systems.

In PAT imaging systems, US pulse-echo imaging is often performed alongside PA imaging.^41,49,72 The detection of the US waves produced from PA and pulse-echo are performed in an interleaved fashion, achieving frame rates that are limited by laser pulse repetition (typically 20 Hz). The US imaging electronics architectures are robust and mature enough to be used for acquiring and processing PA signals. Synchronization with the laser Q-switch trigger is required to ensure that the acquisition timing of the PA signals does not interfere with the pulse-echo emission. Lately, it is possible to perform this using commercially available US imaging platforms that provide access to the acquisition architecture. Other modes of US imaging acquisitions (e.g., Doppler or contrast enhanced modes) can be used in conjunction with the PA acquisitions for a variety of biomedical applications.¹⁰² 

In most PAM systems, US imaging is not performed as the resolution of the US transducers typically used is insufficient for detecting the fine features in the sample that could potentially be of interest. But with an UHF transducer, US imaging resolution can range from 10 μm at 100 MHz to nearly 1 μm at 1 GHz with very wide bandwidth signals, ideally suited for signal analysis techniques.¹⁰³ At these frequencies, the lateral and axial resolution of PA and US are similar and ideal for imaging single cells. To enable US imaging in a PAM system, pulse echo circuitry must be added, increasing the hardware complexity.^82,104–106 This requires a pulse generator capable of 50–100 Vpp monocycle pulses at the desired center frequencies. The transducer must now transmit as well as receive, and a fast switch or electronic-matching circuitry is required to control the sending and receiving signals to the transducer. A low jitter (<1 ns) trigger source should be used to coordinate the pulse generator, switch, laser, and digitizer. These components can cost >$5000; however, they add an additional complementary imaging modality to the system.

Once a PA wave is generated, the propagation and interaction of the wave the surrounding media are defined by US physics.¹⁰⁷ The effects of frequency-dependent attenuation; diffraction; and the transducer focus, geometry, center frequency, and bandwidth will affect PA imaging in the same manner as with US imaging.¹⁰⁸ These physical effects, along with the illumination scheme, alter the PA signal and play a key role when selecting the optimal image reconstruction algorithm for a given application.¹⁰⁹ A detailed treatment of PA imaging reconstruction approaches can be found in prior tutorials and reviews that describe universal backprojection and multiple inversion schemes^5,110 or time reversal reconstruction approaches.^111,112

In most PA imaging applications, only the amplitude of the time-domain signals is utilized to generate images. In doing so, quantitative information that can be obtained from the time domain signal by using signal processing methods is lost. Signal processing methods can be used to extract information on the structure (size, shape, and orientation) of the underlying absorber structure.¹¹³ Waveforms corresponding to the analytical solutions of the time-domain PA wave equation alongside with their frequency spectra (obtained by the Fourier Transform) are shown in Fig. 3 for simple geometries. Unlike the backscattered US signals that predominantly contain frequencies present in the initial transmit pulse, PA signals are inherently broadband.¹¹³ However, ultrasonic transducers typically have limited bandwidth. Attenuation during the propagation of the PA wave suppresses the high-frequency components, altering the PA waveforms. Tissues, thus, behave as a low pass filter for the propagating US waves.

In Fig. 5, the PA waveform emitted by an infinite layer, cylinder, and sphere are shown for unlimited bandwidth (solid line) and band-limited (dashed line) waveforms. The band-limited signals for each structure appear similar, which makes it difficult to differentiate between objects with different morphologies.¹¹⁴ However, the spectra of these signals are distinctly different, with minima and maxima occurring at specific frequencies that depend on the shape, size, and acoustic properties of the source object that produced them.^4,113 In AR PA imaging, numerous PA sources of various sizes, shapes, and orientations may be present within the illumination/resolution volume. Broad illumination combined with sparse source concentration results in the constructive and destructive interference of the bandlimited PA signals, which can lead to the formation of PA speckle as observed in Fig. 2(e).^108,115 A comparison of measured spectral patterns to theoretical models allows for studying the size and/or morphological properties of single cells and structures in small animals using both PAM and PAT.²⁸ 

Next, we cover signal analysis techniques that can be applied to PAM and PAT acquisitions to extract details that cannot be obtained from the maximum-amplitude generated images alone. Rather, these approaches exploit the additional information encoded in these signals to probe the morphological and functional properties of the underlying PA source. It is important to distinguish at this point the various PA sources that these approaches are capable of interrogating. A chromophore refers to a light-absorbing molecule such as Hb that can produce a PA signal. An absorber refers to a collection of chromophores that can act, for all practical purposes, as “one” PA signal source, such as a red blood cell (RBC) consisting of a collection of Hb molecules. An absorbing structure refers to a collection of absorbers that are localized in space, thus acting as “one” absorbing structure insofar as the detection device can sense. An example of this are blood vessels consisting of numerous red blood cells (RBCs). Even though each RBC is a source of the PA waves, a detection device (transducer) may not be able to resolve the contribution of each RBC within a blood vessel.¹¹⁶ 

Once the US and PA signals are acquired, signal processing techniques can be used to extract information that cannot be acquired through resolving anatomical structures. To improve SNR, multiple signals from the same spatial location are averaged to remove background noise. A background signal (acquired before or after the main sample acquisition) can be subtracted to remove persistent artifacts. A bandpass filter is applied to reduce system noise outside of the transducer bandwidth, and then the desired signal within the acquisition is gated. Finally, the spectrum is calculated using the FFT introduced in Sec. II D. In this section, we will describe three frequency-based approaches for analyzing signals in PAM. These include (i) frequency spectrum fitting for estimating the absorber size, (ii) frequency band equalization for reconstructing images at various length scales through frequency filtering, and (iii) F-mode for introducing a new source of frequency-based contrast. Following each technique, we will discuss PAM applications where these analysis techniques can extract information about the cells and tissues that cannot be obtained through conventional imaging approaches.

In PAM, measured PA or US signals are compared to established theoretical models (such as those in Fig. 5) to determine information about the size or shape of single cells. Cells are inhomogeneous, viscoelastic, nucleated, and composed of organelles scattered throughout; however, over the length scales commonly used in UHF microscopy (8–4 μm at 200–500 MHz), they can be considered roughly homogeneous. Most cells have negligible endogenous absorption in the visible wavelengths; however, melanoma cells contain melanin, which strongly absorbs light through the UV-near IR spectrum. Figures 6(a) and 6(b) show the PA signal and spectrum, respectively, from a single melanoma cell measured using a 375 MHz transducer and a 532 nm laser. Spectral minima at specific frequencies result in a unique pattern that occurs due to the size and acoustic properties of the cell. Using the known properties of cells (speed of sound of roughly 1570 m/s and a density of 1050 kg/m³),^117,118 the PA signal and spectrum were calculated using a theoretical model,^119,120 and compared to the acquired signal as shown in Fig. 6(b). By matching the frequencies of the minima and maxima to the theoretical models, the cell size was determined to be 17.2 μm. Using pulse echo US, the same phenomena occur, where frequency minima exist in the US backscatter from the cell.^121–123 As with the PA spectrum, the frequency minima in the US spectrum depend on the cell size and acoustic properties. Figures 6(c) and 6(d) show the US signal and spectrum from the same melanoma cell, respectively. By comparing the spectrum to the US scattering model,¹²¹ the size was determined to be 16.9 μm. In this work, these spectral US and PA techniques were used to measure 47 melanoma cells. Excellent agreement between the average diameter predicted using either modality was observed (17.1 ± 4.4 μm using US, and 17.0 ± 4.6 μm using PA), which also compared well with optical microscopy size estimates. The same technique was applied to measure the size of perfluorocarbon droplets containing gold nanoparticles in liquid and on a solid boundary.^116,124

The frequency spectrum fitting method can be extended to measure the nucleus-to-cytoplasmic (N:C) ratio of single cells. Cellular features, such as the N:C ratio, are an important metric in tumor grading and are typically evaluated using stained histological samples.¹²⁵ In one example, MCF-7, PC-3, and MDA-MB-231 cancer cells were stained with DRAQ5, a supravital fluorescent dye that stained cell nuclei. Using the same methodology as above, the cell diameter was measured using US, and the nuclear diameter was measured using PA with a 375 MHz transducer and a 532 nm laser. An US signal backscattered from an MDA-MB-231 cell and a PA signal emitted from its dyed nucleus are shown in Figs. 6(e) and 6(g), respectively. Figures 6(f) and 6(h) show the corresponding frequency spectra and the best fit to analytical solutions that can be used to derive cell/nucleus diameter. Once the cell and nucleus diameters are known, the N:C ratio can be calculated using a ratio of the cross-sectional area of the nucleus divided by that of cytoplasm.¹²⁶ Using this approach, the average N:C ratio of MDA-MB-231 cells was 0.54 ± 0.19. This has potential applications in differentiating cell populations and the detection of circulating tumor cells in blood samples.

The cells and nuclei are ideal spherical morphologies with good agreement to theoretical models. It is also possible to employ such techniques to measure other morphologies without spherical symmetry, such as RBCs. RBCs strongly absorb through the UV-near IR spectrum and are homogeneous but with an irregular bi-concave shape. The morphology of RBCs is important and enables their efficient transport throughout the vascular system. RBC shape irregularities, such as swelling, spherical changes, or protrusions, can affect their function.¹³⁰ The RBCs are typically about 7.8 μm in diameter but only 1–2 μm thick. PA signals were measured from RBCs in horizontal and vertical positions relative to the 375 MHz transducer and the 532 nm laser. The PA spectra for the two RBC are very different, and the orientation can be obtained by examining the PA spectra [Fig. 6(i)]. Then, RBCs in a horizontal position were measured as their thickness increased due to swelling. The PA spectrum for a normal, healthy RBC is nearly flat [Fig. 6(j), cyan]; as the RBC swells, the higher frequencies become increasingly attenuated, resulting in an increase in the spectral slope measured at 275–300 MHz. A similar phenomenon occurs with echinocytes, which are spiculated RBCs; the average changes in the spectral slope can be used to evaluate the health of RBCs during blood storage.^131–133 The findings of this work were then extended to measure the quality of blood bags during their 42-day storage period.^132,134 These measurements and simulations show that the PA spectra for RBCs are similar below 100 MHz regardless of their size, shape, or orientation. UHF frequency-based analysis is required to differentiate these parameters at these length scales.

The analysis methods described until now all used a PAM system and manual measurements of each individual cell, which is a long and laborious process. Creating a high throughput method to measure the US and PA signals from single cells would significantly increase the number of cells to enable a more accurate determination of the cellular properties and increase the probability of detecting rare pathophysiology changes. An acoustic flow cytometer was developed by using a polydimethylsiloxane (PDMS) microfluidic device to flow focus single cells through a target area, where the US and PA signals of the passing cells were measured using a 375 MHz US transducer and a 532 nm laser [Fig. 6(k)].^135,136 Using this system, thousands of cells can be measured in a short period of time, enabling a distribution of cell size like the one produced using a Coulter counter. This approach was used to size melanoma cells and dyed AML cells using US and PA¹³⁶ and AML and HT29 cells using US¹²⁹ and also to count the number of RBCs and WBCs in a sample of human blood.¹³⁷ 

Up until this point, we have focused on quantitative analysis of PA signals using spectral analysis techniques of single RF lines. However, frequency-based PA techniques need not be restricted in this way, and it is possible to generate images of biological samples composed of multiple RF lines and spectra. A simple yet powerful technique that demonstrates the capabilities of frequency-based imaging is the frequency band equalization method utilized in wideband optoacoustic mesoscopy systems (10–180 MHz), which bridge PAM and PAT in terms of both imaging depth and resolution.^138–140 The raw acquired PA signals are split into a low band (e.g., 10–60 MHz) and a high band (60–180 MHz) through bandpass filtering. Beamforming is then applied to each of these datasets separately, generating two PA images. The high frequency image is weighted by a factor which is selected to optimize the contrast between the low frequency and high frequency image. A composite image can then be formed by overlaying the weighted high frequency image on the low frequency image. This technique allows for the visualization of features with high-frequency content that are usually masked by higher amplitude low-frequency components originating from larger structures.¹³⁸ 

Figure 7(a) shows how frequency band equalization can be used to image the entire volume of a preclinical mouse tumor and generate high resolution images of the tumor vasculature at different scales [Fig. 7(b)].¹⁴⁰ In general, the PA signals from large vessels have a majority of their energy at lower frequencies and dominate the PA image. By using an appropriate choice of the weighing factor, the signal contribution from the smaller vessels, which have high frequency content,¹¹³ can be boosted so that these structures have equal representation in the reconstructed image, when they would otherwise not be seen.¹³⁸ The practical usefulness of this approach becomes obvious when one examines the effect of the transducer receiving bandwidth on the resolution of PA images of the skin [Fig. 7(c)]. As shown in Fig. 7(d), the capillary loops present in the human skin can be individually resolved when higher frequencies are utilized in the image reconstruction process.¹⁴¹ It is also possible to utilize this approach for monitoring the effects of vascular targeted cancer therapies in preclinical models.^139,140

Another frequency-based PA imaging technique is F-mode imaging, introduced by Moore et al.^142,143 The spectrum is calculated and then divided into several narrow (on the order of 1 MHz) frequency bands. For a PA dataset represented by a 3D array indexed using two spatial coordinates (corresponding to the scan location) and one time coordinate (corresponding to the signal arrival time), this is most easily accomplished by applying a 1D FFT only along the time domain dimension. Taken alone, each frequency bin might only have a small difference in energy compared to adjacent frequency bins, and the fluctuations between bins may seem stochastic. However, if several bins from the same power spectrum are combined into a larger band, then size-dependent spectral features, such as those shown in Fig. 5(b), can be captured. A separate F-mode image is then generated for each spectral band by setting the intensity of a pixel to be proportional to the relative power in the spectral band at that scan location compared to all other locations in the image. As previously discussed, the frequency spectrum of a PA signal for a given absorbing structure will typically contain both spectral maxima and minima at frequencies predetermined by the structure's size and shape. Thus, for a given object, specific spectral bands can be selected that either enhance or diminish the visibility of the object compared to other features in the image. For example, if an F-mode image is generated for a frequency where the spectral signal is strong compared to other regions in the image, then the structure will be readily apparent in the resultant image. Otherwise, if the spectral signal is weak, then the structure will not be observed in the resulting F-mode image. Strategically choosing to generate F-mode images in bands that contain a spectral minimum, such as those shown in Fig. 5, allows for the generation of an image for which the visibility of absorbers of a certain size or shape is negligible.

PA images of blood vessels lend themselves well to F-mode enhancement, as the vessels span a wide range of sizes in the body, and it has been demonstrated that the spectral content of a PA wave emitted by a blood vessel is dependent upon the vessel diameter.¹⁴⁴ Zebrafish larvae, which are ideally suited for UHF PAM imaging due to their small size and high optical transparency, are ideal models for the demonstration of this technique in vivo. Figure 8(a) depicts a schematic rendering of a region within a zebrafish larva that contains the posterior cardinal vein (PCV), intersegmental vessel (ISV), and dorsal aorta (DA), as well as a conventional PA image acquired using a 532 nm laser and a 80 MHz transducer [Fig. 8(b)].¹⁴⁵ In Fig. 8(b), the entire spectral range was used to generate the image, and all three vessels are equally visible. The PA frequency spectra of the three individual vessels as given by the solution to the PA wave equation are shown in Fig. 8(c) and allow for the visualization of the spectral bands that should be selected to optimize contrast between the different vessels in F-mode images. F-mode images created using spectral bands at 12.5, 31.5, and 97.5 MHz are shown in Figs. 8(d)–8(f), respectively. At 12.5 MHz, only the largest vessel, the PCV, is visible. At 31.5 MHz, the PCV and DA are visible. At 97.5 MHz, the PCV is minimized, while the smaller ISV and DA vessels are visible. These results are in good agreement with theoretical models,^4,115,146 which predict that the 22 μm vessel should have the greatest spectral magnitude at lower frequencies. As the frequency is increased, the relative difference between the curves representing the 22 and 15 μm vessels decreases. Finally, at high frequencies, the signal from the smallest vessel should dominate. The F-mode imaging technique has been utilized to analyze blood vessel mimicking phantoms at clinical frequencies and perform high-resolution label-free imaging of single biological cells and their nucleoli.¹⁴³ 

PAT allows for the illumination of a larger field view compared to PAM, providing flexibility in targeting multiple anatomical sites at the expense of resolution.¹⁴⁷ In preclinical imaging, single pulse panoramic photoacoustic computed tomography and spiral volumetric optoacoustic tomography have been developed to achieve real-time, whole-body small animal imaging for tracking circulating tumor cells¹⁴⁸ and the mapping of exogenous contrast agent distribution.²² However, clinically translatable PA devices that can achieve more clinically relevant penetration depths typically consist of a linear array transducer and a limited field of view illumination.⁴¹ Here, one cannot resolve individual absorbers within the illuminated volume but only capture the image arising from the superpositions of signals at multiple biological length scales.¹⁴⁹ This typically occurs when more than 10 photoacoustic sources per acoustic resolution volume (which cannot be individually resolved) are simultaneously illuminated by a laser pulse.¹⁰⁸ Such limited view PAT reconstructions result in the appearance of speckle as evidenced from the acoustic resolution images of tumors,¹⁵⁰ livers,⁷⁴ kidneys,¹⁵¹ intestines,¹⁵² thyroid,⁴⁸ and flowing blood.¹⁵³ 

Although acoustic resolution PA images typically cannot resolve sub-resolution chromophores, the PA signals that are used to reconstruct such images can be analyzed to obtain information that is not readily apparent in the images. Signal analysis techniques can be used to extract information related to the shape/size/orientation of the underlying scatterers (for US) and absorbers (for PA). This is a direct consequence of the physics of ultrasonic scattering¹⁰⁷ and the physics of PA wave generation.⁴ In US imaging, QUS techniques have emerged as signal-based analysis methods for characterizing biological tissue.¹⁵⁴ Given the parallels between US and PA imaging, these methods can be adapted to retrieve physical properties of unresolved optical absorbers in PA as well. As multiwavelength PA imaging can provide functional information (e.g., blood sO₂), it is possible to combine this information with the QUS analysis to provide insight into the tissue microstructure.¹⁵⁵ 

The radial artery in the wrist is easily accessible and plays an important role in understanding the overall systemic circulation.¹⁶⁸ PAT is ideally suited for studying blood's hemodynamic properties at the radial artery due to its superficiality.¹⁶⁹ Algorithms based on the analysis of PA data to quantify the changes in the radial artery sO₂ as a function of flow rate have been developed.^153,169,170 Figure 10(a) shows the sO₂ maps obtained by spectrally unmixing the contributions of oxy- and deoxyhemoglobin at various pulsation rates. The sO₂ variation with time is shown in Fig. 10(b). Using blood flow phantoms, it was demonstrated that the aggregation of RBCs during the pulsatile flow impacts oxygen release to the surrounding medium.¹⁵⁸ As shown in Fig. 10(b), the sO₂ decreases during the peak-systolic phase where the high shear rates in flowing blood disaggregate the RBCs. This is also observed in Fig. 10(c) where the sO₂ is plotted as a function of the number of cells that comprise each RBC aggregate. The number of cells is estimated from the frequency analysis of the simultaneously ultrasound backscatter RF signals. The number of cells comprising the aggregate increases alongside the sO₂ when the flow is in the end-diastolic luminal contraction phase.¹⁷⁰ 

The variation in sO₂ (ΔsO₂) can be examined as a function of the healthy participant's age [Fig. 10(d)].¹⁶⁹ As the oxygen release is a function of RBC aggregation, higher aggregability in older individuals suggests potential impairment in nearby tissue oxygenation. These findings suggest that PAT RF techniques can play an important role in monitoring physiologically relevant parameters that offer insights into underlying pathologies, such as diabetes-induced hyperaggregability.¹⁷¹ 

Given the large optical absorption cross section of RBCs,²⁷ and since biological cells devoid of melanin have low endogenous absorption in the visible optical spectrum, it is thought that the vascular structures filled with RBCs are the dominant sources of the tissue PA signal. The PA analysis techniques described in Sec. II C have been used to develop PA biomarkers of treatment response.^{20,108,134,155,162,172,173} Signal-based biomarkers were correlated to the histopathology of tumors to identify the vascular structures giving rise to the PA signals, and how these change with treatment. Figure 11 provides a snapshot of some of these techniques. For instance, simulations of tightly packed, non-resolvable PA absorbers show that the generalized gamma envelope statistic fit parameter a [Eq. (5)] can quantify changes in effective absorber size and concentration [Fig. 11(a)]. Similar approaches can potentially be used to monitor effective changes in the vessel structure in tumors.

It is possible to combine time and frequency domain techniques to study the functional and structural changes in the tumor vasculature following treatment [Fig. 11(b)]. The drop in tumor sO₂ as early as a few hours post-treatment is accompanied by a decrease in the PA spectral slope.¹⁶² Histology reveals that these effective changes are correlated to the vascular damage that the treatments cause. These changes in spectral parameters can be used to differentiate treatment responders and non-responders as early as a few hours after treatment. These findings have also been extended to other types of cancer treatments such as vascular disrupting agents¹⁷³ and combinations of nanobubbles with radiation.¹⁵⁵ 

It is also possible to utilize the cepstral analysis to study the spatial organization of non-resolvable tumor vasculature [Fig. 11(c)]. Using vascular tree models that mimic the tumor vasculature,¹⁷² it is possible to simulate PAT images representative of the chaotic tumor vasculature. By simulating increased vasculature branching that occurs during tumor angiogenesis,¹⁷⁴ the cepstrum can potentially be used to estimate the (non-resolvable) effective vessel separation. As the number of vessels increases, their spatial separation decreases. Taken as a whole, these functional metrics of tumor metabolism and spatial organization metrics of the tumor vasculature could be useful in quantifying the biophysical changes that occur during cancer treatments. This is particularly important for facilitating the clinical translation of PA imaging so that it can enter mainstream radiology.

Spectral unmixing of PAT signals can also provide useful information with regards to non-resolvable absorbers such as collagen. Collagen plays an important role in fibrosis, a chronic disorder where the buildup of extracellular matrix can cause irreversible damage to kidneys and livers.¹⁷⁵ A collagen unmixing algorithm has been developed for the detection of fibrosis in donated kidneys available for transplant.^151,176 As the demand for kidney transplant increases, physicians are being forced to accept donations from sicker and older donors which contain significantly higher degrees of fibrosis. High levels of donor fibrosis reduces the kidney longevity and it is important to quantify the fibrotic burden at the time of transplant.¹⁷⁷ The signal-based unmixing approach can map the distribution of collagen in both anatomical orientations of whole human kidneys [Fig. 12(a)]. Moreover, 3D scanning allows for quantification of the spatial distribution of collagen and facilitates the comparison with the distribution determined using the gold standard PSR staining [Fig. 12(b)]. In this way, PAT can be used to assess the spatial heterogeneity of fibrosis within human kidneys. This imaging approach can guide the location where biopsies are performed. The accuracy of our unmixing technique is assessed through comparisons with the gold standard histology [Fig. 12(c)].

The same PAT technique can be used for the quantification of liver fibrosis in preclinical animal models. It is possible to quantify the collagen deposition in different levels of liver fibrosis induced in mice [Fig. 12(d)]. The PA spectral unmixing (top row) demonstrates the increasing levels of fibrosis as confirmed in the PSR histology images (bottom row). There is excellent agreement between the PA derived collagen score and the collagen score assessed with the gold standard PSR staining [Fig. 12(e)]. These data suggest that it might be possible to utilize the PAT signal analysis to guide both the progression of preclinical liver fibrosis and to monitor the efficacy of antifibrotic drugs.¹⁷⁸ 

PAM and PAT have traditionally been used for imaging, resolving structures with length scales from micrometers to centimeters. When the structures cannot be resolved, signal analysis techniques can be used to make inferences about the underlying tissue structure. US signal analysis techniques from QUS are largely translatable to PA imaging. This Tutorial has highlighted signal-based analysis techniques to extract structural, functional, and metabolic information from cells and tissues using PAM and PAT. This includes the frequency analysis of the PA signal which can used in PAM to extract information about the effective absorber size and shape and to enhance PAM images by relying on the additional contrast encoded in the PA power spectra. In PAT, traditional QUS techniques such as envelope statistics, power spectrum regression analysis, and cepstral analysis can be used to extract structural and functional information about the tissue microstructure, while spectral unmixing, commonly used in optical imaging, is used to estimate the chromophore concentration and enhance the PA images.

US is an accessible imaging modality commonly used in the clinic. Quantitative ultrasound signal analysis techniques provide information about the unresolvable tissue microstructure and can be used to analyze PA signals. As PA imaging moves from pre-clinical to clinical applications, we expect quantitative analysis techniques to become more commonly used. Since US and PA imaging use complementary hardware, this enables applications that combine the US structural images with PA functional and metabolic images. The additional benefits of using quantitative US and PA signal analysis techniques will make US/PA clinical imaging a highly used tool for diagnostic imaging applications.

E.H., M.J.M., and E.M.S. contributed equally to this work.

The authors would like to acknowledge the financial support of the following funding agencies that have supported various aspects of work through the years: Natural Sciences and Engineering Research Council of Canada, Canadian Institutes for Health Research (CIHR), Canada Research Chairs Program, Canada Foundation for Innovation, Ontario Ministry and Research and Innovation, Canadian Cancer Society, Terry Fox Foundation, Collaborative Health Research Partnerships, Canadian Blood Services, and the Federal Economic Development Agency for Southern Ontario Investing in Commercialization Partnerships Initiative. E. Hysi is supported through CIHR Banting Fellowship and a Kidney Foundation of Canada Fellowship through the Kidney Research Scientist Core Education and National Training Program. The following individuals are also gratefully acknowledged for their contributions to this work: Muhannad Fadhel, Tae-Hoon Bok, Vaskar Gnyawali, Darren Yuen, Elizabeth Berndl, and Xiaolin He.

The data that support the findings of this study are available from the corresponding author upon reasonable request.