Biaxial Tensile Prestress and Waveguide Effects on Estimates of the Complex Shear Modulus Using Optical-Based Dynamic Elastography in Plate-Like Soft Tissue Phantoms Free

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. Most initial studies have 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 elastography 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. Here, transverse wave motion on the cornea is modeled as Rayleigh–Lamb waves [1]. Blood vessels, as well, have been modeled using cylindrical shell equations considering fluid-structure interaction and Lamb wave dispersion [2–4]. Limited studies on cardiac elastography have also acknowledged the frequency-dependent (i.e., wavelength-dependent) waveguide behavior of the heart wall [5,6].

Often, when elastography studies are done under varying nonzero quasi-static prestress 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. A recent article provides a summary of the literature relevant to this issue, in particular for uni-axially prestressed cylindrically-shaped structures, as well as biaxially prestressed plate-like structures [7]. In another article, the impact of neglecting the prestress effect in cornea elastography is quantified [8].

In the present study, focused on biaxially prestressed plate-like structures, the confounding effects of finite dimensions and prestress are further explored, analytically and experimentally. Furthermore, we articulate and evaluate a strategy for decoupling prestress and waveguide effects from estimates of material shear viscoelastic properties.

Here, $ksh[ω]$ is the complex-valued, frequency-dependent shear wave number and $μ[ω]$ is the complex-valued, frequency-dependent shear modulus, comprised of the shear storage modulus, $μR[ω]$⁠, and the shear loss modulus, $μI[ω]$⁠, such that $μ[ω]=μR[ω]+jμI[ω]$ where $j=−1$⁠. The attenuation rate of the wave as it propagates is governed by the imaginary part of the wavenumber: $Imag[ksh[ω]]$⁠. In a viscoelastic material, both the shear storage and loss moduli affect both the phase speed and attenuation rate. In a purely elastic material, $μI=0$⁠, there is no attenuation and the phase speed is independent of frequency (nondispersive) and reduces to $c=μ/ρ$⁠. While some studies have assumed pure elasticity (no viscosity), often their analyses are generalizable to the linear viscoelastic problem for harmonic motion by adding the imaginary shear loss modulus to form the complex shear modulus.

Here u, v, and w refer to the displacement component in the x, y, and z direction, respectively, and subscripted $x, y, z, and t$ after a comma refer to partial derivatives with respect to that spatial or time dimension. The term $κ$ is the bulk modulus of the material, which will affect compression wave behavior, but not shear wave behavior. In biological soft tissue or other “nearly incompressible” materials, $κ$ is more than a few orders of magnitude greater than $μ$ such that the Poisson's ratio for the isotropic material approaches, but does not equal 0.5.

Here, $θ$ is the angle between the direction of propagation and the z-axis.

where $I¯1=J−23I1$⁠, with $I1=tr[C]$⁠, $C=FTF$⁠, $J=Det[F]$⁠, and $Jm$ is a limiting parameter for $I¯1−3$⁠. Here, $F$ is the deformation gradient, $C$ is the right Cauchy-Green tensor, and $J=1$ if the material is incompressible. The first term of the strain energy function is based on the isochoric deformation of the isotropic material, and the second term only exists if the material is compressible. The shear modulus is given by $2∂ψ∂I1$⁠, which in the incompressible limit is equal to $μ$ and is consistent with infinitesimal strain theory. Consequently, Eqs. (5) and (6), governing planar shear wave phase speed should remain valid here for nearly incompressible materials. Utilization of other strain energy functions is left for future study. It is also possible to incorporate the prestress field directly into the strain energy function [19,20]. If it is assumed that the superimposed wave motion on top of the finite deformation is small, one may still recover Eqs. (5) and (6) for planar shear wave phase speed.

Applying the nonhomogeneous boundary condition, Eq. (9), we expect only the lowest order solution from Eq. (14) to be mechanically driven, and the corresponding wavenumber $kRL$ will be close in value to, but slightly larger than $ksh$⁠, associated with a slightly slower wave propagation speed and shorter wavelength.

Note, in this equation the only “r” dependence is within $J0(kRLr)$⁠.

A novel approach is proposed to account for biaxial prestress aligned in the plane of the plate. In previous studies on transversely isotropic materials not under prestress [23–26] the last author has shown that, by distorting the geometry based on direction and polarization-dependent planar phase speeds one can then solve an equivalent isotropic problem. This approach to the transversely isotropic problem was called “Transformation Elastography”. Biaxial prestress causes a similar, though not identical, direction and polarization dependence of the planar shear wave phase speed, as exhibited in particular in Eq. (6). The same approach is adapted to the acoustoelastic problem here. Note, different geometric distortions are needed depending on whether the wave motion of interest is a slow or a fast shear wave, based on propagation direction and polarization. The out-of-plane axisymmetric wave motion detailed in the previous subjection matches the criteria for a fast shear wave.

When the formulas are expressed in terms of $ϵR$⁠, it is clear that one does not need to know $σ$ or $μ$ a priori in order to calculate the equivalent dimensions. One only needs know the radial strain. If the compressibility is non-negligible, Eqs. ((16) and (17)) will need some adjustment.

A numerical case study was conducted to understand the interactions between biaxial prestress and waveguide behavior, as well as to validate the proposed transformation acousto-elastography (TAE) approach proposed above. Case study geometric and material property values are provided in Table 1. These material property and geometry values are chosen to match those for the material and fixture used in experimental studies described later in this article.

An axisymmetric finite element (FE) nonlinear static deformation study followed by a time-harmonic analysis based on linearization of the stiffness matrix after deformation was conducted using ansys Mechanical APDL 2020 revision 2 (ansys, Canonsburg, PA). A plate-shaped structure with radius R = 25.4 mm and thickness H = 8.5 mm was defined using Plane183 eight-node elements with individual element side dimensions of ∼0.1 mm. This resulted in 16,341 nodes and 5,334 elements. (Use of different element dimensions confirmed resolution was sufficient for these studies.) A mixed u-P (displacement-pressure) solution approach was used due to the near-incompressibility of the element material properties, given in Table 1, in order to avoid locking, which can occur when using a displacement formulation that is dependent on volumetric strain calculations that are difficult to accurately calculate when Poisson's ratio is near 0.5. In this formulation, as opposed to a more efficient formulation based only on displacements, the hydrostatic pressure or volume change rate is interpolated on the element level and solved on the global level independently in the same way as displacements. The material elasticity was defined using a Gent model, which only requires specification of the shear storage $μR$ and bulk $κ$ moduli (specified in ansys as a compressibility factor $2/κ$⁠), as well as the parameter $Jm$⁠. With the nonlinear deformation analysis enabled, a radial displacement boundary condition was applied to the nodes on the outer edge of the model, with displacements in the r direction of 0, 0.635, 1.27, 2.54, and 5.08 mm in five separate simulations, resulting in percentage strains in the radial direction of 0%, 2.5%, 5%, 10%, and 20%, respectively. After application of each deformation, a time-harmonic perturbation analysis was performed by applying a displacement to the same nodes on the outer edge of the phantom, but in the out-of-plane z direction with amplitude of 10 microns and at a frequency f = 300 Hz. For the harmonic (dynamic) analysis, in addition to shear storage and bulk moduli, it is also necessary to specify density $ρ$ and shear loss modulus $μI$⁠, which is specified via beta damping in ansys, where the beta damping term will equal the ratio of shear loss to storage moduli, $η$⁠, divided by the frequency in radians/second, $2πf$⁠. The static displacement and stress conditions due to the initial radial deformation, and the dynamic displacement, due to the harmonic excitation, were output for every node into text files, which were imported into matlab for analysis.

In summary, we have articulated here a strategy to estimate the material's shear storage and loss moduli based on elastography measurements acquired in the presence of a known static deformation, but potentially unknown prestress condition and unknown material moduli. Equations ((16) and (17)) are expressed in terms of the prestrain, independent of the prestress and material properties.

Predictions based on the TAE concept to estimate material shear viscoelastic property values at 300 Hz, $μR+jμI$⁠, independent of prestress conditions are in Table 2. The case of no prestress (0% prestrain) serves as a check on the accuracy of the finite element model and matlab -based least squares curve fit in the process given above in this section. In other words, this level of error in estimate is independent of the TAE approach. Whereas the actual value of the complex shear modulus is $μR+jμI=39.840+j6.542 kPa$⁠, the prediction at 0% prestrain is $40.1+j6.78 kPa$⁠. The “Estimate Error” shown in the table is the difference of the calculated and actual complex moduli, divided by the actual complex modulus. It is 0.78% and 0.47% for estimates of the shear storage and lost moduli, respectively. For 2.5%, 5%, 10%, and 20% radial prestrain cases, the TAE approach yields calculations of the complex shear modulus with the error in the shear storage and loss moduli predictions usually within 2%, except reaching 2.15% in the case of 10% prestrain.

This is compared to the biaxial stress value $(σFEA)$ determined from the finite element simulation, which is the averaged radial stress in all elements of the model. For the given loading condition, the stress was uniform throughout the model. Here, the percent error in prestress prediction stays within a few percent reaching a maximum of 2.03% when the radial prestrain reaches 20%.

To experimentally evaluate TAE-based estimates of material viscoelastic properties, a test configuration similar to, though not the same as the numerical study of Sec. 3 was constructed. This configuration was inspired by earlier studies of our group that identified the confounding effect of intraocular pressure (IOP) on estimates of cornea shear elasticity properties based on optical coherence elastography measurements [8]. It was shown that in-plane biaxial tensile stress in the cornea due to nonzero IOP affected measured Rayleigh–Lamb wave motion induced on the cornea, independent of the cornea shear elastic properties.

The experimental configuration, shown in Fig. 2, includes a circular plate-like phantom 8.5 mm thick and with 25.4 mm radius made of a moldable polymer, Ecoflex-30™. The polymer plate is fixed around its perimeter to a cylindrical container that is vibrated axially to create a radially converging, ideally axisymmetric, flexural wave pattern on the phantom material, similar to the FEA study in Sec. 2. Simultaneously, the cylindrical container is statically pressurized, which both creates some curvature of the phantom and places it under a nonzero axisymmetric multi-axial and nonhomogeneous prestress field. A scanning laser Doppler vibrometer (SLDV; PSV-400, Polytec, Waldbronn, Germany), is used to measure out-of-plane flexural motion on the phantom, as described in previous studies [21,27].

From the Ref. [28], it was found that for Ecoflex-10, the following values accurately describe shear viscoelastic properties over the range of 200 Hz to 7.75 kHz: $μ0=13.3 kPa$⁠, $μα=2 kPa·sα$⁠, and $α=1/3$⁠. For Ecoflex-30, measurements conducted over the range of 200 Hz–1 kHz yielded: $μ0=27 kPa$⁠, $μα=1.5 kPa·sα$⁠, and $α=0.3$⁠. Thus, while Ecoflex-30 is “stiffer” (higher $μ0$⁠) under static conditions, due to its lower viscosity the magnitude of its complex shear modulus increases at a slower rate with frequency (⁠$α=0.3$⁠), as compared to that of Ecoflex-10 (⁠$α=1/3$⁠).

The measured flexural wave motion on the phantom material using the SLDV system is exported and read into matlab for further analysis. Typical experimental measurements of the Rayleigh–Lamb wave motion are shown in Fig. 3. As in the numerical study of Sec. 3, the “lsqcurvefit” command was used to optimize the fit of Eq. (15) to the experimentally-measured out-of-plane harmonic displacement along a radial line of the phantom (averaged from measurements along 40 radial lines spaced 9 degrees apart). Unlike in the previous numerical study, it is recognized that in the experimental case, the curvature of the plate means that the SLDV is not measuring motion exactly normal to the plate surface at all locations. The SLDV “imposes” a planar geometry with a radius R that is constant, not changing with increasing pressure-induced prestrain. To account for this in the TAE analysis the initial least squares curve fit to determine $kRL$ is multiplied by $(1+ϵR)$⁠. In the experimental case the radial strain, $ϵR$⁠, is estimated based on its correlation to the out-of-plane deformation caused by the static pressure, which is measured using a micrometer. This relationship is established via finite element studies using the phantom geometry. (The experimental measurement of vertical displacement with the micrometer was judged to be accurate to the nearest 0.1 mm with confidence, though the micrometer precision was 0.01 mm. We believe this translates to estimates of radial deformation of the order of 0.1 mm, which in turn provides strain precision and accuracy to 0.4%.)

Next, dynamic small amplitude axisymmetric excitation is applied to the phantom along its outer boundary and dynamic displacement profiles along 20 equally angularly spaced diameters were averaged, and then used to identify best-fit values for $kRL$ and $wza$⁠.

Once $kRL$ is determined, the Rayleigh–Lamb wavelength $λRL=2π/Real[kRL].$ Next, as in the numerical study both $λRL$ and $kRL$ are scaled to estimate the undeformed case via TAE using Eqs. ((18) and (19)). The complex-valued shear wave number $ksh/TAE$ and the shear wavelength $λsh/TAE$ in the undeformed state then are estimated via solution of Eq. (14) with $He$ from Eq. (17) in place of H, again using a least squares optimization algorithm. From $ksh/TAE$ an estimate of the material's complex shear modulus, $μR+jμI$⁠, is found using Eq. (20).

In summary, this is a strategy to estimate the material's shear storage and loss moduli based on elastography measurements acquired in the presence of a known static deformation, but unknown prestress conditions and unknown material moduli. Once the material moduli are estimated, then this information can be used to estimate the averaged biaxial stress in the phantom, based on the measured estimate of the radial strain, $ϵR$⁠.

The approach described in this section to estimate the complex shear modulus and biaxial stress is first tested using a finite element simulation that mimics the experiment by applying a pressure load to one side of the plate that is statically fixed at its perimeter. It is then applied to experimental data.

Predictions based on the TAE concept applied to the “curved surface” numerical study to estimate material shear viscoelastic property values at 300 Hz, $μR+jμI$⁠, independent of prestress conditions are in Table 3. (As in the “flat surface” studies summarized in Table 2, the “ground truth” value for the complex shear modulus is $μR+jμI=39.840+j6.542 kPa).$ As compared to the flat surface numerical study of Sec. 3, the TAE approach does not do quite as well in enabling estimates of the complex shear modulus independent of the preload condition, with the error reaching 6% when the preload pressure reaches 100 cm H₂O, which induces a radial prestrain of $ϵR=7.03%$⁠. Also, it is found that estimates of the radial biaxial prestress, $σTAE$⁠, based on Eq. (21), are not directly relatable to the actual prestress condition, as the actual condition is not simply biaxial and is not homogenous; it varies with both radial and axial location. In addition, the prestress field is not simply in the radial r direction or in the local tangential plane of the phantom; there is a significant through-plane predominantly z component, as well as a significant x-z shear stress component. Figure 4 is a plot of the principle stress obtained from the FE static simulation that is predominantly aligned with the radial direction but follows the local tangential plane of the phantom.

Predictions based on the TAE concept applied to experimental data to estimate material shear viscoelastic property values at 300 Hz, $μR+jμI$⁠, independent of prestress condition are in Table 4. For pressure-induced radial prestrain cases, the TAE approach yields calculations of the complex shear modulus with the error in the shear storage and loss moduli predictions of no more than 6% and 16%, respectively, over a range of pressure-induced prestrains up to 11.3%. And, as in the numerical studies, the biaxial radial stress $(σTAE)$ is estimated using Eq. (21), though in light of the curved surface numerical study, the unknown actual prestress field for this case is expected to be more complex and not simply biaxial.

The numerical finite element and experimental studies of the previous two sections have quantified the confounding effects of finite dimensions, surface curvature and nonzero prestress on the elastography approach to estimating material viscoelastic properties in an isotropic plate-like structures under biaxial normal stress within the plane of the plate and under more complex but still axisymmetric stress conditions. Additionally, a novel coordinate transformation approach, TAE, has been introduced that, in theory, enables one to estimate material viscoelastic properties independent of the nonhomogeneous prestress conditions without requiring a priori knowledge of those stress conditions. Rather, only the amount of deformation, or strain, from the unstressed condition is required. If the prestress field is uniform, the method can then estimate the prestress.

The numerical and experimental studies show the promise of the TAE approach. In the flat surface numerical study, the material shear storage $μR$ and loss $μI$ moduli were determined with accuracy from less than 0.5% up to 2.15% error as prestrain levels increased to 20%. Subsequent estimates of the uniform radial stress field $σ$ were comparably accurate over this prestrain range, reaching 2.03% error at 20% prestrain. The TAE approach does not account for geometric nonlinear effects that increase with deformation and thus, this may be a significant source of the error in estimates of both the shear moduli and the prestress. Another source of error or discrepancy is an incompressibility assumption used to formulate the approach. Because the bulk modulus is about five orders of magnitude greater than the shear moduli, we believe the assumption of incompressibility is justified. Though, the approach could be reformulated to account for compressibility. Accounting for geometric and material nonlinearity, which increase as prestrain increases, is left for future development of the TAE approach. Also left for future numerical simulation studies, is the robustness of the TAE approach when confronted with simulated noise in the measurements.

In the curved surface numerical study, which more closely approximates the experimental study, estimates of the complex shear modulus independent of prestress condition reach 6% error at the highest pressure considered—100 cm H₂O, which generated a 7% radial strain. The estimate of a biaxial prestress value is less meaningful in this case, as the finite element study shows that the stress field is not simply biaxial and is not homogeneous. Rather, there is variation as a function of both radial r and axial z position (Fig. 4).

Finally, in the experimental study, one must also consider that the static shear modulus governing the initial deformation due to the prestress, is typically different from the frequency-dependent shear storage modulus that will govern response to harmonic excitation. This, and the fact that there may be some residual stress even in the unstressed condition, is why the deformation induced by the pressure loading is significantly greater than that in the curved surface numerical study. Despite this, estimates of the complex shear modulus under nonzero prestress conditions up to 11% prestrain, generally stayed close to each other and within 10% of the expected value. The value predicted in the unstressed case had the greatest difference and it is thought this may be due to unmodeled residual stresses caused by the assembly of the experimental fixture. Based on finite element simulations, we do not believe that the fact that the SLDV measured motion of the surface of the plate as opposed to the centerline was a major source of discrepancy in estimating wave numbers.

Logical next steps for the strategy introduced here for decoupling prestress from material shear stiffness estimates include consideration of more complex geometry and stress conditions, as well as anisotropic and nonuniform material properties. Conducting studies at multiple frequencies would enable identification of a rheological model that could help discern static (constant) versus frequency-dependent components of the shear storage and loss moduli. Use of torsional excitation of waves polarized in the plane of the plate, would excite the slow shear waves (Eq. (5)), which would provide an additional means of estimating the shear modulus [29]. Finite element models based on medical images that can provide detailed geometry [30], deformation under varying loading conditions and, if needed, measures of anisotropy and inhomogeneity, may provide a way to advance the TAE technique beyond simple geometries and assumptions of isotropy and homogeneity.

• National Science Foundation (NSF) (Grant No. 1852691; Funder ID: 10.13039/100000084).

• National Institutes of Health (NIH) (Grant No. AR071162; Funder ID: 10.13039/100000069).

There is no conflict of interest.