1. INTRODUCTION

For noninvasive elasticity imaging of soft-tissue organs, local application of a known force through focused ultrasound beams has been suggested and implemented in a number of papers [1,2]; this process has come to be known as remote palpation. Local sinusoidal forcing, the frequency of which can be varied, is applied to a selected region (also called the region-of-interest, ROI) using focusing ultrasound transducers, which results in oscillatory strain in the vicinity of the ROI [3,4] that is measured. Transient ultrasound forcing also creates shear wave in the neighborhood of the ROI, which propagates through the object with its amplitude tracked by a Doppler ultrasound scanner [5]. An ultrasound pulse-echo correlator has also been employed to make quantitative measurements in the context of remote palpation imaging [6]. A major turning point in this line of research came with the discovery of the vibro-acoustic (VA) wave, a secondary acoustic emission caused by the local ultrasound forcing [3,4]. The VA wave, detected outside the object, is shown to carry information on the mechanical stiffness of the object in so far as its amplitude and phase can be used to create a qualitative image of the stiffness [7–10]. Attempts were also made to obtain quantitative images (of, say, shear modulus distribution) from VA amplitude and phase [7]. Since an accurate recovery of the elastic property from VA data would involve accounting for the effect of the inhomogeneous visco-elastic medium on the VA wave, quantitative inversion of VA data involves solving a difficult inverse problem [11].

In our recent work, a simplified inversion was proposed by measuring, instead of amplitude and phase, the spectral peaks of the VA wave [12,13]. In this, we have taken a clue from resonant ultrasound spectroscopy (RUS), proposed for ascertaining mechanical properties of objects such as granite samples from the measured ultrasound (or, VA) spectral peaks. So long as one is able to gather the vibrational modes of the ROI from the VA spectrum, these modal frequencies can be inverted for visco-elastic properties of the ROI. In Ref. [12], it has been shown, both theoretically and through experiments, that most of the spectral peaks of the VA signal are indeed the resonant modes of the ROI. In fact, it has also been shown that the VA pressure waves, though primarily influenced by the bulk modulus of the material, also carry information on the local shear modulus through their coupling with the shear waves generated in the ROI. Consequently, a method of elasticity recovery from frequency data, insensitive to noise, is demonstrated in Ref. [12]. In addition, a frequency-based ultrasound remote palpation technique is demonstrated in Ref. [14] for estimating the elastic modulus of arterial vessels.

Measurement of VA spectral peaks in Ref. [12] is done through a brute-force approach, by scanning the ultrasound drive frequency and looking for frequencies at which the VA amplitude peaks. A limitation of this method is owing to its inability to detect higher-order modes, which are characteristically weak. The same limitation is observed in Ref. [15] in the application of RUS for characterizing rock samples: a dissipating medium will quickly smother the weak high-frequency modes present in the VA signal as it propagates through the object being studied. This paper demonstrates an optical probe for examining the ROI to map these VA spectral peaks: a coherent light beam sending photons to meander through the object, including the ultrasound focal volume, with the facility to measure the decay in photon correlation orchestrated by its scattering from moving particles in the ROI. This is the well-known diffusing-wave spectroscopy (DWS), used to study temperature-induced Brownian motion in a multiple-scattering object. A small modification is effected by the presence of ultrasound-induced periodic motion of the scattering centers in the ROI, a tag with which we can identify photons that have intercepted the ROI and thus recover Brownian dynamics pertaining to the ROI [16].

In DWS, the temporal variation in the mean-squared displacement (MSD) of Brownian particles is estimated. In Ref. [17], an important observation is made highlighting the difference in the MSD plots between cases when ultrasound-induced dynamics superimposing the Brownian motion is present or withdrawn. The plateau of the MSD-versus-time plot contains noise-like fluctuations when ultrasound is present. In Ref. [17], a model for simulating the behavior of MSD is put forth in which the dynamics of the ROI (considered a continuum) is modeled through that of a representative particle, the system particle. The system particle is modeled to obey a generalized Langevin equation (GLE). The major contribution of Ref. [17] is the means to arrive at a multiplicative noise term, which when used in the GLE captures the fluctuations in the MSD plot. The multiplicative noise represents the overall effect on the system particle of a distribution of “bath particles” in the ROI, thus capturing the nonlocal interparticle interactions. If material length scales are comparable to deformation and forcing space–time length scales, terms to account for nonlocal interactions must be built in to the model for it to accurately represent the system’s macro behavior. In Ref. [17], this is done by accounting for the micro-rotations the particles undergo (modeled using the standard micropolar theory [18]) in addition to translation. A salient feature of these micro-rotations is the inherent randomness, which arises owing to the uncertainty relation that exists between two strain operators, one containing micro-rotational information, and the other, translational information. In Ref. [19], it has been shown that the dynamics of a representative particle, in the context of ultrasound-assisted tomography, is indeed the dynamics of the ROI. In the present case, they have one Brownian component owing to the thermal fluctuations and one induced by ultrasound forcing, the latter containing both translational and rotational modes. In the absence of ultrasound, the weak thermal stress alone is insufficient to excite these additional modes, enabling the conventional GLE without the multiplicative noise to represent the dynamics and leading to MSD plots without fluctuations. But the fluctuation, owing its origin to the ultrasound exciting additional modes within the ROI, in fact captures the extra dynamics. In this work, we exploit this information present in the noise-like fluctuations in MSD plots. In fact, as seen in Section 4.B, the Fourier transform of the fluctuations gives us the frequencies of modes excited by the ultrasound forcing within the ROI. From the experimentally gathered modal frequencies, a complete recovery of the elasticity tensor for an orthotropic medium, as shown in Section 4, is possible. A complete recovery of the tensor components for an orthotropic (or, more generally, anisotropic) material, such as a malignant tumor, should be quite helpful to track the progress of the disease [20].

A brief summary of the rest of the paper is as follows. In Section 2, the basic theoretical foundation of the present work, which claims that the modal frequencies of the ROI are, in fact, contained in the fluctuations in the MSD plateau, is set forth. Here, the eigenvalue equation linking the mechanical properties of the material of the ROI to the modal frequencies is also described. This is our forward equation, connecting material properties to what is measured, viz., the natural frequencies. We give explicit expressions when the material is orthotropic, as in pork fat, or isotropic, as in agar [21]. Section 3 gives a description of the experiments, a setup for carrying out ultrasound-assisted DWS. Here, provision is made to gather photons that intercept the insonified region for computing intensity autocorrelation. Section 4 contains the details of the inversion algorithm used to recover orthotropic and isotropic elastic constants from the experimentally gathered modal frequency data. Section 5 includes a discussion of the major findings of this work. Our concluding remarks are in Section 6.

2. DYNAMICS OF THE ROI

We present in this section a brief description of the overall dynamics of the ROI, which has a temperature-induced Brownian component and an ultrasound-induced deterministic component. There are two major routes to model this dynamics: (1) treat the object as a continuum and use balance equations; (2) a random distribution of particles makes up the object with a single so-called “system particle” modeled to capture the dynamics of the whole object, with a distribution of “bath particles” used to capture the effect of nonlocal interactions of the system particle with the rest of the body. The basic work highlighting the physics of the problem along these lines has already been done, e.g., the continuum mechanics approach in Ref. [12] and the particle approach in Ref. [17]. Briefly drawing upon this basic theory, for the sake of completeness, we give an outline of these two approaches.

A. Dynamics Using the Continuum Model

In Ref. [12], generation of an acoustic field and an attendant shear wave, consequent to the sinusoidal forcing in the ROI by ultrasound, is discussed. It is shown that the low-frequency acoustic wave, also known as the VA wave, propagates carrying signatures of the shear wave (containing shear modulus information in the ROI), which is itself attenuated in the vicinity of the ROI. The cardinal equation governing the propagation of the VA wave with the pressure source confined to the ROI is Eq. (4) of Ref. [12], which is reproduced here:

B. Dynamics from the System Particle Approach

Here the continuum model is replaced by a collection of particles (harmonic oscillators), and the motion of the ROI is assigned to that of a single particle. The single particle, called the system particle, is generally assumed to have a single, predominant translational degree of freedom (DOF). After incorporating the effect of neighboring particles, a GLE is derived and used to evolve this DOF [22]. However, it is well known that the standard form of GLE may not be able to capture the macro system response, which is influenced by the material micro-structure. As discussed in Ref. [17], this inability is obvious when the forcing and response length scales are comparable to the micro length scales of the material. The present case of ultrasound forcing in a tissue-like object is an example where nonlocal interactions caused by micro-structural length-scale effects greatly influence the macro behavior of the object. In Ref. [17], a modified GLE is derived that has a multiplicative noise term with which the experimentally measured temporal evolution of the MSD of the system particle is more accurately reproduced. For this, recourse is taken to the micro-polar continuum theory, thereby incorporating nonlocal micro-structural effects in the model. The micro-polar theory additionally incorporates the rotational DOF for a particle. Using an uncertainty relation between the rotational and translational strain operators, a Gaussian, history-dependent noise term (referred to as the internal noise) is seen to modulate the stiffness term in the GLE [17]. With the internal noise, the solution to the modified GLE captures the growth, and more so, the fluctuations in the plateau of the MSD curve from an ultrasound-assisted local DWS experiment (see Fig. 3 in Ref. [17]). The absence of noise-like fluctuations when the ultrasound is withdrawn is because micro-rotations, from thermal movements alone, are negligibly weak. In addition, assuming the initial response to be temperature-driven fluctuations, a correlation structure for the additive noise term, similar to the second fluctuation dissipation theorem of Kubo, is also derived in Ref. [17].

The final form of the new GLE is (Eq. (20) of Ref. [17])

Equation (3) is solved for the evolution of the MSD, which can be verified against the data from DWS experiment. The MSD should be measured for the Brownian particles within the ROI, which is the intersection of the focal volumes of the two ultrasound transducers (see Section 3 dealing with experiments). Just as decay of intensity autocorrelation allows one to compute the MSD averaged over the entire object, the decay of the modulation depth overriding the intensity autocorrelation helps us arrive at the MSD for Brownian particles in the ROI [16]. This measurement is further discussed in Section 3.B.

Two major observations from the MSD-versus-time curve (Figs. 2 and 3 of Ref. [17]) are worth underlining, especially as one of these forms the basis for measurement of data for the inverse problem of elasticity recovery considered in this work. The first is that the transients in both the curves, one with and the other without ultrasound, are identical, confirming that the visco-elastic spectra recovered from both the curves match. This is, after all, anticipated, for the experiments in Ref. [17] were done on homogenous object samples. The second is, as stated earlier, the appearance of a noise-like fluctuation in the plateau when the ultrasound is on. As discussed earlier, this noise is accounted for by considering that the particles have rotational degrees of freedom in addition to translational. It is the ultrasound forcing that triggers the rotational DOFs, which in turn are coupled with the translational ones. It is thus no surprise that the prominent peaks in the Fourier transform of the noisy plateau indeed correspond to some of the natural frequencies of vibration of the ROI, with the weight decreasing with the mode number. As shown in the next subsection, we can also compute the resonant modes of vibration of the ROI given its mass density, elastic constants, and shape. Identification of spectral components of these fluctuations with the modal frequencies of the vibrating ROI is, we believe, the novelty of the work presented here. It is the enrichment of data through the DWS-orchestrated measurements that allowed characterization of a non-isotropic object. With experimental means of measuring these modal frequencies, and the availability of an easily implemented routine to compute them, it was easy to set up a scheme to invert the measured modal frequencies for the elastic parameters of the material in the ROI. [We note, in passing, that the only equation required in the forward model to compute the modal frequency data repeatedly is the eigenvalue equation (Eq. (4)]) to be discussed in the next section. However, to complete the forward model, i.e., to show that the DWS experiment indeed gives a noisy plateau in the MSD that contains the modal frequencies, we need the GLE to simulate the dynamics of the Brownian particles and their MSD at least once to verify.]

C. Computation of Natural Frequencies of Vibration of the ROI

This section describes how natural frequencies of a freely vibrating body can be computed from the shape, size, density, and elastic properties of the object, the details of which are already given in Ref. [12]. We repeat a summary of it here for the sake of completeness. The principle used in arriving at a computable equation is the well-known variational approximation for undamped, unforced mechanical systems: for a body with a stress-free boundary condition and with displacement vector $\mathbf{u}$ satisfying the elastic wave equation, the Lagrangian $\mathcal{L}={\int }_V(E_K-E_P)\mathrm{d}V$ of the system is stationary. Here, $E_K$ defines the kinetic energy equal to $(1/2)\rho {\omega }^2{u}_i{u}_i$, and $E_P$ the potential energy equal to $(1/2)C_{ijkl}{\partial }_j{u}_i{\partial }_l{u}_k$, where the Einstein summation convention over repeated indices is followed. Moreover, $u_i$, etc., are components of $\mathbf{u}$, which have harmonic time variations. $\rho$ and $C_{ijkl}$ are the material density and the elasticity tensor, respectively. With the expansion of $u_i$ in terms of a basis function set, $\{B_{\lambda }\}$, defined over $V$, $u_i={\tilde{u}}_{i,\lambda }{B}_{\lambda },\lambda =1,2.\dots N$, and using the principle that $\delta \mathcal{L}=0$, we reach the following symmetric eigenvalue problem [23]:

The introduction of a basis $B_{\lambda }$ [23] using monomials of the form $B_{\lambda }=x_1^{l(\lambda )}{x}_2^{m(\lambda )}{x}_3^{n(\lambda )}$ paved the way for the development of an easily computable forward model, which can be applied for any arbitrary-shaped (and sized) object. The $\{B_{\lambda }\}$-s, being a non-orthonormal set, makes $\mathbf{M}$ non-diagonal.

In this work, resonant modes were extracted starting from the eigenvalue problem of Eq. (4), making use of a standard commercial software ABAQUS. The inputs to ABAQUS modal analysis are (i) geometry of the tissue sample, (ii) material parameters such as density $\rho$ and elasticity tensor $C_{ijkl}$, (iii) appropriate boundary conditions. (We apply the Dirichlet condition of zero amplitude of vibration on the already identified boundary of the ROI. The procedure used to identify this Dirichlet boundary is given in Ref. [12] and is not repeated here.) The elasticity 4-tensor $\mathbf{C}=\{C_{ijkl}\}$ has 21 independent elastic constants for a general anisotropic material. Presently, we assume our material to be orthotropic. This assumption helps by reducing the number of independent constants to nine [24,25], and also allows us to visualize the constants along the three orthogonal axes.

The stress-strain relation for this case is given by $\mathbf{\sigma }=\mathbf{C}\mathbf{ϵ}$, where $\mathbf{ϵ}$ is the small strain tensor. Writing this in a matrix-vector form, one has

 

Fig. 1. Schematic representation of the object (slab) as an orthotropic material showing axes of symmetry.

Download Full Size | PDF

Our objective here is to employ the forward model of Eq. (4) to invert the measured natural frequencies of the ROI for its Young’s and shear moduli. In the next section, experiments to gather the resonant modes through ultrasound-assisted DWS are described. Both homogeneous phantoms and slabs cut from pork fat are made use of.

3. EXPERIMENTS

A. Preparation of the Phantom

The phantom used is made of agar gel molded into slabs of dimensions $65\mathrm{mm}\times 17\mathrm{mm}\times 18\mathrm{mm}$. The stiffness of the gel fabricated is increased by increasing the concentration of agar powder from 0.5% to 1.0% to 1.5% (by weight) in distilled water. The basic recipe of making the gel is given in Ref. [26]. Turbidity is created by adding Intralipid to the solution. Transport mean-free-path ($l^{*}$) of approximately 1.4 mm was obtained by adding Intralipid stock solution to make up 0.75% by volume [26]. From the same fabricated gel, discs of diameter 25 mm were also cut for independent measurement in a rheometer. The measured storage and loss moduli obtained from the rheometer (Discovery HR-3 rheometer from TA Instruments) using the frequency-sweep method, are shown in Figs. 2(a) and 2(b).

Freshly harvested pork fat was procured and preserved in a phosphate-buffered saline (PBS) solution of pH 7.4, maintained at $27\pm {0.5}^{0}C$. The PBS solution provided the fat samples an isotonic environment. The fat pieces were trimmed to make slabs of dimensions $65\mathrm{mm}\times 17\mathrm{mm}\times 18\mathrm{mm}$. The slab, once prepared, was immediately used in the experiment.

 

Fig. 2. (a) Variation of the storage moduli versus frequency of the agar slabs as measured by rheometer. (b) Same as (a) but showing loss moduli.

Download Full Size | PDF

B. Diffusing-Wave Spectroscopy with Ultrasound Tagging

The experimental setup is schematically shown in Fig. 3. This is basically a DWS setup for gathering intensity autocorrelation data [27]. It has a laser source illuminating the object under study and a detector, a photon-counting PMT, for gathering the multiple-scattered photons from the object. The current signal from the PMT is converted to a voltage signal and is sent through a signal conditioner and amplifier before being given to a digital autocorrelator [28], which time shifts and integrates the received signal. The output, the intensity autocorrelation, is stored in a computer for further processing.

 

Fig. 3. Schematics of the experimental setup. The ergodic medium-object slab composite object, immersed in a water bath (maintained at $27\pm 0.5$°C) is insonified by the intersecting focal volume (ROI) of two ultrasound transducers. The ROI is illuminated by a coherent beam from a He–Ne laser. The scattered light is detected by a single-mode fiber and given to a PMT. The output from the PMT after signal conditioning is input to the correlator and then to the personal computer. PA, power amplifier; US, ultrasound; ROI, region of interest; DCFG, dual channel function generator; CORR, correlator; COM, personal computer.

Download Full Size | PDF

The present setup has, in addition, provision for insonifying a region in the object with two ultrasound beams intersecting at their focal region. Dual-beam insonification allows forcing of the ROI at a low frequency, which is the difference in the drive frequencies of the two synchronously driven transducers, at which the acoustic signal is produced by “mixing” the two ultrasound pressure signals. With the average power from the two ultrasound transducers restricted to 15 mW, we made sure the temperature in the ROI did not increase appreciably from the background temperature. This low-frequency forcing excites the natural frequencies of vibration of the ROI, which are characteristically of low frequency. The transducers and object are immersed in a water bath for ensuring acoustic impedance matching. There is also provision for alignment of the transducers and object with the help of translation-cum-rotation stages.

In the experiments we take temporal measurements and time averages. One can replace the ensemble averages required for connecting decay in autocorrelation to the MSD of Brownian particles by time averages only if ergodicity can be assumed for the stochastic process governing the measured quantities [29,30]. Variation in photon number reaching the detector from either the agar gel or pork fat is not the realization from an ergodic process. As suggested in Ref. [29], we enforced ergodicity by sandwiching a turbid medium of suspended polystyrene balls in water immediately behind the object slab. From separate measurement of intensity autocorrelations of the sandwich and liquid medium, the intensity autocorrelation from the object slab alone is extracted using the method described in Ref. [29] as the “two-cell technique.” In our experiments, the ergodic medium is prepared to have an $l^{*}$ around 3 mm. For this, 63 μL of glycerol is first diluted with 6100 μL of distilled water that forms the base to which is added polystyrene balls of 2 μm diameter. The final solution has $\sim 0.02\%$ by volume of polystyrene balls.

In Fig. 3, the sample is a sandwich of the slab followed by the “ergodic” medium. The ultrasound transducers (V394- SU-F50MM-PTF, Olympus NDT NE, Inc) are aligned in such a way that their focal regions intersect at their thinnest mid-region to fall approximately in the middle of the slab. The transducers are driven synchronously using signals derived from a dual-channel function generator and power amplifiers, one at a frequency of 1.00 MHz, and the other at $1.00\mathrm{MHz}+\mathrm{\Delta }\varpi \mathrm{Hz}$ with provision to tune $\mathrm{\Delta }\varpi$ from a few hundred hertz to a few kilohertz. (In our experiments, $\mathrm{\Delta }\varpi$ is kept at 1 kHz.) As stated earlier, this arrangement helps us produce a low-frequency driving force from a low-frequency acoustic wave through mixing of the ultrasound waves in the common focal volume (the ROI), which would excite the inherently low-frequency natural modes of the vibrating ROI. The photons from a He–Ne laser (HRR 170 from Thorlabs, power 15 mW) illuminating the sandwich are captured using a single-mode fiber and given to the PMT (Hamamatsu, H7360-03). The fiber is carefully aligned to capture a single speckle formed through interference of the emerging photons. The autocorrelator to which the signal from the PMT is given after pre-processing can sample up to 16.7 ns (DAC, Flex 021d [28]).

What we measure in the experiments is the intensity autocorrelation measured at a point, $\mathbf{r}$, in the object, denoted by $g_2(\mathbf{r},\tau )$, where $\tau$ is the delay time. The effect of the low-frequency forcing confined to the ROI is a periodic modulation in $g_2(\mathbf{r},\tau )$ as it decays to small values. As one closely examines this modulation, it is observed that the modulation itself decays with $\tau$, as revealed by the strength of, say, the first fundamental frequency of the modulation [16]. This decay in modulation depth is in fact owing to the Brownian motion in the location tagged by the intersecting ultrasound beams, i.e., the ROI. This is because the ultrasound forcing at low frequency and the periodic movement of scattering centers introduced by it are available only in the ROI. Therefore, if one recovers the MSD from the decay of modulation depth versus $\tau$, then it must belong only to the particles in the ROI. This way of carrying out DWS, called local DWS pertaining to a region marked by ultrasound tagging, was established in Ref. [16].

We have conducted experiments with the three prepared samples of agar slabs of increasing stiffness, and also a slab of pork fat obtained from the hind side. We have recorded $g_2(\mathbf{r},\tau )$ from sandwiches of these samples with a cuvette containing the ergodic scattering medium. Prior to collecting these data sets, experiments were carried out to record $g_2(\mathbf{r},\tau )$ for the ergodic medium alone. From these sets of measurements, following the procedure outlined in Ref. [29], $g_2(\mathbf{r},\tau )$ corresponding to the three agar slabs and the pork fat are obtained. From $g_2(\mathbf{r},\tau )$, the modulus of amplitude autocorrelation, $|g_1(\mathbf{r},\tau )|$, is extracted using Siegert’s relation, which is $g_2(\mathbf{r},\tau )=1+f|g_1(\mathbf{r},\tau )|$ [31]. The coupling constant, $f$, is obtained from the value of $g_2(\mathbf{r},\tau )$ at $\tau =0$. As mentioned earlier, these $g_2{(\mathbf{r},\tau )}^{\prime }s$, as well as $g_1{(\mathbf{r},\tau )}^{\prime }s$ have a modulation overriding them. We used short-time Fourier transform with a window of suitable length such that a minimum of a few cycles of the periodic modulation is within the window to extract the strength of the fundamental frequency in the modulation in $g_1(\mathbf{r},\tau )$, $M(\mathbf{r},{\tau }_i)$ at ${\tau }_i$ (also called the evolutionary power spectrum [32]), the center of the window.

In a standard DWS experiment, $g_2(\mathbf{r},\tau )$, obtained from the hard-wired autocorrelator, has to cover a large range of delay time, starting from a few microseconds to tens of seconds; a (quasi) logarithmic scale is automatically applied in the $\tau$ axis, selecting the so-called multi-tau scheme to process the intensity data. However, to obtain $M(\mathbf{r},\tau )$ from $g_1(\mathbf{r},\tau )$, we need uniform sampling of the input beyond the Nyquist rate set by the fundamental frequency of modulation in the autocorrelation. To obtain uniform sampling, the multi-tau scheme is replaced by the “photon-mode” scheme available with the autocorrelator (DAC, Flex 021d [28]). This scheme can acquire the intensity data at any (uniform sampling time as small as 16.7 ns. A Matlab 1-d autocorrelation routine computes the autocorrelation of acquired intensity data. We note that the 1 kHz frequency of modulation superimposing $g_2(\mathbf{r},\tau )$ is large enough to extract the low-frequency resonant modes we are expecting without violating the sampling theorem.

From the agar slabs and pork fat, $g_2(\mathbf{r},\tau )$ was experimentally measured, from which $g_1(\mathbf{r},\tau )$ and $M(\mathbf{r},\tau )$ are computed as described above. $M(\mathbf{r},\tau )$ is related to the MSD ($⟨\mathrm{\Delta }{r}^2(\mathbf{r},\tau )⟩$) through [33]

In Eq. (6) $L$ is the average thickness of the ROI, and $k_0$ is the modulus of the illumination wave vector. This equation is inverted for $⟨\mathrm{\Delta }{r}^2(\mathbf{r},\tau )⟩$.

4. RESULTS

A. Verification of GLE with Experimental MSD Data

Representative plots of $⟨\mathrm{\Delta }{r}^2(\mathbf{r},\tau )⟩$ versus $\tau$, for pork fat, are shown in Fig. 4. Also shown is the corresponding plot obtained from $g_1(\mathbf{r},\tau )$, i.e., without the ultrasound. It is seen that the plateau region is without fluctuations when the ultrasound forcing is removed. To verify our theoretical model given in Eq. (3), we have solved it for $⟨\mathrm{\Delta }{r}^2(\mathbf{r},\tau )⟩$ within a Monte Carlo setup [17]. Using a Prony series approximation for the damping memory kernel, the unknown coefficients of the series are approximated, using an inverse procedure, with the help of the experimentally available MSD plot. With this, we have been able to convert a set of non-Markovian stochastic differential equations (SDEs) into a set of Markovian ones. A nonlinear filtering scheme is employed, as described in Ref. [34], to extract the unknown Prony series coefficients and the strength of the multiplicative noise term. This exercise leads to a solution of $⟨\mathrm{\Delta }{r}^2(\mathbf{r},\tau )⟩$, closely matching its experimental counterpart (see Fig. 5). The strategy of the filtering scheme was to drive the error between simulated and experimental MSDs to a zero-mean, pure Brownian noise (which in a more general technical setting is called a martingale) by manipulating the unknowns, and using the above SDE as the process equation.

 

Fig. 4. Variation of MSD with $\tau$ from the measured $g_1(\tau )$ from the pork slab. The thin curve corresponds to when the ultrasound is on, and the thick curve to when it is off.

Download Full Size | PDF

 

Fig. 5. Comparison of MSD versus $\tau$ plot from a pork slab, obtained experimentally (thick curve) with the corresponding one obtained through solving the GLE [Eq. (3)] (thin curve).

Download Full Size | PDF

B. Experimental Measurement of Resonant Modes and Inversion for Elastic Constants Using an Ensemble Kalman Filter

The MSD plots are computed from experimentally obtained $M(\tau )$ plots corresponding to the three agar slabs and pork fat (see Fig. 4). The noisy plateau regions are Fourier transformed, which revealed spectral peaks, $\mathrm{\Delta }{\varpi }_i,i=1,\dots ,4$, at various frequency locations [Figs. 6(a)–6(d)]. Our present objective is to invert these peak frequencies for the elastic properties of the material of the ROI, namely, the three Young’s moduli and the three shear moduli, as appearing in the forward equation, Eq. (4). Before doing the inversion, for which the ensemble Kalman filter (EnKF) is made use of, we verify the experimental data, for which all the elastic constants are assumed to be known, by computing the modal frequencies using Eq. (4). The comparison is in Table 1. It is observed that the experimentally measured spectral peaks correspond to some of the modal frequencies. We notice that the ultrasound forcing is not able to excite all the computed spectral peaks. We now attempt an inversion of the excited, experimentally measured spectral peaks for the elastic parameters.

Table 1. Verification of the Experimentally Measured Resonant Modes of the Objects through Those Computed Using a Modal Analysis [23]^a

View Table | View all tables in this article

 

Fig. 6. (a) Power spectrum obtained from the plateau of the experimentally obtained MSD versus $\tau$ curve from 0.5% agar slab, displaying the peaks (resonant modes) at 75.5, 78.7, 84.2, and 90.4 Hz. (b) Same as (a), but for 1% agar slab. The peaks are at 108.0, 110.0, 135.0, and 142.0 Hz. (c) Same as (a), but for 1.5% agar slab. The peaks are at 135.8, 151.0, 165.0, and 191.4 Hz. (d) Same as (a), but for the pork slab. The peaks are at 214.8, 228.9, 253.9, and 269.5 Hz.

Download Full Size | PDF

1. Inverse Problem

To invert the measured eigenvalues (natural frequencies) for the mass and stiffness matrices ($\mathbf{M}$ and $\mathbf{\Gamma }$), and through them the elastic constants, we go back to Eq. (4), the eigenvalue problem. We rearrange Eq. (4) as

For a nontrivial solution of Eq. (7), the following condition has to be met:

The left-hand side of Eq. (8) is a polynomial in $\omega$ of degree $N$, where $N$ is the number of DOFs associated with the eigenvalue problem of Eq. (4). Assuming that the entries in $\mathbf{M}$ are known, the (ill-posed) inverse problem we want to tackle is the recovery of elastic constants that make up the entries in $\mathbf{\Gamma }$ from the four natural modes of vibration we measured through DWS. We pose this as a stochastic filtering problem and solve it using the EnKF [35]. The motivation for and basics of the approach are briefly given in the next section, and for details the reader is referred to [35]. We reiterate that the material is considered orthotropic with nine unknowns, out of which three, the Poisson’s ratios, are considered not as important sets of parameters as the other six, on which we concentrate. Moreover, the interpretation of the elastic constants is dependent on the selection of the coordinate axes.

2. Reconstruction Based on Stochastic Evolution

The inversion technique based on the evolutionary stochastic method involves defining the unknown nine elastic parameters and the four measured resonant frequencies as diffusive stochastic processes in an appropriate probability space. The motivation to select stochastic filtering for this job can be made clear from the following arguments. (1) Recovery of elastic constants from the sparse data of just a few natural frequencies is deterministically ill-posed. For instance, use of a regularized Gauss–Newton scheme may possibly result in a number of equally admissible solutions, depending on the regularization, selection of which in itself is elaborate and time-consuming [36]. (2) Deterministic methods are unable to rationally account for the model and measurement error/noise in the estimated elastic constants. (3) In the stochastic filtering approach, noise is an integral part of the model, and the goal of filtering is to guide the evolving parameters to the “best” estimate, conditioned on the noisy measurements up to the current instant. In the present context, the goal of the filtering problem is as follows. First, we define the unknown elastic constants by a vector stochastic process (say, ${\mathbf{X}}_t$) evolving over pseudo-time, which is a monotonically increasing function of the iteration number in the algorithm. The elastic constants are, in fact, stationary, but unknown. In order to estimate them, we allow them to evolve as diffusive stochastic processes in pseudo-time,

5. DISCUSSION

We have used the above algorithm to reconstruct the nine independent orthotropic elastic constants that characterize the material. It is noticed that out of the nine, three parameters—viz., Poisson’s ratio in the three directions—remained unchanged around the initial guess corresponding to the nearly incompressible material in our experiments, beyond a few initial recursions. Therefore, even though we ran the algorithms with nine unknowns, for all practical purposes, we had only six: three each of the Young’s and shear moduli. If we assumed full anisotropy and needed to recover 21 independent elastic constants, our experience was that, with just four measurements, the posterior or filtering density had too many peaks to enable meaningfully converged estimates of the unknown parameters. Even for the orthotropic case, in order to ensure proper convergence, it was necessary to impose a priori (and physically admissible) constraints on set of parameters. For physical realizability of orthotropic parameters, the following specific constraints need to be met [25]:

A simple strategy to always be within the constraints specified above is to start the algorithm with a “good” guess of the parameters. The good initial guess was obtained by solving the problem, initially assuming isotropic properties, leading to good guesses of the Young’s and shear moduli and Poisson’s ratio. Figure 7 shows the evolution of the Young’s moduli along the three directions with the recursion number, and Fig. 8 the same for shear moduli. Insets in these figures show the initial evolution of the elastic parameters, assuming isotropic behavior. The evolution of the elastic parameters, shown in Figs. 7(d) and 8(d), shows that for the harvested pork back fat, these are dependent on direction. The same for agar phantoms shows independence with respect to direction, revealing the material to be isotropic [Figs. 7(a)–7(c) and 8(a)–8(c)]. The coordinate system, showing axes 1, 2, 3 against a photograph of the fat sample, is shown in Fig. 9. The appearance of the end of a 1D rod-like structure running perpendicular to the surface and into the slab is clearly seen. This is a reinforcement providing extra protection in the direction marked 1. The Young’s modulus along this direction, noted as $E_1$, has a higher value of 205.8 kPa, as reconstructed from the measured natural frequencies, compared to $E_2$ and $E_3$, which are 146.6 kPa and 85.2 kPa, respectively. The values of shear moduli are 33.1 kPa ($G_{23}$), 38.2 kPa ($G_{13}$), and 48.1 kPa ($G_{12}$) (notice that $G_{12}$ is higher). We also note, in passing, that pork fat harvested from other places away from the skin, such as from near internal organs, manifests isotropic behavior. Table 2 gives a summary of the reconstructed elastic parameters for all the slabs: three for agar and one for pork fat. In the case of the agar, we have cross-verification with measurements using a rheometer. The comparison for pork fat is from results published in Ref. [40], wherein the material is assumed to be isotropic, giving a single average value just for the Young’s modulus, which is 140  kPa, close to one of the values reconstructed (146 kPa).

Table 2. Comparison of Elastic Parameters Estimated by the Present Method with Measurements from Rheometer or Published Results, as the Case May Be

View Table | View all tables in this article

 

Fig. 7. Estimated Young’s moduli, $E_1$, $E_2$, and $E_3$ from the measured resonant modes for slabs (a) with 0.5% agar, (b) with 1% agar, (c) with 1.5% agar, and (d) pork fat. Notice that for pork, being orthotropic, $E_1$, $E_2$, and $E_3$ are different at convergence. For all four cases, the object was modeled to be linear, orthotropic material. Inset shows the initial convergence with assumption of isotropy.

Download Full Size | PDF

 

Fig. 8. Same as Fig. 7, but showing shear moduli (a) with 0.5% agar, (b) with 1% agar, (c) with 1.5% agar, and (d) for pork fat. Notice that for pork, being orthotropic, $G_{23}$, $G_{13}$, and $G_{12}$ are different at convergence. For all four cases, the object was modeled to be linear, orthotropic material. Inset shows the initial convergence with assumption of isotropy.

Download Full Size | PDF

 

Fig. 9. Pork back fat used during the experiment with the coordinate axis (inset), representing the material directions. The surface visible in this figure was just below the skin, as can be confirmed from its grainy appearance.

Download Full Size | PDF

Figure 10 shows the root mean square error (RMSE) with iteration, the error as computed in the measurement domain, i.e., from the modal frequencies. The RMSE ($f_k$) for each step $k$ was calculated as ${\mathrm{RMSE}}_k=\sqrt{\sum _{i=1}^4{(f_i^{\mathrm{experiment}}-{\widehat{f}}_i^k)}^2}$. Note that in the first 20 iterations when the material was assumed to be isotropic, the error stagnated from around recursion 3 onwards. (See also the insets in Figs. 7 and 8.) From recursion 20 onwards, we had assumed orthotropy, and the error began falling until it reached acceptably low values beyond 100 recursions. For agar, which is isotropic, the behavior of the error curve is quite different, showing complete stagnation en route to convergence in the first three iterations, when modeled as isotropic, and no further improvements when assumed to be orthotropic (Fig. 10). Modeling the material as truly anisotropic, the convergence of error in the case of pork fat would have been faster and to lower values, for the tissue samples generally do exhibit complete anisotropic behavior.

 

Fig. 10. Behavior of root-mean-square error with iteration number during the recovery of the orthotropic elastic parameters of pork fat. During the first 20 iterations, the sample was modeled isotropic, and in the next 100, orthotropic.

Download Full Size | PDF

Another important source of error in the parameter recovery is from the possible unreliability of the measured data, which are the modal frequencies. Unlike in our previous work [12], our measurement here is not from the VA wave propagated to the boundary through the background. Even if it is, the location of modal frequencies will be influenced by the “transfer function” of the background if the VA wave is modulated by a multiplicative term containing $e^{j{\omega }_{0}x}$, which will shift the VA spectrum by ${\omega }_0$. This is unlikely from a biological tissue forming the background. However, in our measurement we use coherent photons to transport the ROI dynamics to the boundary. As proven in Ref. [16], the ultrasound-assisted DWS provides a “local” measurement from the modulation depth decay, and the modulation itself created by the location specific interference of two ultrasound beams. Therefore, the present DWS-based measurement gives us variation of $\mathrm{\Delta }{x}^2(\mathbf{r},\tau )$ pertaining to the ROI. However, there can be errors in measurement of modulation depth, which are passed on in greater measure to the retrieved $\mathrm{\Delta }{x}^2(\mathbf{r},\tau )$ and also the modal frequencies computed. This is reflected in the comparison of simulated and measured natural frequencies given in Table 1, where in places the match is not quite close. An analysis of error of elastic property recovery owing to variation in the measured resonant frequencies is given in Table 3. It is seen that the error introduced in elastic parameter due to errors in modal frequency measurement is more pronounced for lower-order modes than higher-order modes. The maximum error is in $G_{12}$ (10 kPa), owing to an error of 5 Hz in the measurement of the first modal frequency, and the minimum is in $E_1$ (3 kPa) for the same error in the measurement of the fourth modal frequency.

Table 3. Error in Recovered Young’s and Shear Moduli of Pork Fat Owing to Error in the Measured Modal Frequencies

View Table | View all tables in this article

6. CONCLUSION

We have demonstrated a method to noninvasively measure elastic constants in tissue and tissue-like phantoms when they are assumed to have orthotropic material behavior. The method uses focused ultrasound beams for remote palpation, which induces vibro-acoustic waves from its common focal volume (the ROI). Measurement of the resonant modes of vibration of the ROI has enabled us to recover the Young’s and shear moduli from the eigenvalue equation of the freely vibrating ROI. In this work, we have successfully recovered some of these natural frequencies using DWS-based measurement of MSD growth with respect to time. The Fourier transform of the fluctuations in the plateau of the MSD curve contained spectral peaks corresponding to ROI modes excited by the ultrasound forcing. With the measured modal frequencies, the eigenvalue equation is inverted for the elastic constants using a stochastic filtering method, square-root filter implementation of the ensemble Kalman filter. The recovered Young’s and shear moduli, the result of convergence reached at the end of approximately 120 recursions, demonstrate how the present scheme works. For the nearly isotropic agar gel, the three Young’s moduli converged to a single value representing the isotropic material constant; this was also true of the three shear moduli. These isotropic values were verified through independent rheometer measurement. For pork fat, which has orthotropic behavior, we reconstructed six elastic constants. One of the Young’s modulus values is seen to be close to that reported elsewhere, assuming pork fat to be isotropic.

To the best of our knowledge, this is the first instance where resonant ultrasound spectroscopy has been employed for anisotropic characterization of soft tissues. Estimation of anisotropic elastic constants provides the physician with extra information, such as which of the constants are sensitively affected by progress of tumor or interventive treatment methods. However, many more investigative trials on real tissue, healthy or diseased, must be undertaken before this method can become a potentially useful diagnostic imaging tool.