Model-based elastography: a survey of approaches to the inverse elasticity problem

3.2.1 Weak form

The weak form is a weighted-integral statement that is equivalent to the governing partial differential equation and the natural boundary conditions of the problem. The weak form of the governing elasticity equation for each approach to elastographic imaging can be derived as described in (Reddy, 1993). For simplicity, we will describe only the finite element implementation for the harmonic case (i.e., equation (8)), since the procedure is identical for the quasi-static and transient cases. The weak form of equation (8) with a scalar weighting function ϕi(x, y, z) is given by:

where ${\oint }_{\Gamma }$ represents an integration over the boundary ( Γ) of the element, ϕ_i is a scalar basis associated Γ with the element, n̂ represents the outward-pointing normal vector, and t̂ represents the Nuemann boundary condition. If the scalar components of the displacement vector, u, are u(x, y, z),v,(x, y, z),w(x, y, z) in the x,y, and z directions, respectively, then the Galerkin approximation of the displacement may be derived by expanding the scalar components in the ϕ basis to give:

where N is the number of nodes associated with each element in the finite element mesh and u_j ,v_j , w_j are the axial, lateral, and elevational displacement components at each element node. The Galerkin weak-form finite element model is obtained by substituting equation (18) into equation (17):

where f is the 3N × 1 global force vector and K(μ,λ,ρ ) that is given by

where the coefficients of the global stiffness matrix are determined by the basis function and the material properties (i.e., μ and λ ).

3.2.2 Dimensionality reduction using the plane-strain and plane-stress approximation

Despite the technological advances in 3D ultrasound imaging, ultrasound is predominately a two-dimensional imaging modality. Thus, most investigators in ultrasonic imaging typically reduce the 3D elasticity problem to a two-dimensional problem, using either a plane-strain or plane-stress approximation. To illustrate the differences between the two approximations, let’s consider a linear elastic 1 1 solid, Ω, of uniform thickness, h , bounded by two parallel planes ($z=-\frac{1}{2}h$ and $z=\frac{1}{2}h$) and a closed 2 2 boundary, Γ. If the thickness of h≫Ω, then the problem can be considered a plane-strain problem. In this case, we assume thatε₁₃ = ε₂₃ = ε₃₃ = 0. The plane-strain approximation is typically used when structures above and below the ultrasound scan plane have motion confined to the elevation (z) direction (Kallel et al., 2001; Doyley et al., 2000; Skovoroda et al., 1995; Barbone and Bamber, 2002). However, when h≪Ω, the problem is considered a plane-stress problem. In this case, we assume that σ₁₃ = σ₂₃ = σ₃₃ = 0). The plane-stress approximation generally applies to thin plates; however, (Sumi et al., 1995b) used a plane-stress approximation to reduce the 3D elasticity problem to 2D. In reality, the plane-strain and plane-stress assumptions are valid only for special cases such as phantoms with cylindrical inclusions, or when elastography is performed using the constrained imaging method described in (Kallel and Ophir, 1997). Therefore, errors are typically incurred when the 3D elasticity problem is modeled using either approximation.

4.1.1 Direct inversion

Raghavan and Yagle (1994) proposed a direct inversion scheme for recovering shear modulus. They derived a linear system of equations by re-arranging the PDEs that describe the direct problem for the plane-strain condition. The PDEs that (Raghavan and Yagle, 1994) derived are given by:

where ${\partial }_1=\frac{\partial }{\partial x}$, the unknowns are shear modulus (μ ) and hydrostatic pressure ( p)¹, and the coefficients are functions of the internal tissue strains that are related to the measured displacements (see equation (6)). The weakness of this approach is that both the shear modulus (μ ) and the hydrostatic pressure ( p ) on the boundary must be known to solve equation (21).

Although Rhagavan and Yagle’s computer simulation demonstrated that this inversion scheme could produce very encouraging displacements, they observed that the performance of the reconstructed elastograms degraded rapidly with increasing measurement noise. Another weakness of this inversion scheme is that no imaging system to date can measure hydrostatic pressure. To overcome this limitation, (Skovoroda et al., 1995) used an analytical method to eliminate the pressure term from equation (21). Eliminating the pressure term gives:

As in equation (21), the shear modulus must be known on the boundary of the region of interest (ROI) to solve equation (22). Skovoroda et al. (Skovoroda et al., 1995) were the first to demonstrate that this information could be obtained by exploiting the stress-continuity properties of soft tissues. They demonstrated that contours of small shear modulus variations could easily be defined, and that the modulus along the boundary could be computed by processing axial strain elastograms as follows:

where _μr is the relative shear modulus, and ${\epsilon }_{22}^{\mathit{int}}$ and ${\epsilon }_{22}^{\mathit{ext}}$ are the axial strains within and outside the boundary of ROI. Equation (22) contains high-order derivatives that will amplify measurement noise, which could compromise the quality of ensuing modulus elastograms as demonstrated in figure 4. Despite this potential limitation,, the viability of this inversion scheme has been demonstrated in gelatin and ex vivo kidney phantoms (Skovoroda et al., 1999; Skovoroda and Aglyamovv, 1995; Chenevert et al., 1998). Skovoroda et al. (1999) cast the inverse problem as an integral rather than a differential form — an approach that was first described in (Sumi et al., 1995b) — to make the technique less susceptible to measurement noise. To bolster performance further, they computed lateral displacement using the incompressible method described in (Lubinski et al., 1996).

Bishop et al. (2000) proposed to eliminate the pressure term appearing in equation (21) by partitioning the matrix, rather than by doing this analytically as described in (Skovoroda et al., 1995). But the resulting formulation proved to be ill-conditioned. Consequently, (Bishop et al. a , 2000) constrained the solution by employing the Tikhonov regularization method.

Sumi et al. (1995b) also proposed a direct inversion scheme; however, they solved the inverse problem for the plane-stress case, using the following PDEs:

where the unknowns are spatial derivatives of relative Young’s modulus, and the coefficients are strains and their spatial derivatives. They computed the shear modulus at a given point (x, y) within the tissue relative to a reference point —let’s say: (A, B) — by employing a line integral. The feasibility of this inversion scheme has been demonstrated in phantoms and excised tissues (Sumi, 2007; Sumi and Nakayama, 1998); nevertheless, the plane-stress condition is typically not relevant for most clinical applications. To make the technique more clinically relevant, (Le Floc’h et al., 2009) extended the concept to the plane-strain case. The equations they derived for solving the inverse problem for the plane-strain condition are given in compact form as follows:

Although theoretically feasible, equation (28) is difficult to solve, because $\frac{\nabla p}{E}$ cannot be measured in practice. However, (Le Floc’h et al., 2009) demonstrated that the second term could be used to highlight the boundaries of different tissue types.

4.1.2 Iterative Inversion

The inverse problem can also be viewed as a parameter-optimization problem, where the goal is to find the shear modulus that minimizes the error between measured displacement or strain fields, and those computed by solving the direct problem. This inversion approach has been successfully used in several emerging medical imaging modalities, as discussed in (Paulsen et al., 2005; Yorkey et al., 1987), and was applied to elastography by (Kallel and Bertrand, 1996; Doyley et al., 1996). The objective function that is minimized typically has the following form:

where u is the solution to the direct problem computed from the shear modulus distribution, μ , by solving equation (9). Minimizing equation (29) with respect to shear modulus variations is a nonlinear process; however, we can do this iteratively using a variety of techniques as discussed in (Vogel, 2002). Proposed optimization techniques can be classified into three broad categories: (1) those that require only the value of the functional π for different values of the parameters μ (i.e., non-gradient approaches); (2) those that require the value of the functional and its derivative with respect to the material parameters (called the gradient vector); and (3) those that require the first and second derivative of the functional with respect to the material parameters (called the Hessian matrix). All three optimization methods are illustrated below.

4.2.1 Direct inversion

Skinus et al. (2000) solved the inverse harmonic elastography problem using a linear system of PDE that they derived by solving the wave propagation model described in (Landau et al., 1986) in the frequency domain, and expressed Young’s modulus as a symmetric tensor. Lorenzen et al. (2003) used this inversion scheme to demonstrate that MRE can detect changes in breast tissue elasticity during the monthly hormonal cycle. Using the medical standard, the first day of menses is counted as day 1 in the woman’s cycle. Using MRE, (Lorenzen et al., 2003) showed that on day 5, the median value of elasticity for fibroglandular adipose tissue declined significantly, but at day 14, the same tissue’s elasticity increased noticeably.

Manduca et al. (2001) proposed to solve the inverse harmonic elastography problem using a direct inversion scheme, which they referred to as algebraic inversion of differential equation (AIDE). By assuming local homogeneity, the equations were solved separately at each pixel using only data from a local neighborhood to estimate local derivatives as described in (Oliphant et al., 2001). Very encouraging phantoms and patient results have been obtained using this technique. However, the large difference in magnitude of shear modulus and the second Lamé coefficient (i.e., kPa for shear modulus and GPa)prevents simultaneous estimation ofμ and λ , Manduca et al. (2001) assumed that the divergence of the displacement vector in equation (8) was negligible (i.e., ∇ · u ≈ 0). Using this assumption, μ was estimated from a single component of motion as follows:

Romano et al. (1998) developed a direct-inversion method by using the variational or weak-form of equation (8), and appropriately chosen test functions to estimate both Lamé constants. The advantage of using the weak form of equation (8) is that the derivative is calculated not from the measured data, but rather from a smooth test function.

Park and Maniatty (2006) also described a direct-inversion scheme for reconstructing shear modulus from displacements measured during time-harmonic excitation. Their approach assumed that the time-harmonic displacements were divergence-free, which reduces the governing equation to:

If the gradient of the shear modulus is also neglected, equation (38) reduces to the Helmholtz equation.

4.2.2 Iterative inversion method

Like its quasi-static counterpart, shear modulus can also be computed by treating the image-reconstruction problem as a parameter-optimization problem. However, to accosmplish this all component of the internal tissue displacement must be known. Using computer simulation, Van Houten et al. (2001) demonstrated that the displacement fields of an oscillating, 3D, isotropic, linear elastic body was not acurately characterized using 2D plane approximation. Consequently, no ultrasonic methods have been proposed to solve the inverse-harmonic-elastographic problem because without symmetries, the three dimensional case cannot be approximated.

The objective function to be minimized is identical to that used in the quasi-static elastography equation (29). However, the computational overhead required to solve the full 3D elasticity problem at the resolution MR data set with the Hessian method would make the computations infeasible on contemporary processors. For example, if a typical MR displacement data set were discretized to 16 × 256 × 256 image slices, the corresponding parameter set matching the MR resolution would have over a million elements (assuming a description based on single parameter of elasticity). Each parameter update for this large property description would require over ≈ 1e¹⁸ floating-point operations to invert the Hessian matrix, and an additional 1e¹⁸ operations to generate the matrix beforehand. Van Houten et al. (2001) demonstrated that this issue can be circumvented by dividing the reconstruction field-of-view into a series of overlapping subzones, and expanding equation (29) as a sum over all the subzones as follows:

where _uz (_μz ) represents the displacements on the zth subzone computed by solving the direct problem from the shear modulus (_μz ), and ${\mathit{u}}_z^m$ represents the corresponding MR measured tissue displacements. They assumed that minimizing the sum in equation (39) was equivalent to the sum of minimization of the individual subzones:

which involves equating derivatives of displacements with respect to subzone shear modulus to zero, and solving the resulting set of nonlinear equations with the Gauss-Newton method. The resulting matrix solution at iteration ( i+1 ) has the form:

where $\left[J\right({\mu }_z^i\left)\right]$represents the Jacobian matrix for a given subzone, and α_z is the regularization parameter on the subzone level. Van Houten et al. (2001) deployed the subzones in a random and overlapping manner from a diminishing list of seed points that corresponded to locations in the computational domain that had n not been previously updated; however the manner in which the subzones are deployed is not essential as demonstrated in (Doyley et al., 2007). When the shear modulus at all locations is updated, a global sweep is concluded and the process is repeated. Doyley et al. (2007) demonstrated that the computatiional load of this inversion technique increased linearly with increasing subzones — a limitation that they circumvented by reducing the size of the reconstruction field of view and the overlap between subzones. Doyley et al. (2005) also demonstrated that the spatial resolution of elastograms computed using this technique was on the order of 5 mm. Figure 7 shows a representative example of an elastogram obtained from a health volunteer using the subzone-inversion technique, which illustrates that the resolution of the elastograms was sufficiently high to visualize fibroglandular tissue from the adipose tissue. The advantage of this inversion scheme is that it is ideally suited for a parallel-computing platform, because the sub-domain discretizations are computationally independent as demonstrated in (Doyley et al., 2005).

Van Houten et al. (2003) reported results of a pre-clinical study that they performed on five healthy volunteers. Their results demonstrated that the elastic properties of the breast fibro-glandular and fatty breast tissues, measured in in vivo with the subzone inversion technique, were comparable to those reported in literature.

5.1.1 Quasi-static methods

Sridhar et al. (2007a, b) developed an elastographic imaging approach for characterizing viscoelastic materials. More specifically, they acquired time-varying axial strain elastograms when a viscoelastic material was subjected to a constant stress, and constructed a creep curve at each pixel from the time-varying strain elastograms. They computed the creep compliance from the the ratio of the time-varying strain, ε (t) , to the applied stress; and computed the complex compliance, D^*(s) , from the Laplace transform of the creep compliance:

where D’(s) and D” (s) represent the storage and loss compliance, respectively. Sridhar et al. (2007a) used a three-parameter Kelvin-Vogit rheological model to predict the strain response at each pixel within the time-varying strain elastogram. In the time domain, the resulting strain response is given by:

where ε_o is the instantaneous strain, i.e., the strain incurred immediately after compression; ε₁ is the viscoelastic strain amplitude; and t₁ is the retardation time, i.e., the time required for the tissue to become fully deformed. They estimated the three model pameters by fitting the rheological model to the measured creep response at each pixel. Using this tecnique, this group demonstated that t₁ is the most useful parameter for discriminating between malignant and benign breast tumors (Qiu et al., 2008). More specifically, t₁ values measured from malignant tissue were smaller than those measured from the surrounding healthy breast tissue; whereas, the converse was observed for benign tumors.

5.1.2 Harmonic methods

Hoyt et al. (2008) proposed a viscoelastic approach based on sonoelastography imaging. They measured the visoelastic properties of gelatin and muscle samples, by fitting a Voigt model to the dispersive shear speed obtained using the crawling-wave method described in (Wu et al., 2006). These authors validated their technique by conducting simulations, phantom studies, and human studies. The simulation studies revealed a 2.3 % difference in the computed and true shear speeds under ideal measurement conditions. The phantom studies revealed a 1 % error in shear modulus that had been computed using the proposed technique, relative to that measured using a mechanical testing system. They also observed a slight increase in shear speed with increased frequency. Statistically significant differences were observed in the shear and loss moduli of relaxed and volunteered contracted muscles. The authors also observed a slight increase in shear speed with frequency. To minimize any anisotropic effects, they obtained data parallel to the muscle fiber.

Asbach et al. (2008) used their multifrequency MR elastography method to measure the viscoelastic properties of normal liver tissue versus diseased liver tissue taken from patients with grade 3 and 4 liver fibrosis. Like Hoyt et al. (2008), they measured shear-wave and shear-attenuation dispersion; but, rather than using ultrasound, they computed shear wave and shear attenuation from the Fourier transform of the complex MR displacement fields. They computed the shear modulus and viscosity variations within the tissue by fitting a Maxwell rheological model to the measured data, by solving the linear viscoelastic wave equation in the frequency domain. They observed that fibrotic liver tissue had a higher viscosity (14.4 ± 6.6 Pa s) and elastic modulus (μ₁ = 2.91 ± 0.84kPa and μ₂ = 4.83 ± 1.77kPa) than normal liver tissue. Their results revealed that although liver tissue is dispersive, it appeared as non-dispersive between the frequency range of 25 — 50 Hz. This research group has also measured the viscoelastic properties of brain tissues (Hamhaber et al., 2010; Sack et al., 2009), but in that case, they characterized the viscous properties of the brain by fitting a Voigts model to measured complex modulus. They also observed that the viscosity didn’t agree with the predictions of the Kramers-Kronig relation (Klatt et al. , 2010; Madsen et al., 2008; Urban and Greenleaf, 2009). Klatt et al. (2010) also measured the viscoelastic properties of the liver. In this case, the spring-pot model was used to study the dispersive behavior of the viscoelastic properties between frequency ranges of 25, 37.5, 50 and 62.5 Hz. Like Hoyt et al. (2008), they observed that the stiffness and viscosity of muscle increased with voluntary contraction.

5.1.3 Transient methods

Catheline et al. (2004) were the first to propose a method to visualize the viscoelastic properties of soft tissues. Using transient elastography, they measured the spatial variation of the time-harmonic displacement field, and used a plane-wave approximation to compute the complex wave number ( k = k’ + iα_T ) of the propagating waves, as follows:

where u_z(x) is the measured displacement, and FT is the Fourier Transform. The wave speed and attenuation coefficients of the propagating shear wave relate to the complex wave number as follows:

where 𝕽 and 𝕵 represents the real and imaginary component of the complex wave number, respectively.

Catheline et al. (2004) used Agar-gelatin phantoms to demonstrate that the shear wave speed ( c₂ ) and attenuation coefficient (_αT ) computed using equation (45) were comparable to those measured independently from phase and amplitude measurements. They computed the shear modulus (μ ) and viscosity (η ) by fitting the measured speed of sound and attenuation to Voigt’s and Maxwell’s rheological models, whose wave speed and attenuation were given by:

and similarly, the attenuation was given by:

where the superscripts V and M represent the Voigt and Maxwell models, respectively. The recovered shear modulus values were independent of the rheological model employed, but researchers observed an order of magnitude difference in viscosity values recovered with the two models. This was not surprising, because the goodness of fit of the wave-speed dispersion data to both rheological models were comparable, but Catheline et al. (2004) had demonstrated that the Voigt model provided a better fit to the attenuation-dispersion data than the Maxwell model, which proved that attenuation coefficents were independent of frequency despite the strong frequency dependency that was apparent in the data. Catheline et al. (2004) also had conducted studies on excised bovine muscles, which revealed two apparent wave speeds: either a fast or a slow wave was observed, depending on the polarization. They obtained a good estimate of both wave speeds using equation (46), but there was an overestimation of attenuation, because the model assumed the displacement field arose solely from transverse waves.

Sinkus et al. (2005b) developed a direct-inversion scheme to visualize the mechanical properties of visocelastic materials, in which a curl operation was performed on the time-harmonic displacement field u(x, y,z,t) = u(x, y,z)e^jωt to remove the displacement contribution of the compressional wave. They dervived the governing equation that describes the motion incurred in an isotropic, viscoelastic medium by computing the curl of the PDEs that describe the motion incurred by both transverse and compressional shear waves. The resulting PDEs for transverse waves are given in compact form by:

Sinkus et al. (2005b) developed a direct-inversion scheme from equation (48), in which μ and η (viscosity) were the unknowns. They evaluated the inversion scheme using (a) computer simulations, (b) phantom studies, and (c) patient studies. Their simulation studies revealed that the proposed algorithm could accurately recover shear modulus and viscosity from ideal displacement data. However, with noisy displacements, a good estimate of shear modulus was obtained only when the shear modulus of the simulated tissue was < 8 kPa; the inversion scheme overestimated the shear modulus values when actual stiffness of the tissue was larger than 8 kPa. A similar effect was observed when estimating viscosity, albeit much earlier (i.e., the algorithm provided good estimates of of viscosity when μ < 5 kPa). Although the shear modulus affected the bias in the viscosity measurement, the authors demonstrated that the converse did not occur; i.e., the viscosity did not affect the bias in shear modulus. Despite these issues, their phantom studies revealed that inclusions were discernible in both μ and η -elastograms, and the viscosity values agreed with previously reported values for gelatin (0.21 Pa s). The patient studies revealed that the shear modulus values of malignant breast tumors were noticeably higher than those of benign fibroadenomas, but there was no significant difference observed in the viscosity of the tumor types, a result that would appear to contradict those reported in (Qiu et al., 2008).

Vappou et al. (2009) proposed a two-step approach for quantifying the viscoelastic properties of tissue. They measured the real component of the wave number (k’) from the phase ( ϕ) of the Fourier transform of the time-varying displacement at the excitation frequency:

They also measured the ratio of the real to the imaginary component of the complex shear modulus, from the phase between shear strain images and radiation force:

From R and k’ , they computed the k” (the imaginary component of the complex wave number) using the following relation:

Using this relation, they computed the storage ( G’(ω ) ) and the loss ( G”(ω )) as follows:

They validated their technique by performing phantom studies, and observed good agreement between the storage moduli that had been computed with the proposed technique, and those that had been measured with a commercially available rheometer. However, there was large discrepancy between the loss moduli.

Schmitt et al. (2010) used a similar approach to characterize the viscoelastic properties of vascular tissues. In their approach, they computed the real and imaginary components of the complex wave number from the Fourier transform of the real component of the complex time-harmonic displacement (i.e., u(x, y,t) = u e^i(k’+jαT) o e^iϕ ). They computed the attenuation coefficient by fitting a line to the natural log of the absolute value of the complex time-harmonic displacement. Storage and loss moduli were computed using using equation (60) as described in (Vappou et al., 2009), where k” = α . For materials with vascular geometries, they fitted the complex modulus that was computed using an analytic model to the measured data, using MATLAB’s nonlinear solver. They measured the dispersion of G’ and G” in aortic samples between 540-670 Hz, and conducted experiments on rat liver where G’(ω ) and G”(ω ) were 119.24 ± 61.6 Pa and 96.7 ± 7.9 Pa, respectively, at 250 Hz. The group also demonstrated that the viscoelastic properties of blood clots could be characterized using the proposed technique. However, in these studies of complex moduli, they estimated the storage and loss moduli by fitting a rheological model to the measured data. More specifically, they compared the goodness of fit of five rheological models—Maxwell, Kelvin-Voigt, Jeffery, Zener, and a third-order Maxwell—and observed that the Zener model gave the best fit to the data. G’ had the maximum chang at the beginning, and stabilized after 120 minutes; whereas G” was constant, with a spike between 38 and 81 Hz, followed by a gradual decrease in amplitude.

5.2.1 Quasi-static methods

Konofagou et al. (2001) introduced poroelastography. Using the finite-element method, they studied the temporal behavior of the radial-to-axial strain ratio within a homogenous, cylindrically shaped poroelastic material during stress relaxation. They predicted that three mechanical responses would be seen in poroelastic materials; specifically, that:

• (a) The radial-to-axial strain ratio elastograms would have a uniform value of 0.5 immediately after compression, which was consistent with the theoretical predictions made earlier by Armstrong et al. (1984) using an analytical model.

• (b) The radial-to-axial strain ratio within the sample would tend towards an equilibrium value that was equal to the Poisson’s ratio of the solid matrix (ν_s ) when the compression was sustained.

• (c) The time taken to reach equilibrium (the time constant) would depend on permeability of the solid matrix and the length of the fluid path — i.e., that materials with a high permeability would reach equilibrium more quickly than materials with a low permeability.

Righetti et al. (2004) demonstrated the feasibility of poroelastography experimentally. By studying the transient mechanical behavior of drained and undrained tofu samples (Wu, 2001), Righetti and colleagues demonstrated the following transient mechanical behaviors: (1) the mean radial-to-axial strain ratio of undrained tofu samples was approximately 0.5 immediately after compression, which was consistent with the numerical predictions of Konofagou et al. (2001); (2) the mean radial-to-axial strain ratios of drained and undrained tofu samples were noticeably different, but no significant difference was observed in rehydrated and undrained tofu samples; (3) the mean radial-to-axial strain ratio of undrained tofu samples decayed from 0.5 towards an equilibrium value that was slightly higher than Poisson’s ratio measured for the drained state; however, the rate of decay of the measured radial-to-axial strain ratio was noticeably slower than that predicted by the analytical model described in (Armstrong et al., 1984); and finally, (4) the mean radial-to-axial strain ratio of drained tofu samples displayed some transient mechanical behavior, albeit only slight, which suggests that besides being porous, the tofu samples were also viscoelastic.

Berry et al. (2006a) derived analytical expressions to predict the internal strain field of cylindrically shaped poroelastic materials that are undergoing stress relaxation. By extending the analytical model described in (Armstrong et al., 1984), they derived analytical expressions for predicting the spatio-temporal behavior of radial strain, axial strain, and radial-to-axial strain ratio elastograms. Using the derived equations, they predicted that immediately after compression, the radial-to-axial strain ratio would be equivalent to 0.5, and both the radial-to-axial ratio would decay to towards an equilibrium value that is equal to the Poisson’s ratio of the solid matrix when compression is sustained. This theoretical prediction was consistent with that predicted by Konofagou et al. (2001) using the finite-element method. They also demonstrated that there was good agreement between the mean radial-to-axial strain ratio predicted using the derived equation and those computed using either the analytical model described in (Armstrong et al., 1984), or the finite-element method. Using the modified analytical model, Berry et al. (2006a) predicted two mechanical behaviors that were not previously reported. More specifically, they predicted that (1) a plateau would be present in the radial-to-axial strain ratio profile in the early stages of stress relaxation, because fluid flow does not begin simultaneously throughout poroelastic materials; and that (2) near the surface, the radial-to-axial strain ratio would overshoot the Poisson’s ratio of the matrix analytical model to the time-dependent strain elastograms.

In a companion paper, Berry et al. (2006b) reported the results of experiments that they had performed to corroborate both the theoretical predictions reported in Berry et al. (2006a), and predictions made using the finite-element method concerning the effect of non-slip boundary conditions on the strain response within a poroelastic material undergoing stress relaxation. Using the finite-element method, they predicted that the transient strain response of cylindrically shaped poroelastic materials with slip boundary conditions would be noticeably different from those observed in samples with non-slip boundary conditions. More specifically, they predicted that: (1) immediately after compression and at equilibrium, region variations would be present in the radial, axial, and radial-to-axial strain elastograms obtained from poroelastic materials with non-slip boundary conditions; (2) axial strain elastograms would be time-dependent and spatially varying, because during stress relaxation, the axial strain in the central region would transfer to regions near non-slip boundaries; and (3) radial-to-axial strain profiles acquired from regions far away from the nonslip boundaries would be similar to those acquired for samples with slip boundary conditions. Berry et al. (2006b) also introduced a new elastogram known as the volumetric strain elastogram, which was the most useful elastogram for studying the transient mechanical behavior of poroelastic materials. They showed that, immediately after compression, this parameter was zero (incompressible) for samples with slip and non-slip boundary conditions. Although their experimental results departed from the theoretical predictions, they demonstrated that useful parametric images could be obtained experimentally using their model-based reconstruction approach. However, this approach was not without problems; local variations were observed in the parametric images obtained from a homogeneous sample, but the results were sufficiently encouraging to warrant further investigation.

The main limitations of stress-relaxation studies are that temporal variation is not typically observed in the axial direction during unconfined testing, and that the poor lateral resolution of current ultrasound systems compromises the quality of radial-to-axial strain elastograms. To address this issue, Righetti et al. (2007b) studied the temporal behavior of porous materials during creep. They demonstrated that the quality of transient axial strain elastograms was sufficiently high to differentiate between various grades of tofu samples. They generated time-constant elastograms by fitting an exponential function to each pixel in the time-sequence axial-strain elastogram. Using this technique, they demonstrated that they could detect focal poroelastic inclusions embedded in a poroelastic background with a contrast-to-noise ratio.

5.2.2 Harmonic Methods

Poroelasticity has also been explored using steady-state harmonic magnetic resonance elastography (Perrinez et al., 2010). More specifically, Perrinez et al. (2010) described a novel, 3D inversion scheme for recovering the mechanical parameters from volumetric displacement data obtained with MR in a porous media. They formulated the inverse-elasticity problem by treating it as a parameter-optimization problem, in which the goal was to minimize the difference between (1) measured displacements and (2) those computed by solving the partial differential equations that describe time-harmonic behavior of a poroelastic medium, which is composed of a porous, compressible, linear elastic solid matrix and a viscous, incompressible fluid. The resulting system of PDEs that describes the forward problem is given by:

The parameter β is defined as:

where p is the time-harmonic pore pressure; κ is the hydraulic conductivity; ϕ is the matrix porosity; and _ρf and _ρa represent the pore-fluid density and the apparent mass density, respectively.

Minimizing the objective function with respect to Lamé constants is a nonlinear process that can be solved using the Gauss-Newton iterative scheme. At each iteration, updates to the mechanical parameters are computed as follows:

where {δδ } = {δμ,δλ, δp}, ^J*T J is the self-adjoint Hessian matrix, α is the regularization parameter, J^* is the complex Jacobian or sensitivity matrix, and u^m and u^c are the measured and computed displacement fields, respectively.

Perrinez et al. (2010) solved equation (55) using the sub-zone inversion approach described in (Van Houten et al., 2001; Van Houten et al., 1999), which they implemented on a parallel-computing platform as described in (Doyley et al., 2004). They performed simulation and phantom studies to assess the performance of the proposed algorithm. The simulation study revealed that the algorithm could recover good elastograms (i.e., shear modulus, the second Lamé constant, and pore-pressure amplitude) in the presence of 5% additive noise, and that hydraulic conductivity influenced the performance of the reconstruction method. More specifically, the reconstruction method overestimated shear modulus when image reconstructions were performed using hydraulic conductivity less than 1×10⁻⁹, and considerable variability was observed in recovered images when hydraulic conductivities were greater than the true values. Perrinez et al. also performed a phantom study to assess the performance of the poroelastic reconstruction method relative to the linear-elastic reconstruction method described in (Van Houten et al., 2001; Van Houten et al., 1999; Doyley et al., 2004) . In this study, they perfoormed MR elastographic imaging on a phantom that contained a cylindrically shaped gelatin inclusion that was embedded in a block of tofu. Elastograms recovered with the linear -elastic inversion scheme showed high variations of the recovered shear modulus in both the surrounding poroelastic background and the linear elastic inclusion; whereas, the poroelastic-inversion scheme produced noticeably better images.

The advantage of the poroelastic inversion scheme as a reconstruction method is that it produces images of the pore-pressure that could prove useful in studying edema or hydrocephalus. Figure 8 shows examples of poroelastic elastograms obtained from the brain of a healthy volunteer using the poroelastic inversion scheme described in (Perrinez et al., 2010).

5.3.1 Geometric nonlinearity

When soft tissues deform by a small amount (an infinitesimal deformation), their geometry in the undeformed and deformed states is similar, and thus the deformation is characterized using engineering strain. However, when soft tissues experience a finite deformation, their geometries are noiceably different in the undeformed and deformed states. In such cases, errors are incurred when deformations are significantly different; thus, engineering strain provides an accurate measure of the deformation. To characterize finite deformation, we first have to define a reference configuration, which is the geometry of the tissue under investigation in either the deformed or undeformed state. The Green-Lagrangian strain tensor can be used to characterize the deformation incurred during finite deformations, which is defined as:

The nonlinear term is neglected when the magnitude of the spatial derivative is small, to produce the linear strain tensor as defined in equation (6). The relationship between stress and strain is nonlinear even for a linearly elastic material when it is undergoing finite deformations. Conseqently, Skovoroda et al. (1999) proposed an iterative technique to compute shear modulus for materials undergoing finite deformations.They performed studies using a linear elastic phantom that was undergoing finite deformation, to evaluate the quality of the ensuring elastograms relative to those produced using equation (22). The results of this investigation revealed that a smaller standard deviation was incurred in elastograms computed using the nonlinear reconstruction method for large deformations, versus elastograms computed using the linear-elastic reconstruction method (i.e., based on equation (22)).

5.3.2 Material nonlinearity

Some materials exhibit nonlinear material properties that are typically described using a strain energy density function. Among the strain engery functions proposed in the literature, the most widely used for modeling tissues are (a) the Neo-Hookean hyperelastic model, and (b) the Neo-Hookean model with anexponential term. Oberari et al. (2009) used a different model, the Veronoda-Westman strain energy density function, to describe the finite displacement of a hyperelastic solid that is undergoing finite deformation, which is defined by:

where the terms I₁ and I₂ are the first and second invariants of the Cauchy-Green strain tensor, μ₀ is the shear modulus, andΓ denotes the nonlinearity. The authors proposed to produce a secant modulus elastogram, by minimizing the following function:

Where ϑ is the total variation diminishing (TVD) regularization functional. For the nonlinear case, they proposed to reconstruct a nonlinear parameter, and the shear modulus at zero strain by minimizing the following functional:

In this case, there were two displacements: one obtained at small strain value, denoted by ${\mathbf{u}}_1^m$ and the other obtained at larger strain, denoted by ${\mathbf{u}}_2^m$; the corresponding predicated axial displacements were denoted by u₁ and u₂, respectively. The weighting factors w₁ and w₂ were selected to ensure that both the large deformation data and the small deformation data contribute in roughly equal measures.

Equations (58) and (59) were minimized using the quasi-Newton methods (.i.e., the Broyden-Fletcher-Goldfarg-Shanno method). They compared the performance of modulus elastograms computed using the secant and nonlinear iteriative inversion techniques, using data obtained from voluteer breast-cancer patients, one with a benign fibroadenoma tumor and the other with an invasive ductal carcinoma(IDC). For the fibroadenoma case, the tumor was visible in modulus elastograms that had been computed using small strain and large strain (12 %), although the contrast of the elastograms computed at large strain (7:1) was lower than that computed at smaller strain (10:1). The fibroadenoma tumor was not visible in nonlinear-parameter elastograms. The inclusion in the patient with IDC was discerible in shear modulus elastograms recovered using small and larger strains; however, the stiffness contrast of the modulus elastograms recovered at both small and high strains were comparable, and the IDC tumor was visible in nonliear-parameter elastograms. This result is one of several that have demonstrated the clinical value of nonlinear elastographic imaging. Specifically, model-based elastography can characterize the nonlinear behavior of soft tissues and may be used to differentiate between benign and malignant tumors. Figure 9 shows examples of shear modulus and nonlinear parameter images of the breast that were computed using the nonlinear reconstruction technique.

(a) Regularization

The ill-pose issue occurs when pertinent information about the solution is not available; therefore, the goal of regularization is to introduce a priori information, such as smoothness, in the reconstruction process. The two challenges we may encounter when using regularization are (i) selecting the most appropriate method, and (ii) deciding the optimum value of the regularization parameter. Too little regularization produces unusable modulus elastograms, while excessive regularization typically produces low-resolution modulus elastograms.

Discovering the most appropriate regularization technique for elastography is still an open research question, but the Tikhonov regularization is currently the most commonly used method. From our experience, elastograms produced with the total variation diminishing (TVD) regularization method usually possess better (contrast recovery and contrast-to-noise ratio) than elastograms computed using the Tikhonov regularization method, although the TVD regularization method typically does produces blotchy images (Richards and Doyley, 2011). The H1-seminorm regularization method could prove to be a better choice, since this regularization method typically produces elastograms that do not contain “blotchy” artifacts, its performance is comparable to modulus elastograms produced with the TVD method. However, before settling on any given regularization method more detailed studies must be conducted.

Developing objective methods to select the optimum value of the regularization parameter is another concern. To avoid the temptation of “tweaking”, the regularization parameter should be selected objectively using either the L-curve or the generalized cross-validation method (Vogel, 2002). However, since the L-curve and the generalized cross valiation methods are not appropriate for clinical applications, the author recommends that a statistical approach to the image reconstuction problem be employed, since statistical-based reconstruction methods provide a precise description of the regularization parameter (Van Houten et al. 2003).

(b) Spatial priors

In addition to regularization, other methods could transform the ill posed inverse elasticity problem to a well-posed one. For example, Barbone and Bamber (2002) suggested that model-based elastography be performed with multiple, independent displacement fields. Although feasible for quasi-static breast elastography, this method would be difficult to implement in other approaches to elastography. An alternate approach is to incorporate structural information in the image reconstruction process (Richards et al., 2010; Baldewsing et al., 2006; Baldewsing et al., 2005b; Le Floc’h et al., 2009; Le Floc’h et al., 2010). Structural information can be obtained by segmenting images obtained from ultrasound, MR, or other sources. In quasi-static elastography, spatial priors are typically used to impose hard constraint on the reconstruction process through a procedure known as parameter reduction (referred to here as “hard prior reconstruction”). More specifically, the shear moduli of all pixels in a given region as defined by the spatial prior are lumped together. Hard-prior reconstruction methods do not require regularization because the reconstruction problem is well conditioned; however, the technique is prone to errors because the segmentation of congruent features may contain classification errors. The author recommends that spatial prior should be used to impose soft constraints but not hard constraints, since this should minimize the effect of segmentation errors on the reconstruction process.