1.

Introduction

Soft-tissue organ imaging with near infrared (NIR) light is pursued vigorously because of its application in early diagnosis of cancer based on the measured spectral variation of ${\mu }_a$. This modality known as diffuse optical tomography (DOT) however produces only low resolution images limited by the diffusion propagation of photons through soft tissue with large ${\mu }_s^{\prime }$s (reduced scattering coefficient). Ultrasound assisted optical tomography (UAOT) is developed as a possible remedy to this poor resolution, which takes advantage of tight focusing of ultrasound (US) waves in tissues and restricting imaging of optical as well as mechanical properties to the focal volume (the region of interest, ROI) using the so-called US-tagged photons. UAOT began by constructing qualitative images of ${\mu }_a$ in the ROI from the measured $M$ of light autocorrelation $G_1(\tau )$. There were several attempts to move onto quantitative ${\mu }_a$ and ${\mu }_s$ recovery^1–4 which began with modeling of US-induced effects in the ROI which are oscillation of scattering centers and refraction index modulation ($\mathrm{\Delta }n$).^5–7 Theoretical expressions are derived for $M$ in terms of the phase modulation ($\varphi$) and its fluctuation ($\mathrm{\Delta }\varphi$) picked-up by light from the ROI which connected $M$ to ${\mu }_a$, ${\mu }_s^{\prime }$ and also to $E$. In spite of the above there were no attempts of quantitative recovery of these properties from $M$. This work is intended to primarily fill this gap by quantitative recovery of ${\mu }_a$, ${\mu }_s$ from $\langle \mathrm{\Delta }\varphi \rangle$ and $E$ from $M$. $\mathrm{\Delta }\varphi$, as demonstrated here, has a non-zero mean $\langle \mathrm{\Delta }\varphi \rangle$ when the US frequency ($\omega$) is small ($<\mathrm{kHz}$) and the anisotropy factor $g$ is large. The ultrasound modulation can be decreased when the optical scattering increases, thereby reducing the modulation depth of UAOT signals.³

Both $M$ and $\langle \mathrm{\Delta }\varphi \rangle$ can be theoretically computed from $G_1(\tau )$ the amplitude autocorrelation of detected photons arrived at using the expressions given in Refs. 8 and 9, also Monte Carlo simulation of detected scattered photons. As demonstrated earlier⁹ $\langle \mathrm{\Delta }\varphi \rangle$ computed using MC simulations is shown to be non-zero for the phantom with large $g$ (here $g=0.89$) when $\omega <1\mathrm{kHz}$ and the insonified volume ($V_{\text{in}}$) does not exceed 3 to 4 times ${({\ell }^{*})}^3$ (${\ell }^{*}$ is the transport mean path given by $1/{\mu }_s^{\prime }$). In this work, because we use $\langle \varphi \rangle$ to compute ${\mu }_a$ and ${\mu }_s$ changes we would like to keep $\langle \varphi \rangle$ large so that with ${\mu }_a$ and ${\mu }_s$ there is a large range of variation in $\langle \varphi \rangle$. For this we keep in our simulations and experiments $V_{\text{in}}\approx {({\ell }^{*})}^3$. With reference to Fig. 1 the phase fluctuation ($\mathrm{\Delta }\varphi$) picked up by light in a typical scattering event owing to a displacement of the scatterer by $\mathbf{u}$ is $\mathbf{u}({\widehat{k}}_s-{\widehat{k}}_i)=|\mathbf{u}||{\widehat{k}}_s-{\widehat{k}}_i|\cos \frac{\theta }{2}$. Therefore decorrelation in $\mathrm{\Delta }\varphi$ will be dictated by $\langle \cos \frac{\theta }{2}\rangle$ (and not by $\langle \cos \theta \rangle$) and therefore the anisotropy factor, when considering this decorrelation, should be replaced by $g_n=\langle \cos \frac{\theta }{2}\rangle$ and ${\ell }^{*}$ used above is arrived from this new $g_n$ (denoted it by ${\ell }_n^{*}$) which is larger than the usual ${\ell }^{*}$.

Fig. 1

A schematic diagram of scattering of light by a scattering center which undergoes vibration with an amplitude $\mathbf{u}$. The incident light with wave vector, ${\mathbf{k}}_{\mathbf{i}}$ is scattered with an associated wavelength along ${\mathbf{k}}_{\mathbf{s}}$.

First we verify that $\langle \varphi \rangle$ indeed can be non-zero. Theoretically we arrive at a distribution of $\varphi$ at an array of detector points on one face of a slab of dimension $40\times 40\times 50{\mathrm{mm}}^3$ and optical properties, $g=0.89$, ${\mu }_s=20{\mathrm{cm}}^{-1}$, and ${\mu }_a=0.18{\mathrm{cm}}^{-1}$ by launching ten million photons from a point on the opposite side.⁹ From the random distribution of $\varphi$ a histogram of its distribution is computed. It is seen that when $\omega =1\mathrm{MHz}$ the distribution is uniform with zero mean and when $\omega =100\mathrm{Hz}$ it has a mean of 130 deg. The figure is not included for want of space. The insonified region for the computation is taken as a cube of volume $4\times 7\times 4{\mathrm{mm}}^3$ which is $\sim {({\ell }_n^{*})}^3$. The histogram arrived at from the experimentally measured $\varphi$ is shown in Fig. 2(a) and 2(b) which also verifies almost the same non-zero $\langle \varphi \rangle$ obtained from computation (the experiments are described further down).

Fig. 2

(a) Experimental plot of histogram of $\mathrm{\Delta }\varphi$ at US beat frequency of 100 Hz. (b) Experimental plot of histogram of $\mathrm{\Delta }\varphi$ at US frequency of 1 MHz.

We use the histogram to compute $\langle \varphi \rangle$. To compute $M$, $G_1(\tau )$ is computed from the arrived photons at a detector which is photon-path probability density weighted (only those paths which intersect ROI are considered) and its Fourier transform amplitudes at $\omega ={\omega }_o$ the acoustic frequency is computed. The details are in Ref. 10. The computed as well as experimentally measured $M$ and $\langle \varphi \rangle$ are used to recover the average ${\mu }_a$, ${\mu }_s$ and $E$ ($\langle {\mu }_a\rangle$, $\langle {\mu }_s\rangle$ and $\langle E\rangle$, respectively) of the insonified ROI in the object, as detailed below.

We calibrate our measurements with respect to ${\mu }_a$, ${\mu }_s$ (against $\langle \varphi \rangle$) and $E$ (against $M$). Theoretically, to get plots of $M$ versus $E$ one has to go through a number of steps for a given $E$: 1. compute the US radiation force in the ROI, 2. set-up and solve the force-balance equation to find $\mathbf{u}(\mathbf{r})$ the distribution of US-induced displacement of the scatterers, and 3. transport photons using MC simulation and arrive at $G_1(\tau )$ and $M$. Details of the above steps are given in Ref. 10. While doing the above we can also compute the distribution of $\varphi$ and $\langle \varphi \rangle$. However, in regard to $\varphi$ we have an easier route which employs the semi-empirical relation first proposed in Ref. 11 in the context of frequency-modulation diffuse optical tomography (DOT) connecting $\langle \varphi \rangle$ and ${\mu }_a$, ${\mu }_s$, and $g$. It is given by

The corresponding ${\langle \mathrm{\Delta }\varphi \rangle }_{{\mu }_s}$ when ${\mu }_s$ changes is

Equation (1) is verified using data from MC simulation and found to be fairly close representation of the actual variation of $\langle \varphi \rangle$. The results are not given for want of space.

For extracting $E$ from $M$ we made use of the frequency at which a dominant natural mode of vibration of the ROI coincides with the US frequency. The resonance for a jelly like soft tissue mimicking object is at low frequencies and therefore to insonify the object at low frequencies we mix two US beams oscillating at slightly different (and adjustable) frequencies. The mixed difference frequency pressure wave drives the scattering centers in the ROI. The resonance can be observed in the $M$ versus $\omega$ plot as a peak, which is first obtained from computed $M$ and the experiments. At these low frequencies, as observed in Ref. 8, one of the effects of US forcing, namely $\mathrm{\Delta }n$ is small. Therefore, in this study, we neglected $\mathrm{\Delta }n$. The resonance frequency (${\omega }_o$) is used to compute $\langle E\rangle$ of the ROI using the method of bisection as detailed in Ref. 10.

2.

Methods and Materials

The experimental set-up in Fig. 3 is a parallel speckle detection system which uses a source-locked detection. The illumination is from a laser diode with wavelength at 785 nm, emitting 25% duty-cycle square waves which is synchronized to the oscillators driving the two US transducers and can be phase shifted with respect to the square-wave. The focusing US transducers operate at frequencies of 1 MHz and $1\mathrm{MHz}+\mathrm{\Delta }f$ where $\mathrm{\Delta }f$ can be varied from $\sim 50\mathrm{Hz}$ to 1  kHz and are driven by ultrastable oscillators. The object (phantom), which is either poly-vinly alcohol (PVA) or a milky liquid in a cuvette, is immersed in a water bath for acoustic impedance matching. The two US transducers are adjusted to have a bisecting focal region within the phantom. The region of intersection is the ROI which is driven by an acoustic radiation force at a frequency $\mathrm{\Delta }f$. To obtain $M$ and $\langle \mathrm{\Delta }\varphi \rangle$ one records four sequential speckle intensities using the CCD camera focused in the exiting light, one each at phase difference between the acoustic and laser drive signals of 0, 90, 180, and 270 deg. From the four intensities $M$ and $\varphi$ corresponding to each pixel in the CCD camera are computed.^8,9 The average $\langle \varphi \rangle$ is computed from the histogram of $\varphi$ and also the average of $M$. The histograms of the experimental measured $\varphi$ (Fig. 2) verifies the fact that $\langle \mathrm{\Delta }\varphi \rangle$ is non-zero when $\mathrm{\Delta }f(=\mathrm{\Delta }\omega /2\pi )<1\mathrm{kHz}$.

Fig. 3

(a) The schematic diagram of the experimental setup. (b) A typical cross section of the experimental set-up in which US transducers, diode laser and CCD camera are fixed. The CCD camera is kept opposite to the laser diode to capture the transmitted light through the sample. The axis through the diode laser, and the CCD detector bisects the angle subtended by the transducer axes (i.e., 60 deg) at the insonified ROI.

We have used two phantoms in the experiments. First is a PVA phantom whose optical and mechanical properties can be tailored to match the breast tissue and a liquid phantom. PVA slabs of dimensions $40\times 40\times 50{\mathrm{mm}}^3$ are made with ${\mu }_a$, ${\mu }_s$ tailored to be $0.18$, $20{\mathrm{cm}}^{-1}$, respectively, and $E$ to be 11.39 , 23.42, and 40.35  kPa in three different specimens. Experiments are conducted using them and $M$ versus $\mathrm{\Delta }f$ plots obtained are shown in Fig. 4 which are compared to the eigen mode distribution (natural frequency response) of the vibrating ROI computed using ANSYS package.¹² The prominent resonant mode matches well with the peak in $M$ versus $\mathrm{\Delta }f$ curve in all the three cases considered. From these measured resonant frequencies the value of $E$ of the ROI are computed as explained in Ref. 10 which are found to be 11.42 , 23.79 , and 40.78  kPa for the three slabs.

Fig. 4

(a) The plot of frequency response (amplitude of vibration) of the insonified region against the ultrasound beat frequency using ANSYS package. The elasticity co-efficients used are 11.39, 23.42, and 40.35 kPa. Poisson’s ratio of 0.499 and mass density of $1000\mathrm{kg}/{\mathrm{m}}^{-3}$ are used in all the simulations. (b) The experimental plot of the modulation depth (M) against the ultrasound beat frequency. The experimental error is in the range 0.46% to 7.80%. PVA phantoms with two, three, and four freeze-thaw cycles which correspond to the elastic co-efficients of 11.39, 23.42, and 40.35 kPa are used.

The second is a liquid inhomogeneity obtained by mixing milk, water, and drops of India ink. By proper mixing of the above we were able to vary ${\mu }_s$ from 2.85 to $8.24{\mathrm{cm}}^{-1}$ and ${\mu }_a$ from $3.2\times {10}^{-3}$ to $792.8\times {10}^{-3}{\mathrm{cm}}^{-1}$. These were independently verified using a transmitted intensity measurement as in Ref. 8. With the liquid in a cuvette we have conducted experiments to find mainly $\langle \varphi \rangle$. In the first set, fixing ${\mu }_s$ at $2.85{\mathrm{cm}}^{-1}$ we varied ${\mu }_a$ from $3.2\times {10}^{-3}$ to $792.8{\mathrm{cm}}^{-1}$ and ${\langle \varphi \rangle }_o$ was computed in each case. With water as reference we have also found ${\langle \varphi \rangle }_r$ (at ${\mu }_a$ and ${\mu }_s$ at $0.025{\mathrm{cm}}^{-1}$ and $14.41{\mathrm{cm}}^{-1}$, respectively) and then their ratio ${\langle \mathrm{\Delta }\varphi \rangle }_{{\mu }_a}$. For the ${\mu }_a$ and ${\mu }_s$ values used in the experiment we have also computed $\langle \mathrm{\Delta }\varphi \rangle$’s using MC simulation. The results are given in Fig. 5. These more or less follow the variation predicted by Eq. (2). Therefore from the experimentally measured ${\langle \mathrm{\Delta }\varphi \rangle }_{{\mu }_a}$ one can simply deduce ${\mu }_a^{\text{in}}$ if ${\mu }_a^b$ is known. In the second set we fix ${\mu }_a$ at $3.2\times {10}^{-3}{\mathrm{cm}}^{-1}$ and varied ${\mu }_s$ from 2.85 to $8.24{\mathrm{cm}}^{-1}$. The measured and computed ${\langle \mathrm{\Delta }\varphi \rangle }_{{\mu }_s}$’s are plotted in Fig. 5 which also verify the variation predicted by Eq. (4) leading to a simple read-out of ${\mu }_s^{\text{in}}$ from the measured ${\langle \mathrm{\Delta }\varphi \rangle }_{{\mu }_s}$.

Fig. 5

The plot of $\text{phase}={\langle \mathrm{\Delta }\varphi \rangle }_o/{\langle \mathrm{\Delta }\varphi \rangle }_r$against (a) absorption coefficient ${\mu }_a$ and (b) scattering coefficient ${\mu }_s$. In these plots ultrasound driving voltage is used as experimental parameter with $40V_{pp}$ and $60V_{pp}$. The operating ultrasound modulation beat frequency is 750 Hz. The experimental error is 0.10% to 13.82%. Liquid phantoms are used as the tissue-samples.

3.

Conclusion

In conclusion we have demonstrated the possibility of recovering the average ${\mu }_a$, ${\mu }_s$, and $E$ corresponding to the insonified ROI from the UAOT measurements of $M$ and $\langle \varphi \rangle$. The hitherto neglected $\varphi$ considering it as a zero-mean random variable is shown to have a non-zero mean when there is a large scattering anisotropy and acoustic frequency is small. This means that $\langle \varphi \rangle$ which carries useful information can be employed for the recovery of optical properties.

The other two advantages of the present study are: 1. Since we have not used a diffusion model in our inversion, the method should work well when the ROI is in a low scattering region like water-filled cyst. 2. Since we concentrate on the ROI, the reconstruction from these are not effected by bulk movement of body (like coming through respiration).