Spring–damper equivalents of the fractional, poroelastic, and poroviscoelastic models for elastography

In MR elastography, it is common to use an elastic model for the tissue's response in order to interpret the results properly. More complex models, such as viscoelastic, fractional viscoelastic, poroelastic, or poroviscoelastic ones, are also used. These models appear at first sight to be very different, but here it is shown that they may all be expressed in terms of elementary viscoelastic models. For a medium expressed with fractional models, many elementary spring–damper combinations are added, each of them weighted according to a long‐tailed distribution of time constants or relaxation frequencies. This may open up a more physical interpretation of fractional models. The shear‐wave component of the poroelastic model is shown to be modeled exactly by a three‐component Zener model. The extended poroviscoelastic model is found to be equivalent to what is called a non‐standard four‐parameter model. Accordingly, the large number of parameters in the porous models can be reduced to the same number as in their viscoelastic equivalents. While the individual displacements from the solid and fluid parts cannot be measured individually, the main use of the poro(visco)elastic models is therefore as a physics‐based method for determining parameters in a viscoelastic model.


INTRODUCTION
In modeling of data from MR elastography, it is common to use a simple elastic model for the medium. This is the case even for ultrasound elastography. For more accurate modeling, three families of models are used: these are the linear viscoelastic models, such as the Kelvin-Voigt and Zener models, the fractional extensions of these models, and poroelastic models based on the theory of Biot.
The linear viscoelastic models are among those that have been used for fitting the frequency dependence of shear-wave data from MR elastography in the brain. 1 The fractional Kelvin-Voigt model was fitted to breast MR elastography data by Sinkus et al 2

and also analyzed by Holm and
Sinkus 3 and compared with other models for elastography by Zhang and Holm. 4 It may be argued that these single-phase models are too simplistic and that in tissue a biphasic model, which distinguishes between solid and liquid phases, would be more accurate. This potential has already been demonstrated in models of the quasi-static biomechanics of hydrocephalus 5 and infusion-induced swelling in the brain. 6 The poroelastic model has also been used for elastography. In the work of Konofagou et al, 7 ultrasound elastography was modeled with either an elastic model (i.e. without viscosity) or a simplified poroelastic model, which depended on porosity, composite density, and fluid density. Similarly, Perriñez et al 8 evaluated the full poroelastic model versus an elastic one for MR elastography. McGarry et al 9 went one step further and compared the poroleastic model with a viscoelastic one with a complex shear modulus. They also tried to use these models for inversion and observed that the viscoelastic model produced better reconstructions at 50 Hz, while the poroelastic model was superior at 1 Hz. HOLM pores but only sees a homogeneous effective medium. Further, since it is sensitive to signals from hydrogen atoms in the voxel, one cannot separate the signal from the solid and the fluid, even if the solid should contain fluid that is considered not to be in the pores. Likewise, ultrasound is scattered from the soft tissue, which is mainly part of the solid matrix. Ultrasound will therefore produce a strain estimate in the solid matrix that is only indirectly influenced by the wave in the fluid. 7 In this article, however, the goal is not to develop a poroelastic model for reconstruction of the MR elastography image but only to discuss it in the framework of explaining variations in parameters due to physiological and physical parameters.
The claim of this article is that these models are much more similar than they appear to be. The fractional models can be developed as sums of ordinary viscoelastic elements weighted in a particular way. The fractional model has its strength, in that it gives a parsimonious description of the phenomenon, i.e. one with a minimal number of parameters, especially when power-law variation in frequency is observed, but it is not fundamentally different. Likewise, the poroelastic model for the wave mode of interest for elastography, the shear wave, can be described in terms of standard viscoelastic models. In this case, it is the viscoelastic model that requires the smallest number of parameters. The strength of the poroelastic formulation is that it provides a way to find out how these parameters depend on physical and even physiological parameters. The linear viscoelastic model is therefore first generalized in Section 2 from simple two-and three-component models to chains of spring-damper elements described by time-and frequency-spectral functions. These functions are used to show that the fractional viscoelastic models in Section 3 can be described as sums of ordinary viscoelastic models. A surprising result is that the weighting in the sum follows a long-tailed distribution, reminiscent of that produced by a fractal geometry.
The shear-wave solution of the poroelastic theory is then developed in Section 4, from both the original formulation of Biot 11 and also that of Stoll. 12 Somewhat unexpectedly, it is found that there is a one-to-one correspondence between the poroelastic shear-wave response and that of a simple spring-dashpot network.

LINEAR VISCOELASTIC MODELS
The linear viscoelastic model is expressed in several different ways. These ways are needed in order to generalize from simple two-and three-component models (e.g. Kelvin-Voigt and Zener) to higher order models.
In order to illustrate the similarities between models, it is sufficient here to express the models in one dimension, although three dimensions are really needed for a complete shear-wave description. This section, as well as Section 3, builds to a large degree on the work of Mainardi. 13

Three descriptions
In linear viscoelasticity, there are three different ways of describing the medium's response. The first is the hereditary model of Boltzmann,14,15 where the constitutive equation is a convolution integral 13 : The kernel G(t) is called the relaxation modulus and represents a fading memory, i.e. one where changes in the past have less effect now than more recent changes. In order to ensure causality, the kernel, G(t), has to be zero for a non-negative time. The relaxation modulus is the strain response of the system to a step excitation in strain.
The second description is a linear differential equation between stress and strain with constant coefficients: [ The third description is in the form of a relaxation spectrum or Prony expansion (Mainardi's Equation 2.28) 13 : A physically realizable model is obtained if G(t) is modeled by a parallel network of N − 1 pairs of springs and dashpots in series where the spring constants are E n and the viscosities of the dashpots are n = n E n . Furthermore, both the spring constants and the viscosities are non-negative. This is the Maxwell-Wiechert model of Figure 1. The additional left-hand spring leads to the constant G e = E e , the equilibrium modulus, and if there is a dashpot directly across the terminals (e.g. if E 1 = 0) then the response will have an impulse at time zero as well, given by G − . This is the case in the Kelvin-Voigt model, which will soon be discussed.
The model of Equation 2 can also be expressed as a convolution, as in Equation 1. However, here we are interested in that subset of convolution kernels that can be characterized as fading memory kernels and not all the coefficient sets of Equation 2 will lead to such a model. Therefore the linear differential equation is the most general model. Also, a model described as a relaxation spectrum (Equation (3)) can always be described with a linear differential equation, but not vice versa. In addition, a relaxation spectrum model, Equation 3, always results in a fading memory model, as it consists of sums of positively weighted falling exponentials.

of 12
The frequency-domain equivalent of Equation 2 is also useful: The dynamic modulus, E( ), is often called G( ) in elastography, but, since G means the relaxation modulus in this article, E is used instead.

Elementary linear viscoelastic models
The two simplest viscoelastic models are the Kelvin-Voigt and the Zener models, as shown in Figure 2. The Kelvin-Voigt model is found by just keeping E e and 1 in Figure 1. The linear differential equation is and the relaxation modulus is Thus p = 0, q = 1, a 1 = 0, and b 1 = in Equations 2 and 4. Further, G e = E e , G − = , and G (t) = 0 in Equation 3. The dynamic modulus is where = ∕E e . The even simpler elastic model, which is often the reference in elastography, as noted in the Introduction, is obtained by setting the viscosity, , to 0.
The Zener model adds one more term to the linear differential equation of the Kelvin-Voigt model: and the relaxation modulus is For this model, the dynamic modulus is The addition of the extra spring leads to a frequency-dependent denominator in the dynamic modulus, which gives more degrees of freedom in fitting this model to data.
The dynamic modulus can also be included in a dispersion relation when propagating waves are involved. In Equations 13-15 of Holm and Näsholm, 16 it is shown that the complex wave number, k, in that case is where is the density and ( ) is the dynamic compliance or compressibility, the inverse of the dynamic modulus. This result is important in the subsequent analysis of the poroelastic model. Note that the roles of and are reversed here, compared with the work of Holm and Näsholm. 16 Here the convention of Mainardi 13 is followed instead.
A slightly rewritten version of the compliance based on factoring of Equation 11 will be needed when analyzing the poroelastic model: Judging from Figure 1, a Zener model is characterized by two springs and a damper. These three parameters may be expressed in different ways and one common way is via the low-frequency asymptote of the propagation speed, c 0 , and the two crossover frequencies, = 1∕ and = 1∕ .
The two frequencies express approximately where the phase velocity starts rising and where it approaches its asymptotic value, 16 Taking into consideration that c 2 0 = E e ∕ and thus depends on both the shear modulus and the density, then a Zener medium will depend on four independent parameters: E e , , , and .

Time-and frequency-spectral functions
The continuous generalization of Equation 3 is given by 13,17,18 where b = G (0) is a non-negative constant, which for the Zener model is There is a corresponding frequency-spectral function, which is By substituting Ω = 1∕ so that dΩ = −d ∕ 2 , this gives The relaxation modulus is, as noted, the stress response of the system to a step excitation in strain. The relaxation spectrum could also have been expressed with the creep response, which is the strain response to a step excitation in strain. It is denoted by J(t) and plays the same role as G(t) in where a = J (∞) is a non-negative constant, which for the Zener model is a = (1∕E)(1 − ∕ ). 18 This corresponds to a decomposition into a sum of elementary models in series, as shown in Figure 3. This is the conjugate of the Maxwell-Wiechert model of Figure 1, i.e. a different configuration of springs and dashpots that has the same characteristics. 19 It follows that the conjugate of the Zener model in the right-hand part of Figure 2 is also the top three elements of Figure 3.
The frequency-spectral function for the creep is found by a transformation analogous to that for relaxation: and the frequency-spectral function is This frequency spectral function will be used to find the spring-damper equivalent of the fractional models.   (2)). For the Zener model, this is The dynamic modulus is now while, in the simpler Kelvin-Voigt model shown in Figure 4, the constant = 0: The dynamic modulus is In the limit as the frequency and/or viscosity becomes very large, the dynamic modulus approaches E( ) → (i ) . This is equivalent to the case where the spring of Figure 4 can be neglected. That seems often to be the case for elastography data. 20 compliance of both models also follows a Mittag-Leffler function. That means that the creep time-and frequency-spectral functions of the two models will also be similar.

Long-tailed time-and frequency-spectral function
In previous works, 13,17,18 it has been shown that, for the fractional Zener model, the creep time-spectral function of Equation 17 is It was plotted by Caputo and Mainardi 13,18 with linear axes and is a decreasing function of for small , with a more and more pronounced peak as approaches 1. For the non-fractional case, = 1, it is a delta function at ∕ = 1, showing that it is equivalent to a single relaxation process in that case.
Here the function is plotted on a log-log scale in Figure 5. That brings out the properties of this particular function in a different way than if it were plotted in a linear plot. The log-log plot fixes attention on the asymptotes, rather than the peak, and they are The distribution of elementary spring-damper models in the medium, as given in Figures 1 and 3, has different asymptotes for small and large time and both of them are power laws. As long as the exponent of the high-frequency tail falls off more slowly than −2 , i.e. for < 1, the variance of the distribution will not exist. Such long-tailed distributions are scale-invariant.
In Mainardi's works, 13,23 the frequency-spectral function is also given. It can be found from Equation 18 using the normalizing constant a = J (∞) = (1∕E)(1 − ( ∕ ) ) for the fractional Zener model. The creep frequency-spectral function of Equation 19 is then This result was given by Näsholm and Holm 24 in the form above. There, it was derived with a starting point in the multiple relaxation theory of Nachman et al. 25 It should be noted that the frequency-spectral function is exactly the same as the time-spectral function, except for a scaling factor.
The plot for = 0.8 and = 0.01 in Figure 6 shares many properties with the plot of the time-spectral function. The asymptotes are also very similar 26 : Interestingly, the relaxation spectrum approaches a single power law for the limiting case → 0, where both the low-and high-frequency parts of the relaxation spectrum will approach Ω −1 . See the dash-dotted line for = = 0.01 in Figure 6, with asymptotes Ω −1+ for low frequencies

Relationship with fractals
Much effort has been expended on finding the relationship between fractional models and fractal geometries, beyond the semantic similarity. A general connection has not been found, but a geometrical and physical interpretation of fractional integration and differentiation was given by Podlubny. 28 There are also studies that are more specific to wave propagation and attenuation. The results are then often for the alternative attenuation mechanism of scattering rather than absorption. This has been shown theoretically by Fouque et al 29

Relationship with hysteresis
The results found here have some resemblance to those derived from hysteresis. Hysteresis is a hypothetical loss element with a constant phase lag between stress and strain at all frequencies, 35 resulting in attenuation that increases linearly with frequency. 36 The dynamic modulus for hysteresis is E( ) = K + iH (28) where K and H are constants. The model leads to a noncausal model and therefore 35 is concerned with making a band-limited approximation, where the constants are allowed to vary with frequency, K( ) and H( ), in order to ensure causality.
By comparing the dynamic moduli, Equations 28 and 23, it is evident that hysteresis can be viewed as the limiting case as → 0 for the fractional Kelvin-Voigt model. Thus an alternative way to ensure causality is to employ a fractional Kelvin-Voigt model with parameter = = 0 + .
Hysteresis is particularly interesting in the context of fractional operators, because the result for the limiting case → 0 of the previous section is similar to that recently found for hysteresis by Parker. 37 In Parker's work, it was shown that the hysteresis model is the same as a sum of relaxation processes weighted with a long-tailed power law Ω −1+ . This parallels our discussion of the asymptotic result of Equation 27 for = and Parker's result can be interpreted as a special case of Equation 27.

POROELASTIC MODEL
The Biot poroelastic model deals with a saturated porous medium with a solid phase and a fluid phase. Wave propagation in such a medium is described by a set of coupled vector wave equations, as given in Equation 4.2 of Biot. 11 The variables are the displacement vectors u and U for the displacement in the solid and the fluid, respectively. That article assumes that, as the fluid moves in the pores, the flow is laminar and that losses are given by Darcy's law and are proportional to the relative velocity, (u − U)∕ t. The theory predicts three solutions: two compressional waves and one shear wave. In a biological porous medium like cancellous bone (bone with a low volume fraction of solid, less than 70%), all three waves have been detected. 38 This model is much more complex than the viscoelastic models of the previous sections. It also seems more appropriate for a complex medium like brain, liver, or bone. Here, a surprising exact relationship between the poroelastic and viscoelastic models will be shown for the shear-wave solution.

Biot's original formulation
As elastography is only concerned with shear waves, we restrict the analysis here to that mode. The following dispersion relation was derived by Biot 11 (Equations 7.5-7.6 therein): where r = 1 1 + 22 This expression depends on three normalized densities: 11 = 11 , 22 = 22 , 12 = 12 (31) and the aggregate or composite density, which is These formulae depend on the porosity, , and the fluid and solid densities. The parameter 12 represents a negative mass coupling density between fluid and solid, 11 represents the total effective mass of the solid moving in the fluid, and 22 is the total effective mass of the fluid. There is also a characteristic frequency: where is the fluid viscosity and B is the permeability. The effective characteristic frequency of Equation 30 is, however, a lower frequency: 39 where is a structure constant related to tortuosity. Comparison with Equations 12 and 13 shows that Equations 29 and 30 are exactly the same as those of the Zener model. This is a remarkable result, which shows that, for shear waves, the poroelastic model is also a linear viscoelastic model. What distinguishes it from ordinary linear viscoelastic models is that it provides a sophisticated way of determining the parameters from physical considerations.

Biot-Stoll formulation
The parameters of the original Biot theory are often considered to be hard to find in practice and the theory has been rewritten in terms of the relative displacement between fluid and solid, u−U, or the volume of fluid that has flowed in or out of an element, = ∇(u − U), in combination with the solid displacement u. The material parameters are then transformed to a new set of parameters.
In that case, Hovem's 40 Equation 16.77 gives the shear dispersion of the low-frequency Biot model. It can also be found from Equation 12 of where the sign of has been reversed compared with the original articles, due to a different definition of the Fourier transform here (see also Appendix A of Holm and Näsholm 41 ). The mass coupling density is where is the tortuosity. Comparing Equation 35 with Equation 12 shows that the dynamic modulus is As expected, this is also equivalent to that of a Zener model, as can be seen by comparison with Equation 11. The constants when = 1∕ and = 1∕ are Insertion of the expression for the mass coupling density, Equation 36, in the equation for the effective characteristic frequency, Equation 34, shows that this frequency is also the same as ∕(2 ).

Redundant parameters in the Biot shear-wave model
The low-frequency Biot model is characterized by the ten parameters shown in Table 1. Seven of them affect shear-wave propagation, as indicated in the right-hand column, but several of them are connected, as Equation 38 indicates.
The parameter depends on the mass coupling density c . In addition it depends on the ratio of the viscosity, , and the permeability, B. They do not influence any of the other parameters, and c 0 , so it is clear that it is their ratio that matters. Likewise, the sound velocity, c 0 , depends on the ratio of the shear modulus, r , and the aggregate density, .
Furthermore, for the third parameter, , the denominator in the expression is usually larger than 0.8, so the ratio of and is not very sensitive to changes in the densities , c , or f . That means that neither is it very sensitive to changes in tortuosity and it will mainly be that depends on in the ratio ∕( B).
In this way, three of the seven parameters can be said to be redundant and the seven parameters from Table 1 ( , B, , f , s , , r ) may be reduced to four by combining the first three into one, ∕( B). Alternatively, the four parameters may be stated as ∕ , c , , and r . In case one wants to compute exactly, a fifth parameter, f , also needs to be included.
In this way, the number of parameters in the Biot shear model with some approximation is four, as in the Zener medium model, and five in the exact case.

Extension to poroviscoelasticity
In the sediment acoustics field, the Biot model has been amended in order to extend the model from a rigid porous frame, like a porous rock or bone, to one where the grains are allowed to move. This is a model for a saturated sediment and it is not unlikely that it may be more appropriate for tissue than the rigid frame implied in the Biot model.
In this case, viscosity is introduced in the solid frame in addition to the flowing liquid. The model is called the Biot squirt flow and viscous drag (BICSQS) model. 42 This viscosity is added by allowing a relaxation model for the shear modulus, thus allowing for a frequency-dependent complex shear modulus or dynamic shear modulus, described by a characteristic frequency : When the relaxation is included in Equations 37 and 11, the result is In the limit, as the frequency tends to zero the dynamic modulus approaches E(0) = and as the frequency approaches infinity it approaches E(∞) = ∞. This model is called the non-standard four-parameter model. 19 Its spring-damper realization is shown in Figure 7, with the equivalent conjugate model to the right. Therefore, even in this case an equivalent viscoelastic model may be found.
This model has been further developed to cover water-saturated sand as well. In that case, there is also coupling between the relaxation model for the bulk modulus and that for the shear modulus; see Equation 15 of Chotiros and Isakson. 43 The coupling term is a complex function of frequency that involves Bessel functions, but it can often be approximated with a simpler polynomial expression. 44 In that case, the model may also be expressed by springs and dampers, but an even more complex system than that of Figure 7. However, for biological tissue, the simpler model of Equation 39 should be adequate.

DISCUSSION AND CONCLUSION
The concept of time-and frequency-spectral decompositions of viscoelastic systems has been developed and applied to the fractional Zener and Kelvin-Voigt models. This summation of multiple elementary models is an idea that is independent of the fractional models. In fact, it is the idea behind the relaxation-spectrum models of Kelvin and Wiechert dating from 1888 and 1893, respectively. 15 In addition to the coverage by Mainardi 13 , Tschoegl's book 19 devotes several sections to it (Chapters 3.5-3.6). Similar ideas of model fitting have also been used in biomechanics and elastography. In the white matter model of Cheng and Bilston, 45 the shear relaxation modulus was, for instance, modeled with three very slow relaxation terms (0.01 Hz and slower) and a similar model with two terms was used by Caenen et al. 46 The idea behind the particular weightings implied in the time-spectral and frequency-spectral functions found here is, however, to make the sum approximate the behaviour of the fractional models, i.e. in terms of power laws in the frequency domain. Therefore it parallels the modeling of arbitrary power-law attenuation in medical ultrasound over a limited frequency range with a few terms by Yang and Cleveland. 47 The surprising result is that the weighting has a long-tailed distribution, with different power laws for low and high arguments. Here it characterizes the distribution of individual relaxation processes, in both time and frequency domains.
The poroelastic model depends on ten independent parameters, of which seven influence the shear-wave solution, and seems to build on an entirely different foundation from viscoelastic models. Despite this, it has been shown to be equivalent to a Zener model. This was done by comparing dispersion relations and, not unsurprisingly, the result is the same whether the dispersion relations from the original Biot theory or those from the Biot-Stoll theory are analyzed. This result is an extension of the work of Bardet,48 where it was shown that, in the low-loss/low-frequency approximation, an equivalent can be found between the poroelastic model and a Kelvin-Voigt model. This was done by comparing approximate expressions for the velocity and attenuation. As the Zener model in the low-loss/low-frequency case is equivalent to the Kelvin-Voigt model, 16 the result found here also agrees with that result.
It is also shown that the seven parameters of the poroelastic model can be reduced to five, and even four as in the Zener model, if a small approximation is allowed. When the Biot model is extended to include viscosity in the frame, as in the work of Chotiros and Isakson, 42 an extra damper has to be added to the Zener model, making it into what is called the non-standard four-parameter model. 19