Measurement and imaging of viscoelasticity, fat and collagen content of soft tissues

Abstract

Various examples are provided related to viscoelasticity, fat and collagen determination in soft tissues. In one example, a method for evaluation of soft tissue includes obtaining measured particle velocities from generated mechanical waves in tissue of a subject; determining viscoelastic parameters of soft tissue at a point or in a region of interest (ROI) from wave propagation characteristics derived from the measured particle velocities; and determining fat or collagen content at the point or in the ROI from the wave propagation characteristics. In another example, a system includes an ultrasound scanner for SWE; and a computing device that can obtain measured particle velocities from generated mechanical waves in tissue of a subject; determine viscoelastic parameters of soft tissue at a point or ROI from wave propagation characteristics derived from the measured particle velocities; and determine fat or collagen content at the point or in the ROI.

Description

Docket: 221407-2140 MEASUREMENT AND IMAGING OF VISCOELASTICITY, FAT AND COLLAGEN CONTENT OF SOFT TISSUES CROSS REFERENCE TO RELATED APPLICATIONS

[0001] This application claims priority to, and the benefit of, co-pending U.S. provisional application entitled “Measurement and Imaging of Viscoelasticity, Fat and Collagen Content of Soft Tissues” having serial no.63 / 652,389, filed May 28, 2024, which is hereby incorporated by reference in its entirety. STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH OR DEVELOPMENT

[0002] This invention was made with government support under HL145268 awarded by the National Institutes of Health, and DMS2111234 awarded by the National Science Foundation. The government has certain rights in the invention. BACKGROUND

[0003] There has been growing interest in estimating tissue viscoelasticity using ultrasound and magnetic resonance imaging given its role as a non-invasive biomarker for many diseases. Different elastography techniques such as transient elastography (TE, FibroScan), magnetic resonance elastography (MRE), and ultrasound shear wave elastography (SWE) have been used to characterize disease progression. SUMMARY

[0004] Aspects of the present disclosure are related to viscoelasticity, fat and collagen determination in soft tissues. In one aspect, among others, a method for evaluation of soft tissue comprises obtaining measured particle velocities from generated mechanical waves in tissue of a subject; determining viscoelastic parameters of soft tissue at a point or in a region of interest (ROI) from wave propagation characteristics derived from the measured particle velocities; and determining fat or collagen content at the point or in the ROI from the wave propagation characteristics. The viscoelastic parameters can be determined using a twinpeaks method by minimizing mismatch between the measured ^^^^^^and ^^^^^^peaks and thepredicted ^^^^^^ and ^^^^^^ peaks. The predicted ^^^^^^ and ^^^^^^ peaks can be obtained as analyticalpredicted ^^^^^^ and ^^^^^^ peaks can be generated from particle velocities obtained from simulating SWE experiments. A complex-valued viscoelastic modulusDocket: 221407-2140 can be estimated directly as an arbitrary function of frequency. Viscoelasticity can be parametrized in a functional form, with the viscoelastic parameters estimated through inverse optimization. The parametrization can capture direct mechanical behavior comprising elasticity and viscosity. Viscoelasticity can be parametrized in a functional form, where the viscoelasticity is linked to the fat or collagen content, which is directly estimated through inverse optimization. Spatial maps of the fat or collagen content can be estimated through inverse optimization and relationships linking fat and collagen contents to mechanical properties. The viscoelastic parameters or the fat or collagen content can be determined using a convolutional neural network (CNN) based at least in part upon the measured particle velocities, the CNN trained over multiple acoustic radiation force (ARF) widths. The CNN can be selected from a plurality of trained CNNs based upon a minimum misfit between synthetic and measured responses. The training can be performed either by a single CNN model across multiple ARF widths or by multiple CNN models, with one CNN model associated with each ARF width.

[0005] In another aspect, a system for evaluation of soft tissue comprises an ultrasound scanner configured for shear wave elastography (SWE); and a computing device comprising a processor and memory. The computing device can be configured to at least: obtain measured particle velocities from generated mechanical waves in tissue of a subject; determine viscoelastic parameters of soft tissue at a point or in a region of interest (ROI) from wave propagation characteristics derived from the measured particle velocities; and determine fat or collagen content at the point or in the ROI from the wave propagation characteristics. The viscoelastic parameters can be determined using a twin peaks method by minimizing mismatch between the measured ^^^^^^ and ^^^^^^ peaks and the predicted ^^^^^^ and ^^^^^^ peaks. The predicted ^^^^^^ and ^^^^^^ peaks can be obtained as analytical relations. The predicted ^^^^^^ and ^^^^^^ peaks can be generated from particle velocities obtained from simulating SWE experiments. A complex-valued viscoelastic modulus can be estimated directly as an arbitrary function of frequency. Viscoelasticity can be parametrized in a functional form, with the viscoelastic parameters estimated through inverse optimization. The parametrization can capture direct mechanical behavior comprising elasticity and viscosity. Viscoelasticity can be parametrized in a functional form, where the viscoelasticity is linked to the fat or collagen content, which is directly estimated through inverse optimization. Spatial maps of the fat or collagen content can be estimated through inverse optimization and relationships linking fat and collagen contents to mechanical properties. The viscoelastic parameters or the fat or collagen content can be determined using a convolutional neural network (CNN) based at least in part upon the measured particle velocities, the CNN trained over multiple acoustic radiation force (ARF) widths. The CNN can be selected from a plurality of trained CNNs based upon a minimum misfit between synthetic and measured responses. The training can beDocket: 221407-2140 performed either by a single CNN model across multiple ARF widths or by multiple CNN models, with one CNN model associated with each ARF width.

[0006] Other systems, methods, features, and advantages of the present disclosure will be or become apparent to one with skill in the art upon examination of the following drawings and detailed description. It is intended that all such additional systems, methods, features, and advantages be included within this description, be within the scope of the present disclosure, and be protected by the accompanying claims. In addition, all optional and preferred features and modifications of the described embodiments are usable in all aspects of the disclosure taught herein. Furthermore, the individual features of the dependent claims, as well as all optional and preferred features and modifications of the described embodiments are combinable and interchangeable with one another. BRIEF DESCRIPTION OF THE DRAWINGS

[0007] Many aspects of the present disclosure can be better understood with reference to the following drawings. The components in the drawings are not necessarily to scale, emphasis instead being placed upon clearly illustrating the principles of the present disclosure. Moreover, in the drawings, like reference numerals designate corresponding parts throughout the several views.

[0008] FIGS.1A-1C illustrates examples of a homogeneous bulk, ultrasound shear wave elastography (SWE) experimental setup, and measured particle velocity response, in accordance with various embodiments of the present disclosure.

[0009] FIGS.2A-2D illustrates an example of visualization of ^^^^^^ and ^^^^^^ peaks, in accordance with various embodiments of the present disclosure.

[0010] FIGS.3A-3D illustrate examples of twin peaks for low and high viscosity mediums, in accordance with various embodiments of the present disclosure.

[0011] FIGS.4A-4C illustrate examples of twin peaks simulation results, in accordance with various embodiments of the present disclosure.

[0012] FIGS.5A and 5B illustrate examples of twin peaks for noise-free in silico data, in accordance with various embodiments of the present disclosure.

[0013] FIGS.6A-6D illustrate examples of twin peaks for noisy in silico data, in accordance with various embodiments of the present disclosure.

[0014] FIG.7 illustrates an example of spatiotemporal response, in accordance with various embodiments of the present disclosure.

[0015] FIGS.8A-8C illustrate examples of twin peaks for in silico data from spring pot model, in accordance with various embodiments of the present disclosure.Docket: 221407-2140

[0016] FIGS.9A-9E illustrate an example of ex vivo validation of twin peaks method, in accordance with various embodiments of the present disclosure.

[0017] FIGS.10A-10C illustrates an example of problematic data sets for ex vivo acquisition, in accordance with various embodiments of the present disclosure.

[0018] FIGS.11A-11E illustrate an example of ex vivo inversion of the last two data sets, in accordance with various embodiments of the present disclosure.

[0019] FIGS.12A-12D illustrate an example of an in vivo application, in accordance with various embodiments of the present disclosure.

[0020] FIGS.13A and 13B illustrate examples of the effect of acoustic radiation force (ARF) width on shear wave velocity response, in accordance with various embodiments of the present disclosure.

[0021] FIG.14 illustrates an example of a convolutional neural network (CNN) that can be used for viscoelastic parameter estimation, in accordance with various embodiments of the present disclosure.

[0022] FIGS.15A and 15B illustrate examples of mean absolute percent error (MAPE) confusion matrices for shear modulus (^^) and relaxation time (^^), in accordance with various embodiments of the present disclosure.

[0023] FIG.16 illustrates an example of a confusion matrix showing ARF width recovery accuracy, in accordance with various embodiments of the present disclosure.

[0024] FIGS.17A-17B and 18A-18B illustrate examples of prediction error (signed and absolute) across test set, in accordance with various embodiments of the present disclosure.

[0025] FIGS.19A and 19B illustrate a comparison of mean absolute precent error using a combined CNN with unknown ARF width versus proposed selection approach, in accordance with various embodiments of the present disclosure.

[0026] FIG.20 is a schematic block diagram illustrating an example of a system employed for estimation of viscoelasticity, fat and collagen content of soft tissue, in accordance with various embodiments of the present disclosure. DETAILED DESCRIPTION

[0027] Disclosed herein are various examples related to viscoelasticity, fat and collagen determination in soft tissues. Reference will now be made in detail to the description of the embodiments as illustrated in the drawings, wherein like reference numbers indicate like parts throughout the several views.

[0028] Transient elastography (TE) primarily measures only elasticity which has lead researchers to focus on magnetic resonance elastography (MRE) and ultrasound shear wave elastography (SWE) for viscoelasticity inversion, and despite the accuracy of MRE, theDocket: 221407-2140 cost associated with MRE is significantly higher than TE or SWE. Thus, the focus is on SWE given that it is widely available, inexpensive, and can estimate both elasticity and viscosity.

[0029] Many initial efforts using ultrasound-based methods to characterize tissue viscoelasticity are based on evaluating the phase velocity and assessing the dispersion, i.e., the variation of the phase velocity with frequency. The dispersion curves (phase velocity as a function of frequency) can be measured with different signal processing approaches such as the phase gradient and two-dimensional (2D) Fourier transform, as well as some additional recently proposed methods; the resulting curves can be fit with different rheological models of viscoelasticity. The idea of dispersion matching was extended and the difference between the shear wave speeds from particle displacements, velocities and accelerations utilized to invert for tissue viscoelasticity.

[0030] As an alternative to shear wave dispersion, shear wave velocity and attenuation can be used together to characterize the complex-valued viscoelastic modulus. Compared to shear wave velocity, shear wave attenuation can be difficult to estimate because of the need to incorporate information about the shear wave source. Many acoustic radiation force (ARF) push beams can be approximated as a cylindrical source and the attenuation due to the diffraction of the source needs to be accounted for in the estimation process. A method for reconstructing images of shear wave velocity and attenuation has been proposed by accounting for these diffraction effects. The AMUSE method was developed to estimate the phase velocity and attenuation, after incorporating diffraction correction. In subsequent work, the reduced accuracy of the AMUSE method was noted due to truncation of signals and proposed finite window, 2D Fourier transform (fw2DFT) to remove this limitation. Later, the two-point frequency shift approach was proposed, building on an initial contribution utilizing the ideas of amplitude spectra frequency shift to measure attenuation. This study was further extended for better characterization of attenuation using power spectra frequency shift of shear waves measured at two spatial positions. A modified frequency shift approach was also used to prevent the outlier attenuation values in the presence of noise.

[0031] The other related contributions include the use of convolutional neural networks and full wave inversion. Despite the many contributions, viscoelasticity estimation is inherently difficult given that the adverse effects of noise are amplified due to the signal decay associated with tissue viscosity. Focusing on approaches based on 2D Fourier transform, as discussed before, one can consider various quantities to estimate elasticity and / or viscosity: (a) signal spread, (b) peaks (dispersion curves), or (c) the entire signal as done for arterial viscoelasticity. Of these, peaks, i.e. dispersion curves, are most robust against noise and thus have been routinely used for estimating elasticity. Unfortunately, dispersion curves alone may not be sufficiently sensitive for accurate viscoelasticity inversion. In this disclosure, it is observed that there exists another peak, different from the traditional dispersion curve,Docket: 221407-2140 providing the necessary sensitivity to viscosity while retaining the robustness of peaks with respect to noise. Specifically, first note that in the presence of viscoelasticity, the peaks of the particle velocity response (in the Fourier domain) can be defined in two different ways: ^^^^^^ peaks, which are wavenumbers associated with maximum response while sweeping through frequencies, and ^^^^^^ peaks which are correspondingly the frequencies at the maximum response while sweeping through wavenumbers. Importantly, these two peaks diverge, and the extent of the divergence depends on tissue viscosity, indicating that matching both these peaks would lead to inversion of both elasticity and viscosity. Thus, this can be called the twin peak method (TPM), which is the main subject of the disclosure.

[0032] A new model-based inversion approach is presented using SWE that can incorporate any viscoelastic model and predict the viscoelasticity accurately. The “^^^^^^ peaks” and “^^^^^^ peaks” or twin peak matching can be used to invert for the viscosity and elasticity. The approach was validated with in silico and ex vivo porcine liver data and it was applied to in vivo human liver to demonstrate the precision and efficacy of the method for estimating viscoelasticity and thereby contributing to the enhancement of biomarker. After defining the problem setup and theoretical background, the mathematical rationale behind the TPM method is summarized. The details of the methodology describe a forward model that takes in viscoelastic parameters as input and computes simulated ^^^^^^ and ^^^^^^ peaks. The forward model is then used in an iterative inversion framework that translates experimental ^^^^^^ and ^^^^^^ peaks to estimates of viscoelasticity parameters. Finally, the TPM’s accuracy is examined using in silico and ex vivo experiments, followed by in vivo application. Methods

[0033] Shear Wave Elastography: Problem Setup: A model-based optimization approach is proposed to estimate viscoelastic material properties using Shear Wave Elastography (SWE). Shear waves are generated using acoustic radiation force (ARF), applied by a focused ultrasound beam at a specific region of the tissue. FIG.1A is a graphical representation of a homogeneous bulk. ARF generates shear waves that propagate within the tissue, and the resulting particle velocity is measured in the ^^ directionon the ^^ െ ^^ plane. The objective is then to process these measurements to estimate thetissue viscoelasticity. FIG.1B illustrates an experimental setup of SWE. The resulting waveis measured as particle velocity in the ^^ direction in the ^^ െ ^^ plane as shown. Consideringthe noise in the signal, the particle velocity is often averaged in the ^^ direction over a small 3 mm deep horizontal strip at the center depth of the push, resulting in the particle velocityrepresentation in the ^^ െ ^^ domain. FIG. 1C shows an ^^ െ ^^ representation of a measuredparticle velocity response. The objective is then to process the ^^ െ ^^ representation ofparticle velocity to estimate the tissue viscoelasticity.Docket: 221407-2140

[0034] In this study, the tissue is considered to be a homogenous bulk medium, and ARF is essentially a time-dependent body force. The resulting displacement can be obtained by the incompressible elastodynamics equation with the field variable being the three- dimensional (3D) displacement vector. Given the dominance of shear waves and the relatively long shape of the ARF along the ^^ axis, the equation is often simplified to a wave equation with the scalar field variable ^^ representing the displacement in the ^^ direction: ^^^^^ െ ^^ ∗ ^^ଶ^^ ൌ ^^. (1)^^ is the density, the over ^^^^^^ modulus in shear (inverse, ^^ଶis the Laplacian, and ^^ represents the ARF. The ∗ symbol is the convolution operator capturing the memory effects of viscoelasticity. The ARF is characterized by its shape in space ^^^, multiplied with a rectangular step function in time ^^௧, i.e., ^^^^^,^^, ^^, ^^^ ൌ ^^^^^^,^^, ^^^^^௧^^^^, (2)^^ ൌ ^1 0 ^ ^^ ^ ^^,௧0 ^^^^ℎ^^^^^^^^^^^^,(3)where ^^ is the duration of the ARF, which is chosen to be 400 microseconds in this study.The velocity ^^ ൌ ^^^ , depends on the unknown viscoelastic modulus ^^. The goal is to estimatethe modulus ^^ from v , measured along the ^^ axis.

[0035] Basic Idea: Twin Peak Method (TPM): Important to the approach is theobservation that the ^^ െ ^^ particle velocity, examined in the frequency ^^^^ – wavenumber^^^௫^ domain has a peak with wide spread compared to that from an elastic medium. This spread, quantified using the full-width at half-maximum, was used to estimate viscosity but is prone to noise. The spread can instead be quantified by examining two peaks, so called ^^^^^^ and ^^^^^^ peaks, which are distinct and are expected to be robust against noise. FIG.2D illustrates the two diverging peaks. The ^^^^^^ peaks are obtained from examining the ^^ െ^^ transform of the particle velocity and, for each wavenumber ^^௫, finding ^^, or equivalentlycyclic frequency ^^ ൌ 2^^^^, associated with the maximum absolute value, i.e.,^^^^^௫^ ൌ ^^^^^^^^^^^^|^^^^^௫,^^^|, (4) where ^^ is the Fourierplotting |^^^^^௫,^^^| after normalizing for each ^^௫by themaximum value. This will make the maximum value for each ^^௫to be exactly one, visually highlighting the peak. This is shown in FIG.2B. Note that such normalization and plotting are used only to illustrate the idea and not in the actual methodology.

[0036] The ^^^^^^ peaks are similarly obtained by locating the wavenumber ^^ corresponding to maximum response for each ^^ (or ^^), visualized again with appropriateDocket: 221407-2140 normalization as shown in FIG.2C. The two peaks are shown together in FIG.2D, which clearly illustrates the spread, indicating that viscoelasticity can be estimated by fitting both peaks. While at the end, a standard optimization-based inversion technique is utilized to estimate viscoelasticity from these peaks, in the remainder of the section, a basic mathematical insight is provided into the method with the help of a simplified one- dimensional (1D) wave propagation model.

[0037] Specializing Equation (1) to an infinitely long axisymmetric push leads to a 1Daxisymmetric wave equation written in polar coordinates. Applying √^^ modulation, i.e. ^^^ ൌ√^^^^ leads to a standard 1D wave equation that governs the far-field: ^^ డమ௨^ డ௧మെ ^^ ∗డమ௨^ డ௫మൌ ^^௧^^^^^^^^^^, (5)where the ARF shapehere that the ^^coordinate coincides with the radial direction and hence used interchangeably. Fourier transforming Equation (5) in ^^ and ^^ results in:  ^^^ 2U ^ Gk 2x U ^ F ( ^ ) ,  (6) where ^^ and ^^௫ areF ^ F ( kx , ^ ) are the Fourier transforms of the displacement and ARF respectively, G isthus the frequency-dependent complex shear modulus, which is the Fourier transform of therelaxation modulus. Focusing on the simpler Kelvin-Voigt rheological model, G is given by,G ^ G 0 (1 ^ i^^ ) ^ ^c2 s(1 ^ i ^^ ) , (7)where ^^^ is the(the ratio between viscosity and elastic modulus). Further, idealizing ARF to be an impulseresult in F ( ^ ) ^ 1, leading to the frequency-wavenumber representation of thedisplacement: U( k ,^) ^ 1^ 1(8)The particleas, ^^^^^ ൌ ^^^^^^ ൌ^ఠ . (9)

[0038] ^^௫, obtain the frequency at which the amplitude of the response is maximum, leading to ^^^^^^Docket: 221407-2140 peak, and (b) obtain ^^^^^^ peak by traversing along ^^ (or ^^), and locating ^^௫associated with maximum response: ^^^^^^ peak: ^^^^^௫^ ൌ ^^^^^^^^^^^^|^^^^^௫,^^^|, (10) ఠ ^^^^^^ peak: ^^௫^^^^ ൌ ^^^^^^^^^^^^|^^^^^௫ ,^^^|. (11)^ Starting with Equation (10), the ^^^^^^ peak can be obtained by finding ^^ that maximizes the absolute value of the right hand side in Equation (9). This is equivalent to minimizing the sum of squares of the real and imaginary parts of its inverse. Such minimization is done bytaking the partial derivative with respect to ^ and equating it to zero, i.e.,డ ^ ଶ ೞమ^^మఠ ^^ ఠെ ^^^ ^ ^^ସ^ ସ ଶడ^ ^௫^^ ൨ ൌ 0 , (12)which results in the^^^^^௫^ ൌ ^^^^^௫. (13)Similarly, the ^^^^^^ peaks can be obtained setting the partial derivative with respect to ^^ as zero: డ ^ೞమ^^మଶ ସସ ଶ^ ^^െ ^^^ ^ ^^^^^௫^^ ൨ ൌ 0, (1)resulting in the^^௫^^^^ ൌ. (15) ^ೞ√^ାఛమఠమ

[0039] The ^^^^^^ anddifferent, except for the purely elastic case, i.e. ^^ ൌ 0. Furthermore, the deviation betweenthe two peaks is a function of ^^, indicating that ^^, thus the viscosity, can be estimated from the peak deviation. Elasticity can be obtained from the slope of the dispersion curve represented by Equation (13), which is the classical approach for estimating bulk tissue elasticity. It can thus be concluded that both elasticity and viscosity can be estimated by matching both ^^^^^^ and ^^^^^^ peaks. Moreover, peaks are expected to be independent of the ARF magnitude and less sensitive to noise compared to the width of the peak.

[0040] Equations (13) and (15) can be used to match the experimental peaks in FIGS. 2B and 2C, respectively, but have several shortcomings: (a) the analytical peaks in Equations (13) and (15) are limited to the Kelvin-Voigt model, which may not be accurate for real tissues; (b) Equations (13) and (15) assume an idealized impulse for ARF, which may not be the case in reality; (c) diffraction correction through√^^ modulation is not accurate in the near-field and in the presence of viscosity; and importantly, (d) signal processing detailsDocket: 221407-2140 such as sampling and truncation in space and time affect the spread and thus the location of the experimental peaks but not the analytical peaks.

[0041] To overcome the above limitations, the approach to viscoelasticity estimation does not directly use Equations (13) and (15), but builds on the idea that the experimental peaks (FIG.2D) are sensitive to unknown viscoelasticity parameters. Specifically, an inverse optimization approach is used where the particle velocity is obtained from a forward modelon the ^^ െ ^^ grid that is identical to the experimental measurement. The simulated peaks areobtained from the simulated response using a procedure identical to that for computing the experimental peaks. The mismatch between the experimental and simulated peaks are quantified using least-squares error, which is iteratively minimized to estimate the viscoelastic parameters. Details of (a) the forward model to obtain the simulated response, (b) the approach to compute ^^^^^^ and ^^^^^^ peaks, followed by (c) the inversion approach used to solve for the unknown viscoelasticity parameters are provided. Details are provided of the procedure and parameters used for in silico, ex vivo and in vivo studies performed in this study.

[0042] Forward modeling for response computation: Given the homogeneous, unbounded approximation for the tissue, the modeling approach is analytical in nature. Specifically, the Fourier transform is first used to convert the problem into a set of decoupled (pseudo) differential equations in time. Depending on the rheological models, these equations are solved either directly in the time domain, or in the frequency domain. Both approaches are presented, followed by a more efficient 2D spatial approximation, which makes the inversion process more practical with respect to the computational effort.

[0043] Frequency-Wavenumber Approach for General Viscoelasticity. Fourier transforming Equation (1) in both space and time results in, ^^^ 2U ^ Gk 2 U ^ F , (16)where ^^ ൌ ^^^^^,^^^ and F ^ARF respectively. ^^ ൌ ^^^௫ , ^^௬, ^^் ௭^ is the wavenumber vector, where ^^௫ ,^^௬ , ^^௭ represent thespatial frequencies in ^^, ^^, ^^ directions. ^^ ൌ ^^^௫ଶ ^ ^^௬ଶ ^ ^^௭ଶ is the magnitude of thewavenunber vector. Equation (16) immediately results in the particle displacement and velocity fields: ^F^^i^ FFinally,domain velocity response:Docket: 221407-2140 ^^^^^, ^^, ^^, ^^^ ൌ ^^^^^^^^^^^^,^^^^. (18)Note that, for this frequency-domain simulation, there are no restrictions on the underlying viscoelastic model and one can consider spring-pot, Kelvin-Voigt, fractional Voigt, or a completely general viscoelastic model (along with appropriate parametrization needed for inversion). The procedure is implemented in MATLAB (Mathworks, Natick, MA), which is also used for all the other computations performed in this disclosure. Other platforms such as Julia or Python can be utilized as well.

[0044] Time-Wavenumber Approach for Kelvin-Voigt Model: While the above frequency-wavenumber approach has general applicability, time-wavenumber approach would be more efficient for the special case of the Kelvin-Voigt model and is presented here. For the Kelvin-Voigt model, the viscoelastic modulus can be written in an operator form as, ^^ ൌ ^^^^ଶడ ^^1 ^ ^^. (19) డ௧ ^ Substituting the above in Equation (1), and applying the Fourier transform in space results in, ^^ డమ௨ഥ డ௧మ^ ^^^^ଶ^ ଶ ଶ ଶడ௨ഥ ^^ ^ത^ ^ ^^^^^ ^^ ^^డ௧ൌ ^^, (20)where, ^̄^ ൌ ^̄^^^^, ^^^ and ^^ ൌ ^^^^^, ^^^ are the spatial Fourier transforms of the ^^ and ^^respectively. The above essentially represents a damped vibration problem, i.e.,^^^ ^డ௨ഥ డ௧మ^డ௧^ ^^^ത^ ൌ ^^ , (21)with mass ^^ ൌ ^^, damping ^^ ൌ ^^^^^ଶ^^ଶ^^ and stiffness ^^ ൌ ^^^^^ଶ^^ଶ. The natural frequency ^^^ ൌ^^^⁄ ^^ ൌ ^^^^^ and the damping ratio ^^ ൌ ^^⁄ 2√^^^^ ൌ ^^^^^^^⁄ 2. Considering that ARF is arectangular pulse in time as mentioned in Equation (3), the resulting particle velocity can be derived as, v F kt ^ ^ ^1t ^ ^2t^ t ^ t^^^,ଶ ൌ ^^^^െ^^ േ ^^^ଶ െ 1^. (23)We obtain the finalFourier transform in space: ^^^^^, ^^, ^^, ^^^ ൌ ^^^^^^^^̄^^^^, ^^^^ (24)Docket: 221407-2140

[0045] Two-Dimensional (2D) Approximation: While the above subsections address 3D simulation in space, noting that the ARF push is often long in the ^^ direction and the response is often used only at the mid-depth, the simulation can be simplified from 3D to 2D. Essentially, the push, thus the response, are assumed to be independent of ^^ and vary onlyin the ^^ െ ^^ plane. This leads to a 2D version of Equation (1) in the ^^ െ ^^ plane. Whentransformed into the Fourier domain, only wavenumbers ^^௫ and ^^௬remain, with ^^௭ ൌ 0 (sinceno variation in ^^ direction). The remaining details of the formulation in the previous section stay the same. Such a simplification leads to significant savings in the computational cost (e.g.5.5 seconds for 2D analysis as opposed to 1250 seconds for 3D, on a 12௧^Gen IntelCore i9-12900k computer with 64 GB RAM). These computational cost savings appear tocome with minimal effects on accuracy as discussed in the Results section. As illustrated in FIG.4B in the Results section, 2D inversion of synthetic data generated from 3D forward model does not result in any significant errors, leading us to advocate the use of 2D forward models for TPM inversion. Identifying Peaks

[0046] The TPM is based on matching the experimental peaks with simulated peaks. Thus, the measured particle velocities from SWE data as well as from the computed response using the forward modeling procedure described above, must be converted to appropriate pair of ^^^^^^ and ^^^^^^ peaks. Care is exercised in (a) reducing the effect of noise in the experimental data, and (b) avoiding any unintended discrepancies in the definition of peaks that could lead to errors in the eventual estimation of the viscoelasticity parameters.

[0047] Focusing first on the experimental peaks and associated noise effects, start withthe ^^ െ ^^ data as discretely sampled in the SWE motion data, with appropriate values of ^^^^,^^^^, ^^^^௫and ^^^^௫. A standard procedure in, e.g., dispersion-based inversion is to correct for geometric spreading with√^^ modulation. Based on experience with experimental data, such modulation is not used for TPM due to its effect on amplifying noise away from the ARF. More importantly, since TPM is being used by matching with the simulated peaks from a consistent forward model and not the simplified analytical peaks in Equations (13) and (15), such modulation is no longer necessary. At the end, the signal is appropriately padded andFourier transformed in space and time to get the ^^ െ ^^ motion data. ^^^^^^ peak is then thefrequency ^^ corresponding to the maximum absolute ^^ െ ^^ response for any givenwavenumber ^^௫. To visualize this peak, one can normalize absolute ^^ െ ^^ data for each ^^௫,by the maximum for that ^^௫. See e.g. FIG.2B. A similar procedure can be followed for computing and visualizing ^^^^^^ peaks (see FIG.2C). Note again that normalization is used only if visualization is desired; it is not a part of the TPM method.Docket: 221407-2140

[0048] For computing simulated peaks, it is important to keep in mind that many of the processing steps described for experimental data processing has effects on the peak locations, e.g., signal truncation in experimental data translates to convolution with the Sincfunction in ^^ െ ^^ domain leading to increased spread and potential overestimation ofviscoelasticity. Such complexities can be avoided by obtaining simulated peaks in the same way as the experimental peaks, e.g., the response obtained from the forward model is sampled and truncated exactly the same way as the SWE measurements, followed by usingthe same windowing and padding parameters to obtain the ^^ െ ^^ data that is used forobtaining the simulated peaks.

[0049] Frequency and wavenumber ranges: Peak matching in TPM is performed in a particular frequency and wavenumber range that are dependent on acquisition parameters and expected range of material properties. While specific examples are provided in the results section, some guidelines are provided to determine these ranges.

[0050] The resolutions Δ^^ and Δ^^ determine the upper limits of the wavenumber and frequency ranges, respectively. Thus, Δ^^ and Δ^^ can be required to be not very large inorder to have sufficient range of ^^ െ ^^ signal for peak matching. At the outset, it may appearthat the lower limit of the ^^ െ ^^ range should be governed by ^^^^௫ and ^^^^௫ due to spectralleakage effects, but it turns out that the ranges can go to much lower. This is aided by the fact that the approach utilizes the exact same signal processing parameters for both measurement and simulated data, indicating that all the signal distortions would be identical between the two and can be confidently compared. Noting however that the signal distortions can dominate the peak deviation resulting in reduced sensitivity, as an ad hoc criterion, fix the lower limit to approximately the points where the two peaks cross each other.

[0051] Once a good enough ^^ െ ^^ range is determined for the particle velocity signalbased on acquisition parameters, the range of reliable twin peaks need would depend on the expected material properties, especially viscosity as measured through relaxation time. For amaterial with low viscosity, higher ^^ െ ^^ range is needed to see measurable deviation (seean example of twin peaks for low viscosity medium ^^^  ൌ 0.2 ^^^^^ in FIG. 3A). On the otherhand, for a material with higher viscosity, the deviation can be observed in lower ^^ െ ^^ range(see an exampe of twin peaks for high viscosity medium ^^^  ൌ 0.8 ^^^^^ in FIG. 3B). This ismathematically evident through Equations (13) and (15), where peak deviation is a function of the product of frequency ^^ and relaxation time ^^.

[0052] FIGS.3A and 3B also indicate that the deviating peaks have significant artifacts once they exceed a particular frequency and wavenumber, well below the range determined by Δ^^ and Δ^^. This phenomenon is attributed to spatiotemporal window effects (GibbsDocket: 221407-2140ringing) and can be explained by examining the absolute value of ^^ െ ^^ signals which areprocessed to obtain the peaks. These are shown in FIGS. 3C, ^^ െ ^^ response for lowviscosity ^^^  ൌ 0.2 ^^^^^, and 3D, ^^ െ ^^ response for high viscosity ^^^  ൌ 0.8 ^^^^^. Essentially,Gibbs ringing from low ^^ and ^^ signal pollutes the signals in high ^^ – low ^^ and low ^^ – high ^^ regions. For high viscosity materials, the signal decays faster with frequency and wavenumber, compared to low viscosity matrials. This allows the polluted signal to dominate, leading to peak jumping and other artifacts seen in FIGS.3A and 3B.

[0053] In summary, the ^^ െ ^^ range of peaks are governed by Δ^^, Δ^^, windowingapproach and viscoelastic material properties. In addition, to address the effects of noise in real data, the repeatability of the two peaks from across various replicates of SWE experiments is also used to determine the ranges for reliable experimental peaks. Theprocess of considering the above effects to determine the ^^ െ ^^ range for peak matching iscurrently performed in an ad hoc manner. Automation can be performed after more robustly quantifying uncertainties arising in real signals, where various windowing and other signalprocessing ideas can be used to optimize the ^^ െ ^^ range for reliable TPM inversion.Inversion using Twin Peaks

[0054] The unknown viscoelastic parameters of the tissues can be estimated by matching simulated peaks with measured peaks. This can be done by minimizing the objective function ^^^^^: ^^^^௩ ൌ ^^^^^^^^^^^^^^^^^, (25) ^^ where ^^ represents the vector of viscoelastic parameters and ^^^^௩is the final estimate of the parameter vector. Define the objective function to be the sum of the relative least-squares differences between measured and simulated twin peaks: ^N1N^ ^2 2 ^^2^where ^^^^^^^ and ^^^are simulated peaks, and ^^^and ^^^^^^are measured peaks. The minimization in Equation (25) is carried out using the BFGS method, through MATLAB’s "fminunc" function. Given the limited number of parameters, a finite difference derivative was utilized with a relative step size of 10ିହ. In silico StudiesDocket: 221407-2140

[0055] In silico verification was performed not only as a natural first step before validation with real experimental data, but also to further examine the accuracy of inversion based in 2D models when applied to in silico data obtained from 3D models, which is more representative of real experimental data. Start with a detailed inversion for the Kelvin-Voigtmodel with ^^^ ൌ 2 ^^^^^^ and ^^ ൌ 0.5 ^^^^, then expand to other levels of viscosity by changingthe relaxation time per Table 2 in the Results section. For 3D models, the following are usedΔ^^௫ ൌ 39 ^^^^^^ / ^^, Δ^^௬ ൌ 40 ^^^^^^ / ^^, Δ^^௭ ൌ 35.5 ^^^^^^ / ^^, Δ^^ ൌ 33.6 ^^^^^^ / ^^, ^^௫,^^௫ ൌ12,600 ^^^^^^ / ^^, ^^௬,^^௫ ൌ ^^௭,^^௫ ൌ 6,280 ^^^^^^ / ^^, ൌ 31,416 ^^^^^^ / ^^ for simulation. For 2Dmodels, the same parameters were used (Δ^^௭, ,^^௫no longer applicable). At the end,the spatio-temporal response is obtained with Δ^^ ൌ 0.25 ^^^^ and Δ^^ ൌ 0.1 ^^^^, and istruncated to ^^^^௫ ൌ 30 ^^^^ and ^^^^௫ ൌ 15 ^^^^, representative of typical experimenal data.Also consistent with experimental data processing, the data were extended using zero padding in ^^ and ^^ to result in 2^ସpoints in both directions; this resulted in Fourierresolutions of ^^^^ ൌ 0.61 ^^^^ and Δ^^ ൌ 1.53 ^^^^^^ / ^^.

[0056] A synthetic ARF with ^^ / ^^ ൌ 1.25 was used as the forcing function. The ARFprofile was computed using Field II, which was then used in the forward model using 3D time-wavenumber formulation. The parameters for the L7-4 linear array transducer (Philips Healthcare, Andover, MA) were used for the modeling of ARF. This is a 128-element array transducer with element width of 0.283 ^^^^, pitch of 0.308 ^^^^, element height of 7 ^^^^, andelevation focus of 25 ^^^^. A beam was focused at ^^ ൌ 25 ^^^^ using 64 active elements in amedium with longitudinal wave speed of 1540 ^^^^ି^and an attenuation of 0.5 ^^^^ / ^^^^ / ^^^^^^. The pressure and intensity were simulated for a three-dimensional region to be used as an input to the forward models. The intensity was used as a scaled version of the ARF.

[0057] The response at ^^ ൌ 0 , i.e., at the center of the ARF, is used to generate the insilico peaks, which can then be inverted using the TPM with the 2D forward model. The verification can then be performed with varying viscosity, i.e., with relaxation time (^^) ranging from 0.1 ^^^^ to 1.0 ^^^^. The procedure discussed in the Methods section was used to pick thefrequency and wavenumber ranges. For lower viscosity, i.e., 0.1  ^ ^^ ^ 0.5 ^^^^, peakmatching can be performed for 100 ^ ^^ ^ 500 ^^^^,  500 ^ ^^௫ ^ 1500 ^^ି^. For higherviscosity, i.e., 0.5 ^ ^^ ^ 1 ^^^^, peak matching can be performed for 100 ^ ^^ ^ 200 ^^^^,300 ^ ^^௫ ^ 800 ^^ି^.

[0058] Next, the inversion was repeated with in silico data polluted with noise: ^^^^^^௬ ൌ ^^^௬^௧^^௧^^^1 ^ ^^^^^^ ^ ^^^^^ ൈ ^^^^^^^ ^^^௬^௧^^௧^^^, (27)wheregenerated from the standard normal distribution, i.e., with zero mean and unit standardDocket: 221407-2140 deviation. A wide range of ^^^and ^^^were chosen as represented in Table 1 to perform the in silico verification with wide range of noise levels. Table 1. Noise levels used for in silico verification Noise level Noise Quantity Additive, ^^^Multiplicative, ^^^(%) (%)0 0 Free Low 0.35 15 Medium 0.70 30 High 1.00 40

[0059] The compatibility of the TPM was also checked with the spring-pot viscoelastic model, i.e., ^ ఈ ^^ ൌ ^^ఠ ^^^ , (28)where ^^^ is the reference(7) and (28) are not the same). The ^^^is an arbitrary reference frequency introduced tomaintain dimensional consistency. The sampling parameters Δ^^, Δ^^, ^^^^௫,^^^^௫ and the ^^ െ^^ bin sizes are considered identical to Kelvin-Voigt model discussed previously. For thisstudy, choose ^^^ ൌ 2 ^^^^^^ and ^^ ൌ 0.4 and ^^^ ൌ 2^^ ൈ 200 ^^^^^^ / ^^. In silico data wereobtained by detailed 3D simulation of wave propagation using the frequency domain forwardmodel, followed by adding noise per Equation (27), again with ^^^ ൌ 0.30 and ^^^ ൌ 0.007(medium noise level), as informed by visual comparisons with real experimental data. Theresponse at ^^ ൌ 0, i.e., at the center of the ARF, is used to generate the in silico peaks.These peaks were then inverted for 100 ^^^^ ^ ^^ ^ 500 ^^^^ and 500 ^^ି^ ^ ^^௫ ^ 1500 ^^ି^,using the 2D version of the forward model. Ex vivo Studies

[0060] The TPM was validated with ex vivo testing on a porcine liver. Liver samples were acquired from pigs after euthanasia. The pigs were dedicated for medical education or cardiovascular research on protocols approved by the Mayo Clinic Institutional Animal Care and Use Committee. SWE experiments were performed using a Verasonics V1 system (Verasonics, Inc., Kirkland, WA) equipped with a linear array transducer (L7-4, PhilipsHealthcare, Andover, MA). For the ARF push, a ^^ / ^^ ൌ 1.25 was used with a 400 ^^^^ at4.09 ^^^^^^. The detection was performed using plane wave compounding ^െ4∘, 0∘,^4∘^ witha pulse repetition period of 80 ^^^^ between each plane wave transmission. AfterDocket: 221407-2140 compounding, the effective pulse repetition period was 240 ^^^^ and an effective frame rate of 4166 ^^^^. The ^^ component of the particle velocity for both ex vivo resulting from ARF push was estimated from the acquired in-phase / quadrature (IQ) data using an autocorrelation method. Four SWE experiments were performed at each of three apparently homogenous sites of the liver. A total of twelve datasets were generated considering left and right propagating waves. To reduce the noise, the particle velocity responses are averaged over a3 mm strip along the ^^ direction, centered at the focal depth of the ARF, resulting in the ^^ െ ^^data. The data has a resolution of ^^^^ ൌ 0.154 ^^^^, ^^^^ ൌ 0.24 ^^^^, ^^^^௫ ൌ 19.2 ^^^^ and^^^^௫ ൌ 23.5 ^^^^. The signal is further truncated by a window, 3 ^ ^^ ^ 13 ^^^^ and 3 ^ ^^ ^13 ^^^^, to highlight the main signal by eliminating the noise away from the signal and the nearfield uncertainties associated with unknown and potentially unfocused ARF. For the Fourier transform, we perform zero padding so that the number of points is 2^ସin bothspatial and temporal directions, leading to ^^^^ ൌ 0.254 ^^^^ and ^^^^ ൌ 2.49 ^^^^^^ / ^^.

[0061] For validation purposes, the liver was cored at four locations away from large blood vessels for independent characterization. The resulting 9.5 ^^^^ cylindrical samples were mechanically tested with hyper-frequency viscoelastic spectroscopy (Rheospectris C500+, Rheolution, Inc., Montreal, Quebec, Canada), to obtain the storage ^^^^^ and loss ^^^^^ moduli. In vivo Studies

[0062] To examine in vivo applicability of the TPM, analysis was performed on data from SWE experiments with a curved array transducer (C5-2v, Verasonics, Inc., Kirkland, WA) on in vivo human liver in a healthy individual. The data were collected based on a protocol approved by the Mayo Clinic Institutional Review Board and each subject provided written informed consent. The ^^ component of the particle velocity for in vivo resulting from ARF push was estimated in the same manner of ex vivo. The in vivo measured responseconsists of Δ^^ ൌ 0.24 ^^^^,Δ^^ ൌ 0.36 ^^^^, ^^^^௫ ൌ 65.9 ^^^^ and ^^^^௫ ൌ 32 ^^^^. Similar to theex vivo measurements, due to far field noise and ARF uncertainties, the response istruncated to a smaller ^^ െ ^^ window, 10 ^ ^^ ^ 20 ^^^^ and 2.52 ^ ^^ ^ 19.8 ^^^^, to eliminatenoise artifacts in the far field at early times. For the Fourier transform, zero padding was performed so that the number of points is 2^ସin both spatial and temporal directions, leadingto Δ^^ ൌ 0.17 ^^^^ and Δ^^ ൌ 1.60 ^^^^^^ / ^^.Results and Discussion

[0063] Comparison between Frequency Domain and Time Domain Forward Models: In this section, the results from frequency and time-domain forward models (2D versions) for the Kelvin-Voigt model are briefly compared. The simulation results are truncated as done for TPM inversion, and the two sets of peaks are compared in FIGS.4A-4C, which showsDocket: 221407-2140 twin peaks from time domain simulation (dashed lines) and frequency domain (solid lines) simulation. To shed light on the truncation effect on the peaks, three gradually reducingwindow sizes are considered: ^^^^௫ ൌ 90 ^^^^, ^^^^௫ ൌ 45 ^^^^ (FIG. 4A), ^^^^௫ ൌ 30 ^^^^,^^^^௫ ൌ 15 ^^^^ (FIG. 4B), and ^^^^௫ ൌ 15 ^^^^, ^^^^௫ ൌ 7.5 ^^^^ (FIG. 4C). While the peaks fromdifferent windows are different, they match between frequency and time-domain analyses verifying their equivalence. Given the efficiency of time-domain simulation (taking only 0.82 second compared to 5.5 seconds for frequency domain simulation), the use of time domain simulation is advocated whenever the Kelvin-Voigt model is used for inversion.

[0064] The plots also shed light on the effects of the windowing as discussed in the Methods section. At lower frequencies and wavenumbers, truncation introduces the artifact of peaks crossing each other, with crossing point occurring at larger frequncies and wavenumbers as the window is made smaller. At higher frequencies, the spread spuriously increases which is associated with the convolution effect with the Sinc function. Importantly, these artifacts have predictable and identical efffects, justifying the use of peak mathching as long as the windowing is done in the same way for both experimental and simulated data. In silico Verification

[0065] Noise-free Inversion with Kelvin-Voigt Model: FIG.5A illustrates an example of twin peaks from the data (dashed) along with inverted peaks (solid) and FIG.5B shows an example of an objective function for noise-free in silico data. FIG.5A shows the match between noise-free in silico peaks (obtained from 3D simulation) and inverted peaks(obtained from 2D simulation). The inverted parameters are ^^^ ൌ 2.04 ^^^^^^ and ^^ ൌ 0.51 ^^^^,which are close to the true parameters ^^^ ൌ 2.00 ^^^^^^ and ^^ ൌ 0.50 ^^^^. The objectivefunction defined in Equation (26) is plotted in FIG.5B, showing a bowl shape with a clear minimum. The eigenvalues of the Hessian of the objective function near the optimum parameter values were examined and they were found to be 2.07 and 15.3, which are essentially measures of the principal curvatures of the bowl. Instead of dealing with uncertainty quantification, focus on parameter identifiability and thus the ratio of large to small eigenvalues being small (and positive). Given that they are of the same order, both modulus and relaxation time are independently identifiable, which is consistent with visual observations. These observations illustrate the robustness and accuracy of TPM, as well as the underlying 2D approximation in the forward model.

[0066] Table 2 summarizes the inversion results for varying levels of relaxation time, illustrating the effectiveness of the TPM method even with a simpler 2D forward model. InTable 2 the relative errors are very low, except for the low viscosity case of ^^ ൌ 0.1 ^^^^. Thehigh relative error for ^^ ൌ 0.1 ^^^^ is attributed to the twin peaks not diverging significantlyenough in the chosen ^^ െ ^^ range, consistent with Equations (13) and (15) and theDocket: 221407-2140 discussion in the subsection “Frequency and Wavenumber Ranges”. TPM inversion wasrepeated with an expanded range of 200 ^ ^^ ^ 1000 ^^^^,  800 ^ ^^௫ ^ 3500 ^^ି^, whichresulted in a significant improvement in accuracy with 10% and 0.3% errors in ^^ and ^^^, respectively, confirming the hypothesis. Table 2. Inversion results for noise-free in silico data Actual Parameter Inverted Parameter Absolute Error Error Percentage (%) ^^^^^^^^^^^^ ^^^^^^^^ ^^^^^^^^^^^^ ^^^^^^^^ ^^^^^^^^^^^^ ^^^^^^^^ ^^^^ ^^2.00 0.7 1.97 0.68 0.03 0.02 1.5 2.85 2.00 1.0 2.03 0.96 0.03 0.04 1.5 4.00

[0067] Inversion using Synthetic Data with Noise: FIG.6A illustrates an example of Twin peaks from the data (dashed) and final inverted peaks (solid) and FIG.6B shows an example of an objective function for noisy in silico data. FIG.6A shows a typical match between noisy in silico peaks (medium noise, informed by visual comparisons with real experimental data), and inverted peaks. The corresponding objective function, shown in FIG. 6B, has a bowl shape with a clear minimum, indicating the robustness of TPM against noise.The errors in inversion of G 0 are shown in FIG. 6C (errors in elastic modulus for tenrealizations), indicating ^ 5% error for various levels of relaxation time (the mean andstandard deviations are shown for these errors with the thick lines representing the mean and the thin lines representing the standard deviation). The errors in viscosity are larger, approximately at 10%, for relaxation times less than 0.45 ^^^^ and approximately 5% for higher relaxation times (shown in FIG.6D, errors in viscosity). As expected, these errors are larger than those from noise-free data but are still acceptable.

[0068] Table 3 summarizes the results for different noise cases which informs us about the robustness of the TPM with wide range of noise where actual parameter being ^^^ൌ2 ^^^^^^ and ^^ ൌ 0.5 ^^^^. FIG. 7 provides visualization of the spatiotemporal representation ofthe noisy data (top row), corresponding objective functions (middle row) and peak matches for different noise levels (bottom row - dashed from data and solid is inverted). Columns represent (a) no noise, (b) low noise, (c) medium noise, and (d) high noise. Similar to FIG. 5B, the eigenvalues of the Hessian at the optimal point are computed and are given by (2.07, 15.3), (1.96, 13.6), (1.52, 10.3) and (0.6, 4.7) respectively, indicating that both shear modulus and relaxation time are independently identifiable. Note that, even for high noise,Docket: 221407-2140 the objective functions are bowl shaped with clear minima, the match is generally good, and the relative errors shown in Table 3 are acceptable. Table 3. Inversion results for noisy in silico data Noise Inverted Absolute Error Error Percentage level Parameter (%) ^^^^^^^^^^^^ ^^^^^^^^ ^^^^^^^^^^^^ ^^^^^^^^ ^^^^^^ Noise Free 2.04 0.51 0.04 0.01 2.0 2

[0069] inverting the noisy synthetic data obtained from spring-pot model, shown in FIG.8A, toestimate both the reference modulus ^^^ and the power ^^. FIG. 8A shows in silico x^ t datafrom spring-pot model (with noise). The match between in silico and inverted peaks for spring-pot model is shown in FIG.8B, with twin peaks obtained from data (dashed) andinverted peaks (solid). The inverted parameters, ^^^ ൌ 1.98 ^^^^^^and ^^ ൌ 0.407, are close tothe true parameters of ^^^ ൌ 2 ^^^^^^ and ^^ ൌ 0.4, i.e., with errors ^ 2%. The objective functiondefined in Equation (26) is plotted in FIG.8C which has a bowl shape with clear minimum, and the eigenvalues of the Hessian are (1.92, 26.2), illustrating the robustness of the inversion for both modulus ^^^and power ^^. Ex vivo Validation

[0070] FIGS.9A-9E illustrate ex vivo validation of TPM. FIG.9A shows an example of particle velocity data (for a single SWE sample); FIG.9B shows an example of truncated particle velocity data; FIG.9C shows an example of twin peaks for three SWE samples (dashed lines) along with inverted peaks (solid lines); FIG.9D shows an example of estimated moduli (solid lines) along with Rheospectris measurements (dashed lines); and FIG.9E shows an example of an objective function. Out of eleven data sets, the peaks from three are very closely clustered together, while the scatter within the remaining peaks is more significant. Given this, the three clustered peaks were inverted together, which are shown as the dashed lines in FIG.9C. The peaks are observed to be repeatable for thefrequency range of 70 ^ ^^ ^ 140 ^^^^ and the wavenumber range of 450 ^ ^^௫ ^ 900 ^^ି^;these ranges are thus used for inversion using TPM. The Kelvin-Voigt model was chosen for inversion and the inverted peaks are shown as solid lines in FIG.9C, illustrating a close match. The resulting storage and loss moduli are shown in FIG.9D, along with corresponding Rheospectris measurements. Interestingly, the Rheospectris measurements do not conform to the Kelvin-Voigt model, but the inverted moduli from TPM are stillDocket: 221407-2140 reasonably close to the Rheospectris measurements. The objective function is plotted in FIG.9E, which is bowl shaped; the eigenvalues of Hessian are (0.9, 57.9). Unlike in the in silico studies, more sensitivity to elastic modulus was observed compared to relaxation time (viscosity). Closer examination of Hessian can be performed along with analysis of the noise characteristics, leading to uncertainty quantification (UQ) of the two parameters. Notwithstanding this, given that the eigenvalues are not different by many orders of magnitude, both parameters are independently identifiable, and in fact a clear minimum is observed in the objective function in both directions.

[0071] The validation using the remaining SWE datasets has not been successful due to various reasons of the data quality and model mismatch. FIGS.10A-10B illustrate examples of problematic data sets for ex vivo acquisition. FIG.10A shows an example data set with noise artifacts rendering it unusable; FIG.10B shows an example data set with apparent wavefront bending, reflecting either heterogeneity or more complicated viscoelasticity; and FIG.10C shows an example of twin peaks for all the problematic data sets. Examination of seven of the remaining nine data sets contained either noise artifacts or wavefront bending. Noise artifacts, as shown in FIG.10A, significantly pollute the signal rendering it unusable in the 2D Fourier transform framework. Wavefront “bending” in space- time, shown in FIG.10B (compared to FIG.9A), indicates changing group velocity with ^^, which is an indicator of localized spatial heterogeneity. The hypothesis of heterogeneity is further supported by the variability of spatiotemporal patterns observed across various data sets (note that all these are obtained within the same liver at different locations); three examples are provided in FIGS.9A, 10B and 11A (these varying patterns result in a large scatter of twin-peaks for the seven data sets, presented in FIG.10C, compared to FIG.9C). At the end, the observation of spatial heterogeneity renders TPM inapplicable for these data sets as the method is based on the assumption of homogeneity.

[0072] The final two data sets are similar in quality as those in FIGS.9A and 9B, but different in spatiotemporal distribution (they are measured at the same location with different pushes and are repeatable). TPM was applied to these data sets and the results are presented in FIGS.11A-11E, which shows ex vivo inversion of the last two data sets (analogous to FIGS.9A-9E). FIG.11A shows an example of particle velocity data (for a single SWE sample); FIG.11B shows an example of truncated particle velocity data; FIG. 11C shows an example of twin peaks for three SWE samples (dashed lines) along with inverted peaks (solid lines); FIG.11D shows an example of estimated moduli (solid lines) along with Rheospectris measurements (dashed lines); and FIG.11E shows an example of an objective function. The peak matching, presented in FIG.11C, is not satisfactory, unlike the good match observed in FIG.9C. Moreover, the validation in FIG.11D is not good, and the objective function shape in FIG.11E does not have a clear, strong minimum. TheseDocket: 221407-2140 indicate the existence of model error either with respect to homogeneity or viscoelasticity. While heterogeneity cannot be handled by TPM, the method can be expanded to the more complicated three-parameter fractional Voigt model.

[0073] The limited success of TPM validation with ex vivo data may be attributed to the specimen quality. Whenever the SWE data are repeatable across different locations, i.e., in FIG.9C, TPM not only provides good accuracy but also is robust as evident by the bowl shape of the objective function in FIG.9E. The data were acquired previously for a different study, where no attention was paid to the homogeneity of the specimen, i.e., repeatability of the data across different locations. A more systematic and elaborate ex vivo experimentation would help increase the robustness of TPM validation. In vivo Application

[0074] FIGS.12A-12D illustrate an in vivo application. FIG.12A shows an example of particle velocity; FIG.12B shows an example of truncated x-t data; FIG.12C shows an example of twin peaks (SWE is dashed line and simulated is solid line); and FIG.12D shows an example of an objective function. Given the poor signal-to-noise ratio in the in vivo data(FIG. 12A), the response is truncated further in smaller ^^ െ ^^ window as shown in FIG. 12bfor better characterization of SWE twin peaks (dashed lines in FIG.12C). The inversion wasperformed for frequency range of 60 ^ ^^ ^ 170 ^^^^ and wavenumber range of 300 ^ ^^௫ ^800 ^^ି^selected based on observed repeatability and spread between the twin peaks; these windows are narrower than ex vivo validation due to higher noise in the in vivo data. The inversion resulted in a storage shear modulus of 1.45 kPa and a relaxation time of 0.63 ms, which are in line with a healthy human liver. While these numbers cannot be validated, the reliability and robustness of the inversion process can be examined through plotting the objective value ^^^^^in Equation (26) as a function of the inversion parameters. This is shown in FIG.12D, which has a bowl shape, with the eigenvalues of Hessian being 1.61 and 3.25, illustrating the robustness of inverting for both modulus and relaxation time. This robustness needs to be confirmed with more elaborate in vivo studies with multiple datasets and subjects.

[0075] In an attempt to provide reliable point measurements of not just elasticity but viscoelasticity of soft tissues using shear wave elastography (SWE), the so-called twin peak method (TPM) was developed. The TPM is based on the observation that viscosity results inspreading of the particle velocity in the frequency-wavenumber ^^^ െ ^^^ domain, which in turnresults in the separation of ^^^^^^ peaks obtained by stepping through ^^, and ^^^^^^ peaks obtained by stepping through ^^. Based on the premise that the peaks are typically less sensitive to noise compared to other amplitude-based measures, the TPM estimates the viscoelasticity parameters by matching these peaks from experimental measurements toDocket: 221407-2140 those from simulation. The TPM is shown to be effective through verification using in silico data and validation using ex vivo data. The method is also applied to in vivo data, where the examination of the objective function leads to the observation that the TPM is robust and can be used even with in vivo data with significant noise.

[0076] While largely focused on using a Kelvin-Voigt model for viscoelasticity, the developed approach is applicable to general viscoelasticity, although the underlying model needs to be appropriately parameterized for the sake of inversion (such an approach for the spring-pot model was illustrated, while more complicated models can be considered depending on the context). The method can be further refined to expand and automate the range of the frequencies and wavenumbers used for inversion. Finally, TPM provides point measurements, but it does so by matching the ^^^^^^ and ^^^^^^ peaks that assume homogeneity within the measurement line (plane). This has an averaging effect, and the associated approximation can be investigated when significant heterogeneities are involved; such cases may need imaging algorithms, where the TPM results can potentially serve as initial estimates. Overall, TPM can provide a useful tool to measure tissue viscoelasticity that is accurate and robust against measurement noise. CNN-Based Estimation of Bulk Viscoelasticity Using SWE

[0077] Shear Wave Elastography (SWE) is a powerful imaging modality that measures the mechanical response of tissue to excitations by the acoustic radiation force (ARF). By analyzing the resulting shear wave propagation, it is possible to estimate viscoelastic parameters such as shear modulus (^^) and relaxation time (^^), both of which provide important insights into the structural and pathological state of soft tissues. While elasticity has long been the focus of elastography-based diagnostics, viscosity is increasingly recognized as a complementary biomarker, improving sensitivity to conditions such as fibrosis, brain disease, and tumor development.

[0078] Conventional SWE-based inversion techniques, such as dispersion fitting, twin- peak methods (TPM), or full waveform inversion (FWI), typically assume either sharp or known spatial distribution of ARF, characterized by the spatial width (^^) of the force distribution. However, this assumption rarely holds in clinical practice, where system- dependent variability, patient anatomy, and push implementation can all alter the effective ARF width. This leads to the problem: Inversion techniques that are highly sensitive to ^^ often perform poorly when ARF parameters are uncertain, resulting in inaccurate viscoelastic estimates or misinterpretation of tissue state.

[0079] Different inversion methods show varying degrees of sensitivity to ARF width. For example, while TPM and dispersion methods emphasize peak tracking in the frequency- wavenumber domain and may tolerate certain forms of noise, they still use fundamentallyDocket: 221407-2140 small ARF width. Similarly, FWI attempts to reconstruct tissue properties by matching full wavefields and thus requires even more precise characteristics of the ARF.

[0080] CNNs as a Model-Free Alternative. Convolutional Neural Networks (CNNs) offer a flexible and data-driven alternative that avoids explicit reliance on ARF parameter assumptions. When trained across multiple ARF widths, CNNs can, in principle, learn to generalize across varying ^^, enabling parameter inversion even when the underlying push characteristics are unknown. However, relying on CNN predictions alone may not be enough, and without a physical constraint or selection mechanism, CNNs may misattribute features of the data to incorrect parameter combinations. This is especially risky when dealing with noise or ambiguous inputs.

[0081] To address this, a hybrid strategy is proposed that combines CNN’s flexibility with physics-based selection. The method trains CNNs on synthetic data from the ^^–^^ or ^^–^^domain generated over a range of ^^^, ^^,^^^ values. During testing, CNN predictions are notaccepted at face value. Instead, for each model (or each discrete ^^), a synthetic ^^–^^ or ^^–^^response can be regeneratedmodel with the predicted ^^^, ^^^ values. Thissynthetic plot is compared with the original measured data using the misfit norm in or ^^–^^ domains.

[0082] The model whose prediction yields the lowest loss is selected as the best estimate, thus indirectly identifying the correct ARF width and corresponding viscoelastic parameters. This approach allows the tissue properties to be robustly recovered even under unknown ARF conditions while also enabling evaluation of model accuracy and ^^ inference.

[0083] Theoretical Framework. Shear wave propagation is simulated using a forward model for homogeneous bulk tissue and utilizes a 2D scalar viscoelastic wave model governed by the Kelvin–Voigt constitutive law. The governing equation captures both elastic and viscous behavior and is given by: ^^ பమ௨ൌ ^^∇ଶ^^ ^ ^^∇ଶ ^ப௨^ ^ ^^^^^, ^^^, (29)where ^^^^^, ^^^ isshear modulus, ^^ ൌ ^^^^ is the shear viscosity (with ^^ denoting the relaxation time), and ^^^^^, ^^^is an acoustic radiation force (ARF) modeled as a spatial Gaussian profile with equal width ^^ in both the ^^- and ^^-directions, applied for a fixed duration.

[0084] Applying 2D spatial and temporal Fourier transforms simplify the equation to an algebraic form and the velocity response for each ^^^௫,^^௬^ takes the form, ^^^^^^,^^^^ఠ ^^^^^^^^,^^^, (30)Docket: 221407-2140where ^^ ൌ ^^^௫ଶ ^ ^^௬ଶ and ^^^ ൌ ^^^ / ^^ is the shear wave speed (^^௫ , ^^௬ are the wavenumbers in^^ and ^^ directions). To generate ^^–^^ plots, we compute the inverse 2D Fourier transform of^^^^^^, ^^^ at each timestep to obtain:^^^^^, ^^^ ൌ ℱି^^ ^^^^^^^, ^^^^. (31)Extract ^^^^^,^^ ൌ 0, ^^^ as atraining and evaluation.

[0085] For k–^^ plots, apply a 2D Fourier transform to the windowed ^^–^^ response: ^^^^^,^^^ ൌ ℱ௫,௧^^^^^^^^^^^^^^^^^^, ^^^^, (32)with Gaussianforward model is used for generating synthetic data and validating CNN outputs via regenerated ^^–^^ or ^^–^^ plots, ensuring physical consistency throughout the training and evaluation stages.

[0086] Data Generation. Training data are generated in MATLAB by simulating thousands of responses, evenly distributed across ARF push widths, using the following parameter ranges: ^Shear modulus ^^ ∈ ^1.5,2.5^ kPa^ Relaxation time ^^ ∈ ^0.1,1.0^ ms^ ARF push width ^^ ∈ ^0.10,0.25,0.50,0.75,0.90^ mmEach simulation yields an ^^–^^ or ^^–^^ velocity plot. To replicate realistic measurement conditions, Gaussian noise is added at a signal-to-noise ratio (SNR) of 20 dB. After noise addition, both ^^–^^ and ^^–^^ representations are normalized to unit energy. In addition to windowing mentioned above, the ^^–^^ plots are down-sampled by retaining every fourth data point in both frequency and wavenumber dimensions. This reduction step minimizes data size and enables efficient training. The resulting plots serve as input samples for CNN training or evaluation.

[0087] FIGS.13A and 13B display two of the generated responses using different push widths while holding viscoelastic parameters fixed. FIGS.13A and 13B illustrate the effect of ARF width ^^ on the shear wave velocity response in the space–time (^^–^^) domain. Bothsimulations use ^^ ൌ 1.94 kPa and ^^ ൌ 0.632 ms, with ^^ ൌ 0.25 mm in FIG. 13A and ^^ ൌ 0.90mm in FIG.13B. Despite possessing the same tissue properties, the wave shape and velocity values differ noticeably. Larger ^^ values produce broader, more dispersed wavefronts, highlighting the influence of ARF geometry on shear wave propagation. This highlights the importance of accounting for ^^ during inversion.Docket: 221407-2140

[0088] CNN Training and Architecture. Two model configurations were evaluated: (1) five separate CNNs, each trained exclusively on data from a single ARF width ^^, from which we later employ loss-based selection; and (2) a combined CNN trained on samples spanning all ^^ values. All models operate on normalized ^^–^^ or ^^–^^ representations as input.

[0089] The CNN architecture was developed through manual hyperparameter tuning to balance validation accuracy with training efficiency. FIG.14 illustrates the CNN architecture used for viscoelastic parameter estimation. Convolutional layers extract spatial-temporal features from ^^–^^ or ^^–^^ inputs, followed by dense layers for regression of ^^ and ^^. In the example of FIG.14, it comprises three convolutional layers with increasing filter depths (32,64, 128), each followed by ReLU activation and 2 ൈ 2 max pooling. These layersprogressively extract abstract spatial-temporal features from the generated plots. The resulting feature maps are flattened and passed through two fully connected layers (128 and 64 units), followed by a final dense layer that evaluates the predicted viscoelastic parameters ^^ and ^^, implemented using the Keras API on a TensorFlow backend.

[0090] All models are trained using the Adam optimizer, which provides adaptive learning rates and performs well in practice for small- to medium-scale regression tasks such as the current task. A mean squared error (MSE) loss function is used to penalize larger deviations more heavily, making it well-suited for continuous parameter estimation. Mean absolute error (MAE) is tracked as a secondary metric to provide a more interpretable measure of average prediction error.

[0091] To prevent overfitting, early stopping is employed with a patience of 10 epochs, allowing the model to stop training once validation performance plateaus. Training is performed for a maximum of 20 epochs with a batch size of 32, which was selected for convergence speed and gradient stability. A 70 / 30 balance between training and validation ensures sufficient data for learning while preserving a reliable holdout set for model selection, a common split in CNNs.

[0092] Forward Model-Based Selection Workflow. The central innovation of the method is not in the CNN model itself or the training as the procedures used are quite standard, but in how CNN predictions are evaluated and selected. During testing, CNN outputs are not accepted directly. Instead, an ^^–^^ response was regenerated using the forward model based on the predicted (^^, ^^), and the misfit computed between this synthetic response and the measured response. This process can be repeated for each trained CNN (for each candidate ^^). The misfit is calculated using the ^^ଶnorm between the synthetic and measured ^^–^^ or ^^–^^ fields. The final prediction is chosen based on the one whose reconstructed wavefield produces the lowest misfit. This hybrid approach helps ensure that the selected parameters are internally consistent and physically meaningful.Docket: 221407-2140 Results

[0093] The following results are based on models trained using ^^–^^ domain inputs only. All evaluations are conducted on held-out test data not seen during training.

[0094] Motivation for Model Selection via ^^–^^ Misfit. The sensitivity of CNN predictions to ARF width mismatch can be quantified. FIGS.15A and 15B present the mean absolute percent error (MAPE) confusion matrices in predicted ^^ and ^^ across all combinations of test ^^ (rows) and model ^^ (columns) values, respectively. Errors increase substantially under ^^ mismatch.

[0095] In both matrices, the lowest errors almost always occur when the model ^^ matches the test ^^ (diagonal entries). Off-diagonal entries, corresponding to ^^ mismatch, show substantial increases in error, especially for ^^, where MAPE exceeds 100% in several cases. These findings underscore the need for a robust selection mechanism to handle unknown or variable ARF conditions, as fixed-CNN approaches degrade significantly when applied under mismatched push configurations.

[0096] ARF Width Inference via Model Selection. To address this sensitivity, a model selection strategy can be employed based on ^^–^^ misfit. For each test sample, estimate viscoelastic parameters using each ^^-specific CNN, regenerate the predicted ^^–^^ response via the forward model, and select the ^^-specific CNN that minimizes the misfit.

[0097] FIG.16 shows a confusion matrix of selection accuracy. In FIG.16, the confusion matrix shows ARF width ^^ recovery accuracy via ^^–^^ misfit-based model selection. Each row corresponds to the true ^^ (from the test sample), and each column to the selected model ^^. Perfect diagonal structure reflects accurate inference. This validates that the misfit-based strategy enables reliable and implicit ARF width inference without needing explicit knowledge of the source.

[0098] Viscoelastic Parameter Recovery. Using the model selected via ^^–^^ misfit, the final accuracy in estimating the viscoelastic parameters ^^ and ^^ can be assessed. FIGS. 17A-17B and 18A-18B summarize prediction error across the entire test set. FIGS.17A and 17B illustrate examples of signed percent error in predicted viscoelastic parameters across the test domain. FIG.17A shows the error in shear modulus ^^ and FIG.17B shows the error in relaxation time ^^. Each dot represents a test sample, colored by signed percent error relative to the true value. FIGS.18A and 18B illustrate examples of absolute error in predicted viscoelastic parameters across the test domain. FIG.18A shows the absolute error in shear modulus ^^ (kPa) and FIG.18B shows the absolute error in relaxation time ^^ (ms). Most estimates fall within േ4% error, with average absolute percent errors of 1.47% for ^^ and 3.31% for ^^. The absolute errors average 0.0292 kPa and 0.0148 ms, respectively.Docket: 221407-2140 Prediction accuracy is highest near the center of the sampling space, where dense coverage and CNN generalization bias improve robustness.

[0099] Effectiveness of ^^–^^ misfit-based Selection. To demonstrate the practical value of the ^^–^^ misfit-based selection strategy, it can be compared against a baseline approach in which viscoelastic parameters are predicted directly from a single combined CNN trained across all ^^ values, without model-specific selection. Both approaches operate under the assumption that ARF width is unknown at initial inference.

[0100] FIGS.19A and 19B show the mean absolute percent error (MAPE) across test ^^ values for both methods (^^ and ^^, respectively). FIGS.19A and 19B illustrate the comparison of MAPE using a combined CNN with unknown ARF width (left bar) versus the proposed ^^–^^ loss-guided model selection approach (right bar). ^^–^^ loss guidance improves ^^ accuracy but results in higher error for ^^. For shear modulus ^^ in FIG.19A, the ^^–^^ misfit- guided approach outperforms the unknown-width baseline across nearly all ^^ values, withespecially large improvements at ^^ ൌ 0.90 mm, where error drops from 4.97% to 1.91%. Aslight under performance occurs only at ^^ ൌ 0.50 mm.

[0101] In contrast, for relaxation time ^^ in FIG.19B, the ^^–^^ loss-guided method yields consistently higher errors across all ^^, suggesting that while misfit-based model selection enhances ^^ estimation, it may somewhat reduce the accuracy in ^^ prediction. This reduction in accuracy might suggest that ^^–^^ plots from the forward model are less sensitive to viscosity differences than previously expected, or that the smaller training set size for each ^^-specific model (5,000 samples versus 25,000 for the combined model) may limit the network’s ability to delineate subtle viscosity effects.

[0102] While CNNs are powerful tools for learning complex mappings from velocity response images to physical parameters, they are sometimes not sufficient on their own for robust viscoelastic estimation when ARF push width is unknown. Direct CNN predictions are often susceptible to bias or error when trained and tested on mismatched ARF configurations (see FIGS.15A-15B). Without a misfit-based correction, the resulting viscoelastic estimates may be unreliable.

[0103] The proposed approach addresses this issue by introducing a model selection step that evaluates each prediction via forward simulation and compares the reconstructed ^^–^^ response to the observed one using ^^ଶmisfit. This misfit-guided selection process enables implicit identification of the correct ARF width (see FIG.16) and leads to more accurate and physically consistent viscoelastic estimates for shear modulus ^^ (see FIGS. 19A-19B). However, this approach yields worse performance for relaxation time ^^ compared to the combined CNN baseline, suggesting that model selection may introduce added variability for parameters less strongly expressed in the wavefield. Potential strategies toDocket: 221407-2140 address this limitation include employing a dedicated CNN model for ^^, using the combined model specifically for ^^ estimation (where it demonstrates superior performance), or implementing architectural modifications to the CNN models that enhance sensitivity to viscous effects while preserving the accuracy of ^^ estimation.

[0104] A novel framework for viscoelastic parameter estimation and ARF width inference has been introduced using convolutional neural networks (CNNs) trained on synthetic shear wave elastography (SWE) data. The method is designed to address a key limitation in current inversion techniques which is their sensitivity to the spatial profile of the ARF push and typically unknown in clinical settings. The core innovation lies in the integration of forward modeling as a physics-based selection mechanism. Rather than accepting CNN predictions directly, simulated responses can be regenerated using each candidate prediction and the one that minimizes the misfit with the measured data selected. This post-inference validation enables indirect yet effective inference of the ARF width ^^, and improves the accuracy of the recovered shear modulus ^^. Through this hybrid strategy, CNN-based inversion can remain robust and accurate even when source characteristics are not known a priori. Fat and Collagen Content Estimation

[0105] Inversion, either using the TPM method or using CNN approach, can also be performed to estimate not the viscoelastic parameters but more of the fat and collagen contents. This is performed by linking percent fat and collagen percent area (CPA) to the viscoelastic parameters either using empirical models or micromechanical simulation. This can simply be done by prefixing the forward models in with a model linking histological and mechanical properties, i.e., through a functional relation informed by empirical models or micromechanical simulation.

[0106] With reference to FIG.20, shown is a schematic block diagram of a computing (or processing) device 1000 that can be utilized for simulation, measurement, and imaging of viscoelasticity, fat and collagen content of soft tissues or other applications using the described techniques. In some embodiments, among others, the computing device 1000 may represent a mobile device (e.g., a smartphone, tablet, computer, etc.) or other processing device. Each computing device 1000 includes processing circuitry comprising at least one processor circuit, for example, having a processor 1003 and a memory 1006, both of which are coupled to a local interface 1009. To this end, each computing device 1000 may comprise, for example, at least one server computer or like device. The local interface 1009 may comprise, for example, a data bus with an accompanying address / control bus or other bus structure as can be appreciated.

[0107] In some embodiments, the computing (or processing) device 1000 can include one or more network interfaces 1012. The network interface 1012 may comprise, for example, a wireless transmitter, a wireless transceiver, and a wireless receiver. AsDocket: 221407-2140 discussed above, the network interface 1012 can communicate to a remote computing device using a Bluetooth protocol. As one skilled in the art can appreciate, other wireless protocols may be used in the various embodiments of the present disclosure.

[0108] Stored in the memory 1006 are both data and several components that are executable by the processor 1003. In particular, stored in the memory 1006 and executable by the processor 1003 are viscoelasticity, fat and collagen content process program 1015, application program 1018, and potentially other applications. Also stored in the memory 1006 may be a data store 1021 and other data. In addition, an operating system may be stored in the memory 1006 and executable by the processor 1003.

[0109] It is understood that there may be other applications that are stored in the memory 1006 and are executable by the processor 1003 as can be appreciated. Where any component discussed herein is implemented in the form of software, any one of a number of programming languages may be employed such as, for example, C, C++, C#, Objective C, Java®, JavaScript®, Perl, PHP, Visual Basic®, Python®, Ruby, Flash®, or other programming languages.

[0110] A number of software components are stored in the memory 1006 and are executable by the processor 1003. In this respect, the term "executable" means a program file that is in a form that can ultimately be run by the processor 1003. Examples of executable programs may be, for example, a compiled program that can be translated into machine code in a format that can be loaded into a random access portion of the memory 1006 and run by the processor 1003, source code that may be expressed in proper format such as object code that is capable of being loaded into a random access portion of the memory 1006 and executed by the processor 1003, or source code that may be interpreted by another executable program to generate instructions in a random access portion of the memory 1006 to be executed by the processor 1003, etc. An executable program may be stored in any portion or component of the memory 1006 including, for example, random access memory (RAM), read-only memory (ROM), hard drive, solid-state drive, USB flash drive, memory card, optical disc such as compact disc (CD) or digital versatile disc (DVD), floppy disk, magnetic tape, or other memory components.

[0111] The memory 1006 is defined herein as including both volatile and nonvolatile memory and data storage components. Volatile components are those that do not retain data values upon loss of power. Nonvolatile components are those that retain data upon a loss of power. Thus, the memory 1006 may comprise, for example, random access memory (RAM), read-only memory (ROM), hard disk drives, solid-state drives, USB flash drives, memory cards accessed via a memory card reader, floppy disks accessed via an associated floppy disk drive, optical discs accessed via an optical disc drive, magnetic tapes accessed via an appropriate tape drive, and / or other memory components, or a combination of any twoDocket: 221407-2140 or more of these memory components. In addition, the RAM may comprise, for example, static random access memory (SRAM), dynamic random access memory (DRAM), or magnetic random access memory (MRAM) and other such devices. The ROM may comprise, for example, a programmable read-only memory (PROM), an erasable programmable read-only memory (EPROM), an electrically erasable programmable read- only memory (EEPROM), or other like memory device.

[0112] Also, the processor 1003 may represent multiple processors 1003 and / or multiple processor cores and the memory 1006 may represent multiple memories 1006 that operate in parallel processing circuits, respectively. In such a case, the local interface 1009 may be an appropriate network that facilitates communication between any two of the multiple processors 1003, between any processor 1003 and any of the memories 1006, or between any two of the memories 1006, etc. The local interface 1009 may comprise additional systems designed to coordinate this communication, including, for example, performing load balancing. The processor 1003 may be of electrical or of some other available construction.

[0113] Although the viscoelasticity, fat and collagen content process program 1015 and the application program 1018, and other various systems described herein may be embodied in software or code executed by general purpose hardware as discussed above, as an alternative the same may also be embodied in dedicated hardware or a combination of software / general purpose hardware and dedicated hardware. If embodied in dedicated hardware, each can be implemented as a circuit or state machine that employs any one of or a combination of a number of technologies. These technologies may include, but are not limited to, discrete logic circuits having logic gates for implementing various logic functions upon an application of one or more data signals, application specific integrated circuits (ASICs) having appropriate logic gates, field-programmable gate arrays (FPGAs), or other components, etc. Such technologies are generally well known by those skilled in the art and, consequently, are not described in detail herein.

[0114] Also, any logic or application described herein, including the viscoelasticity, fat and collagen content process program 1015 and the application program 1018, that comprises software or code can be embodied in any non-transitory computer-readable medium for use by or in connection with an instruction execution system such as, for example, a processor 1003 in a computer system or other system. In this sense, the logic may comprise, for example, statements including instructions and declarations that can be fetched from the computer-readable medium and executed by the instruction execution system. In the context of the present disclosure, a "computer-readable medium" can be any medium that can contain, store, or maintain the logic or application described herein for use by or in connection with the instruction execution system.Docket: 221407-2140

[0115] The computer-readable medium can comprise any one of many physical media such as, for example, magnetic, optical, or semiconductor media. More specific examples of a suitable computer-readable medium would include, but are not limited to, magnetic tapes, magnetic floppy diskettes, magnetic hard drives, memory cards, solid-state drives, USB flash drives, or optical discs. Also, the computer-readable medium may be a random access memory (RAM) including, for example, static random access memory (SRAM) and dynamic random access memory (DRAM), or magnetic random access memory (MRAM). In addition, the computer-readable medium may be a read-only memory (ROM), a programmable read- only memory (PROM), an erasable programmable read-only memory (EPROM), an electrically erasable programmable read-only memory (EEPROM), or other type of memory device.

[0116] Further, any logic or application described herein, including the viscoelasticity, fat and collagen content process program 1015 and the application program 1018, may be implemented and structured in a variety of ways. For example, one or more applications described may be implemented as modules or components of a single application. For example, separate applications can be executed for the viscoelasticity estimation workflows. Further, one or more applications described herein may be executed in shared or separate computing devices or a combination thereof. For example, a plurality of the applications described herein may execute in the same computing device 1300, or in multiple computing devices in the same computing environment. Additionally, it is understood that terms such as “application,” “service,” “system,” “engine,” “module,” and so on may be interchangeable and are not intended to be limiting.

[0117] It should be emphasized that the above-described embodiments of the present disclosure are merely possible examples of implementations set forth for a clear understanding of the principles of the disclosure. Many variations and modifications may be made to the above-described embodiment(s) without departing substantially from the spirit and principles of the disclosure. All such modifications and variations are intended to be included herein within the scope of this disclosure and protected by the following claims.

[0118] Although specific terms are employed herein, they are used in a generic and descriptive sense only and not for purposes of limitation.

[0119] As will be apparent to those of skill in the art upon reading this disclosure, each of the individual embodiments described and illustrated herein has discrete components and features which may be readily separated from or combined with the features of any of the other several embodiments without departing from the scope or spirit of the present disclosure.

[0120] Any recited method can be carried out in the order of events recited or in any other order that is logically possible. That is, unless otherwise expressly stated, it is in no way intended that any method or aspect set forth herein be construed as requiring that its steps beDocket: 221407-2140 performed in a specific order. Accordingly, where a method claim does not specifically state in the claims or descriptions that the steps are to be limited to a specific order, it is in no way intended that an order be inferred, in any respect. This holds for any possible non-express basis for interpretation, including matters of logic with respect to arrangement of steps or operational flow, plain meaning derived from grammatical organization or punctuation, or the number or type of aspects described in the specification.

[0121] While aspects of the present disclosure can be described and claimed in a particular statutory class, such as the system statutory class, this is for convenience only and one of skill in the art will understand that each aspect of the present disclosure can be described and claimed in any statutory class.

[0122] It is also to be understood that the terminology used herein is for the purpose of describing particular aspects only and is not intended to be limiting. Unless defined otherwise, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which the disclosed compositions and methods belong. It will be further understood that terms, such as those defined in commonly used dictionaries, should be interpreted as having a meaning that is consistent with their meaning in the context of the specification and relevant art and should not be interpreted in an idealized or overly formal sense unless expressly defined herein.

[0123] Prior to describing the various aspects of the present disclosure, the following definitions are provided and should be used unless otherwise indicated. Additional terms may be defined elsewhere in the present disclosure.

[0124] As used herein, “comprising” is to be interpreted as specifying the presence of the stated features, integers, steps, or components as referred to, but does not preclude the presence or addition of one or more features, integers, steps, or components, or groups thereof. Moreover, each of the terms “by”, “comprising,” “comprises”, “comprised of,” “including,” “includes,” “included,” “involving,” “involves,” “involved,” and “such as” are used in their open, non-limiting sense and may be used interchangeably. Further, the term “comprising” is intended to include examples and aspects encompassed by the terms “consisting essentially of” and “consisting of.” Similarly, the term “consisting essentially of” is intended to include examples encompassed by the term “consisting of.

[0125] As used in the specification and the appended claims, the singular forms “a,” “an” and “the” include plural referents unless the context clearly dictates otherwise.

[0126] It should be noted that ratios, concentrations, amounts, and other numerical data can be expressed herein in a range format. It will be further understood that the endpoints of each of the ranges are significant both in relation to the other endpoint, and independently of the other endpoint. It is also understood that there are a number of values disclosed herein, and that each value is also herein disclosed as “about” that particular value in addition to theDocket: 221407-2140 value itself. For example, if the value “10” is disclosed, then “about 10” is also disclosed. Ranges can be expressed herein as from “about” one particular value, and / or to “about” another particular value. Similarly, when values are expressed as approximations, by use of the antecedent “about,” it will be understood that the particular value forms a further aspect. For example, if the value “about 10” is disclosed, then “10” is also disclosed.

[0127] When a range is expressed, a further aspect includes from the one particular value and / or to the other particular value. For example, where the stated range includes one or both of the limits, ranges excluding either or both of those included limits are also included in the disclosure, e.g. the phrase “x to y” includes the range from ‘x’ to ‘y’ as well as the range greater than ‘x’ and less than ‘y’. The range can also be expressed as an upper limit, e.g. ‘about x, y, z, or less’ and should be interpreted to include the specific ranges of ‘about x’, ‘about y’, and ‘about z’ as well as the ranges of ‘less than x’, less than y’, and ‘less than z’. Likewise, the phrase ‘about x, y, z, or greater’ should be interpreted to include the specific ranges of ‘about x’, ‘about y’, and ‘about z’ as well as the ranges of ‘greater than x’, greater than y’, and ‘greater than z’. In addition, the phrase “about ‘x’ to ‘y’”, where ‘x’ and ‘y’ are numerical values, includes “about ‘x’ to about ‘y’”.

[0128] It is to be understood that such a range format is used for convenience and brevity, and thus, should be interpreted in a flexible manner to include not only the numerical values explicitly recited as the limits of the range, but also to include all the individual numerical values or sub-ranges encompassed within that range as if each numerical value and sub-range is explicitly recited. To illustrate, a numerical range of “about 0.1% to 5%” should be interpreted to include not only the explicitly recited values of about 0.1% to about 5%, but also include individual values (e.g., about 1%, about 2%, about 3%, and about 4%) and the sub-ranges (e.g., about 0.5% to about 1.1%; about 0.5% to about 2.4%; about 0.5% to about 3.2%, and about 0.5% to about 4.4%, and other possible sub-ranges) within the indicated range.

[0129] As used herein, the terms “about,” “approximate,” “at or about,” and “substantially” mean that the amount or value in question can be the exact value or a value that provides equivalent results or effects as recited in the claims or taught herein. That is, it is understood that amounts, sizes, formulations, parameters, and other quantities and characteristics are not and need not be exact, but may be approximate and / or larger or smaller, as desired, reflecting tolerances, conversion factors, rounding off, measurement error and the like, and other factors known to those of skill in the art such that equivalent results or effects are obtained. In some circumstances, the value that provides equivalent results or effects cannot be reasonably determined. In such cases, it is generally understood, as used herein, that “about” and “at or about” mean the nominal value indicated ±10% variation unless otherwise indicated or inferred. In general, an amount, size, formulation, parameter or other quantity or characteristic is “about,” “approximate,” or “at or about” whether or not expressly stated to be such. It isDocket: 221407-2140 understood that where “about,” “approximate,” or “at or about” is used before a quantitative value, the parameter also includes the specific quantitative value itself, unless specifically stated otherwise.

[0130] As used herein, the terms “optional” or “optionally” means that the subsequently described event or circumstance can or cannot occur, and that the description includes instances where said event or circumstance occurs and instances where it does not.

[0131] Unless otherwise specified, temperatures referred to herein are based on atmospheric pressure (i.e., one atmosphere).

Claims

Docket: 221407-2140 CLAIMS Therefore, at least the following is claimed:

1. A method for evaluation of soft tissue, comprising: obtaining measured particle velocities from generated mechanical waves in tissue of a subject; determining viscoelastic parameters of soft tissue at a point or in a region of interest (ROI) from wave propagation characteristics derived from the measured particle velocities; and determining fat or collagen content at the point or in the ROI from the wave propagation characteristics.

2. The method of claim 1, wherein the viscoelastic parameters are determined using a twin peaks method by minimizing mismatch between the measured ^^^^^^ and ^^^^^^ peaks and the predicted ^^^^^^ and ^^^^^^ peaks.

3. The method of claim 2, wherein the predicted ^^^^^^ and ^^^^^^ peaks are obtained as analytical relations.

4. The method of any of claims 2 and 3, wherein the predicted ^^^^^^ and ^^^^^^ peaks are generated from particle velocities obtained from simulating SWE experiments.

5. The method of any of claims 2-4, wherein a complex-valued viscoelastic modulus is estimated directly as an arbitrary function of frequency.

6. The method of any of claims 2-5, wherein viscoelasticity is parametrized in a functional form, with the viscoelastic parameters estimated through inverse optimization.

7. The method of claim 6, wherein the parametrization captures direct mechanical behavior comprising elasticity and viscosity.

8. The method in any of claims 2-7, wherein viscoelasticity is parametrized in a functional form, where the viscoelasticity is linked to the fat or collagen content, which is directly estimated through inverse optimization.Docket: 221407-2140 9. The method in any of claims 1-8, where spatial maps of the fat or collagen content are estimated through inverse optimization and relationships linking fat and collagen contents to mechanical properties.

10. The method of claim 1, wherein the viscoelastic parameters or the fat or collagen content are determined using a convolutional neural network (CNN) based at least in part upon the measured particle velocities, the CNN trained over multiple acoustic radiation force (ARF) widths.

11. The method of claim 10, wherein the CNN is selected from a plurality of trained CNNs based upon a minimum misfit between synthetic and measured responses.

12. The method of claim 11, wherein the training is performed either by a single CNN model across multiple ARF widths or by multiple CNN models, with one CNN model associated with each ARF width.

13. A system for evaluation of soft tissue, comprising: an ultrasound scanner configured for shear wave elastography (SWE); and a computing device comprising a processor and memory, the computing device configured to at least: obtain measured particle velocities from generated mechanical waves in tissue of a subject; determine viscoelastic parameters of soft tissue at a point or in a region of interest (ROI) from wave propagation characteristics derived from the measured particle velocities; and determine fat or collagen content at the point or in the ROI from the wave propagation characteristics.

14. The system of claim 13, wherein the viscoelastic parameters are determined using a twin peaks method by minimizing mismatch between the measured ^^^^^^ and ^^^^^^ peaks and the predicted ^^^^^^ and ^^^^^^ peaks.

15. The system of claim 14, wherein the predicted ^^^^^^ and ^^^^^^ peaks are obtained as analytical relations.

16. The system of any of claims 14 and 15, wherein the predicted ^^^^^^ and ^^^^^^ peaks are generated from particle velocities obtained from simulating SWE experiments.Docket: 221407-2140 17. The system in any of claims 13-16, where spatial maps of the fat or collagen content are estimated through inverse optimization and relationships linking fat and collagen contents to mechanical properties.

18. The system of claim 13, wherein the viscoelastic parameters or the fat or collagen content are determined using a convolutional neural network (CNN) based at least in part upon the measured particle velocities, the CNN trained over multiple acoustic radiation force (ARF) widths.

19. The system of claim 18, wherein the CNN is selected from a plurality of trained CNNs based upon a minimum misfit between synthetic and measured responses.

20. The system of claim 19, wherein the training is performed either by a single CNN model across multiple ARF widths or by multiple CNN models, with one CNN model associate with each ARF width.

Images