Skip to main content

ORIGINAL RESEARCH article

Front. Phys., 21 January 2021
Sec. Medical Physics and Imaging
This article is part of the Research Topic Innovative Developments in Multi-Modality Elastography View all 20 articles

Poroelasticity as a Model of Soft Tissue Structure: Hydraulic Permeability Reconstruction for Magnetic Resonance Elastography in Silico

Damian R. Sowinski
Damian R. Sowinski1*Matthew D. J. McGarryMatthew D. J. McGarry1Elijah E. W. Van HoutenElijah E. W. Van Houten2Scott Gordon-WylieScott Gordon-Wylie1John B Weaver,,John B Weaver1,3,4Keith D. Paulsen,,Keith D. Paulsen1,3,5
  • 1Thayer School of Engineering, Dartmouth College, Hanover, NH, United States
  • 2Département de génie mécanique, Universite de Sherbrooke, Sherbrooke, QC, Canada
  • 3Geisel School of Medicine, Dartmouth College, Hanover, NH, United States
  • 4Dartmouth-Hitchcock Medical Center, Department of Radiology, Lebanon, NH, United States
  • 5Dartmouth-Hitchcock Medical Center, Center for Surgical Innovation, Lebanon, NH, United States

Magnetic Resonance Elastography allows noninvasive visualization of tissue mechanical properties by measuring the displacements resulting from applied stresses, and fitting a mechanical model. Poroelasticity naturally lends itself to describing tissue - a biphasic medium, consisting of both solid and fluid components. This article reviews the theory of poroelasticity, and shows that the spatial distribution of hydraulic permeability, the ease with which the solid matrix permits the flow of fluid under a pressure gradient, can be faithfully reconstructed without spatial priors in simulated environments. The paper describes an in-house MRE computational platform - a multi-mesh, finite element poroelastic solver coupled to an artificial epistemic agent capable of running Bayesian inference to reconstruct inhomogenous model mechanical property images from measured displacement fields. Building on prior work, the domain of convergence for inference is explored, showing that hydraulic permeabilities over several orders of magnitude can be reconstructed given very little prior knowledge of the true spatial distribution.

1 Introduction

In the past two decades magnetic resonance elastography (MRE) has emerged as an imaging modality that extends the capabilities of standard MRI, giving clinicians and researchers unprecedented access to tissue properties in vivo. As it develops rapidly, MRE promises clinicians information akin to the data surgeons acquire through manual palpation - material properties such as stiffness - but noninvasively and quantitatively, with applications in diagnosis and monitoring of diseases. Connecting macroscopic mechanical response to microscopic cellular organization of tissue, MRE opens up new avenues for assessing the etiology of tissue pathology. The modality has been implemented to varying degrees of success on skeletal muscle [19], breast [1017], liver [1824], and brain [25, 26].

Elastography dates back to Gierke et al’s 1952 paper, Physics of Vibrations in Living Tissue - one can discern tissue properties by observing a tissue’s response to vibration [27]. MRE is multidisciplinary at its core, requiring expertise in magnetic resonance imaging, biophysical modelling/theoretical mechanics, and numerical simulations/machine learning. Several reviews are available [2833].

The purpose of the current paper is to explore a poroelastic model of tissue, and convey novel results on the simultaneous imaging of both shear stiffness and hydraulic permeability in silico. These results extend prior studies on poroelastic parameter estimation [3437], and quantify the domain of valid inference. We show that convergence to the correct in-silico property values is possible even when prior knowledge is off by several orders of magnitude. These results are a foundation for ongoing in vitro and in vivo work, giving reasonable bounds for noiseless data.

This paper is organized as follows: Section 2 covers the theoretical details of our tissue model and how inference is accomplished. It also details some of the numerical considerations that have been implemented in our in-house MRE computational platform. Section 3 details the setup of our simulations - inference of one and two properties from in-silico data generated by solving Newton’s Laws for a poroelastic medium. Single and dual property contrast conditions necessary for fidelitous inference are examined. In Section 4 we discuss the results from our simulations and conclude with planned future directions for in vitro and in vivo application.

1.1 Notational Conventions

Throughout, we use coordinate independent tensor notation to emphasize the physical underpinnings of expressions, but, when required to use coordinates, follow the Einstein summation convention where repeated indices in terms are summed over. Indices are drawn from the middle of the Latin alphabet, i,j,k,, and run from 1 to 3. Scalar quantities are written in standard font with no typographical emphasis. Vectors and second rank tensors appear in lower case bolded font, the former in Latin and latter in Greek. Fourth rank tensors are bolded, and in capitalized Latin. When discussing either solid or fluid phases similar symbols are used to denote physical quantities, however the addition of a subscript s or f, respectively, is added to ensure physical clarity.

We denote the standard tensor product with a . An orthonormal vector basis is denoted {e^i}i=13, with e^ie^j=δij where the RHS is the Kronecker delta. At times we use the Cartesian representation of the basis, e^1=x^,e^2=y^, and e^3=z^. We construct the basis for second and fourth rank tensors from these, {e^ij=e^ie^j}i,j=13 and {e^ijk=e^ie^je^ke^}i,j,k,=13 respectively. The second rank identity tensor is denoted 1=e^ii. The dot product, , and double dot product, : , act on neighboring basis tensors via e^ie^j=δij and e^ij:e^k=δikδj, implying contraction on the component indices. For example, given an arbitrary vector, a, second rank tensor, β, and fourth rank tensor, C, we can construct a vector as follows:

aC:β=(aie^i)(Cjkme^jkm):(βnoe^no)=aiCjkmβnoe^ie^je^ke^e^m:e^ne^o=aiCjkmβnoδije^kδnδmo=aiCikmβme^k

The index gymnastics in evaluating such expressions can get pretty messy, which is why we stick to coordinate independent notation as much as possible. The transpose acts on second rank basis tensors as e^ijT=e^ji. The gradient operator is symbolized by a nabla, . Temporal Fourier transforms of quantities appear with the same typographical emphasis as above, having the added decor of a tilde. For temporal partial derivatives we employ fluxion notation, so that f˙ft, f¨2ft2, etc.

2 Modelling Tissue

Biological tissue is a complex medium, the result of almost four billion years of continual evolutionary pressure and the whims of chance. At the microscopic level, the variety of cellular structures in animals would require a vastly more complicated set of modelling assumptions than continuum mechanics can provide. However, tracing over these microscopic degrees of freedom, an effective theory of tissue should have more in common with a coupled solid-fluid description than one consisting of a single continuum. Such a model already exists, developed extensively during the last century by the geophysics community - poroelasticity.

2.1 Poroelasticity

Poroelasticity is an effective mechanical model describing the continuum as a coarse grained biphasic medium consisting of an elastic matrix coupled to a fluid. The former is not simply-connected, but has a connected pore space with a tortuous topology - it is this pore space that the fluid moves through. The coupled dynamics of solid and fluid result in a much richer phenomenology than single phase continuum models. Figure 1 shows a schematic of a porous structure.

FIGURE 1
www.frontiersin.org

FIGURE 1. The model of a porous material consists of an elastic matrix - the solid - and a connected pore space. Pores whose connected boundaries are homeomorphic (topologically identical) to spheres are not a part of the connected pore space, as nothing can flow through them. The connected pore space is filled with a fluid, whose interactions with the walls of the elastic skeleton provide a much richer phenomenology than elastic/viscoelastic models of materials. When a porous substance is intersected with a plane, the cross section, whose area is A, will have holes in it. The effective cross-sectional area, Af, is the area of all of these holes. For an isotropic material, the orientation of the intersecting plane does not change the ratio Af/A. The Delesse-Rosiwal law states that this ratio is equal to the ratio of the pore space volume to the total volume i.e. the porosity [40].

The theory of poroelasticity has its roots in the empirical work of Henry Philibert Gaspard Darcy, a water engineer tasked in 1856 with evaluating the many elaborate public fountains within the city of Dijon, France. Darcy, who had been working as a water engineer for nearly three decades, was pivoting towards pure research, and used this as an opportunity to write his 600 page magnum opus Les Fontaines Publiques de la Ville de Dijon, wherein he describes the sand column experiments pictured in Figure 2 [38]. By examining the flow of water through sand in cylindrical pipes, Darcy discovered a linear proportionality between the rate of fluid volume leaving the bottom of a pipe, ΔVf/Δt, and the difference in piezometric head at the ends of the pipe, Δh,

1AΔVfΔt=KΔhL,(1)

where L is the length of the pipe, and A its cross-sectional area. Hydraulic conductivity is the name given to the constant of proportionality, K, and for water flowing through sand was measured to range from 2.073.22×104 m/s [39]. This name is kept alive today by geologists, though we will see shortly that other fields’ naming conventions have created an aura of ambiguity. Note that K depends on the properties of the fluid used in the piezometers, as well as the porous medium and the fluid flowing through said medium - since geology is primarily concerned with the filtration of water through soils, rocks, and clays, they think of it as intrinsic to the medium. For us, biologically relevant fluids can differ substantially in density and viscosity, so we mustn’t fall into the same mindset. The empirical formula Eq. (1) came to be known as Darcy’s Law, a starting point for understanding how the flow of a fluid is affected by its coupling to a porous material.

FIGURE 2
www.frontiersin.org

FIGURE 2. The setup of one of the sand column experiments, from plate 24, Figure 3 of Darcy’s 1856 Les Fontaines Publiques de la Ville de Dijon [38]. Annotations have been added to clarify the origin of the physical quantities appearing in Eq. (1). The piezometric heads h1 and h2 in the two piezometers on the right are measured using mercury. L is the length of the sand column, and A its cross-sectional area. In this setup the diameter of the column is 0.35 m, and the height of the column is 3.5 m. ΔVf/Δt is the rate at which fluid volume is leaving the column, measured by observing the rise of the water level in the container. One can quickly see that this particular setup is problematic for applying Eq. (1) directly unless the sand column reaches up to near the top piezometer. If not, then the piezometric head must be replaced by the total head so as to includes the pressure head coming from the column of water.

A modern view of Eq. (1) utilizes the instantaneous rate of fluid volume moving through a cylinder at constant velocity, vf, replacing ΔVf/ΔtvfAf, where Af is the effective area seen by the fluid due to the microscopic geometry of the porous medium. The effective area is related to the porosity (volume fraction of fluid, assuming saturation) by the Delesse-Rosiwal law (second equality) [40], ϕ=Vf/V=Af/A. Identifying the pressure differential between the ends of the pipe with the drop in head via ΔP=ρfgΔh (ρf is the fluid density and g the gravitational acceleration), shrinking the pipe to infinitesimal length, L0, and restoring the tensorial nature of the physical quantities, we find

q=1ηκP(2)

where two clarifications are in order. (1) The LHS is the Darcy velocity (or filtration velocity), q=ϕvf, the magnitude of which is called the specific discharge. (2) The RHS is rewritten using the hydraulic permeability tensor, κ - a property belonging solely to the porous medium - and η - the dynamic viscosity of the fluid. The full justification of this separation is possible through a microscopic treatment using homogenization theory [4144], extending the validity of Darcy’s Law to pressure gradients not necessarily generated by gravitational forces.

For an isotropic medium, such as Darcy’s sand column, κ=k1, and we call k the hydraulic permeability. Unfortunately, as mentioned earlier, there is a bit of ambiguity within the literature depending on the scientific field one is enamoured with - geologists measure piezometric head and refer to Darcy’s original K [LT1] as hydraulic conductivity, while most other fields measure pressure gradients and refer to K¯=K/ρfg [L3TM1] as such. Furthermore, the medical and bio-engineering literature has, on occasion, used the terms hydraulic conductivity and permeability interchangeably; it is best to always check units to make sure that you and the author are on the same page. Given the dynamic viscosity of water is103 Pa s, this implies Darcy’s original measurements of sand’s hydraulic permeability on the order of 1011 m2. In honor of this, hydraulic permeability is now measured in units of Darcys, D, so that clean sand in aquifers comes out as 101 D. The unit is defined so that 1D permits a filtration velocity of a cm/s for a fluid with dynamic viscosity 103 Pa s under a pressure gradient of 1 atm/cm, or in SI, 1D=9.86923×1013 m2. The strangeness of this value comes from there being 101,325 Pa in 1atm, rather than a round 105, but who are we to quibble with history?

Karl von Terzaghi, a mechanical engineer, laid out the dynamics of structure consolidation on soil in a series of papers [4547] between 1922 and 1925, extending Darcy’s insights to become the father of soil mechanics. He incorporated the dynamics of soil in response to the fluid by introducing an effective stress tensor, σ, to describe the bulk medium,

σ=σsP1.(3)

Here σs is the stress tensor of the solid, and P is the pore pressure. The effective stress describes the net stress on the solid skeleton, and it’s definition is referred to as Terzaghi’s Principle - the total stress on the elastic matrix is lessened by the pore pressure. Terzaghi’s analysis describes a static, highly symmetric, and fully saturated system in which an incompressible fluid is confined to a one dimensional column of incompressible solid grains. Rearrangements of these grains drive deformation of this soil when a load is applied. The decomposition of the stress into an effective bulk and fluid is a bit like looking at the soil through poor lenses that smear the distinction between microscopic grains and pore fluid - a technique central to homogenization theory in engineering, and coarse graining in physics and chemistry [4850]. One should be a tad careful in their mind’s eye in all that follows - in going to this effective picture we now imagine solid and fluid components existing simultaneously at every spatial point, their dynamics coupled as a result of coarse graining the microscopic structure.

Though revolutionary in application, Terzaghi’s assumptions were ill suited for generality. Thermodynamically speaking, his model contains the pair of dynamical variables {σ,P}, the effective stress and the pore pressure, respectively; the elastic matrix strain,

ε=12(us+(us)T),(4)

is conjugate to the former, where us is the displacement field of the solid matrix. A good guess for the kinematic variable corresponding to the conjugate of the pressure would seemingly be the fluid strain rate, which for a perfect fluid with velocity field vf=u˙f is the single degree of freedom tr(ε˙f)=vf. Though this works to describe a fluid alone and results in the Navier-Stokes equations, this term cannot be coupled to the strain in a linear theory, nor to the solid displacement as that ruins translational invariance. The only option to preserve these necessary ingredients would be the fluid strain, uf. It turns out that this choice almost works; unfortunately it couples the solid and fluid a tad too weakly by predicting a vanishing pore pressure when the fluid moves in sync with the solid. With the obvious choices for a kinematic conjugate used up, it appears as if we must continue our search in darkness.

Maurice Antony Biot, an applied physicist, shed light on the correct coupling in 1941 by introducing [51] the variation in relative fluid content,

ζ=ϕ(ufus)+Γ,(5)

where γ=Γ˙ is referred to as a fluid source. His treatment generalized the Darcy velocity to a relative motion of the fluid with respect to the solid, q=ϕ(vfvs), thus allowing for an analysis of situations in which both the fluid and porous medium are in motion. Note that the relative fluid content decreases when net flow out of a dilating volume of the solid matrix occurs, and increases if the flow is reversed - together with the filtration velocity this is summarized by a continuity equation

ζ˙+q=γ.(6)

Here, we see the fluid source, γ=Γ˙, in action either creating fluid content, γ>0, or annihilating it, γ<0 - a rather useful quantity when considering tissue near veins and arteries.

In Terzaghi’s analysis, bulk volume change is driven by rearrangements of the pore space and flow of the fluid. Biot realized that this physics may be true for nearly incompressible sand grains, but for a softer elastic medium one should also take into account the contribution of solid compressibility to the bulk volume change. An applied isotropic load can now change the volume of the bulk by altering the volume of the solid skeleton due to the fluid pushing back - for reasonably small volume changes the effective support of the fluid with regards to the bulk is then less than the fluid pressure. Biot argued that the correct effective stress should be a modification of Eq. (3),

σ=σsαP1(7)

with the Biot effective stress coefficient, α, a function of the drained effective bulk modulus and the unjacketed solid bulk modulus, K and Ks respectively [52],

α=1KKs.(8)

To clarify, experimentally the unjacketed bulk modulus is measured in undrained conditions wherein a sample is submerged in fluid; the fluid pressure is then allowed to vary and the volume change of the sample is measured. A jacketed test, on the other hand, is done under drained conditions wherein a dry sample is wrapped in a membrane and connected to the atmosphere via a small hole; a load is then applied under constant atmospheric pressure, and the volume change measured. In aerated soil K/Ks1, validating Terzaghi’s treatment [53, 54]. The theory predicts a second type of longitudinal wave within a porous medium, which was observed in 1979 [55]. Biot refined his model by incorporating anisotropy, analyzing dispersion in wave propagation, and examining non-linear extensions of the constitutive relations, earning him the title of father of poroelasticity [5659].

Examining the effective potential can further illuminate the role of α - it is the dimensionless coupling between the volumetric strain of the skeleton and the variation in relative fluid content. Since it and ζ are both dimensionless, the Biot modulus, M, must be introduced to describe the energy density coming from the relative fluid content and its coupling to the elastic matrix, allowing the quadratic potential energy density for a linear poroelastic medium to be written as

U=12ε:E:ε+12Mζ2αMζtr(ε)(9)

where E=E0+α2M11 is an effective elastic modulus, a fourth rank tensor describing the effective coupling between the components of the strain. For an isotropic Hookean material described by Eq. (20), note that the Biot effective stress coefficient and modulus serve to renormalize the second Lamé paramter, λ. The conjugate variables can now be computed from this, giving the poroelastic constitutive relations

σ=Uε=E:εαMζ1(10)
P=Uζ=αMtr(ε)+Mζ(11)

Combining Eq. (10) and Eq. (11) recovers the Biot decomposition, Eq. (7). These expressions clarify the interpretation of M - it is the change in pore pressure caused by a variation in water content under the constraint of no solid dilation, M=P/ζ|tr(ε)=0. It can be written in terms of the Biot effective stress coefficient, as well as the fluid and jacketed solid moduli, Kf and Ks′′ respectively [44, 52],

1M=1K(α(1α)ϕKKs′′+ϕKKf).(12)

Biot, and others, explored the physical meaning of both α and M in a number of manuscripts [53, 60, 61]. The effective potential, Eq. (9), can also be constructed beginning at microscopic scales to justify the expression through homogenization theory, coarse graining, or mixed micro-macro formulations [4143, 6264].

We now use the effective stress and pressure to write the coupled equations of motion for the effective medium and relative fluid flow components in the absence of external stresses. The equations of motion can be derived via the principle of least action [42, 44, 62], however that would take us too far afield, so we introduce them qualitatively. Momentum conservation for the effective medium reads

(1ϕ)ρsu¨s+ϕρfu¨f=σ+(1ϕ)ρsg+ϕρfg,(13)

where the LHS is the total change in inertia of both solid and fluid. If we define the effective density as ρ=(1ϕ)ρs+ϕρf, the LHS can be written in a more illuminating form as ρu¨s+ρfϕ(u¨fu¨s) - the first term is the effective inertia following the acceleration of the solid component, while the second term is the relative acceleration of the fluid with respect to the solid. The RHS is simple to interpret - the divergence of the bulk stress tensor describes the internal forces, while the external forces consist of a gravitational field.

Similarly, momentum conservation for the fluid reads

ϕρfu¨f+ρa(u¨fu¨s)=ϕPϕ2ηk(u˙fu˙s)+ϕρfg.(14)

The LHS here is the change in fluid inertia; however a mass density, ρa, is added for the fluid moving relative to the solid. This addition can be derived from a micro-macro approach, and captures the effect of the fluid being slowed down through interactions with the structure of the pores [43, 44, 65]. It is typically written as ρa=(a1)ϕρf, with a for packed spheres being related to the porosity as a=1+ϕ2ϕ - for glass beads at low frequencies it has been measured as a1.66±0.13, or ϕ0.43±0.05, which coincides nicely with the porosity of loosely packed spheres, 0.41±0.04 [66, 67]. As in Eq. (13), the LHS of Eq. (14) is rewritten using the relative acceleration of the fluid to the solid as ϕρfu¨s+(ϕρf+ρa)(u¨fu¨s). The RHS has the pressure gradient driving the fluid acceleration, but also includes a viscous dissipation term coming from the relative motion of the fluid with respect to the solid, serving to synchronize the former to the latter. The porosity factor for the former is due to the fact that the pore pressure is acting on the effective area seen by the fluid. The ϕ2 factor on the latter can be broken down as follows - one factor comes from the dynamic viscosity, which is proportional to the fluid mass, while the second comes from the hydraulic permeability, which is proportional to the effective area seen by the fluid mass.

Note the inverse dependence on hydraulic permeability - the more easily the solid permits fluid flow the less dissipation occurs. Using the definition of the filtration velocity, Eq. (14), can be written as a generalization of Darcy’s Law, Eq. (2),

q=kη(P+ρfu¨s+ρa+ϕρfϕ2q˙ρfg),(15)

which now includes the inertial forces coming from the motion of the solid skeleton, as well as the fluid’s interactions with the pore space. Accordingly, we drop the subscript on the solid displacement, u=us.

The system of Eqs. (6, 13, 14) together with the constitutive relations Eqs. 10, 11 are a more suitable continuum model of tissues containing a mobile interstitial fluid component than standard viscoelastic continua. Since MRE subjects tissue to oscillatory driving conditions at constant frequency, the system can be simplified by Fourier transforming into frequency space - linearity assures us that distinct temporal modes evolve independently. The driving frequency, ω, and all of its harmonics, ωn=nω n+, contribute to the spatial signal. For each of the d.o.f./hyperparameters x(u,ζ,q,σ,P;g,γ), the Fourier decomposition and its inverse read

x(r,t)=nx˜n(r)eiωntx˜n(r)=1T0Tdtx(r,t)eiωnt.

The decomposition acts formally by replacing xx˜n and tiωn. Decomposing each d.o.f. in the equations of motion, we find the mechanical response at each frequency ωn is described by

iωnζ˜n=q˜nγ˜n(16)
ωn2ρu˜niωnρfq˜n=σ˜n+ρg˜n(17)
ωn2ρfu˜niωnρfβn1q˜n=P˜n+ρfg˜n(18)
 whereβn1=ρa+ϕρfρfϕ2+iηωnρfk.(19)

together with the transformed constitutive relations. Here, we have introduced the poroelastic β factor, a key ingredient in differentiating poroelastic and viscoelastic dynamics. In a static and constant gravitational field we of course have that g˜n=gδ0nz^.

We want to connect the MR images of tissue to this model, so we assume that the MRI phase contrast displacement measurements acquired by MRE trace the motion of the elastic skeleton. Since fluid gives a strong MR signal, this assumption relies on there being significant amounts of fluid in the disconnected pore space, acting as part of the bulk. Since fluid pressure within tissue is critical for homeostasis, we should reduce our system of equations further by eliminating all variables except the pair (u˜,P˜). For analytical simplicity, we choose an isotropic elastic stress tensor parameterized by the standard Lamé parameters μ (shear modulus) and λ,

E0=μe^ije^ij+μe^ije^ji+λ11(20)

After some lengthy algebra, one finds a vector and scalar family of equations for each harmonic,

ωn2(ρβnρf)u˜n=[μu˜n+μ(u˜n)T]+(λu˜nαP˜n)+βnP˜n+(ρβnρf)g˜n(21)
αu˜n[βnu˜n]=1ρfωn2[βnP˜n]1MPn˜+1ωn2[βng˜n]+iωnγ˜n.(22)

Up to this point we have not imposed assumptions on the spatial distribution of the elastic moduli in the constitutive relations, or the hydraulic permeability. We do so now by assuming that the underlying solid and fluid are nearly incompressible relative to the bulk, namely that K/Ks1, K/Ks′′1, and K/Kf1 [68]. Plugging the first of these in Eq. (8) shows that in this regime α1, while using all three in Eq. (12) shows that M10. In this physical regime, the above reduce to [52]

ωn2(ρβnρf)u˜n=[μu˜n+μ(u˜n)T]+(λu˜n)(1βn)P˜n+(ρβnρf)g˜n(23)
iωnρfγ˜n=[(1βn)ρfωn2u˜n+βnP˜nβnρfg˜n].(24)

The vector equation shares some resemblance to viscoelasticity - the main differences being the presence of the complex βn and displacement driving pressure gradient. Both Eq. (13) and Eq. (14) can be rewritten using the filtration velocity, getting rid of all mentions of the fluid displacement. The scalar equation is a divergence condition, implying that the given combination of solid displacement and pressure gradient can only contribute a net curl to the system unless a fluid source is present. Solving this system with given boundary conditions is referred to as the forward problem.

2.1.1 Numerical Considerations I

Analytical solutions of the equations of motion exist for highly idealized geometries [41, 44, 62]; in the context of classical clinical MRE such idealizations do not exist, hence computational methods are required to find approximate numerical solutions [69]. We use a finite element forward solver, described in detail elsewhere [35, 37]. It exploits multiple meshes - displacement fields can be supported on either tetrahedral or hexahedral element meshes, while material properties, such as the shear modulus or hydraulic permeability, are supported on separate hexahedral meshes [70]. Supporting each property on its own mesh provides versatility. It decouples the resolution of the unknown properties being inferred from the mesh resolution of the (u,P) solutions generated by the forward problem, which balances discretization error and computational load without impacting the spatial scale of parameter inference. As we move forward into work on phantoms and human subjects, this versatility will allow for mitigating the effects of noise through coarser property meshes.

2.2 Inference

In MRE we observe the displacement field of the tissue being examined with little or no knowledge of the underlying mechanical properties of the tissue. This inverse problem of using the observation to infer the properties, is ideal for the application of a powerful epistemological tool known as Bayesian inference.

The Bayesian interpretation of probability brings with it many tools, some of which have started to be used by the MRE community [71, 72]. The fundamental idea behind it is simple and intuitive - observation carries with it information, and information leads to a reduction in uncertainty [7377]. Beliefs evolve in order to reduce uncertainty, thereby maximizing the fitness of a cognizant being within an environment that is sufficiently complex to be indistinguishable from one behaving pseudorandomly. Note that, in a Bayesian setting, belief/plausibility are synonymous with probability, and in what follows the terms will be used interchangeably [78, 79].

Applying these ideas to MRE, let us fix a particular mechanical model of tissue so that the forward problem is specified. The space of all possible material property fields is denoted ={θ}. An epsitemic agent has a preexisting, or prior, belief distribution over , namely p(θ)[80]. If they know absolutely nothing about what to expect then the prior is uniform over all possible property fields - a uniform distribution is maximally uncertain [81]. Of course, if they have information that may reduce uncertainty in some way, the prior could be peaked over a region of that is deemed more likely to occur. The space of all possible observations of the tissue displacement field is denoted D={u}. Finally, the posterior over is the informed belief distribution that the epistemic agent has after observing the data, denoted p(θ|u). The information content of a particular belief is defined uniquely [73, 82] in order to satisfy a natural sub-additive constraint on dependent beliefs, I[θ]=logp(θ). Note that strong beliefs entail little information content, while implausible beliefs the opposite. Observing a rare event endows one with more information than observing a common event. It is for this reason that sometimes information content is referred to as surprise - we will be content in interpreting information content as a measure of the implausibility of a belief.

The posterior and prior are connected via bayes’ theorem

p(θ|u)=p(u|θ)p(u)p(θ)(u,θ)p(θ)(25)

where the quantity multiplying the prior is the likelihood of the observation given a model. Note that observations that are more likely in the context of a set of material properties, θ, have (u,θ)>1 and enhance the posterior over the prior, whereas unlikely observations, (u,θ)<1, reduce the belief in that particular θ.

Specifying a poroelastic model with parameters θ={μ(r),k(r)}, along with the hyperparameters {ϕ,ρs,ρf,ρa,η,λ} and boundary conditions, one solves the equations of motion Eqs. (23, 24) for u˜n(r). As will be explained shortly, we need only do this for the fundamental frequency, so we define the model displacement as uM(r,θ)=u˜1(r). We have left the model as an argument here to remind us that the solution for the fundamental mode depends on the parameters that have been chosen.

MRE data acquisition images steady state vibration fields using phase sensitive MRI sequences [28, 8387]. Displacements are measured at a number of times across the harmonic motion cycle, which are processed into a sequence of displacements at each voxel and the fundamental mode of the time series is computed to give uD(r)eiωt for each voxel [88, 89]. Vibrations are often in the range of 25100 Hz and supplied using an external actuator, or alternatively, our intrinsic actuation MRE uses retrospective gating to measure cardiac induced motions at 1 Hz without the need for actuation hardware [90].

With uD and uM in hand, we are much closer to understanding the likelihood. The latter are encoding the material property degrees of freedom, while the former provide constraints. Both are discretized so one must always make sure to interpolate to the same voxel domain, taking into account that there are enough constraints to fully fix the degrees of freedom. Once we have fixed the voxel domain, we can treat the scalar difference between the two as a random variable Δui(θ)=|uD(ri)uM(ri,θ)|. The differences are independent across voxels, though this assumption can be relaxed when the noise structure of the data collection can be properly described. The Δui are at the voxel level, so each one is actually the mean of all the degrees of freedom at the subvoxel level, of which there are a great many. Invoking the Central Limit Theorem allows us to claim that the distribution of Δui is Gaussian for each voxel such that

(u,θ)p(u|θ)=ip(Δui)exp(12iΔui2σi2),(26)

with the σi quantifying the root mean square fluctuations in the subvoxel degrees of freedom of the ith voxel. Since no estimate is available for these values, they are set to unity. Using Bayes’ Rule, Eq. (25)

I[θ|u˜]=log(u˜|θ)+I[θ]12i|uD(ri)uM(ri,θ)|2+I[θ],(27)

where the in the second line means equality modulo a constant independent of θ, leads to Eq. (27) as a cost function; the first term on the RHS is an error function, while the second term, is a regularizer. Finding the set θ that minimizes cost is equivalent to having found the most believable set of parameters - an epistemological interpretation grounded in the tools of information theory.

2.2.1 Numerical Considerations II

There are, of course, multiple ways to attack this particular question, but they all have in common the theme of sliding down the cost surface. If the cost surface is convex, such a slide will eventually get us to the most plausible set of parameters. The idea can be broken down as follows: Let’s say that we are at some arbitrary model θ. Our goal will be to move to a new model θ+δθ so that the information content at this new model is: I[θ+δθ|u]I[θ|u]. We are content to slide either down or horizontally, but never back up the error surface. By an appropriate sequence of choices for the δθ’s one has an iterative procedure for finding the minimum. Different methods rely on expanding the above information content inequality to different orders in δθ, and then making appropriate approximations with the resulting gradients.

Our MRE platform has several of these methods implemented - descent direction can be computed via Gauss-Newton, Conjugate Gradient, or Quasi-Newton updates [37]. Travel distance is modulated by a binary line search once the descent direction is found. Computation of gradients depends not only on the log likelihood, but also on prior information content - the regularizer. Multiple regularization terms have been implemented to take into account prior knowledge on material property limits, their spatial distribution, and smoothness [35, 91, 92].

To manage the computational resources required for multiple solutions of the 3D finite element forward problem required for inference, an efficient zoning strategy has been developed [93, 94] that breaks up the domain of interest into smaller overlapping zones. Measured displacements on the boundary voxels of zones provide Dirichlet boundary conditions, leading to a similar inference problem on a smaller scale. Computational complexity in our MRE-suite for the finite element solution of Eqs. (23) and (24) on a cubic zone of side length L scales as O(L4.6), which is less than the O(L6) complexity of LU factorization due to the sparsity resulting from using local operators [95]. Each zone is treated in parallel, computing as many gradient descent iterations as desired. Property solutions from each zone are collected and stitched together to complete a global iteration [70]. Each global iteration begins with a reseeding of the zone structure to limit propagation of bias or boundary noise inherent in a particular zone throughout the inference.

Convergence of any combination of the above methods is tricky to gauge due to the immense dimensionality of . Furthermore, the use of zoning means each global iteration is performing inference on a slightly different cost surface. To understand convergence we have therefore opted towards a coarse grained approach, projecting the inferred d.o.f. to a lower dimensional subspace consisting of the mean value(s) of the inferred property descriptors. Inference is judged by cost function descent on this subspace. A geometric interpretation of the method is given in Figure 3.

FIGURE 3
www.frontiersin.org

FIGURE 3. A geometric view of the cost surface’s lower dimensional representation. The trajectories in D and (red, green, and blue) are meant to represent how difficult it is to visualize the paths in these high dimensional spaces. In MRE, the dimension of both is proportional to the number of voxels. To visualize the process of inference, we map our trajectories to a lower dimensional space made of the negative log posterior (root mean square of the displacement errors), and the mean property value. The RMS is proportional to the L2 distance between the data displacements, uD, and model displacements, uM. The mean property value is proportional to the L1 distance from the “no property” origin to the model point, θ. During single property inference the resulting lower dimensional representation is 1-dimensional, while for two property inference it is 2-dimensional. Note how the red and green trajectories eventually leads to the minimal RMS. The blue trajectory begins at a high enough error (near vanishing posterior), that our epistemic agent is incapable of finding the correct direction to move in . The dimensionality of is so large that nearly every direction leads to an increase in RMS displacement error, and the lower dimensional representation interprets this as rolling up the cost surface, leading to divergence.

3 Simulation

We use our MRE computational platform - a finite element code capable of simulating displacement fields in poroelastic materials, and an artificial epistemic agent running Bayesian inference on measured displacement fields to extract poroelastic properties. The former situation is referred to as the forward problem, MREf, while the latter the inverse problem, MREi.

The following tests examine MREi’s ability to reconstruct spatial maps of both the hydraulic permeability, k, and the real shear modulus, μ. First we examined our code’s baseline ability to infer distributions of a single property at a time. Next we considered the more complex requirement of inferring both material property distributions simultaneously. Simulated data were generated via MREf, followed by an exploration of the parameter space to see how initial conditions affect material property inference and find a range of starting conditions which converge to the known ground truth.

3.1 Single Property Inference

3.1.1 Forward Problem

Reconstructing the shear modulus or the hydraulic permeability one at a time was tested through configurations depicted in Figure 4(10 Hz) and Figure 7(1 Hz) - a long rectangular block containing two conical inclusions. For both frequencies we tested two distinct scenarios. In the first scenario, inclusions were assigned shear modulus contrast. Relative to the background shear modulus μ0, the inclusion on the right has a shear modulus of μ=χμ0, where χ>1, whereas the one on the left has a shear modulus of μ=χ1μ0. The shear modulus contrast of the right inclusion is positive (μμ0)/μ0=χ1>0, while the one on the left is negative, (μμ0)/μ0=(χ1)/χ<0. In the second scenario the inclusions were given hydraulic permeability contrast. We ran simulations with χ ranging from 1.01 to 10, with similar results, and here, we report on the case χ=2.0, in which the inclusions have material property contrasts 1/2 and 1 - in the range reported for in vivo soft tissue MRE.

FIGURE 4
www.frontiersin.org

FIGURE 4. (A) Boundary conditions on the top/sides(I) and bottom(II) surfaces for the 10 Hz simulation. n^ are outward pointing normals. Pressure has Dirichlet BCs on I and Neumann BCs on II. Displacement has the opposite. The bottom plate oscillates in the y direction at frequency ω. (B) Orthographic projection showing geometry of the conical inclusions. The left inclusion has positive property contrast, and the right one has a negative property contrast. (C) The dynamics of both scenarios - the left figure has inclusions with shear modulus contrast, while the right figure has inclusions with hydraulic permeability contrast. The MREf displacement field over one period is shown for both, the red dots representing in-phase displacements at the start of a cycle. Below these, the motion of the boundary nodes is displayed at four equally spaced times within the cycle of period T.

Referencing Figure 4 and Figure 7, we chose r=4.25 mm, making the block dimensions 1.7cm×3.2cm×1.7cm. The displacements were supported on a tetrahedral mesh with 10,081 vertices and 55,593 elements. Although this is smaller than typical organs imaged with MRE in practice, results were similar to larger values of r, and the computational cost of running a large number of simulations favored smaller geometries for detailed numerical investigations. Material properties were supported on a 0.5 mm cubic lattice mesh with 71,672 vertices and 66,429 elements. The conical inclusions each have a relative volume of Vc/V=π(233)/324.56%. The displacement mesh intersects the inclusions at 453 and 461 vertices for the left and right cones, respectively.

We report results the results for both frequencies, f=ω/2π=1 Hz and 10 Hz. For the shear modulus frequencies between 1100 Hz were tested with similar results. For hydraulic permeability we evaluated a smaller range of frequencies, finding that the reconstruction worsened at frequencies above 10 Hz. At 10 Hz, see Figure 4, fixed displacement boundary conditions were imposed on the bottom face of the block to generate oscillations parallel to its long axis, u=u0eiωty^. At 1 Hz, see Figure 7, oscillations perpendicular to the bottom face were imposed, u=u0eiωtz^ The displacement amplitude in both cases was chosen to be 0.5mm, justifying a linear treatment, u0/L0.031.

Material property values used in MREf are summarized in Table 1. Bulk density was chosen to be similar to brain tissue, while fluid density and dynamic viscosity were selected to fall in the range of cerebrospinal fluid [96100]. The apparent density is a coarse grained effective quantity coming from the interactions between the fluid and the pore space, and we chose a value based on prior studies, noting that a low sensitivity of results to ρa was demonstrated [3436, 68]. Lamé constants have been well studied in MRE, our choices reflecting reasonable values [92]. Porosity was chosen based on mean values for white and grey matter using DTI, as was the background hydraulic permeability [101]. It should be noted that these last two properties are not well measured, with the literature quoting values in a wide range [102105].

3.1.2 Inverse Problem

Inference uses a coarser material property mesh in order to ensure enough constraints to solve the equations of motion. A side effect of the forward and inverse meshes not being identical is that there is no zero error solution. Here we have a 2.0 mm cubic lattice mesh, with 1700 vertices and 1296 volume elements.

We explore the question of how close to the true mean our initial values have to be in order to ensure convergence by examining ranges spanning 2 orders of magnitude above and below the true mean. The background shear modulus is set to μ0=3.0 kPa, and the background hydraulic permeability is k0=101.3 D. Thus, the ground truth mean values for χ=2 are:

μGT=(1+VcV(χ1)2χ)μ03068PakGT=(1+VcV(χ1)2χ)k0103.6D

These expressions are, of course, only approximate for a discrete mesh - deviations on the order of 1% were observed for our constructed meshes. Since these values are so close to the background fields due to the tiny volume occupied by the inclusions, we considered an inversion successful if the measured mean property fell near the region between the background and mean property values.

We initialize inference with a homogeneous material property distribution of either μ or k. Twenty (20) values of shear modulus were chosen ranging an order of magnitude above and below the true mean, 15030,000 Pa. Twenty (20) values of hydraulic permeability were chosen ranging nearly two orders of magnitude above and below the true mean, 100104 D.

The inversion used zones with 10% overlap. Zones were large enough that approximately 500 non-boundary nodes existed in each zone, with about 60 nodes in each zone overlapping with other zones. Gaussian smoothing regularization (spatial filtering) was applied, with the smoothing radius shrinking from 5 to 0.5 mm over the first 50 iterations and held constant for the remaining iterations. A conjugate gradient descent method was used for iterative inference on each zone. During the first 20 global iterations, each zone finds a single CG direction, and the distance moved is computed through a single iteration of a secant linesearch. The next 30 iterations use two of each, and the remaining use three.

Results of inference are summarized for the 10 Hz data in Figure 5 for μ and Figure 6 for k. For the 1 Hz data, the shear modulus results are shown in Figure 8 and the hydraulic permeability in Figure 9.

FIGURE 5
www.frontiersin.org

FIGURE 5. Reconstruction evolution of the shear modulus in the 10 Hz simulation. (A) shows the reconstructed images of the shear modulus distribution inside the simulated phantom. The column images are averages of the numbered (xy)-slices (top row), each of which is 2 mm thick. The rows are labeled by the iteration number of the inference. (B) shows the lower dimensional trajectory of the shear modulus. As above, the shear modulus is measured in kilo-Pascals. The smaller diagram is a magnified region near the converged values. Initial estimates are noted by an asterisk and changes in descent strategy are indicated by small and large red markers, as explained in the text. The light red region represents the strip between the background value of the shear modulus and the ground truth mean value. The thicker trajectory is the one used to make the images in the above panel. Due to a rather long convergence time, best fit linear extrapolations of the trajectories are plotted to show that the convergence region is (nearly) the same for all converging runs.

FIGURE 6
www.frontiersin.org

FIGURE 6. Reconstruction evolution of the hydraulic permeability, k, at 10 Hz. (A) shows the reconstructed images of the hydraulic permeability distribution. The column images are the averages of the numbered (xy)-planes (top row). The rows are labeled by the iteration number of the inference. (B) shows the lower dimensional trajectory of the shear modulus. As above, the hydraulic permeability is measured in Darcy. The smaller diagram is a magnified region near the converged values. Initial estimates are noted by an asterisk and changes in descent strategy are indicated by small and large red markers, as explained in the text. The light red region represents the strip between the background value of the shear modulus and the ground truth mean value. The thicker trajectory is the one used to make the images in the above panel. One can see that the convergence of the permeability occurs well before the final change in descent strategy. It appears that proper inference is robust to initial conditions 2 orders of magnitude greater than the ground truth, while only 1 order of magnitude below.

FIGURE 7
www.frontiersin.org

FIGURE 7. (A) Boundary conditions on the top(I), side(II), and bottom(III) surfaces for the 1 Hz simulation. n^ are outward pointing normals. Pressure has Neumann BCs on I & III, and Dirichlet BCs on II. Displacement has the opposite. Top plate is held fixed, while the bottom oscillates in the z direction at frequency ω. (B) Orthographic projection showing geometry of the conical inclusions. The left inclusion has positive property contrast, and the right one has a negative property contrast. (C) The dynamics of both scenarios - the left figure has inclusions with shear modulus contrast, while the right figure has inclusions with hydraulic permeability contrast. MREf displacement field over one period is shown for both, the red dots representing in-phase displacements at the start of a cycle. Below these, the motion of the boundary nodes is displayed at four equally spaced times within the cycle of period T.

FIGURE 8
www.frontiersin.org

FIGURE 8. Reconstruction evolution of the shear modulus, μ, at 1 Hz. (A) shows the reconstructed images of the shear modulus distribution inside the simulated phantom. The column images are the averages of the numbered (xy)-planes (top row). The rows are labeled by the iteration number of the inference. (B) shows the lower dimensional trajectory of the shear modulus. As above, the shear modulus is measured in kilo-Pascals. The smaller diagram is a magnified region near the converged values. Initial estimates are noted by an asterisk and changes in descent strategy are indicated by small and large red markers, as explained in the text. The light red region represents the strip between the background value of the shear modulus and the ground truth mean value. The thicker trajectory is the one used to make the images in the above panel. Once again we see that convergence of shear modulus is much slower than hydraulic permeability, yet images for a wide range of these simulated phantoms all show the correct structure of inclusions. Best fit lines are added to the zoomed in image to show that all the trajectories are approaching the same value.

FIGURE 9
www.frontiersin.org

FIGURE 9. Reconstruction evolution of the hydraulic permeability, k, at 1 Hz. (A) shows the reconstructed images of the hydraulic permeability distribution. The column images are the averages of the numbered (xy)-planes (top row). The rows are labeled by the iteration number of the inference. (B) shows the lower dimensional trajectory of the shear modulus. As above, the hydraulic permeability is measured in Darcy. The smaller diagram is a magnified region near the converged values. Initial estimates are noted by an asterisk and changes in descent strategy are indicated by small and large red markers, as explained in the text. The light red region represents the strip between the background value of the shear modulus and the ground truth mean value. The thicker trajectory is the one used to make the images in the above panel. Once again we witness that hydraulic permeability inference is robust to initial conditions over two orders of magnitude above the ground truth. With this simulation we see that the lower bound has gotten better as well, robust out to 2 orders of magnitude below the ground truth.

3.2 Double Property Inference

3.2.1 Forward Problem

MREf simulation of combined shear modulus and hydraulic permeability contrast involved the configuration depicted in Figure 10 - a square block containing eight inclusions. The inclusions are spherically symmetric, and have material property multipliers of either 1,χ, or χ1, resulting in either null, positive, or negative material property contrasts, respectively. Corner inclusions had both properties contrasted, while the center inclusions have contrast in one. This setup tests all possible pairs of positive, negative, or null contrasts over all the inclusions. Multiple values of χ were tested ranging from 1.13.0, with similar results - we report the case χ=2.0, corresponding to inclusions with property contrasts of either 1 or −1/2.

FIGURE 10
www.frontiersin.org

FIGURE 10. (A) Boundary conditions on the top(I), side(II), and bottom(III) surfaces for the two property inference simulation at 1 Hz. n^ are outward pointing normals. Pressure has Neumann BCs on I & III, and Dirichlet BCs on II. Displacement has the opposite. Top plate is held fixed, while the bottom oscillates in the z direction at frequency ω. (B) The motion of the boundary nodes is displayed at four equally spaced times within the cycle of period T. One can see the Poisson effect during the apexes of the squeeze/stretch cycles. (C) Orthographic projection showing the geometry of the eight spherical inclusions. The left inclusions have positive shear modulus contrast, while the right ones has negative contrast. The top inclusions have positive hydraulic permeability contrast, while the bottom ones have negative. The MREf displacement field over one period is superimposed, the red dots representing in-phase displacements at the start of a cycle.

Referring to Figure 10, We chose r=3.5 mm, making the block dimensions 3.5cm×3.5cm×1.4cm. MRPEf used a 1 mm tetrahedral displacement mesh (23,275 vertices and 132,084 volume elements) and a 1.0 mm cubic lattice material property mesh (23,104 vertices and 20,535 volume elements). Each inclusion took up π/3001% of the total volume, containing240 vertices and 215 elements - deviations from inclusion smoothness occur at 0.35 steradians.

Boundary conditions were imposed to simulate the top face of the block being held fixed, while the bottom face was vibrated in the vertical direction, u=u0eiωtz^. The amplitude was chosen to be u0=0.5 mm, ensuring that the size of vibrations relative to block geometry is small, u0/H=1/280.041, so that a linear treatment is justified. A frequency of f=ω/2π=1 Hz was applied. As in the case of single property inference, material properties were chosen to reflect known values associated with brain tissue. Our choices are summarized in Table 1.

TABLE 1
www.frontiersin.org

TABLE 1. The first six material properties were used in both MREf and MREi. The last two are background values used only in MREf. Inclusions have multipliers relative to background values as discussed in the text. In MREi, the last two properties, shear modulus and hydraulic permeability, are inferred not specified.

3.2.2 Inverse Problem

MREi deployed a coarser material property mesh than MREf - a 2.5 mm cubic mesh with 1792 vertices and 1350 elements. Each inclusion sampled 12–13 material mesh nodes.

We examined the effects of initial conditions on reconstruction fidelity by applying 20 different homogeneous initial conditions around the mean values, with the shear modulus ranging between 9.13×1025.15×103 Pa, and hydraulic permeability ranging between 8.9×1022.7×105 D. Background shear modulus and hydraulic permeability were μ0=3kPa and k0=101.2 D. Ground truth mean values were

μGT=(1+π100(χ1)2χ)μ03,047PakGT=(1+π100(χ1)2χ)k0102.9D

Once again, these were approximate values on the discrete mesh, but fell within 2% of the observed values. Inference was considered successful if both mean property values fell within or near the region [k0,kGT]×[μ0,μGT].

The first 10 global iterations used a single CG direction and a single secant line search for computing total distance traversed in parameter space while the next 30 applied two of each, and the final 250 invoked three of each. The first 50 iterations used Gaussian smoothing regularization (spatial filtering) on both properties, with the Gaussian standard deviation shrinking linearly from 5 to 0.5 mm. If the change in the cost fell below 1012, inference was considered complete, and further iterations ceased.

Evolution of image reconstruction is summarized in Figure 11 and the lower dimensional representation of the inference trajectory is in Figure 12.

FIGURE 11
www.frontiersin.org

FIGURE 11. Simultaneous inference evolution of both (A) shear modulus, μ, and (B) hydraulic permeability, k, in the 1 Hz experiment setup in Figure 10. The former is measured in kilo-Pascal, while the latter in Darcy. Columns are averages of the numbered (xy)-planes as in the other simulations. Rows are labeled by which iteration they are in MREi.

FIGURE 12
www.frontiersin.org

FIGURE 12. Lower dimensional representation of the inference trajectories. Note that because two properties are being reconstructed, the lower dimensional representation is 2-dimensional. This plot conveys the trajectories through an orthographic projection. (A) shows a large region of (k,μ)-space, including starting points of the trajectories represented by asterisks. Right and top panels represent root mean square values for converging trajectories. Some diverging trajectories are also shown - they are dashed and lighter. The single thick trajectory is the one used for Figure 11. Small and large red markers represent points where the descent strategy changed. The light red horizontal and vertical regions are the critical strips laying between the background property values and the ground truth mean values. Note that, as before, scales are all logarithmic. (B) is zoomed into the blue squares located in the figures described above. The shaded red areas represent the region between the background property value, μ0 or k0, and the mean value of the ground truth distributions, μGT or kGT. Note that in these plots the scales have changed to be linear. Discussion of trends is in the text.

4 Discussion and Conclusion

The dynamical evolution of our MREi code for single property inference and image reconstruction is displayed in Figure 6 and Figure 7 for f=10 Hz data. The top panels of both display reconstructed images at four distinct iterations. By the 20th iteration, both hydraulic permeability and shear modulus are not very distinct, though the former is clearly partitioning the domain into positive and negative contrast regions. By the 50th iteration, both material properties show distinct localization that is similar to the actual inclusions. At this stage, spatial filtering has reached sub-voxel width; hence, fuzziness of the hydraulic permeability image reconstruction relative to the sharper shear modulus image reconstruction is solely due to the geometry of the information content surface.

The f=1.0 Hz reconstructions are summarized in Figure 8 and Figure 9. Hydraulic permeability starts to settle to the correct values by the 50th iteration, though nearly 100 more iterations are needed for the correct shape of the inclusions to appear in part. This is due to the coarseness of the MREi mesh. The k value enters the critical strip typically within the first 20 iterations. We do note that when the initial values of the hydraulic permeability fall two orders of magnitude below ground truth values, MREi begins to diverge. Similar difficulties were not found if initial values overestimated the true hydraulic permeability. The shear modulus did not show this behavior, with image reconstruction recovering inclusion geometry and relative distribution very well by the 50th iteration, although different initial conditions are seen to lead to different final absolute values - a peculiarity noted in prior work [106]. The lower panels of Figures 5, 6, 8, 9 show the inference trajectories, with initial estimates noted by an asterisk and changes in descent strategy indicated by red markers.

Projecting onto the mean field error surface, hydraulic permeability approaches the correct mean field much earlier than shear modulus. Initial values do not affect the late stage reconstruction of the former when they are 2 order of magnitude too small or up to 3 orders of magnitude too high. The CG descent method does appear to have pathological behavior when initial values are more than 2 order of magnitude too small, resulting in an ascent of the error surface. A possible reason for this outcome is that low values of hydraulic permeability induce incompressible elastic behavior known to have stability issues in viscoelastic simulations [37]. The geometry of shear modulus reconstruction occurs with high fidelity within very few iterations, but the values converge much slower to the ground truth, an observation reported in prior work [106]. Projected trajectories of the error evolution in both Figure 5 and Figure 8, however, converge at or near the critical strip for all initial conditions.

The results of the simultaneous property convergence appear in Figure 11 and Figure 12. The former shows the fidelity of image reconstruction, while the latter inference trajectories as material properties are updated at each global iteration. Interestingly, the hydraulic permeability enters its critical strip well before the shear modulus. Mean value convergence, however, is not a strong indicator of image fidelity. The shear modulus inclusion geometry is discernible by the 75th iteration, whereas nearly twice as many iterations are required before hydraulic permeability inclusion geometry begins to resemble ground truth. Furthermore, dependency between the two properties can be observed near the edges of the simulation - hydraulic permeability reconstruction has some vertical patterning that resembles the shear inclusions.

Interestingly, we find that the k-reconstruction trajectories always approach the critical strip from above. Inference tends to overestimate hydraulic permeability with several iterations, and then drives it down. As in the case of single property inference, k initial conditions have a lower limit of about 1 order of magnitude below the ground truth in order for MREi to converge. Shear modulus is underestimated in almost all cases within about 50 iterations, and is then slowly driven upwards to the critical strip, not unlike the inference dynamics for the single property (μ) case.

Prior work [36, 107] has performed similar dual property reconstructions. These results, however, had MREi initial conditions very close to ground truth values, utilized the same property mesh for MREi as used for MREf, and required large amount of regularized smoothing. Our work shows that the first two assumptions can be discarded, and the third relaxed significantly, while maintaining high fidelity for image reconstruction.

There seems to be some confusion in the literature between the terms hydraulic conductivity and hydraulic permeability, with many authors using the two interchangeably. We focused early in Section 2 on making a clear distinction between the two in order to emphasize the specific contributions to fluid flow coming from the viscosity of the fluid and the geometry of the solid. Cerebrospinal fluid has a viscosity ranging between 0.71.0 mPa s [99], blood between 3.067.7 mPa s [108], ascitic fluid between 0.51.5 mPa s [109] - the range covers two orders of magnitude so tissues with similar hydraulic permeabilities can have vastly different hydraulic conductivities. Hydraulic permeability measurements correlated with pathology can therefore be connected to changes in the elastic matrix of the tissue (via its effective continuum).

Our motivation for simulated values of porosity and hydraulic permeability rests on a general lack of knowledge by the scientific community of in vivo values of these parameters. In vitro experiments tend to find small values for both, which we hypothesize is due to pore collapse in sampled tissue. It is easy to imagine that once tissue is removed from the body, the pressure change and unjacketing of the sample results in fluid being expelled. The drop in pore pressure is such that the stresses in the elastic matrix are incapable of supporting the pore space under the influence of gravity, leading to collapse. Once pores are closed, surface to surface cellular interactions may prevent the process from being reversible. Even simple models like the Carman-Kozeny [110, 111] relationship between hydraulic permeability and porosity, kϕ3, would predict drastic changes in the former under reasonable changes in the latter. Furthermore, preliminary work shows that higher values of hydraulic permeability result in higher likelihood, giving us confidence that in vitro measurements of permeability are being underestimated by up to several orders of magnitude. A full analysis is forthcoming.

The shear modulus, on the other hand, is much better understood, with far less variance in the literature. In viscoelastic models, white matter typically falls in the range of 1.52.5 kPa, depending on frequency, while gray matter is between 0.751.5 kPa [107, 112]. Viscoelastic estimates lump together resistance to deformation from both the solid and fluid constituents. The poroelastic shear modulus is the value for the drained solid, so will likely be somewhat lower than viscoelastic approximations, although a similar order of magnitude for most reasonably solid tissues. The frequency dependence of both is interesting, since it means that there is a memory effect in the tissue dynamics that has not been accounted for in the equations of motion. In order for this memory effect to be causal, in the sense that only past dynamics and not future dynamics affect the present state, certain constraints known as Kramers-Kronig relations must be satisfied by real and imaginary components of the moduli [113]. There is still some controversy over whether these conditions are satisfied by moduli measured in vivo. Finding a frequency dependence of the poroelastic parameters would put similar constraints on the equations of motion.

The introduction of βn in the poroelastic equations of motion creates a link between poroelastic and viscoelastic dynamics by means of the imaginary component of the parameter. Past numerical studies attempted to bridge the connection between the two types of behavior [107]. A full analytical treatment is in the works - understanding where the dynamics are most similar and where they are most different is a way of gaining insight into the microstructure of tissue. It also points towards better experimental design, since understanding the parameter spaces of the two can be harnessed to accentuate, and hence test for, behavior intrinsic to only one type of model. Forthcoming work will focus on understanding this connection.

One might wonder what happens when we attempt to infer more properties. Unfortunately, as the number of properties increases, the coarseness of the material meshes must also increase. This is simply due to matching displacement constraints coming from MR measurements to the property degrees of freedom in the equations of motion. Theoretical considerations aside, our numerical suite has a much harder time converging with the addition of an extra property to infer - a significant contraction of the domain of convergence being observed. We believe that this is a numerical issue, not a theoretical one, and an active research direction within our group.

As we move forward with in vitro and in vivo testing, our in silico results give us a good picture at how ignorant we can be when guessing our background values. Hydraulic permeability reconstruction for noiseless data is robust to initial conditions ranging from an order of magnitude below, to 2+ orders of magnitude above the ground truth background. Meanwhile, shear modulus alone has not been observed to diverge, albeit its convergence rate can get quite slow. Adding noise to data introduces fluctuations that will affect reconstruction, as our group has shown in prior work [114]. We found our results to be robust up to contrasts of an order of magnitude (it was small contrasts that were problematic) and believe that as long as signal to noise ratios are well above unity, noise should have only a minor effect on the domain of convergence.

4.1 Concluding Thoughts

Being able to image multiple poroelastic parameters will be a boon for preclinical diagnosis. Getting our computational suite to work in vivo will allow us to apply poroelasticity to the brain, where both cerebro-spinal fluid and blood play an active roll in homeostasis. The robustness under a wide range of initial conditions at low frequencies bodes well for using intrinsic actuation via the cardiac cycle. Tumors tend to have shear moduli between one and two orders of magnitude higher than surrounding tissue - these will be as visible to poroelastic parameter reconstruction as they are for viscoelasticity. The movement of fluid in the brain is critical in several venues. It is a mixture blood movement and interstitial fluid movement over the scale of a fraction of a voxel. We expect the dominant scale of fluid flow observed in MRE will be between that of perfusion which observes macroscopic bulk blood flow over a fraction of a meter and that of diffusion which observes microscopic fluid flow on a cellular scale. Observing fluid flow in this scale has potential to complement diffusion in characterizing edema and dementia. We postulate this scale of fluid flow will also illuminate tissue viability before diffusion becomes relevant allowing us to better identify and characterize reversible tissue damage. It might also provide additional information about functional changes secondary to brain activity.

It is an exciting time to watch as MRE enters it’s third decade of application. As interdisciplinary as its roots are, it is no surprise that with maturity the field is digging deeper into the fertile soil of its disciplinary constituents. Poroelasticity has been one of the crowning theoretical achievements of the geophysics community, and like every good idea, is ripe for cross-pollination. Perhaps it is no surprise that a model developed to describe rocks and soils is finding use in describing living tissue when in 1510 Leonardo da Vinci penned [115], into his Codex Leicester, that The Earth is a living body. Its flesh is the soil, its bones are the strata of rock, its cartilage is the tufa, its blood is the underground streams, the reservoir of blood around its heart is the ocean, the systole and diastole of the blood in the arteries and veins appear on the Earth as the rising and sinking of the oceans.

4.2 Methods and Data Availability

Simulations were performed using an in-house FORTRAN MRE Computational Platform written and developed over the past decade by group members. Both forward and inverse problems were computed on the Discovery Cluster, a 3000+ core Linux cluster at Dartmouth College administered by Research Computing. The data that support the findings of this study are available from the corresponding author upon reasonable request.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available under reasonable request by the corresponding author, without undue reservation.

Author Contributions

DS ran the simulations, did the data analysis, and wrote the manuscript with input and guidance from MM, SG, JW, and KP. MM and EH wrote the FEM code used for the simulations.

Funding

This work was supported by the National Institute for Health Research Grant (No. R01EB027577-02).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We would like to thank Research Computing at Dartmouth for their ever present ability with troubleshooting our cluster usage, and wish Susan Schwarz an adventurous retirement.

References

1. Dresner MA, Rose GH, Rossman PJ, Muthupillai R, Manduca A, Ehman RL. Magnetic resonance elastography of skeletal muscle. J Magn Reson Imaging (2001) 13:269–76. doi:10.1002/1522-2586(200102)13:2<269::AID-JMRI1039>3.0.CO;2-1

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Sack I, Bernarding J, Braun J. Analysis of wave patterns in mr elastography of skeletal muscle using coupled harmonic oscillator simulations. Magn Reson Imaging (2002) 20:95–104. doi:10.1016/S0730-725X(02)00474-5

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Bensamoun SF, Ringleb SI, Littrell L, Chen Q, Brennan M, Ehman RL, et al. Determination of thigh muscle stiffness using magnetic resonance elastography. J Magn Reson Imaging (2006) 23:242–7. doi:10.1002/jmri.20487

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Ringleb SI, Bensamoun SF, Chen Q, Manduca A, An KN, Ehman RL. Applications of magnetic resonance elastography to healthy and pathologic skeletal muscle. J Magn Reson Imaging (2007) 25:301–9. doi:10.1002/jmri.20817

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Chen Q, Manduca A, An K. Characterization of skeletal muscle elasticity using magnetic resonance elastography. Biomedical Applications of Vibration and Acoustics in Imaging and characterizations. New York: ASME Press (2008). doi:10.1115/1.802731.ch8

CrossRef Full Text | Google Scholar

6. Klatt D, Papazoglou S, Braun J, Sack I. Viscoelasticity-based MR elastography of skeletal muscle. Phys Med Biol (2010) 55:6445. doi:10.1088/0031-9155/55/21/007

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Aarabon N, Kozinc A, Podrekar N. Using shear-wave elastography in skeletal muscle: a repeatability and reproducibility study on biceps femoris muscle. PLoS One (2019) 14:1–13. doi:10.1371/journal.pone.0222008

CrossRef Full Text | Google Scholar

8. Yin L, Lu R, Cao W, Zhang L, Li W, Sun H, et al. Three-dimensional shear wave elastography of skeletal muscle: preliminary study. J Ultrasound Med (2018) 37:2053–62. doi:10.1002/jum.14559

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Nelissen JL, Sinkus R, Nicolay K, Nederveen AJ, Oomens CWJ, Strijkers GJ. Magnetic resonance elastography of skeletal muscle deep tissue injury. NMR Biomed (2019) 32:e4087. doi:10.1002/nbm.4087. E4087 NBM-18-0096.R2

PubMed Abstract | CrossRef Full Text | Google Scholar

10. McKnight AL, Kugel JL, Rossman PJ, Manduca A, Hartmann LC, Ehman RL. MR elastography of breast cancer: preliminary results. AJR Am J Roentgenol (2002) 178:1411–7. doi:10.2214/ajr.178.6.1781411

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Lorenzen J, Sinkus R, Lorenzen M, Dargatz M, Leussler C, Röschmann P, et al. MR elastography of the breast:preliminary clinical results. Rofo (2002) 174:830–4. doi:10.1055/s-2002-32690

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Van Houten EE, Doyley MM, Kennedy FE, Weaver JB, Paulsen KD. Initial in vivo experience with steady-state subzone-based MR elastography of the human breast. J Magn Reson Imaging (2003) 17:72–85. doi:10.1002/jmri.10232

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Xydeas T, Siegmann K, Sinkus R, Krainick-Strobel U, Miller S, Claussen CD. Magnetic resonance elastography of the breast: correlation of signal intensity data with viscoelastic properties. Invest Radiol (2005) 40:412–20. doi:10.1097/01.rli.0000166940.72971.4a

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Asbach P, Klatt D, Hamhaber U, Braun J, Somasundaram R, Hamm B, et al. Assessment of liver viscoelasticity using multifrequency mr elastography. Magn Reson Med (2008) 60:373–9. doi:10.1002/mrm.21636

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Hawley JR, Kalra P, Mo X, Raterman B, Yee LD, Kolipaka A. Quantification of breast stiffness using MR elastography at 3 Tesla with a soft sternal driver: a reproducibility study. J Magn Reson Imaging (2017) 45:1379–84. doi:10.1002/jmri.25511

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Bohte AE, Nelissen JL, Runge JH, Holub O, Lambert SA, de Graaf L, et al. Breast magnetic resonance elastography: a review of clinical work and future perspectives. NMR Biomed (2018) 31:e3932. doi:10.1002/nbm.3932

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Balleyguier C, Lakhdar A, Dunant A, Mathieu M, Delaloge S, Sinkus R. Value of whole breast magnetic resonance elastography added to mri for lesion characterization. NMR Biomed (2018) 31:e3795. doi:10.1002/nbm.3795

CrossRef Full Text | Google Scholar

18. Rouvière O, Yin M, Dresner MA, Rossman PJ, Burgart LJ, Fidler JL, et al. MR elastography of the liver: preliminary results. Radiology (2006) 240:440–8. doi:10.1148/radiol.2402050606

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Yin M, Woollard J, Wang X, Torres VE, Harris PC, Ward CJ, et al. Quantitative assessment of hepatic fibrosis in an animal model with magnetic resonance elastography. Magn Reson Med (2007) 58:346–53. doi:10.1002/mrm.21286

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Motosugi U, Ichikawa T, Sano K, Sou H, Muhi A, Koshiishi T, et al. Magnetic resonance elastography of the liver: preliminary results and estimation of inter-rater reliability. Jpn J Radiol (2010) 28:623–7. doi:10.1007/s11604-010-0478-1

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Venkatesh SK, Yin M, Ehman RL. Magnetic resonance elastography of liver: technique, analysis, and clinical applications. J Magn Reson Imaging (2013) 37:544–55. doi:10.1002/jmri.23731

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Chen J, Yin M, Glaser KJ, Talwalkar JA, Ehman RL. MR elastography of liver disease: state of the art. Appl Radiol (2013) 42:5–12.

PubMed Abstract | Google Scholar

23. Kennedy P, Wagner M, Castéra L, Hong CW, Johnson CL, Sirlin CB, et al. Quantitative elastography methods in liver disease: current evidence and future directions. Radiology (2018) 286:738–63. doi:10.1148/radiol.2018170601

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Hoodeshenas S, Yin M, Venkatesh SK. Magnetic resonance elastography of liver: current update. Top Magn Reson Imaging (2018) 27:319–33. doi:10.1097/RMR.0000000000000177

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Xu L, Lin Y, Xi ZN, Shen H, Gao PY. Magnetic resonance elastography of the human brain: a preliminary study. Acta Radiol (2007) 48:112–5. doi:10.1080/02841850601026401

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Kruse SA, Rose GH, Glaser KJ, Manduca A, Felmlee JP, Jack CR, et al. Magnetic resonance elastography of the brain. Neuroimage (2008) 39:231–7. doi:10.1016/j.neuroimage.2007.08.030

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Gierke HE, Oestreicher HL, Franke EK, Parrack HO, Wittern WW. Physics of vibrations in living tissues. J Appl Physiol (1952) 4:886–900. doi:10.1152/jappl.1952.4.12.886

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Manduca A, Oliphant TE, Dresner MA, Mahowald JL, Kruse SA, Amromin E, et al. Magnetic resonance elastography: non-invasive mapping of tissue elasticity. Med Image Anal (2001) 5:237–54. doi:10.1016/s1361-8415(00)00039-6

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Mariappan YK, Glaser KJ, Ehman RL. Magnetic resonance elastography: a review. Clin Anat (2010) 23:497–511. doi:10.1002/ca.21006

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Sarvazyan A, Hall TJ, Urban MW, Fatemi M, Aglyamov SR, Garra BS. AN overview OF elastography - an emerging branch OF medical imaging. Curr Med Imaging Rev (2011) 7:255–82. doi:10.2174/157340511798038684

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Tang A, Cloutier G, Szeverenyi NM, Sirlin CB. Ultrasound elastography and mr elastography for assessing liver fibrosis: Part 1, principles and techniques. AJR Am J Roentgenol (2015) 205:22–32. doi:10.2214/AJR.15.14552

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Low G, Kruse SA, Lomas DJ. General review of magnetic resonance elastography. World J Radiol (2016) 8:59–72. doi:10.4329/wjr.v8.i1.59

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Bayly PV, Garbow JR. Pre-clinical mr elastography: principles, techniques, and applications. J Magn Reson (2018) 291:73–83. doi:10.1016/j.jmr.2018.01.004

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Perriñez P, Marra S, Kennedy F, Paulsen K. “3D finite element solution to the dynamic poroelasticity problem for use in MR elastography,” in Medical imaging 2007: physiology, function, and structure from medical images. Editors A Manduca, and XP Hu, Washington: International Society for Optics and Photonics (SPIE) (2007), 6511, 402–12. doi:10.1117/12.709786

CrossRef Full Text | Google Scholar

35. Perriñez PR, Kennedy FE, Van Houten EE, Weaver JB, Paulsen KD. Magnetic resonance poroelastography: an algorithm for estimating the mechanical properties of fluid-saturated soft tissues. IEEE Trans Med Imaging (2010) 29:746–55. doi:10.1109/TMI.2009.2035309

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Pattison AJ, McGarry M, Weaver JB, Paulsen KD. Spatially-resolved hydraulic conductivity estimation via poroelastic magnetic resonance elastography. IEEE Trans Med Imaging (2014) 33:1373–80. doi:10.1109/TMI.2014.2311456

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Tan L, McGarry MD, Van Houten EE, Ji M, Solamen L, Weaver JB, et al. Gradient-based optimization for poroelastic and viscoelastic MR elastography. IEEE Trans Med Imaging (2017) 36:236–50. doi:10.1109/TMI.2016.2604568

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Darcy H. Les Fontaines publiques de la ville de Dijon. Exposition et application des principes à suivre et des formules à employer dans les questions de distribution d’eau, etc (Paris: V. Dalamont) (1856).

Google Scholar

39. Brown G. Henry Darcy and the making of a law. Water Resour Res (2002) 38 11.1–11.12. doi:10.1029/2001WR000727

CrossRef Full Text | Google Scholar

40. Delesse A. Pour déterminer la copasition des roches. Annales des Mines Paris (1849) 4:379–88.

Google Scholar

41. Coussy O. Poromechanics. New Jersey: John Wiley & Sons, LTD (2005).

Google Scholar

42. Lopatnikov S, Cheng A. Macroscopic Lagrangian formulation of poroelasticity with porosity dynamics. J Mech Phys Solid (2004) 52:2801–39. 10.1016/j.jmps.2004.05.005

CrossRef Full Text | Google Scholar

43. Bear J, Cheng A. Modeling groundwater flow and contaminant transport. Dordrecht: Springer, Vol. 23 (2010), 834. doi:10.1007/978-1-4020-6682-5

CrossRef Full Text | Google Scholar

44. Cheng A. Poroelasticity. Theory and applications of transport in porous media. New York: Springer International Publishing (2016).

Google Scholar

45. Tergazhi K. Der grundbruch an stauwerken and seine verhiltung. Die Wasserkraft 17 (1922) 445–9.

Google Scholar

46. Terzaghi K. Erdbaumechanik aufbodenphysikalischer grundlage. Leipzig: Franz Deuticke (1925).

Google Scholar

47. Terzaghi K. Theoretical soil mechanics. New Jersey: John Wiley & Sons Incorporated (1943). doi:10.1002/9780470172766

CrossRef Full Text | Google Scholar

48. Kadanoff L. Scaling laws for Ising models near Tc. Physics Physique Fizika (1966) 2:263.

CrossRef Full Text | Google Scholar

49. Wilson K. Renormalization group and critical phenomena. I. Renormalization group and the Kadanoff scaling picture. Physical review B (1971) 4:3174.

CrossRef Full Text | Google Scholar

50. Marrink SJ, Risselada HJ, Yefimov S, Tieleman DP, De Vries AH. The MARTINI force field: coarse grained model for biomolecular simulations. J Phys Chem B (2007) 111:7812–24. doi:10.1021/jp071097f

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Biot M. General theory of three-dimensional consolidation. J Appl Phys (1941) 12:155–64. doi:10.1063/1.1712886

CrossRef Full Text | Google Scholar

52. Detournay E, Cheng A. Fundamentals of poroelasticity. Analysis and design methods. Netherlands: Elsevier (1993), 113–71.

CrossRef Full Text | Google Scholar

53. Skempton A. The pore-pressure coefficients a and b. Geotechnique (1954) 4:143–7.

CrossRef Full Text | Google Scholar

54. Bishop A. The influence of an undrained change in stress on the pore pressure in porous media of low compressibility. Geotechnique (1973) 23:435–42.

CrossRef Full Text | Google Scholar

55. Plona T. Observation of a second bulk compressional wave in a porous medium at ultrasonic frequencies. Appl Phys Lett (1980) 36:259–61. doi:10.1063/1.91445

CrossRef Full Text | Google Scholar

56. Biot M. Theory of elasticity and consolidation for a porous anisotropic solid. J Appl Phys (1955) 26:182–5. doi:10.1063/1.1721956

CrossRef Full Text | Google Scholar

57. Biot M. Theory of propagation of elastic waves in a fluid-saturated porous solid. i. low-frequency range. J Acoust Soc Am (1956a) 28:168–78. doi:10.1121/1.1908239

CrossRef Full Text | Google Scholar

58. Biot M. Theory of propagation of elastic waves in a fluid-saturated porous solid. ii. higher frequency range. J Acoust Soc Am (1956b) 28:179–91. doi:10.1121/1.1908241

CrossRef Full Text | Google Scholar

59. Biot M. Nonlinear and semilinear rheology of porous solids. J Geophys Res (1973) 78:4924–37. doi:10.1029/JB078i023p04924

CrossRef Full Text | Google Scholar

60. Biot M. General solutions of the equations of elasticity and consolidation for a porous material. J Appl Mech (1956c) 78:91–6.

Google Scholar

61. Biot M, Willis D. The elastic coefficients of the theory of consolidation. J Appl Mech (1957) 24:594–601.

Google Scholar

62. Santos J. Introduction to the theory of poroelasticity. Indiana: Purdue Mathematics Department (2014).

Google Scholar

63. Nordbotten J, Celia M, Dahle H, Hassanizadeh S. Interpretation of macroscale variables in Darcy’s law. Water resources research (2007) 43:18. doi:10.1029/2006WR005018

CrossRef Full Text | Google Scholar

64. Apostolakis G, Dargush G. Mixed variational principles for dynamic response of thermoelastic and poroelastic continua. Int J Solid Struct (2013) 50:642–50. doi:10.1016/j.ijsolstr.2012.10.021

CrossRef Full Text | Google Scholar

65. Bear J, Bachmat Y. Introduction to modeling of transport phenomena in porous media. Berlin/Heidelberg: Springer Science & Business Media, Vol. 4 (1990).

Google Scholar

66. Berryman J. Confirmation of biot’s theory. Appl Phys Lett (1980) 37:382–4.

CrossRef Full Text | Google Scholar

67. Bonnet G, Auriault J. Dynamics of saturated and deformable porous media. Physics of finely divided matter. Berlin/Heidelberg: Springer (1985), 306–16.

CrossRef Full Text | Google Scholar

68. Perriñez PR, Kennedy FE, Van Houten EE, Weaver JB, aulsen KD. Modeling of soft poroelastic tissue in time-harmonic mr elastography. IEEE Trans Biomed Eng (2008) 56:598–608. doi:10.1109/TBME.2008.2009928

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Lynch D. Numerical partial differential equations for environmental scientists and engineers: a first practical course. Berlin/Heidelberg: Springer Science & Business Media (2004).

Google Scholar

70. McGarry MD, Van Houten EE, Johnson CL, Georgiadis JG, Sutton BP, Weaver JB, et al. Multiresolution mr elastography using nonlinear inversion. Med Phys (2012) 39:6388–96. doi:10.1118/1.4754649

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Franck IM, Koutsourelakis P. Sparse variational bayesian approximations for nonlinear inverse problems: applications in nonlinear elastography. Comput Methods Appl Mech Eng (2016) 299:215–44. doi:10.1016/j.cma.2015.10.015

CrossRef Full Text | Google Scholar

72. Jiang Y, Qian S. Bayesian approach for recovering piecewise constant viscoelasticity from mre data. Acta Mathematicae Applicatae Sinica, English Series (2020) 36:223–36. doi:10.1007/s10255-020-0922-7

CrossRef Full Text | Google Scholar

73. Cox R. The algebra of probable inference. Am J Phys (1963) 31:66–7.

CrossRef Full Text | Google Scholar

74. Jaynes E. Probability theory: the logic of science. Cambridge: Cambridge University Press (2003).

Google Scholar

75. Terenin A, Draper D. Cox’s theorem and the jaynesian interpretation of probability. arXiv preprint arXiv:1507.06597 (2015).

Google Scholar

76. Sowinski D. Complexity and stability for epistemic agents: the foundations and phenomenology of configurational entropy. New Hampshire: Dartmouth College (2016).

Google Scholar

77. Gleiser M, Sowinski D. How we make sense of the world: information, map-making, and the scientific narrative. The Map and the Territory. New York: Springer (2018), 141–63.

CrossRef Full Text | Google Scholar

78. Van Horn KS. Constructing a logic of plausible inference: a guide to cox’s theorem. Int J Approx Reason (2003) 34:3–24. doi:10.1016/S0888-613X(03)00051-3

CrossRef Full Text | Google Scholar

79. Caticha A. Lectures on probability, entropy, and statistical physics. arXiv preprint arXiv:0808.0012 (2008).

Google Scholar

80. Gelman A, Simpson D, Betancourt M. The prior can often only be understood in the context of the likelihood. Entropy (2017) 19:555. doi:10.3390/e19100555

CrossRef Full Text | Google Scholar

81. Caticha A, Preuss R. Maximum entropy and bayesian data analysis: entropic prior distributions. Phys Rev E Stat Nonlin Soft Matter Phys (2004) 70:046127. doi:10.1103/PhysRevE.70.046127

PubMed Abstract | CrossRef Full Text | Google Scholar

82. Shannon C. A mathematical theory of communication. Bell system technical journal (1948) 27:379–423.

CrossRef Full Text | Google Scholar

83. Muthupillai R, Lomas DJ, Rossman PJ, Greenleaf JF, Manduca A, Ehman RL. Magnetic resonance elastography by direct visualization of propagating acoustic strain waves. Science (1995) 269:1854–7. doi:10.1126/science.7569924

PubMed Abstract | CrossRef Full Text | Google Scholar

84. Sinkus R, Tanter M, Catheline S, Lorenzen J, Kuhl C, Sondermann E, et al. Imaging anisotropic and viscous properties of breast tissue by magnetic resonance-elastography. Magn Reson Med (2005) 53:372–87. doi:10.1002/mrm.20355

PubMed Abstract | CrossRef Full Text | Google Scholar

85. Sack I, Beierbach B, Hamhaber U, Klatt D, Braun J. Non-invasive measurement of brain viscoelasticity using magnetic resonance elastography. NMR Biomed (2008) 21:265–71. doi:10.1002/nbm.1189

PubMed Abstract | CrossRef Full Text | Google Scholar

86. Johnson CL, McGarry MD, Van Houten EE, Weaver JB, Paulsen KD, Sutton BP, et al. Magnetic resonance elastography of the brain using multishot spiral readouts with self-navigated motion correction. Magn Reson Med (2013a) 70:404–12. doi:10.1002/mrm.24473

PubMed Abstract | CrossRef Full Text | Google Scholar

87. Johnson CL, Holtrop JL, McGarry MD, Weaver JB, Paulsen KD, Georgiadis JG, et al. 3d multislab, multishot acquisition for fast, whole-brain mr elastography with high signal-to-noise efficiency. Magn Reson Med (2014) 71:477–85. doi:10.1002/mrm.25065

PubMed Abstract | CrossRef Full Text | Google Scholar

88. Wang H, Weaver JB, Doyley MM, Kennedy FE, Paulsen KD. Optimized motion estimation for mre data with reduced motion encodes. Phys Med Biol (2008) 53:2181. doi:10.1088/0031-9155/53/8/012

PubMed Abstract | CrossRef Full Text | Google Scholar

89. Wang H. Optimizing motion encoding and reconstruction in magnetic resonance elastography. New Hampshire: Dartmouth College (2008).

Google Scholar

90. Weaver JB, Pattison AJ, McGarry MD, Perreard IM, Swienckowski JG, Eskey CJ, et al. Brain mechanical property measurement using MRE with intrinsic activation. Phys Med Biol (2012) 57:7275–87. doi:10.1088/0031-9155/57/22/7275

PubMed Abstract | CrossRef Full Text | Google Scholar

91. McGarry M, Johnson CL, Sutton BP, Van Houten EE, Georgiadis JG, Weaver JB, et al. Including spatial information in nonlinear inversion MR elastography using soft prior regularization. IEEE Trans Med Imaging (2013) 32:1901–9. doi:10.1109/TMI.2013.2268978

PubMed Abstract | CrossRef Full Text | Google Scholar

92. Johnson CL, McGarry MD, Gharibans AA, Weaver JB, Paulsen KD, Wang H, et al. Local mechanical properties of white matter structures in the human brain. Neuroimage (2013b) 79:145–52. doi:10.1016/j.neuroimage.2013.04.089

PubMed Abstract | CrossRef Full Text | Google Scholar

93. Van Houten EE, Paulsen KD, Miga MI, Kennedy FE, Weaver JB. An overlapping subzone technique for mr-based elastic property reconstruction. Magn Reson Med 42 (1999) 779–86. doi:10.1002/(sici)1522-2594(199910)42:4<779::aid-mrm21>3.0.co;2-z

PubMed Abstract | CrossRef Full Text | Google Scholar

94. Van Houten EE, Miga MI, Weaver JB, Kennedy FE, Paulsen KD. Three-dimensional subzone-based reconstruction algorithm for mr elastography. Magn Reson Med (2001) 45:827–37. doi:10.1002/mrm.1111

PubMed Abstract | CrossRef Full Text | Google Scholar

95. McGarry M. Rayleigh damped magnetic resonance elastograpy. New Zealand: University of Canterbury. Mechanical Engineering (2008).

Google Scholar

96. Barber TW, Brockway JA, Higgins LS. The density of tissues in and about the head. Acta Neurol Scand (1970) 46:85–92. doi:10.1111/j.1600-0404.1970.tb05606.x

PubMed Abstract | CrossRef Full Text | Google Scholar

97. Cala LA, Thickbroom GW, Black JL, Collins DW, Mastaglia FL. Brain density and cerebrospinal fluid space size: ct of normal volunteers. AJNR Am J Neuroradiol (1981) 2:41–7.

PubMed Abstract | Google Scholar

98. DiResta GR, Lee J, Lau N, Ali F, Galicich JH, Arbit E. Measurement of brain tissue density using pycnometry. Acta Neurochir Suppl (1990) 51:34–6. doi:10.1007/978-3-7091-9115-6_12

PubMed Abstract | CrossRef Full Text | Google Scholar

99. Bloomfield IG, Johnston IH, Bilston LE. Effects of proteins, blood cells and glucose on the viscosity of cerebrospinal fluid. Pediatr Neurosurg (1998) 28:246–51. doi:10.1159/000028659

PubMed Abstract | CrossRef Full Text | Google Scholar

100. Yatsushiro S, Sunohara S, Atsumi H, Matsumae M, Kuroda K. Visualization and characterization of cerebrospinal fluid motion based on magnetic resonance imaging. Hydrocephalus: Water on the Brain (2018) 9:73302. doi:10.5772/intechopen.73302

CrossRef Full Text | Google Scholar

101. Venton J, Bouyagoub S, Harris P, Phillips G. “Deriving spinal cord permeability and porosity using diffusion-weighted MRI data,” in Poromechanics VI (Sixth Biot Conference on Poromechanics), Virginia: American Society of Civil Engineers, Paris, France, July 9‐13, 2017. Sponsored by École des Ponts ParisTech; Instutut Français des Sciences et Technologies des Transports, de l'Aménagement, et des Réseaux (IFSTTAR); Centre National de la Recherche Scientifique (CNRS); Agence Nationale de la Recherche (ANR); and the Engineering Mechanics Institute of ASCE. (2017), 1451–7. doi:10.1061/9780784480779.180

CrossRef Full Text | Google Scholar

102. Cheng S, Bilston LE. Unconfined compression of white matter. J Biomech (2007) 40:117–24. doi:10.1016/j.jbiomech.2005.11.004

PubMed Abstract | CrossRef Full Text | Google Scholar

103. Chen X, Sarntinoranont M. Biphasic finite element model of solute transport for direct infusion into nervous tissue. Ann Biomed Eng (2007) 35:2145–58. doi:10.1007/s10439-007-9371-1

PubMed Abstract | CrossRef Full Text | Google Scholar

104. Støverud K. Modeling convection-enhanced delivery into brain tissue using information from magnetic resonance imaging. Master’s thesis. Netherlands: Utrecht University (2009).

Google Scholar

105. Ray L, Heys J. Fluid flow and mass transport in brain tissue. Fluids (2019) 4:196. doi:10.3390/fluids4040196

CrossRef Full Text | Google Scholar

106. McGarry M, Houten EV, Solamen L, Gordon-Wylie S, Weaver J, Paulsen K. Uniqueness of poroelastic and viscoelastic nonlinear inversion MR elastography at low frequencies. Phys Med Biol (2019) 64:075006. doi:10.1088/1361-6560/ab0a7d

PubMed Abstract | CrossRef Full Text | Google Scholar

107. McGarry MD, Johnson CL, Sutton BP, Georgiadis JG, Van Houten EE, Pattison AJ, et al. Suitability of poroelastic and viscoelastic mechanical models for high and low frequency MR elastography. Med Phys (2015) 42:947–57. doi:10.1118/1.4905048

PubMed Abstract | CrossRef Full Text | Google Scholar

108. Schmid-Schönbein H, Wells RE. Rheological properties of human erythrocytes and their influence upon the “anomalous” viscosity of blood. Ergeb Physiol (1971) 63:146–219. doi:10.1007/BFb0047743

PubMed Abstract | CrossRef Full Text | Google Scholar

109. Gokturk HS, Demir M, Ozturk NA, Unler GK, Kulaksizoglu S, Kozanoglu I, et al. The role of ascitic fluid viscosity in the differential diagnosis of ascites. Can J Gastroenterol (2010) 24:255–9. doi:10.1155/2010/896786

PubMed Abstract | CrossRef Full Text | Google Scholar

110. Carman P. Fluid flow through granular beds. Chem Eng Res Des (1937) 15:S32–S48. doi:10.1016/S0263-8762(97)80003-2

CrossRef Full Text | Google Scholar

111. Kozeny J. Uber kapillare leitung der wasser in boden. Royal Academy of Science, Vienna, Proc. Class I (1927) 136:271–306.

Google Scholar

112. Budday S, Nay R, de Rooij R, Steinmann P, Wyrobek T, Ovaert TC, et al. Mechanical properties of gray and white matter brain tissue by indentation. J Mech Behav Biomed Mater (2015) 46:318–30. doi:10.1016/j.jmbbm.2015.02.024

PubMed Abstract | CrossRef Full Text | Google Scholar

113. Sinkus R, Siegmann K, Xydeas T, Tanter M, Claussen C, Fink M. Mr elastography of breast lesions: understanding the solid/liquid duality can improve the specificity of contrast-enhanced mr mammography. Magn Reson Med (2007) 58:1135–44. doi:10.1002/mrm.21404

PubMed Abstract | CrossRef Full Text | Google Scholar

114. Tan L, McGarry MDJ, Van Houten EEW, Ji M, Solamen L, Zeng W, et al. A numerical framework for interstitial fluid pressure imaging in poroelastic mre. Plos one (2017) 12:e0178521. doi:10.1371/journal.pone.0178521

PubMed Abstract | CrossRef Full Text | Google Scholar

115. da Vinci L, Gates B, Desmond M. Leonardo da Vinci : the Codex leicester: notebook of a genius. Sydney: Powerhouse Publishing (2000).

Google Scholar

Keywords: Elastography, MRI, Biophysics, Biomechanics, Bayesian Inference, Effective Field Theory, Continuum Mechanics, Porous Materials

Citation: Sowinski DR, McGarry MDJ, Van Houten EEW, Gordon-Wylie S, Weaver JB and Paulsen KD (2021) Poroelasticity as a Model of Soft Tissue Structure: Hydraulic Permeability Reconstruction for Magnetic Resonance Elastography in Silico. Front. Phys. 8:617582. doi: 10.3389/fphy.2020.617582

Received: 15 October 2020; Accepted: 04 December 2020;
Published: 21 January 2021.

Edited by:

Jing Guo, Charité – Universitätsmedizin Berlin, Germany

Reviewed by:

Sebastian Hirsch, Charité – Universitätsmedizin Berlin, Germany
Martina Guidetti, Rush University Medical Center, United States

Copyright © 2021 Sowinski, McGarry, Van Houten, Gordon-Wylie, Weaver and Paulsen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Damian R. Sowinski, RGFtaWFuLlNvd2luc2tpQERhcnRtb3V0aC5lZHU=

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.