Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
IEEE Trans Ultrason Ferroelectr Freq Control. Author manuscript; available in PMC 2010 April 1.
Published in final edited form as:
PMCID: PMC2704493

Fractional Derivative Models for Ultrasonic Characterization of Polymer and Breast Tissue Viscoelasticity


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.

I. Introduction

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.

II. Modeling Relaxation and Creep

A. Power-Law Model of Relaxation

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(t), strain εn(t), and strain rate ε˙n(t) to the elastic modulus En(t), each as a function of time, t [35],


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 τn, via


where E′n is the modulus and τ n = η n/E′n is the characteristic relaxation time related to the viscous coefficient ηn for unit n.

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


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


Combining (2)–(4) we obtain


where F(n,t) = n ln bant/τ0.

It is observed with polymers that the time dependence of the elastic modulus E(t) is a power law. Let n0 indicate the viscoelastic unit where e F(n0,t) is maximum and let us assume that the material constants E′n and τn from sequential units vary slowly and monotonically. Consequently, the unit responses are highly correlated, and F(n,t) is well represented by a second-order Taylor series expansion about n0.


since ([partial differential]F/[partial differential]n)n0 = 0. Combining the definition of F(n,t) in (5) with its first derivative being zero, we find that


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 n; therefore (4) gives a > 0. E′n must be positive and remain finite; consequently, 0 < b < 1. In order for e F(n0,t) to exist and be a maximum for all t, E′n and τn must vary in opposite ways. For example, experiments show that E(t) decreases with time. If E′n decrease with n, then τn must increase with n to obtain a power-law form. From this relation, we can derive that 0 < a < 1. These bounds on model constants have been experimentally verified [36]. The power-law assumption also does not permit a to approach 1; that is, from (5),


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,


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


where E1 is the scaled modulus for unit n0, τ is a characteristic relaxation time constant, and Γ is the Gamma function for non-integer arguments.

B. Fractional Derivatives

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


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


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


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


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 E1, 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.

C. KVFD Models for Relaxation and Creep

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


and consisting of elastic and viscoelastic terms. E0 and E1 may have different values; however, τ can be modified so that


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)


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


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 α from 0.2 to 0.8 in steps of 0.2. Lower values of α induce quick responses of low amplitude while higher values are responsible for slower reactions exhibiting greater net strain. These 2 types of behavior respectively describe more elastic and more viscous materials. The more complex expression for compliance in the time domain, as compared with relaxation, is a consequence of the following relation between D(t) and E(t) [35]:

Fig. 1
Modeled strain versus time from (17) shown for different values of α and constant E0 = 1 Pa and τ′ = 10 s. Since a unit step stress is assumed, D(t) = ε(t).

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


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


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.

III. Data Acquisition and Processing

The equations above show that the viscoelastic parameters E0, τ′, 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 (~4N) 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].

Fig. 2
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(τ′αE0sin(απ/2)). The remaining necessary information to approximate the 3 unknowns was extracted from the time-domain strain curve using ε/σ0 ≈ 1/E0. If the acquisition time is long enough, this value may be reliably estimated, which once again justifies the need for long acquisition times when processing data in the Fourier domain. The acquisition time of the in vivo breast data was too short (20 s) to apply Fourier techniques; instead, strain curves were directly fitted to (17) to estimate the model parameters in the time domain.

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.

A. RF Echo Simulations with Deformation

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 E0, τ′, 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 g1 and g2 and associated time constants τ1 and τ2. The normalized shear modulus of the matrix is modeled as


Other FEA model parameters are the elastic modulus, Em, and Poisson's ratio, νm, of the porous solid matrix phase. Finally, the fluid-flow-dependent viscoelastic behavior is specified by the hydraulic permeability, k, of the fluid movement through the porous matrix.

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.

Table I
FEA Simulation Parameters for the Two Simulations*

B. Gelatin Imaging Phantoms

1) Contrast from Gelatin Concentration

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.

2) Contrast from pH

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.

3) Interpreting Gelatin Phantom Results

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].

IV. Results

A. Phantom Simulations

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.

Fig. 3
KVFD model images obtained for the strain image sequences generated from derivation of the interpolated FEA displacements: Object parametric images. Two simulated phantoms with VE properties given in Table I are compared using the KVFD approach. Elastic ...
Fig. 4
KVFD model images obtained for the strain image sequences generated from simulated echo data of 2 deformed FEA phantoms: Ultrasonic parametric images. These images give results similar to those observed in Fig. 3. The parameter images show significant ...
Table II
Average KVFD Model Parameter Values (α, E0, and τ′) Estimated in the Inclusion and Background for the Object (top) and Ultrasonic (bottom) Parametric Images Displayed in Figs. 3 and and44.

Notice from Tables I and andIIII that, though the matrix moduli Em 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 k, Poisson's ratio νm, and time constants τ1 and τ2, differ. In the second simulation, however, these parameters are similar inside and outside the inclusion. The α and τ′ images in Fig. 3 consequently present a very high contrast in the first simulation and very little in the second. Although contrast polarity is preserved in the α and E0 ultrasonic parametric images of Fig. 4, contrast magnitude is reduced relative to that of the object parametric images of Fig. 3 because of the spatial filtering applied to reduce estimation noise. The inverted contrast observed in the τ′ ultrasonic image is due to the correlation between the values of E0 and τ′ given by the Fourier domain fit. The large averaging window used here causes an apparent decrease of the creep amplitude in the inclusion, hence a decrease of E0 and a compensating increase of τ′. The value of α on the other hand is uncorrelated to the other parameter values and appears less affected by the averaging window.

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.

B. Gelatin Phantoms

1) Contrast from Gelatin Concentration

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 E0 is observed (due to the inverse proportional relation between E0 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.

Fig. 5
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 E0 show that ...
Table III
Average KVFD Model Parameter Values (α, E0, and τ′) in the Inclusion and Background for the Experimental Phantom Presented in Section IV-B-1.

2) Contrast from pH

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, E0 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].

Fig. 6
Typical strain images obtained for 2 type-B gels respectively injected in the inclusion with an acidic and basic solution during polymerization, and corresponding average KVFD parameter values in 2 regions. In (a), a strain frame of the acidic inclusion ...

C. In Vivo Breast Tissue Data Study

1) Breast Data Acquisition

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.

2) Average Gelatin Phantom and Breast Tissue Responses

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.

Fig. 7
Average creep response observed in an experimental phantom (a) and in normal breast tissue measured in vivo (b). Data were fit in the time domain according to (17). The corresponding KVFD model parameters are given in Table IV.
Table IV
KVFD Approach Parameter Values (α, E0, and τ′) Obtained from Gelatin Phantom and Normal Glandular Breast Tissue Responses.

3) Cancer Patient Data

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 (α, E0, 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 E0 values in the benign lesions are much larger than those in the corresponding background regions, while lesion and background values for E0 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.

4) Assumptions

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].

V. Conclusion

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.


An external file that holds a picture, illustration, etc.
Object name is nihms122020b1.gif

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.

An external file that holds a picture, illustration, etc.
Object name is nihms122020b2.gif

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.

An external file that holds a picture, illustration, etc.
Object name is nihms122020b3.gif

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.

An external file that holds a picture, illustration, etc.
Object name is nihms122020b4.gif

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):

σ(t)=E1ταΓ(1α)0tdt(tt)αdε(t)dt=E1ταΓ(γ)0tdt(tt)γ1dε(t)dtforγ=1α=E1τα0Dtγ(dε(t)dt)by definition ofDatγ(f)=E1ταD0tα1(dε(t)dt)=E1ταdα1dtα1(dε(t)dt)σ(t)=E1ταdαεdtα.

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



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]