Multiple irradiation sensing of the optical effective attenuation coefficient for spectral correction in handheld OA imaging.

Spectral optoacoustic (OA) imaging enables spatially-resolved measurement of blood oxygenation levels, based on the distinct optical absorption spectra of oxygenated and de-oxygenated blood. Wavelength-dependent optical attenuation in the bulk tissue, however, distorts the acquired OA spectrum and thus makes quantitative oxygenation measurements challenging. We demonstrate a correction for this spectral distortion without requiring a priori knowledge of the tissue optical properties, using the concept of multiple irradiation sensing: recording the OA signal amplitude of an absorbing structure (e.g. blood vessel), which serves as an intrinsic fluence detector, as function of irradiation position. This permits the reconstruction of the bulk effective optical attenuation coefficient μeff,λ . If performed at various irradiation wavelengths, a correction for the wavelength-dependent fluence attenuation is achieved, revealing accurate spectral information on the absorbing structures. Phantom studies were performed to show the potential of this technique for handheld clinical combined OA and ultrasound imaging.

Introduction

In recent years, there has been a growing interest in the development of quantitative spectral optoacoustic (OA) imaging for medical diagnostics [1]. This technique exploits the wavelength-dependent optical absorption properties of specific chromophores in tissue (e.g. hemoglobin), to provide quantitative estimates of their spatially varying concentrations [2], [3], [4], [5], [6], [7]. A physiologically important example is the determination of local blood oxygen saturation, based on the distinct absorption spectra of oxy- and deoxyhemoglobin in the near-infrared range [8], [9]. This is of particular relevance for the study of oxygenation heterogeneity in tumors, the early detection and monitoring of cerebral ischemia in brain, and of other abnormalities characterized by a change in tissue oxygenation or perfusion status [10].

Quantitative spectroscopic OA imaging, however, faces important challenges, as well-described by Cox et al. [11]; and among them is the spectral distortion due to wavelength-dependent fluence attenuation. The retrieval of accurate functional information will only be possible if this wavelength-dependency of the optical attenuation in the tissue is taken into account (which we shall refer to as “spectral correction”) [12], [13]. Many groups have been working on this problem, proposing a number of different correction approaches, each of which has certain advantages and disadvantages [4], [13], [14], [15], [16], [17], [18]. An ideal method to estimate the fluence distribution would be non-invasive, not rely on any additional modality, or any a priori knowledge.

In this study, we present a technique that aims at fulfilling all the aforementioned criteria. Our approach ties in with the concept introduced by Zemp et al. [19], [20], which employs multiple irradiation positions to estimate optical properties of highly scattering media, a concept similar to optical tomography. However, as opposed to superficial optical detection, it employs optically absorbing structures inside the tissue as intrinsic “fluence detectors”. The practical feasibility of this concept was shown in the context of OA microscopy, demonstrating for a single-wavelength a depth-dependent correction for the optical attenuation [21], [22]. The work described in the present study constitutes a natural continuation of the aforesaid study, in that the multiple-irradiation technique is extended to a spectral correction, which is essential to retrieve quantitative functional information. Further, to emphasize the suitability of this technique for clinical OA imaging with a handheld probe, we make use of a linear array transducer for 2D image acquisition.

We conducted phantom experiments with a tissue-mimicking bulk material showing a spectrally varying optical attenuation in which vessel-like cylindrical absorbers characterized by a distinct absorption peak were embedded. Our results demonstrate that the multiple optical irradiation approach allows a successful correction of the measured OA spectrum of the absorbers, which was initially distorted by the spectral attenuation of the surrounding medium.

Theory

Spectral optoacoustic signal amplitude

The initial pressure, p(r, t = 0) = p0,(r), generated at position r after tissue irradiation, acts as a source of acoustic energy emitting a broadband transient wavefield (OA signal). The initial pressure can be described bywhere the Grueneisen parameter Γ(r) characterizes the thermo-elastic expansion of the tissue, μ(r) is the wavelength-dependent absorption coefficient of the chromophore of interest and Φ(r) refers to the fluence at a position r. Under ideal conditions (i.e., full information, thus no detection bandwidth limitation, sample-encompassing detection area), p0,(r) can be reconstructed from the OA signal p(r′, t) measured at detection points r′, by solving an acoustic inverse problem. The reconstructed signal, S(r), can be expressed as

Eq. (2) highlights the proportionality between the reconstructed signal, S(r), at any point r = (x, y, z) inside the investigated tissue and the absorbed energy density, α(r), which corresponds to the product of μ(r) and the fluence Φ(r). This holds true if the product of μ(r) and diameter of the absorbers is significantly smaller than 1 (see [23]). For quantitative spectroscopic OA imaging it is therefore inevitable to have knowledge on the light fluence distribution Φ(r), which itself is not self-evident to formulate, as it depends in a complex way on the tissue geometry and composition. To describe Φ(r) in a simplified situation, the following significant assumptions were made:

First, we presume that the tissue can be modeled as an optically homogeneous bulk material, which is strongly scattering and weakly absorbing. In the diffusion approximation for light propagation in a semi-infinite 3D geometry (infinite along x − y plane and with a boundary at z = 0, infinite along z > 0 – see next section for more details)1 the fluence Φ0,(r|μ) at a virtual intrinsic detection point r = (x, y, z) inside the bulk material is fully determined by the optical effective attenuation coefficient μ. For a point-like irradiation, it can be analytically expressed as:where z denotes the depth of the virtual detection point, and d the distance, from the said point to the source position.

Secondly, we presume that absorbing structures (e.g., blood vessels containing specific chromophores, ) which are optically more absorbing than the background, are sparsely distributed within the bulk medium. The assumption is made that these absorbing structures do not affect the light propagation, so that the actual fluence inside the tissue is undisturbed and can be taken equal to Φ0,(r|μ) (This zero-order assumption is akin to the first Born approximation commonly used in light scattering problems). As such, S(r) can be written as:

Optically absorbing structures as intrinsic fluence detectors

If Φ0,(r|μ) is known for all wavelengths of interest, the undistorted absorption spectrum of the absorbing structures can be retrieved by inverting Eq. (4):

As Φ0,(r|μ) is fully determined by μ, the bulk effective attenuation coefficient can be derived by “measuring” the fluence for various source-detection distances d and fitting Eq. (3) to the outcome. A conventional way to measure the fluence would be to use an external optical detector (e.g., a photodiode) [25], but this requires additional equipment and can be rather cumbersome in a clinical environment. Instead, we make use of the judicious approach proposed in [19], which suggests using the investigated absorbing structures themselves as point-like “intrinsic fluence detectors”. As the OA signal amplitude is proportional to Φ0,(r|μ), the idea is that successively acquired OA images at different λ and with varying irradiation positions effectively sample the fluence for different source-detection distances d.

Materials and methods

Experimental model

To demonstrate the feasibility of the proposed method, we have conducted phantom experiments in a well-defined and controllable setting, where the phantom was chosen to fulfill the assumptions made so far. A liquid fat emulsion was taken as a homogeneous and spectrally attenuating bulk medium because (i) it is commonly considered as a fairly realistic model for optical tissue properties [26], (ii) it is strongly scattering and weakly absorbing, and (iii) it is well-studied, easy to manipulate and perfectly reproducible. The fat emulsion was placed in a large cuboid container, so that boundary effects can be neglected. To mimic sparsely distributed absorbing structures (e.g. blood vessels), cylindrical tubings were embedded in the liquid bulk medium, aligned on top of each other, in parallel to the y-axis (see Fig. 1). The total number of absorbers was limited to three, to avoid excessive reverberations of OA transients between the different absorbers. The tubings were filled with a dye solution, whose absorption spectrum displays a significant wavelength dependency. The tubings were chosen thin to minimize the influence on Φ0,(r|μ), while still providing a good signal-to-noise ratio of the OA signal.

Fig. 1

(a) Experimental model including a spectrally-attenuating bulk medium and three embedded dye-filled tubings; the light propagation distance d is defined as the distance between the ith irradiation spot position and the intersection of the jth cylindrical absorber with the transducer imaging plane. (b) Experimental setup, showing the mounting of the linear transducer and the multimode fiber attached to the translation stage. (c) 2D cross-section of the experimental model, indicating the locations of the three absorbers at depths z laterally centered in the imaging plane. (d) 2D reconstructed OA image showing the vertical alignment of the absorbers inside the spectral attenuating bulk medium. (For interpretation of the references to color in text, the reader is referred to the web version of the article.)

Fig. 1a and c shows a sketch of the experimental model, illustrating the position and orientation of the absorbers inside the bulk medium, the irradiation optics and the linear transducer array, which is placed on top of the phantom surface for OA signal detection. The x–z plane underneath the transducer, represented by a red-dashed rectangle in Fig. 1c, defines the so-called “imaging-plane”, i.e., the plane of maximum sensitivity of the linear US transducer.

To fulfill the requirements needed for the validity of Eq. (3) in terms of irradiation geometry, an optical fiber was used to approximate a point-like irradiation. The irradiation spot was translated along the phantom surface, parallel to the linear transducer in the x-direction, leading to different distances d between the irradiation spot and the virtual fluence detection point r, i.e., the intersection of the cylindrical tubings with the imaging plane (Fig. 1a). The ith irradiation spot position is described by x, the depth of the jth dye-filled tubing of interest by z and the fixed spacing distance between the transducer imaging plane and the fiber translation axis by y. This distance was set to 16 mm, dictated by the finite width of the linear transducer array. The orientation of the linear transducer array, aligned perpendicularly to the embedded absorbers, enabled a 2D reconstruction that represents a cross-sectional view of the dye-filled tubings within the imaging plane. Fig. 1d exemplarily shows an OA image acquired with one irradiation spot position x at one irradiation wavelength λ. It displays the spatially-resolved OA signal, S(r), of the embedded absorbers at various depths z.

In Section 2.1 we made the assumption of a calibrated OA detection system providing full information. Yet in practice, a typical linear array transducer has a limited aperture and detection bandwidth, thus not strictly fulfilling the made assumption. Despite that, S(r) can be regarded as an entity, which is, in case of a moderate optical absorption, proportional to the absorbed energy density α(r), at the cross-sectional position r of the absorbers in the imaging plane and is thus proportional to the fluence, Φ0,(r|μ), “detected” at those positions. The cross-section points of the imaging plane with the vessels act in this manner as intrinsic fluence detectors.

Phantom preparation

The background medium was prepared by diluting an optically scattering stock fat emulsion (SMOFlipid 20%, Fresenius Kabi, Switzerland) with de-ionized water.2 To produce a total volume of 4100 ml, 6 wt% of the stock fat emulsion was mixed with de-ionized water, resulting in a total fat concentration of 1.2 wt%. Optical reference measurements of the fat emulsion were performed in the spectral range from 700 nm to 830 nm, using various optical characterization techniques (see supplementary material). The spectrally-absorbing dye Indocyanine green (ICG) (IR-125, Laser Grade, ACROS Organics, NJ, USA), diluted in ethanol (Ethanol, p.a., Sigma–Aldrich), was used to fill the embedded cylindrical tubings. The absorption coefficient of the ICG–ethanol solution was assessed from extinction measurements using a transmission photo-spectrometer. The measured absorption spectrum agrees with literature values [27]. The peak absorption coefficient value ( 0.35 mm−1 at 785 nm) has the same order of magnitude as the absorption of hemoglobin for the spectral range of interest (see Section 5).

Measurement setup

An overview of the experimental setup is given in Fig. 1b. The cuboid container for the bulk medium was made of polymethylmethacrylate (PMMA) and had geometrical dimensions of 23 cm (length), 17 cm (width) and 14 cm (height). Cylindrical polyurethane (PU) tubings (Harvard Apparatus GmbH, Germany), with an inner-diameter of 0.28 mm and an outer-diameter of 0.6 mm, were used as embedded absorbers filled with the dye solution. The tubings were placed at depths of 14.8 mm, 25.4 mm and 28 mm, respectively. For OA signal generation, a diode-pumped Q-switched Nd:YAG laser (Spitlight DPSS OPO, InnoLas Laser GmbH, Germany) with an integrated optical parametric oscillator (OPO) was used, delivering laser pulses with a tuneable wavelength in the range of 680–900 nm, a pulse duration of 10 ns and a repetition rate of 100 Hz. The laser light was coupled into a 1000 μm core, multimode fiber (Thorlabs Inc., NJ, USA), delivering a wavelength-dependent average output pulse energy at the fiber tip of around 4 mJ. The energy of each single pulse was measured before the fiber coupling with a pyroelectric energy sensor. The fiber terminal was fixed on a motorized translation stage (T-series, Zaber Inc., Canada), allowing an automated, stepwise translation of the irradiation spot parallel to the x-axis. The OA pressure transients, generated from the absorbing inclusions, were detected by a linear array transducer (ATL L7-4, Philips N.V., NL) containing 128 elements at a pitch of 0.298 mm with a center frequency of 5 MHz and fractional bandwidth of 80%. OA signal acquisition and reconstruction was carried out with an ultrasound research system (V 1-64, Verasonics Inc., Redmond, WA), which allowed parallel signal readout of 64 channels and continuous data transfer for real-time processing on a host PC. This system facilitates the acquisition and reconstruction of a full OA image with each single laser pulse, at a frame rate corresponding to the laser pulse repetition rate. The software of the ultrasound research system was programmed in order to alternately acquire US and OA images, which offer a direct one-to-one comparison between the two modalities.

For image acquisition, the transducer was mounted with its imaging plane perpendicular to the cylindrical PU tubings, permitting a cross-sectional image with a high signal-to-noise ratio (SNR) [28]. The fiber tip was dipped into the liquid at a distance of 16 mm from the imaging plane (y-direction). The fiber tip and the transducer aperture were aligned to be at the same depth z relative to the liquid surface. The initial irradiation location x0 corresponded to the minimal light propagation distance to all absorbers. The fiber tip was then laterally translated over a total distance of 25 mm, with a step size of 0.5 mm. For each irradiation position an OA image was acquired. To reduce any non-systematic background, such as thermal noise, 100 OA images were averaged for each irradiation position. The OA images (Fig. 1c) were reconstructed using a frequency domain algorithm [29].

Correction procedure

The implementation of the proposed method for the spectral correction can be summarized in the following steps:

Positioning of the transducer array, as described in the previous section.

Acquiring N * M OA images for N irradiation wavelengths, λ, and M different irradiation spot positions, x.

Determining the OA signal amplitudes, S(r, x), of the jth detection point for each λ and each x by fitting a 2D Gaussian profile to the reconstructed OA signal envelope and taking the resulting prefactor as a measure for the amplitudes.

Estimating μ separately for each detection point and each λ by fitting the diffusion model to the OA amplitude values (as function of d, instead of r as in Eq. (3)), , where d denotes the distance between the ith irradiation spot position and the jth detection point.

Calculating the undisturbed fluence, {Φ0,(r)|μ}, for each λ in the tissue, using the previously estimated μ coefficient values as input to Eq. (3), which is now taken as a forward model.

Finally, knowing Φ0(r|μ), the correction of the spectral optical attenuation in the tissue is possible, according to Eq. (5).

Results

Estimation of the optical effective attenuation

In Fig. 2, the OA signal amplitude, S(r), is exemplarily shown for one wavelength, λ = 760 nm, as function of source-detection distance, d. For visualization purposes, the signals yielded by the three absorbers are represented by different symbols and normalized to the maximum amplitude for each absorber. The fluence model Φ0,(r|μ) (dotted lines) was fitted to S(r), to retrieve the bulk medium's effective attenuation coefficient μ. It is apparent that the fit accurately describes the decay of the measurements for all three absorbers.

Fig. 2

OA signal amplitudes of the three absorbers as function of source-detection distances at λ = 760 nm. For better visual comparison, all amplitude values were normalized to the maximum values of the respective absorber and plotted on a logarithmic scale. The dotted lines are the result of the fit of the diffusion model to the measured amplitude values, whereas the solid lines show the results of MC simulations based on the μ values determined by the reference optical measurements.

The estimated μ values, using the different absorbers as fluence detectors, are in good agreement with the results of the independently performed optical characterization described in the supplementary material. Table 1 summarizes the results with an error estimation based on the 95% confidence interval (CI) of the least-square fit. The estimated μ values based on the first and second absorber lie within the range of the reference measurements. Yet, using the third absorber, μ shows an underestimation relative to the reference by around −0.02 mm−1: This is visible in Fig. 2, where the fitted slope for the third absorber differs from those obtained with the first and second absorbers. This is further discussed in Section 5.

Table 1

Estimated μ values of the bulk medium for one exemplarily selected wavelength (λ = 760 nm), together with the results obtained from the independently performed reference measurements using frequency-domain diffuse optical spectroscopy (ImagentTM).

jth absorber	μeff,λ=760 nm
j = 1	0.092 mm−1 ± 0.0041
j = 2	0.089 mm−1 ± 0.0072
j = 3	0.071 mm−1 ± 0.01
Reference
ImagentTM (760 nm)	0.099 mm−1 ± 0.009

To validate the applicability of the diffusion approximation to our experimental model, we performed Monte Carlo (MC) simulations by thoroughly reproducing the phantom and irradiation geometry. All simulations were conducted by generating 1.5 × 109 photon paths with our general purpose MC software “jaMCp3” [30], [31], [32], [33]. The input optical parameters (μ(r), μ(r), and g(r)) were taken from the reference measurements. The MC simulations of the OA amplitudes agree very well with the decay of the experimental data and the diffusion model, supporting the assumption that the applied model accurately describes the light propagation in the considered experimental situation, where boundaries other than the surface could be neglected.

Fig. 3 presents the estimated μ values retrieved from a single experiment over the full investigated spectral range, from 700 nm to 830 nm. These values are based on fitting the fluence model Φ0,(r|μ) to N = 27 sets of measurements, one for each λ. A single set is in turn composed of M = 50 measurements (each obtained by averaging 100 OA signal acquisitions) corresponding to different irradiation spot positions, x. The error estimation for each set is based on the 95% CI of the single least-square fit. (Note that this entire experiment was repeated five times in total, with a reproducibility error ≤6%, which lies within these 95% CI values.)

Fig. 3

μ values of the background medium as determined using the different absorbers, displayed separately for each absorber (a) j = 1, (b) j = 2, (c) j = 3, over the spectral range of interest (700–830 nm). The results of the independent optical characterization, are presented as reference (gray dashed line, black circles, orange squares). (For interpretation of the references to color in this legend, the reader is referred to the web version of the article.)

For comparison, the results obtained from each absorber are separately plotted along with the reference μ values. The experimental results based on the first and second absorber (Fig. 3a and b) show a very good quantitative agreement with the reference values, except for the slight overestimation between 700 and 720 nm. In contrast, the results obtained with the third absorber show a good agreement with the reference values between 700 and 730 nm, but a strong underestimation of μ for λ > 730 nm.

Spectral correction

The results in Fig. 3 underline the spectrally varying attenuation of the fluence within the bulk medium. This leads to a distorted OA spectrum with respect to the actual absorption spectrum of the investigated ICG solution, as shown in Fig. 4a, c and e. The observed distortion is getting more pronounced for deeper absorber positions. It is particularly accentuated at the spectral range from 700 nm to 790 nm, where μ of the bulk medium strongly changes from around 0.04 mm−1 to 0.09 mm−1 (see Fig. 3).

Fig. 4

OA signal spectra of the 3 embedded absorbers as function of wavelength, {S(r)} and the reference absorption spectrum, (red solid line). For each absorber, (a), (c) and (e) show the amplitude values before correction. Additionally, the expected OA amplitude values (dashed line), calculated according to Eq. (4) and using the μ values retrieved from the optical reference measurements, are plotted to verify the reliability of our measurements. (b), (d) and (f) present the amplitude spectra after spectral correction. All spectra were normalized to their respective mean value for better comparability with the reference. (For interpretation of the references to color in this legend, the reader is referred to the web version of the article.)

To verify the reliability of our measurements, the expected OA amplitude values were calculated according to Eq. (4), using the μ values retrieved from the optical reference measurements and plotted over the actual measurements: see Fig. 4a, c and e. The expected values coincide well with the measured amplitudes, with a slight underestimation for λ ≤ 720 nm. In this range, the OA measurements provide a better estimation of μ than the optical reference measurements, which is apparent in the excellent spectral correction for λ ≤ 720 nm in Fig. 4b, d and f. The underestimation of the μ values in the reference measurements (see Fig. 3) is due to the higher scattering, that reduces the accuracy of the spectrometric analysis presented in the supplementary material.

Our goal is to correct for this distortion, using the relation given in Eq. (5). Having determined μ at all wavelengths, we calculate the corresponding fluence Φ0(r|μ), following the steps (v–vi) described in Section 3.4. The corrected OA spectra are shown in Fig. 4b, d and f for all three absorbers, respectively. The shapes of the corrected spectra of the first and second absorbers (Fig. 4b and d) show a very good agreement with the reference absorption spectrum of the dye-solution. As for the third absorber, despite the observed systematic underestimation of the bulk effective attenuation (see Fig. 3c), the spectral correction is still apparent.

In practical applications such as the determination of local blood oxygen saturation levels, we are interested in the ratio of OA amplitudes at different wavelengths [13], [6], rather than in absolute values of . The reason is the incapability to measure absolute pressure amplitudes owing to the a priori unknown influence of the limitations of the imaging system (e.g. bandwidth limitation, restricted detection aperture) on the OA signal amplitude [28]. Even though these system inherent parameters are not known a priori, they can be considered to be independent of the irradiation wavelength.

Therefore, to provide quantitative evidence of the efficiency of our spectral correction, we are interested in the ratio ROA(r)/R(r). is the OA signal amplitude ratio and is the absorption coefficient ratio, both calculated for different wavelength pairs . The amplitude ratios are affected by two types of error: (i) the distortion of ROA relative to the reference spectrum, caused by the wavelength-dependent attenuation (akin to a systematic error) and (ii) a random error owing to measurement noise, which leads to an uncertainty in the experimental data. ROA(r)/R(r) can be considered as a random variable, which, in the case of an ideal spectral correction, would be uniformly distributed, whereby each pair yields a ratio equal to 1. To illustrate our purpose, Fig. 5 displays the histogram of realizations of the ratio ROA(r)/R(r) obtained from our experimental data at r, before and after spectral correction. These histograms are restricted to spectral regions where the distortion is noteworthy, e.g., 700 nm ≤ λ1 ≤ 725 nm and 740 nm ≤ λ2 ≤ 810 nm. Before spectral correction, the ROA(r)/R(r) values are significantly larger than 1 and the realizations are strongly scattered, as a result of the systematic error, which prevails over the random error. Yet, after spectral correction, the systematic error vanishes and the random error becomes dominant. The ensuing histogram follows a normal distribution with a mean equal to 1 and a standard deviation of 0.15. This, together with the fact that the two histograms do not show any overlap, underlines the relevance of the spectral correction. For our following analysis in the next section, we consider the ratio ROA(r)/R(r) for one particular wavelength pair . The corresponding ratios are shown for all three absorbers in Table 2.

Fig. 5

Histogram of realizations of the ratio of ROA(r)/R(r) for various wavelength pairs , before (blue) and after (red) spectral correction; were restricted to a spectral region with significant spectral distortion, e.g., 700 nm ≤ λ1 ≤ 725 nm and 740 nm ≤ λ2 ≤ 810 nm; after spectral correction, the histogram is dominated by the random error and can be approximated by a normal distribution with a mean value of 1 and a standard deviation of 0.15. (For interpretation of the references to color in this legend, the reader is referred to the web version of the article.)

Table 2

Quantitative comparison of ROA(r)/R(r) of each individual absorber at one wavelength pair, , before and after spectral correction.

jth absorber	Correction [before/after]	ROA(rj) [S800(rj)/S715(rj)]	R_ref(r_j)[μa,800chr(r_j)/μa,715chr(r_j)]	R_OA(r_j)R_ref(r_j)±σ
j = 1	Before	1.73 ± 0.15	2.65 ± 0.03	0.65 ± 0.06
	After	2.51 ± 0.27		0.94 ± 0.1
j = 2	Before	1.52 ± 0.18	2.65 ± 0.03	0.57 ± 0.07
	After	2.57 ± 0.57		0.97 ± 0.19
j = 3	Before	1.55 ± 0.38	2.65 ± 0.03	0.58 ± 0.14
	After	2.2 ± 0.88		0.84 ± 0.33

Discussion and conclusion

Validity of the tissue model

Overall, our results demonstrate that we have successfully applied the concept of multiple optical irradiations to correct for the spectral attenuation of the bulk medium. The optoacoustically determined μ values, based on the first and second absorbers, agree very well with the reference measurements, showing a variation in μ from 0.04 mm−1 at 700 nm to 0.11 mm−1 at 830 nm. These values lay within the physiological range of superficial human tissue, such as breast tissue (∼0.1 mm−1) or muscle tissue (∼0.2 mm−1) [34], which are at the site of potential target applications.3

Yet, for λ> 730 nm, the μ values obtained from the third absorber are consistently underestimated with deviations reaching −0.02 mm−1 with respect to the reference measurements. An explanation for this underestimation is the small spacing in the z-direction (∼2.5 mm) between the second and third absorbers. This geometry was purposefully chosen to demonstrate the fact that the fluence at the location of deeper absorbers is perturbed by the presence of absorbing structures located above (this models an in vivo setting, where, for instance, the trajectories of two or more blood vessels/arteries are bound to intertwine): see preliminary studies reported in [35]. In other terms, the absorbers do disturb to some extent the fluence Φ0,(r|μ), on the contrary to what has been assumed in Section 2. This influence can be seen in the different slopes of the OA amplitudes of the three absorbers in Fig. 2. For the exemplarily chosen wavelength, the slope of the third absorber is less steep than that of the first and second absorbers. Thus, the fitting of the diffusion model to the OA amplitudes generated by the third absorber yields μ values, which differ from the expected ones. We have further validated this explanation with additional MC simulations: see Fig. 6, showing a “close-up” of Fig. 2. It is apparent that the light diffusion model based on the optical reference measurements is appropriate to describe the fluence at the depth of the third absorber (z = 28 mm), in case of a single embedded absorber. Yet, in the presence of two additional absorbers located at z < 28 mm, the fluence is significantly reduced at r, depending on the number of absorbers present. This effect is particularly pronounced for 32 mm ≤d ≤ 35 mm (see Fig. 6, squares and dots), when the irradiation location is positioned close to the center, i.e. when the source-detector distance is the shortest. Under such conditions, the proposed diffusion model fails at correctly describing the fluence at the position r.

Fig. 6

MC simulations were performed at 760 nm to demonstrate the effect of the presence of additional absorbers on the fluence at the position of the third absorber r. It can be seen that the light diffusion model (solid red line) based on the optical reference measurements is inappropriate to describe the fluence at r. (For interpretation of the references to color in this legend, the reader is referred to the web version of the article.)

Ultimately, a few more general remarks can be made regarding the quality of the proposed method to optically characterize the bulk medium. We have quantitatively assessed the accuracy of our results by comparison to independently performed reference measurements of the absorption coefficient μ(r), the scattering coefficient μ(r) and the anisotropy factor g(r), needed to express μ. Naturally, these reference values are to be considered with precaution, insofar as the factor g(r) was not directly measured but estimated based on indirect measurements and different literature sources. Nonetheless, we have noticed that all estimations agree on the spectral trend of g(r). This means that the spectral behavior of μ remains unaffected by the absolute value of g(r), as the latter simply adds a systematic offset, which was taken into account in our error propagation calculations. In comparison, the source-detection distance d is of higher influence on the accuracy of the μ determination. For instance, a change in d by 2 mm results in a variation in the optoacoustically determined μ of up to 20%. This error, however, can be avoided with an appropriate geometrical calibration of the position of the optical fiber relative to the imaging plane. In conclusion, the retrieved μ values show a very good agreement with the independently performed optical reference measurements and thus, underline the general ability of the proposed technique to determine absolute μ values of strongly scattering bulk materials.

Influence of the number of irradiation spot positions and number of acquisitions on the estimation of the optical effective attenuation

The high number (M = 50) of irradiation spot positions, x, for this proof-of-principle study was chosen so as to fully understand the behavior of the light fluence for a large number of d and over a broad spectral range. To provide a good SNR, particularly for x far off the absorbers, 100 OA acquisitions per x were used for averaging. The resulting measurement time per x of around 3 s (including 100 OA acquisitions, image reconstruction and mechanical translation of the irradiation spot position) led to long acquisition times for a single experiment.

With the aim of shortening the acquisition time, we investigated how the accuracy of the estimated μ values is affected by reducing (i) the number of x and (ii) the number of acquisitions per x. Exemplarily, for one absorber (j = 2) and one set of measurements (λ = 760 nm), we estimated μ based on smaller subsets with n number of x values, randomly sampled (×100) out of M. The number of acquisitions for each x was kept constant. The resulting mean and 95% CI of μ for different n values are summarized in Fig. 7a. The results show that for this specific absorber and wavelength, the estimated μ values converge to the reference value for n ≥ 10, accompanied by a decrease in the fit uncertainty.

Fig. 7

(a) Mean and 95% CI of μ, estimated based on (a) randomly sampled subsets of irradiation spot positions with size n out of M, (b) by averaging over k acquisitions for M irradiation spot positions; all results are based on the second absorber (j = 2) acquired at λ = 760 nm.

For the same absorber (j = 2) and set of measurements (λ = 760 nm), we subsequently estimated μ by averaging k acquisitions for M irradiation spot positions. The resulting mean and 95% CI of μ for different k are presented in Fig. 7b. The results show that for the specific absorber and wavelength, the estimated μ values converge to the reference value for k ≥ 15, accompanied by a slightly decreasing fit uncertainty. The underestimation for small k results from the noise contribution to the fit, observed for x positions far off the absorbers. Excluding those noise dominated measurements (i.e., effectively reducing the number of x values) from the fit improves the μ estimation, even for k ≤ 15.

Certainly, these results suggest that n ≥ 10 positions and k ≥ 15 acquisitions would suffice to accurately determine μ, leading in practice to a significant reduction in acquisition time. Yet, it is difficult to draw a universal conclusion with respect to clinical applications, since the slope will be perturbed in an in vivo setting, where not only geometrical boundaries, but also optical heterogeneities play an influential role. Actually, the effect of the optical heterogeneities can already be observed here in the case of the third absorber: see Fig. 6. As such, we can only claim that there is trade-off between the number of irradiation spot positions and the number of acquisitions, and that this needs to be further investigated with in vivo studies.

Accuracy of the spectral correction and ensuing blood oxygen saturation level calculations

We shall recall that the main purpose of this “proof-of-principle” study was to correct for the spectral distortion of the measured OA spectra. The key point was to calculate Φ0,(r|μ), based on the spectrally determined μ using multiple-irradiation sensing, and to subsequently assess the accuracy of the spectral correction. The latter is of particular relevance if quantitative OA imaging is applied to measure physiological parameters such as the blood oxygen saturation level [1].

Compared to the distorted OA spectra, the corrected ones presented in Fig. 4(b), (d), and (f), indicate a significantly better agreement with the reference absorption spectrum of the ICG solution. This qualitative observation can be quantitatively judged by the ratio of ROA(r)/R(r) at 800 nm and 715 nm given in Table 2. After spectral correction, ROA(r)/R(r) of the first (0.93 ± 0.11) and second (0.94 ± 0.2) absorbers agree, within the standard deviation, with the target value (i.e. ROA(r)/R(r) = 1). Yet, the correction of the third absorber (0.84 ± 0.4) does not yield the intended ratio, despite a significant correction compared to the initial value (0.58 ± 0.087). This is related to the observed systematic underestimation of μ, as thoroughly explained above.

The choice of ROA(r)/R(r) was motivated by the target application example, i.e. quantification of the relative concentrations of deoxyhemoglobin (HHb) and oxyhemoglobin (HbO2) based on the distinct optical absorption ratios of HHb and HbO2 at two different wavelengths. Assuming HHb and HbO2 to be the main optically absorbing blood constituents, their relative concentrations at r can be retrieved using the following relation between the known molar absorption coefficient and the resulting OA amplitudes at two different wavelengths:where the product of the total hemoglobin (Hb) concentration (∼150 g Hb/L) and the wavelength-dependent molar extinction coefficient is given as ϵ. The relative concentrations of HHb and HbO2 is denoted by cHHb and cHbO, respectively and the factors κ depend on the wavelength-dependent fluence, as well as on various system parameters. (Although, S, μ, cHbO, and cHHb are all a function of r, this dependency is omitted in Eqs. (6), (7), (8), for sake of simplicity.) The inversion of Eq. (6) yieldsin which det is the determinant of the molar extinction coefficient matrix and refers to the ratio between OA signal amplitudes at two different wavelengths. The resulting blood oxygen saturation (SO2) level can be expressed as a function of the relative concentrations of HHb and HbO2 [36], where the ratio between κ and κ accounts for the wavelength-dependency of the fluence at the two wavelengths of interest. The accuracy of the SO2 level quantification depends on the accurate determination of this ratio, i.e. the spectral correction of the amplitude ratio ROA(r):

To discuss the accuracy of the achieved spectral correction in the context of SO2 level quantification, we extrapolate the measured OA spectra (see Fig. 4) to a literature-based blood absorption spectrum, corresponding to a nominal SO2 level of 90%, which is reasonable for a healthy adult [36] (black solid line in Fig. 8). The extrapolated spectra were obtained by multiplying the ratios of the measured OA spectra (see Fig. 4) with the said nominal blood spectrum. Fig. 8 shows the extrapolated spectra for the second absorber (blue dots), before and after spectral correction, together with the literature-based nominal blood absorption spectrum. The SO2 levels were quantified for the particular wavelength pair using Eq. (8). It should be noted that this wavelength pair, which was already considered in Section 4, was chosen so that the efficiency of our approach could be assessed in the case of a strong initial distortion of the OA spectrum. To underline the sensitivity of the SO2 level determination to small deviations from the extrapolated spectrum, the absorption spectra corresponding to the two “extreme” blood oxygenation states (1% O2 and 99% O2 saturation) are added in Fig. 8.

Fig. 8

OA spectrum of the second absorber (blue dots), before (a) and after (b) spectral correction, extrapolated to a literature-based nominal blood absorption spectrum of 90% SO2 (black solid line); the SO2 level estimation is based on the extrapolated OA spectra at the two wavelengths, λ1 = 800 nm and λ2 = 715 nm; the absorption spectra of oxygenated (99% SO2) and deoxygenated (1% SO2) blood illustrate the full possible absorption range. (For interpretation of the references to color in this legend and in the text, the reader is referred to the web version of the article.)

Table 3 presents the nominal and the estimated SO2 levels, before and after spectral correction, together with the experimentally obtained 95% CI values. The estimated SO2 values before spectral correction show for all three absorbers a significant deviation from the reference value, which leads in practical applications to an erroneous diagnostic outcome. Yet, after spectral correction, the estimated SO2 values approach the nominal value. The corresponding CI values increase after spectral correction, as a consequence of the error in the experimental μ determination. It is apparent that despite the increased CI values, the improvement in SO2 level determination is significant, underlining the power of this technique. Even in the case of a strong initial distortion of the OA spectrum, the accuracy of our results with the first absorber lies within the CI, which is comparable not only to other OA based methods [6], but also to commercially available pulse oximeters (95% CI <±6%) [37], [38]. Moreover, our approach does not require any a priori knowledge about the tissue geometrical properties, unlike the study of Laufer et al. [6], and provides spatially resolved information, not possible with pulse oximetry. Furthermore, we would like to point out that the accuracy of our results could further be increased when using a larger wavelength range.

Table 3

The nominal and estimated SO2 levels based on all three absorbers, before and after spectral correction, together with the corresponding 95% CI.

jth absorber	Correction [before/after]	Estimated SO2 [%]	95% CI [%]	Nominal SO2 [%]
j=1	Before	69.5	±6.4	90
	After	87.2	± 7.3
j=2	Before	62.1	±8.6	90
	After	88.9	± 16.4
j=3	Before	62.8	±18.9	90
	After	82.2	± 26.8

In conclusion, we presented a method to correct for the wavelength-dependent optical attenuation in highly scattering media such as biological tissues, which is of great relevance in the frame of quantitative OA imaging. The method is based on the concept of multiple irradiation sensing which was implemented on a system designed for handheld, clinical combined OA and US imaging. The fact that the proposed technique is non-invasive, does not require a priori knowledge or any other additional detection modalities and has low computational cost, makes it promising for real-time clinical application. The results of this feasibility study emphasize its applicability also for deep optoacoustic imaging. Its benefit for the estimation of blood oxygen saturation levels was underlined. More importantly, this proof-of-principle study helped identifying the challenges and the limitations of this approach. These include the applicability limits of the tissue model in the presence of several strong optical absorbers, as well as the influence of the noise on the fit of the diffusion model. Future steps will involve more realistic situations, i.e., simulation and experimental work on strongly heterogeneous phantoms (multi-layered, non-trivial vasculature model, etc.) with influential boundaries, and the corresponding adaptation of the light propagation model. This will pave the way to quantitative OA imaging in a clinical settings.

Conflict of interest

The authors declare that there are no conflicts of interest.