Nicolas Grenier

Iosif Mendichovszky

Baudouin Denis De Senneville

Sébastien Roujol

Pascal Desbarat

Michael Pedersen

Kevin Wells

Jorgen Frokiaer

Isky Gordon

Measurement of Glomerular Filtration Rate With Magnetic Resonance Imaging: Principles, Limitations, and Expectations

come   L'archive ouverte pluridisciplinaire

U N C O R R E C T E D P R O O F

reproducible information from the entire kidneys because of respiratory movements and to magnetic susceptibility artifacts arising from the surrounding bowel; second, the complex relationship between the observed MR signal changes and the concentration of the contrast agent makes quantitative analysis difficult; third, the different kinetic models that are proposed to calculate GFR from MRI data have to be validated; and fourth, the reported accuracy of these techniques is still unsatisfactory. Therefore, nuclear medicinebased techniques remain at present the reference method for quantification of GFR but MRI should be able to compete in the near future.

Several methods are available for measurement of GFR using MRI and freely filtered Gd-chelates: (1) measurement of the clearance of the agent using blood samplings; (2) measurement of the clearance of the agent using signal intensity (SI) changes within abdominal organs; (3) measurement of the extraction fraction of the agent; and (4) monitoring of tracer intrarenal kinetics.

Determination of Gd-Chelate Concentration

Quantification of physiological parameters such as GFR using contrast-enhanced MRI requires the conversion of SI changes into concentrations changes. The relationship between SI measured on MR images and Gd-chelate concentration is not linear due to T2 and T2* effects occurring at high concentrations. However, there is a linear relationship between and the T1 relaxation rate (R1 ϭ 1/T1) and Gd-chelate concentration. From R1 changes, quantification of the concentration of a contrast agent is possible according to the following linear relationship:

C ϭ ͑ R1 Ϫ R1 0 ͒ ⁄ r
where R1 0 is the bulk R1 relaxivity in the medium or the tissue without contrast agent and r is the specific relaxivity of the agent.

Dynamic R1 change values can be obtained either indirectly using conversion methods from SI, or directly by measuring R1 instead of SI with specific sequences. The indirect methods include the following:

• calibration using a phantom of several tubes, filled with Gd-solutions at various concentrations, and imaged with the same pulse-sequence as the one used for the in vivo study. The acquired SI values are plotted against measured R1 values, and a polynomial fit is made to obtain a calibration curve. A recent study showed that this approach provides results that correlate well with those obtained by a reference technique. [START_REF] Bokacheva | Quantitative determination of Gd-DTPA concentration in T1-weighted MR renography studies[END_REF] However, this approach assumes that the relaxivity of the contrast agent is equivalent in solution and in tissues and that the difference between room and body temperature will have little effect. • an analytical calibration based on the equation driven by the sequence used, [START_REF] Pedersen | Quantitation of differential renal blood flow and renal function using dynamic contrast-enhanced MRI in rats[END_REF] which unfortunately is not straightforward. In addition, some sequences may not even have a an analytical formula that can be used for this method. [START_REF] Rusinek | Optimal dose of Gd-DTPA in dynamic MR studies[END_REF] • R1 of a tissue also can be measured directly using an appropriate pulse-sequence, such as a gradient-echo one with several flip angle values, [START_REF] Buckley | Measurement of single kidney function using dynamic contrast-enhanced MRI: comparison of two models in human subjects[END_REF] or using a very fast inversion recovery pulse-sequence with several inversion times based on the Look-Locker method. 6

Measurement of Gd-Pasma Clearance

Because there is a linear relationship between Gd-chelate concentration and the T1 relaxation rate (R1 ϭ 1/T1), changes in R1 can be used for absolute quantification of the concentration of a contrast agent. Choyke and coworkers [START_REF] Choyke | Hydrated clearance of gadolinium-DTPA as a measurement of glomerular filtration rate[END_REF] showed that measurement of R1 within serial blood samples after intravenous injection of Gd-DTPA was feasible, safe for the kidney, and reliable when compared with 99 m Tc-DTPA clearance (y ϭ 1.04 ϫ 2.2, r ϭ 0.94; Fig. 1). This method presents little advantage compared with other methods because it is time-consuming and requires several blood samples.

Measurement of Gd-Renal Clearance

The renal clearance of a freely filtered Gd-chelate also can be calculated by measuring the slope clearance of the tracer from the extracellular fluid volume (ECFV) after an intravenous injection of the tracer. Using this principle, Boss and coworkers [START_REF] Boss | Quantitative assessment of glomerular filtration rate with MR gadolinium slope clearance measurements: a phase I trial[END_REF] measured, in 16 healthy volunteers, the SI changes, induced by the contrast agent (0.05 mmol/kg of gadobutrol),

U N C O R R E C T E D P R O O F
within several abdominal organs (liver, spleen, and renal cortex) during the course of 80 minutes. The results were compared with simultaneous iopromide clearance measurement as reference standard. The sequence used was a navigatorgated FLASH sequence (TR ϭ 498, TE ϭ 1.25, FA ϭ 8°) with a nonselective magnetization preparation and a spoiler for heavy T1-weighting. The GFR was calculated as the product of the ECFV (ECFV ϭ 0.02154 • weight 0.6469 • height 0.3964 ) and the time constant of the second exponential phase (␣ 2 ). The best concordance between gadobutrol clearance and io-promide clearance was observed within the liver, the exponential fit being performed between 40 and 55 minutes after injection (with a mean paired difference of Ϫ5.9 mL/min/ 1.73 m 2 Ϯ 14.6; Fig. 2). All measurement points were within Ϯ 2 standard deviation values but the maximum deviation form the reference GFR was 29%. With such a wide deviation, this technique cannot be applied to the individual patient. The main drawbacks of this method are the length of MR acquisition (approximately 70 minutes), the sensitivity to body movements, and to the selected time intervals when the analysis is performed. More experience is required with this technique.

Measurement of Extraction Fraction

This method allows a quantification of a single kidney (SK) extraction fraction (EF), based on the measurement of T1 within flowing arterial and venous blood (Look and Locker method) during a continuous Gd infusion [START_REF] Dumoulin | Noninvasive measurement of renal hemodynamic functions using gadolinium enhanced magnetic resonance imaging[END_REF][START_REF] Niendorf | Measurement of Gd-DTPA dialysis clearance rates by using a look-locker imaging technique[END_REF][START_REF] Niendorf | Rapid measurement of Gd-DTPA extraction fraction in a dialysis system using echo-planar imaging[END_REF][START_REF] Niendorf | Rapid in vivo measurement of single-kidney extraction fraction and glomerular filtration rate with MR imaging[END_REF] :

EF ϭ Ca Ϫ Cv/Ca
where Ca and Cv are the arterial and venous Gd concentration respectively. Because of the linear relationship between relaxation rates and concentration of Gd, EF can be expressed alternatively as:

EF ϭ ͑͑ 1 ⁄ T1a ͒ Ϫ ͑ 1 ⁄ T1v ͒͒ ⁄ ͑͑ 1 ⁄ T1a ͒ Ϫ ͑ 1 ⁄ T1 0 ͒͒
where T1a is T1 in the renal artery, T1v in the vein and T1 0 in the blood without Gd. Once EF is calculated for each kidney, values of SKGFR can be calculated if RBF is known:

SKGFR ϭ EF ϫ RBF ϫ ͑ 1 Ϫ Hct ͒
where RBF is the renal blood flow, usually measured with the cine-phase contrast method on the renal artery, and Hct is the level of hematocrit in the blood.

In practice, the T1 of the blood is measured in the infrarenal inferior vena cava for the arterial input (because the renal Figure 2 Measurement of renal clearance. (A) Exponential decay on the logarithmically scaled time-signal intensity curve (displayed for liver parenchyma) indicates linear behavior. Two exponential phases can be distinguished owing to different slopes: the fast mixing phase immediately after contrast medium injection and the more prolonged phase, in which the signal intensity decrease is caused solely by glomerular filtration. The negative slope of this phase is equal to the GFR divided by the ECFV. (B and C) Normalized GFR measured by using gadobutrol clearance for MR evaluation of an ROI in the liver at 40-65 minutes after contrast medium injection is compared with iopromide clearance. (B) Graph shows a plot of the GFR values obtained by using both methods. (C) Bland-Altman plot relates the difference between each measurement pair to the mean value. A mean difference in measurement pairs of Ϫ5.9 mL/min Ϯ 14.6 was calculated. All measurement points are within Ϯ 2 standard deviations from the mean value. (Reprinted with permission from Boss et al. [START_REF] Boss | Quantitative assessment of glomerular filtration rate with MR gadolinium slope clearance measurements: a phase I trial[END_REF] )

Measurement
U N C O R R E C T E D P R O O F
artery is too small) and on the renal vein for the venous output (Fig. 3A). Preliminary animal studies have shown concordant results compared with those calculated from inulin clearance, but the correlation was 0.77, which is inadequate for clinical use at this stage (Fig. 3B). [START_REF] Coulam | Noninvasive measurement of extraction fraction and single-kidney glomerular filtration rate with MR imaging in swine with surgically created renal arterial stenoses[END_REF] No clinical application of this method has been published until now.

Tracer Intrarenal Kinetics With Dynamic MRI

This method is based on modeling of an MR renogram into 2 or more compartments obtained after intravenous injection of a bolus of Gd-chelate.

Dynamic MR Acquisition

A great heterogeneity for parameters of pulse-sequences, doses of injected Gd contrast, methods for conversion of signal intensity into concentration, postprocessing methods, and compartment models is still noted in the literature between the different groups and no consensus exists.

Sequences

All the pulse-sequences used must have a heavy T1-weighting and be fast enough to characterize the vascular phase of the tracer kinetic, which is necessary for assessment of the arterial input function. A reasonable T1 weighting can be accomplished using spoiled gradient-echo (type FLASH, i.e., fast low-angle shot) and RF-flip angles between 45°and 90°, which provide a relative linear relationship between SI and Gd concentration. [START_REF] Rohrschneider | Functional and morphologic evaluation of congenital urinary tract dilatation by using combined static-dynamic MR urography: Findings in kidneys with a single collecting system[END_REF] However, to obtain very short acquisition times and a heavy T1-weighting, a nonselective magnetization preparation has to be used, combined with very short TR and TE values and low flip angles. The preparation pulse is either realized as a saturation pulse (90°) or as an inversion pulse (180°), or alternatively, a more complex pulse scheme (90°-180°) to get rid of the effects of arrhythmia on flow artifacts. [START_REF] Montet | Noninvasive measurement of absolute renal perfusion by contrast medium-enhanced magnetic resonance imaging[END_REF] Therefore, according to selected geometric parameters and the performance of the MRI system used, a compromise may be necessary.

Dose of Contrast Agent

The concentration of Gd within the kidney can be very high because ofwater reasbsorption in the proximal convoluted tubule and within the medulla. High concentrations also can be reached in the abdominal aorta during the bolus. Therefore, to avoid the T2* contribution of the agent occurring at high concentration, the injected dose must be lower than the regular clinical dose of 0.1 mmol/kg, and the patient must be well hydrated. In the literature, the dose of gadolinium used varied from 2 to 16 mL (for adults), with the Gd being diluted to different volumes (in 2 studies, the Gd was diluted to a final injection volume of 60 mL). The rate of injection in humans varied, being either automatic (1 or 2 mL/s) or manual. In the animal studies the dose was between 0.05 and 0.075 mmol/kg body weight with a manual injection technique. There is no consensus about the optimal dose to be used for that purpose. The choice between 0.025 mmol/kg and 0.05 mmol/kg depends on the level of signal-to-noise ratio obtained with the sequence and the system used.

Acquisition Volume

Some groups acquired only one median slice through both kidneys, and, calculating a mean GFR value, extrapolated this value to the entire renal volume to get a total GFR value. This approach is simpler, allowing a high temporal resolution, but neglects the probable heterogeneity of filtration

U N C O R R E C T E D P R O O F
function within the cortex. Also when dealing with diseased kidneys the assumption that the single slice accurately represents the entire kidney may be untrue. Therefore, others prefer to cover the entire parenchyma of each kidney, using multislice or 3D acquisitions, to be able to sum the GFR values of each slice, or even of each voxel in each slice. Comparison of these 2 approaches has never been performed.

Acquisition Plane

In the literature, the dynamic renal images are acquired either in an oblique-coronal plane passing through the long axis of the kidneys or in the transverse plane. Drawbacks of the later are the higher number of slices required for covering the entire kidneys, with a negative impact on the time resolution, the impossibility for movement correction, and the risk of severe inflow effects within the aorta altering the accuracy of the arterial input function.

Arterial Input Function (AIF)

Measurement of renal filtration requires an accurate sampling of the vascular phase of the enhancement to measure the AIF. The AIF is characterized by the SI changes within the suprarenal abdominal aorta. It is used for the different kinetic models to compensate for the noninstantaneous bolus injected into the blood. Not taking this into account will produce an overestimation of GFR because of recirculation of the agent within the vascular space. Sampling the AIF requires a high temporal resolution. Recent advancements in MR instrumentation and hardware facilitate coverage of the entire kidneys with several slices using 3D pulse-sequences while maintaining a temporal resolution around 1.5 seconds. [START_REF] Pedersen | Quantitative measurement of renal function using contrast enhanced MRI: an initial experience[END_REF] Coronal acquisition, coupled with a nonselective saturation or inversion pulse, has the advantage to get rid of inflow effects if sampling of the aorta is far enough from the top of the imaging volume (Fig. 4). These effects are more difficult to control in axial acquisitions.

Postprocessing of Data

If the time resolution of the pulse-sequence is sufficiently high, normal SI-time curve after Gd injection can be characterized by 3 phases (Fig. 4): a first abrupt ascending segment followed by a first peak, corresponding to the "vascular-toglomerular first-pass" or cortical vascular phase; a second slowly ascending segment, ended by a second peak, corresponding to the glomerulo-tubular phase; and a slowly descending segment, corresponding to the predominant excretory function and so-called "excretory phase." If the region of interest (ROI) is restricted to the cortex, the slope of the second ascending segment is reduced. If the ROI is restricted to the medulla, the curve is characterized by a biphasic ascending slope, corresponding to the vascular and the uptake phases, the later segment reflecting the tubular transit and progressive concentration, and is followed by the excretory phase.

These renal curves constitute the basis for functional measurements. Therefore, they have to be optimized to remove The cortical SI time curve can be characterized by 3 phases: a first abrupt ascending segment followed by a first peak, corresponding to the "vascular-to-glomerular first-pass" or cortical vascular phase; a second slowly ascending segment, ended by a second peak, corresponding to the glomerulo-tubular phase; and a slowly descending segment, corresponding to the predominant excretory function and so-called "excretory phase."

Measurement of GFR with

U N C O R R E C T E D P R O O F
the artifacts introduced by respiratory movements during data acquisition, partial volume effects, separation of segments with different functional status and application of appropriate pharmacokinetic models. These tasks will require the development of specific postprocessing techniques.

Movement Correction

Because dynamic MR imaging of the kidneys is performed during free breathing, the main problem with postprocessing these images is correction of respiratory movements. Because gated sequences [START_REF] Moricawa | Feasibility of simple respiratory triggering in MR-guided interventional procedures for liver tumors under general anesthesia[END_REF] would lead to severe penalties in temporal resolution, most of the groups studying renal function either repositioned images manually, on a visual basis before sampling, or applied no correction at all. Absence of correction produces artifacts in time-intensity evolution which can lead to incorrect quantification.

Although problems related to movement can in principle be addressed with image registration algorithms [START_REF] Maintz | A survey of medical image registration[END_REF] applied in a postprocessing step, additional challenges arise from the fact that image amplitude and contrast change with time due to the transit of the contrast agent during bolus passage. A movement correction method developed in our group [START_REF] De Senneville | Improvement of MRI-functional measurement with semi-automatic movement correction in native and transplanted kidneys[END_REF] is derived from a method proposed by Sun and coworkers, [START_REF] Sun | Integrated registration of dynamic renal perfusion MR images[END_REF] allowing the estimation of a rigid transformation and correction for movement in renal MR data with subpixel accuracy. The principle is the following: although the relative intensities between tissues vary with time, edge orientation along tissue boundary is always parallel across the image sequence. A ROI restricted to the renal cortex is manually drawn by the radiologist on an early enhanced image (during the vascular phase of renal enhancement which enhances strongly and specifically the cortex). A binary mask (noted m) was constructed: pixels of the mask inside the cortex had a value of one, and outside a value of zero. The reference image set for motion estimation of all images was the image on which the ROI had been manually set. The displacement of the kidney T was estimated by a 3 parameters rigid body model. Let ref (x, y) and M ref (x, y) be respectively the edge orientation and the edge magnitude for the pixel (x,y) obtained using a Sobel edge detector on the reference image. Identically, let cur (x, y) and M cur (x, y) be respectively the edge orientation and the edge magnitude for the pixel (x,y) obtained on the image to register. An exhaustive search is performed to determine the parameters of the spatial transformation T that maximize the following edge-based consistency metric:

␥ ͑ I ref ͑ x, y ͒ , I cur ͑ T ͑ x, y ͒͒͒ ϭ ͚ x, yʦm M ref ͑ x, y ͒ . M cur ͑ T ͑ x, y ͒͒ . cos͑2͑ ref ͑ x, y ͒ Ϫ cur ͑ T ͑ x, y ͒͒ ͒͒ ͚ x, yʦm M ref ͑ x, y ͒ . M cur ͑ T ͑ x, y ͒͒ (1)
A multistep scheme is used to perform a subpixel registration while reducing enumeration of possible solutions. With help of this accurate estimate, the ROI could be propagated to all images of the time series.

As kidney displacement is assumed to be rigid, the estimated 2D motion is used to register automatically all other slices. If the motion is correctly determined, the kidney border would coincide with the contour of the original region of interest manually drawn by the radiologist around the kidney. An experienced operator could also review the re-aligned data set and manually correct these if necessary.

This algorithm could estimate successfully small motion and more than 20% of reduction on GFR uncertainty on transplanted kidneys

U N C O R R E C T E D P R O O F
In the future,

Renal and Cortex Segmentation

In the clinical studies published in the literature, segmentation and ROI generation were generally done manually and the difficulties in this area were poorly addressed. Applying an adequate model to the renographic data requires an effective segmentation of renal compartments, cortex and medulla, to separate their respective function. Such a manual segmentation is time-consuming and therefore requires automatic methods. Furthermore, to avoid interoperator variability, several segmentations by hand must be done for each renographic data set. A recent review has developed and discussed extensively the principles and the difficulties on segmentation and on renal ROI generation. [START_REF] Michoux | Analysis of contrastenhanced MR images to assess renal function[END_REF] Several groups of methods are offered: simple thresholding, boundary-based techniques, region-based techniques, and other types of semiautomated methods. These algorithms are only geometry-based and, therefore, their efficiency is correlated with image contrast and SNR. A first approach to develop a high-level MR renography segmentation algorithm using graphs as been introduced. [START_REF] Rusinek | Performance of an automated segmentation algorithm for 3D MR renography[END_REF] If the result of the segmentation of the cortex and the medulla is encouraging, it could be further improved by taking into account the topology of the kidney. To do this, the graph representation of the kidney has to be oriented, resulting in a combinatorial map. [START_REF] Braquelaire | Image segmentation with topological maps and interpixel representation[END_REF] This approach should allow an accurate segmentation by considering the cortex as a connected set of voxels and medulla as "holes" (lack of signal). In the future, a real 3D segmentation has to be performed. Computation of a complete 3D renal reconstruction will be necessary to follow the cortex boundary in 3D and to avoid partial volume artifacts.

Partial Volume Effects

The partial volume effect occurs where signals from 2 or more tissues combine to produce a single image intensity value within a particular voxel. This effect becomes evident mostly at the boundaries between tissues, and is present to a greater or lesser degree in all imaging modalities. In functional MRI acquisitions, where spatial resolution is killed for better temporal resolution (such as magnetic resonance renography) the partial volume effect is likely to be significant. This potentially impacts on ROI selection, and its interpretation for subsequent GFR estimation. Work recently presented suggested that voxels of the right kidney cortex adjacent to the liver may have as much as 20% contribution from the liver. (Rodriguez Gutierrez D, Diaz Montesdeoca O, Wells K, et al: Partial volume correction for magnetic resonance renography studies. British Chapter of the ISMRM, 12th Annual Meeting, Guilford, August 23-25, 2006).

Compartment Models

The ideal model for the analysis of the renal dynamic contrast-enhanced MRI data in the kidney remains to be worked out. Most of the reviewed MRI methods used a 2-compartment model, although the more recent publication by Lee and coworkers [START_REF] Lee | Renal function measurements from MR renography and a simplified multicompartmental model[END_REF] proposed a more complex compartment model.

The Rutland-Patlak plot technique is a graphical method based on a 2-compartmental model with the assumptions that the rate of change of concentration in the kidney during the clearance phase is constant if the amount of contrast agent is taken into account during this period. [START_REF] Patlak | Graphical evaluation of blood-to-brain transfer constants from multiple-time uptake data[END_REF] This method has been applied to MRI by several groups. [START_REF] Annet | Glomerular filtration rate: assessment with dynamic contrast-enhanced MRI and a cortical-compartment model in the rabbit kidney[END_REF][START_REF] Hermoye | Calculation of the renal perfusion and glomerular filtration rate from the renal impulse response obtained with MRI[END_REF][START_REF] Pedersen | MR renography based on contrast-enhanced T1-mapping[END_REF][START_REF] Hackstein | Contrast media clearance in a single kidney measured on multiphasic helical CT: Results in 50 patients without acute renal disorder[END_REF][START_REF] Hackstein | Measurement of single-kidney glomerular filtration rate using a contrast-enhanced dynamic gradientecho sequence and the Rutland-Patlak plot technique[END_REF] Generally, the interstitial space is neglected, and Gd is considered to mix freely between the capillary and the glomerular compartments. This assumption may, however, fail in situations where the interstitial space is expanded (pyelonephritis and some examples of acute renal failure) which may lead to overestimation of SKGFR. The assumption that no contrast agent leaves the ROI during the sampling period may justify the use of ROIs encompassing both cortex and medulla. The model is realized as an x-y plot (the Rutland-Patlak plot) using the ratio of the renal concentration/aortic concentration plotted against the ratio of the integral of aortic concentration/aortic concentration. When applied to the second phase of the signal intensity-time curve (glomerulo-tubular uptake), this plot leads to a straight line, with a slope proportional to the renal clearance, and an intercept with the y-axis proportional to the cortical blood volume.

The cortical compartment model [START_REF] Annet | Glomerular filtration rate: assessment with dynamic contrast-enhanced MRI and a cortical-compartment model in the rabbit kidney[END_REF][START_REF] Hermoye | Calculation of the renal perfusion and glomerular filtration rate from the renal impulse response obtained with MRI[END_REF] is another 2-compartment model confined to the cortex. It derives from the Rutland-Patlak model but differs by the constraints that the outflow from the tubules, during the sampling period, is taken into account. This makes it possible to draw ROIs strictly limited to the cortex. The aortic curve is shifted and dispersed taking into account the transit time of the contrast agent between the aorta and the kidney and the tracer dispersion in the glomeruli. [START_REF] Annet | Glomerular filtration rate: assessment with dynamic contrast-enhanced MRI and a cortical-compartment model in the rabbit kidney[END_REF] Calculation of the residue function, (ie, the deconvolution of C(t) with AIF(t)) then exhibits 3 sequential peaks (successively glomerular, proximal tubule, distal tubule) and 2 rate constants, kin and kout that describe the flow into and out of the proximal tubule, meaning that kin represents SKGFR. SKGFR can be calculated according to the following equation: SKGFR ϭ maximal slope of proximal tubule peak/ C ͑ vasc ͒ max To calculate C(vasc)max, the maximum of the vascular peak was divided by the RBV (see chapter on perfusion) to take into account that the contrast agent remains in the extravascular space. This method has provided more accurate results than Rutland-Patlak model in an experimental study in rabbits, using 51Cr-EDTA as a gold standard. [START_REF] Annet | Glomerular filtration rate: assessment with dynamic contrast-enhanced MRI and a cortical-compartment model in the rabbit kidney[END_REF] The multicompartmental model 24 includes 3 cortical compartments (glomerular, capillary, and proximal convoluted tubules, as in previous model), 3 medullary compartments (loops of Henle, distal convoluted tubules, collecting ducts), and the collecting system. Despite being complex, this model has the advantage of assessing some important tubular physiological parameters. Applied to the passage of Gd into the first 2 compartments, the model allows calculation of GFR.

Calculation of Single Kidney GFR

The total SKGFR (in mL/min) can be obtained either by multiplying the mean GFR value, measured on a cortical region of

U N C O R R E C T E D P R O O F
interest, by the functional renal volume being either the renal cortex [START_REF] Hermoye | Calculation of the renal perfusion and glomerular filtration rate from the renal impulse response obtained with MRI[END_REF] or the total renal volume. [START_REF] Buckley | Measurement of single kidney function using dynamic contrast-enhanced MRI: comparison of two models in human subjects[END_REF] In the normal kidney, the cortical volume can also be considered as 70% of the total kidney volume. [START_REF] Hegedus | Cortical volume of the normal human kidney. Correlated angiographic and morphologic investigations[END_REF] Another method for single kidney GFR calculation consists in modeling, on voxel-by-voxel basis, each MR-slice of the kidney (providing renal functional maps) and to add the values from each voxel of each slice (Fig. 6). This method has the advantage it takes into account the heterogeneity of filtration within the organ. To our knowledge, these 2 approaches have still to be compared.

Accuracy and Reproducibility

After multiple blood sampling (2-4 hours) after a single intravenous bolus injection of 51 Cr-EDTA and a single compartment analysis for the measurements of GFR, Chantler and coworkers [START_REF] Chantler | Estimation of glomerular filtration rate from plasma clearance of 51-chromium edetic acid[END_REF] showed a coefficient of variation (CV) of 3.9%. Brøchner-Mortensen and Rødbro 33 showed a precision of 5.5% and a reproducibility (CV) of 4.1% in patients with a plasma clearance Ͼ30 mL/min, whereas with a plasma clearance Ͻ30 mL/min, a precision of 11.6%, and a reproducibility of 11.5% were obtained. Blake and coworkers [START_REF] Blake | Long-term precision of glomerular filtration rate measurements using 51Cr-EDTA plasma clearance[END_REF] showed that the long-term precision of GFR based on the 51 Cr-EDTA clearance had a CV of 9.8%.

Studies using dynamic contrast-enhanced MRI should be accurate and reproducible when compared with a gold standard technique. Several published papers conclude that Gdenhanced MRI provides a reliable estimate of GFR and is suitable for use in both clinical and experimental settings. [START_REF] Hackstein | Contrast media clearance in a single kidney measured on multiphasic helical CT: Results in 50 patients without acute renal disorder[END_REF][START_REF] Hackstein | Measurement of single-kidney glomerular filtration rate using a contrast-enhanced dynamic gradientecho sequence and the Rutland-Patlak plot technique[END_REF][START_REF] Laurent | Effect of essential hypertension on kidney function as measured in rat by dynamic MRI[END_REF][START_REF] Lee | Dynamic three-dimensional MR renography for the measurement of single kidney function: Initial experience[END_REF] However, a recent analysis of the literature showed a great heterogeneity of protocols (e.g., in acquisition mode, dose of contrast, postprocessing techniques) and that these studies presented regression coefficient values that are generally considered inadequate for replacing one technique with an accepted reference method in biomedical sciences. [START_REF] Mendichovszky | How accurate is dynamic contrast-enhanced MRI in the assessment of renal glomerular filtration rate? A critical appraisal[END_REF] The claims in the current literature that using MRI "Gd DTPA clearance provides a good estimate of glomerular filtration [START_REF] Laurent | Effect of essential hypertension on kidney function as measured in rat by dynamic MRI[END_REF] " are not supported and the conclusion that no "approach used provides a wholly satisfactory measure of renal function [START_REF] Buckley | Measurement of single kidney function using dynamic contrast-enhanced MRI: comparison of two models in human subjects[END_REF] " is more accurate. The single publication on reproducibility is not encouraging because the absolute values have a wide % reproducibility, which is not clinically useful. 38 One possible explanation for the "poor" published results on both reproducibility and accuracy is that there is a biological variation in GFR but, obviously, it is unacceptable to introduce methods that have a lower accuracy and reproducibility than existing methods. Only one study undertook the 99m Tc-DTPA clearance during the Gd-MRI [START_REF] Lee | Renal function measurements from MR renography and a simplified multicompartmental model[END_REF] but even these results with a complex mathematical model-based SKGFR the 95% confidence interval varied between 5.8 and 17.9 mL/min compared with the reference technique.

In summary, although MRI has a great potential to provide accurate and reproducible results for measurement of glomerular filtration, none of the proposed methods has been sufficiently validated to be recommended as a routine clinical/ research tool. The high heterogeneity of protocols reflects the number of technical challenges that will have to be solved before to reach a consensus.

Figure 3

 3 Figure 3 Measurement of extraction fraction. (A) Coronal maximum intensity projection image from contrast-enhanced 3-dimensional spoiled gradient-echo MR imaging (3.6/1.0; flip angle, 25°). Lines indicate imaging planes for spiral T1 relaxation measurements. (B) Plot of single-kidney EF determined with inulin measurement and MR imaging. A statistically significant linear correlation between EF measured with MR imaging versus EF IN was found (r 2 ϭ 0.77). (OE) ϭ value observed in kidney without RAS (right kidneys), () ϭ value observed in kidney with RAS (left kidneys). (Reprinted with permission from Coulam et al. 13 )

Figure 4

 4 Figure4SI-time curves obtained on a renal transplant. The aortic curve, for characterization of the arterial input function, is altered by inflow effects only when sampling of SI is done at the top of the volume. The cortical SI time curve can be characterized by 3 phases: a first abrupt ascending segment followed by a first peak, corresponding to the "vascular-to-glomerular first-pass" or cortical vascular phase; a second slowly ascending segment, ended by a second peak, corresponding to the glomerulo-tubular phase; and a slowly descending segment, corresponding to the predominant excretory function and so-called "excretory phase."

  Fig 5 and up to 60% on native kidneys, using Rutland-Patlak method.

Figure 5

 5 Figure 5 Color map of GFR, calculated on a voxel-by-voxel basis, using the Rutland-Patlak plot technique, without (A) and with (B) automatic motion correction.

Figure 6

 6 Figure 6 Color maps of GFR, calculated on a voxel-by-voxel basis, using the Rutland-Patlak plot technique, on the 5 slices of a renal transplant. The cortex was segmented manually.

of GFR with MRI 3

	110
	111
	112
	113
	114
	115
	116
	117
	118
	119
	120
	121
	122
	123
	124
	125
	126
	127
	128
	129
	130
	131
	132
	133
	134
	135
	136
	137
	138
	139
	140
	141
	142
	143
	144
	145
	146
	147
	148
	149
	150
	151
	152
	153
	154
	155
	156
	157
	158
	159
	160
	161
	162
	163
	164
	165
	166
	167

F2

tapraid2/zs1-snuc/zs1-snuc/zs100108/zs12015d08z PossaM Sϭ1 10/9/07 Art: 50172 ce: jm

  -snuc/zs1-snuc/zs100108/zs12015d08z PossaM Sϭ1 10/9/07 Art: 50172 ce: jm tapraid2/zs1-snuc/zs1-snuc/zs100108/zs12015d08z PossaM Sϭ1 10/9/07 Art: 50172 ce: jm

	8	N. Grenier et al
	400 458
	401 459
	402 460
	403 461
	404 462
	405 463
	406 464
	407 465
	408 466
	409 467
	410 411 412 413 414 415 416 417 418 419 420 421 468 469 470 471 472 474 475 476 477 478 479 473 AQ: 7	P R O O F
	422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 480 481 482 483 484 485 486 487 488 489 492 493 494 495 496 497 498 499 500 501 502 503 504 505 AQ: 3 490 491	U N C O R R E C T E D	AQ: 4
	448 506
	449 507
	450 508
	451 509
	452 510
	453 511
	454 512
	455 513
	456 514
	457 515

F6

tapraid2/zs1

tapraid2/zs1-snuc/zs1-snuc/zs100108/zs12015d08z PossaM Sϭ1


Iulius Dragonu

Baudouin Denis De Senneville

Bruno Quesson

Chrit Moonen

Mario Ries
email: mario.ries@imf.u-bordeaux2.fr

Real-Time Geometric Distortion Correction for Interventional Imaging with Echo-Planar Imaging (EPI)

Many MR-guided interventional procedures rely on fast imaging sequences for providing images in real-time with a precise relation between the target position in the image and its true position. Echo-planar imaging (EPI) methods are very fast but prone to geometric distortions. Here, we propose a correction method designed for real-time conditions, adapting existing approaches based on dual EPI acquisition with varying echo times. The method is demonstrated with MR-thermometry for guiding thermal therapies. The proposed approach imposes a small penalty in acquisition speed but adds negligible latency to data processing, an important element for interventions of mobile organs.

In typical interventional procedures such as catheter tracking, MR-guided biopsies, or MR-guided thermal ablation, rapid MR-imaging is typically used to provide functional or positional information in real-time, which is used to retroactively control the interventional procedure. In general, this requires rapid MR-sequences coupled with fast data processing strategies to achieve high frame-rates and low image latencies, allowing guidance of the interventional process and resolution of physiological motion.

Although a large variety of sequences, such as multiecho balanced steady-state free precession (TrueFISP) [START_REF] Deimling | Magnetization prepared true FISP imaging[END_REF] or single-shot echo-planar-imaging (EPI) [START_REF] Mansfield | Multi-planar image formation using NMR spin echoes[END_REF], have been used for this purpose, many rely on EPI readout trains for rapid data acquisition [START_REF] Pruessmann | SENSE: sensitivity encoding for fast MRI[END_REF][START_REF] Sodickson | Simultaneous acquisition of spatial harmonics (SMASH): fast imaging with radiofrequency coil arrays[END_REF]. However, EPI images suffer from geometric distortions due to off-resonance effects resulting from static magnetic field (B 0 ) inhomogeneities, from susceptibility changes between different anatomical regions especially near tissue/air and tissue/bone interfaces and from eddy-current effects. These distortions may lead to positional errors if the images are used for guiding external interventional devices.

Several correction techniques for geometric distortions have been proposed in the past, such as field inhomoge-neity mapping [START_REF] Jezzard | Correction for geometric distortion in echo planar images from B 0 field variations[END_REF][START_REF] Reber | Correction of off resonancerelated distortion in echo-planar imaging using EPI-based field maps[END_REF], point-spread function calculation [START_REF] Zeng | Image distortion correction in EPI: comparison of field mapping with point spread function mapping[END_REF], or dynamically switched phase ramp in kt-space [START_REF] Xiang | Correction for geometric distortion and N/2 ghosting in EPI by phase labeling for additional coordinate encoding (PLACE)[END_REF]. The main limitations of these methods arise either from the inability to rapidly update the distortion information during the acquisition or from the introduced additional computational overhead, which may render them unsuitable for real-time image processing. Another interesting B 0 mapping-based approach was recently proposed by Priest et al. [START_REF] Priest | EPI distortion correction from a simultaneously acquired distortion map using TRAIL[END_REF] using two reduced acquisitions interleaved (TRAIL). This method collects correction data by acquiring dual-echo images with halved image resolution. However, this method requires for the reconstruction the combination of nonsequentially-acquired half-resolution images, which can be a limiting factor in resolving rapid motion. Furthermore, TRAIL suffers from an intrinsic signal loss of ͱ2, which is particularly limiting for interventional applications in regions with short T* 2 such as abdomen or thorax.

We propose a distortion-correction based on field maps, as suggested by Reber et al. [START_REF] Reber | Correction of off resonancerelated distortion in echo-planar imaging using EPI-based field maps[END_REF] in combination with a dynamic update strategy for the correction data and timeefficient data processing, as a simple and efficient method for the real-time correction of geometric distortions in gradient-recalled EPI images. The field maps are temporally updated to compensate for dynamic changes induced by physiological motion or the displacement of interventional devices, by continuously acquiring EPI images with alternating echo-time (TE) values. The proposed approach does not require additional acquisitions, and all the acquired images are used to correct for geometric distortions on the fly.

The proposed method may be particularly useful for MR-guided thermal therapies that induce local hyperthermia using laser, microwave, radio frequency (RF), or highintensity focused ultrasound [START_REF] Kennedy | High intensity focused ultrasound: surgery of the future?[END_REF][START_REF] Levy | Biophysical basis and cardiac lesions caused by different techniques of cardiac arrhythmia ablation[END_REF]. In these applications MRI provides positional and temperature information [START_REF] Hindman | Proton resonance shift of water in the gas and liquid states[END_REF][START_REF] Peters | Ex-vivo tissue type independence in proton-resonance frequency shift MR thermometry[END_REF], and relies on gradient-echo imaging with long TE, often leading to severe distortion artifacts. This study demonstrates the possibility to remove severe image distortions dynamically with MR-thermometry in an ex-vivo RF-heating experiment, and its ability to remove EPI distortions in real-time in an in-vivo experiment imaging on the abdomen of a healthy volunteer under free-breathing conditions.

THEORY

For the case of a typical EPI readout train, the distortion in the read direction is small compared to the distortion in the phase-encoding direction: y pe ϭ y 0pe ϩ ⌬f͑y 0pe ͒ BW pe [START_REF] Deimling | Magnetization prepared true FISP imaging[END_REF] where y pe is the distorted position in the EPI image, y 0pe is the correct position, ⌬f(y 0pe ) is the magnetic field inhomogeneity, and BW pe is the bandwidth per pixel in the phase encode direction.

A field map can be determined by collecting two gradient echo images with different TE values. After calculation of the phase maps from each image followed by phase unwrapping, the field map can be calculated from the phase difference ⌬:

⌬B 0 ϭ ⌬ 2 ⅐ ␥ -⅐ ⌬TE [2]
where ␥ -is the gyromagnetic ratio in Hz/Tesla, and ⌬TE ϭ (TE 2 -TE 1 ) with TE 1,2 representing the alternating echotimes, respectively. Furthermore, the B 0 map can be converted into a pixel shift in the phase-encoding direction (⌬y) according to the following equation:

⌬y ϭ ␥ -⅐ ⌬B 0
BW pe [START_REF] Pruessmann | SENSE: sensitivity encoding for fast MRI[END_REF] Making use of Eqs. [START_REF] Mansfield | Multi-planar image formation using NMR spin echoes[END_REF] and [START_REF] Pruessmann | SENSE: sensitivity encoding for fast MRI[END_REF], one can derive the value of the pixel shift directly from the unwrapped phase difference according to the following relation:

⌬y ϭ ⌬ 2 ⅐ BW pe ⅐ ⌬TE [4]
Two subsequent images are combined in a sliding window to update the field map dynamically for each new acquired image. This leads under the Nyquist criteria to a temporal resolution of 2T acq for the correction, where T acq represents the acquisition time of the EPI image. Note that the modifications of the magnetic susceptibility field induced by organ motion during the acquisition of the two subsequent images have not been taken into account in Eq. [START_REF] Sodickson | Simultaneous acquisition of spatial harmonics (SMASH): fast imaging with radiofrequency coil arrays[END_REF]. Alternatively, a linear interpolation is proposed between the phases ( TE1(n-1) and TE1(nϩ1) ) of two consecutive images with the same TE. This linearly interpolated value gives a more accurate representation of the pixel shift at instant t n according to the following relation:

⌬y ϭ 1 2 ⅐ BW pe ⅐ ͑TE 2 Ϫ TE 1 ͒ ϫ ͩ TE2͑n͒ Ϫ TE1͑nϪ1͒ ϩ TE1͑nϩ1͒ 2
ͪ [5] However, this approach has the inconvenience to introduce a latency of one image to the correction and to reduce the temporal resolution from 2T acq to 3T acq .

For the use of the correction in interventional procedures, it is important to evaluate the upper bound of the precision using Eq. [START_REF] Sodickson | Simultaneous acquisition of spatial harmonics (SMASH): fast imaging with radiofrequency coil arrays[END_REF]:

͑⌬y͒ ϭ 1 2 ⅐ BW pe ⅐ ⌬TE ͱ 2 ͑ TE1 ͒ ϩ 2 ͑ TE2 ͒ [6]
Conturo et al. [START_REF] Conturo | Signal-to-noise in phase angle reconstruction: dynamic range extension using phase reference offsets[END_REF] have demonstrated that the phase error is inversely proportional to the signal-to-noise ratio (SNR) of the complex MR image, i.e. () ϭ (|I|)/|I| ϭ 1/SNR, where I and (|I|) represent the intensity and the noise level of the magnitude image, respectively. Assuming the same noise level in both EPI images and taking into account the T* 2 signal decrease, we derive the following equation for estimating the precision of the displacement value:

͑⌬y͒ ϭ 1 2 ⅐ BW pe ⅐ 1 SNR͑TE short ͒ ⅐ 1 ⌬TE ⅐ ͱ 1 ϩ exp ͩ 2⌬TE
T* 2 ͪ [7] where, SNR(TE short ) represents the SNR of the EPI image acquired with the short echo-time (TE short ).

MATERIALS AND METHODS

Data Acquisition and Processing

All experiments were performed on a 1.5T Philips Achieva clinical scanner using a modified gradient-recalled EPI sequence allowing the acquisition of dynamic images with alternating TE values. In a first processing step, the pixelby-pixel phase difference of two subsequent images is calculated and a temporal phase unwrapping step is performed to remove 2 aliasing. Then, the displacement of each pixel in the phase encoding direction is directly calculated using Eq. [START_REF] Jezzard | Correction for geometric distortion in echo planar images from B 0 field variations[END_REF]. Finally, this information is used to unwrap the current EPI image. The distorted EPI images can be considered as irregularly sampled data, which is mapped on a regularly-sampled Cartesian grid. Adjacent pixels on the irregular grid were remapped and linearly interpolated to obtain the pixel intensities on the Cartesian grid. To avoid interpolation between regions with very low signal intensity and object boundaries, an adaptive thresholding was applied previous to this step. For the in-vivo experiments with low SNR, an additional filtering step was introduced: the thresholding mask was obtained from a filtered version of the magnitude image (Gaussian filter, 3-pixel kernel diameter) and the distortion map was calculated from phase data derived from smoothed complex data (Gaussian filter, 3-pixel kernel diameter), whereas all geometric corrections were applied to the unfiltered complex data.

The unwrapping algorithm was applied to the complex pixel values. This allows undistorted magnitude and phase images to be obtained, which can be used for MRthermometry using the proton resonance frequency (PRF) shift technique [START_REF] Hindman | Proton resonance shift of water in the gas and liquid states[END_REF][START_REF] Peters | Ex-vivo tissue type independence in proton-resonance frequency shift MR thermometry[END_REF].

Heating Study

A fresh ex-vivo pig kidney was installed on a motorized platform performing a periodical linear displacement of 10 mm in 6.3 s. Heating was performed with two copper electrodes using a clinical RF device (Tyco Healthcare/ Radionics, Burlington, MA, USA) operating at 480 kHz frequency with 100 W of maximum output power. RF power at 10 W was applied 15 s from the start of the image sequence with a total duration of 22 s. MRI was performed using the following single-shot EPI sequence: 1024 dynamic coronal images, TR ϭ 57 ms, TE ϭ 27 ms, ⌬TE ϭ 3 ms, flip angle ϭ 25°, FOV ϭ 240 ϫ 128 mm 2 , slice thickness ϭ 6 mm, matrix size ϭ 112 ϫ 60, BW pe ϭ 22.4 Hz, and BW rd ϭ 2593.4 Hz. A T 2 -weighted turbo spin-echo (TSE) sequence with the same resolution was used for comparison.

Motion-compensated MR-thermometry was performed as presented [START_REF] De Senneville | Real time adaptive methods for treatment of mobile organs by MRI controlled high intensity focused ultrasound[END_REF]. This PRF-based method acquires complete set of reference magnitude and phase images for each organ position during a pretreatment step. During the intervention, the current magnitude image is compared to the prerecorded set and the closest match selected by evaluating the intercorrelation coefficient. The corresponding prerecorded phase image serves as reference for current temperature computation. This strategy allows correcting online periodical interscan motion-related errors in PRF temperature methods. To be compatible with the proposed distortion correction, the two separate lookup tables (N ϭ 100) for each of the two TE values were used.

Real-Time In-Vivo Study

Dynamic MRI was performed under free-breathing conditions of the abdomen of a healthy volunteer and reconstructed in real-time. The MR sequence employed the following parameters: 500 dynamic sagittal single-shot EPI images, TR ϭ 78 ms, TE ϭ 39 ms, ⌬TE ϭ 7 ms, FOV ϭ 350 ϫ 227.5 mm 2 , slice thickness ϭ 7 mm, matrix size ϭ 128 ϫ 84, BW pe ϭ 14.1 Hz, and BW rd ϭ 2435 Hz. To avoid chemical shift artifacts in the field map due to fat resonance frequency, a 121 spectral-spatial [START_REF] Meyer | Simultaneous spatial and spectral selective excitation[END_REF] water selective pulse was used. A four-channel abdominal phased-array coil was used for signal detection. The experiment was repeated with a fixed echo-time in order to evaluate the influence of modifications of the magnetic susceptibility field induced by organ motion during the acquisition of two subsequent images. This correction error can be calculated by inserting the phase differences into Eq. [START_REF] Pruessmann | SENSE: sensitivity encoding for fast MRI[END_REF].

k-Space data was transferred directly from the MR acquisition system to an external in-house-developed realtime reconstructor using The Ace Orb (TAO) a real-time implementation of the Common Object Request Broker Architecture (CORBA) communication protocol [START_REF] Pyrali | Techniques for enhancing real-time CORBA quality of service[END_REF]. The entire reconstruction was implemented as a multithreaded processing chain in Cϩϩ/OpenMP on a four-processor AMD Opteron (2.4 GHz). The bulk reconstruction time per image was 4 ms and the image latency was 55 ms post-TE.

A steady-state free precession (SSFP) image with the same image geometry and a two-fold increased image resolution was used for the evaluation of the efficacy of the correction (TR ϭ 3.3 ms, TE ϭ 1.66 ms).

RESULTS

Theoretical Calculations

Figure 1a shows the theoretical precision of the correction given by Eq. [START_REF] Zeng | Image distortion correction in EPI: comparison of field mapping with point spread function mapping[END_REF], evaluated for the typical case of an EPI sequence, as a function of ⌬TE, in pixels per unitary bandwidth. The function was evaluated for two SNR values (10 and 15), which typically represents the lower limit for MR-thermometry or interventional guidance and for the T* 2 values for liver and kidney tissue at 1.5T [START_REF] De Senneville | MR thermometry for monitoring tumor ablation[END_REF]. The minima of this function can be found for a ⌬TE which is equal to 1.108 times the T* 2 of the observed tissue. However, in practice a shorter ⌬TE is preferable for the correction to avoid phase aliasing due to strong local field deviation and to limit the time penalty on the acquisition speed of the sequence. Assuming a precision of less than one-third of a pixel as the acceptance criteria for the correction, and considering the most unfavorable case: a SNR value of 10, the lowest EPI phase-encoding bandwidth of the presented experiments of 14 Hz and the T* 2 value of 15 ms corresponding to the liver tissue, a lower bound for ⌬TE of 6 ms can be estimated as sufficient.

Strong local inhomogeneities can lead to phase differences larger than 2 between subsequent phase maps and thus an incorrect displacement vector is calculated. To avoid this effect, ⌬TE has to be chosen according to the largest field deviation that is present in the area of interest with the help of Eq. [START_REF] Mansfield | Multi-planar image formation using NMR spin echoes[END_REF], as shown in Fig. 1b. For ⌬TE ϭ 5 ms to 7 ms the maximum inhomogeneity that can be corrected is 3.1 ppm to 2.2 ppm at 1.5T. For a typical BW pe ϭ 22 Hz and ⌬TE ϭ 5 ms, distortions up to 9.1 pixels (i.e., Ϯ4.55 pixel from the undistorted position) can be corrected. Considering the voxel size for the phantom study (1.72 ϫ 1.72 ϫ 5 mm 3 ) the proposed method will correct distortions in the phase-encoding direction up to Ϯ7.82 mm with a precision of Ϯ0.57 mm.

Heating Study

The accuracy of the correction is shown in Fig. 2. An edge detection algorithm was used to highlight the borders of the pig kidney in the TSE image (Fig. 2a). This information was overlaid on the distorted image (Fig. 2b) and the corrected EPI datasets. B 0 inhomogeneities Ͼ2.0 ppm were found (Fig. 2d), leading to shifts of up to 5.7 pixels. Figure 2e andf show the standard deviation (SD) of the temperature maps before and after the unwrapping algorithm. The high-intensity spots on both images correspond to the heated region in the vicinity of the copper electrodes. Low values of SD (Ͻ0.9°C) of temperature maps in the nonheated regions show that the precision of the MR-thermometry is not altered by the spatial remapping.

The advantage of a geometry correction for multimodal imaging is illustrated by Fig. 2g andh, which show the single-voxel temperature evolutions obtained by MR thermometry compared to the readings of a fiber optic probe. Since the position of the fiber optic probe was determined by an undistorted high-resolution spin-echo (SE) image (at the position indicated by the arrow in Fig. 2a), the MR thermometry based on the uncorrected EPI shows a temperature difference due to a geometry mismatch. The readings obtained from the distortion-corrected EPI images taken at the probe position lead to a temperature evolution that closely matches that of the readings of the fiber optic probe. Note that the readings remain a good match throughout the entire motion cycle despite the fact that the distortions at the position of the probe vary temporally by 0.5 ppm.

Real-Time In-Vivo Study

Figure 3b shows the amplitude of the distortions of the uncorrected sagittal EPI image displayed in Fig. 3c extracted from the dynamic acquisition. The map shows large local fluctuations especially in the vicinity of tissue boundaries. While the distortions present within the kidney and the back muscle are below the resolution of the EPI images for this experiment, distortions in the phaseencoding direction (in this case anterior-posterior) within the liver exceed 1 cm. Since the lung volume and the location of the organs has a direct impact on the local magnetic field homogeneity, the distortion field changes dynamically throughout the respiratory cycle. In the vicinity of the area indicated by arrow II in Fig. 3c the amplitude of the distortion oscillates by up to 8 mm during respiration. The overlaid anatomical reference points of the uncorrected image (Fig. 3c) and those of the corrected image (Fig. 3d) demonstrate the efficacy of the correction.

The error induced by field modifications during the acquisition of two subsequent images for a pixel within the liver is shown for direct phase subtraction and phase averaging in Fig. 3e andf, respectively. Direct phase subtraction results in a peak-to-peak error of 1.9 pixels and an SD of 0.35 pixels. Note that the error is two-fold periodic: It is modulated with the respiratory frequency and oscillates in average by 0.4 pixels between two subsequent images. Phase averaging removes this effect and reduces the error to 1.1 pixels peak-to-peak and an SD of 0.2 pixels.

DISCUSSION

For typical gradient-recalled EPI images of abdominal organs such as liver and kidney, with typical SNR Ͼ 15, the simulations show that the proposed distortion correction results in a precision of the correction within 0.16 pixels while introducing an echo-time penalty of 5 to 7 ms. The small echo-time penalty and the high precision correspond well with the real-time requirements of MR-guided interventional procedures.

The choice of ⌬TE of 5 ms in the proposed implementation restricts the correction to areas with field variations FIG. 3. a: A high-resolution SSFP image with a short TE was used as a reference image. b: A map of the distortions in millimeters (mm) of the uncorrected EPI image displayed in (c), which is an example extracted from the temporal series. This distortion map is dynamically updated for each acquired image. The kidney shows moderate distortions of ϳ1.5 mm in this example (below the voxel size), whereas image distortions in the liver can be larger than 1 cm, especially in the vicinity of the diaphragm. The main anatomy of the liver and the kidney obtained from the SSFP image was transposed to the distorted EPI image (c) and the distortion-corrected EPI image (d) for comparison. Note the mismatch between the mask and the position of the vessels in the center of the liver and at the interface lung/liver (see arrow), which are removed by the distortion correction. The correction error at position indicated by arrow II for the direct phase subtraction (e) and for the phase-averaged technique (f).

AQ: 4 F3

tapraid5/z3t-mrm/z3t-mrm/z3t00209/z3t4334d09g gockleyj Sϭ7 11/25/08 9:00 Art: 08-9929 less than Ϯ1.56 ppm as shown in Fig. 1b. Areas with more important inhomogeneities therefore show an aliasing artifact of the correction, and thus need to be masked. For a typical single-shot EPI sequence with a BW pe of 14 -22 Hz and a resolution of 2 mm. this corresponds to shifts of Ϯ1.4 cm to Ϯ0.9 cm. Since the proposed method is designed for real-time applications requiring low image latency and thus low computational overhead, this parameter choice represents a design compromise. This limitation can be circumvented by employing a 2D phase-unwrapping step to the phase images prior to the calculation of the field maps [START_REF] Dias | InSAR phase unwrapping: a Bayesian approach[END_REF] as is required for TRAIL [START_REF] Priest | EPI distortion correction from a simultaneously acquired distortion map using TRAIL[END_REF]. However, this remains a computationally-intense processing step, which is often unstable in areas with signal discontinuities and with the strong susceptibility changes frequently encountered in abdominal imaging. In addition, these areas with strong susceptibility changes are subject to rapid intravoxel dephasing, which in turn leads to a low local SNR (Ӷ10) for gradient-recalled imaging and are thus of limited value for MR-thermometry or MR-guidance [START_REF] De Senneville | MR thermometry for monitoring tumor ablation[END_REF].

Another limitation of the proposed correction method arises from the fact that the field maps are calculated using a sliding-window technique. When direct phase subtraction of subsequent images is employed during evolution of the field map, slight position and thus phase changes lead to small errors in the correction map. It can be shown with Eq. [START_REF] Sodickson | Simultaneous acquisition of spatial harmonics (SMASH): fast imaging with radiofrequency coil arrays[END_REF] that the sign of this residual error changes with every new correction map update. Consequently, the relative positioning error between two subsequently corrected images is effectively doubled. The in-vivo data reveals that this problem is in practice frequently encountered on fast-moving objects such as liver and kidney, as shown in Fig. 3e. The resulting dynamic correction artifact can be effectively removed by widening of the sliding window to three images in conjunction with phase averaging as shown in Fig. 3f. However, this reduces the temporal resolution of the correction to one-third of the bandwidth of the employed sequence and introduces a temporal latency of 1 TR. Consequently, this type of correction is only suitable in conjunction with very fast sequences that observe the physiological motion of the target area with a very high temporal resolution.

The heating experiments demonstrate the compatibility of the proposed distortion correction with MR-thermometry during thermoablation. Furthermore, the distortions in the EPI images, introduced by the susceptibility difference between the RF electrodes and the surrounding tissue, are removed. In addition, the slight local frequency changes introduced by temperature variations in the observed tissue theoretically also lead to local image distortions of the EPI phase-encoded images. In the experiment presented, the total temperature difference in the vicinity of the RF needles was 28°C. The resulting frequency shift can be calculated to -0.28 ppm [START_REF] Hindman | Proton resonance shift of water in the gas and liquid states[END_REF]. Using Eq. [START_REF] Sodickson | Simultaneous acquisition of spatial harmonics (SMASH): fast imaging with radiofrequency coil arrays[END_REF] with a BW pe ϭ 22.4 Hz, this leads to a shift of -0.81 pixel, which is superimposed on the distortions caused by the intrinsic magnetic field inhomogeneities. Since the proposed distortion correction dynamically updates the correction maps with each new arriving image, the spatial mismatch due to the temperature evolution can also be removed.

The in-vivo experiments presented demonstrate the ability to correct spatial distortions of EPI sequences for both organs in real-time with a dynamic correction scheme and allows a precise spatial match of thermometry with undistorted anatomical images. The absence of spatial fitting of the field map to polynomial or other model functions is advantageous for a real-time corrections, since this renders the method insensitive to signal discontinuities in regions where field maps cannot be determined, represents the local fluctuations well along organ boundaries compared to a globally-fitted function [START_REF] Jezzard | Correction for geometric distortion in echo planar images from B 0 field variations[END_REF], and finally limits the required processing steps and thus processing time. The resulting low computational overhead leads to complete execution time of 1.5 ms for the entire correction (including the field calculation and the regridding of the complex data) and thus to reconstruction times below 5 ms, making the method compatible with rapid acquisition sequences of more than 20 images per second. However, since the method collects correction data with one-half the temporal resolution of the base sequence, high frame rates are required to sample the changes during the motion cycle with sufficient resolution.

It should be noted that, if areas with large B 0 inhomogeneities lead to larger phase difference than 2 during the chosen ⌬TE, the correction method will not only fail, but can also lead to increased distortions (see arrow III in Fig. 3d). It is thus indispensable that an initial experiment with a short ⌬TE establishes a validity map for the subsequent correction with longer ⌬TE. In case the region of interest is not included in this validity map, either ⌬TE has to be shortened accordingly with concomitant lower precision of the correction or 2D phase unwrapping has to be employed.

CONCLUSIONS

The proposed method is an extension of existing approaches based on dual-TE gradient-recalled images with long readout trains, and especially adapted for the use under real-time conditions for interventional imaging on moving targets. It is compatible with MR thermometry and interventional devices such as needles and catheters, which typically cause severe local image distortions.

The method proposed imposes a small speed penalty on the acquisition sequence but adds a negligible latency to the data processing, which is of great importance for realtime imaging. As a drawback, the method introduces a fluctuating signal magnitude and is therefore not suitable for applications that require the magnetization to remain in a steady-state, such as the observation of bolus passages or functional MRI (fMRI) experiments.

FIG. 1 .

 1 FIG. 1. a: Theoretical values of the uncertainty of the displacement vectors normalized to an acquisition bandwidth in the phase encoding direction of 1 Hz as a function of ⌬TE. The error of the displacement can be computed by dividing the vertical axis number from the graph by the bandwidth in the phase-encoding direction. The error is shown for two different T* 2 values and two different SNRs. The solid line and dotted line represent the correction uncertainty for liver tissue (T* 2 ϭ 15 ms at 1.5T) for an SNR of 10 (solid) and 15 (dotted), respectively. The dashed line and the dashed-dotted line correspond to the values typically found in kidney tissue (T* 2 ϭ 25 ms at 1.5T) for an SNR of 10 (dashed) and 15 (dash-dotted), respectively. b: Upper limit for ⌬TE before phase aliasing occurs as a function of the magnetic field inhomogeneity for a B 0 field of 1.5T in a given region of interest.

FIG. 2 .

 2 FIG. 2. Demonstration of the effectiveness of the proposed method for the ex-vivo heating study. a: TSE image used as a reference. b: Uncorrected EPI image. c: Corrected EPI image. d: The field inhomogeneity conditions of the experiment. Standard deviation of the distorted (e) and corrected (f) temperature maps during the heating experiment. Comparison between the PRF-based MR thermometry (dashed line) and the reference temperature measured with the optical fiber thermometer (solid line) at the position determined on the TSE image (black arrow in (a)) without (g) and with (h) the proposed distortion correction method.

tapraid5/z3t-mrm/z3t-mrm/z3t00209/z3t4334d09g gockleyj Sϭ7 11/25/08 9:00 Art: 08-9929

Grant number: EC-FP6-project LSHB-