Towards an acoustic model-based poroelastic imaging method: I. Theoretical Foundation

Ultrasound in Med. & Biol., Vol. 32, No. 4, pp. 547–567, 2006 Copyright © 2006 World Federation for Ultrasound in Medicine & Biology Printed in the USA. All rights reserved 0301-5629/06/$–see front matter

doi:10.1016/j.ultrasmedbio.2006.01.003

● Original Contribution TOWARDS AN ACOUSTIC MODEL-BASED POROELASTIC IMAGING METHOD: I. THEORETICAL FOUNDATION GEARÓID P. BERRY,* JEFFREY C. BAMBER,* CECIL G. ARMSTRONG,† NAOMI R. MILLER,* and PAUL E. BARBONE‡ *Institute of Cancer Research, Sutton, Surrey, UK; †School of Mechanical and Aerospace Engineering, QUB, Ashby Building, Belfast, UK; and ‡Aerospace and Mechanical Engineering, Boston University, Boston, MA, USA (Received 18 October 2004; revised 27 December 2005; in final form 5 January 2006)

Abstract—The ultrasonic measurement and imaging of tissue elasticity is currently under wide investigation and development as a clinical tool for the assessment of a broad range of diseases, but little account in this field has yet been taken of the fact that soft tissue is porous and contains mobile fluid. The ability to squeeze fluid out of tissue may have implications for conventional elasticity imaging, and may present opportunities for new investigative tools. When a homogeneous, isotropic, fluid-saturated poroelastic material with a linearly elastic solid phase and incompressible solid and fluid constituents is subjected to stress, the behaviour of the induced internal strain field is influenced by three material constants: the Young’s modulus (Es) and Poisson’s ratio (␯s) of the solid matrix and the permeability (k) of the solid matrix to the pore fluid. New analytical expressions were derived and used to model the time-dependent behaviour of the strain field inside simulated homogeneous cylindrical samples of such a poroelastic material undergoing sustained unconfined compression. A model-based reconstruction technique was developed to produce images of parameters related to the poroelastic material constants (Es, ␯s, k) from a comparison of the measured and predicted time-dependent spatially varying radial strain. Tests of the method using simulated noisy strain data showed that it is capable of producing three unique parametric images: an image of the Poisson’s ratio of the solid matrix, an image of the axial strain (which was not time-dependent subsequent to the application of the compression) and an image representing the product of the aggregate modulus Es(1⫺␯s) ⁄ (1⫹␯s)(1⫺2␯s) of the solid matrix and the permeability of the solid matrix to the pore fluid. The analytical expressions were further used to numerically validate a finite element model and to clarify previous work on poroelastography. (E-mail: [email protected]) © 2006 World Federation for Ultrasound in Medicine & Biology. Key Words: Poroelastography, Elastography, Ultrasound, Poroelastic, Biphasic, Strain, Stress relaxation, Analytical, Finite element, Porosity, Permeability, Young’s modulus, Poisson’s ratio, Oedema, Cartilage, Cancer.

INTRODUCTION

To measure the mechanical properties of tissue ultrasonically, work toward analysing rather than displaying ultrasonic echoes began some 25 years ago (see for example, Dickinson and Hill 1980). The development of acoustic elastography (Ophir et al. 1991) in which the echo analysis was aimed at producing an image of tissue strain represented a further important milestone in this subject. By scanning the tissue ultrasonically before and after an applied uniaxial compression, and performing cross-correlation-based RF echo-tracking, Ophir et al. (1991) imaged the induced strain along the direction of beam propagation. Konofagou et al. (1998) later extended the technique to enable strain imaging in a direction perpendicular to the acoustic axis, but in the scan plane. Treating the tissue as a linearly elastic, isotropic material allowed these two imaged components of the

Changes in the mechanical properties of soft biological tissue have long been associated with the onset and progression of disease (Wellman et al. 1999, 2001). This is the basis for the routine clinical assessment of tissue using manual palpation. By pressing gently with the hand or fingers, the physician may subjectively assess the local tissue stiffness and mobility, and relate these to the pathological state of the tissue (Dixon and Mensel 1994; Harrison 1991; Talley and O’Connor 1988).

Address correspondence to: Gearóid P. Berry, Joint Department of Physics, Institute of Cancer Research and Royal Marsden NHS Foundation Trust, Downs Road, Sutton, Surrey SM2 5PT, UK. E-mail: [email protected] 547

548

Ultrasound in Medicine and Biology

strain field, as a first approximation, to be related to the two underlying material constants that fully characterise the deformation of such a material: the Young’s modulus (E) and the Poisson’s ratio (␯). However, in general, a purely elastic model does not accurately describe the mechanical behaviour of biological tissues during compression. Many tissues are known to exhibit time-dependent compression-induced stress and strain (Duck 1990; Fung 1993; Lakes 1998). One possible mechanism for such transient behaviour in materials is the compression-induced flow of unbound mobile fluid through the material (Armstrong et al. 1984; Lakes 1999; Mow et al. 1980a). If this mechanism is the source of all or part of the transient stress-strain behaviour, the material is said to be at least in part poroelastic. Thus, factors such as the vascularity (Dickinson and Hill 1980), the volume of mobile interstitial fluid (Levick 2003) and the ease with which fluid can flow through various fluid compartments such as the interstitium (Levick 2003), may have an important influence on the response of the tissue to an applied compression. Therefore, it may be possible to extract information about these quantities by studying the response of the tissue to the compression. At the very least, if simple elastograms (strain images) are to be used to visually interpret variations in tissue stiffness, or if inverse methods (e.g., Barbone and Bamber 2002; Doyley et al. 2000; Kallel and Betrand 1996; Oberai et al. 2004) are to be successful at reconstructing quantitative images of parameters such as the Young’s modulus, it will be important to understand the implications of neglecting poroelastic behaviour. Both the volume of free interstitial fluid and the ease with which this fluid can move through the interstitium can be important diagnostic indicators. For example, the severity of lymphoedema is often graded by the girth of the affected limb (Moseley et al. 2002; Stanton et al. 2000), which, in turn, is related to the volume of interstitial fluid present. Furthermore, oedema is classified as “pitting” or as “brawny” based on how readily interstitial fluid may be displaced following sustained compression with the thumb or calipers (Roberts et al. 1995; Stanton et al. 2000). The behaviour under compression of materials rich in mobile fluid has previously been mathematically investigated. Poroelastic theory (Biot 1941, 1962; Wang 2000) treats such materials as being composed of two phases: a sponge-like, elastic, porous, permeable, solid phase called the matrix, and a liquid phase that inhabits the pores of the solid matrix. A special case of the Biot theory, developed independently in the biomechanical context, is known as the KLM biphasic theory (Armstrong et al. 1984; Mow et al. 1980a, 1980b, 1984, 1990). The KLM biphasic theory, and more recent generalisa-

Volume 32, Number 4, 2006

tions thereof, considers a poroelastic material to be a binary mixture of a solid and a fluid phase. The KLM biphasic theory has already been applied to the deformation of cartilage tissue (Armstrong et al. 1984; Mow et al. 1980a, 1980b, 1984, 1990) and so is more readily employed for the purposes of this paper. The theory assumes that a poroelastic material consists of a linearly elastic, isotropic, incompressible solid phase saturated with an incompressible fluid phase. It has been shown that a homogeneous sample of such a material is characterised by the following material constants: the Young’s modulus (Es) of the solid matrix, the Poisson’s ratio (␯s) of the solid matrix and the permeability (k) of the solid matrix to the pore fluid (Armstrong et al. 1984). It is important to note that a sample of such a material consisting of incompressible solid and fluid constituents may itself be compressible if the compression causes fluid to leave the sample. Es and ␯s are equivalent to the values of the Young’s modulus and the Poisson’s ratio of the sample under drained conditions i.e., when all compression-induced fluid flow has ceased. The permeability, k, is defined by Darcy’s Law, which states that the flux of fluid through a particular region of a material is proportional to the pressure gradient across the region and is given as follows (Biot 1941): q ⫽ ⫺k ⵜ P

(1)

where q is the fluid volume flowing through a region of unit cross-sectional area per unit time and ⵜP is the pressure gradient across the region. In Cartesian (x,y,z) coordinates, the operator ⵜ is given by: ⵜ⫽i

⭸ ⭸ ⭸ ⫹ j ⫹k ⭸x ⭸y ⭸z

(2)

where i, j and k are the corresponding unit vectors along each of the coordinate axes. The units of q are m s⫺1, those of ⵜP are Pa m⫺1 and, thus, the units of permeability are m4 N⫺1 s⫺1. Alternatively, Darcy’s Law for the flow of fluid through a vertical column of porous material may be written as follows (Bear 1988; Bowles 1992): q ⫽ Ki

(3)

where K is the hydraulic conductivity and i is the hydraulic gradient (equal to the difference in hydraulic head across the column divided by the length of the column); i is dimensionless and thus K is measured in units of m s⫺1. The permeability and the hydraulic conductivity are related by the following expression: K ⫽ k␥w

(4)

where ␥w is the specific weight of the pore fluid measured in units of N m⫺3. ␥w is related to the acceleration

Acoustic model-based poroelastic imaging ● G. P. BERRY et al.

549

due to gravity (g) and the density (␳) of the pore fluid as follows:

␥w ⫽ ␳g

(5)

The units of g are m s⫺2 and those of ␳ are kg m⫺3. Unfortunately, the terms “permeability” and “hydraulic conductivity” are frequently used loosely or interchangeably in the soil mechanics literature. To further confuse issues, these same terms are sometimes used to refer to similar, but not identical, quantities in the physiology literature. It has previously been suggested that the methods of elastography could be extended to image the time-dependent strain field induced during the compression of tissues rich in mobile fluid (Konofagou et al. 1998). The extended technique has been called poroelastography (Konofagou at al. 2001), and preliminary investigations have previously been carried out to examine its potential as a technique for imaging time-dependent strain in cylindrical samples. It is convenient to introduce a cylindrical polar coordinate system that is independent of ultrasonic transducer orientation where the axial, radial and circumferential directions are identified in Fig. 1a. By convention, expansive and compressive strains are treated as positive and negative quantities respectively. For materials with 0 ⬍ ␯ ⱕ 0.5, an axial compression causes both radial and circumferential expansion and thus both the ratio of the radial strain to the axial strain (radial-to-axial strain ratio) and the ratio of the circumferential strain to the axial strain (circumferential-to-axial strain ratio) are negative quantities. When ␯ ⫽ 0, both strain ratios are zero. To facilitate a comparison between this and previous work, when we refer to either of these strain ratios in this paper, we are in fact referring to the quantity obtained when the strain ratio is multiplied by ⫺1, which for materials with 0 ⱕ ␯ ⱕ 0.5 is equivalent to taking the absolute value of the strain ratio. According to this coordinate system and sign convention, Konofagou et al. (2001) used finite element modeling (FEM) to predict the behaviour of the radialto-axial strain ratio induced inside unconfined homogeneous cylindrical samples of a poroelastic material undergoing stress relaxation. The term “stress relaxation” is used to describe the compression of a viscoelastic or poroelastic material by a sustained, constant, applied strain and the term is derived from the reducing stress required to sustain the applied strain when such a material is being compressed (Armstrong et al. 1984; Fung 1993; Lakes 1998; Mow et al. 1980a). A time-dependent, spatially varying radial-to-axial strain ratio that eventually reached equilibrium was predicted by Konofagou’s simulation. By simulating the ultrasonic scanning of the

Fig. 1. (a) A schematic of a cylindrical sample of a poroelastic material of radius a. The axial and radial directions are along the z-axis and r-axis, respectively. The circumferential direction is in the direction of increasing angle ␪. The axial and radial directions used here are, respectively, the axial and lateral directions referred to in Konofagou et al. (2001) and Righetti et al. (2004), and the radial direction used here is equivalent to the radial direction described in Fortin et al. (2003). The objective of this analysis is to determine the compression-induced, time-dependent strain field inside the cylinder. Because of the axisymmetric nature of the problem, the solution throughout the entire cylinder is obtained by determining the strain field on the solution plane, and then rotating the solution plane about the axis of symmetry of the cylinder. (b) The sample is compressed between two compressor plates. The compression is applied in the ⫺z direction. The initial elastic response causes expansion of the sample in the direction of increasing r. Subsequent fluid flow in the direction of increasing r leads to radial contraction of the sample toward the axis of symmetry.

material with a transducer pointed in the direction along which the compression was applied (i.e., in the axial direction), it was demonstrated that, if the echo signalto-noise ratio were high enough, elastographic techniques might be capable of accurately imaging the temporal variation of the radial-to-axial strain ratio. However, the accuracy of the finite element model used by Konofagou et al. (2001) was not validated by comparing its predictions, for a given radial position inside the cylindrical sample, either with those of a corresponding analytical model or with experimental results. Furthermore, no specific proposal was given for using measurements of the spatio-temporal variation of the strain ratio to reconstruct images of the poroelastic material constants described above. Few attempts to measure or image experimentally the internal deformation (i.e., strain) of poroelastic ma-

550

Ultrasound in Medicine and Biology

terials using elastographic techniques have been reported. Fortin et al. (2003) used A-mode elastographic techniques to measure the time-varying radial displacement profile across the diameter of a cylindrical sample of cartilage undergoing stress relaxation. For each time instant, the relative radial displacement at the curved cylindrical boundary of the sample was then used to calculate the time-evolution of the average radial-toaxial strain ratio for the sample as a whole. However, neither local radial strain nor local radial-to-axial strain ratio profiles inside the sample were obtained, nor were they considered in this paper. Righetti et al. (2003, 2004), on the other hand, have recently used elastographic methods to produce time-sequenced images of the radial-to-axial strain ratio inside poroelastic cylindrical samples undergoing stress relaxation. Cylinders of tofu were used for their experiments (tofu is poroelastic and relatively homogeneous samples can be easily obtained). It was concluded that the feasibility of generating time-sequences of strain images to map and quantify the time-dependent mechanical behaviour of porous media had been demonstrated (Righetti et al. 2004). However, although not noted by the authors, on visual inspection of the published images and graphs, it may be seen that these experimental images do not agree fully with those predicted by the simulations of Konofagou et al. (2001), particularly in relation to the spatial distribution of radial-to-axial strain ratio values as a function of time. The specific nature of this discrepancy is discussed below. Previously, only limited attention has been given to the interpretation of the average radial-to-axial strain ratio and the spatial distribution of the local radial-toaxial strain ratio within the sample, in terms of the poroelastic material constants of the sample. According to the KLM biphasic theory, during unconfined compression, the average radial-to-axial strain ratio at equilibrium is equal to ␯s (Armstrong et al. 1984). Konofagou et al. (2001) further showed by simulation that, at equilibrium, the local radial-to-axial strain ratio is uniform throughout such a homogeneous poroelastic material, and is therefore also equal to ␯s. Righetti et al. (2004) reasoned that, at equilibrium, images of the radial-toaxial strain ratio could be interpreted as images of ␯s, and any inhomogeneities in the strain ratio image at equilibrium could indicate local variability in ␯s. Konofagou et al. (2001) suggested that the permeability, k, could be determined from the time constant of the variation in the radial-to-axial strain ratio, but did not discuss how this could be achieved. Previous work, therefore, was extremely important in defining the new field of poroelastography, but left substantial gaps in understanding the nature of the spatio-temporal variation of the strain field within the sam-

Volume 32, Number 4, 2006

ple, and how this behaviour could be interpreted in terms of the material constants of the sample. The purpose of this paper is to increase the level of this understanding and, through so doing, take a first step toward imaging parameters related to each of the material constants (Es, ␯s and k) of a homogeneous cylindrical sample of a poroelastic material. By applying the KLM biphasic theory, also invoked in previous work on poroelastography (Konofagou et al. 2001; Righetti et al. 2004), new analytical expressions are derived that describe, for the first time, how the internal strain field during stress-relaxation of such a material is influenced by the three material constants. The predictions made by these expressions permit the behaviour of the internal strain field to be investigated, and enable previous work by others to be clarified. A model-based parametric imaging technique is then described that works by fitting elastographically measured values for the radial strain as a function of time and radial position to the analytical expression describing the time-evolution of the spatially dependent radial strain. This allows reconstruction of quantities related to the underlying poroelastic material constants. The behaviour of this imaging method, particularly with regard to the uniqueness, accuracy and precision of the solutions that it provides, is studied in simulation, with and without strain measurement noise. Since the analytical expressions are valid only when specific accompanying boundary conditions are satisfied, in preparation for future work that may involve a variety of boundary conditions, one of the analytical expressions is also used to numerically validate the predictions of a finite element model. The mutual agreement obtained lends further confidence to both predictions. Furthermore, experimental evidence of strain field behaviour consistent with predictions made herein has been reported orally, together with an investigation of the experimental feasibility of the poroelastic image reconstruction technique (Berry et al. 2004a; Berry at al. 2004b). A detailed description of this experimental work is provided in a companion paper. METHODS Two methods were employed to solve the biphasic equations that predict the time-evolution of the strain field within a poroelastic cylinder undergoing stress relaxation, viz., calculation from analytical expressions and finite element modeling. In the results section below, the results from these two methods are compared to demonstrate their level of agreement. In this section, we begin by deriving the new analytical expressions required for the first method, the starting point for which is a restatement of an expression obtained by Armstrong et al. (1984). Following this is a description of the finite ele-

Acoustic model-based poroelastic imaging ● G. P. BERRY et al.

ment simulation tool and the procedures for employing it. Finally, we describe the method used for multi-parameter fitting of the theoretical model’s predictions to simulated experimental data, which it is hoped will enable the estimation of quantities related to unknown poroelastic material properties from experimental time-dependent strain field measurements.

Analytical theory Armstrong et al. (1984) considered the problem of the unconfined compression of a cylindrical sample of a homogeneous poroelastic material. The compression was modeled by placing the sample between two rigid, parallel, impermeable plates, Fig. 1b. Perfect slip (i.e., frictionless) boundary conditions were assumed to exist between the plates and the sample. One of the plates was then instantaneously and infinitesimally displaced in the axial direction, subjecting the sample to an externally applied compressive axial strain, ␧zz, which was held constant thereafter. The axial strain history was thus given by: ␧zz共t兲 ⫽

再

t⬍0

0, ⫺␧0,

tⱖ0

⌬V ⫽ ␧ ⫽ ␧zz ⫹ ␧rr ⫹ ␧␪␪ V

(7)

where ␧ is the volumetric strain (also known as the dilatation) and ␧zz, ␧rr and ␧␪␪ are the normal strains in the axial, radial and circumferential directions, respectively. For stress-relaxation, ␧zz is given by eqn 6 and, in the absence of shear strain, the radial and circumferential strains may be written as follows: ␧rr ⫽

⭸u ⭸r

(8)

u r

(9)

␧␪␪ ⫽

erning equation for the radial displacement, u, of the solid phase:

where u is the radial displacement at radial position r. (The reader is referred to Brady and Brown (1985) for a rigorous mathematical derivation of eqn 9 above, or to Appendix A for a geometrical validation). Armstrong et al. (1984) obtained the following gov-

冉

⭸2u 1 ⭸u u 1 ⫹ ⫺ ⫽ ⭸r2 r ⭸r r2 HAk

⭸u r ⭸␧zz ⫹ ⭸t 2 ⭸t

冊

(10)

where HA is the aggregate modulus of the solid matrix and is defined in terms of the other material constants as follows: Es共1 ⫺ ␯s兲 共1 ⫹ ␯s兲共1 ⫺ 2␯s兲

HA ⫽

(11)

The total stress in the sample is defined to be equal to the sum of the stresses in the solid matrix and in the fluid (Armstrong et al. 1984). By assuming that both the radial component of the total tissue stress and the fluid pressure were equal to zero at r ⫽ a (where a is the sample radius), the following solution to eqn 10 was obtained (see eqn 31 in Armstrong et al. 1984):

冤

u a, t ⫽ ␧0 ␯s ⫹ 共1 ⫺ 2␯s兲共1 ⫺ ␯s兲 a共 兲

(6)

where t represents time and ␧0 is the magnitude of the applied axial strain. When a cylindrical sample is compressed, for infin⌬V itesimal strains, the fractional change in volume, , of V a cylindrical volume element is given as follows (Brady and Brown 1985; Jaeger and Cook 1976; Wang 2000):

551

⬁

⫻

兺

n⫽1

冉

冊

⫺␣n2HAkt a2 2 2 兵␣n共1 ⫺ ␯s兲 ⫺ 共1 ⫺ 2␯s兲其 exp

冥

(12)

where u(a,t) is the radial displacement at r ⫽ a as a function of time t and ␣n are the roots of the following characteristic equation: J1共x兲 ⫺

共1 ⫺ ␯s兲 共1 ⫺ 2␯s兲

xJ0共x兲 ⫽ 0

(13)

where J0(x) is a Bessel function of the first kind of order zero and J1(x) is a Bessel function of the first kind of order one. The quantity u/a describes the average radial strain along the entire radius of the cylindrical sample and, thus, eqn 12 describes the time-evolution of the average radial strain in the sample. Referring to eqn 9, in the absence of shear strains, u/a also represents the circumferential or “hoop” strain at the curved cylindrical surface; it does not describe the value of the local radial strain on the curved cylindrical surface. From eqn 12, the average radial-to-axial strain ratio in the sample as a function of time is given as follows: u a, t ⫽ ␯s ⫹ 共1 ⫺ 2␯s兲共1 ⫺ ␯s兲 a␧0 共 兲 ⬁

⫻

兺

n⫽1

冉

冊

⫺␣n2HAkt a2 兵␣n2共1 ⫺ ␯s兲2 ⫺ 共1 ⫺ 2␯s兲其 exp

(14)

552

Ultrasound in Medicine and Biology

Volume 32, Number 4, 2006

冤

Thus, as t¡⬁, the average radial-to-axial strain ratio in the sample tends toward the value of the Poisson’s ratio of the solid matrix with a characteristic time governed by the Poisson’s ratio and aggregate modulus of the solid matrix, the permeability of the solid matrix to the fluid and the sample radius. Note that for non-cylindrical samples, the characteristic time would be influenced by the characteristic length of the sample and the sample geometry. This completes the restatement of the expressions provided or evaluated by Armstrong et al. (1984). We now derive new expressions, for the time dependence of the internal strain field inside the sample. It is trivial to show from eqns 12 or 14 that the radial displacement of the curved cylindrical surface as a function of time is given as follows:

冤

u共r, t兲 ⫽ r␧0

J0共␣n兲 ⫺ ⬁

␯s ⫹

兺

冉

1 ⫺ 2␯s 1 ⫺ ␯s

u共a, t兲 ⫽ a␧0 ␯s ⫹ 共1 ⫺ 2␯s兲共1 ⫺ ␯s兲

⬁

⫻

冊

冥

(15)

To evaluate the spatial dependence of the time-evolution of the strain field, it is first necessary to generalise eqn 15 to obtain an expression that describes the radial displacement as a function both of radial position, r, and time, t. The method that we have used to derive such an expression employs the Heaviside expansion theorem, and was described in detail by Armstrong et al. (1984). When this method is applied, the following new expression is obtained:

冉 冊 冊冦 J1

␣nr J1共␣n兲 a ⫹ ␣nr ␣n a

1 J ␣ ⫺␣ J ␣ 共1 ⫺ ␯s兲 0共 n兲 n 1共 n兲

n⫽1

兺

n⫽1

冉

⫺␣n2HAkt a2 2 2 兵␣n共1 ⫺ ␯s兲 ⫺ 共1 ⫺ 2␯s兲其 exp

冧

冉

exp

⫺␣n2HAkt a2

冊

冥

(16)

for 0 ⬍ r/a ⱕ 1. Note that it can be shown that this new expression reduces to the equation published by Armstrong et al. (1984), viz., eqn 12, at r ⫽ a (this is demonstrated in Appendix B). Differentiating this expression with respect to r and using eqn 8 delivers the following expression for the radial strain as a function of r and t:

⭸u ⫽ ␧0 ␧rr共r, t兲 ⫽ ⭸r

冤

⬁

␯s ⫹

冉

1 ⫺ 2␯s J0共␣n兲 ⫺ 1 ⫺ ␯s

兺

冉 冊 冊冦冉 冊 ␣nr J0 ⫺ a

J1

␣nr a J1共␣n兲 ⫹ ␣nr ␣n a

1 J ␣ ⫺␣ J ␣ 1 ⫺ 共 ␯s兲 0共 n兲 n 1共 n兲

n⫽1

冧

冉

exp

⫺␣n2HAkt a2

冊

冥

(17)

for 0 ⬍ r/a ⱕ 1. From eqns 9 and 16, the circumferential strain as a function of r and t is given as follows:

冤

u ␧␪␪共r, t兲 ⫽ ⫽ ␧0 r

J0共␣n兲 ⫺ ⬁

␯s ⫹

兺

n⫽1

冉

1 ⫺ 2␯S 1 ⫺ ␯S

冊冦 冉 冊 J1

␣nr J1共␣n兲 a ⫹ ␣nr ␣n a

1 J ␣ ⫺␣ J ␣ 共1 ⫺ ␯s兲 0共 n兲 n 1共 n兲

冧

冉

exp

⫺␣n2HAkt a2

冊

冥

(18)

By dividing eqns 17 and 18 by ␧0, the following expressions for the radial-to-axial and the circumferential-to-axial strain ratios as a function of r and t are easily obtained:

Acoustic model-based poroelastic imaging ● G. P. BERRY et al.

J0共␣n兲 ⫺ ␧rr r, t ⫽ ␯s ⫹ ␧0 共 兲

⬁

冉

兺

J0

␣nr ⫺ a

␣nr a J1共␣n兲 ⫹ ␣nr ␣n a

1 J ␣ ⫺␣ J ␣ 共1 ⫺ ␯s兲 0共 n兲 n 1共 n兲

n⫽1

␧␪␪ r, t ⫽ ␯s ⫹ ␧0 共 兲

1 ⫺ 2␯s 1 ⫺ ␯s

冊冦冉 冊 冉 冊 J1

⬁

冉

1 ⫺ 2␯s J0共␣n兲 ⫺ 1 ⫺ ␯s

兺

冊 冉 冊冦 J1

␣nr a J1共␣n兲 ⫹ ␣nr ␣n a

1 J ␣ ⫺␣ J ␣ 共1 ⫺ ␯s兲 0共 n兲 n 1共 n兲

n⫽1

for 0 ⬍ r/a ⱕ 1. Eqns 16 to 20 are only given for the half plane that passes through the axis of symmetry of the cylinder. However, referring to Fig.1a, because of the axisymmetric nature of the problem, the solution throughout the entire cylinder may be obtained by rotating this planar solution about the z axis. Results are presented below for the strain ratios, evaluated from eqns 19 and 20 and plotted as a function of r and t, for three hypothetical cylindrical poroelastic samples with the material constants described in Table 1. The Es and k values were fixed across the three samples, and chosen to match published values for tofu, an experimentally useful poroelastic phantom material (Righetti et al. 2004). The Poisson’s ratio of the solid matrix was varied between the samples over a wide range. The radius, a, of each sample was given by a ⫽ 25 mm and the magnitude of the applied axial strain was given by ␧0 ⫽ 1%. Finite element simulation A commercial finite-element modeling tool (ABAQUS, Abaqus Inc., Providence, RI, USA) was also used to predict the spatial and temporal behaviour of the radial-to-axial strain ratio during the sustained compression of each of the poroelastic samples described in Table 1, under the same boundary conditions as those modeled analytically. ABAQUS uses an “effective stress” princi-

Table 1. Values used for the material constants of the simulated samples of poroelastic materials in Figures 2 to Figure 12. HA may be determined from Es and ␯s Sample name

Es (kPa)

k (m4 N⫺1 s⫺1)

␯s

HA (kPa)

A B C

1.29 1.29 1.29

1.373 ⫻ 10⫺11 1.373 ⫻ 10⫺11 1.373 ⫻ 10⫺11

0 0.2 0.3

1.29 1.43 1.74

冧

553

冧 冉

冉

exp

⫺␣n2HAkt a2

⫺␣n2HAkt exp a2

冊

冊

(19)

(20)

ple to describe the mechanical behaviour of poroelastic elements (Hibbitt et al. 2003). The homogeneous samples were modeled as consisting of a linearly elastic, isotropic, incompressible, permeable solid phase saturated with an incompressible fluid phase. Squeezing the sample between two plates was simulated by displacing the upper surface of the cylinder a distance equivalent to 1% of the sample height in the ⫺z direction while preventing the lower surface of the cylinder from moving in the z direction. Once applied, this axial strain was held constant thereafter. Perfect slip (i.e., frictionless) boundary conditions were enforced on the upper and lower surfaces. The permeability throughout the sample was assumed to be strain-independent. The finite element mesh used to predict the behaviour on the solution plane in Fig. 1a was composed of 400 axisymmetric elements (CAX4RP) arranged in a single row. The left hand side of the model corresponded to the axis of symmetry of the cylinder and the right hand side corresponded to the curved cylindrical surface. The model dimensions were 25 mm (radial) ⫻ 50 mm (axial). Fluid flow was prevented through the upper and lower edges and the left hand side of the model, but was enabled in the positive radial direction by setting the pore pressure equal to zero at the right hand side of the model. The time-step was automatically adjusted by the software during a given sustained compression and was allowed to vary between a minimum of 0.1 s and a maximum of 100 s. As far as can be determined from the information provided in their paper, the only certain differences between this finite element model and that constructed by Konofagou et al. (2001), who also employed ABAQUS, are the values of the material constants of the sample and the cylinder dimensions. Although, we note, that the mesh density was not specified in Konofagou et al. (2001), and the values of the time-steps used were not explicitly stated. The units of permeability used in eqn 19 (m4 N⫺1

554

Ultrasound in Medicine and Biology

Volume 32, Number 4, 2006

Table 2. Fit parameter values used to generate the radial strain data for Sample C in order to test the fitting routine (true values), and the mean pixel values in the corresponding parametric images obtained from the noiseless and noisy radial strain data (imaged values). Error estimates represent one standard deviation Fit parameter

True value

Imaged value (noiseless)

Imaged value (5% noise)

HAk (m2 s⫺1) ␯s ␧0 (%)

2.38902 ⫻ 10⫺8 0.30000 1.00000

2.38912 (⫾0.00066) ⫻ 10⫺8 0.30000 (⫾0.00002) 1.00007 (⫾0.00008)

2.39810 (⫾0.19464) ⫻ 10⫺8 0.30071 (⫾0.00655) 0.99777 (⫾0.02185)

s⫺1) are not the same as those used for permeability in ABAQUS (m s⫺1). This is because hydraulic conductivity is referred to as permeability in ABAQUS, and this alternative use of terminology is noted in the software documentation (Hibbitt et al. 2003). Referring to eqn 4, hydraulic conductivity and permeability are related by the specific weight of the fluid and are numerically equal when the specific weight of the fluid is equal to 1 N m-3. To enable a valid comparison with the predictions of eqn 19, the specific weight of the fluid in the finite element model was set equal to 1 N m⫺3. Curve-fitting algorithm A program was written in Matlab (The MathWorks, Inc., Natick, MA, USA), using an implementation of the Levenberg-Marquardt algorithm, to find quantities related to the poroelastic material constants by obtaining a best fit of eqn 17 to time-dependent radial strain data. The sum of the squared differences was used as the fitting cost function to be minimised by iteration. Eqn 17 contains a sum of exponential terms. The difficulties associated with fitting such an equation to data, particularly in relation to the uniqueness with which the constituent exponents may be separated, are well documented (Atkins 1969; Lanczos 1957; Worseley 1964). In this case, however, the multi-exponential fit is greatly simplified by the fact that the exponents in eqn 17 are related through the roots of the characteristic equation (eqn 13). A three-parameter fitting routine was used to fit to the radial strain data with the three fit parameters p1, p2 and p3 given by: p1 ⫽ HAk p2 ⫽ ␯S

(21)

p3 ⫽ ␧0 As discussed later, this choice of parameters was made to avoid non-uniqueness in the set of parameters that delivered the best fit. While the sample radius influenced the behaviour of the strain field, the radius was treated as a known quantity a priori and was therefore not a fit parameter.

To test the curve-fitting routine, a time sequence of images of the radial strain as a function of time during the stress-relaxation of Sample C (see Table 1) was obtained using the analytical model, i.e., eqn 17. It was first verified that the radial strain data produced by the parameter set [p1, p2, p3] and radius of Sample C were unique to that parameter set and radius combination. This was achieved by sampling parameter space at regular intervals to obtain a large number of parameter combinations other than those belonging to Sample C. These parameter combinations were then used to generate more radial strain data to simulate the behaviour of many different samples with different parameter values but with the same radius as Sample C, and the goodnessof-fit between these new radial strain data sets and the original data generated for Sample C was plotted. The goodness-of-fit measure chosen was the reciprocal of the sum of the squared differences between the data sets. The ability of the curve-fitting algorithm to correctly identify the parameters when searching parameter space was then investigated by applying it to radial strain data that had been generated using the parameter values and radius of Sample C, and producing parametric images of quantities related to the material constants of the sample. To simulate the structure of future experimental data sets, the test data consisted of a 10 ⫻ 10 pixel matrix of radial strain values produced at each time instant by evaluating eqn 17, at 10 equidistant radial positions, and replicating these data in the axial direction to simulate 10 equidistant height positions. When supplied with an initial guess for the parameters [p1, p2, p3], the curve-fitting routine attempted to find the set of three values that best described the time-variation of the radial strain at each pixel in turn given the known value of the sample radius. The result of the fit for the previous pixel became the initial guess for the next pixel. To demonstrate the ability of the curve-fitting routine to converge accurately in spite of a poor initial guess for the fit parameters, the initial guess for the first pixel was chosen to be [p1, p2, p3] ⫽ [3 ⫻ 10⫺7 m2 s⫺1, 0.2, 0.7%] (the true values of the fit parameters of Sample C are given in Table 2). White noise was added to the strain data to simulate experimental conditions. The noise was created by al-

Acoustic model-based poroelastic imaging ● G. P. BERRY et al.

Fig. 2. Predictions of eqn 19 for the time-evolution of the radial-to-axial strain ratio on the solution plane for the three different poroelastic samples described in Table 1. The three samples differ only in their values of ␯s (nu) where ␯s ⫽ 0, 0.2 or 0.3. The time in seconds since the compression was applied is indicated, where t ⫽ 0⫹ s corresponds to the time instant immediately after the compression was applied. After 60,000 s of sustained compression, all three samples are close to equilibrium. The figure consists of a matrix of 15 images. The images on a given row of the figure describe the behaviour of a particular sample at various times; the images of a given column of the figure describe the behaviour of all three samples at a particular time. All images share the same colour scale. The leftmost edge of each image corresponds to r/a ⫽ 0.01 and the rightmost edge corresponds to r/a ⫽ 1.

lowing each radial strain value to fluctuate randomly over a total range equivalent to ⫾5% of the corresponding strain value.

555

time-variation of these strain ratios arises completely from the time-variation of the radial and circumferential strains, since the post-compression axial strain is not time-dependent. Immediately after the compression is applied, at time t ⫽ 0⫹ s, inside all three samples, the value of the circumferential-to-axial strain ratio is equal to 0.5 everywhere. However, although not obvious from Fig. 2, the value of the radial-to-axial strain ratio is equal to 0.5 everywhere except on the curved cylindrical surface (at r ⫽ a), where it is less than 0.5 for all three samples. On the curved cylindrical surface of Samples A, B and C, the value of the radial-to-axial strain ratio is 0, 0.125 and 0.214, respectively. As the compression is sustained, everywhere inside each sample both strain ratios tend toward the value of the Poisson’s ratio of the solid matrix. However, for a given sample, the timedependent behaviour of the radial-to-axial strain ratio differs from that of the circumferential-to-axial strain ratio, and this is especially obvious in the region near the curved cylindrical surface. Figure 4 shows the timeevolution of the average radial-to-axial strain ratio inside each sample, calculated in two ways: firstly, by using Armstrong’s original eqn 14 and, secondly, by obtaining the mean of the pixel values in the images of the radialto-axial strain ratio, which were produced using the new eqn 19 (some of these images are displayed in Fig. 2). Figure 4 proves that the predicted exotic (spatially varying) behaviour of the local radial-to-axial strain ratio inside the samples is entirely consistent with the previously predicted average behaviour of the sample as a whole (Armstrong et al. 1984).

RESULTS Results are provided below for the three sets of computational experiments, in the form of (1) visualisations, computed using the analytical model, of the timedependent strain field inside cylindrical samples of poroelastic materials undergoing stress relaxation, (2) equivalent visualisations computed using the finite element model, and (3) a preliminary demonstration of the feasibility of a method, based on the analytical model, to reconstruct parameters related to poroelastic material constants from noiseless and noisy strain data. Analytical theory Figure 2 shows half-plane images of the time-evolution of the radial-to-axial strain ratio, as predicted by eqn 19, for each of the three samples whose properties are described in Table 1. Figure 3 presents the corresponding images for the circumferential-to-axial strain ratio as predicted by eqn 20. Clearly, the 3D strain field predicted by the analytical model is both time-dependent and spatially varying. It is important to note that the

Fig. 3. Predictions of eqn 20 for the time-evolution of the circumferential-to-axial strain ratio on the solution plane for the three different poroelastic samples described in Table 1. The leftmost edge of each image corresponds to r/a ⫽ 0.01 and the rightmost edge corresponds to r/a ⫽ 1.

556

Ultrasound in Medicine and Biology

Fig. 4. Comparison of the average radial-to-axial strain ratio predicted by eqn 14 (lines) and the mean pixel values in the radial-to-axial strain images predicted by evaluating eqn 19 (points), for the samples in Table 1.

Finite element simulation Figure 5 shows images of the time-evolution of the radial-to-axial strain ratio, as predicted by the finite element model, for each of the three samples whose properties are described in Table 1. Note the good agreement with Fig. 2. In further agreement with the analytical model, the time-variation in this predicted strain ratio arose completely from the time-variation of the radial strain; the post-compression axial strain was predicted to be time-independent. In Fig. 6, the time-evolution of the radial-to-axial strain ratio in Sample A is plotted at three radial positions. The predictions of the analytical theory are superimposed on the plot to illustrate the degree of

Fig. 5. Predictions of FEM for the time-evolution of the radialto-axial strain ratio on the solution plane for the three different poroelastic samples described in Table 1. The leftmost edge of each image corresponds to r/a ⫽ 0.0013, and the rightmost edge corresponds to r/a ⫽ 0.9988. The height of each image is 50 mm. Note the agreement with Fig. 2.

Volume 32, Number 4, 2006

Fig. 6. Comparison of the predictions of the FEM and eqn 19 for a selection of radial positions for Sample A.

agreement between the analytical and finite element predictions. New structure predicted by both analytical theory and finite element simulation Figures 7 and 8 display new structure, not previously anticipated, in the time-dependent behaviour of the radial-to-axial strain ratio. In these figures, the predictions of the analytical theory and finite element model are superimposed for comparison. Figure 7 displays a plateau in the early time response of the radial-to-axial strain ratio for Sample A at radial positions close to the axis of symmetry of the cylinder (qualitatively similar early time responses were also obtained for Samples B and C). Figure 8 illustrates that the time-dependent ra-

Fig. 7. The early time plateau in the radial-to-axial strain ratio for internal points in Sample A as predicted by FEM and eqn 19.

Acoustic model-based poroelastic imaging ● G. P. BERRY et al.

Fig. 8. The radial-to-axial strain ratio overshoot for radial positions near the curved cylindrical surface in Sample C as predicted by FEM and eqn 19.

dial-to-axial strain ratio near the curved cylindrical surface for Sample C temporarily overshoots its equilibrium value (qualitatively similar behaviour was also observed for Sample B but regions near the edge of Sample A tended toward equilibrium without any overshoot). To our knowledge, these features of the time-dependent radial strain response of cylindrical poroelastic samples have not previously been described. Figure 9, obtained solely from the analytical model (eqns 6, 7, 17 and 18), illustrates that, even in an outer region of Sample C, in spite of the radial strain overshoot, the volumetric strain still decreases monotonically toward an equilibrium value. Curve-fitting algorithm Figure 10 illustrates the goodness-of-fit (Q) achieved when comparing the radial strain data generated for Sample C with the radial strain data generated for a large number of other samples with different fit parameter values to those of Sample C but with the same radius. It can be seen that the best fit occurs at a point in parameter space corresponding to the parameter set for Sample C. Figure 11 displays images of parameters p1, p2 and p3, produced by the curve fitting method described above, using noiseless strain data from Sample C and the same strain data mixed with ⫾5% additive white noise. Table 2 lists the known values of the fit parameters used to generate the radial strain data from eqn 17, and the mean and standard deviation of the pixel values in each image. DISCUSSION Analytical predictions The temporal and spatial variation of the strain field induced in three cylindrical samples of poroelas-

557

tic materials, when subjected to an instantaneously applied, sustained, unconfined axial compression, was investigated. The samples differed only in the values of the Poisson’s ratio of the solid matrix. The postcompression exudation of fluid, and the associated radial and circumferential contraction of the sample, caused a time-dependent, spatially varying volumetric strain (which is not shown, but may be inferred from Fig. 2 and Fig. 3) inside each sample, which was dependent on radial position and was independent of the axial position. Ultimately, the volumetric strain is directly related to fluid flow, since for a saturated sample of a poroelastic material with incompressible solid and fluid constituents, the compression-induced reduction in volume of the sample is numerically equal to the volume of the exudate. In the future, the availability of real-time 3D ultrasound scanners may enable direct measurement of this volumetric strain as a function of time and space. With current 2D scanners, only strains in a 2D plane can be imaged. It was, thus, logical in this paper to decompose the analysis of the volumetric strain behaviour into the behaviour of the radial and circumferential components, and to consider these components separately. Furthermore, given the anticipated difficulties in imaging circumferential strains with 2D ultrasound transducers, the predicted time-evolution of the radial strain and the radial-toaxial strain ratio was given more attention in this analysis than that of the circumferential strain. However, unlike previous work on poroelastography, the circumferential strain was considered in some detail, as its behaviour influences the volumetric strain.

Fig. 9. Volumetric, axial, radial and circumferential strains in an outer region of the cylinder (r/a ⫽ 0.99), in response to a ⫺1 % (compressive) axial strain applied to Sample C and sustained thereafter. The behaviour of the strains during the instantaneous application of the compression is not shown. Note that the volumetric strain decays monotonically with time and ␧ ⫽ ␧zz ⫹ ␧rr ⫹ ␧␪␪ for all times.

558

Ultrasound in Medicine and Biology

Volume 32, Number 4, 2006

Fig. 10. Goodness-of-fit between radial strain data generated for Sample C and data generated for samples with various other parameter combinations but with the same radius as Sample C. Figures (a), (b) and (c) were obtained using noiseless strain data and Figures (d), (e) and (f) were produced using data with added white noise.

At t ⫽ 0⫹ s, immediately after the compression is applied, the values of the strain ratios are uniformly distributed throughout the samples (except on the curved cylindrical surface at r ⫽ a). However, as the compression is sustained, strain ratios everywhere inside each sample tend toward an equilibrium value equal to the value of the Poisson’s ratio of the solid matrix of each sample. While this is shown in Fig. 2 and Fig. 3 for three samples, each with a distinct value of ␯s, it can be shown to be true in general by letting t¡⬁ in eqns 19 and 20 to yield: ␧rr ␧␪␪ r, t ⫽ r, t ⫽ ␯s ␧0 共 兲 ␧0 共 兲

as

t→⬁

(22)

Once relaxation has begun, at a given instant before a

given sample has reached equilibrium, the values of both strain ratios are consistently lower in regions closer to the curved cylindrical surface than in regions closer to the axis of symmetry. This results in an inward propagation of lower radial-to-axial strain ratio values as time progresses that is qualitatively consistent with that found in the finite element simulations of Konofagou et al. (2001). However, this does not agree with the results of the radial-to-axial strain ratio imaging experiments of Righetti et al. (2003, 2004), where the direction of this propagation is either not apparent in the images presented, or is seen to be in the opposite direction. It was, however, noted by Righetti et al. (2004) that their radialto-axial strain ratio data acquired near the lateral edges and in the central region of the scan plane were unreli-

Acoustic model-based poroelastic imaging ● G. P. BERRY et al.

559

Fig. 11. Parametric images of HAk, ␯s and ␧0 from simulated radial strain data during the sustained compression of Sample C. Figures (a), (b) and (c) were obtained using noiseless strain data and Figures (d), (e) and (f) were produced using data mixed with noise. The known parameter values used to generate the data and the mean pixel values in the parametric images are presented in Table 2.

able. Furthermore, any radial-dependence of the radialto-axial strain ratio behaviour was not discussed, and a comparison with the predictions of Konofagou et al. (2001) was not made. Since authors of some previous studies do not appear to have anticipated the radial-to-axial strain ratio behaviour described herein, we now provide a conceptual explanation for the physics of the phenomenon, as well as a comparison with previous work. The instantaneously applied compression causes immediate pressurisation of the pore fluid and radial expansion of the solid matrix. To balance the stress on the solid matrix due to the pore fluid pressure, an equal and opposite stress is generated in the solid matrix, which tends to cause radial contraction of the sample

(Armstrong et al. 1984). Since both the solid and fluid phases are assumed to be incompressible, radial contraction can only occur if an equivalent volume of fluid leaves the sample. However, fluid flow cannot occur instantaneously and, therefore, volume is conserved during this initial expansion. Figure 12 illustrates the time-evolution of the hydrostatic pressure in the pore fluid (pore pressure) as a function of radial position during the sustained compression of Sample A. As described in Appendix C, Fig. 12 was obtained both by evaluating an analytical expression describing the spatio-temporal behaviour of the pore pressure and also by using finite element simulation. In general, fluid flow is driven by a pressure gradient through Darcy’s Law. At t ⫽ 0⫹ s, the pore pressure is

560

Ultrasound in Medicine and Biology

Fig. 12. The behaviour of the hydrostatic fluid pressure as a function of radial position and time during the sustained compression of Sample A by a ⫺1% uniaxial compressive strain according to both eqn C1 (lines) and FEM (squares). No axial pressure gradient was predicted to exist. Note the agreement between the analytical and finite element predictions.

uniformly elevated across the radial extent of the sample, and falls abruptly to zero at the curved cylindrical surface. Thus, fluid flow does not simultaneously begin everywhere inside the sample. In response to the steep pressure gradient at the curved sample surface, flow begins first in regions just inside the curved cylindrical surface, causing the pore pressure in these outer regions to decay. This allows the restoring stress in the solid matrix to cause radial and consequently circumferential contraction of the solid matrix in these outer regions, and explains why the observed early time decay in the strain ratios is initially confined to these outer regions. Gradually, via a diffusion process, outward radial fluid flow induces a pore pressure gradient across a progressively larger proportion of the radial extent of the sample. This gradient encourages fluid flow from regions progressively closer to the axis of symmetry, which, in turn, permits volumetric strain to occur in these interior regions. Eventually, the pore pressure gradient extends across the whole radial extent of the sample and radial fluid flow occurs everywhere inside the sample. Thus, by a process of fluid exudation and fluid redistribution in the radial direction, the local volumetric strain inside the sample decays in a time-dependent, spatially varying manner resulting in the inward propagation of lower radial-to-axial strain ratio values as time progresses. Eventually, equilibrium is reached when the values of both strain ratios everywhere inside the samples equal the Poisson’s ratio of the solid matrix of each sample. The larger the ␯s, the faster equilibrium is reached (Fig. 4). This finding differs from the finite element predictions of Konofagou et al. (2001), who reported that

Volume 32, Number 4, 2006

the time taken to reach equilibrium is the same for two homogeneous poroelastic cylinders with the same values of Es and k, and that differ only in the value of ␯s, where ␯s ⫽ 0 or 0.3. However, by considering eqn 19, it is clear that ␯s affects the relaxation time. The influence of ␯s on the time taken for the radial-to-axial strain ratio to reach equilibrium is felt in two ways: firstly, ␯s affects HA by definition and, secondly, ␯s determines the roots of the characteristic equation (eqn 13). Both HA and the roots of the characteristic equation appear as terms in each exponent of the infinite sum, and the larger the ␯s, the larger the magnitude of a given exponent in the infinite sum, causing equilibrium to be reached more quickly. In fact, all three material constants should affect the relaxation time (Armstrong et al. 1984; Konofagou et al. 2001; Righetti et al. 2004). Thus, the time evolution of the relaxation cannot be used by itself to uniquely differentiate between materials with distinct values of k. Two classes of previously unanticipated structure were found in the time-dependent behaviour of the radial-to-axial strain ratio. Firstly, internal regions of the cylindrical samples exhibited a plateau in the early time response of the radial-to-axial strain ratio, and regions closer to the axis of symmetry maintained this plateau for longer than regions closer to the curved cylindrical surface (Fig. 7). This effect occurs because, as described above, fluid flow does not begin everywhere simultaneously; fluid flow occurs first at the curved cylindrical surface and subsequently occurs in regions progressively deeper inside the sample. Secondly, Fig. 8 predicts that the radial-to-axial strain ratio temporarily overshoots its equilibrium value in regions close to the curved cylindrical surface in Sample C. At t ⫽ 0⫹ s, immediately after the compression has been applied, the radial-toaxial strain ratio everywhere except at the cylindrical boundary is equal to 0.5 because of the instantaneous incompressibility of the poroelastic material. The analysis in Appendix D shows that, due to the boundary conditions of zero total radial stress and zero fluid pressure on the outer wall of the cylinder, the radial-to-axial ratio on the boundary at t ⫽ 0⫹ s is ␯s ⁄ 2共1 ⫺ ␯s兲 ⫽ 0.214. At equilibrium, when fluid flow eventually ceases, the radial-to-axial strain ratio everywhere inside the sample (including on the boundary) acquires the value of the Poisson’s ratio of the solid matrix equal to 0.3. Thus, on the boundary, the radial-to-axial strain ratio may be considered to overshoot the equilibrium value before tending gradually toward the value of the Poisson’s ratio of the solid matrix. For regions close to but not on the cylindrical surface, the overshoot is not instantaneous, has a smaller amplitude, and takes a finite time to occur. The physical mechanism behind this newly predicted behaviour may be related to the manner in which the time-dependent volumetric strain is redis-

Acoustic model-based poroelastic imaging ● G. P. BERRY et al.

tributed between the radial and circumferential strains according to the coupling between them as the sample contracts following the initial expansion. This view is supported by noting that the effect is not predicted to occur for Sample A (where ␯s ⫽ 0) since, in this case, ␯s ⁄ 2共1 ⫺ ␯s兲 ⫽ 0 as described in Appendix D. It is important to note that the behaviour of the compression-induced internal strain field predicted by this analysis is consistent with previous theoretical predictions. Appendix B demonstrates that eqn 12, identical to eqn 31 from Armstrong et al. (1984), may be recovered as a special case of eqn 16. Furthermore, Fig. 4 shows that, at a given instant, the mean of the pixel values in the images of the radial-to-axial strain ratio obtained using eqn 19 is equal to the value of the average radial-to-axial strain ratio predicted by Armstrong et al. (1984) at that same instant (see eqn 14). Thus, the exotic local strain field behaviour predicted by this analysis does not prevent the average radial-to-axial strain ratio in the sample from behaving in a manner identical to previous predictions (see Fig. 2 in Armstrong et al. 1984). The predictions of this analysis are also consistent with the experimental results of Fortin et al. (2003). Fortin et al. (2003) calculated the “effective Poisson’s ratio” (equivalent to the average radial-to-axial strain ratio inside the cartilage samples) by detecting the positions of the sample boundaries as a function of time. These time-dependent “effective Poisson’s ratio” data (see Fig. 2 in Fortin et al. 2003) displayed a monotonic decrease with time that is qualitatively similar to that predicted by Fig. 4. The predicted behaviour of the volumetric strain, obtained by combining eqns 6, 7, 17 and 18, is also consistent with previous theoretical work. Firstly, as demonstrated in Appendix E, the governing equation for the volumetric strain, obtained from eqn 10, agrees with that obtained by Biot (1941). This equation predicts that a process that is similar to heat conduction or diffusion governs the temporal and spatial variation of the volumetric strain. However, the radial strain is not required to obey the same equation. Secondly, as illustrated in Fig. 9, for a given radial position with r/a ⫽ 0.99, the predicted radial strain overshoot does not cause the volumetric strain to overshoot its equilibrium value. Instead, because the radial strain is accompanied by circumferential strain (due to the cylindrical geometry) and the axial strain, once applied, is held constant thereafter, the volumetric strain still decays in the monotonic fashion expected for an analogous diffusion or heat conduction process. Although Fig. 9 illustrates this for a single radial position, it is also true at all of the internal radial positions where radial strain overshoot is apparent. An overshoot in the radial-to-axial strain ratio at the curved cylindrical surface was not noted by Konofagou

561

et al. (2001). Different radial-to-axial strain ratio behaviour at the curved surface was instead predicted (see Figs. 3 and 7 in Konofagou et al. 2001), although this was based on an apparent misinterpretation of the meaning of eqn 14 (see below). However, although not noted by the authors, there is evidence of radial-to-axial strain ratio overshoot in some of their simulated images (see Fig. 9 in Konofagou et al. 2001). Eqn 14 has previously been used erroneously to describe the local radial-to-axial strain ratio as a function of time at the curved cylindrical surface (see Figs. 3 and 7, and eqn 4 in Konofagou et al. 2001). This interpretation is mathematically incorrect, as eqn 14 describes the local circumferential-to-axial strain ratio at the curved cylindrical surface, and is equivalent to the average radial-to-axial strain ratio in the material in the absence of shear strains. The variation of the radial-to-axial strain ratio as a function of time at the curved cylindrical surface is instead given by setting r ⫽ a in eqn 19, to obtain: ␧rr a, t ⫽ ␯s ⫹ ␧0 共 兲

⫻

J0共␣n兲 ⫺

⬁

兺

n⫽1

冉

冊

1 ⫺ 2␯s J ␣ 1 ⫺ ␯s 0共 n兲

1 J ␣ ⫺␣ J ␣ 共1 ⫺ ␯s兲 0共 n兲 n 1共 n兲

冉

exp

⫺␣n2HAkt a2

冊

(23)

Righetti et al. (2004) were the first to attempt to obtain agreement between experimental radial-to-axial strain ratio images and the predictions of the KLM biphasic theory. They approached this by comparing the timevariation of the radial-to-axial strain ratio measured at internal regions of a cylindrical sample of tofu with predictions made by a combination of finite element modeling and eqn 14. The radius of the sample that they used experimentally was approximately 37.5 mm, but the value used for the radius when eqn 14 was evaluated was 15 mm (see Fig. 11 in Righetti et al. 2004). The authors’ reasons for choosing this value of the radius for the purpose of evaluating eqn 14 were not stated. However, it is noted by us that the regions from which the experimental strain data were obtained were approximately symmetrically distributed about the axis of symmetry, and lie approximately 15 mm from this axis. This implies that the authors interpreted eqn 14 as describing the radial-to-axial strain ratio at a particular radial position, whereas in fact it describes the average radial-to-axial strain ratio over the entire radial extent of the cylinder (in the absence of shear strains). We note that a comparison between theory

562

Ultrasound in Medicine and Biology

Fig. 13. The predicted behaviour of (a) the average radial-toaxial strain ratio (using eqn 14) for a cylindrical sample with Es ⫽ 1.29 kPa, k ⫽ 1.373 ⫻ 10⫺11 m4 N⫺1 s⫺1, ␯s ⫽ 0.3 and a ⫽ 15 mm, (b) the average radial-to-axial strain ratio (using eqn 14) for a cylindrical sample with Es ⫽ 1.29 kPa, k ⫽ 1.373 ⫻ 10⫺11 m4 N⫺1 s⫺1, ␯s ⫽ 0.3 and a ⫽ 37.5 mm, (c) the local radial-to-axial strain ratio (using eqn 19) at a radial position 15 mm from the axis of symmetry in a cylindrical sample with Es ⫽ 1.29 kPa, k ⫽ 1.373 ⫻ 10⫺11 m4 N⫺1 s⫺1, ␯s ⫽ 0.3 and a ⫽ 37.5 mm.

and experiment for the time-dependent behaviour of the radial-to-axial strain ratio at a particular radial position is possible but, in any such comparison, eqn 19 should be used instead of eqn 14, and the value of r in the equation should equal the radial position at which the radial-to-axial strain ratio data were obtained. Alternatively, if the required data are available, a comparison between theory and experiment for the behaviour of the average radial-to-axial strain ratio across the entire radial extent of the sample can also be made, but such a comparison is only valid when eqn 14 and the correct value for the sample radius is used. Figure 13 illustrates the importance of using the correct equations for the purpose of comparison as well as the correct input value for the radius of the sample. Note that using the correct equation and the correct radius influences the time taken to reach equilibrium. Note also that using the correct equation for the local variation of the radial-to-axial strain ratio predicts the existence of the early time plateau discussed above. Furthermore, although not shown here, incorrectly using eqn 14 in the manner described above in an attempt to describe the behaviour of the local radial-toaxial strain ratio at radial positions inside the sample predicts that strain ratio decay begins first near the axis of symmetry and that lower values of the strain ratio propagate outward toward the curved cylindrical surface.

Volume 32, Number 4, 2006

Finite element predictions To validate quantitatively the numerical accuracy of the finite element model, the unconfined stress-relaxation behaviour of the three poroelastic samples described in Table 1 was simulated using ABAQUS, and the predictions were compared with those of eqn 19. This numerical validation of finite element modeling predictions of the local radial-to-axial strain ratio behaviour throughout the sample has not previously been carried out, and is especially important given that previous experimental investigations (Righetti et al. 2003, 2004) have failed to confirm some of the predictions provided by a previous finite element simulation (Konofagou et al. 2001). The mesh was chosen to be as dense as possible in the radial direction and the time-steps were chosen to be suitably small in order to deliver good approximations to the analytical solutions in an acceptable time. As shown in Figs 2, 5, 6, 7 and 8, good agreement was found with the predictions of the analytical theory. Both analytical and finite element methods were in agreement with regard to both types of previously unanticipated structure in the radial-to-axial strain ratio behaviour. This suggests that the finite element model is a faithful approximation to the analytical theory under the same boundary conditions. As a result of this numerical validation, this finite element method may be used in the future, with more confidence than hitherto possible, to examine the behaviour of the strain field in more general poroelastic materials undergoing more general deformation. The application of the finite element method to solve such generalised poroelastic deformation problems is currently under investigation. Parametric images of the poroelastic material constants Previous poroelastographic investigations have focused on either predicting or imaging the time-dependent behaviour of the radial displacement or the radial-toaxial strain ratio during stress relaxation. Additionally, it has been suggested that an image of the Poisson’s ratio of the solid matrix may be produced by waiting for a sufficiently long time so that equilibrium is reached and the strain ratio everywhere inside the material acquires the value of the Poisson’s ratio of the solid matrix (Konofagou et al. 2001; Righetti et al. 2004). This approach is practical only when the time taken to reach equilibrium is suitably short. Suggestions as to how to image the permeability or the Young’s modulus of the solid matrix, and more practical methods of imaging the Poisson’s ratio of the solid matrix, have not previously been made. In this paper, a model-based iterative method was developed to produce images of parameters related to the material constants of a poroelastic sample of known radius, using knowledge of the expected time-dependent

Acoustic model-based poroelastic imaging ● G. P. BERRY et al.

radial strain during stress relaxation. The uniqueness with which a particular radial strain behaviour could be generated in a sample with a given set of fit parameters and sample radius was first investigated. As seen in Fig. 10, for the region of parameter space sampled, a sample of known radius with a given set of fit parameters produces a unique strain behaviour that cannot be generated by another sample with the same radius and any other combination of fit parameter values. The inversion method was then tested using simulated radial strain data, which successfully produced a set of three simulated parametric images with mean pixel values that were in close agreement with the known values (see Fig. 11 and Table 2). The technique was also found to be relatively tolerant to noise. It was not possible to separate HA and k by this curve-fitting technique using strain data alone. Attempting to do so would lead to non-uniqueness of the set of best-fit parameters for a given radial strain data set. However, parameter separation may be achieved experimentally if HA were known from the equilibrium behaviour. To the authors’ knowledge, this curve-fitting approach provides the first demonstration of a means of producing parametric images directly related to the permeability, aggregate modulus and Poisson’s ratio of a poroelastic material, and we believe it warrants further study as a component of an experimental poroelastic imaging method. Consequences for future experimental poroelastography The 3D strain field inside a homogeneous cylinder of a poroelastic material is clearly predicted to be temporally and spatially dependent (Fig. 2 and Fig. 3). Thus, even when imaging a homogeneous sample, the spatial dependence of the strain field may require the use of high spatial resolution strain imaging. This is especially true if the radial strain overshoot in regions close to the curved cylindrical surface is to be detected. Because of the nature of the temporal dependence of the strain field, a variable time interval between acquired ultrasound RF frames is to be recommended if computer memory for storing the data is limited—more frames being acquired at early times than at later times. This will enable the early time rapid variation of the strain field to be captured while also allowing the less rapid strain field variation at later times to be efficiently recorded. It has previously been pointed out that when post-compression equilibrium is reached, the image of the radial-to-axial strain ratio is an image of the Poisson’s ratio of the solid matrix ␯s (Konofagou et al. 2001; Righetti et al. 2004). The model-based approach described here enables images of the Poisson’s ratio of the solid matrix to be produced even when the strain data available do not extend over a sufficiently long time interval to include the time when equilibrium is reached everywhere inside

563

the material. This may render practical the poroelastic imaging of materials with relatively long relaxation times, and may eventually assist in vivo imaging where patient scan duration is an important consideration. Elastography can produce high precision strain estimates along the direction of beam propagation, but strain estimates in the direction perpendicular to beam propagation (but in the scan plane) are inferior (Konofagou et al. 1998). Since the curve-fitting algorithm performs well even when fed with a single (radial) component of the strain field, the direction of beam propagation can be experimentally arranged to coincide with the radial direction, thus potentially enabling parametric imaging without the need to utilise inferior strain estimates in vitro. Clearly, this transducer orientation is less practical in vivo. However, if experimental conditions cannot be arranged to exactly match the boundary condition assumptions of this analysis, then both radial and axial strain data may need to be considered. Although not demonstrated here, if radial-to-axial strain ratio data are available, it is possible to perform a two-parameter fit to these data using eqn 19, where the two fit parameters p1 and p2 are given by: p1 ⫽ HAk p2 ⫽ ␯S

(24)

Although not illustrated here, we have found that the uniqueness, accuracy and precision properties of the two-parameter fit are similar to those achieved using the three-parameter fit. While the poroelastic imaging technique described here was discussed from an ultrasound perspective, the technique itself is imaging modality independent. We propose that magnetic resonance imaging (MRI) could prove to be a useful poroelastic imaging modality since magnetic resonance elastography (MRE) has already been used to perform 3D strain imaging (Plewes et al. 2000; Weaver et al. 2001). Clinical relevance of biphasic theory and parametric imaging The three material constants of the biphasic theory are clearly related to material properties of clinical relevance: Es is a measure of the elasticity of the poroelastic solid matrix; the KLM biphasic theory assumption of phase incompressibility at physiological pressures allows ␯s to be related to the volume of fluid displaced by the applied compression, which may be indicative of the volume fraction of mobile fluid initially residing in the material (sometimes known as the porosity of the solid matrix); the permeability k is a measure of the ease with which fluid may pass through the poroelastic material. The model-based imaging technique described here pro-

564

Ultrasound in Medicine and Biology

duces two images related to these three material constants and so may be clinically relevant in the assessment of lymphoedema, as described in the introduction. While this model-based approach has the potential to be applied in vitro in its current form with suitable phantom materials under controlled laboratory conditions, in vivo clinical use will require a much greater degree of sophistication. Potential problems associated with in vivo quantitative poroelastography include boundary condition departures from those assumed in this analysis, tissue inhomogeneities, inherent viscoelasticity, anisotropy and nonlinearity of the solid matrix, strain dependence of the permeability and patient movement. Furthermore, since different fluid flow paths through the tissue are each likely to possess distinct characteristic time constants, it may also be important to account for the different scales of tissue permeability arising from the possibility of fluid flow through several physiological compartments as well as the possibility of inter-compartmental fluid exchange (Netti et al. 1997). The use of finite element modeling may hold promise for dealing with some of these challenges, which is one of the reasons why we chose to numerically validate a finite element model as part of the present study. Finally, as has been seen throughout the history of medical imaging, although such factors may indeed make quantitative poroelastography difficult to achieve, methods that display qualitative differences between tissues in highly complex images can often be of great medical value, even although they may require considerable expertise to interpret them. An ability to incorporate a multi-compartmental structure into a generalised poroelastic imaging method could lead to improved diagnosis and treatment of cancer. Such an imaging method could provide estimates of physiological parameters such as the micro-vessel density (perfusion), the ease with which fluid may cross the micro-vessel wall into or out of the interstitium, the micro-vessel tortuosity and the viscosity of the blood in the microvasculature. Relative to healthy tissue, spatially varying perfusion, enhanced filtration across the vessel wall, increased micro-vessel tortuosity and higher blood viscosity are characteristic of some malignant tissue (Jain 1994; Kuszyk et al. 2001; Leach 2001). Higher microvessel densities have also been linked with an increased metastatic risk (Folkman 1995; Kuszyk et al. 2001). In addition, the microvessel density and associated interstitial pressure gradients are believed to influence the efficiency with which chemotherapeutic agents are dispersed throughout the tumour mass (Folkman 1995; Jain 1994; Kuszyk et al. 2001). Thus, a generalised poroelastic imaging method may have the potential to discriminate between healthy and malignant tissue, may identify tumours with a greater probability of metastatic

Volume 32, Number 4, 2006

spread and may be applied to predict which tumours are likely to respond to antiangiogenic and antivascular treatments (Leach 2001). Furthermore, such an imaging method may eventually be used to identify tumours that are likely to resist penetration by chemotherapeutic agents in general, and so might best be treated by other means. CONCLUSION New analytical expressions were derived that describe, in terms of the poroelastic material constants, the time-evolution of the strain field in homogeneous cylindrical samples of a poroelastic material undergoing stress relaxation. The equations were used for three purposes. Firstly, they were employed to establish the nature of the time-dependent strain field inside poroelastic materials undergoing stress relaxation where the accompanying boundary condition assumptions are satisfied. New structure in the strain field, not previously predicted, was found. Secondly, they were used to numerically validate a finite element model under the same boundary condition assumptions. The finite element model may now be used with more confidence to predict the behaviour of the strain field under less restrictive boundary condition assumptions, where the analytical expressions are not applicable. Finally, the expression for the time-evolution of the radial strain was used in a model-based iterative method for reconstructing images of parameters related to the poroelastic material constants. By testing this routine with simulated noiseless and noisy strain data, it was shown to be capable of solving the inverse problem of using measurements of the time-dependent radial strain to image parameters related to the underlying poroelastic material constants. This algorithm appears to have potential for use as part of an experimental modelbased poroelastic parametric imaging system. A companion paper will demonstrate the extent to which this is possible in practice in vitro. Acknowledgments—This work was supported by grants from the Institute of Cancer Research and the Engineering and Physical Sciences Research Council, UK. P.E.B. gratefully acknowledges financial support from NSF.

REFERENCES Armstrong CG, Lai WM, Mow VC. An analysis of the unconfined compression of articular cartilage. J Biomech Eng 1984;106:165– 173. Atkins GL. Multicompartment models for biological systems. London, United Kingdom: Methuen and Co. Ltd., 1969. Barbone PE, Bamber JC. Quantitative elasticity images: What can and cannot be inferred from strain images. Phys Med Biol 2002;47: 2147–2164. Bear J. Dynamics of fluids in porous media. New York: Dover Publications, 1988. Berry GP, Bamber JC, Armstrong C, Miller NR. Towards a modelbased poroelastic imaging method. Proceedings of the Third Inter-

Acoustic model-based poroelastic imaging ● G. P. BERRY et al. national Conference on the ultrasonic measurement and imaging of tissue elasticity. October 17–20, Cumbria, UK: 2004a:96. Berry GP, Bamber JC, Armstrong C, Miller NR. Ultrasonic strainbased imaging of tissue permeability and mobile fluid volume. Ultrasound 2004b;12(4):232. Biot MA. General theory of three-dimensional consolidation. J Applied Physics 1941;12:155–164. Biot MA. Mechanics of deformation and acoustic propagation in porous media. J Applied Physics 1962;33(4):1482–1498. Boas ML. Mathematical methods in the physical sciences. 2nd ed. New York: Wiley, 1983. Bowles JE. Engineering properties of soils and their measurements. New York: McGraw Hill, 1992. Brady BHG, Brown ET. Rock mechanics for underground mining. Melbourne, Victoria, Australia: Allen and Unwin Ltd., 1985. Dickinson RJ, Hill CR. Analysis of tissue and organ dynamics. In: Hill CR, Alvisi C, eds. Investigative ultrasonology 1. Technical advances. Pitman Medical, Tunbridge Wells, Kent, England: 1980: 110 –114. Dixon JM, Mensel RE. ABC of breast diseases: Symptoms, assessment and guidelines for referral. BMJ 1994;309:722–726. Doyley MM, Meaney PM, Bamber JC. Evaluation of an iterative reconstruction method for quantitative elastography. Phys Med Biol 2000;45:1521–1540. Duck FA. Physical properties of tissue: A comprehensive reference book. London, UK: Academic Press Ltd., 1990. Folkman J. Clinical applications of research on angiogenesis. N Engl J Med 1995;333(26):1757–1763. Fortin M, Buschmann MD, Bertrand MJ, Foster S, Ophir J. Dynamic measurement of internal solid displacement in articular cartilage using ultrasound backscatter. J Biomechanics 2003;36:443– 447. Fung YC. Biomechanics: Mechanical properties of living tissues, 2nd ed. New York: Springer-Verlag, Inc., 1993. Harrison TR (ed-in-chief). Principles of internal medicine, Vol II, 12th ed. 1991. Hibbitt, Karlsson and Sorensen, Inc. ABAQUS Theory Manual (Version 6.3). 2003. Jaeger JC, Cook NGW. Fundamentals of rock mechanics, 2nd ed. New York: John Wiley, 1976. Jain RK. Barriers to drug delivery in solid tumors. Scientific American 1994;July:42– 49. Kallel F, Bertrand M. Tissue elasticity reconstruction using linear perturbation method. IEEE Trans Med Imag 1996;15:299 –313. Konofagou E, Ophir J. A new elastographic method for estimation and imaging of lateral displacements, lateral strains, corrected axial strains, and Poisson’s ratio in tissues. Ultrasound Med Biol 1998; 24:1183–1199. Konofagou E, Harrigan TP, Ophir J, Krouskop TA. Poroelastography: Imaging the poroelastic properties of tissues. Ultrasound Med Biol 2001;27(10):1387–1397. Kuszyk BS, Corl FM, Franano FN, et al. Tumour transport physiology: Implications for imaging and image-guided therapy. Am J Roentgen 2001;177:747–753. Lakes RS. Viscoelastic solids. Boca Raton, FL: CRC Press. 1998. Lanczos C. Applied analysis. London, UK: Sir Isaac Pitman and Sons Ltd., 1957. Leach MO. Application of magnetic resonance imaging to angiogenesis in breast cancer. Breast Cancer Res 2001;3:22–27. Levick JR. An Introduction to cardiovascular physiology, 4th ed. London: Arnold Publishers, 2003. Moseley A, Piller N, Carati C. Combined opto-electronic perometry and bioimpedance to measure objectively the effectiveness of a new treatment intervention for chronic secondary leg lymphedema. Lymphology 2002;35:136 –143. Mow VC, Holmes MH, Lai WM. Fluid Transport and mechanical properties of articular cartilage: A review. J Biomechanics 1984; 17(5):377–394. Mow VC, Kuei SC, Lai WM, Armstrong CG. Biphasic creep and stress relaxation of articular cartilage in compression: theory and experiments. J Biomech Eng 1980a;102:73– 84.

565

Mow VC, Lai WM. Recent developments in synovial joint biomechanics. SIAM Review 1980b;22:275–317. Mow VC, Ratcliffe A, Woo SLY, eds. Biomechanics of Diarthrodial Joints, Vol I. New York: Springer Verlag, 1990. Netti PA, Baxter LT, Boucher Y, Skalak R, Jain RK. Macro- and microscopic fluid transport in living tissues: Application to solid tumours. American Institute of Chemical Engineers Journal 1997; 43(3):818 – 834. Oberai AA, Gokhale NH, Doyley MM, Bamber JC. Evaluation of the adjoint equation based algorithm for elasticity imaging. Phys Med Biol 2004;49:2955–2974. Ophir J, Cespedes EI, Ponnekanti H, Yazdi Y, Li X. Elastography: A quantitative method for imaging the elasticity of biological tissues. Ultrasonic Imaging 1991;13:111–134. Plewes DB, Bishop J, Samani A, Sciarretta J. Visualization and quantification of breast cancer biomechanical properties with magnetic resonance elastography. Phys Med. Biol 2000;45:1591–1610. Roberts CC, Levick JR, Stanton AWB, Mortimer PS. Assessment of truncal edema following breast cancer treatment using modified harpenden skinfold calipers. Lymphology 1995;28:78 – 88. Righetti R, Ophir J, Krouskop TA. The feasibility of using elastography for imaging the local lateral-to-axial strain ratios in homogeneously porous media. Proceedings of the Second International Conference on the Ultrasonic Measurement and Imaging of Tissue Elasticity. October 12–15, 2003. Corpus Christi, TX, 2003:38. Righetti R, Ophir J, Srinivasan S, Krouskop TA. The feasibility of using elastography for imaging the Poisson’s ratio in porous media. Ultrasound Med Biol 2004;30(2):215–228. Stanton AWB, Badger C, Sitzia J. Non-invasive assessment of the lymphedematous arm. Lymphology 2000;33:122–135. Talley NJ, O’Connor S. Clinical examination: A guide to physical diagnosis. NSW Australia: Williams and Wilkins and Associates Pty Ltd., 1988. Wang HF. Theory of linear poroelasticity. Princeton, NJ: Princeton University Press, 2000. Weaver JB, Van Houten EEW, Miga MI, Kennedy FE, Paulsen KD. Magnetic resonance elastography using 3D gradient echo measurements of steady-state motion. Med Phys 2001;28(8):1620 –1628. Wellman P, Howe RH, Dalton E, Kern KA. Breast tissue stiffness in compression is correlated to histological diagnosis. Harvard BioRobotics Laboratory, Technical Report 1999. Wellman PS, Dalton EP, Krag D, Kern KA, Howe RD. Tactile imaging of breast masses: First clinical report. Arch Surgery 2001;136(2): 204 –208. Worseley BH. Analysis of decay-type data. Communications of the Association for Computing Machinery 1964;7(1):39 – 44.

APPENDIX A By looking down the axis of symmetry, the outer circumference of a cylinder projects onto a 2D circle. In the absence of shear strains, the radial expansion of this outer circumference in response to an axial compression may be modelled in terms of the radial expansion of a circle. The circumferential strain experienced by a circle when it expands from a radius a to a radius a ⫹ u may be obtained by considering the initial and expanded circumferences, cinitial and cexpanded, respectively, which are given by: cinitial ⫽ 2␲a cexpanded ⫽ 2␲共a ⫹ u兲 The change in circumference leads to a circumferential strain. At the outer surface, where the radial position r is equal to a, the circumferential strain, ⑀␪␪, is given by: ␧␪␪ ⫽

2␲共a ⫹ u兲 ⫺ 2␲a 2␲a

(A1)

u a

(A2)

␧␪␪ ⫽

Thus, in the absence of shear strains, the circumferential strain at the

566

Ultrasound in Medicine and Biology

curved cylindrical surface is equal to the average radial strain inside the cylinder. Eqn A2 may be generalised by the same geometrical reasoning to obtain an expression that applies at all internal radial positions in the cylinder and not just at the outer cylindrical surface. The generalised form of eqn A2, valid for all values of radial position, r, is: ␧␪␪(r) ⫽

u共r兲 r

Volume 32, Number 4, 2006 al. (1984) according to the material constants and radius of Sample A, by finite element simulation, or by evaluating the following expression for the pore fluid pressure, p, as a function of r and t: p共r, t兲 ⫽

(A3)

In the absence of shear strains, eqn A3 describes the average radial strain in a cylindrical region of the cylinder, where the radial extent of the cylindrical region extends from the axis of symmetry to the radial position r.

APPENDIX B

⫻

␧0Es 2共1 ⫹ ␯s兲

冤

⬁

兺

n⫽1

冉

冉 冊 兲冊 共

␣nr ⫺ J 0 共 ␣n 兲 ␣2nHkt a exp ⫺ 2 a 1 ␣n J0 ␣n兲 ⫺ J1共␣n兲 2 2共1 ⫺ ␯s J0

冤

␴trr ⫽ ␭s␧ ⫹ 2␮s␧rr ⫺ p

冉

冊再 冎 冉 exp

再

J1共␣n兲 ⬁ 共1 ⫺ ␯s兲 ⫺ 2共1 ⫺ 2␯s兲 ␣nJ0共␣n兲 ␯s ⫹ J1共␣n兲 n⫽1 1 ⫺ ␣n共1 ⫺ ␯s兲 J0共␣n兲

兺

⫺␣2nHAkt a2

冎 冉 exp

冊冥

where ␭s and ␮s are the Lamé constants of the solid matrix, which may be expressed in terms of Es and ␯s as follows: (B1)

⫺␣2nHAkt a2

冊冥

(B2)

␭s ⫽

Es␯s

共1 ⫹ ␯s兲共1 ⫺ 2␯s兲

␣nJ1共␣n兲 ⫽ ␣2n J0共␣n兲

共1 ⫺ ␯s兲 共1 ⫺ 2␯s兲

(B3)

Substitution of eqns B3 into eqn B2 gives: ⬁

u共a, t兲 ⫽ a␧0 ␯s ⫹ ⫻

兺

␴rrt共r ⫽ a, t兲 ⫽ 0

p共r ⫽ a, t兲 ⫽ 0

n⫽1

冉

⫺␣ HAkt ⫻ exp a2 2 n

冤

u a, t ⫽ ␧0 ␯s ⫹ 共1 ⫺ 2␯s兲共1 ⫺ ␯s兲 a共 兲

兺

n⫽1

冉

(D4)

⫹

As described earlier, at t ⫽ 0 s, immediately after the compression has been applied, ␧␪␪ ⫽ (1/2)␧0 everywhere inside the sample, including on the curved cylindrical surface, and so we may rewrite eqn D4 as follows:

冉

␧rr ⫹

冊

␧0 ⫺ ␧0 ⫹ 2␮s␧rr ⫽ 0 2

冊册

␧rr ␭s ⫽ ␧0 2共␭s ⫹ 2␮s兲 (B4)

冊

⫺␣2nHAkt a2 2 2 ␣ 1 ⫺ ␯ ⫺ 共1 ⫺ 2␯s兲其 兵 n共 s兲 exp

(D3)

(D5)

Rearranging, we obtain:

From which eqn 12 may be obtained as follows: ⬁

(D2)

At the outer edge of the sample, by substitution of eqns D3 into eqn D1, we have:

␭s

共1 ⫺ ␯s兲共1 ⫺ 2␯s兲 ⫺ 2共1 ⫺ 2␯s兲共1 ⫺ ␯s兲 共1 ⫺ 2␯s兲 ⫺ ␣2n共1 ⫺ ␯s兲2

Es 2共1 ⫹ ␯s兲

␭s␧ ⫹ 2␮s␧rr ⫽ 0

J1共␣n兲 共1 ⫺ ␯s兲 ⫽ ␣nJ0共␣n兲 共1 ⫺ 2␯s兲

and

␮s ⫽

Zero total radial stress and zero pore pressure boundary conditions exist on the curved cylindrical surface of the sample as follows:

But, from eqn 13, which the ␣n roots must satisfy, we have:

冋

(D1)

n⫽1

冉 冊

冤

APPENDIX D

兺

1 ⫺ 2␯s J1共␣n兲 J0共␣n兲 ⫺ 2 1 ⫺ ␯s ␣n ⫻ 1 J ␣ ⫺ ␣nJ1共␣n兲 1 ⫺ ␯s 0共 n兲

⫽a␧0

(C1)

The total radial stress, ␴rrt, in the sample is given as follows (Armstrong et al. 1984):

⬁

␯s ⫹

冊冥

Eqn C1 was obtained by employing the Heaviside expansion theorem, in the manner described by Armstrong et al. (1984)), to perform the Laplace inversion of eqn 26 in Armstrong et al. (1984).

This appendix illustrates how eqn 12 may be derived from eqn 16. By setting r ⫽ a, eqn 16 reduces to:

u共a, t兲 ⫽ a␧0

冉

冥

(B5) Thus, eqn 16 reduces to eqn 12 when r ⫽ a as required.

(D6)

Using eqns D2 in eqn D6, we obtain the following expression for the radial-to-axial strain ratio on the curved cylindrical surface of the sample at t ⫽ 0⫹ s: ␧rr ␯s ⫽ ␧0 2共1 ⫺ ␯s兲

(D7)

The quantity on the right hand side of eqn D7 is always less than ␯s for 0 ⬍ ␯s ⬍ 0.5 and for ␯s ⫽ 0, it is equal to 0. Thus, no overshoot occurs for ␯s ⫽ 0. Note also that, by adapting the derivation above or by using finite element analysis, it may be shown that overshoot does not occur in the case of plane-strain conditions.

APPENDIX C

APPENDIX E

The radial distribution of the pore pressure at various times for a cylindrical sample in non-dimensionalised units has previously been presented (see Fig. 3 in Armstrong et al. 1984). The pore pressure distribution for Sample A may be obtained by using any one of three techniques: by dimensionalising the figure presented by Armstrong et

The governing equation for volumetric strain using poroelastic theory For a poroelastic material composed of incompressible solid and fluid constituents, and with a linear elastic, isotropic, solid matrix saturated with a pore fluid that flows through the porous material according to Darcy’s Law, Biot (1941) reported that the following

Acoustic model-based poroelastic imaging ● G. P. BERRY et al. ⭸ ⭸r

equation governs the behaviour of the volumetric strain in the solid matrix, for arbitrary applied stress or strain histories: ⵜ2␧ ⫽

1 ⭸␧ c ⭸t

(E1)

冉

冊

⭸u u 1 ⫹ ⫽ ⭸r r HAk

⭸␧ 1 ⫽ ⭸r HAk (E2)

1 ⭸ r ⭸r

冉 冊 r

⭸ 1 ⭸2 ⭸2 ⫹ 2 2⫹ ⭸r r ⭸␪ ⭸z2

(E3)

共1 ⫹ ␯s兲共1 ⫺ 2␯s兲 Es共1 ⫺ ␯s兲

1 ⭸␧ HAk ⭸t

The governing equation for volumetric strain using the KLM biphasic theory Eqn 10 may be rewritten as follows:

ⵜ 2␧ ⫽

r

(E5)

(E6)

and the governing equation (eqn E1) then becomes: ⵜ 2␧ ⫽

冉

⭸u r ⭸␧zz ⫹ ⭸t 2 ⭸t

1 ⭸ r ⭸r

From eqn E9:

Substituting eqn 11 into eqn E5 and then using E2, we find that the Biot consolidation constant may be rewritten as follows: c ⫽ H Ak

(E8)

冊

(E9)

(E7)

(E10)

and thus, when operating on the volumetric strain, eqn E4 reduces to:

(E4)

Note that eqn E1 is given here in a different form to that published by Biot (see eqn 6.8 in Biot 1941). We ascribe the difference to a typographical error in Biot (1941), as can be seen by considering eqns 4.4 and 6.5 in Biot (1941) and noting that Biot (1941) himself refers to this equation as the equation of heat conduction. Combining eqns E3 and D2 we see that: a⫽

冊

⭸2␧ ⭸2␧ ⫽ ⫽0 ⭸ ␪2 ⭸z2

and the Laplacian operator, ⵜ2, is given in cylindrical polar coordinates as follows (Boas 1983): ⵜ2 ⫽

⭸u r ⭸␧zz ⫹ ⭸t 2 ⭸t

where ␧ is the volumetric strain. Because of the rotational symmetry of the problem and since the applied strain is independent of z, we have:

a is the Biot “final compressibility” given by: 1 ⫺ 2␯s , a⫽ 2␮s共1 ⫺ ␯s兲

冉

But using eqns 7, 8 and 9, and recognising that the applied axial strain is independent of the radial position, we may write:

where c is the Biot “consolidation constant” given by: k c⫽ , a

567

⭸ ⭸r

⭸␧ 1 ⭸ ⫽ ⭸r HAk ⭸t

冉 冊 冉 冊 冉 冊 r

⭸␧ 1 ⭸ ⫽ ⭸r HAk ⭸t

⭸ ⭸r

r

1 ⭸ r ⭸r

冉 冊 r

冉 冉

⭸␧ ⭸r

ru ⫹

u⫹r

⭸␧ 1 ⭸ ⫽ ⭸r HAk ⭸t r

r2 ␧ 2 zz

(E11)

冊

⭸u ⫹ r␧zz ⭸r

(E12)

冊

共r␧兲

⭸␧ 1 ⭸␧ ⫽ ⭸r HAk ⭸t

(E13) (E14) (E15)

Using eqn E15 in E11, the following equation is obtained: ⵜ2 ␧ ⫽

1 ⭸␧ HAk ⭸t

(E16)

A comparison of eqns E7 and E16 demonstrates that the governing equation for the temporal and spatial evolution of the volumetric strain in the cylindrical sample, derived from eqn 10, is identical to that obtained from the theory of poroelasticity. Note that it is the volumetric strain that is governed by this equation and not the radial strain. The equation is identical in form to the diffusion or heat conduction equation. Note that, although not shown here, the governing equation for the volumetric strain is quite different in the case where the material constants are spatially varying.