Impact of depth-dependent optical attenuation on wavelength selection for spectroscopic photoacoustic imaging.

An optical wavelength selection method based on the stability of the absorption cross-section matrix to improve spectroscopic photoacoustic (sPA) imaging was recently introduced. However, spatially-varying chromophore concentrations cause the wavelength- and depth-dependent variations of the optical fluence, which degrades the accuracy of quantitative sPA imaging. This study introduces a depth-optimized method that determines an optimal wavelength set minimizing an inverse of the multiplication of absorption cross-section matrix and fluence matrix to minimize the errors in concentration estimation. This method assumes that the optical fluence distribution is known or can be attained otherwise. We used a Monte Carlo simulation of light propagation in tissue with various depths and concentrations of deoxy-/oxy-hemoglobin. We quantitatively compared the developed and current approaches, indicating that the choice of wavelength is critical and our approach is effective especially when quantifying deeper imaging targets.

Introduction

Blood oxygen saturation (SO2), the ratio of oxygenated to total hemoglobin concentrations in the blood, provides useful functional information on many pathological conditions including cancer. Since oxygen deficiency (hypoxia) in a tumor is strongly associated with a tumor malignancy, metastasis, and resistance to various therapies, the non-invasive estimation of the SO2 level is essential for tumor characterization and treatment planning [[1], [2], [3]]. Naturally, to explore tumor hypoxia, researchers have investigated various medical imaging and sensing techniques such as needle electrodes, near-infrared spectroscopy (NIRS), positron emission tomography (PET), blood oxygen level-dependent (BOLD) magnetic resonance imaging (MRI), and electron paramagnetic resonance imaging (EPRI) [4,5]. However, none of these techniques offers non-invasive and real-time measurements of SO2 with high sensitivity and fine spatial resolution at clinically relevant depths.

As an effective alternative, photoacoustic (PA) imaging is becoming a promising biomedical imaging modality since it provides high optical absorption contrast from endogenous chromophores and exogenous contrast agents with ultrasonic resolution [[6], [7], [8]]. This ability of PA imaging shows substantial potential in the characterization of tumors, atherosclerotic plaque, and brain functionality [[8], [9], [10], [11], [12]]. Moreover, spectroscopic photoacoustic (sPA) imaging can quantify the concentrations of multiple absorbers. Since each optical absorption spectrum is a unique function of optical wavelength, multiple tissue chromophores can be separated using a multi-wavelength set of PA images [[13], [14], [15]]. Thus, because of its ability to unmix concentrations of multiple absorbers, sPA imaging has been extensively used in many applications including the measurement of SO2 and nanoparticle deposition [11,[16], [17], [18], [19], [20]].

To estimate the concentrations of various absorbers, we have to solve optical and acoustic inverse problems [7,16]. The solution of optical inverse problem requires a multi-wavelength or spectroscopic set of PA images. However, sPA imaging extends data acquisition time and, therefore, degrades temporal resolution and increases motion artifacts. In addition, a tunable laser system is less cost-effective in comparison to a small set of single-wavelength lasers. Thus, the careful choice of wavelengths is imperative. To select an optimal set of wavelengths, a method of capturing the unique spectral features of absorbers under a given number of wavelengths was developed previously [21]. In this approach, we used the stability of the absorption cross-section matrix under a pseudoinverse operation as a criterion and selected a set of wavelengths whose absorption cross-section matrix has the largest minimum singular value as the optimal set. The developed method showed the promising potential in the wavelength selection for in vivo imaging [22].

Because chromophores vary spatially, optical fluence is a function of position and wavelength. Thus, the amount of light fluence reaching the region-of-interest (ROI) needs to be taken into account through analytical or numerical modeling of light transport in tissue. This paper investigates the impact of optical wavelength selection in the presence of light fluence variations in tissue. Under the assumption that the fluence distribution in the tissue is known, we introduce a depth-optimized wavelength selection method that accurately estimates the concentrations of absorbers. In practice, the fluence distribution can be determined by an independent measurement [23,24], or accounted for in an inversion scheme that uses a light transport model [16,[25], [26], [27]]. For verification of our approach, we use a Monte Carlo (MC) simulation of light propagation in tissue with a blood vessel containing various SO2 levels at multiple depths. For a quantitative evaluation, we compare the proposed and previous methods by applying the various levels of additive white Gaussian noise and computing the root-mean-square errors (RMSE) of oxy- and deoxy- hemoglobin concentrations.

Theory of spectroscopic photoacoustic (sPA) imaging

The photoacoustic effect describes the generation of sound waves following light absorption and thermoelastic expansion, resulting in a peak local pressure rise, p0. Here, p0 is a function of optical wavelength λ and determined by the combination of three quantities (Grüneisen coefficient Γ (the PA efficiency), optical absorption coefficient μ, and light fluence Φ) as follows [16]:Three wavelength-dependent quantities and the Grüneisen coefficient in Eq. (1) (i.e., p0, μ, Γ and Φ) are all also spatially variant, but for the simplicity, their position vector is not used in this paper. Here, μ is generally a linear combination of absorber concentrations, and thus, can be represented aswhere I is the number of absorbers of interest, ε is the molar absorption cross section of the ith absorber, and C is the concentration of the ith absorber. We assume that the Grüneisen coefficient is constant and the light fluence distribution is known. In addition, although the reconstruction of the acoustic signals received by an ultrasound transducer, depending on tissue heterogeneity, acoustic properties, and detector geometry, is not trivial, here we assume that p0 can be accurately reconstructed [28]. Then, the local ultrasound pressure (p0 in Eq. (1)) can be converted to μ by assuming that the light fluence distribution over wavelength at every spatial location is known.Then, by combining Eqs. (2) and (3) and using vector notation, a set of linear equations is formulated aswhere μ is the M-by-1 vector with reconstructed absorption coefficients (M is the number of wavelengths used in sPA imaging), ε is the M-by-I absorption cross-section matrix, and C is the I-by-1 vector of absorber concentrations (I is again the number of optical absorbers of interest). Then, when M is equal to or greater than I, the concentrations can be computed bywhere C is the estimates of true concentrations (C) and ε is the pseudoinverse of ε. Here, to prevent concentration estimates from becoming less than 0 or larger than 1, further constraints, including a non-negativity constraint, can be applied in Eq. (5) [17,29].

Materials and methods

A depth-optimized wavelength selection method

Although the use of extra wavelengths typically allows us to obtain more precise estimates of concentrations in the presence of noise (but, in the absence of motion), it is desirable to keep the number of wavelengths to as few as possible to avoid motion artifacts, to improve temporal resolution, and to accelerate image acquisitions. To reduce the redundancy in the choice of wavelengths while preserving accuracy, this study introduces an algorithm that finds an optimal set of wavelengths for accurate concentration estimation. Then, with the consideration of the thermal and electronic noise in the system, local acoustic pressure p0, in Eq. (1) can be redefined by adding white Gaussian random noise n:Dividing both sides of Eq. (6) by the Grüneisen coefficient and the optical fluence yieldsBy presenting Eq. (7) as a simpler vector form, we obtain the following equation:where μ is the M-by-1 vector of the reconstructed absorption coefficients, n is the M-by-1 noise vector, and Φ is the M-by-M diagonal matrix where Φ(j,j) = Φ(j). Like Eq. (5), the noise-corrupted estimates of concentrations C can be obtained by taking the pseudoinverse of ε.

From Eq. (9), C becomes C if and only if the noise-related second term becomes 0. Therefore, to remove errors in all concentrations, the following equation needs to be minimized:Then, since n is a random variable (zero-mean Gaussian noise), Eq. (10) can be modified aswhere 1 is an M-by-1 vector containing all 1′s. Our goal is to find a set of wavelengths that minimizes the condition (11) under the given number of wavelengths (M). In other words, a set of wavelengths minimizing Eq. (11) (i.e., the noise-related second term in Eq. (9)) makes C as close as to C (i.e., the noise-free ideal concentration estimates). In addition, it is clear that the Eq. (11) relies on two quantities: ε (i.e., the pseudoinverse of absorption cross-section matrix) and (i.e., the inverse of diagonal fluence matrix). Given the absorbers of interest, the absorption cross-section matrix can be assumed to be known, but the fluence matrix is generally unknown. Thus, in order to solve the Eq. (11), the fluence matrix should be modeled in practice. In this study, a Monte Carlo model was used to model the fluence as described below in Section 3.2. In the minimization process, because the wavelengths typically used in sPA imaging are finite and discrete (wavelengths from 680 nm to 970 nm with a one nm step size are used in this study), every possible set of wavelengths was applied to find an optimal set of wavelengths where the objective function (11) is globally minimized. Here, note that in (11) can be rewritten as and the term is maximized in the method to minimize the noise-related error in solving a spectral unmixing problem. Here, note that after a set of optimal wavelengths is selected by this approach, any spectral unmixing methods can be applied to obtain the concentrations of chromophores.

Monte Carlo study parameters

To investigate our approach under various conditions, the recently introduced MC program (mcxyz.c) [28] was used, which allows us to obtain precise optical fluence maps, including a blood vessel with various SO2 levels. Fig. 1 shows tissue models that consist of air, epidermis (600-μm-thickness), dermis (with a blood volume fraction of 0.2%), and blood. The blood vessel (500-μm-diameter) was placed at three depths (0.4 cm, 1.0 cm, and 1.6 cm) with five SO2 levels (0.00, 0.25, 0.50, 0.75, and 1.00). Here, the SO2 level of blood in dermis is the same as that in the blood vessel. The laser beam applied was uniform and had a radius of 0.4 cm. Fig. 2 shows the absorption coefficients of oxy-, deoxy-hemoglobin, water, and melanin used in our study. For oxy- and deoxy-hemoglobin, typical values of the concentration ratio of hemoglobin in blood (150 g/L) and molecular weight of hemoglobin (64,500 g/mol) were used. The software and the detailed optical parameters of the tissue model, including anisotropy, absorption, and scattering coefficients, were obtained from the OMLC website (please refer to makeTissueList.m and spectralLIB.mat files at https://omlc.org/software/mc/mcxyz/index.html). Tissue parameters for epidermis, dermis, blood, air, and water used in this study were exactly the same as those listed in the website while SO2 levels were varied.

Fig. 1

Tissue models with air, epidermis, dermis, and 500-μm-diameter blood vessels at three depths: (a) 0.4 cm, (b) 1.0 cm, and (c) 1.6 cm.

Fig. 2

Absorption coefficients of oxy- and deoxy-hemoglobin, water, and melanin (data obtained from the spectralLIB.mat file available at https://omlc.org/software/mc/mcxyz/index.html).

To collect the light fluence from the program, a three-dimensional array with 800 × 640 × 2 bins, each 25 μm was created. The optical wavelengths of light ranged from 680 nm to 970 nm with a one nm step size (thus, 291 wavelengths). Therefore, 291 × 5x3 light fluence maps corresponding to the number of wavelengths, SO2 levels, and vessel positions, respectively were acquired. Each simulation time of the light fluence of an individual wavelength was 90 min (resulting in an average number of photons of 1.5 × 106), so this study required a 291 × 5x3 × 1.5-hour simulation. To expedite the simulation, 15 cores of an Intel Xeon processor (Intel Corp., Santa Clara, CA) were utilized simultaneously, allowing us to simulate all fluence maps in 291 × 1.5 h. Overall, the MC program was used to acquire various optical fluence maps with multiple wavelengths (680 nm to 970 nm), five SO2 levels (0.00, 0.25, 0.50, 0.75, and 1.00), three depths of a blood vessel (0.4 cm, 1.0 cm, and 1.6 cm), which took approximately 18 days.

Data handling

After obtaining optical fluence maps, the fluence signals in each blood vessel at each individual wavelength were first averaged and one-dimensional optical fluence spectra with respect to the wavelength was constructed. Then, the averaged optical fluence signals were multiplied by the known absorption coefficients, which reconstructs the PA spectrum. To build the noise-added PA spectrum as in Eq. (6), various amounts of zero mean white Gaussian random noise were generated and added to the reconstructed PA spectrum. In this study, we used the simple noise model by assuming the frequency response of the transducer is ideally flat as typically adopted in other studies for brevity [30,31]. Here, the standard deviations (STD) of noise ranged from 0.0 to 0.1, which results in a reasonably wide range of the signal-to-noise ratio (SNR) in the PA spectrum (about 10.7 dB to 85.1 dB. See Fig. 7). In addition, to obtain statistical results, adding this random noise onto the PA spectrum and estimating concentrations were repeated 104 times.

Fig. 7

Mean SNRs of PA spectrum in blood vessels at three depths (0.4 cm, 1.0 cm, and 1.6 cm) as a function of noise STD.

Assessment

For the quantitative evaluation, this study used the root-mean-squared error (RMSE) of the concentration estimates (i.e., the concentrations of oxy-hemoglobin (HbO2), deoxy-hemoglobin (Hb)) as follows:where C represents a concentration estimate, C is a ground-truth concentration set in the MC simulation, and T is the number of noise trials (here, T=104). Note that C is any element of C in Eq. (9). The RMSE of two estimates: the concentrations of HbO2 and Hb used to quantify SO2 at various vessel depths and SO2 levels as a function of the noise standard deviation were computed. To investigate the robustness of our method under the various amounts of noise, the white Gaussian random noise was repeatedly applied.

In this study, M was set to two (i.e., the minimum number of wavelengths to identify the distribution of two chromophores). Given this condition, we compared two wavelength selection methods: the proposed noise minimization approach and the previous minimum singular value-based approach [21]. (here, these methods are referred to as nmin, and σmin, respectively.) In addition, we applied two methods to reconstruct concentrations: linear least squares with non-negativity constraints (the lsqnonneg function in MATLAB) and without any constraints (referred to as CLLS and LLS, respectively). Although two linear least squares methods were used in this paper, both wavelength selection methods are not restricted to any spectral unmixing method. In summary, two wavelength selection approaches (nmin and σmin) with and without non-negativity constraints (thus, total four methods) were used to estimate the concentrations of HbO2 and Hb at three depths of a vessel and five SO2 levels as a function of noise levels. For quantitative comparisons, the RMSE of the estimates and the average of corresponding SO2 were computed.

Results and discussion

Fig. 3, Fig. 4, Fig. 5, show the optical fluence maps containing a blood vessel at three depths: 0.4 cm, 1.0 cm, and 1.6 cm. All the fluence maps (in the unit of J/cm2 per J delivered or 1/cm2) are displayed on a logarithmic scale. The rows of each figure represent the fluence maps of three representative optical wavelengths (680 nm, 875 nm, and 970 nm) selected from 291 maps, which range from 680 nm to 970 nm with a 1-nm step size. The columns represent five SO2 levels (0.00, 0.25, 0.50, 0.75, and 1.00) in blood vessels. Fig. 3, Fig. 4, Fig. 5 show that the optical fluence varies according to the optical wavelength and the SO2 level. These variations in the fluence indicate that the choice of wavelengths affects the strength of the fluence, which determines the SNR. Furthermore, because of light attenuation, the deeper the vessel is located in tissue, the more the overall fluence decreases. Therefore, the optical wavelength selection becomes more challenging and crucial, particularly when imaging deeply located targets.

Fig. 3

Representative optical fluence maps of the tissue model containing a blood vessel at a 0.4 cm depth. Rows represent wavelengths of light, 680, 875, and 970 nm, and columns represent five levels of oxygen saturation in the blood vessel, 0.00, 0.25, 0.50, 0.75, and 1.00. The unit of fluence Φ is J/m2 per J delivered or 1/m2.

Fig. 4

Representative optical fluence maps of the tissue model containing a blood vessel at a 1.0 cm depth. Rows represent wavelengths of light, 680, 875, and 970 nm, and columns represent five levels of oxygen saturation in the blood vessel, 0.00, 0.25, 0.50, 0.75, and 1.00. The unit of fluence Φ is J/m2 per J delivered or 1/m2.

Fig. 5

Representative optical fluence maps of the tissue model containing a blood vessel at a 1.6 cm depth. Rows represent wavelengths of light, 680, 875, and 970 nm, and columns represent five levels of oxygen saturation in the blood vessel, 0.00, 0.25, 0.50, 0.75, and 1.00. The unit of fluence Φ is J/m2 per J delivered or 1/m2.

To investigate the wavelength- and depth-dependent characteristics of optical fluence, we first averaged the fluence signals of each individual wavelength within each blood vessel as one-dimensional optical fluence functions with respect to wavelengths (Fig. 6). According to optical wavelength, SO2, and depth, the fluence signals vary substantially. As a notable example, the fluence patterns of 0.4-cm and 1.6-cm vessels clearly differ. In contrast to the fluence at the 0.4-cm vessel, the fluence at the 1.6-cm vessel abruptly decreases and becomes negligible as the wavelength increases. This observation suggests that avoiding using higher wavelengths (> 950 nm) when estimating SO2 in the far region is desirable.

Fig. 6

One-dimensional averaged optical fluence signals in three vessel depths at (a) 0.4 cm, (b) 1.0 cm, and (c) 1.6 cm at 0 cm in the lateral direction as a function of wavelength.

To construct the noise-added PA spectra in Eq. (6) and to investigate the robustness of our approach, we added various amounts of white Gaussian random noise to PA spectra reconstructed from the averaged 1D fluence signals and the known absorption coefficients. The standard deviations (STD) of noise spanning from 0 to 0.1 created the wide range of SNR. We repeated this process 104 times to find the mean SNRs of the PA spectra. Fig. 7 shows the mean SNRs over the noise STD. As the noise STD increases, the mean SNRs decrease. In addition, three mean SNR functions are plotted according to three depths of a vessel. The deeper the vessel, the more SNR decreases. After all, we use these noise-added PA spectra to assess two wavelength selection methods (nmin and σmin).

To estimate two concentrations of HbO2 and Hb to quantify SO2 levels, this study uses two optical wavelengths (M = 2, the minimum number of wavelengths to solve the optical inverse problem.) selected by the proposed and the previous methods: nmin and σmin. nmin finds two wavelengths that minimize the added noise term whereas σmin finds two wavelengths that maximize the minimum singular value of the absorption cross-section matrix. Therefore, the former minimizes the noise-related errors, and the latter allows the selected absorption cross-section matrix to have unique spectral features. As a result, σmin selects two wavelengths: 680 nm and 970 nm. Here, since σmin relies on only the stability of absorption cross-section matrix, we use the same set of wavelengths, regardless of the depth or the concentration of chromophores. However, since nmin exploits the optical fluence distribution and the optical fluence varies according to the imaging depth or the concentration of absorbers, our approach adaptively offers optimal sets of wavelengths as presented in Table 1.

Table 1

Two (M = 2) selected wavelengths (nm) used for estimating SO2 by the proposed approach (nmin) at three Vessel depths under five SO2 levels.

	Vessel depth
SO2	0.4 cm	1.0 cm	1.6 cm
0.00	692, 970	688, 895	681, 882
0.25	694, 969	683, 887	687, 887
0.50	683, 970	684, 909	693, 881
0.75	681, 970	680, 906	684, 878
1.00	681, 909	682, 887	680, 893

Table 1 describes the selected sets of wavelengths via the proposed method (nmin). First, where the vessel depth is 0.4 cm, the PA spectrum has relatively high SNR as shown in Fig. 7, and the optical fluence hardly fluctuates or decreases over the wavelength as shown in Fig. 6. Therefore, most of the sets of wavelengths (first column in Table 1) are similar to those (680 nm and 970 nm) chosen by σmin. Nevertheless, small shifts in the wavelength are observable. As SO2 decreases from 1.00 to 0.00, the optical fluence at 680 nm decreases. Hence, nmin tends to select the wavelength with higher optical fluence to satisfy the Eq. (11). As a result, the first wavelength increases from 681 nm to 692 nm as SO2 changes from 1.00 to 0.00. The second wavelengths used for the 0.4-cm vessel mostly remain at around 970 nm, except when SO2 is 1.00. Since the optical fluence at 970 nm decreases gradually when SO2 increases. Thus, when SO2 is 1.00, the wavelength moves to 909 nm, which indicates that the gain in the optical fluence starts to outweigh the loss in spectral distinction. Similarly, as the depth of a vessel increases, nmin selects shorter wavelengths as the second wavelengths.

Fig. 8, Fig. 9 show the RMSE of estimated HbO2 and Hb concentrations under various SO2 levels and vessel depths. This study uses two wavelength selection methods (nmin and σmin) with two constraints (with and without the non-negativity constraint) (i.e., four methods: nmin+LLS, nmin+CLLS, σmin+LLS, and σmin+CLLS). The RMSE is negligible in both HbO2 and Hb at the 0.4-cm depth vessel, since as illustrated in Fig. 7, the overall SNR is sufficiently large at 0.4 cm over all wavelengths (>40 dB). Furthermore, since the choices of wavelengths by two wavelength selection methods (nmin and σmin) are almost identical when the vessel is located at 0.4 cm, the RMSE differences between the two methods are minor, except when SO2 is 1.00. However, the deeper the vessel is located, the more the overall SNR decreases, and thus, the RMSE becomes larger. In addition, from the 1.6-cm vessel, we can observe the largest gap between the RSME of the two wavelength selection methods, which suggests that the wavelength choice is critical in far-depth imaging. Commonly, the non-negativity constraint (CLLS) produces more precise estimates in terms of RMSE.

Fig. 8

RMSE of estimated HbO2 concentrations from four methods: σmin+LLS, σmin+CLLS, nmin+LLS, and nmin+CLLS with respect to the noise STD. Rows represent three vessel depths (0.4 cm, 1.0 cm, and 1.6 cm). Columns represent five levels of SO2 (0.00, 0.25, 0.50, 0.75, and 1.00).

Fig. 9

RMSE of estimated Hb concentrations from four methods: σmin+LLS, σmin+CLLS, nmin+LLS, and nmin+CLLS with respect to the noise STD. Rows represent three vessel depths (0.4 cm, 1.0 cm, and 1.6 cm). Columns represent five levels of SO2 (0.00, 0.25, 0.50, 0.75, and 1.00).

In Fig. 10, we compute mean SO2 levels from the concentration estimates of HbO2 an Hb. Here, since any negative estimate in the concentration of HbO2 or Hb causes amplified or negative SO2 values, the non-negativity constraint is desirable in the estimation of SO2. However, the non-negativity constraint results in a bias in SO2 when SO2 is 0.00 or 1.00. Overall, as the vessel depth increases, the accuracy in mean SO2 levels degrades, shown in Fig. 10. Particularly, when the noise STD is larger than 0.01 at the 1.6-cm vessel (that is, SNR < 30 dB), the σmin method without the non-negativity constraint (i.e., σmin+LLS) fails to reconstruct reliable SO2 levels. The non-negativity constraint partially improves accuracy but still yields substantial errors. However, nmin stably measures SO2 levels regardless of vessel depths, constraints, and concentrations.

Fig. 10

Estimated mean SO2 by four methods: σmin+LLS, σmin+CLLS, nmin+LLS, and nmin+CLLS with respect to the noise STD. Rows represent three vessel depths (0.4 cm, 1.0 cm, and 1.6 cm). Columns represent five levels of SO2 (0.00, 0.25, 0.50, 0.75, and 1.00).

We investigated the effect of the choice of optical wavelengths in the estimation of SO2 using sPA imaging in the presence of wavelength-dependent, spatially-variant fluence through the simulation. Different amount of additive white random noise was applied to demonstrate the changes in desired optical wavelengths, showing that the current method outperforms the previous method especially when the SNR is low (i.e., the noise STD is high). In addition, the wavelength choices in Table 1 indicate that for deeper vessel (i.e., the SNR becomes lower), the wavelengths selected with our method are increasingly different from those estimated using the previous method, resulting in improved accuracy in SO2 estimation (note that the wavelength sets in Table 1 are only optimal in our simulation settings, not optimal in general). Thus, for sPA images with low signal and/or high noise, our method finds a more optimal set of optical wavelengths. However, the noise model used in our study does not reflect all noise components present in sPA imaging. Thus, a future study should include other noise sources, such as thermal acoustic noise from the medium and thermal noise from the transducer, to make the model more realistic [32].

The results clearly show that the imaging depth and tissue properties play a significant role in the optimal wavelength selection for sPA imaging. As researchers seek to reduce the cost and/or improve the temporal resolution of sPA imaging systems, using only small subset of the optical spectrum is an appealing option. In practice, the specific choice of imaging wavelengths is nontrivial and application-dependent. In this work, we assumed that the optical fluence distribution is known. In practice, because spatially-varying optical properties of tissue are unknown, successful modeling of the optical fluence distribution is often challenging and requires several assumptions, such as the number and a type of absorbers and the scattering coefficients. Therefore, the selected wavelengths found by our method likely to be suboptimal. These results indicate that extensive analysis of a variety of tissue properties must be performed prior to selecting a subset of the available wavelengths. One additional consideration is the maximum safe level of laser exposure [33]. In the wavelength range from 680 nm to 970 nm, longer wavelengths allow us to use higher laser energy and achieve higher SNR, which could result more accurate sPA images.

This study assumed that the acoustic inversion problem to reconstruct the local PA spectrum from the received PA spectrum can be solved. However, to reconstruct acoustic pressure generated in local absorbers accurately, we should consider heterogeneous tissue, bandlimited ultrasound detectors, and other physical limitations [34]. In reality, the heterogeneity of tissue may cause variations in the speed of sound, resulting in phase aberration [35]. This degrades the image quality as well as the accuracy in concentration estimation. Thus, the accurate modeling of non-homogenous tissue properties should be considered when reconstructing PA images. In addition, because the received PA signals are inherently wide band and noisy, the sensitivity of the ultrasound transducer is crucial. The limited bandwidth and SNR of the PA spectrum affect the resolution and the accuracy of SO2 estimation [36]. Furthermore, because of frequency-dependent attenuation of acoustic waves, as the imaging depth increases, the center frequency of the acoustic wave downshifts and its bandwidth becomes narrow, degrading the resolution and the SNR. Moreover, using data acquired with a linear array transducer, the limited-view tomographic reconstruction leads to the artifacts in the local PA reconstruction. Finally, the reverberation artifact (or reflection artifact), mainly caused by strong PA signals generated from superficial optical absorbers, is another issue in sPA imaging because it significantly limits the image contrast and the imaging depth [37].

Conclusion

Through the simulation study, we have shown that the wavelength selection for the SO2 estimation in sPA imaging is affected by the corresponding wavelength-dependent optical fluence distribution. To improve the SO2 estimation capability, this paper has introduced the depth-optimized wavelength selection method, which allows us to obtain the concentration estimates with higher accuracy. Assuming that the wavelength-dependent, spatially-variant optical fluence distribution can be reliably obtained, this approach finds a set of wavelengths, which minimizes the noise-related errors in the concentration estimates. Through the MC simulation of light propagation in tissue, we have investigated the effect of various imaging parameters, including the depth and the concentrations of chromophores. The quantitative assessment using RMSE showed that the proposed approach robustly and accurately estimated SO2 over the previous method. This study indicates that the choice of wavelength in sPA imaging is critical, and our depth-optimized method could be effective, especially when quantifying deeper agents.

Conflict of interest statement

The authors declare no competing financial interests.