Ultrasound modulated optical tomography: Young’s modulus of the insonified region from measurement of natural frequency of vibration

R. Sriram Chandran; Debasish Roy; Rajan Kanhirodan; Ram Mohan Vasu; C. Usha Devi

doi:10.1364/OE.19.022837

1. Introduction

Ultrasound-modulated optical tomography (UMOT) is a non-invasive optical imaging modality with the potential of providing spectral variation in absorption and scattering coefficients (μ_a and μ_s respectively) in soft-tissue organs for medical diagnostic imaging [1, 2]. Focused ultrasound (US) beam helps provide a frequency shift to the photons passing through the insonified region-of-interest (ROI) in the body; and these photons can be used to image the optical property of the ROI. Since in a highly scattering object like tissue, the ultrasound beam can be focused more tightly than light beam, UMOT provides parameter reconstruction with an increased spatial resolution [3] than is possible with the usual stand alone near-infrared optical imaging modality, the diffuse optical tomography (DOT) [4].

Functional parameter estimation requires quantitative information on the optical properties, especially μ_a. However, to extract such information accurate mathematical models describing the interaction of US wave with tissue, and the propagation of light through the insonified ROI are required. In the past, considerable effort has been expended to arrive at and numerically study such interaction models [5–8]. However, to the best of our knowledge, there have only been a few attempts at quantitative recovery of optical parameters from UMOT data [9–12].

The route to a quantitative recovery has its moorings in the computation of force experienced by the object in the region insonified by the US transducer. This can be done, making use of the data on the electrical power input to the piezo transducer and propagating the pressure generated on its surface to the focal region [13]. As explained in [14], force is proportional to the average intensity in the incident pressure waves and the acoustic absorption in the region; and the intensity in turn is proportional to the rate of change of the average energy deposited in the region. The applied force is sinusoidal, being of the same frequency as the US transducer [15]. With the force and knowledge of elastic parameters of the tissue, a momentum-balance equation is set up and solved for the periodic oscillation experienced by the ROI. In addition to the oscillations, the acoustic force results in compression and rarefaction of the object resulting in a refractive index modulation of the ROI which is not considered in this work. Finally, photons are transported through the body and their interaction with the oscillating ROI is captured using a model similar to the one used in diffusing wave spectroscopy (DWS) to transport photons through a turbid medium. Homodyning of the light exiting the object results in a modulated speckle pattern whose depth of modulation (M) can be inferred from the intensity correlation (g₂(τ)) of the detected light. This model allows one to relate M to μ_a and μ_s, the refractive index (n) and the mean-squared displacement (< a² >) of the scattering centres of the ROI. Through < a² >, M is also related to the (visco) elastic properties of the ROI.

The above models are used to derive explicit expressions linking M to μ_a,μ_s, n and < a² > (in the ROI) opening up the possibility of recovery of these parameters from the set of measurements {M} over many ‘views’. A partial differential equation (PDE) describing the propagation of the amplitude correlation of light (G) through a turbid medium insonified by an US beam is derived in [7], which can also be employed for the estimation of these object properties from the boundary measurements of G. Alternatively, the above non-linear parameter estimation problem can be recast as a linear source recovery problem of reconstructing the ultrasound-induced sources of perturbation in G (which are a(r) and Δn(r)) from the boundary measurements of this perturbation (G^δ).

The aim of the present work is to demonstrate quantitative recovery of elastic property of a tissue-like object from variation of M.vs.frequency of US forcing, in particular the frequency at which M has a dominant peak. To facilitate measurement of M at different US forcing frequencies we employ (as detailed further in Section 2) a confocal two-region transducer, one region resonating at a frequency (f) and the other at (f + Δf). The mixing of US pressure amplitudes produces a forcing at the beat frequency Δf [15]. By varying the drive frequency of the transducer the forcing frequency can be varied. The essential strategy of the method here is to match (or tune) the forcing frequency at the ROI to the dominant natural frequency of the vibrating ROI (VROI), and the matching of the frequencies is helped by the fact that the peak in M vs frequency curve is indicative of resonance in the VROI. This natural frequency is a function of the Young’s modulus (E) of the material of the ROI, presently assumed to be linear elastic (Hookean) homogeneous and isotropic, which opens the way for a quantitative assessment of E from a measured resonant frequency. We note here that the natural modes of vibration depend not only on elastic properties of the VROI, but also on its shape/geometry and the boundary conditions. Therefore, as further elaborated in Section 3, for the recovery of E one needs to correctly obtain the shape of the VROI followed by an assignment of the correct local boundary conditions.

A summary of the work presented in the rest of the paper is as follows. Section 2 contains computation of the pressure and force in the focal volume of a two-region PZT transducer, and its experimental verification by measuring this pressure distribution. The solution of the momentum balance equation with the forcing from the US transducer for arriving at the dynamic response of the ROI is described in Section 3. Prior to this, we arrive at the shape of the ROI which undergoes non-trivial vibration. Transport of light through the insonified region and computation of the intensity autocorrelation of light and its modulation are the subjects discussed in Section 4. Section 5 provides details of the experimental set-up for measuring M for various Δf ’s. Apart from this, in this section we also describe a method to arrive at the unknown E from the measured dominant natural frequency using the method of bisection. Section 6 gives our concluding remarks.

2. Perturbation from ultrasound: modelling and verification

2.1. Ultrasound-induced force distribution in the focal volume

The confocal composite US transducer we use is driven by an input voltage just large enough to cause a modulation in the focal region and whose effect in the intervening region can be conveniently neglected. Focusing is normally done either using a concave transducer on which PZT is deposited or a plane one in conjunction with a US focusing element or an array with provision for controlled sequential firing for generating a spherical wave. In the present work the two regions of the confocal US transducer are driven by voltages oscillating at frequencies 1 MHz and 1 MHz + Δf where Δf can be varied from as low as 10 Hz to 1 kHz. The composite transducer helps us send MHz waves which have reasonable depth of penetration in the object and thus apply a sinusoidal force at the ROI oscillating at Δf Hz, the beat frequency.

Beating takes place at the ROI because the force is proportional to the acoustic intensity which is the square of the acoustic pressure amplitude. We use p₁(f) and p₂(f + Δf) to denote the sinusoidal pressure amplitudes at the ROI from the two regions in the transducer.

p_{1} (f) = p_{10} \cos (2 π ft + ϕ_{1})

p_{2} (f + Δ f) = p_{20} \cos (2 π (f + Δ f) t + ϕ_{2})

The average intensity, I(X), of the combined US beam is given by

I (X) = \frac{1}{s} \frac{d < \tilde{ξ} (X) >}{dt}

where ξ(X) = < ξ̃(X) > is the time average of the energy deposited at X in the ROI, and s the cross-sectional area we are considering in the intersection volume of the two pressure waves. This energy ξ(X) can be shown to be equal to [14]

ξ (X) = P_{0} + P_{1} \cos (2 π Δ ft + ϕ_{2} - ϕ_{1})

where P₀ and P₁ are constants dependent on

p_{10}^{2}

and

p_{20}^{2}

and the integration time constant T. The radiation force experienced by the ROI in an area around X is given by

F (X) = \frac{2 α I (X)}{c}

where α is the sound absorption coefficient and c the velocity of sound in the object. This force can be shown to be equal to [14]

F (X) = F_{0} \sin (2 π Δ f t)

where

F_{0} = \frac{4 α P_{0}}{T c s}

The pressure amplitude, p₁ or p₂, at the focal region can be computed by propagating the pressure at the transducer surface. The pressure at the transducer surface can be calculated from the electrical parameters of the transducer, such as the equivalent resistance, capacitance and inductance (R, C and L respectively) of the PZT. The acoustic pressure p₀ on the transducer surface is given by [16]

p_{0}^{2} = \frac{2 ρ c}{A} k_{eff}^{2} \frac{V^{2}}{R}

where ρ is the material density and c the acoustic wave velocity, A the surface area of the transducer, V the applied voltage and

k_{eff}^{2}

the conversion coupling coefficient (electrical to acoustic and specified by the manufacturers to be 0.62 in our case) of the transducer.

2.2. Ultrasound-induced pressure distribution in the focal volume

In order to transport the pressure distribution p₀ on the surface of the transducer to its focal volume, where the pressure can be very large, one needs a model that accounts for material nonlinearity apart from absorption and diffraction in the medium. Here we make use of such a model known as the Westervelt equation [13].

\nabla^{2} p - \frac{1}{c^{2}} \frac{\partial^{2} p}{\partial t^{2}} + \frac{δ}{c^{4}} \frac{\partial^{3} p}{\partial t^{2}} + \frac{β}{ρ c^{4}} \frac{\partial^{2} p^{2}}{\partial t^{2}} = 0

The first two terms of the above equation represent diffraction and linear propagation and the third term thermoviscous losses. The nonlinear effects are governed by the fourth term with a constant β which depends on the average transducer pressure and the density of the medium. The acoustic absorption

α = \frac{δ ω^{2}}{2 c^{3}}

where δ is the sound diffusivity.

The Westervelt equation is solved in [13] using a series expansion in spheroidal coordinates and we have adopted the same technique here. Here the values of α and β used are α = 1.5 dBcm⁻¹ and β = 6.2 to match the values quoted for breast tissue in [17]. We have computed the focal region pressure distribution for (i) transducer operating with just one element and (ii) the composite transducer with both the regions energized. We did this computation for US waves going into water as well as the (poly) Vinyl Alcohol (PVA) phantom used in the experiments. The axial (i.e., along the axis of the transducer) and lateral cross sections of the pressure distributions computed are shown in Figs. 1 and 2. Figure 3 shows the mesh and contour plots of the computed pressure in the x–z plane which contains the transducer axis taken to be along the z-direction.

Fig. 1 Computed pressure distribution along the US transducer axis (i.e., (0,0,z)) in the ROI. Here d is the focal length of the US transducer. The origin is assumed to be at the centre of the ROI.

Download Full Size | PDF

Fig. 2 Computed pressure distribution normal to the US transducer axis computed along (x,0,0). This distribution is observed to have radial symmetry. Here a is the aperture radius of the US transducer. The origin is assumed to be at the centre of the ROI.

Download Full Size | PDF

Fig. 3 Mesh and contour plots of the pressure distribution in the (x,0,z) plane. Here d is the focal length and a is the aperture radius of the US transducer. The origin is assumed to be at the centre of the ROI.

Download Full Size | PDF

2.3. Experimental verification of the shape of the focal volume

The US transducer, immersed in water is powered by two signal generators working at 1 MHz and 1.01 MHz and the force in the common focal volume in water is measured by scanning the region with a PVDF (Polyvinylidene fluoride) transducer. For a schematic diagram of the set-up and other details, see Section 5. The acceptance area of the transducer is ∼ 7mm² which is masked to a smaller area of ∼ 0.38mm². The transducer is mounted on an x-y-z stage and the voltages generated are measured at a number of points in the focal region (at many x-y cross-sections). Assuming that the measured voltages are proportional to the force experienced by points in the ROI, cross-sectional plots of the force experienced by the ROI are generated. The measured voltage distributions, giving axial and radial cross-sections, are shown in Figs. 4 and 5 respectively. Their variations are similar to (and compare well with) the computed pressure variation shown in Figs. 1 and 2. Figure 6 is the contour plot of the experimentally measured voltage in the x–z plane.

Fig. 4 Experimentally measured voltage distribution along the US transducer axis (i.e., (0,0,z)) which is an experimental measure of the computed pressure shown in Fig. 1. d is the focal length of the US transducer. The origin is assumed to be at the centre of the ROI.

Download Full Size | PDF

Fig. 5 Experimentally measured voltage distribution normal to the US transducer axis along both (x,0,0) and (0,y,0) which corresponds to the computed pressure distribution shown in Fig. 2. Here a is the aperture radius of the US transducer. The origin is assumed to be at the centre of the ROI.

Download Full Size | PDF

Fig. 6 Contour plot of the experimentally measured voltage (normalized) in the x–z plane. Here d is the focal length and a is the aperture radius of the US transducer. The origin is assumed to be at the centre of the ROI.

Download Full Size | PDF

3. Response of the periodically excited focal volume via solution of the momentum-balance equation

We characterize the vibration induced in the US focal volume by setting-up and solving the eigenvalue problem associated with the unforced linear momentum balance equation in the small-deformation (geometrically linear) regime. The force delivered by the US to the ROI is computed in Section 2.1 and the shape of the ROI is experimentally ascertained in Section 2.3. (A theoretical-numerical verification of the shape of the US focal region can be had from the solution of the forced momentum balance equation as described further in this section). Apart from the force experienced by the ROI and its shape, we need the material properties such as ρ, E and ν (the mass density, Young’s modulus and Poisson’s ratio respectively) to be able to set-up and solve the linear momentum balance equation, which is

ρ \ddot{u} = \nabla \cdot σ + F (x) \cos (2 π Δ f t) I_{Ω}

where σ = σ(x) is the Cauchy stress tensor and I_Ω is the indicator function for the VROI which suffers oscillations. We also assume that ρ is time-invariant. σ in Eq. (7) is a function of the spatial (deformed) coordinates that are themselves unknown. Here x is a spatial coordinate of a material point X after deformation. For the small deformation regime we are interested in, we take F(X) = F(x). A more convenient representation of this equation, referred to the reference (material) coordinates, is given by [18]

ρ \ddot{U} = \nabla \cdot P + F (X) \cos (2 π Δ ft) I_{Ω}

where P is the first Piola-Kirchhoff stress tensor, related to σ through the following identity:

ℱ = I + \frac{\partial u}{\partial X}

and P= det(ℱ) σ · ℱ^−T, where ℱ is the deformation gradient, and U = U(X) is the displacement vector measured in the reference coordinates.

We have employed ANSYS^® (a commercially available finite element programme) to carry out the eigenanalysis of the vibrating ROI, obtaining as a result the subset of its natural frequencies corresponding to the dominant vibration (due to shear deformation along the transducer axis). Towards this, the initially selected ∂Ω (the closed boundary of Ω) obtained experimentally in Section 2.3 is refined during the eigenanalysis. Theoretically verified/corrected ∂Ω is obtained by solving the momentum-balance equation over one full cycle of forcing, followed by identifying the boundary nodes (of the VROI) that separate the vibrating region from its non-vibrating complement. (A good initial guess of ∂Ω can also be had from the force distribution in the ROI. The boundary nodes may approximately be considered as those where F(X) drops to $\sim {(\frac{1}{100})}^{th}$ of its maximum value). From the experimental measurement of pressure in the VROI, we notice that of the two confocal regions of the transducer, one (the outer one with f = 1 MHz) gives rise to pressure distributions slightly different with an asymmetry compared to the other one (Figs. 2 & 5). This is owing to a defect in the spherical mould on which PZT is coated.

The photons, on interaction with the US focal volume, pick-up a phase modulation owing to the oscillations of the scattering centres. As shown in [19] the detected g₂(τ) contains an oscillation of the same frequency as that of the scattering centres. The oscillations of the VROI under US forcing are characterized by a set of eigenvalues (natural frequencies) and eigenvectors (mode shapes) which contain both compression and (dilatational) shear modes. The US vibration at the low beat frequency, however excites mainly shear modes in the VROI (along with the dilatational modes excited by the high frequency wave at 2f). The eigenanalysis helps us identify these shear modes corresponding decidedly to the lower end of the eigenspectrum. The shear waves attenuate rather quickly in the surrounding region of the ROI as shown in [12], but are present in Ω. The shear wave attenuation is frequency dependent being higher as the frequency increases. However, as verified through ANSYS-based numerical assessment of the shear deformed region, the attenuation even at these low frequencies is fast enough for the shear displacement to be close to zero beyond the near-hyperboloidal VROI.

The shear modes are typically of the order of a few Hz to ∼ 300 Hz for the nearly-incompressible (ν = 0.49) jelly-like material used in our experiments and simulations. The forced vibration at the k^th nodal point of the VROI can be expressed as

x_{k} (t) = \sum_{l = 1}^{n} c_{l}^{k} Ψ_{l} (t) + γ_{k} (t)

where Ψ_l(t), a combination of sinω_lt and cosω_lt, is the l^th modal basis function,

c_{l}^{k}

are the coefficients in the expansion, γ_k(t) is the particular solution under US forcing and n is the number of degrees of freedom in the discretized VROI. If the US forcing is at one of the natural frequencies, say ω_p, then γ_k(t) will be of the form

{\tilde{c}}_{p}^{k} Ψ_{p} (t)

. In such case, owing to resonance, the last term will predominate all other terms in Eq. (9) with

| {\tilde{c}}_{p}^{k} | ≫ | c_{l}^{k} | \forall l \neq p

. Thus the ROI will oscillate with the same frequency as that of the individual scattering centres. As long as the beat frequency remains sufficiently low (as in the present case), the resonance observed will be owing to shear deformation.

4. Propagation of light through a turbid medium insonified by the US beam

Light propagating through a turbid medium picks up a phase fluctuation even in the absence of the US forcing owing to the Brownian motion of the scattering centres induced by thermal excitation. This fluctuation manifests itself as a decay in the measured intensity autocorrelation (g₂(τ)). To such an object a focused US beam brings in an extra forcing, deterministic and localised. Extensive work has gone into the modelling of the effect of US forcing on the object. As pointed out in Section 1, scattering centres in the insonified region undergo periodic oscillations. The region itself, owing to the compression and rarefaction induced by the pressure wave,presents to the light beam a refractive index modulation Δn. This modulation, owing to the near-incompressibility of the tissue material, should occur with very low amplitudes with very high frequencies. A typical optical path length $\sum_{i} n_{i} l_{i}$ (here {n_i} and {l_i} are refractive index and path length corresponding to the i^th scattering event) gets perturbed to $\sum_{i} (n_{i} + Δ n_{i}) (l_{i} + u_{i})$ (here u_i is the amplitude of oscillation of the i^th scattering centre in the direction of l_i) which means that the US-induced path length fluctuation, to a first-order approximation, is Σ(n_iu_i + l_iΔn_i), giving the US-induced phase perturbation as Δϕ = Δϕ_u + Δϕ_n = k₀Σn_iu_i + k₀Σl_iΔn_i where Δϕ_u and Δϕ_n are the perturbations owing to oscillation and Δn respectively.

Expressions are available for Δϕ_u and Δϕ_n when the medium has either scattering isotropy or anisotropy, which contain optical, acoustic and elastic properties of the medium as parameters [5]. These are used in the equation for the US-introduced contribution towards the temporal autocorrelation function (G₁(τ)) of the electric field of light. In the weak scattering limit, assuming that there is no correlation between light arriving through different paths, photons travelling along the same path of length, s, have amplitude correlation given by $G_{1, s} (τ) = < E_{s} (t) E_{s}^{*} (t + τ) >$ . G_1,s(τ) is given by [5]

G_{1, s} (τ) = \exp [- F (τ) / 2]

where

F (τ) = < {[\sum_{j = 1}^{N} Δ ϕ_{u, j} (t, τ)]}^{2} >, neglecting Δ ϕ_{n}

Here (N – 1) is the number of scattering events along the path of length s. To arrive at the overall amplitude autocorrelation G(τ), G_1,s(τ) is ensemble averaged over all photon paths s: i.e.,

G (τ) = \int_{0}^{\infty} p (s) G_{1, s} (τ) d s

.

Here p(s) is the probability density function for s. The normalized amplitude autocorrelation g₁(τ) is $\frac{G (τ)}{G (0)}$ from which g₂(τ) is estimated as 1 + β|g₁(τ)|² where β is a constant dependent on the optics used in the data gathering.

In [5], analytical expressions are derived for G₁(τ); further p(s) evaluated through sending photons, using Monte Carlo (MC) simulations or by solving the diffusion equation for light propagation, is used to compute G(τ) and g₁(τ). The effect of periodic phase fluctuation on G(τ) is the appearance of a modulation of the same frequency as that of the US-induced acoustic vibration frequency, Δf, from which we derive our measurement which is the modulation depth M(Δf) obtained as the modulus of the Fourier transform of g₂(τ) evaluated at Δf. The MC simulation carries the advantage that in transporting photons one can compute both p(s) as well as Δϕ_u and their time averages. Δϕ_u is in fact the perturbation in the accumulated phase ϕ originating from the oscillations of scattering centres confined only to the insonified ROI.

Simulations proceed in the following lines: for a given (guessed) E of the phantom material, and the force at the focal region (evaluated using the pressure distribution in the focal region obtained from the Westervelt equation, in Eqs.(1)–(4) [13]), first the amplitude of vibration of scattering centres in the ROI is computed. Then 1 million photons are launched and tracked for those arriving at the detector. From the phase and weight of the photons G(τ) is calculated, from which g₁(τ), g₂(τ) and M are evaluated [19, 20]. The above are repeated for Δf varying from 40Hz to 130Hz. A typical M vs Δf plot is shown in Fig. 7 for which E= 11 kPa. It is observed from these simulations that M goes through a prominent peak (for a particular Δf, say Δf_r) and this Δf_r depends on E of the material. This indicates that the vibration of the ROI induced by the (US) forcing at the beat frequency Δf goes through a resonant peak at Δf_r, which results in a peak in phase modulation of photons traversed through the ROI, observed as a peak in the modulation depth in g₂(τ). This evidence of a resonant peak and the inferring of the frequency from the M vs Δf plot are also verified through experiments. As is shown in Section 5 this resonant frequency can be used to recover E of the VROI.

Fig. 7 The variation of the experimental modulation depth (M^e) with Δf which closely follows its simulated counterpart (M^c). The modulation depth was not measured at Δf = 50 and 60 Hz.

Download Full Size | PDF

5. Experiments

5.1. Fabrication and characterization of the tissue-mimicking phantom

The details of the fabrication and characterization of the PVA phantom are given in [20] and we only give a brief account here. PVA aqueous solutions with concentration of 20% by weight are prepared by heating PVA and demineralized water over a temperature bath at 90 to 93°C for 2 hours with continuous stirring. This solution is allowed to undergo repeated freeze-thaw cycles of freezing at −20°C for 12 hours, followed by thawing at room temperature for 12 hours [21]. E and μ_s in the PVA slab can be tailored by varying the number of freeze-thaw cycles. As mentioned in reference [20] E can be varied from 11 to 97kPa and correspondingly μ_s vary from 1.4 to 4.5 mm⁻¹. For further increasing μ_s we mix polystyrene beads of diameter 0.1 μm whilst forming the PVA slab. The scattering coefficient of the resulting phantom is determined using an indirect method based on the reverse Monte Carlo procedure [20] and E is measured and verified using a dynamic mechanical analyzer (DMA). The E for the three samples used are found to be 11 kPa, 45 kPa and 58kPa. We mixed 0.3ml of polystyrene bead solution in 100ml of stock solution and increased the μ_s to 8.14 mm⁻¹.

5.2. Measurement of modulation depth

The experimental set-up is shown in Fig. 8. Light from a He-Ne laser (HRR170, Thorlabs, power 17 mW) illuminates the object insonified by the US transducer, which is a two element confocal one with a common focal length of 50mm (manufactured by Roop Telsonic Ultrasonix Ltd, Mumbai, India).The driver for the transducer is an ulra-stable dual-channel function generator (Tektronix Model AFG3022B). The output from function generator is amplified by a power amplifier (Tomco, Model No:BT00100-AlphaD-CW). A continuously variable neutral density filter from Thorlabs (NDC-100C- 4M) is used to control the intensity of light entering the object. The light gains access to the object through a hole in the transducer. The photons exiting the object are picked up by a photon-counting photo-multiplier tube (PC-PMT, Hamamatsu H7360-03) and current generated after signal conditioning is given to a correlator (DAC, Flex 021d, from www.correlator.com,Bridgewater NJ 08807) which can give a minimum delay time of up to 12.5 ns. The correlator output is given to the computer for further processing. The US transducer and the object are immersed in a water bath for acoustic impedance matching. The laser light is aligned to be along the US transducer axis. The composite US transducer with two elements oscillating at 1 MHz and (1 + Δf) MHz is driven by an ultra-stable dual-channel function generator with suitable power amplification. To achieve a frequency scan of the ROI starting at a Δf as low as tens of Hz we drive one of them slightly away from resonance. A slight drop in output because of this deviation is neglected. We note that the ability to generate a stable low frequency acoustic wave critically depends on the frequency stability of the driving function generator. Our ability to select a frequency to match the desired computed resonant modes (see Section 3) is hampered because of the limitation of possible frequencies that can be selected in the function generator. The photons are picked up by a single-mode fibre and delivered to the PC-PMT. For maximizing the signal-to-noise ratio the fibre is aligned carefully to capture a single speckle from the boundary photon fluctuation.

Fig. 8 The Experimental Setup: Light from laser L, illuminates the insonified volume through a hole in the ultrasound transducer (UST). The exiting light is detected by the detector (PC-PMT) and is given to the correlator DAC and then to a computer C. The UST is driven by ultra- stable dual-channel function generator (DCFG) after power amplification

Download Full Size | PDF

In the experiment the objects which are slabs made of PVA of dimensions 6cm×6cm×20cm, designed to have a μ_s equal to that of healthy breast tissue (8 mm⁻¹) and μ_a kept at a nominal low value of 0.00025 mm⁻¹ and E = 11 kPa, 45 kPa and 58 kPa are used. Light is gathered over 1200 s and the intensity autocorrelation g₂(τ) is measured. For a single g₂(τ), 2000 samples are collected and averaged. g₂(τ) gathered for slabs with E = 11 kPa to 58 kPa are analyzed to arrive at M^e (experimentally measured modulation depth) which is measured from the area around the peak at Δf in the Fourier spectrum of g₂(τ).

The variation of M vs Δf is shown in Fig. 7 for a typical slab with E = 11 kPa. The computed M (M^c) is also shown in the same figure. M^c is seen to closely follow M^e at the peak, an evidence of the resonance in vibration of the VROI, observed at ∼ 70 Hz.

Figures 9–11 show variation of M^e with Δf for objects with E = 11, 45 and 58 kPa respectively. For comparison, the locations of the resonant modes of the VROI as computed using ANSYS^® (in Section 3) are also shown. These figures prove that the peak resonant modes can be measured using the peak observed in the M^e.

Fig. 9 Variation of the modulation depth (M^e) with Δf compared with the similar variation of the computed resonant amplitude (obtained through ANSYS^®) for the 11 kPa phantom. Both M^e and the computed amplitudes are normalized using their maximum value observed at Δf_r which is 70 Hz for this phantom

Download Full Size | PDF

Fig. 10 Same as those given in Fig. 9 repeated for 45 kPa phantom

Download Full Size | PDF

Fig. 11 Same as those given in Fig. 9 repeated for 58 kPa phantom

Download Full Size | PDF

We make the following observations in regard to the results of Fig. 7 and Figs. 9–11. In Fig. 7 there is agreement between M^e and M^c only around the resonant frequency. A possible reason for the deviation at other frequencies is the error accrued in the computation of the acoustic radiation force consequent to the uncertainty in the values used to account for coupling loss between transducer and object, losses at the interface between water and phantom and also the variable acoustic absorption in the object. An incorrect magnitude of the forcing used does not however shift the resonant frequency but only its strength. Also, the other errors we may have incurred may have a broad spectrum and therefore hardly contribute to the negligibly small bandwidth around the resonant frequency wherein most of the structural energy is concentrated.

In Figs. 9–11. the location of the computed resonant modes are not the same as where M^e is found experimentally. This mismatch is owing to the spread of the resonant modes as dictated by the mechanical properties and shape of the vibrating ROI (as computed by ANSYS) and the fact that the set of M^e was evaluated at frequencies possible to be selected in the function generator around 1.01 MHz. Limitations on the stability of the function generator used precludes generation of more experimental sample points in these frequency ranges.

Since our recovery of E is from the measured resonant frequency data-model mismatch elsewhere has not affected the accuracy of the recovered E. However in an algorithm that tries to reconcile a ‘data-model misfit’, the difference between computed and experimental ‘measurements’ owing to the incorrect input of certain quantities to the algorithm, say the force in the insonified area, can pose serious problems. One possible way to overcome this difficulty is to consider force also as a parameter to be identified in the inversion algorithm.

From the measured resonance modes in oscillations (Δf_r = 70 Hz,140 Hz and 250 Hz for the three slabs) assuming that the E’s are not known, we use the method of bisection to arrive at the unknown E. Thus, for a guess of two E values, say E₁ and E₂ such that E₁ < E₀ < E₂ where E₀ is the unknown Young’s modulus to be found, the resonant modes, Δf₁ and Δf₂, are computed. If the interval between mean $Δ f_{m} = (\frac{Δ f_{1} + Δ f_{2}}{2})$ and Δf_i (i = j with j = 1 or 2) contains Δf_r, the guess of E is modified to $E_{m} = (\frac{E_{1} + E_{2}}{2})$ and E_j thus reducing the length of the original search interval ΔE = E₂ – E₁. This is continued until the computed resonant frequencies match within a specified tolerance. The E values thus obtained from the measured f_r’s are 11.4, 44.8 and 58.7 kPa for the three objects investigated.

6. Conclusions

We have demonstrated, experimentally and with theoretical backing, quantitative recovery of Young’s modulus from UMOT data, which is the modulation depth introduced in the light intensity autocorrelation by the US induced dynamics. The method depends on measurement of a frequency and therefore is unaffected by phase and amplitude variations in the light source. In other words we have a measurement which is unaffected by the fluctuations in detected intensity. Since the detection of modulation depth in UMOT is a challenging task owing to speckle noise, the enhanced SNR afforded in measurement of M because of resonance will make the usual UMOT based recovery of optical absorption coefficient more accurate, which can be cited as an added benefit resulting from the present work. The fact that the variation of Δf provides different sets of modulation depth measurements, can be further utilized to do tomographic inversion for μ_a, n and E for the ROI, which is left for investigation in the future.

References and links

1. C. Kim and L. V. Wang, “Multi-optical-wavelength ultrasound-modulated optical tomography: a phantom study,” Opt. Lett. 32, 2285–2287 (2007). [CrossRef] [PubMed]

2. L. V. Wang, “Ultrasound-mediated biophotonic imaging: A review of acousto-optical tomography and photo-acoustic tomography,” Dis. Markers 19, 123–138 (2004). [PubMed]

3. S. Leveque, A. C. Boccara, M. Lebec, and H. S. Jalmes, “Ultrasonic tagging of photon paths in scattering media: parallel speckle modulation processing,” Opt. Lett. 24, 181–183 (1999). [CrossRef]

4. A P Gibson, J C Hebden, and S R Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50, R1–R43 (2005). [CrossRef] [PubMed]

5. S. Sakadzic and L. V. Wang, “Ultrasonic modulation of multiply scattered coherent light: An analytical model for aniostropically scattering media,” Phys. Rev. E 66, 026603 (2002). [CrossRef]

6. M. Kempe, M. Larionov, D. Zaslarsky, and A. Z. Genack, “Acoustooptic tomography with multiple scattered light,” J. Opt. Soc. Am. A 14, 1151–1158 (1997). [CrossRef]

7. S. Sakadzic and L. V. Wang, “Correlation transfer and diffusion of ultrasound-modulated multiply scattered light,” Phys. Rev. Lett. 96, 163902 (2006). [CrossRef] [PubMed]

8. W. Leutz and G. Maret, “Ultrasonic modulation of multiply scattered light,” Physica B 204, 14–19 (1995). [CrossRef]

9. C. Kim, R. J. Zemp, and L. V. Wang, “Intense acoustic bursts as a signal-enhancement mechanism in ultrasound-modulated optical tomography,” Opt. Lett. 31, 2423–2425 (2006). [CrossRef] [PubMed]

10. R. J. Zemp, C. Kim, and L.V. Wang, “Ultrasound-modulated optical tomography with intense acoustic bursts,” Appl. Opt. 46, 1615–1623 (2007). [CrossRef] [PubMed]

11. X. Xu, H. Zhang, P. Hemmer, D. Qing, C. Kim, and L. V. Wang, “Photorefractive detection of tissue optical and mechanical properties by ultrasound modulated optical tomography” Opt. Lett. 32, 656–658 (2007). [CrossRef] [PubMed]

12. E. Bossy, A. Funke, K. Daoudi, A. Boccara, M. Tanter, and M. Fink, “Transient optoelastography in optically diffusive media,” Appl. Phys. Lett. 90174111 (2007). [CrossRef]

13. T. Kamakura, T. Ishiwata, and K. Matsuda, “Model equation for strongly focused finite-amplitude sound beams,” J. Acoust. Soc. Am. 107, 3035–3046 (2000). [CrossRef] [PubMed]

14. E. Konofagou, J. Thierman, and K. Hynynen, “A focused ultrasound method for simultaneous diagnostic and therapeutic applicationsa simulation study,” Phys. Med. Biol. 46, 2967–2984 (2001). [CrossRef] [PubMed]

15. M. Fatemi, A. Manduca, and J. Greenleaf, “Imaging Elastic Properties of Biological Tissues by Low-Frequency Harmonic Vibration,” Proc. IEEE 91, 1503–1517 (2003). [CrossRef]

16. J. Huang, J.A. Nissen, and Erik Bodegom, “Diffraction of light by a focused ultrasonic wave,” J. Appl. Phys. 71(1), 70–75 (1992). [CrossRef]

17. F. A. Duck, “Nonlinear acoustics in diagnostic ultrasound,” Ultrasound Med. Biol. 28(1), 1–18 (2002). [CrossRef] [PubMed]

18. J.E. Marsden and T.J.R. Hughes, Mathematical Foundations of Elasticity (Dover Publications, Inc., New York, 1993)

19. S. Sakadzic and L. V. Wang, “High-resolution ultrasound-modulated optical tomography in biological tissues,” Opt. Lett. 29, 2770–2772 (2004). [CrossRef] [PubMed]

20. C. Usha Devi, R. M. Vasu, and A. K. Sood, “Design, fabrication, and characterization of a tissue-equivalent phantom for optical elastography,” J. Biomed. Opt. 10, 044020 1–10 (2005).

21. A. Kharine, S. Manohar, R. Seeton, R. G. M. Kolkman, R. A. Bolt, W. Steenbergen, and F. F. M. de Mul, “Poly Vinyl Alcohol gels for use as tissue phantoms in photoacoustic mammography,” Phys. Med. Biol. 48, 357–370 (2003). [CrossRef] [PubMed]

Ultrasound modulated optical tomography: Young’s modulus of the insonified region from measurement of natural frequency of vibration

Abstract

1. Introduction

2. Perturbation from ultrasound: modelling and verification

2.1. Ultrasound-induced force distribution in the focal volume

2.2. Ultrasound-induced pressure distribution in the focal volume

2.3. Experimental verification of the shape of the focal volume

3. Response of the periodically excited focal volume via solution of the momentum-balance equation

4. Propagation of light through a turbid medium insonified by the US beam

5. Experiments

5.1. Fabrication and characterization of the tissue-mimicking phantom

5.2. Measurement of modulation depth

6. Conclusions

References and links

Cited By

Figures (11)

Equations (11)

Optics Express