Fluids Alter Elasticity Measurements: Porous Wave Propagation Accounts for Shear Wave Dispersion in Elastography

In shear wave elastography, rotational wave speeds are converted to elasticity measures using elastodynamic theory. The method has a wide range of applications and is the gold standard for non-invasive liver fibrosis detection. However, the observed shear wave dispersion of in vivo human liver shows a mismatch with purely elastic and visco-elastic wave propagation theory. In a laboratory phantom experiment we demonstrate that porosity and fluid viscosity need to be considered to properly convert shear wave speeds to elasticity in soft porous materials. We extend this conclusion to the clinical application of liver stiffness characterization by revisiting in vivo studies of liver elastography. To that end we compare Biot’s theory of poro-visco-elastic wave propagation to Voigt’s visco-elastic model. Our results suggest that accounting for dispersion due to fluid viscosity could improve shear wave imaging in the liver and other highly vascularized organs.


INTRODUCTION
Shear wave elastography is a fast, non-invasive method to estimate tissue elasticity in vivo. Its first clinical application has been the estimation of liver fibrosis, and it is more and more used as a replacement for liver puncture, the current gold standard. The method is quantitative through in situ measurement of the shear wave speed by ultrasonic or magnetic resonance imaging. An estimate of tissue elasticity is retrieved through inversion of the wave speed. Clinical applications as of today include ultrasonic 1D elasticity estimation (FibroScan, Echosens, Paris, France), ultrasonic 2D elasticity mapping (e.g., Aixplorer, SuperSonic Imagine, Aix-en-Provence, France; GE Logiq, Boston, United States; Canon Aplio, Canon Medical Systems, Tochigi, Japan) and magnetic resonance 2D and 3D elasticity mapping (e.g., General Electric/Philips/Siemens Healthcare). Linear elastic or viscoelastic theory is commonly assumed for the elasticity inversion. However, the order of magnitude of the shear modulus retrieved from elastography and other measurement methods, such as dynamic mechanical analysis (DMA), differs [1]. Furthermore, ex vivo and in vivo elastography measurements vary, making a standardization difficult [2,3]. Finally, the observed in vivo shear wave velocity dispersion with frequency of the liver is badly predicted by visco-elastic theory, and recent in vivo studies [4,5] challenge the commonly assumed correlation of viscosity with fibrosis [6]. As of today questions remain concerning the physical processes behind the observed dispersion. This ambiguity has led recent studies (e.g., [7]) to drop the rheological model in favor of a model-independent dispersion slope as the "viscous" parameter. Our work attempts to explain the observed discrepancies between theory and in vivo elastography data by considering the liver as a porous biphasic medium, which is described by Biot theory, instead of a homogeneous visco-elastic solid. The introduction of a fluid phase introduces a new loss mechanism, the viscosity of the fluid, as well as a change in the effective density of the solid part. Proper physical modeling of the of the shear wave dispersion is key to a quantitative elasticity reconstruction and ultimately might reduce false positives and negatives in elastography liver characterization. We postulate that to properly account for the observed shear wave dispersion in the liver, porous theory needs to be taken into account. This hypothesis is based on our experimental results in a porous phantom and theoretical considerations which are compared to existing in vivo data. The paper is structured as follows: In Section 1.1 we report various approaches to model the liver as a porous material. Publications in elastography which regard the liver as a poro-elastic Biot medium are scarce and consequently the investigated literature spans includes from anatomy, dynamic mechanical analysis, surgery simulation, tissue engineering and numerical fluid dynamics. In Section 3.1 we confirm the results of [8] who found that the shear wave dispersion in a soft porous material is described by poro-elastic wave propagation theory. This is unambiguously shown by comparing the shear wave dispersion in a soft porous phantom saturated by a Xanthan solution with the dispersion of the same phantom saturated by water. The observed results are exclusively explained by poro-visco-elastic theory, which is introduced as Biot's theory in Section 2.1.1. Furthermore, a more intuitive interpretation of Biot's theory is given in Section 2.1.2 and the experimental results are compared to Biot's model in Section 3.1.1. The experimental results are followed by analytic modeling of the liver as a porous material in Section 3.2. We use literature values for the poro-visco-elastic parameters (Section 2.6.1) to compare the theoretical dispersion of a model liver with previously reported in vivo shear wave dispersion measurements (Section 4.2). The analytic modeling shows that changes in the porous parameters can account for measured variations in shear wave speed that visco-elasticity fails to explain (Section 3.2.1, Section 4.2.1). We conclude by comparing Biot theory with commonly used rheological models in Section 4.3.1. Finally, we present challenges for future research in Section 4.3 and a common conclusion in Section 5.

The Liver-A Sponge?
Recently, we have shown experimentally that Biot's theory of poro-elastic wave propagation well predicts the shear wave dispersion in a melamine sponge [8]. In Biot's theory of poroelastic wave propagation shear wave dispersion is caused by the viscosity of the fluid saturating the connected pore space. Several organs that are of interest to elastography imaging have been described as porous materials, e.g., the lung, spleen, brain and liver. The liver is the foremost organ to be clinically imaged using shear wave elastography. Here, it is usually modeled as an elastic or visco-elastic solid. However, fresh liver grafts clearly show the presence of at least one liquid phase: the blood. In fact, the liver is one of the blood-richest organs and has previously been depicted as a fluid-filled sponge by several research groups [9][10][11]. Together with other organs, such as the spleen, the liver serves as a specific blood reservoir of the body [12]. It is thus justified to regard it as a soft porous material saturated by a viscous fluid. Furthermore, blood saturation or liver perfusion and blood viscosity are not constants. Both vary among individuals as well as with time for each individual. Perfusion and viscosity are known to be dependent on many factors such as hydration, diet, health and physical fitness. Liver diseases do not only alter the elasticity of the cellular matrix of the liver but are known to affect the porous parameters, notably porosity and permeability, as well. These insights led different groups to numerically model the liver as a poro-elastic or poro-visco-elastic material in the last decade [1,9,[13][14][15][16][17][18][19]. However, albeit early realization of interstitial fluid flow as an dissipative mechanism for acoustic waves [20], few attempts to link the shear wave dispersion in the liver to its porous and fluid properties have been made to our knowledge. Some recent studies investigated the influence of pore pressure on shear wave speed in the liver [21][22][23], and [24] describes a 2-parameter pore channel flow model for shear wave dispersion. The first group focuses on pressure effects. The second develops a model that incorporates flow in the constitutive equations, similar to Biot's first porous theory, the theory of consolidation [25], and to earlier work in compression elastography [26]. Recently [27], model the liver as a hyper-porovisco-elastic but focus on hyper-elastic deformation and the fluid pressure effect on the shear wave velocity by assuming low porosities and permeabilities. The effect of the viscous fluid on the dynamic wave propagation however, as developed by [28][29][30][31] and presented in Section 2.1.1, seems to have been overlooked by the ultrasound elastography community. Hence, the application of Biot theory to ultrasound shear wave elastography seems overdue.

Biot Theory
A non-phenomenological theory that models porous wave propagation based on the physics of fluid-solid interaction is Biot's theory of poro-elastic wave propagation [28,29]. Biot's shear wave dispersion is governed by the characteristics of the porous structure and the pore filling fluid. As a consequence, the expression for the wave speed differs significantly from that of a purely elastic monophasic wave propagation model. Instead of a single homogeneous density, the model takes into account the fluid density ρ fl and the density of the solid frame ρ fr , which is expressed in terms of bulk density of the solid material ρ s and the porosity ϕ as ρ fr ρ s (1 − ϕ). The mean density is expressed as ρ ρ s (1 − ϕ) + ρ fl ϕ. Furthermore, the model accounts for inertial and viscous coupling between the displacements of the two phases.
The consideration of viscous effects through the fluid viscosity η f leads to the frequency dependence of the shear wave speed, which is expressed in terms of wavenumber and angular frequency: where μ fr is the shear modulus of the foam. The term q can be expressed in terms of the solid, fluid and porous parameters as: where κ is permeability, α is tortuosity and F(ζ) expresses the frequency dependence of the viscosity, or the deviation from Poiseuille flow [32]. The permeability is measured in m 2 and is a measure of fluid transmission through a sample. It is defined as the proportionality constant in Darcy's law of the volumetric flow rate. The tortuosity takes into account the pore shape. It is generally assumed for cylindrical pores and is evaluated as the total length of the pore divided by the shortest path connecting its ends or in other words: "the ratio of the arc length of a curve relative to its end-to-end distance" [16]. It also accounts for the inertial coupling of the solid and fluid phase, where inertial coupling signifies that an accelerating solid drags parts of the surrounding fluid with it [33]. The solid moving in the fluid thus acquires an additional mass due to the deflection of the surrounding fluid volume, which can be seen as a renormalization process (see [34] 1.3.4). The virtual mass of the solid is thus ρ virt ρ fr + ρ a (1 − ϕ)ρ s + ρ a , where the additional or apparent mass is calculated as ρ a (α − 1)ϕρ fl . The derivation from first principles and a detailed description of Biot's theory are beyond the scope of this work and can be found in the cited references as well as several textbooks, e.g., [35][36][37]. In order to follow our reasoning we find it however helpful to give an intuitive interpretation that is accessible to the reader who was unfamiliar with poro-elastic theory before reading the present manuscript.

Loss Mechanism and Significance of Biot's Theory
When compared to linear-elastic or visco-elastic theory some results of Biot's theory might seem counter-intuitive. We can however understand the complex inertial and viscous effects of the fluid in terms of classic elastic theory by looking at the extreme ends of the dispersion curves when ω → 0 and ω → ∞. If we consider that the shear modulus μ fr stays constant with frequency and we want to find an equivalent to Biot's wave speed in terms of the elastic shear wave speed expression v s μ ρ , we immediately see that the density ρ needs to be replaced by an effective density ρ eff (ω → 0) and ρ eff (ω → ∞). At very low frequency, when ω → 0 Hz, the solid and the fluid phase are locked together through viscous coupling. This signifies that due to the viscous drag they move in phase without relative movement between the two. The wave speed is then defined by the true shear modulus of the solid and an effective density made up of the densities of the solid and the fluid phase (ρ eff ρ ρ s (1 − ϕ) + ρ fl p ϕ). In the case of biological tissues, those densities are approximately equal and the low-frequency limit is thus the velocity which we would have calculated using linear elasticity or visco-elasticity and the true shear modulus. At very high frequencies, when ω → ∞ Hz, which is conceptually similar to the case when the porous material is saturated by a superfluid, the shear wave velocity reaches its maximum value. Because of the lack of viscosity, the fluid is not dragged along by viscous coupling and rests on average 1 . The fluid and solid displacements are thus out-of-phase and the shear wave experiences a different effective density: the solid frame density plus an apparent density due to the remaining added mass effect of the fluid. The added mass effect stems from the fact that the fluid is considered to be at rest on average, while there is relative movement between the accelerating solid and the fluid phase (see [38] I, 11). Ignoring the inertial coupling described in 1 , the expression for the effective density is then: The higher the porosity of the medium, the lower the solid frame density and consequently due to v s μ ρ eff , the higher the shear wave speed in the high-frequency limit. Even for biological tissue, where fluid and solid densities are very similar, the shear wave speed at higher frequencies differs thus strongly from the elastic shear wave speed. While the idealized viscous-less case defines the high-frequency limit of the shear wave velocity and the phase-locked case defines the low frequency limit of the shear wave speed, it is obvious from the gap between the two values that a different physical mechanism is required to explain the evolution from low to high frequency. This is done through fluid viscosity, which introduces dissipation and defines the form of the dispersion curve in between the low and high frequency limits. Eqs 2, 3 account for the frequency dependence:ρ is scaled by the factor q which is a frequency dependent function of the fluid viscosity. To stay within the effective density model, the term q ρ could be interpreted as an effective density which is given through the frequency dependence of q(ω).

Visco-elastic Wave Propagation
In Section 3 we compare the poro-visco-elastic model to a simple visco-elastic model. The shear wave speed of the visco-elastic model is calculated using Voigt's model: with μ 0 being solid elasticity and μ 1 being solid viscosity. In the present work we explicitly make the choice to compare Biot's model only to a simple visco-elastic model. Viscosity is still not widely inferred for in clinical practice [39] and most in vivo studies employ a Voigt model with limited success [7]. As a consequence, a model independent measure of the dispersion slope has been proposed as diagnostic value [7] and more complex visco-elastic models, such as Zener's model are being investigated [7,40]. However, these include more springs and dampers and require complex inversions. One key difference between viscous and porous models is: the porous parameters in Biot's model are physical parameters that can be measured independently or determined from physical considerations [40], while the visco-elastic model is phenomenological. Hence, if we were to measure the porous and fluid parameters independently of shear wave elastography, the Voigt and Biot model would exhibit the same number of free parameters, the solid elasticity and viscosity. Note that [40]

Parameter Inversion
In order to fit the experimental data to the Biot and Voigt model we perform a constraint global minimization using the Basinhopping algorithm [41] of the Python library Scipy. The global method was chosen instead of a local method to avoid local minima in the optimization. In short, the Basinhoppin algorithm performs a multitude of local minimizations though the L-BFGS-B (Broyden-Fletcher-Goldfarb-Shanno) algorithm, which is a quasi-Newton method [42]. For each minimization it uses varying starting values and keeps or rejects the result. In Section 3 all minimizations are performed on the least squares [42] of the experimental and analytic shear wave speed of the investigated frequency range: where N is the number of frequencies at which the shear wave speed was measured, v exp is the experimental phase velocity, v sim is the analytic phase velocity and ω i is a experimentally measured frequency. Eq. 2 is used to calculate the poro-visco-elastic shear wave speed v sim . The shear wave speed v sim of the visco-elastic model is calculated using Voigt's model from Eq. 5. The bounds for the different minimizations performed in Section 3 are listed in Table 3.

Penetration Depth
The penetration depth is comparable to the skin depth in electromagnetics and relates the mathematical model outcome of any attenuation model to physical wave propagation (see [40] 2.3.2.2). For example, a model with an elasticity μ 0 → 0 and viscosity μ 1 ≫ μ 0 will result in a calculated wave speed v s > 0 but a penetration depth δ skin → 0. Hence, the wave is a nonpropagating one. We use the penetration depth as a physical constraint on the feasibility of a proposed model with: where α s is the shear wave attenuation, v s is the shear wave speed and ω is the frequency for which the physicality of wave propagation is investigated. We classify a model as being physical if δ skin > 0.5λ. The theoretical penetration depth separating a diffusion from a propagating wave is lower (see [40] 2.3.2), but we judge time of flight and phase velocity measurements on less than half a wavelength as unreasonable in shear wave elastography because near-field effects are likely to dominate the wave-field [43,44].

5.1 Experiment
Similar to [8] we saturate a rectangular Basotect ® melamine resin foam and analyze the shear wave dispersion through the phase velocity measurement of a frequency sweep. Robust results are achieved through ultrafast shear wave imaging of plane shear waves inside the foam. Figure 1 shows the experimental setup: The foam of dimensions 30 × 18 × 12 cm is immersed in a water tank. A metal rod driven by a piston (ModalShop Inc. K2004E01) excites shear waves in the frequency range from 100 to 250 Hz. The chosen frequency range ensures sufficient wave propagation and avoids guided waves which exist at lower frequencies due to the finite sample size. The resulting wavefield is imaged in the x − z plane using shear wave elastography. Ultrasonic plane waves are emitted orthogonal to the direction of shear wave propagation at 4,000 fps using a Verasonics Vantage ™ scanner and a Philips L7-4 probe. Phase-based motion estimation [45] of the ultrasound reflection images results in retrieval of the z-component of the particle velocity throughout the imaging plane. Because the shear wave speed in the sample is two orders of magnitude lower than the ultrasonic wave speed, the particle displacement induced by the shear wave propagation is well resolved in time. The shear wave excitation is exemplary shown in Figure 1A) and the imaging principle in Figure 1B). Because the wave excitation source is extended in z, the resulting wave front is planar, which is visible in the insets in Figure 1. Hence, we sum the particle velocities along the z-axis in order to increase Signal-to-Noise Ratio. The resulting 2D space-time wavefield is then Fourier transformed using Matlab's ™ FFT in order to measure the phase velocities of the frequency sweep. Because a plane wave evolves according to e i(kx−ωt) , the phase of a single frequency evolves linear in space. The phase velocity can thus be measured from the linear fit of the phase along the spatial dimension for each frequency. The phase is fitted using the Random Sample Consensus (RANSAC) algorithm [46,47] adapted from [48]. The benefit of the RANSAC algorithm lies in defining a set of inlier points which are taken into account for slope estimation and a set of outlier points which are disregarded. We employ a coefficient of determination (R 2 ) larger than 0.98 and 70% inliers for each frequency. With ϕ meas ∈ Inlier, the phase velocity is thus measured as: where ϕ is the phase, β is the linear slope of the phase, x is the coordinate of the spatial dimension, C is a Constant and ω is the Frontiers in Physics | www.frontiersin.org September 2021 | Volume 9 | Article 697990 frequency investigated. In contrast to our previous work reported in [8], we perform two experiments that differ only in the viscosity of the saturating fluid. For the high viscosity fluid we choose a 0.1% solution of Xanthan (E415) in water and measure the dispersion at 40°C, where 40°C ensures a good mixing of Xanthan with water. For the low viscosity fluid we measure the dispersion in water at 40°.

Poro-Visco-Elastic Dispersion Modeling in Melamine Foam
The Biot parameters of the foam were independently measured using the acoustic impedance tube method [49] and a Johnson-Champoux-Allard-Lafarge model [49][50][51][52] for [8]. For the different methods detailed in [49], the porosity ϕ was found between 96.7 and 99.7%, the tortuosity α between 1 and 1.02 and the permeability κ between 1.28 × 10 -9 and 2.85 × 10 -9 m 2 . The viscous length σ is between 11.24 × 10 -5 and 13.02 × 10 -5 m as indicated in the microscopic photo at the top of Figure 1, and the density is 8.8 kg m −3 ± 1 kg m −3 . The viscous length is typically employed in modeling air filled porous materials and not needed in the classical Biot model. The other quantities were fixed to ϕ 99%, κ 1.276 × 10 -9 m 2 and α 1.02. The elastic parameters are unknown and are likely to have been altered since the experiments performed by [8] on the same foam sample due to degassing, compression, immersion in water and aging. To investigate if Biot's theory can explain the observed dispersion we perform a joint least squares inversion [53] on the unknown parameters μ fr (shear modulus of the foam frame) and η f (fluid viscosity). To account for viscosity of the foam itself we introduce the complex shear modulus μ fr μ 0 + iωμ 1 , where the range for μ 1 is taken from the values given for melamine in [54]. This extended Biot theory, including solid relaxation mechanisms, has been described by Biot [31], and shown to be consistent with the Kramers-Kronig relations by Turgut [55]. We jointly invert using Eq. 6 and Python's Scipy package [56] as detailed in Section 2.3. Because we use the same foam sample at equal temperature in both experiments, μ fr is constrained to be identical for both experiments. In contrast, the fluid viscosity is bounded to 0.1-100 mPas for the Xanthan experiment and fixed to 0.6 mPas (viscosity of water at 40°) for the foam in the water experiment.

Porous Parameter Values of the Liver
The liver can be roughly separated into hepatocytes, blood vessels, bile ducts and space of disse or perisinusoidal space. The latter is part of the interstitial space and contains extracellular fluid. The blood vessels follow a hierarchical tree-like structure with the central vein and portal arteriole at the top and the liver sinusoids at the bottom. The meso-and micro-circulation are exemplary show in Figure 2, which is reproduced with permission from [57], and are gained from a liver corrosion cast. The macro-, meso-and microscopic structures exhibit different scales of porosity and permeability. While the portal vein for example can be seen as a tubular conduct with a very high Darcy permeability (1 × 10 -3 -1 × 10 -2 m 2 [58]), the permeability of the sinusoidal space has been modeled with permeabilities as low as 1.56 × 10 -14 m 2 and the interstitial space with 3.12 × 10 -17 m 2 [13,59].
Most studies reporting liver permeability values take a microstructure approach, investigating single liver lobules or the permeability across vessel walls. A more general approach was taken by [17], who modeled a 3D liver using a tissue permeability of 1.44 × 10 -9 -2.5 × 10 -9 m 2 and a vascular permeability of 5.0 × 10 -9 m 2 . Here, tissue permeability is the permeability derived from the sinusoids. A permeability of 1.4 × 10 -9 m 2 for an avascular model was also found by [60] for a poro-visco-elastic liver model. These relatively high values are not surprising, considering that the blood flow into the liver sinusoids equals about 1.35 Lmin −1 at a pressure drop of 1.2 kPa [12]. As is the case of the permeability, care has to be taken when comparing literature values for porosity. Considering only sinusoidal porosity, [61] find 7.4 ± 1.1% while [13,62] state values ranging from 14% for a healthy to 24.6% for a cirrhotic liver. According to [10], total liver porosity values range from 50 vol% to 70 mass% of total free fluid, including blood, water and bile. These porosity and permeability values are higher than those used recently by [27] in a poro-hyper-visco-elastic model. They base their consideration on the values by [13], who modeled microcirculatory or sinusoidal permeability based on hexagonal liver lobules. We believe however, that the sinusoidal permeability is not the relevant quantity for shear wave elastography. The region of interest for any liver elastography method is in the centimeter range and far exceeds the scale of a liver lobule (0.15 3 mm 3 in [13]). Therefore, one should take into account mesocirculation and potentially macrocirculation for proper modeling of the permeability and porosity in shear wave elastography.

Variability of Porous Parameter Values
For the average adult the normal blood volume in the liver equals about one third of its mass [63] or 450 ml. However, experiments in cats have shown that the liver can account for up to 400 ml of additional blood per kg of liver and rapidly expel half of the stored blood from the liver [63,64]. Concordantly, [12] states that the human liver can occasionally store up to 0.5-1 L of additional blood in the hepatic vein and sinusoids. Using ultrasonic measurements [65] found mean diurnal variations of 17% in human liver volume. Recently [66], found that mice liver volume oscillates with the feeding and circadian cycle, gaining and loosing 34% each day. Furthermore, fibrosis and cirrhosis have been shown to change the capillarization and thus the porous microstructure of the liver [67][68][69]. For harmonic ultrasound elastography, [70] found that water uptake leads to a change in the hepatic shear wave velocity. They conclude that the shear wave speed is most likely related to a change in vascular load after water intake. These blood volume variations on short timescales for a single organism as well as pathologic changes due to liver disease influence the shear wave speed as predicted by poro-visco-elastic theory. They notably alter the porosity and permeability parameters in Biot's model. Furthermore, the higher the porosity, the more Biot's model becomes sensitive to the fluid viscosity. Our research revealed few studies reporting the tortuosity of the liver pore space. For rat liver sinusoids it was found to be 1.04-1.11 by [71].
The above cited studies give the range of available literature values for the porous parameters of the liver and will serve as input to the analytic modeling of the liver as a poro-visco-elastic solid in Section 3.2. Table 1 summarizes the porous and elastic parameter values reported for the liver.

Permeability Calculation
The permeability of a porous sample with tubular pores of a single fixed radius r is given as [79]: FIGURE 2 | Visualization of the liver vasculature as presented by [57]. The blood vessels, present at all scales, constitute the main source of porosity in our poro-viscoelastic liver model. (A). Reconstruction of three mesovascular trees (hepatic artery, portal vein and hepatic venous system) from liver vascular corrosion casts. Reproduced with permission from Fig. 2 e) in [57]. (B). Scanning electron microscopic image of the sinusoidal microcirculation. Reproduced with permission from Fig. 1 c) in [57].
Frontiers in Physics | www.frontiersin.org September 2021 | Volume 9 | Article 697990 Following [80], a more realistic permeability can be calculated as: where a μ is the mean radius of capillary tubes and f(a) is the probability distribution density function of the tube radius. In Section 4.3.1, we use Eq. 11 to calculate the permeability of the liver corrosion cast of [57].

Poro-Visco-Elastic Liver Model
In Section 3.2, we analytically model how the shear wave dispersion in the liver is influenced by the porosity and fluid phase in the poro-visco-elastic wave propagation model. The liver is regarded as a homogeneous solid saturated by a Newtonian fluid according to Biot's theory ( [31,55]. We compare these analytic curves of a hypothetical liver with some exemplary literature values of in vivo elastography measurements. The measurements are displayed in Figure 4 and stem from Figure 1 of [81]. The data was read off manually by taking seven and nine points from the experimental curves for the F1 and F3 curve, respectively. Finally, we compare the results to a best fit of the experimental results with Voigt's model. The porovisco-elastic modeling is equivalent to the melamine foam modeling presented in Section 2.5.2. However, the foam parameters are now replaced by liver values from the literature. Regarding blood as a Newtonian fluid is of course a simplification, albeit one that is commonly made. A more accurate description of blood is a non-Newtonian visco-elastic shear thinning fluid or sol where the viscosity depends on shear stress, oscillation frequency and Hematocrit [82]. Shear rate is a function of flow and vessel diameter, which vary widely throughout the liver [83]. We thus explicitly make the choice of reducing the model complexity by assuming an average Newtonian blood viscosity under systolic conditions as for example measured by [75]. The Biot parameters are interrelated, hence the sensitivity of the shear wave speed on one parameter depends on the choice of the others. This leads to a variety of possible dispersion curves within the range of literature values given in Section 2.6.1. Exploiting the parameter space is beyond the scope of this manuscript. Instead, in Section 3.2 we present the dispersion curves for a low and a high porosity with changing viscosity. This shows the variety of dispersion curves which are covered by the porous model with fixed elastic and changing porous parameters.
To reduce complexity and because the Biot model is less sensitive towards tortuosity than to porosity ϕ, permeability κ, and fluid viscosity η f we fix the tortuosity to 1.05 throughout our calculations, which is in the range reported in Section 2.6.1.
In Section 3.2 we use the same modeling method to compare the poro-visco-elastic model to statistics of experimental shear wave dispersion data acquired by [81] and to Voigt's model. Equivalent to how [81] measure the experimental dispersion, the dispersion of the analytic model is evaluated using the slope of the linear fit between 73 and 250 Hz. With this range we account for the experimental conditions of the measurements of Figure 5 [81]: The stated problems at low frequencies and the fact that the upper frequency for the fit of each point in Figure 5 varied depending on the quality of the experimental data. The five parameter optimization of the Biot model exposes the ambiguity every multi-parameter inversion suffers. We do not claim to find the true values explaining the experimental data here, but merely want to evaluate if physical changes in porous parameters as reported in Section 2.6.2 are able to explain the experimental trend observed by [81]. We thus explicitly compare with a two parameter, and thus less ambiguous, Voigt model, only. To facilitate the multi-parameter inversion for Biot's parameters we put further physical a priori constraints. The aim of the multi-parameter inversion here is to see if the observed statistical trend in the dispersion data can be explained by changes in the Biot model. We thus restrict the solid viscosity of the Biot model to not exceed that of the point with the lowest Dispersion, effectively fixing the Voigt viscosity for the entire observed data. Furthermore, we suppose increasing porosity and permeability for increasing dispersion, which is reflected in the optimization bounds. The latter results from Eq. 11: a higher porosity must be accompanied by a higher permeability if no geometry changes happen.

Experimental Shear Wave Dispersion in Melamine Foam
The red and blue crosses in Figure 3 show the shear wave dispersion curves for high and low viscosity fluids. The dispersion for the foam in water (blue) differs significantly from the dispersion of the same foam in Xanthan solution (red). For example, the measured phase velocity in the foam  Table 2. The Voigt model is displayed to highlight the differences with the poro-visco-elastic model. It is not a fit but uses the same shear modulus, solid viscosity and average densityρ ρ s (1 − ϕ) + ρ fl ϕ as the Biot model.
Frontiers in Physics | www.frontiersin.org September 2021 | Volume 9 | Article 697990 saturated by the Xanthan solution at 125 Hz is 2 m s −1 (>20%) lower than that of the foam saturated by water.

Poro-Visco-Elastic Dispersion Modeling in Melamine Foam
An inversion for the free parameters μ 0 and μ 1 of the Biot model, as described in  Figure 4 shows the modeling results for a liver with different values of porosity and permeability and an identical set of visco-elastic parameters. The exact parameter values can be found in Table 2.

Porous Dispersion in the Liver
The blue dashed curve represents a low porosity-low fluid viscosity liver, while the blue straight curve represents a high porosity-low fluid viscosity curve. The superimposed black curves are two exemplary shear wave dispersion curves in the liver acquired through ARFI/SWEI with an Siemens SONOLINE Antares scanner and a CH4-1 transducer (Siemens Healthcare, Ultrasound Business Unit, Mountain View, CA), taken from [81]. The dashed curve represents a liver classified with Fibrosis stage F1 and the straight curve represents a liver classified with Fibrosis stage F3. The least-squares best fit of the visco-elastic Voigt model and the experimental data in the frequency range of 100-400 Hz are plotted in red. Both, the poro-visco-elastic model results and the visco-elastic model results respect a penetration depth superior to 0.5 λ (see Section 2.4). The dropoff of the experimentally measured wave speeds at low frequencies is attributed to a roll-off of the Fourier transform by Ref. [81]. They express confidence in their measurements for values >73 Hz. Within this confidence interval the measurements show an inflection at 100 and 150 Hz, which the Voigt model is incapable of reproducing without boosting μ 1 to large values that lead to a penetration depth <0.5λ. In contrast to that, the observed inflection is an intrinsic property of the Biot model with steeper slope at lower frequencies for smaller fluid viscosity. The dotted blue curve represents the modeling result of using the same set of parameters as for the straight blue line except for the fluid viscosity η f . Fluid viscosity is tripled for the dotted blue curve, leading to a significant decrease of the shear wave speed. The dash-dotted blue line was achieved likewise from the model of the dashed blue line. At 200 Hz, a commonly taken shear wave speed reference frequency, the difference between the low porosity-high fluid viscosity (dash-dotted) and the high porosity-low fluid viscosity (straight) model reaches 0.8 ms −1 or 42% of the lower value.

Experimental Support for the Poro-Visco-Elastic Liver Model
Reference [81] conducted a statistical analysis of in vivo measurements of shear wave dispersion slope at 200 Hz versus shear wave speed at 200 Hz for healthy and fibrotic human livers. In Figure 6 of [81], which is replotted in our Figure 5, one can observe a positive correlation (R 2 0.3) of the steepness of the dispersion slope with the measured speed at 200 Hz [81]. Eyeballing the trend we set three dispersion-speed target pairs for the inversion: 2.0 ms −1 -0.25 ms −1 100Hz , 2.5 ms −1 -0.45 ms −1 100Hz , 3.0 ms −1 -0.65 ms −1 100Hz . Using Eq. 6 and the Basinhopping algorithm [41] (see Section 2.3) the inversion for Voigt's model results in the points V1,V2, and V3 of Figure 5. They represent the least squares best fits of the Voigt model and the speed-dispersion pairs. Using the same speed -dispersion pairs, the least squares fit using the poro-visco-elastic model results in the points B1, B2 and B3. The exact parameter values for V1,V2, V3 and B1,B2, B3 can be found in Table 2 and the bounds for the inversion in Table 3.

Poro-Elastic Loss Mechanism of Melamine Foam
The experiment confirms the results by [8] that the shear wave dispersion in a highly porous soft phantom is governed by the  Figure 1 in [81]. "-" F3 Fibrosis classification, "--" F1 Fibrosis classification; Blue: Poro-visco-elastic model developed in this paper, "-" high porosity and permeability, "--", low porosity and permeability,"··" same as "-" with tripled fluid viscosity, "·-" same as "--" with tripled fluid viscosity; Red: Voigt's visco-elastic model, "-" high shear modulus, "--" low shear modulus. The parameter values for the analytic curves are reported in Table 2. viscosity of the fluid. Since the dispersion of the foam in Xanthan solution and in water differ significantly at the same temperature, viscous effects of the solid matrix cannot be the single reason for the wave dispersion. The experiment clearly shows that the difference in dispersion is due to the changed viscosity of the pore filling fluid. As a consequence, it is impossible for Voigt's model to reproduce the observed difference in the dispersion curves. Thus, modeling of shear wave propagation in soft porous materials by poro-elastic theory is justified. Biot's nonphenomenological theory of poro-elastic wave propagation takes the fluid viscosity into account by introducing a mechanism of fluid-solid interaction [28,29].
It should be noted that the rheology of Xanthan solutions is non-Newtonian and as such Xanthan shares the shear dependency of the viscosity with blood. However, Biot's approximation of a Newtonian fluid with one constant apparent viscosity is sufficient for the small frequency range investigated herein, and in order to show that fluid viscosity is the single parameter which alters the shear wave dispersion. The good fit of the dispersion curves, which was retrieved by restricting the elastic and viscous parameters to lie within the range of the literature values for melamine foam, indicates that Biot's theory is well suited to explain the observed shear wave dispersion.
Our experiments show a fundamental difference between porous theory and other rheological models such as Voigt's model: the effect of a change in viscosity. A higher value of the fluid viscosity (Biot) leads to a smaller slope of the dispersion curve while a higher solid viscosity (Voigt) leads to a steeper slope of the dispersion curve. This can be understood in the following way: A higher fluid viscosity signifies stronger solid-fluid coupling. As a consequence, the shear wave encounters an increased effective density.
The experimental results were gained on a melamine sponge, which is an artificial, highly porous, soft solid. Melamine sponge has already been used to mimic ultrasonic scattering in biological porous tissue, namely the lung [84]. While the results of [8] and the experimental results presented herein show that Biot theory can explain shear wave dispersion in soft porous materials, its applicability on in vivo organs will be discussed in the following.

Porous Liver Modeling
The analytic shear wave dispersion curves in Figure 4 are modeled without a change of the elastic parameters and cover the experimental curves for a F1 and F3 fibrosis stage liver. While the Voigt model requires a near two-fold increase of the elasticity and solid viscosity, the poro-visco-elastic model can reproduce both curves without a change in neither the elasticity nor the solid viscosity. Hence, the measurements above 73 Hz could in 2 | Modeling parameter values for the poro-elastic and visco-elastic models. The fluid density and the mean density are set to 1,050 kgm −3 , which is the reported value of arterial blood [100]. ϕ-porosity, κ-permeability, η f -fluid viscosity, μ 0 -shear modulus, μ 1 -solid viscosity.  principle be explained without a change in solid visco-elasticity by using poro-visco-elastic theory. Since fibrosis classification by elastography relies on changes in shear modulus, this would indicate the possibility of false positives in fibrosis classification due to varying liver saturation. If the poro-viscoelastic model can be confirmed in vivo, shear wave elasticity imaging should thus take the porous parameters into account. Considering that perfusion and blood viscosity have been shown to be highly variable on short timescales (see Section 2.6.2), we think that these modeling results should be considered in future clinical research.
The fluid viscosity in the poro-visco-elastic model is 3.1 mPas, which is below the commonly assumed whole blood viscosity of 3.5 mPas but above the viscosities given for microcirculation and interstitial fluid in the liver by [59]. In a recent study [75], observed an anti-correlation of blood viscosity with liver stiffness as measured by transient elastography. The authors suggested pathologic changes as the source of this anticorrelation. The viscosity induced changes in shear wave speed of our porous modeling indicate the possibility of a different explanation. The correlation could be rooted in the inherent effect of the fluid viscosity on the shear wave speed, and as such be an artifact of the elasticity inversion.

Experimental Support for the Poro-Visco-Elastic Liver Model
Our hypothesis that the poro-visco-elastic model is more appropriate than Voigt's model to explain shear wave dispersion in the liver is supported by the analysis of the experimental results by [81] reported in Figure 5. In Voigt's model, an increase in the elasticity modulus will lead to an increase in measured speed at 200 Hz, but to a decrease of the dispersion slope. A much larger increase in solid viscosity (μ 1 ) than in elasticity (μ 0 ) could reproduce the correlation, but to our knowledge no such physical mechanism is known. On the contrary, the predominance of μ 0 for liver elastic behavior is well documented [78]. Furthermore, such an increase in solid viscosity leads to very low penetration depths, effectively making the modeled wave a non-propagating one. Unlike Voigt's model, Biot's poro-visco-elastic model is able to reproduce a correlation of increased shear wave speed with increased dispersion slope for the reported in vivo measurements.
Ex vivo measurements do not reproduce this correlation of dispersion slope with shear wave speed [85]. This could be owed to the fact that ex vivo livers are partially drained. The drainage changes the porous properties, which adds to other factors, such as fluid pressure, that are inhibiting ex vivo experiments from properly reproducing the physics that govern shear wave dispersion in vivo.
While assuming the purely visco-elastic model, the porous and liquid properties of the liver might thus have largely been overlooked previously. Studies supporting this statement were conducted by [21][22][23]86], who compared confined and unconfined ex vivo liver samples and found that the pressure on the fluid in the liver influences the shear wave speed. Reference [24] employed a microchannel flow model equivalent to Darcy's law and noted that it is superior to Kelvin-Voigt models, because it is based on tissue architecture instead of a priori information on springs and dampers. Both findings strongly support our deduction that porosity plays a crucial role in shear wave propagation in the liver. However, they do not discuss the main loss mechanism shown in Section 3.1: the viscosity of the fluid in the Biot model. While these works deduce elasticity changes from increased pressure on the tissue, we introduce the direct role of the fluid phase through porous physics into the calculation of the shear wave speed.
First, it is an entirely non-phenomenological model which is derived from first principles and the porous parameters can in principle all be measured independently. This is the main reason why it is widely applied in geophysical applications, e.g., in oil and gas exploration. Second, both the liver porosity and permeability depend heavily on blood vessel size and are not independent variables. The porosity, i.e., the space occupied by free fluids, can be estimated as the total volume of blood and extracellular fluid. Tubes are a good approximation for the average blood vessel shape and as a consequence permeability and porosity can be calculated from the knowledge of the distribution of vessel sizes with Eq. 11.
Reference [69] analyzed the 3D micro and mesoscopic structure of the hepatic blood vessel tree using high resolution Computer Tomography of a liver cast. As a final check for the Biot model we create a normally distributed probability function for each vessel generation based on their values of the mean and standard deviation of the pore size and length per vessel generation. For the interstitial fluid we assume a pore radius of 0.25 µm, a value derived by [59] for the space of disse, and a Volume percentage of 12% [77,88]. Integrating over the permeabilities for each vessel generation we retrieve κ ≈ 2e −9 m 2 and ϕ ≈ 0. 35. Theses values are in the lower range of the literature values [17], which is in accordance with the assumption that the cast represents an un-or normally perfused liver. Reference [57] also evoke a possible shrinking of the blood vessel diameters due to the casting process.
We assume that a change in liver perfusion in vivo is reached through vaso-dilatation or vaso-contraction [63,64,89]. We can then calculate the necessary dilatational factor to reach the maximum porosity of 70% for the poro-visco-elastic model: For example, increasing the radius of the smallest 8 vessel generations by two and of the remaining generations by 1.2 as well as their lengths by 1.5 for the smallest 8 and by 1.2 for the remaining generations gives a porosity of 70% and a permeability of 5.6 × 10 -9 m 2 . This rough estimation is in accordance with the upper limits of the literature 2 . Furthermore, the dilatational factors are not purely speculative, since [68] found a dilatational factor for the sinusoids of nearly two in cirrhotic rat livers.
As of now, these theoretical porous parameter values are estimated from the liver cast and the literature, but recent imaging advances in ultrasound [90][91][92] and MRI imaging [93] have shown that porosity and blood vessel size estimation in vivo are becoming feasible. Another approach that might advance the imaging of porous properties with shear wave elastography was recently presented by [94]. They suggest to use Bayesian inference to deduce microstructure properties from effective wave parameters. Hence, poro-visco-elastic shear wave speed estimation could in the future replace the phenomenological springs and dampers of the currently applied rheological models with a physical loss mechanism that is determined by measurable porous parameters. From a physicist point of view, the added model complexity is then well justified by a reduction in the model ambiguity. If the poro-viscoelastic model will also result in improved disease classification compared to more complex visco-elastic models will have to be evaluated in the future.

Low-Frequency Dispersion
References [95,96] show, that in the case of magnetic resonance shear wave elastography of the brain, a poroelastic finite element model following [97][98][99] is more accurate than a visco-elastic model at very low frequencies.
Physical reasoning for the observed low speeds at low frequencies requires much lower elastic moduli than are commonly inverted for in transient elastography. An elastography independent measurement of such low elastic moduli is given by [9,78] for in vitro porcine liver. They find storage moduli as low as 300 Pa, using dynamic mechanical analysis (DMA). Similar values are reported for ex vivo bovine livers by [3] on fresh porcine liver grafts, using a rheometer and a fractional Kelvin-Voigt model, and by [60], who employed a poro-hyper-elastic model. Even lower elastic moduli are reported for in-vivo mimicking perfusion tests on ex-vivo liver samples. Reference [10] report linear elastic shear modulus values of less than 100 Pa for liver parenchyma.
These results in combination with poro-visco-elastic modeling indicate that the entire observable dispersion range, including the low frequencies, needs to be taken into account in order to accurately retrieve the elasticity of the liver. Future in vivo studies of the liver should thus focus on acquiring shear wave speeds over a large bandwidth. At very low frequencies, frequency dependent fluid viscosity and permeability should be taken into account due to the pronounced visco-elastic properties of blood.

CONCLUSION
Our results show, that for fluid saturated, soft porous materials, the porous parameters need to be taken into account to properly model shear wave dispersion. The shear moduli retrieved from simple visco-elastic models are prone to overestimations in the presence of porosity. Fluid viscosity, permeability, and porosity of the diseased as well as the healthy liver can vary inter-and intraindividual on short timescales. It is thus justified to assume that the wide range of observed shear wave dispersion curves in the in vivo liver is partly owed to changes in the porous properties of the liver. Our comparison of poro-visco-elastic and simple viscoelastic modeling shows, that porous parameter changes can easily be the reason for false positives in visco-elastic shear modulus inversion. In the here presented porous liver model porosity dominates at low frequencies, and solid viscosity is dominant at high frequencies. Furthermore, the stiffer the matrix of the liver, the more sensitive the shear wave dispersion becomes to the porous parameters in the Biot model.
Compared to more complex visco-elastic models, the porovisco-elastic model does not require further phenomenological parameters. Instead it relies on directly measurable quantities. Furthermore, it gives a physical explanation for the observed shear wave dispersion in the liver. If our results are confirmed in vivo, they are likely to have implications for elasticity imaging of other porous organs, such as the spleen, the kidney and the brain.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
JA and SC contributed equally to conception and design of the study. JA conducted the experiments, literature search and analytical modeling and SC participated therein. JA wrote the first draft of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

FUNDING
This work was supported in part by the Plan Cancer ARC-2018-PhysiCancer-BPALP.