Reflection-mode multiple-illumination photoacoustic sensing to estimate optical properties.

OBJECTIVES: We analyze a reflection-mode multiple-illumination photoacoustic method which allows us to estimate optical scattering properties of turbid media based on fitting light-transport models and explore its limits in optical property estimation and depth-dependent fluence compensation.
BACKGROUND: Recent simulation results show significant promise for a technique called multiple-illumination photoacoustic tomography (MI-PAT) to quantitatively reconstruct both absorption and scattering heterogeneities in turbid medium. Prior to experiments, it is essential to develop and analyze a measurement technique and probe capabilities of quantitative measurements that focus on sensing rather than imaging.
METHODS: This technique involved translation of a 532 nm pulsed-laser light spot while focusing an ultrasound receiver on a sub-surface optical absorber immersed in a scattering medium at 3, 4 and 5 mm below the surface. Measured photoacoustic amplitudes for media with different reduced scattering coefficients are fitted with a light propagation model to estimate optical properties.
RESULTS: When the absorber was located at 5 mm below the membrane in media with a reduced scattering coefficient of 4.4 and 5.5 cm(-1), the true values were predicted with an error of 5.7% and 12.7%, respectively. We observe accuracy and the ability of estimating optical scattering properties decreased with the increased reduced scattering coefficient. Nevertheless, the estimated parameters were sufficient for demonstrating depth-dependent fluence compensation for improved quantitation in photoacoustic imaging.

Introduction

Recently, our group proposed the concept of multiple-illumination photoacoustic tomography (MI-PAT), where optical illumination patterns on the surface serve as sources and signals from sub-surface optical absorbers act as virtual sub-surface detectors [1,2]. Recent simulations show significant promise for using data from combinations of optical source and photoacoustic virtual detector pairs to reconstruct both absorption and scattering heterogeneities [3]. However, no clear technique for experimentally obtaining and fitting such data was presented. In particular, reflection-mode imaging modalities are desirable due to ease of accessibility in thick tissues, and for implementing multiple illumination photoacoustic microscopy (MI-pan class="Gene">PAM) where penetration depths are within the radiative transport regime. Part of the experimental challenge is to design a probe capable of flexible illumination-pattern delivery (including point-illuminations) and co-linear reflection-mode ultrasound detection with maximal signal-to-noise ratio. As proposed MI-PAT/MI-PAM reconstruction strategies rely on forward light propagation models, the probe should have the capability to sense accurate fluence as predicted by models of light transport at varying depths and varying illumination locations. Accurate estimation of optical properties is required to calculate accurate fluence distribution.

Traditional optical property estimation techniques are limited to surface measurements [4-7]. The majority of current photoacoustic imaging applications rely on optical absorption contrast and much less work, however, has been done on to estimate optical scattering properties [8]. In photoacoustic imaging, we have an ability to effectively measure the relative local laser fluence at sub-surface locations. As such, photoacoustic methods may potentially provide information in addition to surface detectors. Some groups have tried to obtain optical properties of biological tissues using photoacoustic waveform analysis in simple well-controlled situations [9-11]. We recently developed an experimental photoacoustic method to effectively estimate the Green's function for light transport in a turbid medium [12]. We translated a pulsed light spot while receive focusing on a fixed sub-surface absorber to obtain a curve of photoacoustic signal amplitude as a function of incremental light translocation distance. We demonstrated that this curve, normalized by its peak value was distinct for different levels of optical turbidity.

We felt it essential to analyze the quantitative capabilities of above method in homogeneous turbid media and to this end, we test its capabilities for sensing (rather than imaging) of bulk tissue scattering coefficients and demonstrate highly accurate agreement between forward light-propagation models and experimental data. Recovering optical properties from sensing data and understanding the capabilities and limitations of the method will be critical to the success of future MI-PAT/MI-PAM quantitative imaging strategies, and we feel is worthwhile of a dedicated study. While similar to our previous work (which demonstrated scattering differences could be qualitatively observed from photoacoustic data), the present work aims to be much more quantitative and explores the capabilities and limitations of the technique for optical property estimation and depth-dependent fluence compensation. This study provides strong feasibility data for the utility of the proposed probe for future reflection-mode pan class="Disease">MI-PAT/MI-PAM imaging studies.

Method

Experiment

We used a custom built probe for our study [12]. This probe consists of a 10 mm × 40 mm right angle prism (P), optical index-matching fluid (IML), and an ultrasound transducer (UT) (Fig. 1). The index matching liquid (Cat# 19569, Cargille-Labs, Cedar Grove, NJ, USA) was kept underneath the prism by attaching a ∼25 μm-thick plastic (Saran wrap) membrane to an acrylic holder which was attached to a XYZ positioning stage. The index-matching fluid which has a reflective index similar to the fused silica prism (n = 1.46), enabled light to be delivered down to the sample and collected from the sample without significant loss, refraction, or inter-medium reflection. For the present study, we selected a single-element focused ultrasound transducer (10 MHz, f = 19.05 mm, F#3.1, spherical focus G series immersion type transducer, pan class="Disease">CD International Technology, Santa Clara, CA, USA) which was mounted horizontally to capture photoacoustic signals deflected from index matching fluid and prism interface. In future MI-PAT/MI-PAM imaging studies, we plan to substitute this single element transducer by a linear array transducer. We used a laser (SureLite III, Continuum Inc., Santa Clara, CA, USA) with a pulse repetition rate of 10 Hz. For simplicity, all measurements were done using 532-nm light. We place a thin glass slide in the beam path to reflect a small amount of light onto a high speed photo diode to capture laser pulse intensity variation. A small light spot was focused at the membrane interface by translating the probe vertically relative to lens L. A human hair was used as an optically absorbing target, and was mounted on an acrylic holder submerged in the medium with a 25 mm × 55 mm opening. Initially, the absorber was mounted on a XZ positioning stage and immersed in a water bath for alignment purposes. Our pulser-receiver (5703PR, Panametrics, Waltham, MA, USA) was set to pulse-echo mode and the ultrasound signal of the human hair was captured by a digital oscilloscope (DPO 7054, Tektronix, Beaverton, OR, USA). While observing the ultrasound signal, the hair sample was translated along the z-axis to position the hair at the required depth below the membrane. This depth was a parameter which was varied in our experiments.

Fig. 1

(a) Schematic and (b) photograph measurement configuration of photoacoustic probe and light delivery. The mirror prism (ML) and the focusing lens (L) in the dashed box were translated along the x-axis to move the light spot. PL: pulsed laser input; G: glass slide; PD: photo diode; UT: ultrasound transducer; IML: index matching liquid; M: membrane; IL: intralipid; D: absorbing dye in a thin tube.

To mimic scattering of human tissue, we diluted 20% Intralipid in water to obtain reduced scattering coefficients of . Note that we performed independent measurements of the reduced scattering coefficient and absorption coefficient of our Intralipid stock-solution using the Oblique-Incidence Diffuse-Reflectance (OIR) technique [5]. We set the absorber at depths of 3, 4 and 5 mm below the membrane for each Intralipid concentration. The light spot was translated relative to the absorber by translating top prism (MP) and focusing lens (L) horizontally without moving other parts. The photoacoustic signal by the absorber was recorded by the digital oscilloscope. We plotted photoacoustic signal vs. light spot translation distance. This curve is effectively a measure of a segment of the Green's function of radiative light transport. To minimize the effect of source intensity variation, each photoacoustic signals is normalized by a corresponding photodiode signal. To demonstrate the ability to compensate for depth-dependent fluence for quantitative imaging, we created a phantom consisting of a silicone tube (ID = 0.8 mm; OD = 1.5 mm) filled with varying concentrations of Crystal Violet dye, immersed in a diluted Intralipid bath. For a selected Intralipid concentration, we fixed the absorber depth and obtain photoacoustic signal for different concentrations of dye. Note that the tube was flush with water before and after filling it with the dye. We used a fixed focus transducer for this study and its response varies with the absorber depth. We obtained the amplitude response of the transducer along the axial direction by measuring ultrasound signals of a human hair immersed in the water at different depths. The amplitudes of this depth-dependent ultrasound signal were applied to normalize the photoacoustic signals at different depths. The multiple-illumination photoacoustic sensing technique was used to obtain the optical properties for the Intralipid concentration. The results were applied in a forward Monte Carlo simulation to estimate the laser fluence at the absorber to normalize photoacoustic signals.

Optical property estimation

Our data was fit to light transport models. First, we considered Monte Carlo simulations because it accurately estimates fluence rate at any given location for a medium with any given optical properties. We used a steady state Monte Carlo simulation program [13] to estimate internal fluence in a semi-infinite medium. The incident beam was narrowly focused with divergence an angle of 3°. Prior to the photoacoustic measurements, we conducted a knife edge experiment and found beam diameter at the membrane was 125 μm. Refractive index of incident medium (IML) and semi-infinite medium were set to 1.46 and 1.33 respectively. Simulations were done for various reduced scattering coefficients for a fixed absorption coefficient obtained using OIR. Gaussian distributed 100million photon packets at the surface of the membrane were launched for each case. We use cylindrical coordinates and captured the lateral fluence rate distribution at the absorber depth. The normalized fluence rate was calculated by normalizing above curve by its maximum value.

Next, we applied the diffusion approximation where the transport mean path is much smaller than the absorber depth. The diffusion constant, D is given by,

In the diffusion approximation, the fluence rate distribution at a point generated by an isotropic point source at a point in an infinite homogenous medium is given by,where and . When we consider a narrow pencil beam on a semi-infinite medium, the fluence rate at the point (r, z) is given by the linear combination of positive and negative point sources [14]. Considering a refractive-index-matched boundary condition:where , , is the inverse transport mean free path. The normalized fluence distribution at the absorber layer can be obtained by dividing Eq. (3) by the maximum fluence rate.where and . It is important to note that, the decay of Eq. (4) does not purely depend on .

The amplitude of the photoacoustic signal (PA) at is given by

where , and are Gruneisen parameter, the laser fluence at depth z and the absorption coefficient of the absorber respectively. At constant temperature and wavelength, the Gruneisen parameter is constant. For a fixed absorber, we can write a similar relationship for maximum photoacoustic signal.

When Eq. (5) is divided Eq. (6), we can obtain a direct relationship between normalized photoacoustic signal and normalized fluence. We obtain curves of for various and combinations and compared with the normalized photoacoustic amplitude results. Each data set is compared with the corresponding experimental data and calculated the mean square error.

Results and discussion

Using experimental methods described above, we first calculated the normalized photoacoustic amplitude data from experimental measurements. Fig. 2(a) shows the normalized photoacoustic amplitude as a function of light beam translocation distance for constant . Measurements were taken when the optically absorbing target (a hair) was at 3, 4 and 5 mm below the membrane, where such distances were measured via time-of-flight of acoustic/photoacoustic signals. We next kept the absorber at a predefined depth (4 mm) and changed the Intralipid concentration to achieve reduced scattering coefficients of 4.4, 5.5 and 11 cm−1 (Fig. 2(b)). We fitted above data with Monte-Carlo simulation results. In Monte Carlo simulations, we compute the fluence distribution in the r–z plane by applying input optical parameters ( and ) that were measured from the OIR technique. The experimental results show good agreement with the simulation results.

Fig. 2

Normalized photoacoustic signal vs. lateral translocation distance results for Intralipid. Experimental results (symbols) and corresponding Monte Carlo simulation results (lines) are given for (a) normalized PA signal of the absorber located at various depths (3, 4 and 5 mm below the membrane) for single Intralipid concentration (11 cm−1), (b) normalized PA signal of the absorber located 4 mm below the membrane for various Intralipid concentrations .

To estimate optical parameters from our photoacoustic measurements, we performed an iterative procedure as follows: for a given (measured using OIR) absorption coefficient, we generated Monte Carlo simulation results for different reduced scattering coefficient values. We then calculated the root mean square (RMS) error between experimental and simulation results at each experimental data point. The average RMS error was calculated by averaging all mean square error values at each lateral location on the curves. The value of for which the mean-square error is minimum is taken as the best estimate. For commical">parison purposes, the actual reduced scattering coefficient is marked as a dashed vertical line (Fig. 3). The minimum pan class="Disease">RMS error occurred at for the absorber at 4 and 5 mm below the membrane, respectively for true . It differs by less than 5.7% from the actual value. However, we found a larger variation for true . When the absorber was located at 4 and 5 mm below the membrane, the minimum errors observed at respectively. The true value was predicted with an error of 4.5% and 12.7%, respectively. At the highest-investigated reduced scattering coefficient the minimum in the MSE-curves (Fig. 3(c)) is difficult to see because the normalized fluence profiles become increasingly indistinguishable. It should also be noted that minor variations in the simulation significantly degrade the relative importance of distant data points. This combination is likely to limit the occurrence of global minima in the error profile. For this study, we assumed that the absorption coefficient of the semi-infinite media is a known parameter. But, it will be an unknown parameter in practical situations and we will have to simulate Monte Carlo simulation for different combination of absorption and reduced scattering coefficients. In that case, root mean square error plot will become a three dimensional curve. One independent axis will represent reduced scattering coefficient and the other axis will represent absorption coefficient. Third axis will be the RMS error. The global minimum of the three dimensional curve will provide the estimated optical parameters.

Fig. 3

Root mean square error between experimental results and Monte Carlo simulation for Intralipid with (a) 4.4 cm−1, (b) 5.5 cm−1 and (c) 11 cm−1. The absorber was immersed at d = 4 and 5 mm. The vertical dashed line represents the true .

Even though the Monte Carlo simulation is computationally expensive, it provides an accurate fluence rate at any given depth for any optical parameter combinations. On the other hand, diffusion approximation that gives solutions quickly, is accurate only when a given depth is longer than multiple transport mean free path and is much smaller than [13]. We check the feasibility of applying the diffusion approximation to extract optical properties from normalized fluence profiles for because transport mean free path is much smaller than the absorber depths for that case. Here, the optical properties extracted from the OIR-method were applied in Eq. (4). The experimental results show a good fit with the diffusion approximation results in Fig. 4(a). Similar to Monte Carlo simulations, we calculated RMS error for various reduced scattering coefficient values. We observed the global minima at when the absorber was at 3, 4 and 5 mm below the membrane respectively (Fig. 4(b)). All three values are less than the actual value and maximum error is as high as 25%. Unlike Monte Carlo simulation, diffusion approximation did not show multiple local minima in pan class="Disease">RMS error curves.

Fig. 4

(a) Experimental results (symbols) and diffusion approximation results (lines). Normalized PA signal of the absorber located at various depths (3, 4 and 5 mm below the membrane) for single Intralipid concentration (11 cm−1). (b) Root mean square error between experimental and diffusion approximation results.

Above results show that this method can estimate accurately with minimal error for smaller and accuracy decreases with increasing . To understand this in detail, we generated noisy normalized photoacoustic amplitude data using Eq. (4) and added normal random noise (mean = 0; std. dev. = 0.1%) for an absorber located 5 mm below the membrane. The purpose of this theoretical simulation is to check the capabilities and limitations of this model to recover optical parameters from noisy photoacoustic amplitude data. Fig. 5(a) shows the simulated noisy data and the best fit curve obtained by the non linear fitting function ‘nlinfit’ in MATLAB (Mathworks Inc.) for . Fig. 5(b) shows estimated and for true values of 4, 7, 10, and 13 cm−1 for noisy data. We observed large errors on estimated for larger true . On the other hand, we witnessed large errors on estimated for smaller true . The reason of this behavior is likely explained in simulated curves for different true values (Fig. 5(c)). As increases, the curves are getting too close to each other and estimating accurate become harder for noisy data. The decay of these curves that partially depend on have direct relationship with are almost similar for smaller . That increases the errors of for smaller . We found that depending on the initial guess, the non-linear fitting procedure was quite noise-sensitive and may end up in local minima.

Fig. 5

(a) Simulated noisy normalized photoacoustic amplitude data set (red) and its best fit curve per Eq. (4) (black) for a medium with and . Absorber was located at 5 mm blow the membrane. (b) Estimated and data for true . was kept constant at 0.1 cm−1. (c) Simulated curves for different .

Despite these difficulties, we also examine the ability to compensate for depth-dependent fluence so that effective local absorption can be sensed without appreciable bias due to depth, an important requirement of future MI-PAT/MI-PAM imaging techniques. Fig. 6 shows measured photoacoustic amplitudes from dye-filled tubes at varying dye concentrations and for three different tube depths. Ideally, we would like photoacoustic images to convey information about optical absorption without depth-dependent fluence biasing. Using the multiple-illumination sensing technique described above, we estimated the optical parameters of the Intralipid, and fed this into a Monte Carlo forward model to predict fluence. Additionally accounting for the transducer response along the axial direction as described in the Methods section, we applied compensation factors to the measured data. After compensation, modified photoacoustic amplitudes from all depths show similar relationship with optical absorption as shown in Fig. 6(b). Results are commical">parable to previous studies using OIR to estimate optical parameters [15], however, in this case, we used the multiple-illumination photoacoustic sensing technique.

Fig. 6

(a) Photoacoustic signal amplitude obtained for a tube filled with dye having different absorption coefficients, immersed at different depths. (b) Modified amplitude after normalizing by the local fluence and the transducer response. Fluence was estimated using a Monte Carlo model with optical parameters estimated from the multiple-illumination sensing method described above. (Dashed lines depict best fit for each data set.)

While quantitative estimation of optical parameters with our technique may not be as noise-robust as other purely optical techniques, our data does suggest this technique is more suitable for recovering optical properties of a medium with small to moderate and it could serve as a calibration method for future MI-PAT/MI-PAM imaging studies. In vivo, blood vessels could serve as optically absorbing targets. It remains to be seen how well the proposed method might work in living subjects due to optical heterogeneity. Absorber depth is presently gaged from time-of-flight of photoacosutic signals and presumes sound-speed heterogeneity. Future work will introduce optical absorption and scattering heterogeneities, and apply inverse-algorithms, presently under development, to reconstruct images of these optical heterogeneities, both in the transport and diffusion-regimes. Future work may also need to account for acoustic heterogeneities and account for frequency-dependent ultrasound attenuation.

Conclusion

We describe a reflection-mode multiple-illumination photoacoustic method which permits estimation of optical properties of turbid media based on fitting light-transport models. The method entails measuring the photoacoustic amplitude originating from a sub-surface absorber as a surface light spot is translated laterally. The measured curves effectively represent the Greens function of light transport. We applied Monte Carlo and diffusion approximation models to extract optical properties from experimental curves using fitting procedures. We plotted root mean square error between experimental and simulation results and selected the global minimum. The accuracy and the ability of estimating optical scattering properties decreased with the increased reduced scattering coefficient and estimating optical absorption properties decreased with the decreased reduced scattering coefficient. Nevertheless, the estimated parameters were sufficient for demonstrating depth-dependent fluence compensation for improved quantitation in photoacoustic imaging. Future work incorporating an imaging array transducer into multiple-illumination probe design should prove fruitful opportunities to implement multiple-illumination photoacoustic tomography and microscopy studies.

Conflict of interest statement

The authors declare that there are no conflicts of interest.