Partial Volume Effects in Dynamic Contrast Magnetic Resonance Renal Studies

Partial volume effects in dynamic contrast magnetic resonance renal studies

D. Rodriguez Gutierreza, 1, ,
K. Wellsa, 2, ,
O. Diaz Montesdeocab, 3, ,
A. Moran Santanab, 4,
I.A. Mendichovszkyc, 5, ,
I. Gordon, ,

Abstract

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% inr2fitting 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.

Keywords

Magnetic resonance imaging;
Point spread function;
Partial volume effect;
GFR

Dynamic contrast enhanced magnetic resonance imaging (DCE-MRI) 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 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 using51Cr-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 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]and[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. 1provides 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.5T 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.5mm thickness (no gap), an in-plane voxel dimension of 1.56mm×1.56mm, TR=3.34ms, TE=1.67ms and an acquisition bandwidth of 590Hz/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.63ms, flip angle=17°, acquisition matrix=128×104, 400mm×325mm FOV and 1500Hz/pixel bandwidth. This dynamic image data were acquired every 2.5s over a period of several minutes, with each dynamic dataset consisting of 3D volumes with 18 slices of 7.5mm thickness (no gap) and an in-plane voxel dimension of 3.1mm×3.1mm. 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.

Fig. 1.

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

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 inFig. 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 inFig. 3)[13]. This 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 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.

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 inFig. 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.

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 thez-direction. The image on the right is an example of the gaussian–sinc mixture model (line) fitted to the data.

In order to locate the above mixing map at the correct position within the functional image data, the functional data were up-sampled 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 thezdirection), 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–180s of the data, as shown inFig. 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 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–90s 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:

…1

wherexandyrepresent 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 (seeFig. 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 thez-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ξpandξp+1that minimises the difference between them:

….2

wherexiis denotes theith location along profilep. Oncex′ was foundξp+1was 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 inFig. 3.

The sampling ofk-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 acquiredk-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.

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

…..3

wherebGS,hS,hG,μS,μG,σS, andσGare the fitting parameters