Partial volume effects in dynamic contrast magnetic resonance renal studies

European Journal of Radiology 75 (2010) 221–229

Contents lists available at ScienceDirect

European Journal of Radiology journal homepage: www.elsevier.com/locate/ejrad

Partial volume effects in dynamic contrast magnetic resonance renal studies D. Rodriguez Gutierrez a,1 , K. Wells a,2 , O. Diaz Montesdeoca b,3 , A. Moran Santana b,4 , I.A. Mendichovszky c,5 , I. Gordon c,∗ a

CVSSP, Faculty of Engineering & Physical Sciences, University of Surrey, UK EUITT, Universidad de Las Palmas de Gran Canaria, Spain c Radiology and Physics Unit, UCL Institute of Child Health, London WC1N 1EH, UK b

a r t i c l e

i n f o

Article history: Received 2 February 2009 Received in revised form 3 April 2009 Accepted 21 April 2009 Keywords: Magnetic resonance imaging Point spread function Partial volume effect GFR

a b s t r a c t This is the first study of partial volume effect in quantifying renal function on dynamic contrast enhanced magnetic resonance imaging. Dynamic image data were acquired for a cohort of 10 healthy volunteers. Following respiratory motion correction, each voxel location was assigned a mixing vector representing the ‘overspilling’ contributions of each tissue due to the convolution action of the imaging system’s point spread function. This was used to recover the true intensities associated with each constituent tissue. Thus, non-renal contributions from liver, spleen and other surrounding tissues could be eliminated from the observed time–intensity curves derived from a typical renal cortical region of interest. This analysis produced a change in the early slope of the renal curve, which subsequently resulted in an enhanced glomerular filtration rate estimate. This effect was consistently observed in a Rutland–Patlak analysis of the time–intensity data: the volunteer cohort produced a partial volume effect corrected mean enhancement of 36% in relative glomerular filtration rate with a mean improvement of 7% in r2 fitting of the Rutland–Patlak model compared to the same analysis undertaken without partial volume effect correction. This analysis strongly supports the notion that dynamic contrast enhanced magnetic resonance imaging of kidneys is substantially affected by the partial volume effect, and that this is a significant obfuscating factor in subsequent glomerular filtration rate estimation. © 2009 Elsevier Ireland Ltd. All rights reserved.

Dynamic contrast enhanced magnetic resonance imaging (DCEMRI) of the kidneys has been suggested as a possible alternative to radioisotope renography (RR) for estimating differential renal function (i.e. filtration). In DCE-MRI the injection of DTPA, a pure glomerular filtrate, tagged to gadolinium (Gd) allows MRI renography to be undertaken. Data is acquired as a function of time and the distribution of the Gd-DTPA in the kidney reflects different aspects of renal function depending on the time frame selected for analysis. Whilst DCE-MRI has the same basis as RR, advantages of DCE-MRI compared to RR include the lack of ionising radiation, increased

∗ Corresponding author. Tel.: +44 0 207 905 2202. E-mail addresses: [email protected] (D.R. Gutierrez), [email protected] (K. Wells), [email protected] (O. Diaz Montesdeoca), [email protected] (I.A. Mendichovszky), [email protected] (I. Gordon). 1 Biotronics3D Ltd. The Gatehouse, 64 Orchard Place, E14 0JW London, UK. Tel.: +44 0 207 093 0903. 2 Centre for Vision, Speech and Signal Processing, Faculty of Engineering & Physical Sciences, University of Surrey, Guildford GU2 7XH, Surrey, UK. Tel.: +44 0 1483 686036; fax: +44 0 1483 686031. 3 c./Gomera, 15, 35100 Maspalomas, Las Palmas, Spain. Tel.: +34 650160207. 4 C./Goya, 58, 35200 Telde (Playa del Hombre), Las Palmas, Spain. Tel.: +34 616531388. 5 Tel.: +44 0 207 905 2202. 0720-048X/$ – see front matter © 2009 Elsevier Ireland Ltd. All rights reserved. doi:10.1016/j.ejrad.2009.04.073

spatial resolution and availability of volumetric data that contains both tracer kinetics and anatomical information. To date, DCE-MRI and, even RR techniques have failed to robustly estimate absolute glomerular filtration rate (GFR) when compared to ‘gold standard’ plasma-sampling methods using 51 Cr-EDTA or inulin clearance [1]. Prior work has identified various factors that need to be considered for accurate quantification of glomerular filtration rate on DCE-MRI. These include movement artefacts [2], the selection of a suitable region of interest (ROI) for the analysis [3] non-linear relationship between signal intensity and gadolinium concentration [4] and suitability of the models used to estimate GFR [5]. Within this area, the effects of image distortions due to non-linear phase from magnetic susceptibility have thus far been ignored, and these may be significant at the boundary of the kidney and in particular, for small cortical ROIs. Movement correction is often considered as the first pre-processing step needed for subsequent analysis, as respiratory motion causes kidney displacement mainly in the cranial-caudal direction. The selection of a suitable ROI for analysis has been reported as crucial and this has been comprehensively reviewed in the literature [6]. As glomerular filtration is essentially a cortical function, the common approach is for cortical ROIs to be generated. Whilst all the aforementioned aspects are recognised as confounding factors on GFR quantification, there has, until now, been no consideration of the

222

D.R. Gutierrez et al. / European Journal of Radiology 75 (2010) 221–229

partial volume effect (PVE) in DCE-MRI in the published literature. The PVE occurs where the signals from two or more tissues combine to produce a single image intensity value within a particular voxel. This is a result of the finite bandwidth of the image acquisition system, occurring mostly near the boundaries between tissues, and is present to a greater or lesser degree in all imaging modalities. In particular, functional imaging techniques, such as DCE-MRI, offer relatively poor spatial resolution, often sacrificed for better temporal resolution. The resulting larger point spread function (PSF), and concomitantly larger voxels (compared to the size of the voxel in the anatomical areas of interest), are more likely to present significant PVEs. This potentially impacts on ROI selection, and its interpretation for subsequent GFR estimation. Whether using manual, whole-kidney [7] or semi-automated segmentation of cortical and medullary regions [8], the result of these techniques is the definition of a ROI where every voxel is labelled as either belonging solely to one particular tissue class (inside the ROI) or not (outside the ROI). In contrast, partial volume (PV) analysis assigns a set of mixing values to each voxel, representing the corresponding fractional signal component from various adjacent tissue structures that contribute to the observed voxel signal intensity. This robustly accounts for the convolution action of the system’s point spread function and the signal mixing this produces in adjacent structures. Whilst there is significant interest in PV for image quantification of brain data, PV quantification techniques have yet to be investigated for DCE-MRI renography. The brain analysis methods can be considered to fall into two basic approaches: those that utilise probabilistic approaches to infer the most likely tissue composition of a voxel [9,10], or those that utilise a mixing map estimated from anatomical data registered onto functional data [11]. The aim of the work presented here is to demonstrate the magnitude of the PV effect in DCE-MRI renography using a method based on the mixing map approach. We propose a method for PV correction in dynamic DCE-MRI data using anatomical MRI images and knowledge of the PSF of the particular MRI pulse sequence used to acquire the dynamic data. These are used to calculate mixing vectors for each voxel within a user-defined ROI and to de-compose the observed intensity for each voxel into its constituent parts corresponding to the contributions of each influencing tissue. We assessed the PV effect using the Rutland–Patlak approach to estimate GFR. 1. Materials and methods 1.1. Overview of the methodology This approach employs several pre-processing steps before a relative measure of GFR can be undertaken, whether using PVE analysis or not. Fig. 1 provides a schematic diagram of the flow of data processing used in the methodology, the most significant aspects of which are described in the ensuing sub-sections detailed below. The method was applied to DCE-MRI data obtained from scanning a cohort of 10 healthy volunteers on a 1.5 T Siemens Avanto scanner running a syngo MR 2004V 4VB11A platform. Each human volunteer was first scanned using a true-FISP pulse sequence to acquire breath-hold anatomical data. This 2D data consisted of 18 slices with 7.5 mm thickness (no gap), an in-plane voxel dimension of 1.56 mm × 1.56 mm, TR = 3.34 ms, TE = 1.67 ms and an acquisition bandwidth of 590 Hz/pixel. Oblique-coronal DCE-MRI data volumes were then acquired using a SPGR 3D-FLASH Volumetric Interpolated Breath-hold Examination (VIBE) with TE/TR = 0.53/1.63 ms, flip angle = 17◦ , acquisition matrix = 128 × 104, 400 mm × 325 mm FOV and 1500 Hz/pixel bandwidth. This dynamic image data were acquired every 2.5 s over a period of several minutes, with each

Fig. 1. Flow chart showing the main processing stages used in the PV correction method.

dynamic dataset consisting of 3D volumes with 18 slices of 7.5 mm thickness (no gap) and an in-plane voxel dimension of 3.1 mm × 3.1 mm. The in-plane voxel dimension was deliberately set to be twice that of the anatomical data for reasons given below. The dynamic data was then motion corrected [12]. The anatomical and dynamic scans were optimized (18 slices) to ensure full coverage of the kidneys in all volunteers, regardless of the kidneys’ dimensions, anatomical position or spatial orientation. The anatomical data (from the FISP sequence) were manually segmented, under clinical expert guidance, into a set of binary tissue templates including kidney, liver, spleen and background, as illustrated in Fig. 2. Our definition of ‘background’ represents any voxel signal emanating from a signal that is neither kidney nor spleen/liver. Each tissue template was then convolved with an experimentally determined 3D PSF representing the intrinsic impulse response of the dynamic SPGR 3D-FLASH (VIBE) pulse sequence used for dynamic acquisition (details of the experimental PSF phantom work appear below, illustrated in Fig. 3) [13]. This

D.R. Gutierrez et al. / European Journal of Radiology 75 (2010) 221–229

Fig. 2. (Top row) Tissue templates overlaid on anatomical MR data on left. The right image shows the large rectangular ROIs used in the latter part of the analysis, and a colour-coded representation of the number of tissues contributing to each voxel (with contributions >1%). Voxels containing signal contributions from one tissue are in red, two tissues in green and three tissues in blue. Note this only illustrates the number of tissues being PV mixed – the overall magnitude of these mixing components is shown in Fig. 4. Nonetheless, this demonstrates that the kidney/liver and kidney/spleen boundaries, that corresponds to the region most often used for ROI analysis of kidney function, contains mixtures from three different tissue components. (Bottom row) Functional image data used to manually define the cortical ROIs used in the analysis. Colour-coding corresponds to voxels obtained from the above mixing map. This demonstrates the heterogenic nature of PV mixing in the cortical region. Without PV correction, contamination of the assumed renal signal by the influence of the surrounding tissue signals would produce significant variation in estimated GFR, dependent on size and position of the ROI used for analysis.

produced a blurred representation of the spatial influence of each functional tissue class, with values in the range [0, 1]. Superimposing each convolved template onto the corresponding organ position in the original anatomical image produced a mixing map representing the action of the PSF on the functional data. Where the blurred templates overlap, the relative magnitude of each overlapping contribution in each voxel describes the relative signal contribution

223

from each adjacent tissue class to the observed intensity. Where there is no overlap in the templates, then the voxel may be considered to be a pure voxel, representing functional information from only a single tissue class. Thus mixing vectors were assigned to each voxel, with mixing vector components representing the fractional contribution of the adjacent tissues to the observed voxel intensity. Mixing vector components of less than 1% were truncated to zero. Note that in this work, in common with other PV correction methods and compartmental models of the kidney, we make an assumption of instantaneous mixing as well as a spatially invariant time–intensity distribution for each single tissue type. In order to locate the above mixing map at the correct position within the functional image data, the functional data were upsampled to the same resolution as the template/anatomical data. For registration, a method similar to that previously described [14], based on complex phase-image correlation, was extended to 3D to allow for through plane motion, and then employed here for each functional volume. Acquiring the anatomical data at 2× finer voxel dimensions and in the same slice orientation compared to the functional data allows a simple bilinear interpolation to be used in the up-sampling stage. This also provides greater control over locating the mixing map at the correct position on the functional data, giving the registration potentially sub-voxel functional precision in the coronal plane, and allows for slice motion to be accounted for. Final placement was checked by visual inspection. In about half the cases this required a further manual translation of one to two voxels to improve visually assessed alignment. Given that the registrations were all manually checked, it is unlikely that any of these were misplaced by more than one functional voxel volume, particularly as the data had been up-sampled to allow finer positioning in the final registration. However, to examine this potential source of error, deliberate misplacement of the cortical ROI by up to two voxels in each of the three orthogonal axial directions was applied. This produced changes in the mixing components of typically a few percent in individual voxels, except where the displacement produced a gross mismatch in overlapping regions (in the z direction), in which case the error was clearly visible. For subsequent analysis, a cortical ROI was then manually defined using a central kidney slice taken from the first 160–180 s of the data, as shown in Fig. 2. As the functional images had already been registered onto the aforementioned anatomic data, then this cortical ROI was also, implicitly, correctly registered. Using this ROI, observed time–intensity curves could then be extracted. The robustness of this approach to reproducibly recover PV-corrected time–intensity data has also been tested by using an alternative large rectangular ROI for PV analysis that encompasses the entire

Fig. 3. Phantom used for LSF estimation (left). Horizontal edges in the black squares labelled ‘a’ were used for estimating the LSF at 0◦ , and vertical edges labelled ‘b’ for 90◦ . The phantom was then manually rotated to capture the edge in the z-direction. The image on the right is an example of the gaussian–sinc mixture model (line) fitted to the data.

224

D.R. Gutierrez et al. / European Journal of Radiology 75 (2010) 221–229

kidney and contains significant non-renal contamination. This is described in further detail below. An arterial input function (AIF) for each volunteer was extracted from an aortic ROI in the functional data and used in conjunction with Rutland–Patlak analysis [15] to produce an estimate of relative GFR. The period corresponding to the first 30–90 s after baseline shift of the dynamic data was then used for Rutland–Patlak analysis, where contrast enhanced PV ‘overspill’ contamination from the medullary region was considered negligible. PVE-corrected and non-corrected Rutland–Patlak plots have then been compared for each volunteer. Details of the preliminary methodology used. 1.2. Estimation of the Point Spread Function A linear, shift-invariant imaging system can be completely characterised by its PSF, i.e. the system response to a unit impulse. Since any signal can be described as a series of shifted infinitesimal impulses, the output of a linear, shift-invariant system can be considered a sum of PSFs, each shifted to the location and scaled according to the height of the corresponding impulse. Although the linear shift-invariant imaging system is an idealised concept and is, in practice, never fully realised, many MR sequences can nonetheless be considered sufficiently shift-invariant to be commonly characterised by their PSF [16]. In practice, the PSF can be estimated from a set of Line Spread Functions (LSFs) at different orientations. The LSF is an integrated profile of the PSF along a given direction:



∞

LSF(x) =

PSF(x, y, z) dy dz

(1)

−∞

where x and y represent orthogonal directions. By measuring multiple LSFs at various angles, the PSF can be constructed. However, in the case where the PSF is circularly symmetric, a single LSF orientation may be used to characterise the PSF. To obtain LSFs, an edge phantom containing perspex blocks (see Fig. 3) was imaged using the same acquisition sequence as used for dynamic renal volunteer data acquisition as described above. LSFs were estimated by differentiating edge profiles at 0◦ and 90◦ orientation in-plane, and for the z-direction. To minimise the effects of noise and improve sampling, multiple edge profiles along the edge were co-registered before differentiation. This was achieved by finding the relative displacement, x , between two adjacent edge profiles  p and  p+1 that minimises the difference between them: argmin x

  p (xi ) − p+1 (xi )

To model these effects, a gaussian–sinc mixture, LSFGS (x), was used to model the experimental LSF data: LSFGS (x)=bGS +hS × sinc

x s



−S +hG ×exp



−(x − G )2 2G2



(3)

where bGS , hS , hG, S , G ,  S , and  G are the fitting parameters. This empirical model was fitted to each individual LSF for all three component directions (coronal, sagittal and trans-axial directions) using iterative least squares. These fittings were performed on different edge locations within a slice and on different slices of the phantom data and the average result for each orientation was used to construct a 3D model gaussian–sinc PSF with a resultant Full Width at Half Maximum (FWHM) of 6.8 mm in the x and y (coronal in-plane) directions and 17.2 mm in the z-direction. Fig. 3 illustrates the phantom used and the LSF model fit to the data [13]. 1.3. Human data acquisition Oblique-coronal anatomical and DCE-MRI data volumes from healthy volunteers with assumed normal renal function were acquired on a 1.5 T Siemens Avanto scanner (Siemens Medical Solutions, Erlangen, Germany) as described above. The cohort consisted of 10 volunteers in the age range 22–36 years (mean 27.8 years) with weight varying between 46 and 109 kg (mean 74.3 kg). The use of human volunteers for this study was approved by the local ethics committee. The oblique-coronal orientation along the long axis of the kidney was used to minimise movement artefacts due to respiration. The motion in this plane was considered to be dominated by 2D translation in the coronal plane. However, the motion correction step using a method previously described [12] considered both 3D translation and rotation, using an assumption of rigid body behaviour, consistent with other prior work in this area [2,3,12]. The anatomical scan was acquired as described above. It was followed by the acquisition of contrast enhanced dynamic data using a SPGR 3D-FLASH pulse sequence (VIBE) as described above. This acquisition protocol was deliberately designed to encompass both kidneys and aorta, so that arterial input function (AIF) data could be extracted from a single dynamic renal volume, and moreover, to allow the same acquisition protocol to be used on volunteers with significantly different internal organ locations and/or body mass indices. The Gd-DTPA dose (Magnevist, Schering, Germany) was 0.05 mmol/kg body weight, injected as a bolus at 2 ml/s using an automatic injector (Spectris). The contrast agent bolus was immediately followed by a 15 ml saline flush injected at the same speed.

(2)

i

where xi is denotes the ith location along profile p. Once x was found  p+1 was shifted accordingly and the same process was applied to the next profile. This resulted in a satisfactorily well-sampled edge profile that could be successfully differentiated to produce a reliable estimate of the corresponding LSF profile as shown in Fig. 3. The sampling of k-space during image acquisition results in a finite set of spatial frequencies being captured. Thus, it is common in the literature to consider the PSF as a sinc function; that is, the inverse Fourier Transform of an ideal low-pass filter. However, in addition to this, there are further filtering effects such as relaxation resulting in exponential signal decay during acquisition. Other factors that may also contribute to the observed PSF include not only signal decay during acquisition, but limited bandwidth of acquired k-space data, other low-pass filtering applied on acquisition, exponential signal decay and magnitude operator on complex Fourier-transformed data and mapping into an integer interval for image display.

1.4. Template and ROI definition The anatomic MR data demonstrating major organs (kidneys, liver and spleen) were manually segmented on all slices that contained kidney tissue, shown in Fig. 2. The higher resolution of the anatomic data compared to the functional data minimises any subsequent misrepresentation of the kidney boundary. This might occur, for example, as the kidney is often surrounded by a layer of fat, which may produce phase-cancellation effects resulting in loss of boundary signal. The resulting stack of tissue templates for each organ was then convolved with the aforementioned 3D PSF, to define the spatial extent of PSF-induced blurring, and hence PVE signal mixing, described in further detail in the next section. ROI were difficult to draw on the most ventral and dorsal slices of the kidneys in the oblique projection where there were relatively small amounts of kidney present. However, by selecting the central kidney slice for subsequent analysis, the effects of these most outlier slices on PVE-induced signal mixing was considered negligible.

D.R. Gutierrez et al. / European Journal of Radiology 75 (2010) 221–229

For subsequent PV analysis, two manually defined ROIs taken from a central slice for each kidney were used: a cortical ROI as defined by one of the authors (IG) and a large rectangular ROI. To define a cortical ROI (Fig. 2, bottom row), normally indistinguishable in an anatomic MR image, the first 50–70 s of the dynamic data following baseline shift were used, demonstrating enhancement of the cortex by the assumed dominant effect of glomerular filtration, compared to renal perfusion and tubular enhancement. A cortical ROI (as opposed to a parenchymal one) was chosen for two reasons: the cortex is the site where glomerular filtration takes place, and, we wanted to specifically select ‘cortical’ voxels to illustrate the PV effect on filtration estimation. The large rectangular ROIs considered here consisted of a rectangular abdominal region surrounding the entire kidney and deliberately encompassing adjacent anatomy (Fig. 2, top right). This was done to test the subsequent recovery of the true renal signal using the proposed PV correction method as we recognised that the large ROI contained high levels of signal corruption from PV mixing of adjacent tissue. No attempt was made to investigate the PV effects from the internal (non-cortical) kidney tissue. 1.5. Data analysis Assuming linear mixing of the contributing MR signals within a particular voxel, the problem of partial volumes might be thought of as finding for every voxel ω, a mixing vector, ␣ω : ˛ω = {˛ω1 , ˛ω2 , . . . , ˛ωN }

(4)

where N is the number of tissue signals contributing to the overall observed signal intensity in a particular voxel. Thus, the mixing vector describes the contribution of each tissue class j to the observed intensity of a particular voxel ω, so that:

 N

˛ωj = 1

(5)

j=1

To generate the mixing vectors for each voxel location, the structural tissue templates described above were registered with the movement-corrected dynamic data described above and convolved with the 3D PSF. The superposition of these convolved templates onto the corresponding anatomical locations produces overlap at their boundaries, and thus creates a component mixing map or ˛map of individual mixing components ˛ωj for each voxel location ω. The number of individual components assigned to each voxel (in this case, from one to a maximum of three) is illustrated by colour coding in Fig. 2. In this particular work, N = 3, and Eq. (5) is used to calculate the magnitude of the background mixing component, once kidney and liver/spleen components have been determined from the overlap of their convolved tissue templates. In this context, the background is defined as the complement of the kidney and liver/spleen templates, i.e. the group of voxels excluded from the templates. Having obtained the corresponding mixing maps, the observed intensity oω (t) for each voxel location, ω, at time t, within a ROI with ˝ voxels within or surrounding the kidney can be defined as o␻ (t) =

N 

˛ωj × ij (t)

where o(t) = [o1 (t), o2 (t), . . ., o˝ (t)]T , i(t) = [i1 (t), i2 (t), . . ., iN (t)]T , and A is an ˝ × N matrix with elements ˛ωj . Provided there are more observations per ROI than tissue types, that is ˝ ≥ N, then solving for i(t) is an over determined problem, where in solutions can be iteratively found, in a least-squares sense, which then minimises the effects of measurement errors. Thus solving Eq. (7), produces an estimate of the magnitude of the true signal tissue signal components, ij (t) at time t, that contribute to the observed signal intensity. Repeated over each scanning interval leads to estimates of the unmixed (i.e. PV-corrected) time–intensity curves for each component tissue. Note that within this framework and the aforementioned assumptions of instant mixing, etc., once templates and mixing maps have been defined, then ROI definition becomes a more arbitrary process: ROI definition simply selects those ˝ voxels (with ˝ ≥ N) that the user wishes to include in Eq. (7) to solve for the N (=3 in our case) tissue classes to be recovered. This allows us to produce PV-corrected estimates of the time–intensity behaviour of each of the N tissue components being considered. In order to quantify the effect of PV correction on the dynamic data, we have used the Rutland–Patlak [15] approach that will estimate filtration of a renal curve. This is a well-known model that has been successfully used in brain and kidney functional analysis. It is simple to apply, produces a closed-form solution to estimate (in this case) relative GFR without need for a multi-parameter estimation stage, apart from an analytical straight-line fit. The model only assumes unilateral flow from the arterial compartment to the renal cortex and that no filtered tracer leaves the ROI. Provided that analysis is undertaken within the appropriate timeframe, then the above assumption will remain approximately valid [17]. Thus, glomerular filtration is estimated from plotting the following straight-line equation: RROI (t) = k1 × AROI (t)

t 0

AROI () d AROI (t)

+ k2

(8)

where AROI (t) and RROI (t) represent the relative mean enhancement obtained from aortic arterial and renal ROIs, respectively, defined as {S(t) − S(0)}/S(0), following the method of de Priester [18]. Using Rutland–Patlak analysis between 30 and 90 s after the first rise of the signal from the baseline, S(0), then leakage of contrast agent from the cortex is considered insignificantly small, and thus k1 represents an estimate of GFR, and k2 an estimate of the blood pool contribution. We use a semi-quantitative analysis based on the renal mean enhancement for the ROI which is sufficient for the purpose of demonstrating PVEs in the analysis of DCE-MRI renal imaging. PVEs within the aorta have also been ignored in this study as these represent a simple scaling factor common to both raw and renal PV-corrected dynamic data. It is perhaps also worth explicitly stating that this work is confined to relative estimates of GFR quantification, as absolute quantification would otherwise require conversion from voxel intensity to absolute [Gd] and accurate knowledge of the total 3D volume of the cortex. These are non-trivial issues for absolute GFR estimation. The processing times for PV correction were approximately 20 s per kidney using Matlab (The MathWorks, Natick, MA) implementation on 32 bit dual core 1.86 GHz 3 Gb Windows Vista PC.

(6)

j=1

2. Results

where ij (t) represents the ‘true’ unmixed underlying intensity for tissue class j, at time t, and ˛ωj is the weighting or influence of a particular tissue class, j, at location ω. In matrix form, many individual voxel observations can be combined into a set of equations of the form: o(t) = A × i(t)

225

(7)

The colour coding in Fig. 2 illustrates the number of tissue components that influence the observed individual voxel intensities with the two ROIs used for analysis. This demonstrates there is significant presence of adjacent tissue signals in the boundary region between the liver/spleen. Fig. 4 quantifies this influence, which represents a typical distribution seen within the volunteer datasets,

226

D.R. Gutierrez et al. / European Journal of Radiology 75 (2010) 221–229

Fig. 4. Corresponding bar charts demonstrating tissue signal composition from the ROIs in Fig. 2 (i.e. the sum of all mixing values for each tissue inside an ROI) used in the analysis. Top row—Cortical ROI: demonstrates the contributing tissue class signals from the right and left kidney, liver, spleen and background. Bottom row—Large ROI (encompassing the entire kidney and surrounding anatomy). Note the significant non-renal components in each of these ROIs which need to be accounted for in any subsequent analysis.

showing a breakdown of the mixing components summed across each ROI. In both ROIs there is significant presence from non-renal tissue. Table 1 further illustrates the variation observed in tissue composition for the single-slice cortical ROIs across the cohort of volunteers for the single kidney slice. As can be seen, the perirenal components typically contribute approximately 20–30% to the overall tissue composition within the ROIs, demonstrating significant contamination from adjacent tissues in what was thought to be a ‘pure’ cortical ROI. Observed time–intensity plots were extracted from the cortical and rectangular ROIs for each volunteer dataset. Using Eq. (7), the observed time–intensity curves were then de-composed into their constituent parts, illustrated in Fig. 5 (top row), using the aforementioned technique. The middle row of Fig. 5 illustrates the true assumed (unmixed) temporal signals from spleen/liver ROIs. These were extracted from ROIs taken from the central regions of these large organ templates, where we assume the influence of PV mixing is negligible. These were compared to the extracted component liver/spleen signals (following PV analysis) from the left/right rectangular ROIs. These demonstrate excellent agreement between the

Table 1 Overall tissue contributions (as a percentage) to the cortical ROIs (shown in Fig. 2) from kidney, liver and background. Volunteer

Right kidney cortical ROI

Left kidney cortical ROI

Kidney

Liver

Background

Kidney

Spleen

Background

1 2 3 4 5 6 7 8 9 10

80.1 69.7 65.7 82.9 76.3 75.6 79.3 72.1 81.0 76.3

0.8 3.8 1.3 0.7 1.1 0.9 1.1 3.2 1.5 4.2

19.1 26.5 33.0 16.4 22.6 23.5 19.6 24.7 17.5 19.5

75.3 73.6 64.7 81.2 67.8 76.0 78.3 72.7 78.0 78.5

0.3 0.0 0.4 0.0 2.4 0.1 0.7 0.7 2.8 1.3

24.5 26.4 34.9 18.8 29.8 24.0 21.0 26.6 19.2 20.2

Mean Std. Dev.

75.9 5.4

1.9 1.3

22.2 5.0

74.6 5.1

0.9 1.0

24.5 5.1

Approximately a quarter of the contribution to the MR signal represented within those ROIs originates from non-renal tissues.

pure or unmixed signals and extracted PV-corrected estimates from the large rectangular ROI know to be heavily corrupted (see Fig. 4), supporting the assumptions used in the analysis, in particular the empirical PSF model and its assumed shift-invariance. Removing the effect of these components using the aforementioned PV correction causes a change in the shape of the recovered parenchymal renal curves with a steeper slope in the critical period of filtration as can be seen from the time–intensity plots in the bottom row of Fig. 5. This also illustrates the similarity between PV-corrected time–intensity plots for the cortical ROIs and the large rectangular ROIs, demonstrating the robustness of the approach to the subjective task of renal cortical ROI definition: this is despite the rectangular ROI suffering from >50% non-renal contamination. Rutland–Patlak analysis was applied to the PV-corrected and uncorrected data, with estimated k1 (GFR) produced. Table 2 illustrates the robustness of this approach to ROI selection, taking the data from Fig. 5 (bottom row), showing excellent consistency in estimated GFR, despite widely different initial (uncorrected) estimates of k1 . In both cases the estimated GFR has increased and the r2 fit has improved. Table 3 examines the effect of the PV correction on Rutland–Patlak analysis across the entire volunteer cohort. In all cases, PV correction has produced an enhanced estimate of k1 and a better conformation to the Rutland–Patlak model, as shown by an improved r2 (Fig. 6). This is perhaps surprising, given the relatively simple imaging model presented here. However, the excellent consistency of the results across different tissues, ROIs and different volunteers suggest that, at this region or tissue-based level of analysis, the PVE may be the dominant obfuscating factor in GFR estimation using DCE-MRI methods.

Table 2 Effect of ROI selection on k1 estimation, before and after PV correction. ROI

Left kidney

Cortical ROI Large ROI

Right kidney r2

k1

r2

k1

Before

After

Before

After

Before

After

Before

After

0.015 0.005

0.023 0.023

0.90 0.84

0.96 0.97

0.014 0.005

0.022 0.021

0.93 0.80

0.97 0.98

D.R. Gutierrez et al. / European Journal of Radiology 75 (2010) 221–229

227

Fig. 5. Top row: The observed time–intensity plot (labelled ROI, plotted as continuous line) extracted directly from the large rectangular ROIs in Fig. 2. Using the PV analysis described in the text, PV-corrected responses were extracted from this curve for left and right kidney (labelled LK, RK, respectively) liver (L), spleen (S), and background (BK). The left column corresponds to the right kidney and the right column to the left kidney throughout the figure. Note the dramatic change in shape of the PV-corrected LK/RK plots compared to the flatter attenuated raw ROI plots. Middle row: Using the PV correction method, the recovered (PV-corrected) liver and spleen components () re-plotted from the top row, compared to liver and spleen curves obtained from ROIs containing only ‘pure’ voxels (i.e. only one signal component taken from the corresponding mixing map) within the liver and spleen (continuous line). This illustrates successful recovery of the true tissue component curves by PV analysis. Bottom row: The observed (i.e. uncorrected) time–intensity plot from the cortical ROIs (lower plots), and the PV-corrected time–intensity plots from cortical (+) and rectangular ROIs (upper continuous line). The similarity between both PV-corrected time–intensity plots illustrates the robustness of the approach to the subjective task of renal cortical ROI definition during the period used for Patlak analysis.

Table 3 Summary table showing the effects on k1 (relative GFR) before and after PV correction of cortical ROIs with reported percentage change. The r2 fitting consistently demonstrates better model fitting after the PV correction. This suggests that after the removal of the corrupting enhancement from other tissues, the resulting data is more consistent with the Rutland–Patlak model. Volunteer No.

Left kidney

Right kidney

k1 (×10−3 /s)

k1 (×10−3 /s)

r2

r2

Before

After

% inc

Before

After

Before

After

% inc

Before

After

1 2 3 4 5 6 7 8 9 10

4.8 6.0 5.2 7.6 3.6 7.6 8.4 4.8 6.0 6.0

7.2 9.2 8.4 1.2 7.2 1.4 1.2 7.6 7.6 9.2

35 36 39 33 48 44 32 40 21 34

0.92 0.91 0.91 0.83 0.82 0.75 0.95 0.77 0.88 0.90

0.95 0.96 0.96 0.89 0.95 0.88 0.98 0.88 0.91 0.96

4.4 6.0 5.2 6.0 2.8 7.2 8.0 4.4 6.8 5.6

7.2 1.1 9.2 7.6 6.0 11.2 11.2 7.6 8.8 8.8

37 44 43 24 52 35 29 41 22 39

0.90 0.92 0.95 0.79 0.72 0.80 0.97 0.81 0.93 0.93

0.95 0.96 0.97 0.86 0.90 0.88 0.98 0.90 0.94 0.97

Mean Std. Dev.

6.00 1.50

6.02 3.36

36.2 7.36

0.86 0.07

0.93 0.04

5.64 1.53

36.6 9.37

0.87 0.08

0.93 0.04

7.87 2.90

228

D.R. Gutierrez et al. / European Journal of Radiology 75 (2010) 221–229

Fig. 6. Exemplar Rutland–Patlak analysis showing before (top) and after (bottom) PV correction of the intensity data for the right (left column) and left (right column) kidneys from one of the volunteers. Note the change in slope following PV correction. The Pearson coefficients top left R2 = 0.928 top right R2 = 0.904. Following PV correction bottom left R2 = 0.974 bottom right R2 = 0.955.

3. Discussion The quantification of the renal curve derived from dynamic Gd contrast enhanced MR renography is the subject of many publications. As the majority of the glomeruli are located in the renal cortex, most authors use a cortical ROI when absolute quantification of GFR has been reported. However, the above study has shown that the dynamic signal from the kidney contains significant contributions from non-renal tissue, which may account for some of the variation observed in estimated GFR seen in the prior literature, dependent on ROI location and size, when PVEs are ignored. Although there has been recognition of the importance of the PVE within the brain, there is no quantification of its effect in the kidney on DCE-MRI in the literature to date. The above decomposition of the time–intensity curves from a user-defined renal ROI into its constituent components reveals the importance of considering partial volume effects on a ROI placed on the cortex. This study shows that within a central kidney slice containing a region of interest naively thought to be ‘pure’ cortex of the left and right kidney, there was substantial contribution from the adjacent tissues (∼25%). This effect can be expected to increase dramatically for those volumetric slices at the extreme ends of the extent of the kidney. However, accurate definition of these regions is non-trivial, and is the subject of on-going developments. Our analysis also shows that removing the additional signal contribution due to partial volume effects from surrounding tissues (liver/spleen and peri-renal) has a considerable effect (20–50% change) on the slope of the recovered cortical Rutland–Patlak curve in this phase of the examination, the critical time for the quantification of the GFR. We have also shown that using PV correction produced a better fit to the data in Rutland–Patlak model space,

suggesting the successful elimination of non-renal time dependent signal components. Conventionally, the signal from a cortical ROI, when applied to a model, uses an assumption that all voxels comprise the same tissue type, and so have the same signal; or at least the average signal is somehow characteristic of the region. Yet the outer cortical glomeruli are much denser than those in the juxta-medullary region and thus a lower signal peak will be produced in this latter region compared to the outer cortex. This is a limitation of all current kidney models. With further development the PV model presented here one could address this issue, although in this initial demonstration, the kidney is simplified as a single tissue type, and we have only addressed extra-renal factors. Further analysis considering PVEs from renal medulla and pelvic regions might be expected to produce a better representation of the true, or unmixed, cortical curve. However, during the crucial time necessary for analysis of GFR from DCE-MRI it can be assumed no gadolinium has yet reached the collecting system. Thus overspill from these regions might be considered negligible in this analysis. If the method described here were to be extended to a voxel-byvoxel estimate of functional behaviour (i.e. to produce parametric images), then other well known DCE-MRI effects ignored here (e.g. magnetic field susceptibility variations and potential PSF shift variance due to Gd concentration change), may need to be re-examined. Our paper is the first publication discussing partial volume effects in MR urography and thus, it was intended as a pilot study. In this paper we studied mainly the PV effects arising from the liver and the spleen, as these organs are most likely to contribute the most to the partial volume in the kidneys. However, further studies will be needed to fully understand the PV contributions from the renal medulla, renal collecting system, muscle and bowel

D.R. Gutierrez et al. / European Journal of Radiology 75 (2010) 221–229

in all anatomical directions (antero-posterior, lateral and craniocaudal).

[6]

4. Conclusion In conclusion, this study has demonstrated the importance of the PVE on the observed time–intensity curve and subsequent Rutland–Patlak analysis for DCE-MRI. Its importance when quantifying GFR by means of DCE-MRI was evident and has been widely ignored in the renal MRI literature. Renal cortical ROIs have a contribution from non-renal tissues that may be as high as 20–30% due to the PVE for the typical scanning sequences suggested for DCEMRI. A method has been presented to successfully remove these unwanted components that is robust to the subjective task of ROI definition. Removal of such non-renal PVEs from renal ROIs can significantly change the slope of renal time–intensity curves and so affect the estimation of GFR by DCE-MRI.

[7]

[8]

[9]

[10]

[11]

[12]

Acknowledgements The authors would like to thank Baudouin Dennis de Senneville and Nicolas Grenier from the University of Bordeaux for access and assistance with the motion correction software used in this work. Kidney Research UK have funded the scanning of volunteers.

[13]

[14]

References [1] Chantler C, Barratt TM. Estimation of glomerular filtration rate from plasma clearance of 51-chromium edetic acid. Archives of Disease in Childhood 1972;47(254):613–7. [2] Sun Y, Moura JMF, Ho C. Subpixel registration in renal perfusion MR image sequence. In: Proceedings of the IEEE International Symposium on Biomedical Imaging: Macro to Nano, vol. 1. 2004. p. 700–3. [3] de Priester JA, Kessels AGH, Giele ELW, den Boer JA, Christiaans MHL, Hasman A, van Engelshoven JMA. MR renography by semiautomated image analysis: performance in renal transplant recipients. Journal of Magnetic Resonance Imaging 2001;14(2):134–40. [4] Rusinek H, Lee VS, Johnson G. Optimal dose of Gd-DTPA in dynamic MR studies. Magnetic Resonance in Medicine 2001;46(2):312–6. [5] Annet L, Hermoye L, Peeters F, Jamar F, Dehoux JP, van Beers BE. Glomerular filtration rate: assessment with dynamic contrast-enhanced MRI and a

[15]

[16] [17]

[18]

229

cortical-compartment model in the rabbit kidney. Journal of Magnetic Resonance Imaging 2004;20(5):843–9. Michoux N, Vallee J-P, Pechere-Bertschi A, Montet X, Buehler L, Beers BEV. Analysis of contrast-enhanced MR images to assess renal function. Magnetic Resonance Materials in Physics, Biology and Medicine 2006;19(4): 167–79. Pedersen M, Shi Y, Anderson P, Stodkilde-Jorgensen H, Djurhuus JC, Gordon I, Frokiaer J. Quantitation of differential renal blood flow and renal function using dynamic contrast-enhanced MRI in rats. Magnetic Resonance Medicine 2004;51(3):5F–17. Lee VS, Rusinek H, Noz ME, Lee P, Raghavan M, Kramer EL. Dynamic threedimensional MR renography for the measurement of single kidney function: initial experience. Radiology 2003;227(1):289–94. Kwan RK-S, Evans AC, Pike GB. MRI simulation-based evaluation of imageprocessing and classification methods. IEEE Transactions on Medical Imaging 1999;18(11):1085–97. Van Leemput K, Maes F, Vandermeulen D, Suetens P. A unifying framework for partial volume segmentation of brain MR images. IEEE Transactions on Medical Imaging 2003;22(1):105–19. Quarantelli M, Berkouk K, Prinster A, Landeau B, Svarer C, Balkay L. Integrated software for the analysis of brain PET/SPECT studies with partial-volume-effect correction. Journal of Nuclear Medicine 2004;45(2):192–201. de Senneville BD, Mendichovszky IA, Roujol S, Gordon I, Moonen C, Grenier N. Improvement of MRI-functional measurement with automatic movement correction in native and transplanted kidneys. Journal of Magnetic Resonance Imaging 2008;28(October (4)):970–8. Rodriguez Gutierrez D, Diaz Montesdeoca O, Moran Santana A, Wells K, Mendichovszky I, Gordon I. MR-based renography as a replacement for radionuclide diagnostic renography studies. In: Proceedings of the IEEE Nuclear Science Symposium & Medical Imaging Conference Record, IEEE Volume 6, 3. 2007. p. 4556–63. Giele EL, de Priester JA, Blom JA, den Boer JA, van Engelshoven JM, Hasman A. Evaluation of two cortical fraction estimation algorithms for the calculation of dynamic magnetic resonance renograms. Computer Methods and Programs in Biomedicine 2002;67(3):169–76. Patlak CS, Blasberg RG, Fenstermacher JD. Graphical evaluation of blood-tobrain transfer constants from multiple-time uptake data. Journal of Cerebral Blood Flow Metabolism 1983;3(1):1–7. Steckner MC, Drost DJ, Prato FS. Computing the modulation transfer function of a magnetic resonance imager. Medical Physics 1994;21(3):483–9. Hackstein N, Kooijman H, Tomaselli S, Rau WS. Glomerular filtration rate measured using the Patlak plot technique and contrast-enhanced dynamic MRI with different amounts of gadolinium-DTPA. Journal of Magnetic Resonance Imaging 2005;22(3):406–14. de Priester JA, den Boer JA, Christiaans MH, Kessels AG, Giele EL, Hasman A, van Hooff HP, van Engelshoven JM. Automated quantitative evaluation of diseased and nondiseased renal transplants with MR renography. Journal of Magnetic Resonance Imaging 2003;17(1):95–103.