{'fraction_non_alphanumeric': 0.04463427052101601, 'fraction_numerical': 0.016744922781014877, 'mean_word_length': 4.77130264446621, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonlinear Characterization of Tissue Viscoelasticity with Acoustoelastic Attenuation of Shear-Waves

Shear-wave elastography (SWE) measures shear-wave speed (SWS), which is related to the underlying shear modulus of soft tissue. SWE methods generally assume that soft tissue viscoelasticity is independent of mechanical loading, however, soft tissues are known to have viscoelasticity changing nonlinearly with pre-compression. Hence, characterization by SWS alone is insufficient, where nonlinear properties can be seen as confounders but may also be utilized as additional bio-markers. The viscoelastic nature of a medium is fully characterized by its storage and loss moduli, which are related to SWS and shear-wave attenuation (SWA). In this work, we study SWA characteristics as a function of applied strain to measure nonlinear viscoelastic parameters in soft tissues. For this purpose, we apply incremental quasi-static compression on the samples while measuring SWS and SWA, from which we derive storage and loss moduli to estimate nonlinear viscoelastic parameters as a function of applied strain using acoustoelasticity (AE) theory. Results from tissue-mimicking phantoms with varying oil percentages and ex-vivo porcine liver experiments demonstrate the feasibility of the proposed approach. In both these experiments, SWA was observed to decrease with applied strain. For 10 % compression in ex-vivo livers, shear-wave attenuation decreased on average by 28 % (93 Np/m), while SWS increased on average by 20 % (0.26 m/s).


I. INTRODUCTION
Shear-wave elastography (SWE) is a noninvasive imaging technique that maps shear-wave speed (SWS) in a tissue. Shear-waves are induced by acoustic radiation force and their propagation is observed using ultrasound imaging to capture lateral shear-wave travel speed [1], which is related to the underlying shear modulus of tissue. SWS measurements have been used in many clinical applications including diagnosis of diseases in liver, breast, and kidney [2], [3], [4]. SWS however has shown poor specificity in breast cancer diagnosis [5] and in the early-stage diagnosis of nonalcoholic fatty liver disease (NAFLD). For the latter, a major reason is the fat accumulation with steatosis, which may cause a reduction in SWS, acting as confounder for the expected SWS increase with simultaneously advancing fibrosis [6], [7].
For the further understanding and diagnosis of pathological changes in soft tissues, additional biomechanical markers would be greatly beneficial. Although SWE assumes that soft tissues are linear and elastic, they are inherently nonlinear and viscoelastic [8]. Several bio-markers have been proposed to characterize tissue nonlinearity and viscoelasticity to complement SWS measurements. For instance, in [9]  SWS (i.e., variation with frequency) was shown to correlate with the degree of viscosity induced by steatosis. In general, storage G and loss G moduli can together fully characterize viscoelastic behaviour of a medium. These moduli can be derived from SWS and shear-wave attenuation (SWA) (i.e., decay of shear-wave amplitude with distance) [10]. In [11], it was found for magnetic resonance elastography (MRE) of human liver tumors that the loss modulus provides better differential diagnosis between benign and malignant tumors than the storage modulus. With this motivation, several studies have since been conducted to estimate SWA accurately [12], [13], [14], [15]. In [12], SWA was measured at a particular frequency from its full width at half-maximum (FWHM) in kspace from the 2D-FFT of shear-wave propagation vs time profiles. Dependency of FWHM on length of spatial observation window (i.e. the segment length of the propagation-distance) used to compute 2D-FFT was investigated and accounted for with a model in [16]. SWA as the rate of shear-wave amplitude decay based on a cylindrical wavefront assumption was studied in [13], [17]. In [14], a frequency-shift (FS) method was introduced to measure SWA by fitting a model to the amplitude spectrum of shear-wave displacements at different propagation distances. Similarly, a two-point frequency-shift (2p-FS) method was developed in [15] by fitting a model to the amplitude spectrum of shear-wave displacements at any two propagation distances.
Nonlinear characteristics of soft tissues using SWE have been studied by measuring the change in SWS values at various incremental stresses using the theory of acoustoelasticity (AE) in quasi-incompressible media [18], [19], [20], [21], [22]. In [18], [19] nonlinear properties of liver samples were investigated. In [20], it was shown that benign and malignant lesions in the breast exhibited different nonlinear shear modulus. In [21], changes in the nonlinear shear modulus were studied in bladder samples before and after treatment with formalin. In [22], the effect of compression area and nonlinear characteristics of kidney samples was investigated for progressive and regressive compression. In [23] and [24], the Demiray-Fung hyperelastic model was employed to quantify nonlinear properties of ex-vivo porcine liver tissue and ex-vivo brain tissue, respectively. Frequency variations of the nonlinear shear modulus were studied recently in [25] using SWS dispersion characteristics in ex-vivo porcine liver samples. Note that in all of the above mentioned nonlinear characterization studies, the shear storage modulus is approximated as the shear modulus by neglecting the attenuation term, under the assumption that soft tissues are elastic. The storage modulus with the attenuation term, the loss modulus and SWA characteristics as a function of applied stress/strain have so far not been investigated, which is essential in order to fully quantify the nonlinear viscoelastic characteristics of soft-tissues under the influence of stress (compression). This is highly relevant given that the loss modulus is potentially superior to the storage modulus in differential diagnosis of benign and malignant tumors in the human liver [11]. Furthermore, SWA is has shown to be effective in segregation of transplanted livers with acute rejection vs. no rejection [12]. Such a characterization would not only yield additional bio-markers for diagnosis, but would also help in taking into account any potential confounding effects that may be introduced, e.g. by compression, during the measurement of SWS and SWA.
In this paper, based on AE theory and using acoustic radiation force (ARF) induced shear-waves, we study SWA changes with applied strain levels for characterizing tissue nonlinear viscoelastic properties, in particular storage and loss moduli.

A. Data acquisition
For this study, we employ shear-waves induced via ARF using the supersonic shear-wave imaging technique [26], which generates a cylindrical shear-wave front. We use three subsequent high-intensity ARF pushes of 200 µs duration at three axially separated foci with a separation of 5 mm in depth. Shear-wave propagation is then tracked utilizing ultrafast ultrasound imaging at 10 k frames-per-second and coherent compounding of three angled plane waves at (-8°,0°,8°), using a moving average filter of last three frames. This acquisition sequence was programmed in a research ultrasound machine (Verasonics, Seattle, WA, USA) with a 128-element lineararray transducer (Philips, ATL L7-4) operated at 5 MHz center frequency.
Shear-wave acoustoelasticity experiments were conducted by applying compression using a motorized three-axis linear stage in 1% increments up to 10% strain. The ultrasound probe, attached to a positioning stage, was used for inducing shear-waves via ARF and imaging the resulting shear-waves. To ensure uni-axial compression on phantoms, a custom designed 3D printed compression plate of dimensions 100 mm × 100 mm (with a small window for imaging through) was attached on the front of the ultrasound probe, as shown in Fig. 1a. From the shear-wave in-phase and quadrature (IQ) signal acquired at each compression step, the shear-wave propagation is tracked using a 2D Loupas autocorrelation method [27] after directionally filtering the left and right propagating waves [28] along the depth axis z, leading to displacement profiles u(r, t), where r is the lateral distance to the ARF source, and t the time.

B. Shear-wave propagation in viscoelastic medium
ARF-induced shear-wave in a viscoelastic medium is commonly modeled as a cylindrical wave originating from an infinite source line [29], defined in the frequency domain as: whereû(r, ω) is the Fourier transform of the shear wave displacement field u(r, t) polarized along depth (i.e. z-axis), ω is the angular frequency,k(ω) is the complex wave number, c p (ω) is the phase velocity, and α(ω) is the frequency dependent attenuation. As seen in Eq. (1a), the shear-wave amplitude decreases as a function of the propagation distance r due to medium viscous absorption and diffraction. Indeed, the shear-wave phase velocity and attenuation can be derived based on the Fourier domain spectrum in Eq. (1a). Taking the natural logarithm of the amplitude |û(r, ω)| and separating the phase terms to equate to phase angle û(r, ω) yield, respectively: Where A(ω) is a constant dependent on the frequency ω, and θ 0 is the phase angle at the wave source (ARF push point).
For a linear and isotropic viscoelastic medium, the complex shear modulusĜ(ω) is defined by the storage G (ω) and loss G (ω) moduli and is related to the complex shear-wave numberk through where ρ represents the tissue density. Substituting Eq. (3) in Eq. (1b) allows to derive the shear storage and loss moduli as In conventional SWE, shear storage modulus is often referred simply as shear modulus and it is commonly estimated from Eq. (4a) by neglecting shear-wave attenuation α(ω) as

C. Acoustoelasticity in viscoelastic medium
The influence of uniaxial compression on the speed of acoustic waves in elastic, isotropic, homogenous and nearly compressible media are studied using the theory of AE. For a shear-wave polarized along the axis of deformation, Gennisson et al. [30] approximated a linear dependence of shear modulus on uniaxial strain as follows: where ρ is the tissue density, v g the shear-wave group speed, µ 0 the shear modulus at stress-free (resting) condition (i.e., = 0), and A the third-order elastic constant. However, in a viscoelastic medium, the shear-wave velocity becomes complex and frequency-dependent, and hence the shear modulus also becomes complex and frequencydependent. Substituting Eq. (3) in Eq. (6) results in whereμ 0 (ω) = µ 0 (ω) + iµ 0 (ω) is the complex shear modulus at stress-free condition, andÂ(ω) = A (ω) + iA (ω) is the 3rd-order complex nonlinearity coefficient, composed of an elastic A (ω) and viscous A (ω) nonlinearity components.
In this study, we use Eq. (7) to estimate a complex nonlinearity coefficient from multiple complex shear moduli measurements at different strains. A schematic of our experimental process can be seen in Fig. 2. In practice, at each sample compression step, the phase velocity c p (ω) and attenuation α(ω) are computed based on Eq. (2) using the Fourier transform of tracked shear-waves (the detailed procedure will be explained later in the section. II-F). From c p (ω) and α(ω), the complex shear modulusĜ(ω) is derived using Eq. (4). Finally, the linear model in Eq. (7) is fitted separately to the real and imaginary moduli to estimate the nonlinear elastic A (ω) and viscous A (ω) constants (see Fig. 2c). Note that even the estimation of the elastic component of the 3rd-order nonlinearity constant would be affected whether or not shearwave attenuation is taken into account. In [25], shear-wave attenuation was neglected by estimating an elastic constant A(ω) from G(ω). This is equivalent to substituting shear-wave phase velocity c p (ω) instead of shear wave group speed v g in Eq. (6). We compare this approach with A (ω) computed using actual storage modulus G (ω) based on our derivation  (7) to quantify total nonlinearity of the targeted medium. In summary, we compute A, A , A , and |Â(ω)| using following equations, respectively: At any given frequency ω, from multiple readings of under different compressions, right-hand-side parameters are first found as described earlier. We then solve for any lefthand-side nonlinearity parameter in least-squares sense, i.e. a line-fitting problem.

D. Tissue mimicking phantom experiments
Four tissue mimicking gelatin phantoms were prepared with increasing castor oil content of {0%, 5%, 10%, 15%} to model increasing viscosity similarly to the studies in [9], [17], herein we call these phantoms {#1..#4}. For the phantom preparation, we follow a procedure similar to [31]: We pre-heated distilled water to ≈70°C before mixing in predetermined percentages of gelatin (G2500, Sigma-Aldrich chemical, USA) and castor oil (259853, Sigma-Aldrich chemical, USA). Anionic surfactant (Pril Original, Germany) of concentration 4 cc/L was added to make oil droplets sufficiently small for castor oil to mix homogeneously with gelatin solution. After the mixture is cooled down to ≈50°C, 1% of cellulose (S5504, Sigma-Aldrich chemical, USA) was added to provide scattering, essential for shear-wave tracking. Phantom compositions are tabulated in Tab. I. The mixture was then let to cool down further and, while still in a semi-liquid state, was poured in a 100×40×50 mm 3 mould. Until pouring, the mixture was stirred at a low speed throughout the preparation process, in order to avoid oil separation from gelatin and to ensure a homogeneous scatterer distribution. The phantom mixtures were then refrigerated at 4°C for 18 h. Prior to data acquisition, the phantoms were removed from their moulds and kept at room temperature for 6 h for a uniform temperature distribution. Shear-wave data was acquired by placing the phantoms on a hard flat surface, while applying compression with the compression plate parallel to this surface, as shown in Fig. 1b.
To improve SNR, shear-wave displacement fields over 11 acquisitions at each compression step were averaged prior to Fig. 2. Schematic illustrating the overall procedure of determining the nonlinear viscoelastic parameters. Shear-wave data is acquired at each compression step by applying stress on the targeted medium placed on a rigid surface using an ultrasound transducer connected to a linear motor stage setup (a). Shear-wave speed and attenuation is computed over a range of frequencies at each compression step to first estimate storage (G ) and loss (G ) moduli with Eq. (4) and then extract nonlinear parameters using Eq. (7). For e.g., shear wave speed and attenuation at frequency 300 Hz with applied strain are plotted in (b) and the derived G and G are shown in (c) along with its linear-fit to estimate in nonlinear viscoelastic parameters A and A , respectively. further processing steps. Compression with 0.5 mm steps up to 5 mm were applied, corresponding to 1% increments with the given phantom thickness.
To take care of the bias, that may arise due to changes in experimental conditions, three cycles of compression were performed on each phantom with 15 mins waiting in-between the measurements, such that the phantoms could recover from applied deformation during previous data acquisition. During this waiting period, the phantom were kept in a plastic storage box to prevent drying.

E. Ex-vivo liver experiments
To test the feasibility of the proposed measurements on actual tissue and to study the mechanical characteristics of liver tissue, an ex-vivo study was performed on four porcine livers, fresh on the same day they were brought from the local slaughterhouse. The livers were placed in a flat container as shown in Fig. 1c, and shear-wave data were collected for compression steps of 0.3 mm up to 3 mm for sections of livers that are roughly 30 mm in thickness. Similarly to phantom experiments, 3 sets of measurements were collected for each liver with 20 mins waiting in-between measurements for the liver parenchyma to recover. The liver samples were kept in 0.9% saline water between the experiments to avoid dehydration (also to prevent the intrinsic water squeezing out of the porous liver parenchyma under its own weight). Similarly to [25], to improve SNR, displacements fields over 30 acquisitions were averaged prior to further processing steps.

F. Implementation
Using supersonic approach, we induced shear-waves displacement fields u(z, r, t) by three consecutive ARF pushes with 5 mm separation, as shown in Fig. 3a. To increase SNR, we averaged axial displacements from an axial region-ofinterest (black box in Fig. 3a). To ensure a planar wavefront, we used a 5 mm axial window around the center push. This then reduces the displacement field u(z, r, t) in the ROI to u(r, t) = ROI u(z, r, t).
Calculations were performed for various combinations of start and end propagation distances to derive the best estimate of c p (ω) and α(ω) (see Fig. 5a,c). For this purpose, while computing c p (ω) and α(ω), corresponding coefficient of determination (R-square) of fit was computed to evaluate goodness of the fit and plotted in Fig. 5b,d corresponding to estimates in Fig. 5a,c respectively. Only the estimates with R-square greater than 0.95 were considered for phantoms and greater than 0.90 for ex-vivo porcine livers for further processing. Note that this computational procedure is different from [13] and [32]. In [32] shear-wave attenuation values were computed by fixing first propagation distance and searching for the optimal end propagation distance with minimum standard error of fit, while in [13] full transducer aperture size [0-15 mm] was considered. The selection for start and endpoint is crucial, because near field shear-wave variations and faster decaying higher shearwave frequency components create aberrations in the near field and in the far shear-wave propagation distances respectively (see Fig. 4e). At the in-between propagation distances, the shear-wave amplitudes are uniform due to formation of full shear-wave displacement profiles. Considering these limitations on amplitude profiles, it is difficult to accurately estimate shear-wave attenuation and depending on the chosen start and end point of propagation distances, estimated attenuation values are either over or underestimated. Hence, computing c p (ω) and α(ω) values over varying start and end lateral distances and taking the average of estimates having goodness-of-fit   greater than specific threshold are expected to yield more accurate α(ω). This approach is similar to the FS and 2P-FS methods in [14] and [15], respectively. These limitations are also valid in case of the phase profiles, i.e. the c p (ω) derivation, even though less prominent compared to amplitude profiles, as can be seen in Fig. 4e. Hence, the estimation of c p (ω) is based on the same procedure.
The resulted values are marked in Fig. 4c for c p (ω) and in Fig. 4f for α(ω), the same process was repeated for a range of frequencies, for the SWE data acquired at different strain levels, and for different compression cycles. In the following step, storage G (ω) and loss G (ω) moduli were derived from averaged c p (ω) and α(ω) over different compression cycles using the relationship mentioned in Eq. (7). To fit the functions in Eq. (8) and derive the nonlinear viscoelastic parameters from G (ω) and G (ω), we used the MATLAB least-squares fitting function fit. In order to quantify the goodness of the fit, we have used the metric root-mean-square error (RMSE) similar to [22] and considered only the outcomes having RMSE<0.25 kPa.

A. Phantom results
Normalized amplitude spectra calculated at propagation distances r = {5.1, 6.6, 8.1, 9.6, 11.1, 12.6} mm at stressfree condition are shown for each phantom in Fig. 6. Since the viscosity increases with increasing oil percentage, the amplitude spectrum gets narrower from phantoms #1 to #4, with a meaningful SNR not left above 500 Hz in #4, whereas the peak amplitude not being reached until ≈250 Hz in #1. This common frequency range [250, 500] Hz is thus chosen to comparatively analyze all phantoms, following the similar convention in [33]. Fig. 7 shows our estimated frequency-dependent shear-wave phase velocity and attenuation values at stress-free condition, i.e. ≈ 0. The error bars correspond to the standard deviation of the respective estimates having R-square value greater than 0.95. At any given frequency it is observed that, with increased oil percentage (i.e. from phantom #1 to #4), the shear-wave phase velocity decreases considerably, whereas the attenuation increases. Indeed, phantoms #1-#4 respectively exhibit a shear-wave velocity dispersion of {0.01, 0.02, 0.03, 0.08} [m/s]/100 Hz, and a shear-wave attenuation dispersion of {20.5, 24.5, 43.0, 84.5} [Np/m]/100 Hz, computed from the profiles in Fig. 7 using a linear model. For all the phantoms, the estimated shear-wave phase velocity and attenuation, as well as the storage and loss moduli derived from the former two, are shown in Fig. 8 as a function of frequencies and applied strain. The maps look homogeneous due to the large colorscale ranging across the phantoms as well as the relatively lower nonlinearity of gelatin phantoms (compared to liver tissue presented later).
In order to visualize absolute changes of corresponding parameters as a function of applied strain, in Fig. 9 we present the shear-wave phase velocity, attenuation, storage and loss moduli as averaged with respect to frequency axis and offset by their nominal values, respectively, SWS 0 , SWA 0 , G 0 , and G 0 at stress-free state. It is observed that, having 10% strain applied in phantoms #1-#4, respectively, shear-wave phase velocity changes by +2.9%, +2.5%, +2.4%, and +0.06%, while attenuation changes by -27.7%, -16.1%, -10.0%, and -0.8%, indicating an increase in phase velocity whereas a decrease in attenuation with applied strain. As expected from the relatively small nonlinearity of these phantoms, shear storage G and loss G moduli as functions of applied strain resemble the trends of shear-wave velocity and attenuation in Fig. 9a,b respectively. As a biomechanical marker, the nonlinearity parameter at a selected frequency (range) can be used. Alternatively, we report herein frequency-averaged values of these parameter estimates, for further analysis and interpretation. Frequencyaveraged values together with their standard deviations are presented visually in Fig. 10 and also tabulated in Tab. II. It can be seen that the nonlinear elastic parameters A and A and nonlinear viscous parameter A are decreasing exponentially with increasing castor oil percentage.
The absolute value of the nonlinear elastic parameter A without considering the attenuation term is 5%, 6%, 6% and 13% lower than true nonlinear elastic parameter A measured considering the attenuation term for phantoms #1-#4, suggesting the importance of attenuation term while computing the shear storage modulus. The standard deviation of A and its percentages with respect to its mean values are ±5.3%, ±9.5%, ±11.5% and ±164.4% for phantoms #1-#4 respectively, suggests an increasing variation of A with respect to frequency as oil percentage increases. A similar trend can also be observed in nonlinear parameters A and A .

B. Ex-vivo liver experiments
The feasibility of our method was next tested on four exvivo porcine livers subjected to 10% maximum strain. The normalized amplitude spectra calculated at propagation distances r = {5. 1, 6.3, 7.5, 8.7, 9.9} mm at stress-free condition are shown for each liver in Fig. 11. Porcine liver shear-wave spectra is seen to be much narrower compared to that of phantoms, due to the higher intrinsic shear-wave attenuation. For a   comparative analysis between the livers, a common frequency range [100, 200] Hz was chosen, given the maximum peak and maximum sufficient signal ranges, as described earlier for the phantoms. For these liver samples, estimated shear-wave phase velocity and attenuation, as well as the storage and loss moduli derived from the former two, are shown in Fig. 13 as a function of frequencies and applied strain. These maps are observed to be much less uniform compared to those of phantoms in Fig.8, due to highly complex, nonlinear, and dispersive nature of actual tissues.
Similar to the phantom experiments, shear-wave phase velocity, attenuation, storage and loss moduli were averaged with respect to frequency axis and subtracted from its frequency averaged values SWS 0 , SWA 0 , G 0 and G 0 at stress-free state respectively, and plotted in Fig. 14 to quantify cumulative effects in frequency with respect to applied strain. It is observed that, shear-wave phase velocities exhibit +17.8%, +23.7%, +5.6% and +34.9% while attenuation exhibit -28.1%, -47.0%, -03.7% and -33.7% relative changes per 10% applied strain in the livers #1-#4, respectively. It is interesting to note that, in contrast to phantom experiments, for the livers the changes in storage modulus with respect to applied strain resembles the changes in attenuation while changes in loss modulus with respect to applied strain resembles the changes in phase velocity. Fig. 15 illustrates viscoelastic nonlinear parameters A , A and its magnitude |Â|, and nonlinear elastic parameter A. Fig. 9. Frequency-Averaged shear-wave phase velocity (a), attenuation (b), storage and loss moduli (c) with respect to frequencies on the corresponding maps in Fig. 8 and offset by their nominal values, respectively, SWS 0 , SWA 0 , G 0 , and G 0 at stress-free state are plotted as function of applied strain ( ), the error bars represent the standard deviation of these values. Fig. 10. Frequency-averaged elastic A and viscoelastic nonlinear parameters A , A , and its magnitude |Â| are extracted from phantoms #1-#4 that were subjected to compression of 0.5 mm steps for a total of 5 mm (≈ 10%). The errorbars indicate the standard deviation across frequencies.
These nonlinear mechanical parameters are averaged over the frequencies [100-200 Hz] and its standard deviations are presented as error bars. These values also tabulated in Tab. III. Similar to the observations in the phantom experiments, the absolute value of the nonlinear elastic parameter A without considering the attenuation term is +13.1%, +2.5%, -8.0%, and +62.0% higher than true nonlinear elastic parameter A measured considering the attenuation term for livers #1-#4. The standard deviation of A and its percentages with respect to its mean values are ±33.9%, ±24.0%, ±141.6%, and ±32.9% for livers #1-#4. This suggests a highly dispersive characteristic of actual tissue's non linear constants, similar to the observations in [25]. A related trend can also be observed in nonlinear parameters A and A .

IV. DISCUSSION
In the phantom study, an increasing castor oil percentage gave increased shear-wave velocity and attenuation dispersion profiles similar to the observations in [9]. As the oil percentage increases, speed-of-sound decreases and attenuation of the medium increases [34] which resulted in poor SNR of acquired shear-wave displacement fields. This is evident from the larger error bars of shear-wave phase velocity and attenuation estimates for the phantom #4 (gelatin 08% and castor oil 15%). Note that the error bars of shear-wave attenuation estimates are larger than the phase velocity estimates since the amplitude spectrum profiles (Fig. 4e) of shear-waves are generally noisier than the phase profiles (in Fig. 4b).
Shear storage modulus as a function of applied strain as shown in Fig. 9c suggests that nonlinearity decreases with increased castor oil percentage and is highest in the phantom with only gelatin (i.e. stiffer medium), which is in agreement with the findings in [35], where nonlinearity is reported to decrease with increasing oil percentage in tissue mimicking phantoms during mechanical tests.
It is observed that storage modulus (G ) was increased while loss modulus (G ) was decreased with applied strain as shown in Fig. 9c. Note that such G change with applied strain depends on the operating points of phase velocity and attenuation and their shifts with applied strain in a medium at a frequency. To demonstrate this, we plotted the theoretical storage (G ) and loss modulus (G ) values, respectively, in Fig. 16a,b over a range of different shear-wave phase velocity and attenuation values at 200 Hz frequency. Since shear-wave attenuation values decrease and phase-velocities increase as a function of applied strain, Fig. 16b suggests that a loss modulus (G ) direction of variation (decrease or increase) as a function of applied strain depends on the given shearwave attenuation and phase velocities in the phantom/tissue of interest. For instance, if the applied strain shifts tissue properties in the direction of yellow arrow in Fig. 16b, then G would increase; whereas a change in red arrow directions would yield a reduction in G .
For all the phantoms used herein, A > A is observed as in Fig. 10. Nevertheless, this observation cannot to be generalized and depends on the material parameter change as well. To   illustrate this, the difference between shear storage modulus ( Fig. 16b) and conventional shear modulus (Fig. 16c), i.e. G − G , is plotted in Fig. 16d. Note that the change in G (or G ) with applied strain is the parameter A (or A ) given Eq. (8). For any applied strain causing material changes along the yellow arrow in Fig. 16d, the difference G − G is seen to increases, i.e. change in G being larger than change in G , hence A > A ). The red arrow exemplifies a behaviour opposite to the one above. In Fig. 12b, attenuation is seen to increase linearly as a function of frequency, in agreement with the findings in [15], [14]. As seen in Fig. 14a, Fig. 8a,b. Note that in both the phantom and ex-vivo liver experiments the absolute change in shear-wave phase velocity is correlated with the absolute change in attenuation with applied strain. This characteristic may be due the displacement of extra- Fig. 13. Shear-wave phase velocity (1st row), attenuation (2nd row), storage modulus G (3rd row) and loss modulus G (bottom row) with respect to frequencies and applied strain levels for ex-vivo porcine livers #1 -#4. cellular fluids and reduced gaps in extracellular matrix upon compression. This might in turn lead to a higher tissue density, thereby increasing the shear-wave velocity and decreasing the attenuation. To study general applicability of this observation, future studies with more soft tissues are needed.
Stress-free shear moduli [1.05,1.55] kPa obtained in our exvivo liver experiments are comparable with the values reported in [36], where from multiple harmonic excitations, mean shear moduli of 2.2 kPa and 1.8 kPa are reported, respectively, using Kelvin-Voigt and Zener models. Average G to G ratio from our ex-vivo liver experiments is 1.1, suggesting that the loss and storage moduli are in the same order and both need to be  Fig. 13 and offset by their nominal values, respectively, SWS 0 , SWA 0 , G 0 , and G 0 at stress-free state are plotted as function of applied strain ( ), the error bars represent the standard deviation of these values. Fig. 15. Elastic A and viscoelastic nonlinear parameters A , A , and its magnitude |Â| are extracted from 4 ex-vivo porcine livers that were subjected to compression of 0.3 mm steps for a total of 3 mm (≈ 10%). The errorbars indicate the standard deviation across frequencies.
taken into account to fully characterize the liver or the soft tissues. Hence, loss modulus cannot simply be ignored and is important to consider for the AE effect to fully quantify the nonlinear viscoelastic characteristics of a soft tissue. Thirdorder nonlinear elastic parameter A of porcine livers were reported in [25] in the range of [25,60] kPa at 200 Hz. These findings are corroborated by the frequency-averaged A values over [100,200] Hz range reported in this work.
In addition to A and A , the magnitude of the thirdorder nonlinearity parameter, |Â|, reported in this study could potentially serve as a new bio-marker, along with its loss factor (tan δ = A /A ). For angle δ = 0°, the medium exhibits purely elastic nonlinearity, and for 90°purely viscous. Meanwhile, |Â| can capture and quantify the severity of either nonlinearity. For instance in Fig. 15, for liver #2, A < A , but for liver #4, vice versa A > A . Magnitude of overall nonlinearity |Â| is seen to be similar, while the δ changing drastically. Such reparametrization may then provide an alternative understanding, additional correlative evidence, and potentially simpler interpretation in the clinics.
All liver samples in this study were cut into rough rectangular boxes before measurements for standardization. Measured ex-vivo values are likely different from their in-vivo counterparts due to lack of perfusion and potential dehydration during the experimental process. Note as a limitation that uniaxial compression may not be straight-forward to apply in-vivo due to anatomical complexity and constraints around targeted tissue. Extremities and the breast, nevertheless, can be compressed with controlled handheld setups similar to [37], [38] to allow for the proposed nonlinear viscoelastic parameterization.

V. CONCLUSION
In this work, nonlinear characteristics of tissue viscoelasticity using acoustoelastic attenuation of shear-waves have been studied. We first computed and studied the characteristics of frequency-and strain-dependent velocity and attenuation val-ues of acoustic radiation force induced shear waves. Next, we derived storage and loss moduli from these values to compute nonlinear mechanical parameters by fitting parametric models of acoustoelastic theory in viscoelastic media under uniaxial compression. We have studied these with gelatin phantoms of different oil emulsions, similarly to literature aiming to emulate NFALD. Results indicate the following observations: (i) dispersion of shear-wave attenuation and velocity increased with higher oil percentage; (ii) variation of shear-wave velocity and attenuation with respect to compression decreased with higher oil percentage; and (iii) magnitudes of nonlinear viscoelastic parameters decreased with higher oil percentage.
A study with four ex-vivo porcine livers were carried out to test the feasibility of the proposed methods on actual tissue. Results reveal similar observations that the variation of shear-wave velocity is correlated with the variation of shear-wave attenuation as a function of applied strain. With both the phantom and ex-vivo studies, we show that using conventional shear modulus (i.e., omitting shear-wave attenuation) introduces considerable error in nonlinear parameter estimation. Based on our results, we recommend the use of storage modulus (computed using both shear wave velocity and attenuation) for computing a true nonlinear elastic parameter. Nonlinear biomechanical characterization of tissues can add complementary information to conventional SWE, helping in early diagnosis and staging of different diseases such as NAFLD. In future clinical studies our methods shall be compared with clinical gold-standard such as biopsies, to assess clinical feasibility and value.