The combined importance of finite dimensions, anisotropy, and pre-stress in acoustoelastography.

Dynamic elastography, whether based on magnetic resonance, ultrasound, or optical modalities, attempts to reconstruct quantitative maps of the viscoelastic properties of biological tissue, properties that are altered by disease and injury, by noninvasively measuring mechanical wave motion in the tissue. Most reconstruction strategies that have been developed neglect boundary conditions, including quasistatic tensile or compressive loading resulting in a nonzero prestress. Significant prestress is inherent to the functional role of some biological tissues currently being studied using elastography, such as skeletal and cardiac muscle, arterial walls, and the cornea. In the present article, we review how prestress alters both bulk mechanical wave motion and wave motion in one- and two-dimensional waveguides. Key findings are linked to studies on skeletal muscle and the human cornea, as one- and two-dimensional waveguide examples. This study highlights the underappreciated combined acoustoelastic and waveguide challenge to elastography. Can elastography truly determine viscoelastic properties of a material when what it is measuring is affected by both these material properties and unknown prestress and other boundary conditions?

INTRODUCTION

Background and motivation

Dynamic elastography methods–based on optical, ultrasonic, and magnetic resonance imaging modalities–aim to quantitatively map the shear viscoelastic properties of biological tissue, which are often altered by disease and injury. Optical methods, including early work using stroboscopy to more recent, higher resolution methods using optical coherence tomography or laser Doppler vibrometry, have the advantage of the highest resolution of the three modalities but are constrained to the penetration depth of light, limiting their clinical use to measurements at or near the surface, such as the skin or cornea. Ultrasound (US)-based elastography using continuous or transient excitation has greater penetration depth and maintains high resolution near the surface (though not as high as optical methods), while also being readily available, contributing to its wide use in research since the late 1980s. Magnetic resonance elastography (MRE), introduced in 1995, has the highest depth of penetration, even behind hard tissue obstacles, such as the skull, and the ability to encode oscillatory motion in all three dimensions simultaneously; but it is limited in resolution as compared to optical methods and near-surface US methods and is the most expensive, as compared to all other modalities. Regardless of the imaging modality, dynamic elastography methods share common traits. They typically involve mechanical stimulation followed by measurement and analysis of resulting transverse wave motion in order to estimate or reconstruct the tissue's shear viscoelastic properties.

Most initial studies focused on larger organs, such as the liver or brain, where boundary effects were assumed negligible. But as elastography expands to other anatomical regions where dimensions in at least one direction are smaller or of comparable length to bulk shear wavelengths–such as in slender skeletal muscles, blood vessels, the heart wall, and the cornea–boundary effects become non-negligible and must be considered. Researchers using optical coherence elastography (OCE) to assess the viscoelastic properties of the cornea have long recognized this, adapting models to include waveguides by treating the cornea as a plate-like structure that is fluid-loaded on one side. Here, transverse wave motion on the cornea is modeled as Rayleigh–Lamb waves. Blood vessels, as well, have been modeled using cylindrical shell equations considering fluid–structure interaction. Limited studies on cardiac elastography have also acknowledged the frequency-dependent (i.e., wavelength-dependent) waveguide behavior of the heart wall.

Application of dynamic elastography to tissues with aligned fibrous structure resulting in local transverse isotropic mechanical properties, such as can be found in striated skeletal and cardiac muscle, as well as brain white matter, may benefit from analysis that takes into consideration anisotropy of the tissue. Recognizing this, some groups have pioneered research in this direction over the past few decades, using ultrasound (US)-based elastography, as well as magnetic resonance (MR)-based elastography. Many of these studies have tried to tackle the associated inversion problem. Multiple configurations or a multi-directional shear wave excitation source may be needed in order to generate and measure shear wave motion that will be affected by its displacement polarization direction and propagation direction in an anisotropic material.

In addition to a uniaxially aligned fibrous structure, a layered nonfibrous or fibrous structure, such as in the cornea, can also have a significant effect on shear wave behavior. The cornea can be considered as a transversely isotropic material where the axis of isotropy is perpendicular to the plane of the cornea. Here, it has been shown that both an in-plane and out-of-plane shear moduli are required to fully describe the elastic properties of corneal tissue. This important observation has provided a possible explanation for the differences in magnitude seen with elastic moduli measured using dynamic elastography techniques compared to traditional uniaxial tensile testing of cornea tissue.

Often, when elastography studies are done under varying quasi-static pre-stress conditions, observed changes in mechanical wave behavior are attributed solely to the nonlinear property of the tissue: it has been suggested that its shear and viscous constants are highly dependent on the tensile load and associated deformation. In previous muscle elastography studies, it has been surmised that the shear elastic modulus increased with passive muscle loading or with muscle activation. In other studies, it is simply observed that the shear wave velocity increased with increased passive or active loading of the muscle. While the latter is an indisputable observation–the velocity increased with increased load–the former (change in moduli value) ascribes the wavelength change solely to a change in tissue material properties, which we propose may or may not be responsible for a fraction of the observed changes in measured transverse wave speed.

Some have recognized the influence of compressive forces on elastography measurements of shear waves on phantoms and organs. An early study using MR elastography recognized that both tensile load and material elastic moduli affected transverse wave motion in skeletal muscle under tension. This was an in vivo study of the tibialis anterior (TA) and lateral gastrocnemius (LG), showing that shear wavelength increased when the muscle was stretched and when the muscle was contracted. In the Discussion section, a linear equation was put forth to account for both tensile and shear modulus effects, as an explanation for the observed changes. A more recent MR elastography study made a similar observation but replaced the bulk shear wave expression with one based on a Timoshenko beam under tension formulation; thus, also accounting for waveguide effects, and implemented this in an inversion strategy to assess tensile forces on the individual muscles of the forearm.

Objectives

In the present study, we review the theoretical principles of mechanical wave motion in a normally prestressed transversely isotropic material, and then consider two- and one-dimensional “thin” waveguides subjected to tensile loading, relating fundamental observations to applications in cornea and muscle elastography. Some of the presented concepts, without the supporting analyses detailed here, were previously summarized by the last author in an invited abstract and oral presentation at a recent meeting of the Acoustical Society of America.

PRE-STRESS UNDER SMALL DEFORMATION IN A TRANSVERSE ISOTROPIC MATERIAL

Building upon Tweten et al., we start with a linear elastic nearly incompressible, transversely isotropic (NITI) material as our model for biological tissue with aligned fibrous structure subjected to deformation that is sufficiently small in amplitude to justify the assumption of linearity. A linear elastic NITI material may be fully described using bulk modulus and three additional parameters which can be a combination of two tensile moduli, and two shear moduli, , where the subscripts denote whether the principle direction is perpendicular or parallel to the fiber direction. In other words, are in the direction perpendicular to the fibers (parallel to the plane of isotropy), and are in the direction parallel to the fibers (parallel to the axis of isotropy). We define shear anisotropy and tensile anisotropy . Note also that ; thus, there are only three independent parameters.

The transverse isotropic elasticity matrix using nomenclature from Guidetti et al. and taking the z-axis at the axis of isotropy is Adding only normal (no shear) static pre-stresses aligned with the x, y, and z directions, respectively, leads to the following governing equations of motion, where u, v, and w refer to the displacement component in the x, y, and z direction, respectively, and subscripted and t after a comma refer to partial derivatives with respect to that spatial or time dimension:

Consider plane wave propagation in the x, y, and z directions. Based on Eqs. (2)–(4) phase speeds for propagation in these three directions with two possible polarizations for shear waves, and one polarization for compression waves, is provided in Table I. (Note, this formulation, based inherently on a small strain assumption due to the applied stress and linear elastic theory, can be made to match predictions of phase speed in isotropic and transverse isotropic acoustoelastic models that allow for larger strain values by appropriate choice of third order Landau coefficients. For the isotropic case, this is achieved by setting the third order coefficient , which is within the range of values reported for agar-gelatin and polyvinyl alcohol soft tissue phantom materials based on ultrasound elastography measurements under compressive stress loading. For the transverse isotropic case, three additional third order Landau coefficients are set to specific values to match the small strain assumption used here. See the Appendix for additional explanation.)

TABLE I.

Phase speed squared density of planar waves as function of propagation and polarization directions.

Polarization →	x direction	y direction	z direction
↓ Propagation
x direction	κ+μ_⊥43+49ζ	μ_⊥+σ_x−σ_y2	μ_⊥1+ϕ+σ_x−σ_z2
y direction	μ_⊥+σ_y−σ_x2	κ+μ_⊥43+49ζ	μ_⊥1+ϕ+σ_y−σ_z2
z direction	μ_⊥1+ϕ+σ_z−σ_x2	μ_⊥1+ϕ+σ_z−σ_y2	κ+μ_⊥43+169ζ

ADDING BOUNDARY EFFECTS IN THE 3-DIMENSIONAL PROBLEM

Let us consider a point force or torque acting at the geometric center of a circular cylinder of a transversely isotropic (TI) material under a uniaxial stress aligned with both the axis of the cylinder and the axis of isotropy of the TI material, as shown in Figs. 1(a) and 2(a). For the case of plane wave propagation in this material, taking and taking , slow and fast shear wave speeds squared, respectively, are

FIG. 1.

(Color online) Fast shear waves in the y-z plane generated by a harmonic point force at the origin polarized parallel to the axis of isotropy (z-axis) and direction of uniaxial stress , using material properties from our previous work (Refs. 45 and 51–53) that correspond approximately to passive muscle with fibers along the z-axis (Table II). Simulations are for no pre-stress (b) and (c) or a uniaxial prestress from elongation along fiber direction that is equivalent to the stress from a 10% MVC contraction (d) and (e). All axes in (b) and (e) are in mm. Cylinders with radii of R = 20 mm (b) and (d) or R = 5 mm (c) and (e) are shown, each surrounded by a water-like material (low elastic moduli). Colors correspond to the amplitude of the shear waves, scaled by ( is the distance from the source), to compensate for attenuation away from the source. One quadrant is shown; the other three are mirror images. Here, -motion amplitudes (in the y-z plane) are plotted.

FIG. 2.

(Color online) Slow shear waves in the y-z plane generated by a harmonic point torque at the origin polarized orthogonal to the axis of isotropy (z-axis) and direction of uniaxial stress , using material properties from our previous work (Refs. 45 and 51–53) that correspond approximately to passive muscle with fibers along the z-axis (Table II). Simulations are for no pre-stress (b) and (c) or a uniaxial prestress from elongation along fiber direction that is equivalent to the stress from a 10% MVC contraction (d) and (e). All axes in (b)–(e) are in mm. Cylinders with radii of R = 20 mm (b) and (d) or R = 5 mm (c) and (e) are shown, each surrounded by a water-like material (low elastic moduli). Colors correspond to the amplitude of the shear waves, scaled by ( is the distance from the source), to compensate for attenuation away from the source. One quadrant is shown; the other three are mirror images. Here, x-motion amplitudes are plotted.

TABLE II.

Parameter values for case studies.

Parameter	Nomenclature	Value(s)
Bulk modulus	κ	2.6 GPa
Shear storage modulus in plane of isotropy	μ_⊥R	2.77 kPa
Ratio of shear loss to storage moduli	η=μ_⊥I /μ_⊥R=μ_∥I/μ_∥R	0.15
Shear anisotropy	ϕ	1
Tensile anisotropy	ζ	2
Uniaxial tensile stress/ μ_⊥R	σ/μ_⊥R	0 or 1
Density	ρ	1000 kg/m3
Frequency	f	1 kHz

Here, is the angle between the direction of propagation and the axis of isotropy and uniaxial stress.

A numerical integral solution for the response created for the unbounded case of a harmonically oscillating point force (infinitesimal dipole) in a TI material has been found, as well as analytical approximations to it. The authors are unaware if such a solution exists for the case with uniaxial stress or for the case of a point torque, or for any of these cases when the medium has finite boundaries. Consequently, for the following study, a numerical finite element (FE) approach is taken using ANSYS Mechanical APDL Version 2019 R1 (Ansys, Canonsburg, PA). The FE prediction was validated against the exact integral solution for the unstressed point force and taking the radius of the cylinder to be large enough that waves attenuate before reaching the boundary. The case study parameters shown in Table II are typical of soft biological tissue and match those used previously by the last author, but with the addition of cases with a nonzero uniaxial tensile stress (it is either 0 or equal to ).

First, a cylinder 70 mm in height and 40 mm in diameter was defined using an axisymmetric mixed u-P formulation with Plane183 8-node elements with individual element side dimensions of 0.1 mm. A point force was applied in the vertical “z” direction, parallel to the fiber and tensile load direction at the node located at the geometric center of the cylinder. Given the “z” polarization of the source input, it will predominantly drive fast shear waves. But, since this is a point source, the propagating wave field is not planar and thus, does not perfectly follow Eq. (2). The in phase steady state response for cases without and with a nonzero tensile load are shown in Figs. 1(b) and 1(d). Next, the point force is replaced by a point torque at the same location, oriented in the “z” direction, so that it predominantly drives waves that are polarized in the circumferential direction, which is orthogonal to the axis of isotropy and thus, behave like slow shear waves, approximately governed by Eq. (3). The in phase steady state response for cases without and with a nonzero tensile load are shown in Figs. 2(b) and 2(d).

In the simulations described above in this section, for the chosen material properties, cylinder diameter and excitation frequency, we see more than 8 wavelengths going from the excitation source (cylinder axis) to the cylinder free outer boundary at a radius of 20 mm. The effect of the free boundary is not apparent in this case where wavelength is an order of magnitude less than the characteristic cross-dimension. But consider this exact same case with the cylinder radius R reduced from 20 mm down to 5 mm, in other words, only a few wavelengths across. The simulations are shown in Figs. 1(c), 1(e), 2(c), and 2(e), identical to the cases in Figs. 1(b), 1(d), 2(b), and 2(d), respectively. While the wave pattern in the immediate vicinity of the source is similar, it has been substantially altered beyond a few wavelengths from the source. Differences caused by the presence of the uniaxial stress are still evident, but there are also complex wave patterns further from the source. These patterns are not described at all by Eqs. (5) and (6), although a modified form (to account for ) of the equations in Sec. 8.2.2 of Graff might apply.

TWO-DIMENSIONAL WAVEGUIDE WITH IN-PLANE BIAXIAL NORMAL STRESS

Theory

Next, we will consider a thin plate of thickness in the z-direction, with z still the axis of isotropy, but take: and take , as shown in Fig. 3. We will focus on shear waves in the x-y plane that are polarized in the z direction. In an unbounded three-dimensional medium, based on Eqs. (2)–(4) and Table I, the phase speed of plane waves propagating in the x-y plane with z polarization are

FIG. 3.

(Color online) Two-dimensional waveguide.

However, we will employ the following assumptions based on Mindlin's first contribution that follow from assuming that plate thickness is less than the shear wavelengths of interest (Sec. 8.3 of Ref. 54): . Furthermore, we will assume that , which is consistent with Euler–Bernoulli thin-plate theory. Utilizing Eq. (2.1) of Ref. 47, and inputting the above assumptions: Adapting Sec. 8.3 of Ref. 54, we multiply Eqs. (8) and (9) by z and integrate across the plate thickness from , and we directly integrate Eq. (10) across the plate thickness from . Neglecting higher order terms, this leads to where are bending moments about the plate, and are shear forces, and is an externally applied force per unit area to the plate (Fig. 3). Expressions for the bending moments are Neglecting rotational inertia terms in Eqs. (11) and (12) since h is small, taking partial derivatives with respect to x and y, respectively, solving for and and substituting into Eq. (13) leads to the following:

This can be written as where Neglecting external normal load per unit area and assuming a planar wavefront in the x direction, partial derivatives with respect to lead to zero and the equation simplifies to the following: The general solution to Eq. (20) for harmonic motion is given in the form , leading to the following set of solutions: Thus, we have two propagating waves (α) in the + or – x direction, and two non-propagating (near field or evanescent) waves (β) in the + or – x direction. For the propagating waves, the phase speed will be: . Taking the limit that we see that and thus, for the elastic case . Note, there is dispersion even in the lossless elastic case (neglecting viscosity); this is the classic unstressed thin transverse plate vibration solution. On the other hand, taking the limit of tension , we then drop from Eq. (20) and reformulate the solution to find there are two propagating solutions in the + or – x direction with . This is the classic thin membrane vibration solution. Note, none of these solutions match the expression for shear wave propagation under the same conditions but with h very large (3-dimensional medium) which is provided in Eq. (7).

Application to cornea elastography

The eye can be modeled as an internally pressurized spherical vessel where an increase in the internal pressure strains the walls of the ocular tissue and induces a tensile, circumferential hoop stress . The range of tensile stress can be roughly estimated using Laplace's law for a spherical pressure vessel as , where is the intraocular (gage) pressure (IOP), is the radius of the sphere, and is the thickness of the sphere walls. By assuming a radius of 10 mm and thickness of 0.6 mm, IOP values of 0, 5, 10, 15, and 20 mmHg will induce hoop stresses of 0, 5.56, 11.11, 16.67, and 22.22 kPa, respectively. For comparison, finite element (FE) studies using more precise geometry have reported hoop stresses within the range of roughly 15–25 kPa for around 15–18 mmHg IOP.

To evaluate whether tensile in-plane biaxial pre-stress, being equated to hoop stress, within this range can influence transverse wave motion on the cornea, we use equations from Sec. IV A with , , and , yielding . Pre-stress was incremented between 0 and 20 kPa to estimate the tensile effect of IOP on wavenumber [Fig. 4(a)] and phase speed [Fig. 4(b)]. For comparison, Young's modulus was also separately incremented between 15 and 75 kPa, with zero hoop stress [Figs. 4(c) and 4(d)].

FIG. 4.

(Color online) Thin-plate theory for comparing the effects of tensile prestress and stiffness on wave dispersion. (a) The linear wavenumber vs frequency dispersion curves demonstrate that estimated prestress values incremented between 0 and 2.0 × 104 Pa can alter wave behavior independently of material stiffness (where  = 3.0 × 104 Pa). (b) The corresponding phase velocity curves demonstrate that tension strongly influences low-frequency phase velocities while high-frequency waves are also affected. (c) Dispersion curves demonstrate that stiffness affects wave behavior similar to tension when incremented between 1.5 and 7.5 × 104 Pa (where  = 0 Pa), though its effects at high frequencies are more apparent. Additionally, when comparing (a) and (c), the low-frequency components below ∼1 kHz are affected less by stiffness than tension. (d) This is reflected in the corresponding phase velocity curves where low frequency phase velocities are less influenced by stiffness. Further comparison of (b) and (d) show that both tension and stiffness can change the high-frequency phase velocities that are used to estimate corneal stiffness.

Simulation results shown in Fig. 4 suggest that increasing both the tensile pre-stress (via increasing the IOP) and Young's modulus (a material property) can alter phase speed similarly. The slope of the dispersion curves decreases as both tension and material stiffness increase which corresponds to an increase in phase velocities. Importantly, this indicates that the high frequency asymptotes of the dispersion curves (above ∼3 kHz) that are normally used to estimate corneal shear moduli can be significantly influenced by both pre-stress and Young's modulus. While both effects are similar, differences in the dispersion curves can be observed, particularly in the lower frequency regions (below ∼1 kHz) from the approximate thin-plate model; tensile pre-stress appears to strongly influence the low-frequency dispersion behavior, while the effects of stiffness are more apparent at higher frequencies. These results suggest that tensile pre-stresses within the range predicted for normal and pathological IOP values can indeed influence wave behavior and that differential analysis of low- and high-frequency dispersion components may help to discern the effects of pre-stress and stiffness. It appears that if these pre-stress effects are not accounted for, an overestimation of corneal shear moduli may arise.

In Sun et al. (2021), analysis of the effect of tensile pre-stress and a foundation stiffness is further explored by integrating models for them into the full three-dimensional theory introduced in Sec. II. Trends identified using the thin plate theory are consistent with what is found using the more precise model, in that failure to account for the effects of pre-stress may result in overestimations of the corneal elastic moduli, particularly at high IOPs.

ONE-DIMENSIONAL WAVEGUIDE WITH UNIAXIAL NORMAL STRESS

Starting with the transverse isotropic material model introduced in Sec. II, consider a rod whose axis is aligned with the z-axis, the axis of isotropy, with cross-sectional dimensions in the x and y direction that are significantly less than wavelength. In other words, take the numerical study of Sec. III but with r smaller. Also, as in Sec. III, set and . Undertaking an analysis similar to Sec. IV, we arrive at the pre-tensioned Euler–Bernoulli thin beam described in Section 3.3.4 of Graff for x-polarized transverse wave propagation of the beam along its z-axis, Here, is the area moment of inertia about the y axis ( for a circular cross section of radius r), and is the cross-sectional area in the x-y plane. Note the similarity to the Eq. (20) of Sec. IV.

Consider this an initial approximation for a thin muscle under a tensile load along its axis. As before, the general solution form is: , where , is the circular frequency in rad/s and has four possible solutions, Like in Sec. IV A, we have two propagating waves (α) in the + or – x direction, and two non-propagating (near field or evanescent) waves (β) in the + or – x direction. For the propagating waves, the phase speed will be: . Taking the limit that we see that and thus for the elastic case , which is the classic thin (Euler–Bernoulli) beam transverse vibration solution. On the other hand, taking the limit of tension we then drop from Eq. (22) and reformulate the solution to find that there are two propagating solutions with . This is the classic transverse thin string vibration solution. As noted before in Sec. IV A, a case between the extremes of either neglecting or still does not match the phase speed of bulk shear waves, which again is given by Eq. (7).

Application to muscle elastography

In a recent study of shear wave elastography on an ex vivo cat soleus by Bernabei et al., where the tensile preload could be precisely measured while conducting elastography measurements on a passively tensioned muscle as well as when it is activated via an electrical current, it was found that under passive tensile loading, the wave speed closely followed the string solution with . However, under muscle activation phase speed appeared to fall between the two extreme cases articulated above (Sec. V A). Other studies have shown similar trends, but have ascribed the increase in to be due to an increase in the elastic moduli of the muscle, in other words its nonlinearity under deformation. The cat soleus has a physiological cross sectional area (PCSA) of about 1 cm2, or a diameter of about 1 cm. Human muscles studied using elastography can have a wide range of PCSAs, for example, as small as 0.68 cm2 in the abductor pollicis brevis to 21 cm2 in the medial gastrocnemius.

CONCLUDING REMARKS: NONLINEARITY AND THE GLASS HALF FULL

The experimental method of elastography generally involves small oscillatory deformations and it is reasonable to assume linear theory is valid. However, the static deformations, due to pre-stress in the applications described above (Secs. IV and V), may be significant and not be reasonably modeled based solely on linear systems theory. A logical approach is to apply a more robust nonlinear theory to model the static deformation created by the pre-stress using an appropriate strain energy function, followed by linearization about the new static equilibrium in order to derive an elasticity matrix that can be used for the “linearized” acoustoelastography problem. A few groups in the elastography community have been investigating this and identify higher order elastic constants. This remains an active area of research with a wide range of strain energy functions having been proposed, with it proving difficult to experimentally evaluate which may be most appropriate for a given application.

Even without the added complexity of geometric and material nonlinearities due to large static deformation, the elastography problem, especially in anisotropic materials, is made more challenging by the presence of static or quasi-static stress loads that are inherent to the normal physiological function of many biological soft tissues. They cannot be avoided. So is the increase in phase speed as a tissue is further stretched caused by the stretching (tensile load) itself, or is it caused by the higher order elasticity coefficients needed to fully characterize the material nonlinearity of the tissue over a large deformation range? Can these effects be teased apart through careful experimental measurements that fully investigate all manner of direction and polarization of mechanical wave motion at multiple static deformation levels? Are such experiments feasible? Recognition of these daunting technical challenges may lead one to think that the technique of elastography is, indeed, a half empty glass, unable to confidently determine that which it originally set out to determine, a material's elastic constants. But this would be missing the forest for the trees. Being clear-eyed about the challenges that we face, we believe the glass is half full, in that once we overcome the known complexities of the problem, acoustoelastography will be able to noninvasively quantitatively map not only a material's complex stiffness, but also the complex stress field it is under. Stiffness and stress are both critical to the physiological function of many biological tissues, and both can be uniquely altered by pathology, injury, and response to therapy. A tool that can nondestructively measure both independently would be profound.