Quantitative Analysis of Renal Perfusion by Arterial Spin Labeling

The signal intensity differences measured by an arterial-spin-labelling (ASL) magnetic resonance imaging (MRI) experiment are proportional to the local perfusion, which can be quantified with kinetic modeling. Here we present a step-by-step tutorial for the data post-processing needed to calculate an ASL perfusion map. The process of developing an analysis software is described with the essential program code, which involves nonlinear fitting a tracer kinetic model to the ASL data. Key parameters for the quantification are the arterial transit time (ATT), which is the time the labeled blood takes to flow from the labeling area to the tissue, and the tissue T1. As ATT varies with vasculature, physiology, anesthesia and pathology, it is recommended to measure it using multiple delay times. The tutorial explains how to analyze ASL data with multiple delay times and a T1 map for quantification. This chapter is based upon work from the COST Action PARENCHIMA, a community-driven network funded by the European Cooperation in Science and Technology (COST) program of the European Union, which aims to improve the reproducibility and standardization of renal MRI biomarkers. This analysis protocol chapter is complemented by two separate chapters describing the basic concept and experimental procedure.


Introduction
The signal intensity differences measured by an arterial-spin-labelling (ASL) magnetic resonance imaging (MRI) experiment are proportional to the local perfusion, which can be quantified with kinetic modeling. To calculate the perfusion map, an additional T 1 relaxation map should be obtained from inversion recovery data acquired with the same image acquisition method (eg, Spin-Echo EPI or RARE) and geometry setting.
Creation of perfusion maps involves nonlinear fitting of a tracer kinetic model to the ASL data. Depending on the type of ASL sequence-pulsed or (pseudo)continuous-the model is slightly different. Two key parameters in the quantification are the arterial transit time (ATT), which is the time the labeled blood taken to flow from the labeling area to the tissue, and the tissue T 1 . As ATT varies with vasculature, physiology, anesthesia, and pathology, it is recommended to measure it using multiple delay times (i.e., multiple TI in pulsed ASL, or multiple postlabeling delays in continuous ASL). In the kidney, the cortex and medulla have different T 1 . T 1 could also change depending on the pathological conditions. Therefore acquiring a T 1 map is desirable. This chapter focuses on how to analyze ASL data with multiple delay times and a T 1 map for quantification.
This analysis protocol chapter is complemented by two separate chapters describing the basic concept and experimental procedure, which are part of this book.
This chapter is part of the book Pohlmann A, Niendorf T (eds) (2020) Preclinical MRI of the Kidney-Methods and Protocols. Springer, New York.

Software Requirements
The method described in this chapter provides a detailed description for a solution in MATLAB, but can be adopted to other platforms: If the ASL data of multiple delay times was acquired by separate scans (eg, one delay time for one scan) but not all within the same scan, there could be a difference in the intensity scaling among separate scans due to changes in receiver gain and/or internal scaling factor during image reconstruction. Such a difference could lead to artifactual intensity among delay times leading to bias in the model fitting and inaccurate perfusion quantification. On Bruker MRI systems, the default intensity scaling is set to maximize the dynamic range of the output format. Therefore scans with lower signal would be magnified. The intensity scaling method in the "Reconstruction" class should be changed to "User Scaling" with a scaling factor of 1 for all the scans.
To ensure the same intensity scaling of the ASL data, the receiver gain may be extracted from the header information ("RG" variable in the "acqp" file) and used to rescale the data.

Format Conversion
Convert all data into 4D NifTi format. Consider magnifying the voxel size by 10 in the image header if motion correction tools designed for humans (eg, FSL or SPM) will be used (see Note 1).

Quality Control/ Data Exclusion
During the ASL image acquisition, movement of the body or kidney itself can occur. Artifacts may also arise from ASL labeling pulses. In order to construct accurate perfusion maps, it is important to ensure that the scan series don't show severe movement, spikes, banding artifacts, or sudden changes in SNR. Scans with poor quality or large movement should be discarded.
1. Open each dataset by a NifTi image viewer.
2. Adjust the window/level (or reduce the maximum intensity range) to visualize the image better.
3. Scroll through each time frame to visually check whether there are sudden changes in intensity or movement.

Motion Correction
ASL perfusion imaging, due to its subtraction between label and control images, is very sensitive to movement. Despite respiratory trigger, slight movement between scans could still be present. Motion artifact maybe reduced by rigid-body motion correction tools in many software packages. Since most motion correction algorithms (regardless correlation or square error based cost function) rely on intensity changes for detecting movement, the intensity changes between control and label images and between scans of different delay times (particularly for FAIR ASL) could lead to pseudo motion. Besides, the T 1 mapping data should also be coregistered to the ASL data. Therefore, it is recommended to realign images to the same reference target, (e.g., the M 0 image or an averaged image). Additionally, use "cost functions" that are less dependent on the global intensity difference, such as normalized correlation or mutual information, and avoid using least square error (see Note 2).
1. Use FSL mcflirt to do motion correction on the ASL data (ASL_DATA_TI1.nii) using the mean image as the reference target, normalized correlation as the cost function, and spline for interpolation: mcflirt -in ASL_DATA_TI1.nii -out rASL_DATA_TI1 -cost normcorr -meanvol -spline_final.
3. Repeat the same procedures until all the ASL data are coregistered.
3.3 Quantification of M 0 and T 1

Model Equations
Both M 0 and T 1 are important parameters needed in ASL kinetic model. Besides, the inversion efficiency will be used for FAIR ASL quantification to correct its "labeling efficiency." The most common method to obtain these parameters is to fit the model of an inversion recovery experiment to the signal intensity data of each pixel at various TI, M(TI), using the following equation: where M 0 is the fully relaxed magnetization signal which is a scaling factor that includes many parameters such as the proton density together with the coil sensitivity and the signal gain of the system. α is the inversion efficiency of the inversion pulse (α ¼ 1 for a perfect inversion), and T 1 is the longitudinal relaxation time.

Starting Values
An important step of the curve fitting process is the choice of suitable starting values for each parameter that will be determined by the fitting algorithm (NB: different starting values may lead to different results!). Here we describe how to derive starting values from the SI of each pixel.

1.
Step through all pixels in the images using for loops for the pixel coordinates x and y.
3. As starting value for the parameter M0 use the largest SI (in Matlab: startVal_M0 ¼ max(SI_vector)).
4. For estimating a starting value for T1, one can exploit the fact that the null point of an inversion recovery, where the SI is 0, is TI ¼ ln2 Â T 1 ¼ 0.693 Â T 1 . Estimate the starting T 1 value by finding the TI whose SI is minimal: 5. Find the index of the smallest value in the vector (in Matlab: [min_Value,min_Index] ¼ min(SI_vector)).
6. Get the respective TI using the index from the last step and divide by ln2 to estimate T 1 (in Matlab: startVal_T1 ¼ TI_for-Fitting(min_Index)/0.693).
As for a starting value for α, one can assume that is close to 1 (in Matlab: startVal_alpha ¼ 1).

Fitting Algorithm
Least squares algorithms, such as the Levenberg-Marquardt and Trust-region methods, are the most commonly used curve fitting algorithms for T 1 -mapping. They work by minimizing a cost function, which describes the deviation of the fitted curve from the data points. With starting values near the optimal solution they quickly converge, but with starting values far away from the solution, the Levenberg-Marquardt algorithm will slow down significantly. Also, there is the risk that it may converge to a local minimum (rather than the global minimum) and hence produce a wrong result.
In contrast, the Trust-region method (a further development of the Levenberg-Marquardt algorithm) will quickly converge, even with suboptimal starting values, and it will always find the global minimum. However, it does require the definition of lower and upper limits for the parameters to be fitted. Hence, for T 1 -mapping, where such limits can easily be defined, the Trustregion method is well suited (see Note 3).
1. Create for loops for pixel coordinates x and y to step through all pixels of the image.

Model Equation for FAIR-ASL
The pairwise subtracted perfusion-weighted signals at different TIs, ΔM(TI), can be fitted to the following kinetic function [1] (see Fig. 1): where 1/T 1app ¼ 1/T 1 + f/λ, f is the quantified perfusion in ml/100 g/min, λ is the blood-tissue partition coefficient. In the above equation, M 0 represents the equilibrium magnetization derived from T 1 mapping, T 1 the tissue longitudinal relaxation time, and α (the inversion efficiency) are all determined from the curve fitting of the inversion recovery T 1 mapping data. λ is a value typically derived from literature which ranges between 0.52 and 0.94 have been reported in [2]. Here, λ ¼ 0.9 ml/g can be used based on studies in the brain [3]. However, the suitable value for the kidney remains to be determined. T 1a is the longitudinal relaxation time of arterial blood, which is field strength dependent. A value of 2.21 s can be used at 7 T [4]. And the blood T 1 at other field strength can be calculated using the equation in [4] (see Note 5).

Pairwise Subtraction
The first step is to generate the ΔM image of each TI, and then concatenate the results of different TI together to form a 4D timeseries data: 1. Use the following FSL command to do pairwise subtraction: asl_file --data = rASL_DATA_TI1 --ntis = 1 --iaf = tc --diff --out = rASL_DeltaM_TI1 --mean = rASL_DeltaM_TI1_mean Note that the --iaf option is to specify the order of labelcontrol pair where tc means the label (tag) image is acquired first and then followed by control image. Use --iaf ¼ ct if the data is acquired with the opposite order. where n is the number of TI.
3. If the data is acquired with varying TIs in one scan, the command could be: asl_file --data = rASL_DATA --ntis = n --iaf = tc --diff --out = rASL_DeltaM --mean = rASL_DeltaM_mean The resulted rASL_DeltaM_mean image will then be used in the following fitting process.

Visual Display
1. Display the parameter map, which is a matrix with floating point numbers, as an image (in Matlab: imagesc(mapF);) (see Fig. 2).
2. Remove axis labels and ensure that the axes are scaled equally, so that pixels are square and not rectangular (in Matlab: axis off; axis equal;).

Regional Analysis
To obtain quantitative medulla and cortex perfusion, the kidney can be manually or automatically segmented based on the T 1 map or other kinds of structural images.

Comparison with Tissue Values Under Hypercapnia or Hyperoxia
Further validation of imaging derived results can be performed using datasets acquired under physiological manipulations that alter blood flow, such as hypercapnia or hyperoxia. As CO 2 is a potent vasodilator, the blood flow under hypercapnia condition should be higher than that in normal gas condition. Hyperoxia is expected to lead to the opposite result with O 2 acting like a vasoconstrictor (see Fig. 2).

Comparison with Reference Values from the Literature
Obtained values for both cortex and medulla tissue segmentation can be compared against reference values, shown in Table 1. It should be noted that blood flow is also modulated by anesthesia. For example, the commonly used isoflurane is a vasodilator and can increase blood flow depending on dosage. Therefore, the 3. An approach that is frequently used to increase the speed of the curve fitting, is to fit a linear model to the log of the SI data. This only requires a few simple modifications to the protocol described here. We recommend using exponential fitting (unless speed is very important), because the log-scaling leads to the noise error at each TI being scaled differently and this will impact on the fitted curve. If you do want to use linear fitting it would be good to initially perform both, exponential and linear fitting, on a few data sets and compare the results so as to better appreciate the impact of this choice.
4. Pseudocolor representations can be extremely useful for analyzing T 1 -maps or perfusion maps, since they generally enhance  [5] 550-750 140-230 C57BL/6J mouse, 1.5% isoflurane (gas unknown) [6] 344-575 ND C57BL/6N mouse, isoflurane (dose and gas unknown) [7] 288.4 AE 51.3 ND Spraque-Dawley rat, 1.5% isoflurane (gas unknown) [8] 750 AE 80 ND Fisher 344 rat, thiobutabarbitural, 69-135 mg/kg (gas unknown) [9] the perception of differences within the value range of the parameter. But one needs to be careful when choosing a color map. Parameter maps displayed as images in pseudocolor can potentially be misleading. Small differences in the underlying values may be artificially emphasized by a change in color hue, or a significant parameter difference that can easily be seen in a gray-scale image may be flattened or hidden if there is little change of hue or brightness over a certain range of the colorscale. It is recommended to always use the same color map and scale to improve comparability.
5. The transit time is ignored in the above equation as the delivery of spins is almost instantaneous with a small gap between the imaging and inversion slices used in FAIR ASL. Whether the transit time is negligible can be verified by adding the transit time into the kinetic model. 6. To get the correct unit, the T 1 , T 1a and T 1app should be expressed in "second" and the perfusion value, f, will need to be multiplied by 6000 to convert to per 100 g/min.