Dual-wavelength photoacoustic atlas method to estimate fractional methylene blue and hemoglobin contents

Abstract. Significance Methylene blue (MB) is an exogenous contrast agent that has the potential to assist with visualization and penetration challenges in photoacoustic imaging. However, monitoring the local concentration between MB and endogenous chromophores is critical for avoiding unnecessary MB accumulations that could lead to adverse effects such as hemolysis when exposed to increased dose and photodamage when exposed to high laser energies. Aim We developed a modified version of a previously proposed acoustic-based atlas method to estimate concentration levels from a mixture of two photoacoustic-sensitive materials after two laser wavelength emissions. Approach Photoacoustic data were acquired from mixtures of 100-μM MB and either human or porcine blood (Hb) injected in a plastisol phantom, using laser wavelengths of 710 and 870 nm. An algorithm to perform linear regression of the acoustic frequency response from an atlas composed of pure concentrations was designed to assess the concentration levels from photoacoustic samples obtained from 11 known MB/Hb volume mixtures. The mean absolute error (MAE), coefficient of determination (i.e., R2), and Spearman’s correlation coefficient (i.e., ρ) between the estimated results and ground-truth labels were calculated to assess the algorithm performance, linearity, and monotonicity, respectively. Results The overall MAE, R2, and ρ were 12.68%, 0.80, and 0.89, respectively, for the human Hb dataset and 9.92%, 0.86, and 0.93, respectively, for the porcine Hb dataset. In addition, a similarly linear relationship was observed between the acoustic frequency response at 2.3 MHz and 870-nm laser wavelength and the ground-truth concentrations, with R2 and |ρ| values of 0.76 and 0.88, respectively. Conclusions Contrast agent concentration monitoring is feasible with the proposed approach. The potential for minimal data acquisition times with only two wavelength emissions is advantageous toward real-time implementation in the operating room.

Depending on the dose, exposure time, and number of target cells, adverse effects due to misuse of contrast agents include acute inflammation, 22 apoptosis, 22,23 necrosis, 14 cellular toxicity, 24 allergies, 25 reduction in cellular viability, 26 nephropathy, 27 hemolysis, 28 and photodamage. 29 In addition, maximum absorption and clearance times (i.e., the time for the drug to completely leave the system) range from 5 min to 24 h among the drug delivery approaches reported in the literature, 20,[30][31][32] requiring different dose, delivery, and monitoring protocols.
Currently, most approaches to measuring the concentration levels of chromophores rely on the acquisition of photoacoustic responses from multiple laser wavelength emissions paired with spectral unmixing methods. [33][34][35][36][37][38][39][40] However, as real-time implementations are necessary for monitoring chromophores with short clearance times in applications such as photoacoustic-guided surgery, [41][42][43][44] traditional techniques are typically not feasible because of the lengthy acquisition times associated with transmitting multiple laser wavelengths to achieve a single estimate. 45 In addition, traditional spectral unmixing techniques do not typically consider differences in acoustic spectra, which has the potential to provide additional information for differentiation of exogenous and endogenous chromophores.
We previously proposed an acoustic frequency-based method to discriminate photoacoustic responses from different materials and overcome challenges with traditional spectral unmixing techniques. 46 By measuring the photoacoustic response from only two laser wavelengths, the initial version of our dual-wavelength atlas method achieved comparable sensitivity, specificity, and accuracy to traditional spectral unmixing methods 47,48 and related classifiers. 49,50 Based on the linear relationship between contrast agent concentration and photoacoustic amplitude, [51][52][53] we hypothesize that our method can be extended to locally characterize the volumetric ratio between two photoacoustic-sensitive materials. Specifically, we focus on the exogenous chromophore MB and the endogenous chromophore hemoglobin (Hb), motivated by the potential for photoacoustic-based catheter interventions performed with an optical fiber housed within a cardiac catheter inserted through a major vein. 54 MB would be administered through the catheter and eventually mix with Hb to enhance the visualization of target structures in vascular pathologies (e.g., atherosclerositic plaques, 55 thrombi, 56 tumor neovasculatures, 57,58 and endothelia 59 ), and the local concentration of MB versus Hb would be monitored with our proposed method with either external 54 or intravascular 60 ultrasound sensor placement. This monitoring will enable real-time interventions to avoid the adverse effects produced by the unnecessary accumulation of MB. [61][62][63][64] In addition, the MB used in the present study serves as a surrogate for other potential intravascular contrast agents that would similarly benefit from local concentration monitoring to prevent cell toxicity. 24,60 The remainder of this paper is organized as follows. Section 2 details the phantom experimental setup, a novel dual-wavelength atlas method to estimate MB and Hb concentrations, and quantitative metrics for performance evaluation. Section 3 presents the quantitative evaluation of mixture estimation performance, as well as the evaluation of linear and monotonic trends of estimated concentration versus ground-truth labels, including an assessment of appropriate region sizes to perform the proposed estimation. Section 4 discusses insights from these results, and Sec. 5 summarizes the clinical impact of the proposed methods.
among chambers due to different air bubble distributions, a single chamber was filled with either a 100-μM aqueous solution of MB from Fisher Scientific (Waltham, Massachusetts), blood (Hb) or a combination of MB and Hb. A 1-mm-diameter optical fiber was inserted in the filled chamber, and the fiber tip was positioned ∼20 mm below the top surface. The optical fiber was connected to a Phocus Mobile laser (Opotek Inc., Carlsbad, California), transmitting laser light with wavelengths of 710 and 870 nm and with a laser energy of 3 mJ. The resulting photoacoustic signals were received by an Alpinion L3-8 linear array ultrasound probe (Alpinion Medical Systems, Seoul, South Korea) positioned on the lateral wall of the phantom, ∼20 mm from the hollow chamber cross section. This experimental setup, which is shown in Fig. 1(a), has been described in previous publications. 46,65 Two types of Hb samples were used in this study. First, 23 vials of fresh human Hb were obtained up to two days after blood draw and storage, mixed in a single container, and exposed to air for ∼1 h to minimize the variability of Hb oxygenation levels. Second, experimental whole porcine blood (Innovative Research, Novi, Michigan) was used on the fourth day of its reported 24-day lifetime.
A total of 11 concentration levels were prepared using different 2-mL mixtures of Hb and MB. These concentrations ranged from 0% to 100% in 10% increments of MB volume percentage. For example, for the 2-mL MB-Hb mixture, a concentration of 60% MB consisted of 1.2 mL of 100-μM MB and 0.8 ml of Hb. Each concentration mixture was manually stirred with a syringe in a separate container. Unless otherwise stated, the % concentrations reported in this manuscript refer to the % MB concentration. Five trials were conducted for each concentration level by fixing the light-delivering optical fiber at 0 deg, with 20 photoacoustic frames (10 per wavelength) acquired per trial using a frame rate of 5 Hz. It is worth noting that two blocks of 64channel aperture data were required to obtain an image created from 128 receive elements, which reduced the possible frame rate from 10 Hz (based on the 10-Hz laser pulse repetition frequency) to 5 Hz.

Dual-Wavelength Atlas Method for Mixture Estimation
The framework for mixture estimation is shown in Fig. 1(b), based on the dual-wavelength atlas method previously proposed for the identification of individual chromophores. 46 With this setup, frequency domain information is expected to be useful because the optical fiber is placed inside the PVCP chamber, thus generating fluence maps that diminish radially from the tip of the optical fiber. This fluence distribution generates unequal acoustic frequency responses, as volume regions of different sizes are being excited. Thus, the proposed algorithm is anticipated to leverage frequency domain information to differentiate photoacoustic responses from various concentrations of chromophores.
To implement the proposed approach, conventional delay-and-sum (DAS) beamforming was first employed to create photoacoustic images for each wavelength emission, concentration, and trial. In contrast to preceding work, 66 we used M-weighted short-lag spatial coherence (SLSC) 67 with a cumulative lag M ¼ 20 instead of conventional SLSC with M ¼ 5 to generate the binary masks that segmented signals of interest and provided ground-truth labels. This specific Mweighted SLSC cumulative lag value (i.e., M ¼ 20) was chosen based on the comparison of SLSC 68 images with M-weighted SLSC. First, a representative SLSC 68 photoacoustic image from the blood dataset was generated with M ¼ 5. The same photoacoustic frame was then processed with M-weighted SLSC with M values ranging from 10 to 30. Binary masks generated with a −3-dB threshold were then created for each image, and the similarity was defined with the Dice coefficient: 69 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 6 2 8 where AND is the logical "AND" operator between the binary masks. This process was then repeated for one trial of the human Hb dataset, containing 11 concentrations, 2 wavelengths, and 10 frames. The Dice coefficients calculated between SLSC and M-weighted SLSC masks were then plotted as a function of M, resulting in M ¼ 20 yielding the highest mean similarity to the SLSC masks with the least standard deviation. The modification from SLSC in our preceding work 66 to M-weighted SLSC in this work was implemented to increase the influence of lower lags over higher lags, thus creating less disjointed binary masks. A single binary mask was obtained per trial, resulting from the logical inclusive "OR" operation of 20 masks (obtained from 10 frames per two laser wavelengths).
As noted in Fig. 1(b) and in our previous work, 66 for each image, only those pixels included in the coherence masks were used for feature extraction, training, and classification. In-phase and quadrature (IQ) demodulation was implemented with 2.75 MHz and 85% bandwidth, and the resulting analytic spectra from a sliding window of axial kernels was stacked along the frequency axis. For the last step of feature extraction, principal component analysis (PCA) was applied to the power of the stacked spectra to reduce the complexity of the feature space.
The training dataset consisted of a random trial of 0% MB (i.e., 100% Hb) and 100% MB, and the remaining trials and concentration levels were included in the testing datasets. After the dual-wavelength atlas was constructed [i.e., PCA MB and PCA Hb in Fig. 1(b)], a mixture estimation algorithm was designed to obtain the concentration level of MB from a testing sample. First, N ¼ 1000 samples were randomly selected from the MB and Hb atlas. Then, assuming that mixtures are linear combinations of pure MB and Hb concentrations, a concentration distribution C 0 ∈ R N×1 was obtained using the following equation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 3 1 7 where PCA x is the projected spectra of x% testing concentration, k:k 1 is the Manhattan norm, and each PCA matrix is of size R N×p with p principal components. A histogram filter was then applied to the C 0 vector, where the bin with the highest number of individuals out of 10 bins was extracted. The mean value C 00 was computed from the reduced C 0 vector, and the process was repeated for each valid pixel and frame of the trial. Finally, a median filter was applied to the concentration maps with a kernel of 0.64 × 0.70 mm × 30 in the axial, lateral, and frame dimensions, respectively. The 0.64-mm axial kernel length was selected to reside within the frequency response of the received photoacoustic signals, which ranged from 2 to 2.75 MHz, corresponding to wavelengths ranging from 0.77 to 0.56 mm, respectively. As the lateral and axial resolutions of an image differ, the lateral kernel width (i.e., 0.70 mm) was chosen to be the closest possible match to the axial kernel length. The frame dimension of the kernel reduced the concentration tensor from C 00 ∈ R N i ×N j ×N f to C ∈ R N i ×N j , where N i and N j are the number of rows and columns, respectively, of the reconstructed photoacoustic image. Because Eq. (2) is applied on frequency data, the three most influential parameters of the dual-wavelength atlas method for the estimation of C are the axial kernel size, threshold of spectral log compression, and number of principal components. An increase in axial kernel size improves the accuracy of the estimation of spectral amplitudes in the Fourier space.
The threshold of spectral log compression adds a tolerance to the discrimination between spectral peaks, which affects the variance that is computed in the PCA. The number of principal components controls the amount of information that is reduced to the feature space.
To provide complementary information regarding the location of photoacoustic targets in the imaging plane, the displayed concentration maps were overlayed with an ultrasound image of the plastisol chamber processed with locally weighted SLSC (LW-SLSC), 70 using a regularization factor α ¼ 1 and a 1.25 mm × 1.2 mm kernel. The LW-SLSC parameters were obtained from our previous optimization study. 66 A factor of α ¼ 1 corresponded to an equal weight between cost and penalty function for calculating the optimal coefficients in LW-SLSC. An axial kernel size of 1.25 corresponded to approximately 4.5λ, where λ is the wavelength associated with the center frequency of the L3-8 ultrasound probe.

Quantitative Metrics
To evaluate the overall accuracy of the dual-wavelength atlas method for mixture estimation, the estimated concentration levels C were first ordered as a function of the ground-truth labels k. The function f 1 ∶CðkÞ was then compared with f 2 ∶C ¼ k, representing a 1:1 relationship between estimated and true concentrations. Finally, the coefficient of determination R 2 was used to quantify the performance of the dual-wavelength atlas method: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 5 0 6 where C ki is the estimated concentration of measurement i for a ground-truth concentration k and C is the mean estimated concentration. In addition, the R 2 metric was used to assess the linearity of the photoacoustic spectra obtained from mixtures of MB and Hb by implementing a linear regression of spectral amplitudes as a function of the ground-truth concentrations.
Similarly, monotonicity was evaluated in two cases: (1) between photoacoustic spectral amplitudes of a specific frequency and ground-truth concentrations and (2) between estimated and ground-truth concentrations. The monotonic trend of these data pairs were quantified with Spearman's rank correlation coefficient ρ, described in the following equation: 71 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 3 5 6 ρ ¼ covðRðXÞ; RðYÞÞ σ RðXÞ σ RðYÞ ; where RðXÞ and RðYÞ are the rank variables of the measurements X and Y, respectively, X is either the spectral amplitudes or estimated concentrations, Y is the ground-truth concentration, cov represents covariance, and σ is the standard deviation of the rank variables. Spearman's ρ values of 1 and −1 represent perfect monotonic trends between measurements X and Y that are increasing and decreasing, respectively. It is worth noting that Spearman's ρ measures the level of monotonicity, rather than determining whether or not an evaluated function is monotonic.
In this study, we characterize jρj ≥ 0.8 as a strong monotonic trend. The mean absolute error (MAE) was used to assess the accuracy of the generated concentration maps: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 2 0 2 where C ki is the estimated concentration of pixel i for a ground-truth concentration k and N k is the total number of pixels of the concentration map k.
Finally, processing times were computed for each stage of the dual-wavelength atlas method applied on the human Hb dataset. These stages were (1) generation of coherence masks, (2) generation of acoustic spectra and spectra stacking, (3) PCA, and (4) estimation of MB concentrations, with the last stage including projection to the feature space using the principal component coefficient obtained from the atlas. These processing times were evaluated individually when training our method (i.e., constructing the atlas) and testing our method on a single frame (i.e., estimating a concentration map). The processing times for stages (3) and (4) were only calculated for training and testing, respectively. For training, robustness in the estimation of computation times was achieved by averaging the processing times from 5 different trials of 0% concentration. Similarly, robustness in the estimation of computation times for testing was achieved by averaging the processing times from every frame of the testing dataset detailed in Sec. 2.2. Computation times were measured using the MATLAB (Natick, Massachusetts) environment executed on an Intel Core i5-6600K processor with 32 GB RAM and a TITAN Xp graphical processing unit (GPU). To enhance the computation speed of stage (1) for all results reported in this manuscript, M-weighted SLSC was implemented on the GPU described above, following the architecture of preceding work. 68 In addition, this GPU-M-weighted SLSC approach and the previously-reported GPU-SLSC approach 68 were both implemented to obtain results for Eq. (1) in Sec. 2.2. Figure 2 shows the results from the photoacoustic spectra obtained with mixtures of MB and human Hb. Stacked photoacoustic spectra are shown in Fig. 2(a) with the y axis denoting mixture concentration and the x axis denoting acoustic frequency for optical wavelengths 710 (left) and 870 nm (right). Each concentration block separated by the dashed white lines in Fig. 2(a) shows 20,000 spectra samples that were randomly selected from the total number of training trials, frames, and kernels. Although a subtle frequency shift is observed in the spectra obtained with the 710-nm wavelength when increasing the MB concentration level, the spectral shift across the mixture for the 710-nm response was not strongly monotonic (ρ ¼ −0.63). Similarly, no apparent spectral shift was observed for the 870-nm response (ρ ¼ −0.27). One possible cause of the frequency shift in the spectra obtained with the 710-nm wavelength is the photoacoustic interaction with residual particles of MB in the PVCP chamber. These particles were unable to be removed as they stained the chamber walls, and they likely added frequency components to the overall spectral response of different chromophore concentrations. Figure 2(b) shows the linearity and monotonicity evaluation results as a function of the acoustic frequency for multiple image dynamic ranges. The y axes denote the ρ (left) and R 2 (right) for and linearity and monotonicity evaluations, respectively. The x axis of each plot is divided into acoustic frequencies for optical wavelengths of 710 nm (left) and 870 nm (right). In the acoustic frequency range of 0.5 to 5 MHz for the 870-nm wavelength, as the dynamic range decreases from 40 to 10 dB, the mean ± one standard deviation of improvements in Spearman's ρ and R 2 were −0.47 AE 0.15 and 0.45 AE 0.14, respectively. For dynamic ranges >40 dB, the results primarily overlap, and there is less improvement. This overlap likely occurs because the log compression reaches the resolution limit of the acquired photoacoustic amplitudes. However, for the 710-nm wavelength, only the acoustic frequency range from 2.4-5 MHz shows ρ and R 2 changes when varying the log compression from 10 to 40 dB. In the acoustic frequency range of 0.5 to 2.4 MHz, results primarily overlap regardless of the chosen dynamic range. This overlap likely occurs because of the strong spectral peak observed from 1.5 to 2.5 MHz for nearly every MB concentration in Fig. 2(a) (i.e., the spectral shift pattern described in the preceding paragraph). This strong spectral peak is robust against log compression, which generates minimal changes to Spearman's ρ and variation of R 2 .

Concentration Estimations from MB and Human Hb
Focusing on the R 2 results, multiple strong linear trends were observed when results were obtained with 870-nm laser wavelength, yielding a mean AE one standard deviation R 2 of 0.74 AE 0.14 for frequencies ranging from 0.5 to 5 MHz and thresholds from 10 to 100 dB. In contrast, no strong linearity was observed for the 710-nm region, yielding an R 2 ¼ 0.32 AE 0.21 for the same range of frequencies and thresholds as the 870-nm region. The strong spectral peak observed in Fig. 2 at the 710 nm wavelength yielded a maximum R 2 of 0.55 with a dynamic range of 10 dB, whereas the maximum peak observed at 870 nm was 0.81 with a dynamic range of 50 dB. These results suggest that Hb has a stronger influence on the linear relationship of the acoustic frequency versus concentration compared with MB. Figure 2(c) shows the linear fit for the acoustic frequency of 2.3 MHz and the dynamic range of 60 dB [i.e., parameters that produced the maximum R 2 in Fig. 2(b)]. The y axis denotes the log-compressed spectral amplitudes, and the x axis denotes the ground-truth concentrations. For each violin plot, the shape, horizontal line, white circle, and gray box represent the kernel density estimate of the data points, mean, median, and interquartile range, respectively. R 2 ¼ 0.78, suggesting a strong linearity and monotonicity between the spectral amplitudes of photoacoustic signals and ground-truth concentrations. Spearman's ρ was 0.89, which indicates that there is a high level of monotonicity. Figure 3 shows the optimization of the dual-wavelength atlas method using the human Hb dataset by maximizing the R 2 fit of a 1:1 slope of estimated versus ground-true concentrations. The variation of R 2 values as a function of the three most influential parameters of the dualwavelength atlas method is shown in Fig. 3(a). The z axis denotes the axial kernel size, and the orthogonal axes denotes the number of principal components used and the threshold for spectra log compression, respectively. Although increasing the axial kernel size improved the accuracy of the spectral amplitude generated with fast Fourier transforms, no clear improvement was observed when measuring the overall R 2 values. Similarly, parameter sets using two principal components yielded R 2 ≥ 0.7, suggesting a good estimation performance when using more than one principal component. However, no improvement was observed for PC ≥ 3. Finally, increasing the threshold from 20 to 60 dB improved the R 2 values when other parameters remained fixed, and no additional improvement was observed for dynamic range thresholds ≥60 dB. Figure 3(b) shows an example of concentration estimation versus ground-truth concentration when using the optimal set of parameters that yielded the greatest R 2 value, denoted by the dashed black box in Fig. 3(a). For this optimized estimation result, the R 2 value and Spearman's ρ were 0.80 and 0.89, respectively, indicating a high degree of linearity and monotonicity, respectively. Figure 4 shows the estimated concentrations when using porcine Hb. The sensitivity of estimated R 2 values while varying two parameters of the dual-wavelength atlas method is shown in Fig. 4(a) with the y axis denoting the threshold for spectra log compression and the x axis denoting the number of principal components used. These R 2 values were computed by fixing the axial kernel size to 1.8 mm [i.e., optimal parameter obtained in Fig. 3(a)]. Overall, the effect of varying the number of principal components and the dB threshold implemented prior to estimating the R 2 values showed similar trends to those observed in Fig. 3(a) (i.e., decreasing estimation performance with increasing the number of principal components). However, the optimal set of parameters included a threshold of 40 dB instead of the 60 dB threshold obtained in the experiment using human Hb [i.e., Fig. 3(a)], with R 2 ¼ 0.86 and ρ ¼ 0.93. Figure 4(b) shows detailed estimation results when using the optimal set of parameters from Fig. 4(a). The overall performance of the dual-wavelength atlas method when using the porcine Hb was measured by computing the MAE between the vector constructed with the mean of each violin plot and the ground-truth labels, yielding a value of 6.53%. For comparison, the MAE obtained from the optimal parameter set using the human Hb [i.e., Fig. 3(b)] was 10.49%. Thus, the dual-wavelength atlas method generated more accurate estimates of MB concentration when using porcine Hb instead of human Hb, which can be attributed to the uniformity in chemical composition and oxygenation levels of the porcine Hb samples. Figure 5 shows example DAS photoacoustic images, coherence masks, and an estimated concentration map for a ground-truth concentration of 60% of a single trial. These DAS images were normalized to the maximum amplitude value obtained from both 710-and 870-nm laser wavelength responses and displayed with 35-dB dynamic range to offer a direct comparison across wavelengths. The contours within each image represent the −3 dB boundary of the coherence mask with an area of 5.32 mm 2 . It is worth noting that these masks include the low-amplitude regions of the photoacoustic images that appear if minimal signal is present, which highlights the benefit of the coherence mask that is used to determine the presence of coherent signals, regardless of amplitude. The MAE between the concentration map and the ground truth is 9.68%. Table 1 reports the average processing times for each stage of the dual-wavelength atlas method using the porcine Hb dataset. Without considering acquisition times and transferring times of raw photoacoustic data, the construction of the dual-wavelength atlas was completed in <6 min, and the mean ± one standard deviation processing time of a single concentration map was 9.2 AE 0.5 s. This <10 min processing time is relatively fast in comparison with the total procedure times of cardiovascular interventions (e.g., the average procedure time of perfusion coronary interventions is 76 AE 31 min 72 ), and it can be performed prior to the initiation of the procedure. Figure 6 shows examples of concentration maps obtained from one trial per ground-truth concentration results. The 30% concentration map produced the greatest MAE of 14.80%, whereas the 100% concentration maps produced the lowest MAE of 0.72%. These deviations agree with the overall errors shown in Fig. 4(b), where 30% and 100% ground-truth concentrations reported MAEs of 18.29% and 0.75%, respectively. Overall, the dual-wavelength atlas  Our initial dual-wavelength atlas method used LW-SLSC beamforming to extract meaningful photoacoustic data. 46 The method herein uses M-Weighted SLSC beamforming, rather than LW-SLSC beamforming to reduce the extensive processing times required to generate coherence masks for a total of 2200 frames (i.e., 11 concentrations ×2 wavelengths ×10 frames ×5 trials ×2 datasets) of photoacoustic data. However, considering that only a single ultrasound image is needed to show the structural detail surrounding the concentration maps in Figs. 5 and 6, LW-SLSC was employed to beamform the background ultrasound image and improve the contrast, edges, and interpretation of the structural details surrounding the photoacoustic signals.

Effect of Mask Size on Concentration Estimation Performance
To illustrate the impact of mask size on estimation performance, coherence masks were generated from M-weighted SLSC images with coherence thresholds ranging from 0.3 to 0.9 in increments of 0.02. These masks were used to segment photoacoustic signals from the human blood dataset and generated concentration maps of different sizes. Then, for each coherence threshold, the MAE was calculated for the all 10 frames per 11 concentration levels. Figure 7 shows the performance of the dual-wavelength atlas method when using segmented masks generated with varying coherence thresholds. Overall, the MAE decreases when using higher coherence thresholds to generate the segmentation masks. Therefore, we selected a coherence threshold of 0.7 to generate all coherence masks in the experiment results. Figure 8 shows examples of concentration maps with a mixture of 30% MB (top) and 90% MB (bottom) when using coherence thresholds of 0.32, 0.48, and 0.70 (from left to right, respectively). The measured MAE are reported in the lower right corner of each image. For the 30% MB mixture, when increasing the coherence threshold from 0.32 to 0.70, a MAE decrease from 12.06% to 6.23% (i.e., 5.83% decrease) is observed. The reduced error is attributed to the absence of inaccurate estimates near the bottom of the 0.70 concentration map. In contrast, for the 90% MB mixture, when increasing the coherence threshold from 0.32 to 0.70, a MAE decrease of from 5.81% to 5.06% (i.e., 0.75% decrease) is observed, which can be attributed to a more uniform distribution of the mixture across the chamber. Overall, decreasing the mask size of the region of interest with more appropriate thresholding, increases the performance of the dual-wavelength atlas method, as observed in Figs. 7 and 8.

Discussion
We developed a novel photoacoustic-based, dual-wavelength atlas method that accurately estimates and generates concentration maps of a mixture of endogenous and exogenous chromophores (i.e., MB and Hb, respectively). The method builds on our previous dual-wavelength approach 46 and measures the photoacoustic spectra obtained from two laser wavelength emissions as a linear combination of the chromophore spectra stored in an atlas. Linearity and  monotonicity were confirmed with analyses of acoustic spectra and estimated concentrations as the ground-truth concentration increased, as shown in Figs. 3 and 4, respectively.
Method optimization [see Fig. 3(a)] resulted in requiring only the first principal component for feature reduction, which agrees with previous optimization results of the dual-wavelength atlas method. 46,66 However, the number of principal components as well as the axial kernel size required for accurate estimation may increase with noise and diverse photoacoustic responses obtained from different factors, such as fiber tip geometries, [73][74][75] vessel size, [76][77][78] absorber size, 79-81 ex vivo, 82 and in vivo 83,84 data. Future studies will investigate the effect of these factors in the optimization process of our method.
There are three main factors that contribute to the variability of the estimated concentration maps shown in Fig. 6. First, the energy fluctuations provided by the laser system translated into unequal fluence among frames (e.g., the mean ± one standard deviation of the laser energy at the fiber tip was 3.13 AE 0.4 mJ for the human Hb experiment). Second, decreasing the coherence threshold chosen to segment meaningful photoacoustic signals when generating M-weighted SLSC masks decreased the estimation performance of the dual-wavelength atlas method, as shown in Figs. 7 and 8. Third, although precautions were taken to minimize differences in oxygen saturation among Hb samples, such variations were not negligible. Energy fluctuations among consecutive frames and differences in oxygen-saturation levels can be reduced when using fast-tuning lasers and fresh Hb from a single patient, respectively, in more realistic clinical scenarios.
When extending our approach to other ultrasound systems, implementation with a different transducer bandwidth would require an additional IQ calibration step to maximize chromophore concentration estimation performance. Typically, the IQ compression of radiofrequency signals is conducted using a modulation frequency equal to the central frequency of the ultrasound transducer. However, as observed in the stacked spectra shown in Fig. 2, most of the frequency content resides within 1 to 4 MHz, whereas the corresponding transducer center frequency is 5.5 MHz. Therefore, the modulation frequency and bandwidth parameter of the IQ modulation step should be adjusted 66 to segment meaningful spectra from the total bandwidth of the transducer and thus enhance the estimation performance of the proposed algorithm.
In vivo deployment of the dual-wavelength atlas method for mixture estimation requires two considerations. First, given that constructing a photoacoustic atlas is relatively fast (i.e., <10 min), extracting presurgical Hb samples from a single patient would be beneficial for minimizing the estimation variability due to Hb samples from different patients at different draw times in the training set. Second, smaller vasculature are expected to yield different spectral responses because the frequency components are dependent on the volume of the chamber that is irradiated. 78,85 For reference, the diameter of the PVCP chamber simulated the average size of the inferior vena cava (i.e., 0.46 to 2.26 cm in diameter 86 ). Therefore, vessel-specific parameter optimization is a possible future direction that can be explored with various chamber diameters.

Conclusion
The work contained herein is the first to present an acoustic-based photoacoustic estimator that relies on training sets to estimate concentration levels from mixtures of photoacoustic-sensitive materials. The proposed method consisted of measuring the acoustic spectra obtained from two laser wavelength emissions as a linear combination of the chromophore spectra stored in an atlas. This linear combination assumption was confirmed with phantom experiments. In clinical practice, we envision dual excitation wavelengths illuminating the region of interest with a fast-tuning laser source, providing real-time labeling of photoacoustic-sensitive regions with a parallelized version of the algorithm. The results from the presented experiments are promising for real-time monitoring of the concentration of contrast agents in the operating room.