3.1.

Without Fluence Correction

The effect of fluence can be understood by noting that a change in amplitude at a point in a photoacoustic image, $\delta p_0$, due to a change in absorption coefficient there, $\delta {\mu }_a$ (in this case due to a change in wavelength) is given by

At first glance, Eq. (9) might suggest that a good way to satisfy requirement 1 would be to choose wavelengths for which $|\delta {\mu }_a/{\mu }_a|$ is large. However, as will be demonstrated in the example in Sec. 4.2, this would be wrong because the two terms in Eq. (9) are not independent and a large change in absorption coefficient will result in a large change in fluence. It turns out that the best way to satisfy Eq. (9) is to choose wavelengths at which the absorption coefficient is similar, i.e., $\delta {\mu }_a$ is small. While such a choice minimizes spectral coloring, it introduces another problem. When there is only a small change in the absorption coefficient and a consequent small change in fluence, there will be only a small change in the photoacoustic amplitude. It is therefore important to consider the effect of noise, specifically the noise in the difference image $\delta p(\mathbf{x})$. This leads to requirement 2.

Requirement 2 can be subdivided into two related requirements.

Requirement 2b puts a bound on the condition number of the molar absorption coefficient matrix given the noise level and the accuracy required in the ${\mathrm{sO}}_2$ estimates. In an extreme case, the molar absorption coefficient matrix will be singular (have no inverse). For example, when the wavelengths have been chosen at two isosbestic points for hemoglobin (e.g., 586 nm and 808 nm; see Fig. 4), $\delta p_0$ will be the same for all oxygenation levels, so there is no way to determine, given a measured $\delta p_0$, what the corresponding value of ${\mathrm{sO}}_2$ should be. There is no sensitivity to ${\mathrm{sO}}_2$ for this choice of wavelengths. The requirement that the matrix is not singular is another way of saying that $\delta p_0$ must map uniquely onto ${\mathrm{sO}}_2$. In a less extreme case, but one in which the condition number is large, the inversion will be very sensitive to the noise in the images, so even a low level of noise could lead to uncertainties in the ${\mathrm{sO}}_2$ estimates that are greater than its full range (0% to 100%), rendering them meaningless. An example of such a case would be where the molar absorption coefficient of Hb is almost proportional to the molar absorption coefficient of ${\mathrm{HbO}}_2$ at the chosen wavelengths, i.e., ${\alpha }_{\mathrm{Hb}}({\lambda }_n)\approx K{\alpha }_{{\mathrm{Hb}0}_2}({\lambda }_n)$, for the chosen wavelengths ${\lambda }_n$ and some constant $K$. This could arise, for example, if the two wavelengths were chosen to lie between 875 and 925 nm, where the spectral characteristics of Hb and ${\mathrm{HbO}}_2$ are similar (see Fig. 4). In this wavelength range, the measured $\delta p_0$ is only weakly dependent on blood oxygenation making accurate ${\mathrm{sO}}_2$ estimates very challenging in the presence of noise. In Ref. 42, maximizing the condition number was used as the sole criterion for choosing the optimal wavelengths. The risk of considering only condition number (which is to say the sensitivity to noise) when choosing the wavelengths, and failing also to take into account the effect of spectral coloring, is that the recommended wavelengths could well be those at which spectral coloring has the greatest deleterious effect on the accuracy.

In summary, requirement 1 must be satisfied to avoid the estimate of ${\mathrm{sO}}_2$ suffering from spectral coloring. Requirement 2 must be satisfied to avoid the estimate of ${\mathrm{sO}}_2$ being dominated by noise. Both requirements 1 and 2 apply to every image voxel independently and must be satisfied if ${\mathrm{sO}}_2$ estimates from photoacoustic images are to be considered accurate, when no attempt has been made to correct for the unknown fluence spectrum and a linear inversion [Eq. (6)] used.

3.2.

With Fluence Correction

There may be situations in which an approximate correction for the spectral coloring can be incorporated into the routine to estimate ${\mathrm{sO}}_2$. For example, an estimate of the fluence spectrum from a measurement or a model can be divided out of the measured photoacoustic amplitudes to give a measured photoacoustic spectrum closer to the absorption coefficient spectrum. To obtain accurate estimates of ${\mathrm{sO}}_2$ when using a fluence correction, similar conditions to those above apply.

3.3.

Multiwavelength Case

The two-wavelength case was studied above so that the essential requirements for accurate estimation of ${\mathrm{sO}}_2$ could be clearly identified. When multiple wavelengths are used, equivalent requirements apply: it is still necessary that the fluence variation be minimized between the wavelengths, it is still necessary that the molar absorption coefficient matrix be well-conditioned, and it is still necessary that the noise level in the measurements are sufficiently low.

4.1.

Digital Phantom

The cubic $5\times 5\times 5\mathrm{mm}$ domain had cubic 25 μm voxels and index-matched, optically absorbing boundaries. A blood-filled tube, $250\mu \mathrm{m}$ in diameter, was embedded in the domain at a depth of 1 mm parallel to the illumination surface. The domain was illuminated with a collimated 2.5 mm diameter column of light (top-hat profile). The MC model of light transport used ${10}^{12}$ photons, which was determined to be sufficient to ensure ${\mathrm{sO}}_2$ estimates were not impacted by MC noise (see Sec. 6). The arrangement is shown in Fig. 1.

Fig. 1

Digital vessel phantom: $5\times 5\times 5\mathrm{mm}$ computational domain with a blood vessel 1-mm beneath the surface and collimated wide-field illumination beam of 2.5 mm diameter. The coordinate system used throughout the paper is also shown.

The optical properties were chosen to be representative of soft biological tissue.⁴³ The blood within the tube had a total hemoglobin concentration $c_{\mathrm{HbT}}=150\mathrm{g}{\mathrm{l}}^{-1}$ and oxygen saturation of 90% ($c_{\mathrm{HbO}2}=135\mathrm{g}{\mathrm{l}}^{-1}$ and $c_{\mathrm{Hb}}=15\mathrm{g}{\mathrm{l}}^{-1}$). Two cases were considered for the background optical properties: (a) no spectral dependence and (b) spectral dependence typical of a capillary bed, modeled as 55% water and 45% blood with a total hemoglobin concentration, $c_{\mathrm{HbT}}^{\mathrm{bg}}=5.63\mathrm{g}{\mathrm{l}}^{-1}$ and 60.7% oxygenation ($c_{\mathrm{HbO}2}=3.42\mathrm{g}{\mathrm{l}}^{-1}$ and $c_{\mathrm{Hb}}=2.11\mathrm{g}{\mathrm{l}}^{-1}$). The scattering coefficient throughout the domain followed the relationship $21.3{\mathrm{mm}}^{-1}{\left(\frac{\lambda }{500\mathrm{nm}}\right)}^{-1.2}$, and the Henyey–Greenstein phase function was used, with an anisotropy factor of $g=0.9$ for the background and $g=0.9945$ for blood. The initial pressure distribution was calculated using Eq. (1) with $\kappa$ set to 1 with no loss of generality.

Figure 2 shows examples of the initial acoustic pressure and fluence distributions at wavelengths of 580 and 600 nm with a spectrally varying background [case (b)]. Both are maximum intensity projections in the $y$ direction and are normalized by the maximum value of acoustic pressure or fluence, respectively. Spectral coloring is due to the wavelength dependence of the fluence distribution, which can be clearly seen in this example, as the fluence distributions at 580 and 600 nm are notably different, especially in the region of the tube. This difference is due to the absorption coefficient of blood being an order of magnitude greater at 580 nm than at 600 nm, and it would not be surprising if it were to impact on the ability to accurately estimate oxygen saturation using these wavelengths.

Fig. 2

Spectral coloring. Top row: MIP of initial acoustic pressure in tube phantom with spectrally varying background illuminated with top hat profile at 580 and 600 nm, normalized by maximum acoustic pressure value. Bottom row: MIP of fluence in tube phantom with spectrally varying background illuminated with top hat profile at 580 and 600 nm, normalized by maximum fluence value.

4.2.

Estimates using 26 Wavelengths over the 500- to 1000-nm Range

As a starting point for this investigation, images were obtained at 26 wavelengths evenly spaced at 20-nm intervals between 500 and 1000 nm. These were used simultaneously to estimate ${\mathrm{sO}}_2$ using Eq. (5) with $\mathrm{\Phi }\propto \mathcal{I}$. Both the cases of (a) no spectral dependence in the background and (b) spectral dependence typical of tissue were studied.

Figure 3 shows a slice in the $x-z$ plane, halfway along the $y$ axis, of the error in the resulting ${\mathrm{sO}}_2$ estimates in the tube for both the spectrally flat (upper plot) and spectrally varying (lower plot) backgrounds. (Throughout the paper, error will be calculated as $|{\mathrm{sO}}_2^{\text{true}}-{\mathrm{sO}}_2|$.) The estimates in the tube differ voxel-to-voxel, and the histograms show how the distribution of the oxygenation estimates varies with depth. The voxels included in the histograms are those within the tube in the $x-y$ slices at four depths (a, b, c, d), where each $x-y$ slice is orthogonal to the $x-z$ slice shown. The ${\mathrm{sO}}_2^{\text{mean}}$ values are the means of the associated histograms. The global mean ${\mathrm{sO}}_2^{\mathrm{g}.\text{mean}}$ is the mean of the ${\mathrm{sO}}_2$ estimates in all voxels in the tube.

Fig. 3

${\mathrm{sO}}_2$ estimates using 26 wavelengths between 500 and 1000 nm. Top row: $x-z$ slice of error in oxygenation midway along tube’s length when the background is spectrally flat and the error originates from spectral coloring of the fluence in the blood-filled tube. Dashed lines used to indicate four $x-y$ slices (a–d) from which histograms of ${\mathrm{sO}}_2$ are formed (true value of 90% indicated by black line). Bottom row: $x-z$ slice of error in oxygenation midway along tube’s length when the background optical properties vary with wavelength. Both the cases with and without a depth-dependent fluence correction are shown. Dashed lines used to indicate four $x-y$ slices (a–d) from which histograms of ${\mathrm{sO}}_2$ are formed (true value of 90% indicated by vertical black line). (a–d) ${\mathrm{sO}}_2^{\text{mean}}$ represent the mean of the ${\mathrm{sO}}_2$ estimates of the voxels contained within the tube in slices a, b, c, and d, respectively, for the corrected data. The global mean ${\mathrm{sO}}_2^{\mathrm{g}.\text{mean}}$ of the corrected data, which represents ${\mathrm{sO}}_2$ estimates averaged over all voxels enclosed within the entire tube, is shown to the right of the histograms.

The histograms of the error in the ${\mathrm{sO}}_2$ estimates show that, even in the case of a spectrally flat background, the mean ${\mathrm{sO}}_2$ is at best 75% (compared to the true value of 90%) and the accuracy decreases rapidly with depth through the tube. As there is no spectral coloring from the background absorbers in this case, it is due wholly to the presence of the blood within the tube. These results provide a compelling illustration of the adverse influence of spectral coloring. They suggest that even under these highly benign conditions of a spectrally flat background, reasonable accuracy, without any fluence correction, is only possible if the spatial resolution of the imaging instrument is sufficient to resolve those few voxels that lie in the center of the $x-y$ plane within the uppermost region (slice “a”). Any spatial averaging, particularly in the $z$ direction, due to the resolution limitations of a practical imaging instrument, will rapidly produce very significant errors. For the spectrally varying background (Fig. 3 lower plot), the situation is even worse to the extent that accuracy is significantly compromised even for the “best” voxels in the center of the most superficial slice “a.” The error is larger because of spectral coloring by the background region as well as the blood within the tube itself.

The results in Fig. 3 also show that there is a tendency to underestimate ${\mathrm{sO}}_2$ relative to the true value of 90%. This is the result of the inversion being weighted by wavelengths, where the change in fluence between illumination wavelengths is large. It can be seen in Fig. 4, which shows the spectra of the absorption coefficient in the tube and the fluence averaged over the tube volume that the change in fluence is most significant between 500 and 600 nm. (In a separate set of simulations, not shown here, linear inversions only using wavelengths in the 500- to 600-nm range always resulted in the underestimation of ${\mathrm{sO}}_2$ relative to the true value of 90%, whereas when only wavelengths greater than 600 nm were used, ${\mathrm{sO}}_2$ was overestimated.)

Fig. 4

Plot of absorption spectrum in the tube (dashed blue) for 90% oxygenated blood and average fluence spectrum (green) in the tube. Also shown are the spectra for fully oxygenated blood (red) and deoxygenated blood (black). Fluence simulated using MC with ${10}^{12}$ photons, collimated, illuminated 2.5 mm in diameter at wavelengths in 500- to 1000-nm range, and averaged over tube volume.

Estimating ${\mathrm{sO}}_2$ over this wavelength range without accounting for the spectral coloring of the fluence is clearly not likely to give accurate values. To make a first-order correction for the fluence while retaining the simplicity of the linear inversion, the depth-dependent fluence correction, Eq. (11), was applied to the dataset obtained with a spectrally varying background. The resulting estimates are shown in the lower plots in Fig. 3. Comparing the corrected and uncorrected histograms, it is clear that the ${\mathrm{sO}}_2$ estimates for this case are improved by 20% to 30% and are now comparable to the case with a spectrally flat background. Although the fluence correction provides a clear benefit, the accuracy remains, at best, only comparable to that achieved in the spectrally independent case with acceptable accuracy ($\sim 6\%$) only achieved for those voxels in the center of the $x-y$ plane in the most superficial slice “a.”

In summary, the results in this section show that, when using 26 wavelengths in the range 500 to 1000 nm, spectral coloring has a highly deleterious impact on the accuracy of the ${\mathrm{sO}}_2$ estimates, to the extent that acceptable accuracy is unlikely to be achieved in practice using this wavelength range. However, Fig. 4 indicates that the error is likely to be most significant between 500 and 600 nm, where the wavelength dependence of the fluence is greatest. This suggests that, in contrast to conventional optical transmission spectroscopy, there may be an advantage to selecting wavelengths where the absorption of blood varies slowly with wavelength (which here is also where the absorption coefficient is low) because this will minimize spectral coloring. For this reason, in the next section, consideration is given as to whether more accurate estimates could be obtained by selecting wavelengths beyond 600 nm, where blood absorption and its spectral dependence is relatively low.⁴⁴

4.3.

Optimal Wavelength Pairs in the 500- to 1000-nm Range

In general, when estimating some quantity from experimental measurements, it is often reasonably assumed that the more data available the better the estimate will be.⁴⁵ The implicit intuition is that the effect on the estimate of the errors in the individual instances of the data will be symmetrically distributed about the true value and so “cancel out,” resulting in an estimate whose accuracy increases with the amount of data used. Here, while using more wavelengths may improve the SNR and may improve the conditioning of the multiwavelength inversion, the effect of spectral coloring will remain and could counteract both these advantages. Unfortunately, the coloring that the fluence gives the photoacoustic spectra is not guaranteed to cancel out, and so using images at more wavelengths does not guarantee a better estimate of ${\mathrm{sO}}_2$ when spectral coloring is not corrected for. Following this reasoning, using fewer wavelengths can potentially give more accurate estimates, but only if the spectral coloring at the wavelengths chosen is low. This will be examined in this section, but rather than investigating all possible combinations of wavelengths, pairs of wavelengths will be investigated. This has two advantages: first, the practical advantage that it is a manageable number of sets of wavelengths, and second, the advantage that the simplicity of using two wavelengths allows insights to be gleaned, which might be harder to see were larger groups of wavelengths used.

In order to demonstrate that the wavelengths that give the best estimates of ${\mathrm{sO}}_2$ are those that satisfy requirement (1), ${\mathrm{sO}}_2$ was estimated using Eq. (5) for all possible pairs of the 26 wavelengths in the 500- to 1000-nm range, i.e., ${}^{26}{C}_2=325$ wavelength pairs. Three cases were examined: (a) spectrally flat background and no fluence correction ($\mathrm{\Phi }\propto \mathcal{I}$), (b) spectrally varying background and no fluence correction, and (c) spectrally varying background and the depth-dependent correction [Eq. (11)]. Figure 5 shows histograms of the ${\mathrm{sO}}_2$ estimates for these three cases, for just the pairs of wavelengths that gave a mean ${\mathrm{sO}}_2$ closest to the true value. The voxels that satisfy the condition in Eq. (9) that the relative fluence change be much smaller than the relative absorption change (here interpreted as 2 orders of magnitude smaller) are highlighted in red.

Fig. 5

Histograms of ${\mathrm{sO}}_2$ estimates in each voxel within the tube for (a) spectrally flat with mean ${\mathrm{sO}}_2$ of 90.3%, (b) spectrally varying background with mean ${\mathrm{sO}}_2$ of 91.6%, and (c) spectrally varying with depth-dependent fluence correction with mean ${\mathrm{sO}}_2$ of 88.8% obtained using linear inversion at 620/920 nm, 700/820 nm, and 700/820 nm, respectively; ${\mathrm{sO}}_2$ estimates for which the ratio $(\delta \mathrm{\Phi }/\mathrm{\Phi })/(\delta {\mu }_a/{\mu }_a)<0.01$ shown in red and $(\delta \mathrm{\Phi }/\mathrm{\Phi })/(\delta {\mu }_a/{\mu }_a)\ge 0.01$ in black. Notes: vertical line indicates true oxygenation of 90%; there are no ${\mathrm{sO}}_2$ estimates for which $(\delta \mathrm{\Phi }/\mathrm{\Phi })/(\delta {\mu }_a/{\mu }_a)<0.01$ for the spectrally varying background (uncorrected) due to severity of spectral coloring of fluence; note the scale change on vertical axis between flat and various backgrounds. These wavelength pairs yielded the most accurate ${\mathrm{sO}}_2$ estimates of all pairs examined. ${\mathrm{sO}}_2^{\mathrm{g}.\text{mean}}$ represents the mean of the ${\mathrm{sO}}_2$ estimates all voxels enclosed within the entire tube and is termed the global mean ${\mathrm{sO}}_2$.

Comparing Figs. 3 and 5, it is evident that these wavelength pairs provide significantly improved accuracy compared to using all 26 wavelengths across the 500- to 1000-nm wavelength range. The global mean ${\mathrm{sO}}_2$ values in all cases are within 2% of the true value of 90% in Fig. 5, whereas the global mean ${\mathrm{sO}}_2$ values achieved when using the 500- to 1000-nm spectral range are subject to an error of almost 30% (Fig. 3). When the background is spectrally flat, Fig. 5(a) shows that ${\mathrm{sO}}_2$ is accurately estimated (to within 0.5%) and that estimates made under the condition that the relative change in the absorption coefficient far exceeds the accompanying relative change in fluence are more accurate. These estimates correspond to regions at the surface of the tube and at the center of the $x-y$ slice. When the background has some spectral dependence, ${\mathrm{sO}}_2$ estimates are less accurate with errors of 1% to 2%, and much higher variance, as shown Fig. 5(b). It can be observed that there are no voxels that satisfy Eq. (9) in Fig. 5(b) due to spectral coloring being significant at 700/820 nm. Figure 5(c) shows that the use of a depth-dependent fluence correction shifts estimates to a global mean value of 88.8% in the tube, and again the values closest to the true value are those that satisfy Eq. (9).

Given that in most practical settings the ratio $(\delta \mathrm{\Phi }/\mathrm{\Phi })/(\delta {\mu }_a/{\mu }_a)$ is unknown, and therefore cannot be used directly as an indicator of the accuracy of the ${\mathrm{sO}}_2$ estimate, it is desirable to obtain a result that is of more general applicability. One could, for example, determine the wavelength pairs $({\lambda }_1,{\lambda }_2)$ that yield the most accurate ${\mathrm{sO}}_2$ estimates on average over the tube, thereby providing some indication of which wavelengths are best-suited to recovering ${\mathrm{sO}}_2$. Figure 6 illustrates this by showing the error in oxygenation, averaged over the tube, obtained at a range of wavelengths ${\lambda }_1$ and ${\lambda }_2$ between 500 and 1000 nm when the background was spectrally varying. Four different tube ${\mathrm{sO}}_2$ values were considered. The bottom-left triangles show estimates obtained without the depth-dependent fluence correction, whereas the top-right triangles contain values obtained when the data were corrected using Eq. (11).

Fig. 6

Matrix of average error in oxygenation for 26 wavelength pair combinations in 500- to 1000-nm range for different values of oxygenation in the tube, ${\mathrm{sO}}_2^{\text{tube}}$. In each figure, the bottom-left triangle represents ${\mathrm{sO}}_2$ estimates made with a spectrally varying background and the top-right triangle represent ${\mathrm{sO}}_2$ estimates made with a spectrally varying background using a depth-dependent fluence correction. The white pixels represent estimates with an error $>100\%$. In the regions marked “a,” this error is due to spectral coloring; in the regions marked “b,” it is due to poor conditioning. Top left: ${\mathrm{sO}}_2^{\text{tube}}=70\%$. Top right: ${\mathrm{sO}}_2^{\text{tube}}=80\%$. Bottom left: ${\mathrm{sO}}_2^{\text{tube}}=90\%$. Bottom right: ${\mathrm{sO}}_2^{\text{tube}}=100\%$.

Several observations emerge from the results in Fig. 6. First, it is noticeable that there are two distinct spectral regions “a” and “b,” indicated by white-colored pixels, where the error is particularly high ($>100\%$). The error is large for region “a” because it encompasses the spectral range of 500 to 600 nm, where blood absorption is strongly wavelength-dependent, resulting in significant spectral coloring. Region “b” represents the spectral range of 820 to 900 nm and exhibits poor accuracy for a different reason, namely, poor conditioning of the molar absorption coefficient matrix as the spectra of Hb and ${\mathrm{HbO}}_2$ vary almost in proportion, i.e., the ratio ${\alpha }_{\mathrm{Hb}}/{\alpha }_{{\mathrm{HbO}}_2}$ is approximately constant in this wavelength range. Second, if the wavelengths that give the most accurate estimate of ${\mathrm{sO}}_2$ are selected, a depth-dependent fluence correction may not provide a noticeable improvement over no correction in the accuracy of the ${\mathrm{sO}}_2$ estimates. This is consistent with Figs. 5(b) and 5(c) and is likely to be due to the fact that the spectral coloring is less significant when wavelengths are well-chosen. The exception to this is where one wavelength lies in the 500- to 580-nm range. In this case, the application of the depth-dependent fluence correction provides a benefit as spectral coloring is more severe in this wavelength range. This may be important in applications that necessitate the use of shorter wavelengths, for example, in situations where high SNR is required or the excitation wavelength range of the laser source is limited.

Figure 6 also shows that in each case inversions, whether with or without the fluence correction, yield the highest accuracy at wavelengths between 620 and 1000 nm (excluding region “b”), with few wavelength pairs outside this range producing accurate ${\mathrm{sO}}_2$ values. Moreover, the number of wavelength pairs within this wavelength range yielding an error between 0% and 100% (the range of the color bar in Fig. 6) does not vary significantly for the four different tube ${\mathrm{sO}}_2$ values. However, the number of wavelength pairs yielding reasonable accuracy (an error $<10\%$) is strongly dependent on the value of ${\mathrm{sO}}_2^{\text{tube}}$ and decreases with decreasing ${\mathrm{sO}}_2^{\text{tube}}$. This can be seen most clearly by comparing the maps for ${\mathrm{sO}}_2^{\text{tube}}=70\%$ and ${\mathrm{sO}}_2^{\text{tube}}=100\%$. These show that, while almost any pair of wavelengths between 620 and 1000 nm is likely to yield reasonable accuracy for ${\mathrm{sO}}_2^{\text{tube}}=100\%$, the same does not apply at ${\mathrm{sO}}_2^{\text{tube}}=70\%$, where there are significantly fewer wavelengths yielding accurate estimates. This is because the wavelength dependence of blood absorption increases with decreasing ${\mathrm{sO}}_2$ for wavelengths below the isosbestic point at 800 nm.

At low ${\mathrm{sO}}_2^{\text{tube}}$ values, spectral coloring, therefore, increases for the majority of wavelengths in the 620- to 1000-nm range, so there are fewer wavelengths that provide the low levels of spectral coloring required for accurate estimates. Moreover, at low ${\mathrm{sO}}_2$ values, those wavelength pairs that do yield accurate estimates become increasingly sparsely distributed in wavelength space. Under the conditions of low ${\mathrm{sO}}_2^{\text{tube}}$, the combination of these factors makes it very challenging to select the optimal wavelengths since they depend critically on the tube ${\mathrm{sO}}_2$, which is generally unknown a priori in practice. (Interestingly, additional simulations, not shown, showed that, while the optimal wavelengths were dependent on the tube ${\mathrm{sO}}_2$, they remained relatively independent of changes in the background ${\mathrm{sO}}_2$. This may be due to the relatively low background absorption in this simulation.)

In summary, the results in this section suggest that, unlike the 500- to 600-nm wavelength range, where spectral coloring severely compromises accuracy under almost all conditions, it may be possible to achieve reasonable accuracy in some circumstances by using longer wavelengths between 620 and 1000 nm. For example, for high ${\mathrm{sO}}_2^{\text{tube}}$ values ($>90\%$), a modest level of accuracy may be achieved by selecting any pair of wavelengths in this range. This relaxed wavelength selection criterion also has the important practical advantage that accurate estimates are not then contingent upon accurate a priori knowledge of ${\mathrm{sO}}_2^{\text{tube}}$, just that it is around 90% or higher. This could be useful, e.g., when measuring ${\mathrm{sO}}_2$ in an artery, where the physiological range of ${\mathrm{sO}}_2$ is quite narrow. However, for lower ${\mathrm{sO}}_2^{\text{tube}}$ (e.g., 70%), wavelength selection becomes more problematic as there are only relatively few sparsely distributed wavelength pairs that provide reasonable accuracy and identifying these wavelengths requires a priori knowledge of ${\mathrm{sO}}_2^{\text{tube}}$, which is generally unavailable in vivo. Even our seemingly encouraging conclusion that blind optimal wavelength selection within the 600- to 1000-nm range is feasible when ${\mathrm{sO}}_2^{\text{tube}}$ is high ($>90\%$) may be unduly optimistic since it is based on the assumption of noise-free data. In practice, any imaging system will introduce noise, which may make certain wavelengths, where the absorption coefficient is low, impractical to use as the SNR would be too low, as discussed in Sec. 4.4.

4.4.

Signal-to-Noise Ratio

The suggestion above to select wavelengths above 620 nm, where blood absorption changes less with wavelength than in the 500- to 600-nm range, means SNR is sacrificed in favor of the accuracy of the oxygenation estimates. This trade-off is made evident by considering that we require that the change in the initial pressure, $p_0$, as a result of the change in wavelength from ${\lambda }_1$ to ${\lambda }_2$, exceeds the minimum detectable change in $p_0$, the “image noise equivalent pressure”: $|\delta p_0|>\eta$, or

The inequality in Eq. (12) was applied to the data in Fig. 6 in order to illustrate the effect that taking SNR into account can have on the wavelength selection. This was performed for 325 wavelength pairs in the 500- to 1000-nm range and typical values for the parameters in Eq. (12) were used: $\eta =0.2\mathrm{kPa}$ and $\mathrm{\Gamma }=0.15$.⁴⁷ The fluence was computed using the MC-generated data used in the inversion and was scaled by an incident fluence of $20\mathrm{mJ}{\mathrm{cm}}^{-2}$. Using these parameters yielded mean image SNRs in the range 26.5 to 45.7 dB, depending on wavelength and the value of the tube ${\mathrm{sO}}_2$. These SNRs are significantly higher than those of high quality in vivo photoacoustic images of superficial blood vessels^48,49 and therefore on the optimistic side of what can be achieved in practice. The data in Fig. 6 were then replotted in Fig. 7 but masking the pairs of wavelengths for which the above inequality in Eq. (12) is not satisfied. In other words, those wavelengths for which the change in initial pressure from ${\lambda }_1$ to ${\lambda }_2$ was less than $\eta$.

Fig. 7

Matrix of average error in oxygenation for 26 wavelength pair combinations in 500- to 1000-nm range for different values of oxygenation in the tube, ${\mathrm{sO}}_2^{\text{tube}}$. In each figure, the bottom-left triangle represents ${\mathrm{sO}}_2$ estimates made with a spectrally varying background and the top-right triangle represents ${\mathrm{sO}}_2$ estimates made with a spectrally varying background using a depth-dependent fluence correction. The white pixels represent estimates with an error $>100\%$. In the regions marked “a,” this error is due to spectral coloring; in the regions marked “b,” it is due to poor conditioning. Light blue regions indicate wavelengths that do not satisfy the condition in Eq. (12). (A) ${\mathrm{sO}}_2^{\text{tube}}=70\%$, (B) ${\mathrm{sO}}_2^{\text{tube}}=80\%$, (C) ${\mathrm{sO}}_2^{\text{tube}}=90\%$, and (D) ${\mathrm{sO}}_2^{\text{tube}}=100\%$.

For all the oxygenation levels studied, the number of wavelengths that yield accurate estimates in the presence of noise is considerably reduced compared to the noise-free case (shown in Fig. 6). For the remaining wavelength pairs, the depth-dependent fluence correction can provide an improvement in accuracy, particularly for low values of ${\mathrm{sO}}_2^{\text{tube}}$. This is because for low values of ${\mathrm{sO}}_2^{\text{tube}}$, most of the wavelength pairs that satisfy Eq. (12) comprise at least one wavelength that lies within the spectral range of 500 to 620 nm, where blood absorption is high and varies rapidly with wavelength. The resulting spectral coloring can, therefore, be mitigated by the fluence correction.

In terms of optimal wavelength selection, there are several spectral regions that yield reasonable accuracy. For ${\mathrm{sO}}_2^{\text{tube}}=100\%$ and ${\mathrm{sO}}_2^{\text{tube}}=90\%$, the spectral ranges defined by regions “a” and “b” in Fig. 7 provide accurate estimates with errors predominantly below $\sim 10\%$. However, the number of wavelength pairs in these regions decreases with decreasing ${\mathrm{sO}}_2^{\text{tube}}$ to the extent that these regions disappear entirely for ${\mathrm{sO}}_2^{\text{tube}}=70\%$. This is because $\delta {\mu }_a$ becomes progressively smaller with decreasing ${\mathrm{sO}}_2$ for wavelengths that lie in these spectral ranges until it no longer satisfies Eq. (12). Perhaps, a more useful spectral range is that defined by region “c,” where one wavelength is 600 nm and the other must lie between 800 and 1000 nm. Any wavelength pair lying in region “c” will not only satisfy Eq. (12) but will also yield reasonable accuracy for all the ${\mathrm{sO}}_2^{\text{tube}}$ values represented in Fig. 7. This is because blood absorption at 600 nm is high so $\delta {\mu }_a$ arising from the wavelength change ${\lambda }_1=600\mathrm{nm}$ to ${\lambda }_2>800\mathrm{nm}$ is large enough to satisfy Eq. (12) but not so high as to introduce excessive spectral coloring. However, although accuracy appears to be reasonable for these wavelength pairs, it does decrease significantly with reducing ${\mathrm{sO}}_2^{\text{tube}}$. For example, for the spectral range defined by region “c” and for ${\mathrm{sO}}_2^{\text{tube}}=100\%$, if ${\lambda }_2=600\mathrm{nm}$, then the error varies from 2.8% to 8.1% depending on the choice of ${\lambda }_1$. By contrast for ${\mathrm{sO}}_2^{\text{tube}}=70\%$, the corresponding error varies from 7.7% to 18% for region “c.” The wavelength pair that yielded the best accuracy across all four ${\mathrm{sO}}_2$ values in Fig. 7 was found to be $600/800\mathrm{nm}$ with an error ranging from 9.8% for ${\mathrm{sO}}_2^{\text{tube}}=70\%$ to 3.6% for ${\mathrm{sO}}_2^{\text{tube}}=100\%$.

Although this appears to suggest that reasonable accuracy may be possible using the $600\text{-}/800\text{-}\mathrm{nm}$ wavelength pair, this should be treated with caution because the spectral region for one of the wavelengths is narrow, and a higher noise level or a small mismatch between the model and an experimental scenario could easily shift or eliminate it. In this study, optimistic assumptions of perfect image reconstruction and a homogeneous background have been made, so there are no artifacts, and the typical SNR achieved in vivo is often lower than used here. The latter will have the effect of reducing the number of wavelengths that satisfy Eq. (12) and for those remaining wavelengths that do satisfy it, they will be biased to shorter wavelengths, where spectral coloring is greater and accuracy reduced. Furthermore, the specific wavelength pairs that do yield accurate estimates will be increasingly dependent on the ${\mathrm{sO}}_2^{\text{tube}}$ values, which are unknown, making it hard to identify those wavelength pairs with confidence.

In summary, in the previous section it was observed that it may be possible to achieve accurate estimates for ${\mathrm{sO}}_2$ values by selecting wavelength pairs in the 600- to 1000-nm range. However, noise has the negative impact of not only significantly reducing the number of wavelength pairs that provide accurate estimates but also making the specific values of these wavelengths increasingly dependent on the ${\mathrm{sO}}_2^{\text{tube}}$. Indeed, it is the latter that arguably presents the greatest challenge since ${\mathrm{sO}}_2^{\text{tube}}$ is typically unknown a priori, determining it being the point of the measurement, so there is no certain way of identifying those few wavelength pairs that yield accurate estimates.

6.1.

Level of Uncertainty in ${\mathrm{sO}}_2$ Estimates Due to MC Noise

Although MC models are considered the “gold standard” (Ref. 51), it is important to recognize that due to the stochastic nature of simulations, the estimate of the fluence in each voxel will exhibit some variance. As the number of photons increases, the variance will reduce, and the fluence estimate will converge to the solution given by the RTE. For these simulations, it is necessary to use sufficient photons that the variance (the “MC noise”) falls below an acceptable level. To assess what level this needs to be, we examined the effect of the MC noise on two-wavelength ${\mathrm{sO}}_2$ estimates calculated from initial pressure distributions generated using the MC models. The standard deviation in oxygenation at a particular voxel, ${\sigma }_{{\mathrm{sO}}_2}$, can be estimated from the standard deviation in the initial acoustic pressure simulated for that voxel, ${\sigma }_{p_0({\lambda }_i)}$, by assuming that the initial pressure at the two wavelengths, $p_0({\lambda }_1)$ and $p_0({\lambda }_2)$, is uncorrelated⁵² and using the standard formula: