Towards an acoustic model-based poroelastic imaging method: II. experimental investigation

Ultrasound in Med. & Biol., Vol. 32, No. 12, pp. 1869 –1885, 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.07.013

● Original Contribution TOWARDS AN ACOUSTIC MODEL-BASED POROELASTIC IMAGING METHOD: II. EXPERIMENTAL INVESTIGATION GEARÓID P. BERRY,* JEFFREY C. BAMBER,* NAOMI R. MILLER,* PAUL E. BARBONE,† NIGEL L. BUSH,* and CECIL G. ARMSTRONG‡ *Joint Department of Physics, Institute of Cancer Research and Royal Marsden NHS Foundation Trust, Sutton, Surrey, UK; †Aerospace & Mechanical Engineering, Boston University, Boston, MA, USA; and ‡School of Mechanical and Aerospace Engineering, QUB, Ashby Building, Belfast, UK (Received 6 February 2006; revised 19 June 2006; in final form 13 July 2006)

Abstract—Soft biological tissue contains mobile fluid. The volume fraction of this fluid and the ease with which it may be displaced through the tissue could be of diagnostic significance and may also have consequences for the validity with which strain images can be interpreted according to the traditional idealizations of elastography. In a previous paper, under the assumption of frictionless boundary conditions, the spatio-temporal behavior of the strain field inside a compressed cylindrical poroelastic sample was predicted (Berry et al. 2006). In this current paper, experimental evidence is provided to confirm these predictions. Finite element modeling was first used to extend the previous predictions to allow for the existence of contact friction between the sample and the compressor plates. Elastographic techniques were then applied to image the time-evolution of the strain inside cylindrical samples of tofu (a suitable poroelastic material) during sustained unconfined compression. The observed experimental strain behavior was found to be consistent with the theoretical predictions. In particular, every sample studied confirmed that reduced values of radial strain advance with time from the curved cylindrical surface inwards towards the axis of symmetry. Furthermore, by fitting the predictions of an analytical model to a time sequence of strain images, parametric images of two quantities, each related to one or more of three poroelastic material constants were produced. The two parametric images depicted the Poisson’s ratio (␯s) of the solid matrix and the product of the aggregate modulus (HA) of the solid matrix with the permeability (k) of the solid matrix to the pore fluid. The means of the pixel values in these images, ␯s ⴝ 0.088 (standard deviation 0.023) and HAk ⴝ 1.449 (standard deviation 0.269) ⴛ 10–7 m2 s–1, were in agreement with values derived from previously published data for tofu (Righetti et al. 2005). The results provide the first experimental detection of the fluid-flow-induced characteristic diffusion-like behavior of the strain in a compressed poroelastic material and allow parameters related to the above material constants to be determined. We conclude that it may eventually be possible to use strain data to detect and measure characteristics of diffusely distributed mobile fluid in tissue spaces that are too small to be imaged directly. (E-mail: [email protected]) © 2006 World Federation for Ultrasound in Medicine & Biology. Key Words: Poroelastography, Elastography, Ultrasound, Poroelastic, Biphasic, Strain, Stress relaxation, Finite element, Porosity, Permeability, Young’s modulus, Poisson’s ratio, Edema, Cancer, Cartilage.

INTRODUCTION

al. 2006; Levick 2003; Ottensmeyer et al. 2004). Furthermore, certain pathologies such as cancer (Jain 1994; Kuszyk et al. 2001) and lymphoedema (Levick 2003) are believed to alter the local fluid flow characteristics of the affected tissue. Thus, imaging the mechanical behavior of soft biological tissue, and relating these images to the local fluid flow through the tissue, may be of diagnostic and prognostic value. The mechanical behavior of soft tissue has previously been modeled by poroelastic and biphasic theories, both of which acknowledge the presence of mobile fluid within the tissue. Although the term was not originally

Soft biological tissue is porous and contains mobile fluid distributed between various fluid compartments such as the vasculature, the interstitium and the lymphatic system (Levick 2003). The mechanical behavior of such tissue can be influenced by the compression-induced flow of this mobile fluid (Elmore et al. 1963; Kerdok et

Address correspondence to: Gearóid 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] 1869

1870

Ultrasound in Medicine and Biology

used by Biot, the concept of a poroelastic material was developed to describe the consolidation of moist soils under an applied load (Biot 1941). A poroelastic material consists of a porous elastic solid matrix and a poreresident fluid that is free to move through the pore space in response to an applied pressure gradient such as might be created by compressing the material. The concept of a biphasic material was developed in a biomechanical context independently of poroelastic theory. A biphasic material consists of a mixture of two mechanically interacting phases, one solid and one liquid, and several biphasic models with varying degrees of complexity have been proposed (Armstrong et al. 1984; Holmes and Mow 1990; Huang et al. 2001; Kwan et al. 1990; Mak 1986; Mak et al. 1987; Mow et al. 1980, 1990; Soltz et al. 2000). One of the simplest of these, which has since become known as the KLM biphasic model (Armstrong et al. 1984; Mow et al. 1980, 1990), assumes that the biphasic material is composed of a linear elastic, isotropic, incompressible solid phase saturated with an incompressible fluid phase. The Biot theory of poroelasticity and the KLM biphasic theory are equivalent when used to model the behavior of saturated poroelastic materials composed of incompressible solid and fluid constituents (Simon et al. 1992). By invoking biphasic theories, previous investigators predicted that tissues rich in mobile fluid would respond to a sustained compression in a time-dependent fashion. For a limited number of geometries and boundary conditions, solutions to the governing equations were obtained that related the underlying material constants of a compressed tissue sample to the time-evolution of either the reaction force on the compressor or the average (global) strain induced inside the tissue (Armstrong et al. 1984; Mow et al. 1980, 1990). The time-evolution of these quantities was also determined experimentally in samples of articular cartilage and, by combining experimental measurements with theoretical predictions, average values of the material constants of the tissue samples were extracted by inverse techniques (Ateshian et al. 1997; Brown and Singerman 1986; DiSilvestro et al. 2001a, 2001b; Mow et al. 1980, 1989; Soltz et al. 1998). Historically, experimental investigations of the strain response of tissue samples to an applied compression were conducted without the benefit of medical imaging technologies and were thus limited to measuring average strain behavior inside the sample, which could be determined by external observation. However, the relatively recent development of acoustic elastography (Ophir et al. 1991) has provided a technique for imaging the local strain induced inside a material in response to an applied compression. When applied to a poroelastic material, the compression can be sustained for a period of time to create stress-relaxation loading conditions and,

Volume 32, Number 12, 2006

Fig. 1. (a) A schematic diagram of the cylindrical coordinate system and a cylindrical sample of a poroelastic material of radius a. The axial direction is defined to be along the z-axis, the radial direction is defined to be along the r-axis and the circumferential direction is along the direction of increasing ␪. (b) The sample was compressed axially between two impermeable compressor plates. Perfectly adhesive boundary conditions were assumed to exist on the contact surfaces between the compressor plates and the sample.

as suggested by Konofagou et al. (1998), if the sample is ultrasonically scanned during this time period, the timeevolution of the local strain inside the sample during stress relaxation might be imaged. To assess whether elastographic techniques might be capable of imaging this behavior, it would be advantageous to possess accurate predictions of the expected local strain response during stress relaxation. However, few predictions regarding the nature of this behavior appear in the literature. Before describing such previous work, 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. Although not

Experimental poroelastic strain imaging ● G. P. BERRY et al.

specifically referring to stress relaxation, Biot (1941) showed that, under the assumptions of phase incompressibility and sample saturation, the volumetric strain inside a poroelastic sample is governed by a simple conduction or diffusion equation. Because conduction or diffusion processes are, in general, time dependent and spatially varying, the compression-induced volumetric strain is therefore predicted to behave in a similar manner. Mow et al. (1980, 1990) considered the axial compression and subsequent stress relaxation of a confined KLM biphasic material with a porous compressor. Although an analytical solution describing the spatial variation of the timeevolution of the axial strain was not reported, a schematic representation of the predicted strain distribution inside the confined sample as a function of time was presented (see Fig. 6b in Mow et al. 1980). The strain behavior was predicted to be both time dependent and spatially varying, and the relaxation process, which exhibited diffusion-like behavior, was forecast to begin first at the permeable boundary and subsequently in regions progressively deeper inside the tissue. The same compression problem was considered by Soltz et al. (1998) and, on this occasion, an analytical solution was provided for the solid matrix axial displacement inside the compressed sample during stress relaxation. However, we found this expression to be incorrect and we present the correct solution in Appendix A, accompanied by the corresponding expression for the time-evolution of the axial strain inside the sample. Since the date of submission of the present manuscript, and during the process of its revision, we have become aware that Soltz et al. (2006) have published a corrigendum that demonstrates agreement with the solution that we independently obtained for the axial displacement. The development of elastography prompted further theoretical investigation of the expected local strain behavior in compressed poroelastic materials. Konofagou et al. (2001) used finite element modeling (FEM) to predict the behavior of the induced ratio of the radial strain to the axial strain (“radial-to-axial strain ratio” or simply “strain ratio”) during the sustained unconfined compression of homogeneous cylindrical samples of poroelastic materials under stress relaxation loading conditions. A time-dependent spatially varying strain field was predicted inside the samples such that the strain ratio was found to decay with time toward an equilibrium value equal in magnitude to the Poisson’s ratio of the solid matrix. This diffusion-like relaxation process was predicted to begin first at the permeable curved cylindrical surface of each simulated sample and subsequently in regions progressively closer to the axis of cylindrical symmetry. Elastographic techniques have already been applied to image the time-evolution of the local strain in com-

1871

pressed cylindrical samples of tofu undergoing unconfined stress relaxation (Righetti et al. 2003, 2004). However, although not noted by the authors, the imaged spatial distribution of the strain ratio at a given instant after the compression had been applied may be seen to have departed from that predicted by earlier finite element simulation (Konofagou et al. 2001). In particular, the time-evolution of the spatial distribution of the strain ratio did not demonstrate the predicted advance with time of lower strain ratio values from the curved cylindrical surface inwards toward the axis of symmetry. In some cases, the reverse effect is evident in the images reported (Righetti et al. 2003, 2004). However, Righetti et al. (2004) cautioned that strain ratio values obtained in the central region of the scan plane and in the regions close to the edge of the scan plane might be unreliable. Further experiments found no statistically significant spatial dependence of the strain response (Righetti et al. 2005). With the awareness, as described before, of the absence of direct experimental evidence confirming the predicted diffusion-like behavior of the strain field in poroelastic materials undergoing sustained compression, we decided to investigate the problem both theoretically and experimentally, to attempt to reconcile theory and experiment. The present paper concentrates on the experimental component of the study, but before describing this it is important to review briefly the outcomes of the theoretical investigation, which were published in Berry et al. (2004a, 2004b, 2006). Because the finite element method is essentially an approximation technique, an analytical approach was pursued. We derived new analytical solutions, which predicted the behavior of the strain field inside homogeneous cylindrical samples of poroelastic materials undergoing stress relaxation in unconfined compression, in the absence of contact friction between the compressor plates and the sample. The poroelastic samples were modeled as KLM biphasic materials. Our predictions supported some, although not all, of those made by Konofagou et al. (2001). Like Konofagou et al. (2001), our analysis predicted diffusion-like strain behavior such that, when the strain ratio images were viewed in a time sequence, lower strain ratio values appeared to advance inwards with time from the curved cylindrical surface toward the axis of symmetry. However, we found different behavior at the curved cylindrical surface than that described by Konofagou et al. (2001). Our work also made several other advances. Firstly, new structure in the strain response not previously anticipated was predicted. Secondly, the analytical model was used to test the numerical accuracy of the strain and pore pressure predictions of commercial finite element software (ABAQUS, Abaqus Inc., Providence, RI, USA) and agreement was found when a sufficiently

1872

Ultrasound in Medicine and Biology

dense finite element mesh and a suitably small time-step were used. Thirdly, by considering the behavior of the pore pressure inside the sample, we described the physical mechanism responsible for the predicted strain response. Fourthly, we showed that the predictions of the KLM biphasic model had previously been incorrectly adapted (Righetti et al. 2004) in an attempt to demonstrate agreement with experimental observations of local strain behavior inside a compressed poroelastic sample. Finally, we demonstrated with simulated data that a knowledge of the time-evolution of the strain field may be used to obtain two quantitative parametric images related to the three intrinsic material constants of the poroelastic material: these are an image of the Poisson’s ratio of the solid matrix (␯s) and an image of the product of the aggregate modulus of the solid matrix (HA) and the permeability of the solid matrix to the pore fluid (k). The reader is referred to Berry et al. (2006) and to Appendix B for a mathematical definition of these material constants. Because our analytical investigation of the problem, as described above, predicted diffusion-like strain behavior similar (but not identical) to that predicted by the finite element study of Konofagou et al. (2001), there remained a pressing need to seek the missing evidence from controlled experiments with poroelastic phantoms to determine whether such behavior could be observed. The purpose of this paper is, therefore, to present the results of such experiments. Using elastographic techniques, the strain was imaged as a function of time inside cylindrical samples of a poroelastic material (tofu) undergoing stress relaxation. To assist with the interpretation of the strain images, finite element simulation was used to extend the predictions of a previous analytical and finite element investigation (Berry et al. 2006) by exploring the influence of experimentally realistic boundary conditions on the strain response. In addition, the model-based parametric image reconstruction algorithm mentioned above was used to produce two parametric images, each related to one or more of the three material constants of one of the samples. This paper represents an elaboration on preliminary experimental results presented earlier, which clearly demonstrated strain-field behavior broadly consistent with theoretical predictions (Berry et al. 2004a, 2004b, 2005a, 2005b). MATERIALS AND METHODS Previous theoretical predictions for strain field behavior in compressed poroelastic samples (Berry et al. 2006; Konofagou et al. 2001) and a previously described parametric imaging technique (Berry et al. 2006) assumed the absence of contact friction between the sample and the compressor plates. In this current work,

Volume 32, Number 12, 2006

efforts were made to minimize this contact friction experimentally to facilitate a comparison with previous theoretical predictions and to increase the applicability of the parametric imaging technique. Nevertheless, frictionless conditions could not be experimentally realized and, thus, finite element simulation was used to investigate the expected influence of contact friction on the strain field behavior. As described below, the effect of perfectly adhesive contact friction was simulated because it was expected to represent the greatest departure from the frictionless case. Nevertheless, it was anticipated that the experimentally imaged strain field behavior would lie somewhere between that predicted by previous theoretical work (Berry et al. 2006) and the simulation results described herein. In this section, the finite element model is first described, followed by a description of the physical experiments. Finally, a preliminary application of the previously developed parametric imaging technique is described.

Finite element simulation Commercial finite element modeling software (ABAQUS) was used to predict the time-dependent behavior of the strain field inside a cylindrical sample of a poroelastic material undergoing unconfined stress relaxation compression between two impermeable parallel planar plates, as shown in Fig. 1b, where the interfaces between the sample and plates were perfectly adhesive. A 1 % compressive axial strain was instantaneously applied and sustained thereafter while preventing radial and circumferential expansion of the surface of the sample in contact with the compressor plates. The sample was assumed to obey the assumptions of the KLM biphasic theory described previously. The permeability was defined to be a constant and was thus strain independent. The material constants of the sample were chosen to match reported values for tofu (Righetti et al. 2004) and were given by [HA,␯s, k] ⫽ [1.74 kPa, 0.3, 1.373 ⫻ 10–11 m4 N–1 s–1]. The sample radius and height were 25 mm and 50 mm, respectively. Axisymmetric finite elements (CAX4RP) were used to predict the behavior on the solution plane shown in Fig. 1a. The mesh was chosen to be sufficiently dense and the time steps were chosen to be suitably small to allow accurate predictions to be obtained within a reasonable time. A uniform mesh density was used consisting of 200 elements in the radial direction and 400 elements in the axial direction. An automatic time-stepping regime was employed that allowed the time-step to vary between 1 s and 10000 s over the course of the relaxation. Fluid flow out of the sample was only permitted through the righthand side of the model, which corresponded to the location of the curved cylindrical surface. This compression

Experimental poroelastic strain imaging ● G. P. BERRY et al.

1873

problem is identical to that described previously (Berry et al. 2006), except for the imposition of perfectly adhesive boundary conditions on the upper and lower surfaces of the sample. To facilitate a comparison with experimental data described below, before displaying the finite element predictions, the predictions obtained for the solution plane were reflected through the axis of symmetry and copied onto the opposite side, so that the solution across the entire radial extent of the sample could be displayed. Fig. 2. The experimental compression apparatus.

Experiments The objective of the experiments was to determine whether elastographic techniques could be used to observe strain field behavior qualitatively consistent with the predictions of poroelastic theory in a variety of poroelastic phantoms subjected to sustained unconfined compression between two impermeable compressor plates. Four cylindrical soyabean gel (tofu) phantoms of radius 25 mm to 28 mm were produced from different types of commercially available tofu (Wing Yip, Croydon, Surrey, UK). Four different samples (samples 1 to 4) were used to show that the general pattern of strain behavior observed was not specific to a single type of tofu phantom. The mechanical properties of each sample were not independently measured in this study and these properties may have differed significantly between the samples. Note that, in this paper, as described in Fig. 1a and b, the axial direction is defined to coincide with both the axis of symmetry of the cylinder and the axis of compression. Although this naming convention is consistent with that used by Fortin et al. (2003), it differs from much previous work in acoustic elastography, where the axial direction is often taken to be along the direction of sound propagation. Our naming convention was chosen because, unlike much previous work on elastography, the compressor and the ultrasound transducer were decoupled in this investigation and were oriented perpendicularly to each other, as described below. As seen in Fig. 2, the axial compression was applied using an inverted elevating table (Time and Precision Industries Ltd., Basingstoke, Hampshire, UK) that could be raised or lowered by a stepper motor under computer control. Attached securely to the elevating table was a block of Perspex that formed the upper compressor plate. The samples were rested on a sheet of Perspex, which formed the lower compressor plate. The faces of the compressor plates were either polished (in the case of sample 1) or were Tefloncoated and lubricated with corn oil (for samples 2 through 4) to reduce friction between the compressor

faces and the samples. The experiments were carried out in a Perspex tank of degassed water. A commercial ultrasound scanner (128XP, Acuson, Mountain View, CA, USA) and an L7 linear array transducer operating at a centre frequency of 5 MHz were used to scan the tofu samples before and periodically after the application of the compression. A Mylar membrane on one side of the tank formed an acoustic window, allowing ultrasound B-mode imaging of the tofu samples across their radial dimension. Thus, the direction of sound propagation in this study was defined to be along the radial direction. The ultrasound scanner had been modified to make the frame sync, line sync and postbeam-formed, postdepth-gain-compensated, analogue intermediate frequency (IF) echo signal data (center frequency 2 MHz) available in real time. The IF signal was digitized at 20 MHz and 10 bits directly to personal computer (PC) random access memory (RAM) by an MV-1000 (Mutech Corporation, Billerica, MA, USA) variable scan frame grabber. A typical B-mode image of one of the tofu samples is shown in Fig. 3. A constant-rate ramped compression was applied over an interval of 2 s to 4 s, inducing a global axial strain in the samples of approximately 0.75 % to 1.5 %. This applied strain was sustained thereafter, producing stress relaxation loading conditions. The duration of each compression was varied from sample to sample and between compressions of the same sample. For sample 1 (the first sample investigated), computer hardware limited the time-sequence length to 220 IF frames acquired at a user-defined interframe interval that was held constant for a given compression of the sample. Because even the order of magnitude of the relaxation time of such tofu samples was initially unknown, sample 1 was compressed three times and, although the interframe time interval was fixed for a given compression of this sample, it was varied between compressions so that the total duration of the 220-frame acquisition could be adjusted. In this way, sample 1 was subjected to successive compressions of a progressively longer duration until the majority of the relaxation process was captured. Between

1874

Ultrasound in Medicine and Biology

Fig. 3. B-mode image of a typical soyabean gel phantom. Image dimensions are approximately 49 mm in the horizontal (radial) direction ⫻ 38 mm in the vertical (axial) direction. The compressor plates are not visible because all phantoms used were sufficiently tall to fill completely the 38 mm lateral extent of the field-of-view of the L7 transducer, and ultrasound imaging was carried out in the region approximately midway between both plates.

successive compressions, the sample was allowed to reabsorb exuded fluid by withdrawing the compressor and allowing the sample to equilibrate in the water tank for at least 10 h. Compressions of samples 2 to 4 were carried out after a hardware upgrade enabling the acquisition of as many as 880 IF frames at an interframe interval that could be varied during the compression of a given sam-

Volume 32, Number 12, 2006

ple. During the stress-relaxation of a poroelastic sample, the rate of change of the average strain inside the sample is higher at early times than at late times (Berry et al. 2006). To capture this early-time behavior with sufficient temporal resolution while still efficiently recording the remainder of the time-dependent variation of the strain field, for samples 2 to 4, a high frame rate (6.5 frames/s) was used initially and the interframe interval was allowed to increase gradually as the relaxation proceeded. A variable interframe interval was achieved by running the MV-1000 in external trigger mode, and IF frames were acquired on command when the card was pulsed with an externally produced voltage via a DIO PCI236 card (Amplicon Ltd., Brighton, West Sussex, UK). Synchronization between frame acquisition and the application of the compression was achieved by pulsing the MV-1000 and driving the stepper motor of the compressor with the same PCI236 card hosted on the same PC. Application software written in C⫹⫹ (Borland Corporation, Cupertino, CA, USA) controlled the operation of the PCI236. Using a cubic spline interpolation scheme, the IF data were upsampled to 80 MHz and the number of A-lines was increased by a factor of four. Radiofrequency (RF) data were recovered from the interpolated IF data by mixing with an appropriately-scaled local oscillator frequency and low-pass filtering the mixed output (Doyley et al. 2001). Two-dimensional cross-

Fig. 4. FEM predictions, in the case of perfectly adhesive boundary conditions, for the time evolution of the radial strain, axial strain, absolute value of the radial-to-axial strain ratio and volumetric strain across the entire axial and radial dimensions of a plane through the axis of symmetry of the cylindrical sample. Images in the top row of the radial strain and strain ratio columns have a different color scale to the other images in those columns. Note that thresholds were applied to the color scales in all columns to prevent the behavior at the corners of the images from reducing the contrast in the central region of the images (the predicted behavior of the central region is of greater interest for the purposes of comparison with the available experimental data). The center line of the images in the axial direction corresponds to the axis of symmetry of the simulated sample.

Experimental poroelastic strain imaging ● G. P. BERRY et al.

1875

Fig. 5. FEM predictions analogous to those shown in Fig. 4, but in the absence of contact friction. The images were obtained by reformatting and rescaling data already presented in Berry et al. (2006).

correlation RF echo-tracking was performed between the first pre-compression frame in a time sequence and all subsequent frames in that sequence, and parabolic interpolation was used to find the subpixel peak of the correlation function. In this way, the resulting time sequence of displacement images displayed the timedependent displacement of acoustic scatterers in the samples relative to their original positions before the compression was applied. Radial and axial strain estimates were subsequently obtained by median-filtering the displacement images and applying a leastsquares strain estimator (LSQSE) (Kallel and Ophir 1997) to the filtered displacement images in each direction independently. Each strain ratio image was obtained by dividing each radial strain image by the corresponding axial strain image on a pixel-by-pixel basis. All times were measured relative to the instant at which the compressor reached its final position, i.e., at the end of the ramp. To cope with the computational demands of displacement tracking and strain estimation across a large number of RF frames, calculations were carried out in parallel on a load-balanced sixnode PC cluster. Parametric imaging For one of the tofu samples (sample 4), a parametric imaging technique, described in Berry et al. (2006), was used. Under frictionless boundary conditions, the spatiotemporal evolution of the absolute value of the radial-toaxial strain ratio inside a homogeneous cylindrical poroelastic sample subjected to unconfined stress relaxation compression is given as follows:

冏

冏

␧rr r, t ⫽ vs ␧zz 共 兲

J0共␣n兲 ⫺ ⬁

⫹

兺

n⫽1

冉 冊 兲 兲冦 冉 冊

共1 ⫺ 2vs 共1 ⫺ vs

J0

␣nr ⫺ a

J1

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

1

冧

J ␣ ⫺␣ J ␣ 共1 ⫺ vs兲 0共 n兲 n 1共 n兲 ⫻ exp(

⫺␣n2HAk t) a2

(1)

where ␧rr and ␧zz are the normal strains in the radial and axial directions, respectively, as a function of time t and radial position r, a is the sample radius, J0(x) and J1(x) are Bessel functions of the first kind of order zero and order one, respectively, and ␣n are the roots of the corresponding characteristic equation (Berry et al. 2006). Two parametric images related to the three material constants of the sample were obtained by curve-fitting eqn (1) to the time-evolution of the radial-to-axial strain ratio. The two fit parameters [p1, p2] used were given as follows: p1 ⫽ HAk p2 ⫽ vs

(2)

The technique used an implementation of the LevenbergMarquardt algorithm to find a best fit to the experimental data by minimizing the sum of the squared differences by iteration. This fitting routine has previously been shown to deliver a best fit, which, for a particular set of material

1876

Ultrasound in Medicine and Biology

Volume 32, Number 12, 2006

constants and known sample radius, is accurate, precise, relatively tolerant to noise and unique (Berry et al. 2006). RESULTS Finite element simulation Figure 4 shows the predicted strain (radial, axial and volumetric) and strain ratio distributions inside the simulated sample as a function of time. All quantities in Fig. 4 are time dependent and spatially varying, and equilibrium is eventually reached after approximately 60000 s. It is clear that the mechanical behavior differs from that predicted for the frictionless case shown in Fig. 5 (these predictions are repeated from Berry et al. (2006) after reformatting to facilitate a comparison with present results). Nevertheless, in both the perfectly adhesive and frictionless cases, immediately after the application of the compression (at time t ⫽ 0⫹ s), there is no volumetric strain inside the sample, indicating that the sample is incompressible at this time. Experimental strain images For each sample, 3-D datasets (each consisting of a time-sequence of 2-D strain images) were obtained. Because the LSQSE strain imaging technique necessarily results in image cropping along the direction in which it is applied, neither the axial strain near the compressor plates nor radial strain near the curved cylindrical boundary could be imaged. Furthermore, strain ratio imaging could only be carried out in the central region of the samples where the axial and radial strain images overlapped. In this central region, the imaged quantities were found to show little dependence on the axial coordinate. Thus, the 3-D datasets could be reduced to 2-D by column-wise averaging the imaged quantity in the axial direction at each radial position and, in this way, the time evolution of the radial dependence of the imaged quantity over a large number of RF frames, each acquired at

Fig. 6. The column-averaged time evolution of (a) radial strain, (b) axial strain and (c) absolute value of the radial-to-axial strain ratio for three successive compressions of sample 1. Each image consists of 219 rows, and each row can be considered to be equivalent to a single image at a given instant shown in either Figs. 4 or 5. For each row, the time since the start of the compression is as indicated. Positive strain values indicate expansion and negative strain values indicate compression. Displacement tracking was carried out using a reference kernel with dimensions 2.1 mm (radial) ⫻ 1.5 mm (axial) and a 70% overlap between adjacent displacement estimates. The length of the LSQSE estimators used in the radial and axial directions were 20 and 40 points, respectively. The dimensions of the strain images before column-wise averaging were approximately 44 mm (radial) ⫻ 14 mm (axial).

Experimental poroelastic strain imaging ● G. P. BERRY et al.

1877

Fig. 7. The column-averaged time evolution of radial strain, axial strain and absolute value of the radial-to-axial strain ratio for a single compression of sample 2. Each image consists of 36 rows. Each row was obtained from a single strain image selected from the full time sequence of available strain images such that the inter-row interval was 100 s. The signal processing parameters used were the same as for Fig. 6. The dimensions of the strain images before column-wise averaging were approximately 40 mm (radial) ⫻ 14 mm (axial).

a different time, could be concisely displayed. Figures 6 to 8 show the time-dependent strain behavior inside samples 1 to 3 summarized in this 2-D form, where each row in each of the images may be considered to be equivalent to a single image at a given time in the time sequences displayed in Figs. 4 and 5.

To produce parametric images for sample 4 using eqn (1), strain ratio data at each point on the imaged plane as a function of time were required. Thus, Figs. 9 and 10 show a time sequence of radial and axial displacements, and strain and strain ratio images on the imaged plane through the axis of symmetry of sample 4,

Fig. 8. The column-averaged time evolution of radial strain, axial strain and absolute value of the radial-to-axial strain ratio for a single compression of sample 3. Each image consists of 36 rows where the inter-row interval was 100 s. The signal processing parameters used were the same as for Fig. 6. The dimensions of the strain images before column-wise averaging were approximately 44 mm (radial) ⫻ 14 mm (axial).

1878

Ultrasound in Medicine and Biology

Volume 32, Number 12, 2006

because of the reduced contact friction that resulted from the use of Teflon and lubricant for these other samples. Parametric imaging Figure 12a displays the time dependence of the strain ratio at two internal positions in sample 4, where one position (pixel 1) was located on the axis of symmetry and the other (pixel 2) was situated at the rightmost edge of the strain ratio images. The curve of best fit to the data is also illustrated. Furthermore, the early time behavior of the strain ratio at these two positions is shown in Fig. 12b, demonstrating evidence of an early time plateau for the pixel located on the axis of symmetry as predicted previously (Berry et al. 2006). Parametric images of HAk and ␯s in the central region of sample 4 are shown in Fig. 13a and b, respectively. The mean and standard deviation of the pixel values in these images were: HAk ⫽ 1.449 (⫾ 0.269) ⫻ 10–7 m2 s–1 and ␯s ⫽ 0.088 (⫾ 0.023). Figure 13c displays the corresponding minimized cost–function matrix, where smaller costfunction values indicate better agreement between theory and experiment. DISCUSSION Fig. 9. The time evolution of radial and axial displacement throughout the imaged plane for a single compression of sample 4. Radial displacements to the right and left are represented by positive and negative values of displacement, respectively. Downward axial displacements are represented by positive values of displacement. Image dimensions are approximately 49 mm (radial) ⫻ 30 mm (axial). The dimensions of the displacement tracking reference kernel used were 4.1 mm (radial) ⫻ 1.5 mm (axial), and a 70% overlap between adjacent displacement estimates was applied. A 5 ⫻ 5 point median filter was used to remove outliers.

but only for a limited number of available times during the sustained compression. In Fig. 11, to demonstrate the effect of contact friction on the induced strain field, the radial strain distribution inside sample 1 is displayed over a large field of view for a limited number of available times. Because radial strain can be imaged without the application of the LSQSE in the axial direction, a larger field of view can be achieved in the axial direction than is possible with either axial strain or strain ratio images, thus enabling strain imaging in regions very close to the compressor plates, where the influence of contact friction might be expected to be most obvious. It is seen that the radial strain exhibits a dependence on the axial coordinate. Although a similar effect was seen for the other samples (not shown), the dependence on the axial coordinate was less pronounced than for sample 1, possibly

Finite element simulation The behavior of the strain field under perfectly adhesive boundary conditions is clearly predicted to differ from the frictionless case investigated previously (Berry et al. 2006). Of particular relevance to the experimental component of this study are three classes of behavior that are predicted to occur in the presence of contact friction but not in its absence. First, both immediately after the compression is applied and at equilibrium, regional variations in the radial and axial strains and the strain ratio distributions are evident. Second, the local axial strain is predicted to be both time dependent and spatially varying such that, as the relaxation proceeds, axial strain in central regions of the sample is transferred to regions nearer the compressor plates. Furthermore, both the radial strain and the strain ratio exhibit a dependence on the axial coordinate. Despite the differences between the frictionless and perfectly adhesive cases, radial profiles taken across the central region of the sample, far from both adhesive compressor plates, exhibit behavior that qualitatively resembles that of the frictionless case, insofar as lower values of radial strain and strain ratio propagate inwards toward the axis of symmetry and the dependence on the axial coordinate is weaker than in regions closer to the compressor plates. Experimental strain imaging The behavior of the volumetric strain is of fundamental importance in the study of the mechanical behav-

Experimental poroelastic strain imaging ● G. P. BERRY et al.

1879

Fig. 10. The time evolution of radial strain, axial strain and absolute value of the radial-to-axial strain ratio throughout the scan plane for a single compression of sample 4. LSQSE estimators used were 10 and 20 points in length in the radial and axial directions, respectively. The radial and axial strain images were cropped to be equal in size to the strain ratio image, with approximate dimensions of 37 mm (radial) ⫻ 20 mm (axial).

ior of poroelastic materials (Berry et al. 2006) and, thus, it would be desirable to image this quantity directly. Its fundamental importance is well illustrated in Fig. 4, where the exotic behavior of the radial strain, axial strain and strain ratio immediately after the compression is applied can be seen to contribute to a uniform volumetric strain distribution, which can be easily interpreted to show that no fluid flow out of the sample occurs at this time. However, volumetric strain could not be imaged in this study because of a lack of a 3-D strain imaging capability. Nevertheless, using 2-D ultrasound imaging and acoustic elastographic techniques, solid matrix displacements and strains in the radial and axial directions and the radial-to-axial strain ratio were imaged on a plane through the axis of symmetry of four cylindrical samples of tofu undergoing unconfined stress-relaxation compression. Acoustic elastographic techniques can produce high precision strain estimates along the direction of sound propagation, but displacement and, hence, strain estimates on the scan plane in the direction perpendicular to the acoustic axis are inferior (Hein and O’Brien 1993; Konofagou et al. 1998; Lubinski et al. 1996). Thus,

Fig. 11. The time evolution of the distribution of radial strain on the imaged plane during the third compression of sample 1. Image dimensions are approximately 44 mm (radial) ⫻ 34 mm (axial).

1880

Ultrasound in Medicine and Biology

Volume 32, Number 12, 2006

symmetry. A similar pattern was evident in the timedependent behavior of the strain ratio. The specific nature of this diffusion-like strain response to the applied axial compression arises from the presence of pore fluid inside the sample and its compression-induced flow (Berry et al. 2006). Axial compression causes pore fluid pressurization and radial expansion of the sample beyond the sample’s equilibrium radial dimensions. Subsequent fluid exudation out of the sample facilitates the radial contraction of the solid matrix toward its equilibrium dimensions. The response is spatially varying because fluid flow does not simultaneously begin everywhere inside the sample: it begins first at the curved cylindrical surface and, subsequently, in regions progressively deeper inside the sample. The delay experienced by internal regions in beginning the journey toward equilibrium can manifest itself in the formation of an early-time strain ratio plateau, as seen in Fig. 12b. Qualitatively similar behavior was observed across all samples, indicating that the general pattern of strain behavior was not specific to a single type of tofu phan-

Fig. 12. (a) The time-dependent strain ratio at two internal locations and corresponding curve fits. Both points are at the same height in the sample but one point is on the axis of symmetry (pixel 1) and the other is nearer to (but not on) the curved cylindrical surface of the phantom (pixel 2). (b) A close-up of the early time strain ratio behavior for the same two points.

because of the transducer orientation used in this study, the strain signal-to-noise ratio was greater in the radial direction of the coordinate system described in Fig. 1a than in the axial direction. It can be seen in Figs. 6 to 10 that the application of the compression caused samples to compress axially and, initially, to expand in the radial direction. The initial radial expansion of the samples was followed by a timedependent spatially varying radial contraction (relative to the sample dimensions immediately after the compression was applied). During this contraction phase, lower values of radial strain were seen to propagate inwards from the curved cylindrical surface toward the axis of

Fig. 13. Parametric images of (a) HAk, (b) ␯s and (c) the corresponding cost–function matrix in the central region of sample 4. Image dimensions are approximately 37 mm (radial) ⫻ 20 mm (axial).

Experimental poroelastic strain imaging ● G. P. BERRY et al.

tom. It is unclear why similar behavior was not detected in previous elastographic investigations of compressed cylindrical tofu phantoms (Righetti et al. 2003, 2004, 2005). However, the transducer in these previous studies was oriented such that the direction of sound propagation coincided with the axis of compression, and we note that this configuration would have been suboptimal for imaging changes in radial strain associated with radial fluid flow. Although the observed behavior of the strain response was broadly consistent with previous predictions obtained using analytical (Berry et al. 2006) and finite element techniques (Berry et al. 2006; Konofagou et al. 2001), some departures from previous predictions did occur. First, according to the perfect-slip analysis (Berry et al. 2006; Konofagou et al. 2001), the radial-to-axial strain ratio of the poroelastic sample was predicted to be distributed uniformly inside the sample, both immediately after the compression was applied and at equilibrium. However, in Figs. 6, 7, 8 and 10, it can be seen that, although the strain ratio was close to 0.5 immediately after the compression was applied, it was not uniformly distributed throughout the samples at this time. Instead, strain ratio values were initially lower in regions near the curved cylindrical surface than in regions near the axis of symmetry of the sample. Also, at the end of the compression sequences, strain ratio image values were not uniformly distributed throughout the samples. Second, as seen in Figs. 6, 7, 8 and 10, despite the use of a constant (post-ramp) applied strain, a time-dependent reduction in axial strain was observed in the imaged (central) region of all samples. Third, as shown in Fig. 11, the radial strain in a region of sample 1 of sufficient axial extent to include the regions very close to the compressor plates exhibited a dependence on the axial coordinate. Nevertheless, the strains and strain ratio in the central region of the samples showed little axial variation, which justified our decision, in the case of samples 1 to 3, to average the imaged quantity in the axial direction so that the radial variation over a large number of time instants could be concisely shown. As demonstrated herein using finite element modeling, all of these departures from the earlier perfect-slip predictions are consistent with the existence of contact friction between the samples and the compressor plates. Furthermore, the presence of bow-shaped contours in the radial displacement images of Fig. 9 appears to indicate that regions near the compressor plates experienced some degree of radial confinement and, thus, appear to confirm the existence of contact friction. In addition to the absence of contact friction, previous predictions (Berry et al. 2006) were made under several other simplifying assumptions. For example, the poroelastic sample was previously modeled as a linear

1881

elastic isotropic porous solid matrix composed of an incompressible solid phase and saturated with an incompressible pore fluid. In addition, the applied axial compression was assumed to be infinitesimal and instantaneously applied. It is important to note that failure to satisfy these other assumptions experimentally could also have contributed to the departure of the observed strain behavior from previous predictions. Indeed, these assumptions may only have been satisfied experimentally to a limited extent for several reasons. First, the mechanical properties of the tofu samples were largely unknown. The solid matrix may have been inherently viscoelastic, nonlinear elastic, orthotropic or anisotropic because we did not exclude the possible existence of such behavior by mechanical testing. Furthermore, although water is known to be incompressible at pressures encountered in this experiment, it is not clear whether the solid material from which the solid matrix was constructed was incompressible. Second, in an attempt to select only samples with a high degree of pore saturation, B-mode imaging was used to reject candidate tofu samples containing large gas bubbles. However, the presence of smaller gas bubbles may have gone undetected. Third, to achieve an acceptable strain signal-to-noise ratio in the radial and axial directions, a small applied axial strain of approximately 0.75 % to 1.5 % was used, but this strain may have been sufficiently large to invalidate significantly the infinitesimal strain assumption. Finally, although a good approximation to an instantaneously applied strain was achieved by using a ramp time that was short (2 to 4 s) relative to the relaxation time of the tofu (ⱖ1 h), some relaxation may have occurred before the first post-ramp image was acquired. Parametric imaging The local time-dependent strain inside a sample of a poroelastic material undergoing stress relaxation is influenced by three material constants (HA,␯s, k) and the characteristic length of the sample (Berry et al. 2006). Note that the Young’s modulus of the solid matrix may be used instead of HA as one of the material constants (Berry et al. 2006). Thus, an analysis of the compressioninduced time-dependent strain response could assist in material characterization. A model-based poroelastic reconstruction technique, which shares the same assumptions as the previous frictionless analysis (Berry et al. 2006), has recently been developed for this purpose, enabling the production of parametric images of ␯s and the product of the other two material constants HAk (Berry et al. 2006). It was not possible to separate uniquely the two components of the HAk parametric image using the available strain data alone. To illustrate the applicability of the technique, the same model-based reconstruction algorithm was applied to the time se-

1882

Ultrasound in Medicine and Biology

quence of strain ratio data obtained for one of the tofu samples (sample 4). Although contact friction was present during the experiment, parametric imaging was only carried out in the central region of the sample, where the influence of this friction was less pronounced than in regions near the compressor plates. Because only tofu samples that appeared to be homogeneous were chosen for this study, it was anticipated that the two parametric images obtained would also be homogeneous. However, although the standard deviation of both parametric images was relatively small, as seen in Fig. 13, the parametric images did exhibit some regional variation that may not have been random. It is important to recognize that, although any inhomogeneities in the sample may have contributed to this regional variation, the images themselves should not be interpreted as displaying the spatial variation of the parametric values. This is because of the a priori assumption of poroelastic homogeneity made by the model. Thus, like the nonimaging model-based methods described in the introduction, the important quantities delivered by this model-based poroelastic imaging technique are the mean parametric values in the material, obtained by calculating the mean of the pixel values in the parametric images. However, unlike earlier methods, the use of local strain data means that, in principle, this technique may eventually be extended to obtain meaningful local estimates of these parametric values in layered and other heterogeneous poroelastic media, if it is known how the strain field at a particular location and at a given time is influenced by the spatial distribution of material constant values. This approach is analogous to the production of Young’s modulus images in elastic materials from knowledge of the local strain field and the associated boundary conditions (e.g., Barbone and Bamber 2002; Doyley et al. 2000; Kallel and Betrand 1996; Oberai et al. 2004). However, as described below, even in the absence of a generalized model-based reconstruction technique, the strain data alone are a rich source of information about the mechanical properties of the material. To validate quantitatively the accuracy of this material-constant reconstruction method, independent techniques for accurately and precisely determining the values of the material constants are required for the purposes of comparison. However, it has been suggested that existing soil mechanics techniques for measuring permeability exhibit experimental errors as large as one order of magnitude (Bowles 1988). Furthermore, it is important to note that the required elastic moduli for any such comparison are those of the solid matrix and not those of the solid material from which the matrix is constructed. Such elastic moduli can be determined either when the pore fluid has been removed from the sample or when the sample has reached equilibrium after an applied compression (the sample may remain saturated in this case). However, in the former case, remov-

Volume 32, Number 12, 2006

ing all of the pore fluid without damaging or changing the properties of the solid matrix would be a challenging task and, in the latter, some uncertainty would exist with regard to exactly when equilibrium had been reached. Nevertheless, Righetti et al. (2004, 2005) attempted to perform mechanical measurements of these material constants in several types of tofu and, when the mean pixel values of the parametric images in Fig. 13 are compared with these independent measurements, agreement is obtained to within the large uncertainty of those independent measurements. Significance of the observed diffusion-like strain behavior To our knowledge, the diffusion-like strain behavior observed in all tofu samples in this study represents the first experimental detection of this phenomenon in any poroelastic material by any imaging modality and is a significant observation for several reasons. First, as described in the introduction, such behavior is a key prediction of poroelastic theory. The experimental demonstration of this effect increases our confidence both that elastographic techniques may be used to image the spatio-temporal behavior of the strain in compressed poroelastic media and that, at least in the case of simple phantom materials such as tofu, the Biot poroelastic theory may be used to predict and interpret the strain field behavior. Whether this is also true in the case of soft biological tissue in vivo is currently under investigation. Second, because such behavior is a characteristic signature of a poroelastic response (Berry et al. 2006; Biot et al. 1941), it reveals that mobile fluid is present in the sample. This is especially significant because the presence of this fluid is not obvious in the B-mode images obtained (Fig. 3), suggesting that strain imaging can detect the presence of mobile fluid even when that fluid is accommodated in tiny interconnected pores, the dimensions of which are below the spatial resolution of the ultrasound system. Thus, in the future, this imaging technique might be useful for oedema diagnosis because oedematous tissues are known to possess relatively large volumes of mobile interstitial fluid (Levick 2003), yet a reinterpretation of the results of previous ultrasound investigations suggests that the fluid-filled pores are too small to be imaged directly (Dent 2000; Doldi 1992; Mellor 2004; van der Veen 2001). The technique might also be applicable to characterizing the mechanism(s) behind the transient behavior known to occur in all soft biological tissue (Fung 1993) because the existence or absence of a diffusion-like strain response might help to implicate or exclude fluid flow as a source of the observed transient behavior. In addition, information about the nature of the compression-induced fluid flow can be obtained by studying the diffusion-like behavior of the strain field. As

Experimental poroelastic strain imaging ● G. P. BERRY et al.

predicted by Konofagou et al. (2001), and now demonstrated experimentally herein, the observed strain behavior can identify the route taken by the fluid as it flows through the sample because fluid flow occurs in the opposite direction to which the diffusion-like relaxation process propagates (Berry et al. 2006). Furthermore, although not demonstrated here, it may eventually be possible to image the pore fluid flow velocity by invoking biphasic or poroelastic theory to relate local displacement and strain to the local fluid velocity (see eqn (68) in Mow et al. (1980) and eqn (14) in Armstrong et al. (1984)). Finally, although the correspondence between analytical and finite elemental predictions with regard to this diffusion-like strain behavior has previously been demonstrated (Berry et al. 2006), the experimental discovery of this effect now provides qualitative experimental validation of the finite element predictions. Thus, in the future, the finite element software may be used with more confidence to predict strain field behavior in poroelastic materials with more general geometries and boundary conditions. CONCLUSIONS For the first time, radial and axial displacements and strains, and the radial-to-axial strain ratio, were imaged in cylindrical samples of a poroelastic material undergoing stress relaxation by scanning across the radial extent of the cylinders. As the compression was sustained, lower radial strain and radial-to-axial strain ratio values were observed to advance inwards from the permeable boundaries at the curved cylindrical surface toward the axis of symmetry. The detection of this diffusion-like phenomenon was repeatable within and between samples and provides the first qualitative validation of the predictions of poroelastic theory with regard to this behavior. The observed behavior was found to depart in some respects from previous predictions (where frictionless boundary conditions were assumed), but finite element modeling was used to show that such behavior was consistent with the existence of contact friction between the sample and the compressor. Finally, it was demonstrated that images of parameters related to the material constants of a poroelastic sample may be derived from knowledge of the experimentally imaged time-dependent strain ratio in the sample during stress relaxation. The mean pixel values in these parametric images were found to agree with values derived from previously published data for tofu samples. The results demonstrate that it may eventually be possible to use strain data to detect and measure characteristics of diffusely distributed mobile fluid in tissue spaces that are too small to be imaged directly.

1883

Acknowledgements—This work was supported by grants from the Institute of Cancer Research and the Engineering and Physical Sciences Research Council, UK. P. E. Barbone 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. Ateshian GA, Warden WH, Kim JJ, Grelsamer RP, Mow VC. Finite deformation biphasic material properties of bovine articular cartilage from confined compression experiments. J Biomechanics 1997;30(11/12):1157–1164. Barbone PE, Bamber JC. Quantitative elasticity images: What can and cannot be inferred from strain images. Phys Med Biol 2002;47: 2147–2164. Berry GP, Bamber JC, Armstrong CG, Miller NR, Barbone PE. Towards an acoustic model-based poroelastic imaging method: I. Theoretical foundation. Ultrasound Med Biol 2006;32(4):547–567. Berry GP, Bamber JC, Miller NR, Bush NL. Ultrasonic strain-based imaging of tissue permeability and mobile fluid volume [abstract]. Ultrasound 2004a;12(4):232. Berry GP, Bamber JC, Armstrong C, Miller NR. Towards a modelbased poroelastic imaging method. Proceedings of the Third International Conference on the ultrasonic measurement and imaging of tissue elasticity [abstract]. October 17–20, Cumbria, UK: 2004b; p 96. Berry GP, Bamber JC, Bush NL, Miller NR. The spatio-temporal variation of the strain field inside compressed poroelastic materials. Proceedings of the Fourth International Conference on the ultrasonic measurement and imaging of tissue elasticity [abstract]. October 16 –19, Austin, TX, USA: 2005a; p 97. Berry GP, Bamber JC, Armstrong C, Miller NR. The behaviour of the strain field inside compressed poroelastic materials [abstract]. Ultrasound 2005b;13(4):261–262. Biot MA. General theory of three-dimensional consolidation. J Appl Phys 1941;12:155–164. Bowles JE. Engineering properties of soils and their measurements. New York: McGraw Hill; 1992. Brown TD, Singerman RJ. Experimental determination of the linear biphasic constitutive coefficients of human fetal proximal femoral chondroepiphysis. J Biomech 1986;19(8):597– 605. Dent CL, Scott MJ, Wickline SA, Hall CS. High frequency ultrasound for quantitative characterization of myocardial edema. Ultrasound Med Biol 2000;26(3):375–384. DiSilvestro MR, Zhu Q, Wong M, Jurvelin JS, Suh JKF. Biphasic poroviscoelastic simulation of the unconfined compression of articular cartilage: I. Simultaneous prediction of reaction force and lateral displacement. J Biomech Eng 2001a;123:191–197. DiSilvestro MR, Zhu Q, Suh JKF. Biphasic poroviscoelastic simulation of the unconfined compression of articular cartilage: II. Effect of variable strain rates. J Biomech Eng 2001;123:198 –200. Doldi SB, Lattuada E, Zappa MA, Pieri G, Favara A, Micheletto G. Ultrasonography of extremity lymphedema. Lymphology 1992;25: 129 –133. Doyley MM, Bamber JC, Fuechsel F, Bush NL. A freehand elastographic approach for clinical breast imaging: System development and performance evaluation. Ultrasound Med Biol 2001;27:1347– 1357. Doyley MM, Meaney PM, Bamber JC. Evaluation of an iterative reconstruction method for quantitative elastography. Phys Med Biol 2000;45:1521–1540. Elmore SM, Sokoloff L, Norris G, Carmeci P. Nature of “imperfect” elasticity of articular cartilage. J Appl Physiol 1963;18(2):393– 396. Fortin M, Buschmann MD, Bertrand MJ, Foster S, Ophir J. Dynamic measurement of internal solid displacement in articular cartilage using ultrasound backscatter. J Biomech 2003;36:443– 447. Fung YC. Biomechanics: Mechanical properties of living tissues, 2nd ed. New York: Springer-Verlag, Inc., 1993.

1884

Ultrasound in Medicine and Biology

Hein IA, O’Brien WD. Current time-domain methods for assessing tissue motion by analysis from reflected ultrasound echoes-a review. IEEE Trans Ultrason Ferroelectr Freq Control 1993;40(2): 84 –102. Holmes MH, Mow VC. The nonlinear characteristics of soft gels and hydrated connective tissues in ultrafiltration. J Biomech 1990; 23(11):1145–1156. Huang CY, Mow VC, Ateshian GA. The role of flow-independent viscoelasticity in the biphasic tensile and compressive responses of articular cartilage. J Biomech Eng 2001;123(5):410 – 417. Jain RK. Barriers to drug delivery in solid tumors. Sci Am 1994; 271(1):58 – 65. Kallel F, Bertrand M. Tissue elasticity reconstruction using linear perturbation method. IEEE Trans Med Imag 1996;15:299 –313. Kallel F, Ophir J. A least-squares strain estimator for elastography. Ultrason Imaging 1997;19(3):195–208. Kerdok AE, Ottensmeyer MP, Howe RD. Effects of perfusion on the viscoelastic characteristics of liver. J Biomech 2006;39(12):2221– 2231. 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. Tumor transport physiology: Implications for imaging and image-guided therapy. Am J Roentgen 2001;177:747–753. Kwan MK, Lai WM, Mow VC. A finite deformation theory for cartilage and other soft hydrated connective tissues: I. Equilibrium results. J Biomech 1990;23(2):145–155. Levick JR. An introduction to cardiovascular physiology, 4th ed. London: Arnold Publishers; 2003. Lubinski MA, Emelianov SY, Raghavan KR, Yagle AE, Skovoroda AR, O’Donnell M. Lateral displacement estimation using tissue incompressibility. IEEE Trans Ultrason Ferroelectr Freq Control 1996;43(2):247–256. Mak AF. The apparent viscoelastic behavior of articular cartilage—The contributions from the intrinsic matrix viscoelasticity and interstitial fluid flows. J Biomech Eng 1986;108(2):123–131. Mak AF, Lai WM, Mow VC. Biphasic indentation of articular cartilage: I. Theoretical analysis. J Biomech 1987;20(7):703–714. Mellor RH, Bush NL, Stanton AW, Bamber JC, Levick JR, Mortimer PS. Dual-frequency ultrasound examination of skin and subcutis thickness in breast cancer-related lymphedema. Breast 2004;10(6): 496 –503. Mow VC, Kuei SC, Lai WM, Armstrong CG. Biphasic creep and stress relaxation of articular cartilage in compression: theory and experiments. J Biomech Eng 1980;102:73– 84. Mow VC, Gibbs MC, Lai WM, Zhu WB, Athanasiou KA. Biphasic indentation of articular cartilage – II. A numerical algorithm and an experimental study. J Biomechanics 1989;22(8/9):853– 861. Mow VC, Ratcliffe A, Woo SLY (eds). Biomechanics of diarthrodial joints, vol 1. New York: Springer Verlag, Inc.; 1990: 215–260. 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. Ottensmeyer MP, Kerdok AE, Howe RD, Dawson SL. The effects of testing environment on the viscoelastic properties of soft tissues (page 9 –18). In Cotin S, Metaxas DN (eds): Lecture Notes in Computer Science. Medical Simulation: International Symposium, ISMS 2004, Cambridge MA, USA, June 17-18, 2004. Proceedings. Berlin/Heidelberg, Germany: Springer, 2004. Righetti R, Ophir J, Srinivasan S, Krouskop T. 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

Volume 32, Number 12, 2006 tissue elasticity [abstract], October 12-15, 2003, Corpus Christi, TX, USA. 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. Righetti R, Ophir J, Krouskop TA. A method for generating permeability elastograms and Poisson’s ratio time-constant elastograms. Ultrasound Med Biol 2005;31(6):803– 816. Simon BR. Multiphase poroelastic finite element models for soft tissue structures. Appl Mech Rev 1992;45(6):191–218. Soltz MA, Ateshian GA. Experimental verification and theoretical prediction of cartilage interstitial fluid pressurization at an impermeable contact interface in confined compression. J Biomech 1998; 31:927–934. Soltz MA, Ateshian GA. Corrigendum to “Experimental verification and theoretical prediction of cartilage interstitial fluid pressurization at an impermeable contact interface in confined compression”: [J Biomech 1998;31:927–934]. J Biomech 2006;39(3):594. Soltz MA, Ateshian GA. A conewise linear elasticity mixture model for the analysis of tension-compression nonlinearity in articular cartilage. J Biomech Eng 2000;122(6):576 –586. van der Veen P, Vermeiren K, Von Kemp K, Lamote J, Sacre R, Lievens P. A key to understanding postoperative lymphoedema: A study on the evolution and consistency of oedema of the arm using ultrasound imaging. Breast 2001;10:225–230.

APPENDIX A We wish to determine the nature of the internal strain field induced in the solid matrix of a confined sample of a poroelastic material of height, h, when such a sample is subjected to a stress relaxation compression by a permeable free-draining compressor plate. The problem is schematically represented in Fig. 14. The axial displacement in the solid matrix of the sample is governed by the following differential equation (Mow et al. 1990; Soltz et al. 1998): ⭸ 2u 1 ⭸u ⫽ ⭸z2 HAk ⭸t

(A1)

where u(z,t) is the axial displacement of the sample at height z and time t. The compression is modeled by applying an axial displacement to the upper surface of the sample while the position of the lower compressor is fixed. Thus, the behavior of the lower surface of the sample is described by the following boundary condition:

Fig. 14. The 1-D confined compression problem considered by Soltz et al. (1998). The compressive strain, ␧0, is applied over a time t0 (and sustained thereafter) by displacing the upper surface of the sample at a constant rate in the axial direction with a porous compressor, and the lower surface is fixed in position. The direction of compression-induced fluid flow is indicated.

Experimental poroelastic strain imaging ● G. P. BERRY et al. u共0, t兲 ⫽ 0

The spatial derivative of eqn (A5) yields the induced axial strain, ␧, as a function of depth and time and is given as follows:

(A2)

A time-dependent ramped axial displacement is applied to the upper surface of the material as follows: u共z ⫽ h, t兲 ⫽

再

⫺V0t, 0 ⬍ t ⬍ t0 ⫺V0t0,

⫽

where t0 is the ramp duration and V0 is the descent velocity of the compressor in the axial direction. At t ⫽ t0, the applied displacement is held constant, enabling a stress relaxation compression scheme. The following initial condition is used: u共z, 0兲 ⫽ 0

(A4)

The solution to eqn (A1) subject the boundary and initial conditions described in eqns (A2)-(A4) is given as follows: u共z, t兲

⫽

冦

⫺

⫺

V0tz 2V0h2 ⫹ 3 h ␲ HAk

V0t0z 2V0h2 ⫹ 3 h ␲ HAk

⬁

兺

n⫽1

⬁

兺

共⫺1兲n

n⫽1

n3

共⫺1兲n n3

共e

共1 ⫺ e

n2 ␲ 2 H A k t0 h2

冉 冊 冉 冊

⫺ 1 sin

n␲ z , h

兲e

sin

⫺n2␲2HAk t h2

兲

␧共z, t兲

(A3)

t0 ⱕ t

⫺n2␲2HAk t h2

1885

冦

⫺

⫺

V0t 2V0h ⫹ 2 h ␲ HAk

V0t0 2V0h ⫹ 2 h ␲ HAk

⬁

兺

n⫽1

⬁

兺

共⫺1兲n n2

n⫽1

共⫺1兲n n2

共e

共1 ⫺ e

n2 ␲ 2 H A k t0 h2

冉 冊 冉 冊

⫺ 1 cos

n␲ z , h

兲e

cos

⫺n2␲2HAk t h2

兲

⫺n2␲2HAk t h2

0 ⱕ t ⱕ t0

n␲ z , t0 ⱕ t h (A6)

Although not shown here, when eqn (A6) is evaluated, the strain field behaves in a time-dependent spatially varying manner consistent with the schematic predictions of Mow et al. (1980, 1990).

APPENDIX B There is a production type-setting error in eqn (2) in Berry et al. (2006). The correct equation for the gradient operator is given by:

0 ⱕ t ⱕ t0

n␲ z , t0 ⱕ t h (A5)

ⵜ ⫽ ជı

⭸ ⭸ ⭸ ⫹ ûជ ⫹ kជ ⭸x ⭸y ⭸z

(B1)

where ជı , ûជ and kជ represent the corresponding unit vectors along each of the Cartesian coordinate axes and should not be confused with the permeability, k, and hydraulic gradient, i, introduced previously (Berry et al. 2006).