Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2704493

Formats

Article sections

- Abstract
- I. Introduction
- II. Modeling Relaxation and Creep
- III. Data Acquisition and Processing
- IV. Results
- V. Conclusion
- References

Authors

Related links

IEEE Trans Ultrason Ferroelectr Freq Control. Author manuscript; available in PMC 2010 April 1.

Published in final edited form as:

PMCID: PMC2704493

NIHMSID: NIHMS122020

The authors are with the Departments of Bioengineering and Electrical and Computer Engineering, Beckman Institute for Advanced Science and Technology, University of Illinois at Urbana-Champaign, Urbana, IL 61801 (e-mail: ude.sionilli@2ppayr)

The publisher's final edited version of this article is available at IEEE Trans Ultrason Ferroelectr Freq Control

The viscoelastic response of hydropolymers, which include glandular breast tissues, may be accurately characterized for some applications with as few as 3 rheological parameters by applying the Kelvin-Voigt fractional derivative (KVFD) modeling approach. We describe a technique for ultrasonic imaging of KVFD parameters in media undergoing unconfined, quasi-static, uniaxial compression. We analyze the KVFD parameter values in simulated and experimental echo data acquired from phantoms and show that the KVFD parameters may concisely characterize the viscoelastic properties of hydropolymers. We then interpret the KVFD parameter values for normal and cancerous breast tissues and hypothesize that this modeling approach may ultimately be applied to tumor differentiation.

Mammography and sonography are often the initial imaging techniques applied to patients when breast tumors are to be diagnosed. Both modalities are considered able to better detect the presence of a lesion than to differentiate between benign and malignant lesions [1], [2]. Nearly all focal breast lesions appear as hypoechoic regions in sonograms; diagnosis requires careful visual evaluation of the lesion boundary, image texture, and shadow features. To help improve discriminability, elasticity imaging was proposed because of its ability to reveal the presence of a stiff desmoplastic reaction specific to regions surrounding some malignant tumors [3]. Recent clinical experience in examining diagnostic performance of elasticity imaging for breast disease is very encouraging [4]–[6].

Static elasticity imaging techniques generate strain maps describing deformation patterns resulting from small, quasi-static, uniaxial compressions applied to the tissue surface as the tissue is scanned ultrasonically. However, strain images are just the tip of the diagnostic information iceberg provided by mechanical properties. It is now known from molecular biology studies that cancerous epithelial cells must send and receive molecular signals to and from surrounding stromal cells if they are to develop into malignant tumors [7]. To effect malignant transformation, tumorigenic signaling pathways induce structural modifications to the extracellular matrix (ECM) and alter the viscosity of extracellular fluids, thus inducing changes to viscoelastic properties of the tissue. The most noticeable change is the palpable stiffening of breast stroma. Changes in the cellular mechanical environment can profoundly influence the progression of disease [8], [9]. Therefore, time-varying features of the viscoelastic (VE) response of breast tissues to gentle deforming forces may be an important untapped source of information about both the biology of the malignant processes and the medicine of diagnosis and treatment monitoring.

Radio frequency (RF) echo signals from pulse-echo ultrasonic imaging systems are often used to image the strain response of tissues to applied stress stimuli [10]. Information about the VE properties of tissues is found from spatiotemporal variations in the corresponding strain fields. We form static strain images by processing ultrasonic RF echo frame sequences using a multi-resolution cross-correlation-based displacement algorithm or a regularized optical-flow algorithm [11]; the algorithm choice depends on the amount of applied deformation. Both algorithms were designed to minimize strain noise for relatively large (>1%) deformations. Later generation algorithms [12], [13] provide superior strain estimates for the very small displacements (<0.1%) associated with radiation force stimuli [14] and viscoelastic creep imaging techniques [15]–[17]. For each volume element of tissue imaged, a time sequence of strain estimates is measured from which VE imaging parameters are extracted.

Viscoelastic imaging parameters are selected from the parameters of rheological models applied to time-varying strain estimates. Model parameters summarize material properties of multiphasic polymeric media such as hydrogels and breast stroma. The best rheological models for our application yield just a few imaging parameters that are descriptive of the deformation physics and also are directly connected to the biology of disease. We found that a biphasic poroviscoelastic theory originally proposed by Mak [18], [19] and later modified by Suh and DiSilvestro [20] for articular cartilage could be adapted to reasonably represent the material behavior of hydrogels [17] and breast stroma [21]. The various components of the medium are grouped into 2 phases: a solid phase, consisting of a collagenous ECM with associated glycoproteins [22] and cells, and a fluid phase consisting of an interfibrillar liquid. Gelatin hydrogels are structurally simpler than stromal tissues. Gels are an aggregate matrix of denatured type I collagen with electrically charged molecular side chains. The charged side chains act to “structure” nearby water molecules in a manner similar to the function of glycoproteins in stroma. The fluid viscosity thus varies depending on matrix density and structure.

In both tissues and gels, the solid matrix forms an ensemble of different size channels or pores through which fluid moves when the medium is loaded. However, the solid matrix itself behaves viscoelastically by virtue of the numerous hydrogen-bonded cross links between matrix fibers. The strain response of a hydropolymer medium subjected to a step stress (creep experiment) has been described by a generalized Kelvin-Voigt model consisting of the sum of many real exponential functions [23]. While such a model is easily adaptable to random pore sizes and variations in matrix cross-link strength, it requires at least 2 parameters for each Voigt unit to fully represent the mechanical response. The dimensionality of this straightforward model is very large.

We describe a method for reducing the dimensionality of the generalized Kelvin-Voigt model by applying singular-value decomposition methods [24]. Imaging parameters were selected from the largest eigenvalues of the creep curve and hence represented the dominant features of the polymer response. The complex compliance spectra from breast stroma and hydrogels place them in a class known as amorphous, weakly cross-linked polymers with distinct solid and fluid spectral responses [23]; therefore a 2-component Kelvin-Voigt model was proposed. Parameters were then associated with the 2 phases. However, spectral separability of the phases depends significantly on boundary conditions. For the unconfined, uniaxial compressions usually applied in static strain imaging techniques, overlap of the bimodal response was often observed. Clearly, we need a rheological model that considers the continuous nature of mechanisms governing the slow relaxation of polymers. In experimental situations, where no clearly dominant mechanistic components emerge from creep data, the rheological model should provide a concise set of model parameters for imaging.

Viscoelastic models based on fractional derivatives (FD) were introduced by Sloninsky in 1967 [25] to find a parsimonious representation for complex media. The mechanical responses predicted by such models were found to be consistent with the molecular theory of polymers by Bagley and Torvik [26]. Schiessel and Blumen [27] derived hierarchical mechanical analogs to fractional derivative elements and models by assembling numerous springs and dashpots (elastic and viscous terms, respectively) in series and parallel. The FD models of increasing complexity were proposed [28] to simulate the rheological behavior of synthetic polymers as well as biological cells and tissues [29]–[31]. Recently, others contributed statistical, physical, and mathematical justifications for FD applications [32]–[34].

Our task in this report is to review the extensive literature on Kelvin-Voigt fractional derivative (KVFD) rheological models for the study of hydropolymer dynamics at low applied force frequencies. We then develop methods where the imaging parameter can be interpreted in terms of the biphasic properties.

We begin by reviewing how fractional derivatives arise naturally from slowly varying distributions of exponential relaxations (Kelvin-Voigt units) in the VE response. We then apply a 3-parameter KVFD model to represent the dynamic behaviors of hydrogels and *in vivo* breast tissues. Images of the corresponding FD parameters for simulated hydropolymers, experimental phantoms, and *in vivo* breast tissues are generated and interpreted in terms of the stiffness and fluidity of the material.

In a stress-relaxation experiment, the force response of a sample to a known applied deformation is measured versus time. The FD model assumes the polymer network undergoing viscoelastic relaxation can be modeled by the linear superposition of many elementary viscoelastic (Maxwell) units. The *n*-th viscoelastic unit is given by the classic Boltzmann superposition principle that relates stress *σ _{n}*(

$${\sigma}_{n}(t)={\int}_{0}^{t}d{t}^{\prime}{E}_{n}(t-{t}^{\prime}){\dot{\epsilon}}_{n}({t}^{\prime}),$$

(1)

where we assume each Maxwell unit is linear, time-invariant, and initially at rest. The modulus describes material properties of the *n*-th polymer component through 2 constants, *E′ _{n}* and

$${E}_{n}(t)={{E}^{\prime}}_{n}\phantom{\rule{0.3em}{0ex}}{e}^{-t/{\tau}_{n}},$$

where *E′ _{n}* is the modulus and

The relaxation modulus for the polymer is the sum of contributions from all *N* units [34],

$$E(t)=\sum _{n=0}^{N}{{E}^{\prime}}_{n}\phantom{\rule{0.3em}{0ex}}{e}^{-t/{\tau}_{n}},$$

(2)

where unit constants are assumed to vary slowly with *n*. This interdependency can be modeled recursively with 2 scaling parameters *a* and *b*:

$${{E}^{\prime}}_{n}=b{{E}^{\prime}}_{n-1}={b}^{n}{{E}^{\prime}}_{0}$$

(3)

$${\tau}_{n}=\frac{{\tau}_{n-1}}{a}=\frac{{\tau}_{0}}{{a}^{n}}.$$

(4)

$$E(t)={{E}^{\prime}}_{0}\sum _{n=0}^{N}{b}^{n}{e}^{-{a}^{n}t/{\tau}_{0}}={{E}^{\prime}}_{0}\sum _{n=0}^{N}{e}^{F(n,t)},$$

(5)

where *F*(*n,t*) = *n* ln *b* − *a ^{n}t*/

It is observed with polymers that the time dependence of the elastic modulus *E*(*t*) is a power law. Let *n*_{0} indicate the viscoelastic unit where *e ^{F}*

$$F(n,t)\simeq F({n}_{0},t)+\frac{{(n-{n}_{0})}^{2}}{2}{\left(\frac{{\partial}^{2}F}{\partial {n}^{2}}\right)}_{{n}_{0}},$$

(6)

since (*F/n*)*n*_{0} = 0. Combining the definition of *F*(*n,t*) in (5) with its first derivative being zero, we find that

$${n}_{0}=\frac{\text{ln}(\alpha {\tau}_{0})-\text{ln}\phantom{\rule{0.2em}{0ex}}t}{\text{ln}\phantom{\rule{0.2em}{0ex}}a},\phantom{\rule{0.5em}{0ex}}\text{where}\phantom{\rule{0.2em}{0ex}}\alpha =\frac{\text{ln}\phantom{\rule{0.2em}{0ex}}b}{\text{ln}\phantom{\rule{0.2em}{0ex}}a}.$$

(7)

We restrict the range of acceptable model constants *a* and *b* to those producing physical models of elastic moduli with power-law form. Realistic models require *τ _{n}* > 0 for all

$$\underset{a\to 1,N\to \infty}{\text{lim}}E\left(t\right)=\frac{{E}_{0}^{\prime}}{1-b}{e}^{-t/{\tau}_{0}}$$

is a pure exponential. Regardless, we note that *α* is defined for the range of constants 0 < *a*, *b* < 1.

We now posit that a continuum of viscoelastic units contributes to the elastic modulus, *E*(*t*), i. e., *N* → ∞. A continuum model is consistent with the hypothesis that the viscoelastic response of the collagen matrix results from relaxation of hydrogen bonds having a continuous spectrum of bond strengths [23]. Combining (5)–(7), we follow Kawada [34] to obtain the power-law form,

$$\begin{array}{c}E(\text{t})\simeq {{E}^{\prime}}_{0}\sum _{n=0}^{\infty}{e}^{\phantom{\rule{0.1em}{0ex}}F({n}_{0},\text{t})+[{(n-{n}_{0})}^{2}/2]{[({\partial}^{2}F)/(\partial {n}^{2})]}_{{n}_{0}}}\hfill \\ \phantom{\rule{2em}{0ex}}={{E}^{\prime}}_{0}\sum _{n=0}^{N}{e}^{\phantom{\rule{0.1em}{0ex}}\alpha (\text{ln}(\alpha {\tau}_{0})-\text{ln}(t)-1)-[{(n-{n}_{0})}^{2}/2]\alpha {(\text{In}\phantom{\rule{0.2em}{0ex}}a)}^{2}}\hfill \\ \phantom{\rule{2em}{0ex}}={{E}^{\prime}}_{0}{e}^{\phantom{\rule{0.1em}{0ex}}\alpha (\text{ln}(\alpha {\tau}_{0})-1)}{t}^{-\alpha}\sum _{n=0}^{N}{e}^{-\alpha {(\text{ln}a)}^{2}[{(n-{n}_{0})}^{2}/2]}\hfill \\ \phantom{\rule{2em}{0ex}}\simeq {{E}^{\prime}}_{0}{e}^{\phantom{\rule{0.1em}{0ex}}\alpha (\text{ln}(\alpha {\tau}_{0})-1)}\sqrt{\frac{2\pi}{\text{ln}a\text{ln}b}{t}^{-\alpha}}\hfill \\ \phantom{\rule{2em}{0ex}}\simeq \Phi {t}^{-\alpha},\hfill \end{array}$$

(8)

where all time-independent factors are collected into Φ. On a log-log scale, *E*(*t*) is a straight line with slope −*α* and intercept ln Φ. Applying unitless parameters, we find the alternative form

$$\frac{E(t)}{{E}_{1}}=\frac{1}{\Gamma (1-\alpha )}{\left(\frac{t}{\tau}\right)}^{-\alpha}\phantom{\rule{0.2em}{0ex}},$$

(9)

where *E*_{1} is the scaled modulus for unit *n*_{0}, *τ* is a characteristic relaxation time constant, and Γ is the Gamma function for non-integer arguments.

Power-law viscoelastic responses suggest properties of self-similarity and memory, which are the features described by fractional derivatives. The Riemann-Liouville definition of a fractional derivative operator applied to a function *f* is

$${a}_{}{\mathcal{D}}_{t}^{-\gamma}(f)\triangleq \frac{1}{\Gamma (\gamma )}{\int}_{a}^{t}d{t}^{\prime}f({t}^{\prime}){\phantom{\rule{0.2em}{0ex}}(t-{t}^{\prime})}^{\gamma -1}.$$

(10)

Combining (1), (9), and (10), the net constitutive equation may be written as shown in the Appendix.

$$\sigma (t)={E}_{1}{\tau}^{\alpha}\frac{{d}^{\alpha}\epsilon}{d{t}^{\alpha}}.$$

(11)

Eq. (11) is more intuitive when expressed in the radial frequency domain, *ω.* Applying the Fourier derivative theorem we find

$$\sigma (\omega )={E}_{1}{\tau}^{\alpha}{(i\omega )}^{\alpha}\epsilon (\omega ),$$

(12)

from which the complex modulus and compliance are, respectively, defined as

$${E}^{\ast}(\omega )\triangleq \frac{\sigma (\omega )}{\epsilon (\omega )}={E}_{1}{(i\omega \tau )}^{\alpha}\phantom{\rule{0.2em}{0ex}}\text{and}\phantom{\rule{0.2em}{0ex}}{D}^{\ast}(\omega )=\frac{1}{{E}_{1}{(i\omega \tau )}^{\alpha}}.$$

(13)

Eq. (13) defines FD units referred to as “spring-pots” because of their ability to combine large assemblies of Maxwell or Voigt (exponential) units to predict non-exponential relaxation and creep behavior, respectively. Advantages of this model include just 3 parameters, *α*, *τ*, and *E*_{1}, that have straightforward physical interpretations. When *α* → 0, the fractional unit behaves like a Hookean spring; when *α* → 1, it behaves as a Newtonian dashpot. For intermediate values of *α*, it behaves as a viscoelastic material. Though a spring-pot could more accurately be modeled by a hierarchical assembly of springs and dashpots [27], our goal here is to keep a concise feature space. Unfortunately, a spring-pot alone cannot model the relaxation and creep responses typically observed in polymers, which can plateau over time, and the model must be expanded further.

Adding a spring in parallel with the spring-pot unit allows for the relaxation plateau to be observed experimentally. This is the KVFD model given by

$$\sigma (t)={E}_{0}\epsilon (t)+{E}_{1}{\tau}^{\alpha}\frac{{d}^{\alpha}\epsilon (t)}{\phantom{\rule{0.2em}{0ex}}d{t}^{\alpha}},$$

(14)

and consisting of elastic and viscoelastic terms. *E*_{0} and *E*_{1} may have different values; however, *τ* can be modified so that

$$\sigma (t)={E}_{0}\epsilon (t)+{E}_{0}{{\tau}^{\prime}}^{\alpha}\frac{{d}^{\alpha}\epsilon (t)}{d{t}^{\alpha}}.$$

(15)

This form is convenient when expressing the time-domain behavior in response to a step stimulus.

Considering stress-relaxation data, where a uniaxial compressive step strain is applied to a sample, we have *ε˙*(*t*) = *ε*_{0} *δ*(*t*) and from (15)

$$\frac{\sigma (t)}{{\epsilon}_{0}}=E(t)={E}_{0}+\frac{{E}_{0}}{\Gamma (1-\alpha )}{\left(\frac{t}{{\tau}^{\prime}}\right)}^{-\alpha}.$$

(16)

Because we measure strain ultrasonically and not stress, we usually conduct creep experiments to image viscoelastic properties. In these experiments, a constant step stress is applied to the sample. Schiessel *et al*. [37] showed that the compliance expression for the KVFD model corresponding to the modulus expression in (16) is

$$\frac{\epsilon (t)}{{\sigma}_{0}}=D(t)=\frac{1}{{\text{E}}_{0}}{\left(\frac{t}{{\tau}^{\prime}}\right)}^{\alpha}{E}_{\alpha ,1+\alpha}\left(-{\left(\frac{t}{{\tau}^{\prime}}\right)}^{\alpha}\right),$$

(17)

where ** E** is a type of hypergeometric function known as the generalized Mittag-Leffler function [38], (also, see the Appendix). Fig. 1 displays the modeled time behavior of compliance according to (17) for values of

$${\int}_{-\infty}^{t}E(t-{t}^{\prime})D({t}^{\prime})d{t}^{\prime}=t.$$

(18)

Computation is simplified by transforming (14) into the Fourier domain resulting in the complex modulus

$${\text{E}}^{\ast}(\omega )={E}_{0}+{(i\omega {\tau}^{\prime})}^{\alpha}{E}_{0}={E}_{0}\phantom{\rule{0.2em}{0ex}}(1+{({e}^{i\pi /2}\omega {\tau}^{\prime})}^{\alpha}).$$

(19)

The loss modulus, which describes the deformation energy dissipated as heat, is given by the imaginary part of the complex modulus:

$$\mathfrak{J}{E}^{\ast}(\omega )\triangleq {E}^{\u2033}(\omega )={(\omega {\tau}^{\prime})}^{\alpha}{E}_{0}\text{sin}\left(\frac{\alpha \pi}{2}\right).$$

(20)

Plotting ln *E*″(*ω*) versus ln *ω*, we find a linear function with positive slope *α*. Curve fitting is therefore significantly faster in the frequency domain than in the time domain.

The equations above show that the viscoelastic parameters *E*_{0}, *τ′*, and *α* may be estimated in the frequency domain from either the loss modulus *E″*(*ω*) or the loss compliance *D*″(*ω*), depending on whether stress relaxation or creep experiments are employed. Estimates can be obtained directly from time-domain data at considerably greater computational cost. Time-domain processing is preferred for short-duration strain acquisitions (*in vivo* breast data) because of lower noise and parameter estimation uncertainty.

Typical computational times varied depending on the strain acquisition time and whether time- or frequency-domain analysis is applied. Using a single node on a PC cluster, parametric images processed in the frequency domain were formed in about 60 s for a 100 × 100-sample strain image sequence 20 s in duration recorded at 4 frames per second (fps). The same data processed on the same cluster node but using a time-domain algorithm required more than 3 h of processing time. Increasing the strain acquisition to 2000 s while reducing the frame rate downward from 4 fps as time increases, requires an average of 5 min of processing time in the frequency domain and several days in the time domain. Computational times were the same for processing simulated, phantom, and tissue echo data; the results depended only on the volume of data.

Echo data were acquired during creep experiments as illustrated in Fig. 2. A linear array transducer was used to apply a step force at the top surface and maintain it constant while ultrasonic RF echo frames were recorded. A Siemens Sonoline Antares system (Siemens Medical Solutions USA, Inc., Mountain View, CA) with an ultrasonic research interface was used for all experiments. A VF10-5 linear array transducer transmitted broadband 8 MHz pulses. Echo frames were recorded beginning immediately before the application of a 1 s ramp load (~4*N*) at rates between 1 and 4 fps for up to 2000 s. The frame rate was set on the Antares system by a waveform generator that triggered the ECG module on the scanner. Because scanner memory was limited, there was a few-seconds gap between groups of 60 to 100 RF frames as memory contents were downloaded to disk. For breast scanning, the frame rate was fixed at 13 fps for a 10 to 20 s acquisition time. The frame rate for RF echo simulations (described below) was adjusted between 0.01 and 10 fps to properly sample the creep curves. Sampling issues and deconvolution of the strain series to estimate compliance spectra were detailed in [21].

Ultrasonic data acquisition is illustrated. From a time series of RF echo frames, a strain image sequence is formed and processed to generate KVFD model parameter images, including the *α* image. The initial and final radial frequencies in the analysis **...**

A time series of strain images was formed from the RF echo frames off line by applying a multi-resolution cross-correlation strain algorithm [11]. For each pixel in the spatially registered strain image series pictured in Fig. 2, we plotted the computed strain values versus time to obtain a creep curve. As described previously [24], we applied a preprocessing step to the strain data that eliminated purely elastic and purely viscous contributions, thus isolating the viscoelastic creep response, e.g., (17). Applying (16) and (20) to creep curves, we used Fourier methods to estimate the viscoelastic parameters in the frequency domain for all data except the *in vivo* breast data which were processed in the time domain. The frequency range in the loss spectrum plotted on log-log axes that exhibits a linear response was selected for analysis.

As stated earlier, a numerical linear fit was performed in the log-log domain over a user-defined frequency range on the loss modulus. This fit provided an estimate for the slope *α* and the intercept ln(*τ′ ^{α}E*

In all cases, parametric noise was reduced by averaging creep curves from a local neighborhood of strain pixels before parameter estimation. Spatial averaging thus reduces image noise at the cost of parametric spatial resolution. An example of an *α* image using a 27 × 49 pixels spatial averaging window is shown in Fig. 2, top right.

We simulated a time series of ultrasonic RF echo signal frames from uniformly scattering media with spatially varying viscoelastic properties. The goal was to discover how spatial variations in material properties influence spatial variations in estimates of *E*_{0}, *τ′*, and *α*.

Two-dimensional random scattering fields were generated in Matlab (The MathWorks, Inc, Natick, MA) to facilitate simulation of RF echo frames with fully developed speckle such as those from uniform graphite–gelatin phantoms [17]. The scattering field was numerically generated so that the coordinate of each scatterer was preserved as a floating-point value. Corresponding to the scattering field was a meshed map of viscoelastic material properties that was used in a finite element analysis (FEA). Displacements from the FEA output were used to reposition each scatterer in the scattering field to simulate deformation. The result was convolved with an ultrasonic pulse-echo impulse response typical of the Siemens Antares system data. Details of the RF echo simulations are provided elsewhere [39].

ABAQUS commercial FEA software (Simulia, Providence, RI) was used to simulate deformations of the scatterer field. The software allows the modeling of a biphasic poroviscoelastic medium with spatial heterogeneities in which the strain response to a load is determined by a viscoelastic solid matrix and fluid flow through the porous solid matrix. The viscoelastic nature of the solid matrix phase is specified using the relaxation constants of the discrete shear relaxation spectrum *g*_{1} and *g*_{2} and associated time constants *τ*_{1} and *τ*_{2}. The normalized shear modulus of the matrix is modeled as

$$g(t)=\frac{G(t)}{G(0)}=1-{g}_{1}(1-{e}^{-t/{\tau}_{1}})-{g}_{2}(1-{e}^{-t/{\tau}_{2}}).$$

Other FEA model parameters are the elastic modulus, *E _{m}*, and Poisson's ratio,

We simulated echo data from the central axial plane of the gelatin phantom with a circular inclusion assuming plane strain deformation. The phantom used in the FEA was modeled using an 8-node plane strain porous element (CPE8P) in ABAQUS. A uniform compressive stress of 100 Pa was applied to the top surface; the side surfaces were unbounded and the bottom surface permitted no vertical displacement. The acoustic properties of the medium were statistically homogeneous. The FEA parameter values of the background and inclusion for 2 simulation studies are given in Table I. Derivatives of the interpolated FEA displacements provided us with ideal time series of strain images, referred to here as object strain frames. RF echo data were generated by displacing scatterers according to mesh parameters, generating RF echoes, and processing the results with the MRCC algorithm to obtain ultrasonic strain frame sequences.

A 5-cm cubic graphite-in-gelatin phantom block containing a stiff cylindrical inclusion was imaged in this study [24]. This phantom contained 5.5% w/w type-A animal-hide gelatin. Embedded within the block was an 8-mm-diameter inclusion with a concentration of 8% gelatin. Throughout the phantom, graphite powder (3% w/w) was added to generate ultrasonically tissue-like scattering and absorption. Also 0.1% by volume formaldehyde was added as a chemical cross linker. The balance of the phantom material was distilled water. Although the inclusion provided little ultrasonic contrast, it generated contrast in viscoelastic features. The higher gelatin concentration of the inclusion created a stiffer, more solid polymer, that relaxed slower than its background.

Two additional 5-cm cubic graphite-in-gelatin blocks were constructed. Both were made of 8% type-B gelatin throughout. However, during gelation, a linear track of acidic fluid was injected in one block and a basic fluid track in the other. Acids or bases applied during gelation have been shown to weaken the polymer structure by causing a charge imbalance during polymerization, which respectively causes a softening or stiffening of the material. With these phantoms, we expect to see little ultrasonic contrast and significant viscoelastic contrast from the fluid injection sites. The goal is to use *α* images to see how pH affects the liquid-solid components of the hydrogel. pH effects may play an important role in the elasticity image contrast of malignant tissues.

Gelatin gels as we have prepared them are not intended for use as tissue-mimicking materials. Although both are multiphasic polymers, the solid matrix of gelatin gels has the aggregate structure of denatured collagen instead of the triple helical structure of natured type I collagen that constitutes the ECM of breast stroma. Gelatin gels also lack the ground substance found in all connective tissues of the body. As we discussed previously [17], the mechanical behavior of the 2 media are roughly similar to encourage the use of gelatin gels for developing ultrasonic elasticity imaging techniques. The most useful feature of gelatin is in creating phantoms for imaging viscoelastic parametric contrast. However comparisons of the 2 media show measurable differences in details of the dynamics [21].

Two different software phantom simulations are conducted; in each case a stiff circular inclusion is positioned in a soft background material. Viscoelastic model parameters for these 2 simulations are listed in Table I. From the FEA presented in Section III-A, we generate displacement maps for each phantom as a function of time. FEA displacement fields are processed to generate the object strains and object parametric images in Fig. 3. These maps are also used to model the displacements of random scattering fields, and RF waveforms are generated to simulate echo signals from a linear array. The ultrasonic parametric images obtained from these RF echo simulations are shown in Fig. 4. The average KVFD parameters estimated from the images presented in Figs. 3 and and44 are listed in Table II.

Notice from Tables I and andIIII that, though the matrix moduli *E _{m}* and strains observed for both simulations are similar, most FEA input parameters vary. In the first simulation, the viscoelastic properties of the inclusion and background, modeled by the fluid permeability

Figs. 3 and and44 show that hydropolymers may be differentiated based on their viscoelastic properties through the study of the KVFD parameter images, which differentiate between fluidic and solid-matrix responses of polymers.

A gelatin block with a high-concentration cylindrical inclusion was imaged during a 1460 s creep acquistion. Fig. 5 shows that the inclusion strains less. Hence, a significantly higher elastic modulus *E*_{0} is observed (due to the inverse proportional relation between *E*_{0} and the amplitude of the creep). The lower *α* parameter values in the inclusion indicate a greater solid matrix response. Though this phantom has different viscoelastic properties, the trends of the KVFD features (Table III) observed here are consistent with those observed in the ultrasonic parametric images for simulation (1) of Section IV–A. We previously noticed that a processing artifact was responsible for the contrast observed in the *τ′* ultrasonic parametric image. Here, the *τ′* image presents a similar contrast which may therefore be due to this artifact or may suggest a slow fluidic motion in the inclusion introduced by its higher concentration of gelatin.

KVFD parameter images for an experimental gelatin phantom with a stiff inclusion of diameter 8 mm. The blurring effect is due to the expanded size of the spatially averaging window needed to decrease strain noise. The values observed for *E*_{0} show that **...**

Two gelatin phantoms (isoelectric pH 5.6) were injected, respectively, with acidic and basic solutions during polymerization according to the technique presented in Section III-B-2. The strain image sequences obtained for these phantoms with inclusions polymerized at pH 4.6 [Fig. 6(a)] and pH 6.6 [Fig. 6(b)] show that an acidic fluid injection forms a soft inclusion region while a basic fluid injection locally stiffens the gel. This information is corroborated by the study of the average KVFD parameter values obtained in regions inside and outside the inclusion of these 2 phantoms. As would be expected from the strain images, *E*_{0} is lower in the acidic inclusion and higher in the basic one. The *α* study confirms the matrix changes following the inclusion. Indeed, the average value of *α* is lower in basic than in acidic inclusions, which indicates that basic gels are more elastic. The values of *τ′* suggest that water moves more rapidly in the basic inclusion than in the acidic inclusion, which can be explained by differences in mobile charge density [40].

Breast data were obtained from 3 normal female volunteers between the ages of 23 and 28 years [21]. Five patients between the ages of 42 and 85 years with biopsy-verified breast lesions were also examined [41]. Two to 4 acquisitions were recorded for each subject on the Antares scanner, as described in Section III. Patients were positioned supine and the breast was scanned AP with the chest wall as compression support. Conversely, volunteer subjects were positioned to lie on either side. The bicep of the lower arm supported the adjacent breast from below while data were acquired from above by contact scanning with medio-lateral positioning. Normal volunteers were instructed to breathe with shallow diaphragmatic movements to minimize breast motion during data acquisition that lasted up to 200 s. Breast patients were instructed to hold their breath during echo frame acquisitions up to 20 s. For both patients and volunteers, the hand-held transducer assembly was used to suddenly apply and hold constant a small compressive force (<5N) to the breast surface during acquisition. Motion from the heart beat was ignored because its signal fell outside the measurement bandwidth.

Spatially averaged creep curves were fit in the time domain to (17), as shown in Fig. 7. KVFD model parameters (Table IV) show that *α* tends to be higher in normal breast tissue than in gelatin, which is consistent with previous observations [21] that breast tissues have a less elastic and more fluidic response.

Clinical ultrasonic data were previously obtained for breast cancer patients with nonpalpable and biopsy-confirmed benign and malignant tumors [41]. These data were fit in the time domain to the KVFD model, (17). Table V gives the average parameter values inside and outside the lesion for 2 patients with benign lesions and 3 patients with malignant tumors.

Average KVFD Model Parameter Values (*α*, *E*_{0}, and *τ′*) Obtained from *in Vivo* Measurements in the Lesion and Background of 2 Benign (1–2) and 3 Malignant (3–5) Breast Tumors from Breast Cancer Patients.

Because the data sets were acquired over relatively short times, the fit is quite imprecise. The following are hypothesis formulated based on the observed results. Table V shows that *τ′* values are generally larger in lesions than in the background, and that *E*_{0} values in the benign lesions are much larger than those in the corresponding background regions, while lesion and background values for *E*_{0} tend to be more similar in malignant tumors. This observation is consistent with the histological view of benign lesions having a solid, elastic nature because of their dense concentration of normal collagen. Nonpalpable malignant lesions, on the other hand, retain a more fluidic viscoelastic response even with increases in collagen density because of differences in the ECM proteins. These observations are contradicted by the values of *α*. For all parameters, more data sets with longer acquisition times are needed to draw any conclusions about diagnostic value.

While we have assumed that the medium is loaded uniaxially and by a step function in time, neither assumption is strictly valid. The high geometric symmetry of gelatin blocks, where the scan plane bisects the cube and any inclusions, the freely slipping surfaces, the unrestricted lateral boundaries, and the machine-applied load all lend confidence to the assumption of uniaxial compression in phantoms. Conversely, *in vivo* breast tissues are geometrically asymmetric, highly heterogeneous in mechanical properties, and perhaps even anisotropic. We are certainly less confident in the uniaxial-loading assumption *in vivo*, and have yet to explore the consequences of violations on parameter estimates. The assumption of step loading has been studied to develop data processing guidelines and correction factors where necessary [21].

In this paper, we proposed applications of fractional derivatives to rheological modeling to reduce the feature space for materials characterization and thus provide parameters for imaging. After mathematically justifying the use of fractional derivatives, we imaged and analyzed the 3 parameters of the KVFD model for simulated and experimental hydropolymers as well as *in vivo* breast tissue data.

We showed that regions with significantly different viscoelastic responses can be detected with this method and may be differentiated provided the data acquisition time is long enough and strain noise is limited. Though we were not able to consistently differentiate nonpalpable breast tumors in preliminary patient data, we identified typical parameter values and were able to show, based solely on the 3 parameters of the model, that breast tissues are significantly more fluidic in responding to mechanical forces than gelatin phantoms. We were also able to differentiate an acidic region in gel from a basic region, and thus consider a similar technique, with improved image quality, to be very promising for detecting hypoxic tissues.

The slope parameter *α* is uncorrelated with other rheological model parameters and is largely unaffected by spatially averaging. Provided it can be precisely estimated over the available acquisition time, future work may focus on the analysis of *α* images.

The authors gratefully acknowledge Yusuke Kawada for his assistance. Mallika Sridhar and Karen Lindfors are acknowledged for their parts in the patient studies at University of California, Davis Medical Center.

This report is based upon work supported by the NIH/NCI under Award No. R01CA082497.

**Cecile Coussot** obtained her B.S. degree in engineering at Supelec, Gif-sur-Yvette, France, in 2006, and is currently working toward her M.S. degree in the Department of Bioengineering, University of Illinois at Urbana-Champaign, Urbana, IL. She is the recipient of one of the 2007 Roy S. Carver fellowships. Her research interests include modeling and analyzing polymer and breast tissue viscoelasticity and developing new algorithms for the processing of ultrasonic data.

**Sureshkumar Kalyanam** received a Ph.D. degree in aeronautics and astronautics from Purdue University, West Lafayette, IN, in 2004 for his research work in the area of domain switching and its effects on the fracture behavior of piezoelectric materials. He was a postdoctoral researcher in the Computational Fracture Mechanics Group, Department of Civil Engineering, University of Illinois at Urbana-Champaign (UIUC) between 2004 and 2006, where he studied the effects of delamination cracks in advanced Al-Li alloys using 3D FEA. He is currently a visiting research assistant professor in the Bioengineering Department, UIUC. His research focuses on using poroviscoelastic FEA to study the influence of solid matrix deformation and fluid movement on the response of gelatin hydrogels and the development of indentation testing combined with numerical methodology for characterizing soft biological materials.

**Rebecca Yapp** received her B.S. degree in physics from Illinois State University, Normal, IL, in 2004, and her M.S. degree from the University of Illinois at Urbana-Champaign in 2008. She is currently a Ph.D. candidate in bioengineering at UIUC. Rebecca's research interests include determining contrast parameters for breast cancer viscoelasticity imaging and understanding the responsible mechanisms.

**Michael F. Insana** (M'85) received his Ph.D. degree in medical physics from the University of Wisconsin–Madison in 1983. His dissertation topic was acoustic scattering and absorption in biological tissues. He was a research physicist at the FDA from 1984 to 1987, where he developed methods for describing tissue microstructure from the statistical properties of echo signals. From 1987 to 1999, he was with the Department of Radiology at the University of Kansas Medical Center, Kansas City, KS. There he directed an imaging research lab that applied ultrasonic imaging to the evaluation of progressive renal failure and breast cancer. Between 1999 and 2004, Mike was Professor of Biomedical Engineering at the University of California, Davis, where he directed the graduate program. He is currently Professor of Bioengineering and ECE at the University of Illinois at Urbana-Champaign. His research includes the development of ultrasonic instrumentation and methods for imaging soft tissue viscoelasticity and blood flow. The goal is to understand basic mechanisms of tumor formation and responses to therapy. He is also interested in the principles of imaging system design and performance evaluation, and signal processing, detection, and estimation. He is a member of the Acoustical Society of America, Senior Member of the IEEE, Fellow of the Institute of Physics, Fellow of the American Institute of Medical and Biological Engineering, and Associate Editor for *IEEE Transactions on Medical Imaging*.

Derivation of (11):

$$\begin{array}{c}\sigma (t)=\frac{{E}_{1}{\tau}^{\alpha}}{\Gamma (1-\alpha )}{\int}_{0}^{t}d{t}^{\prime}{(t-{t}^{\prime})}^{-\alpha}\phantom{\rule{0.2em}{0ex}}\frac{d\epsilon ({t}^{\prime})}{d{t}^{\prime}}\hfill \\ \phantom{\rule{2em}{0ex}}=\frac{{E}_{1}{\tau}^{\alpha}}{\Gamma (\gamma )}{\int}_{0}^{t}d{t}^{\prime}{(t-{t}^{\prime})}^{\gamma -1}\frac{d\epsilon ({t}^{\prime})}{d{t}^{\prime}}\text{for}\phantom{\rule{0.3em}{0ex}}\gamma =1-\alpha \hfill \\ \phantom{\rule{2em}{0ex}}={E}_{1}{\tau}^{\alpha}{\phantom{\rule{0.2em}{0ex}}}_{0}{\mathcal{D}}_{t}^{-\gamma}\left(\frac{d\epsilon (t)}{dt}\right)\phantom{\rule{0.3em}{0ex}}\text{by definition of}\phantom{\rule{0.3em}{0ex}}{{}_{a}\mathcal{D}}_{t}^{-\gamma}(f)\hfill \\ \phantom{\rule{2em}{0ex}}={E}_{1}{\tau}^{\alpha}{{}_{0}\mathcal{D}}_{t}^{\alpha -1}\left(\frac{d\epsilon (t)}{dt}\right)\hfill \\ \phantom{\rule{2em}{0ex}}={E}_{1}{\tau}^{\alpha}\frac{{d}^{\alpha -1}}{d{t}^{\alpha -1}}\left(\frac{d\epsilon (t)}{dt}\right)\hfill \\ \sigma (t)={E}_{1}{\tau}^{\alpha}\frac{{d}^{\alpha}\epsilon}{d{t}^{\alpha}}.\hfill \end{array}$$

Definition of the general Mittag-Leffler function (17):

$${E}_{\alpha ,\beta}(z)=\sum _{k=0}^{\infty}\frac{{z}^{k}}{\Gamma (\alpha k+\beta )}\phantom{\rule{1.2em}{0ex}}\text{for}\phantom{\rule{0.5em}{0ex}}\alpha ,\beta >0.$$

1. Kossoff MB. Ultrasound of the breast. World J Surg. 2000;24:143–157. [PubMed]

2. Marini C, Traino C, Cilotti A, Roncella M, Campori G, Bartolozzi C. Differentiation of benign and malignant breast microcalcifications: Mammography versus mammography-sonography combination. Radiol Med (Torino) 2003;105:17–26. [PubMed]

3. Garra BS, Cespedes EI, Ophir J, Spratt SR, Zuurbier RA, Magnant CM, Pennanen MF. Elastography of breast lesions: Initial clinical results. Radiology. 1997:79–86. vol. [PubMed]

4. Thomas A, Kümmel S, Fritzsche F, Warm M, Ebert B, Hamm B, Fischer T. Real-time sonoelastography performed in addition to B-mode ultrasound and mammography: Improved differentiation of breast lesions? Acad Radiol. 2006;13:1496–1504. [PubMed]

5. Thomas A, Warm M, Hoopmann M, Diekmann F, Fischer T. Tissue Doppler and strain imaging for evaluating tissue elasticity of breast lesions. Acad Radiol. 2007;14:522–529. [PubMed]

6. Zhi H, Ou B, Luo BM, Feng X, Wen YL, Yang HY. Comparison of ultrasound elastography, mammography, and sonography in the diagnosis of solid breast lesions. J Ultrasound Med. 2007;26:807–815. [PubMed]

7. Elenbaas B, Weinberg R. Heterotypic signaling between epithelial tumor cells and fibroblasts in carcinoma formation. Exp Cell Res. 2001;264:169–184. [PubMed]

8. Discher DE, Janmey P, Wang YL. Tissue cells feel and respond to the stiffness of their substrate. Science. 2005;310:1139–1143. [PubMed]

9. Suresh S. Biomechanics and biophysics of cancer cells. Acta Biomater. 2007;3:413–438. [PMC free article] [PubMed]

10. Greenleaf JF, Fatemi M, Insana M. Selected methods for imaging elastic properties of biological tissues. In: Yarmush Martin., editor. Annual Review of Biomedical Engineering. Vol. 5. Palo Alto: Annual Reviews; 2003. pp. 57–78. [PubMed]

11. Pellot–Barakat C, Frouin F, Insana MF, Herment A. Ultrasound elastography based on multi-scale estimations of displacement regularized fields. IEEE Trans Med Imag. 2004;23:153–163. [PMC free article] [PubMed]

12. Viola F, Walker WF. A spline-based algorithm for continuous time-delay estimation using sampled data. IEEE Trans Ultrason Ferroelectr Freq Control. 2005;52:80–93. [PubMed]

13. Pinton GF, Trahey GE. Continuous delay estimation with polynomial splines. IEEE Trans Ultrason Ferroelectr Freq Control. 2006;53:2026–2035. [PubMed]

14. Nightingale KR, Palmeri ML, Nightingale RW, Trahey GE. On the feasibility of remote palpation using acoustic radiation force. J Acoust Soc Am. 2001;110:625–634. [PubMed]

15. Barannik EA, Girnyk A, Tovstiak V, Marusenko AI, Emilianov SY, Sarvazyan AP. Doppler ultrasound detection of shear waves remotely induced in tissue phantoms and tissue in vitro. Ultrasonics. 2002;40:849–852. [PubMed]

16. Sinkus R, Siegmann K, Xydeas T, Tanter M, Claussen C, Fink M. MR elastography of breast lesions: Understanding the solid/liquid duality can improve the specificity of contrast-enhanced MR mammography. Magn Reson Med. 2007;58:1135–1144. [PubMed]

17. Sridhar M, Liu J, Insana MF. Elasticity imaging of polymeric media. ASME J Biomech Eng. 2007;129:259–272. [PMC free article] [PubMed]

18. Mak AF. The apparent viscoelastic behavior of articular cartilage: The contributions from the intrinsic matrix viscoelasticity and interstitial fluid flows. ASME J Biomech Eng. 1986;108:123–130. [PubMed]

19. Mak AF. Unconfined compression of hydrated viscoelastic tissues: A biphasic poroviscoelastic model. Biorheology. 1986;23:371–383. [PubMed]

20. DiSilvestro MR, Suh JK. Biphasic poroviscoelastic behavior of hydrated biological soft tissue. ASME J Appl Mech. 1999;66:528–535.

21. Sridhar M, Insana MF. Ultrasonic measurements of breast viscoelasticity. Med Phys. 2007;34:4757–4767. [PMC free article] [PubMed]

22. Muir H. Proteoglycans as organizers of the intercellular matrix. Biochem Soc Trans. 1983;11:613–622. [PubMed]

23. Ferry JD. Viscoelastic Properties of Polymers. New York: Wiley; 1980.

24. Sridhar M, Liu J, Insana MF. Viscoelasticity imaging using ultrasound: Parameters and error analysis. Phys Med Biol. 2007;52:2425–2443. [PMC free article] [PubMed]

25. Sloninsky GL. Laws of mechanical relaxation processes in polymers. J Polym Sci, Part C. 1967;16:1667–1672.

26. Bagley RL, Torvik PJ. A theoretical basis for the application of fractional calculus to viscoelasticity. J Rheol. 1983;27:201–210.

27. Schiessel H, Blumen A. Hierarchical analogues to fractional relaxation equations. J Phys Math Gen. 1993;26:5057–5069.

28. Song DY, Jiang TQ. Study on the constitutive equation with fractional derivative for the viscoelastic fluids—Modified Jeffreys model and its application. Rheol Acta. 1998;37:512–517.

29. Djordjevic VD, Jaric J, Fabry B, Fredberg JJ, Stamenovic D. Fractional derivatives embody essential features of cell rheologi-cal behavior. Ann Biomed Eng. 2003;31:692–699. [PubMed]

30. Heymans N. Fractional calculus description of non-linear viscoelastic behaviour of polymers. Nonlinear Dynam. 2004;38:221–231.

31. Liu H, Oliphant TE, Taylor L. General fractional derivative viscoelastic models applied to vibration elastography. Proc IEEE Ultrason Symp. 2003:933–936.

32. Chatterjee A. Statistical origins of fractional derivatives in viscoelasticity. J Sound Vibrat. 2005;284:1239–1245.

33. Pfitzenreiter T. A physical basis for fractional derivatives in constitutive equations. Z Angew Math Mech. 2004;84:284–287.

34. Kawada Y, Nagahama H, Hara H. Irreversible thermodynamic and viscoelastic model for power-law relaxation and attenuation of rocks. Tectonophysics. 2006;427:255–263.

35. Tschoegl NW. The Phenomenological Theory of Linear Viscoelastic Behavior. Berlin: Springer-Verlag; 1989.

36. Kawada Y, Nagahama H. Viscoelastic behaviour and temporal fractal properties of lherzolite and marble: Possible extrapolation from experimental results to the geological time-scale. Terra Nova. 2004;16:128–132.

37. Schiessel H, Metzler R, Blumen A, Nonnenmacher TF. Generalized viscoelastic models: Their fractional equations with solutions. J Phys Math Gen. 1995;28:6567–6584.

38. Mittag–Leffler MG. Sur la nouvelle fonction *E*_{α}(*x*) Comptes Rendus Acad Sci Paris. 1903;137:554–558.

39. Liu J, Abbey CK, Insana MF. Linear approach to axial resolution in elasticity imaging. IEEE Trans Ultrason Ferroelectr Freq Control. 2004;51:716–725. [PMC free article] [PubMed]

40. Yapp RD, Insana MF. pH-induced contrast in viscoelasticity imaging of biopolymers. Phys Med Biol. 2009;54:1089–1109. [PMC free article] [PubMed]

41. Qiu Y, Sridhar M, Tsou JK, Lindfors KK, Insana MF. Ultrasonic viscoelasticity imaging of nonpalpable breast tumors: Preliminary results. Acad Radiol. 2008;15:1526–1533. [PMC free article] [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |