Analysis Protocols for MRI Mapping of the Blood Oxygenation–Sensitive Parameters T2* and T2 in the Kidney

Renal hypoxia is generally accepted as a key pathophysiologic event in acute kidney injury of various origins and has also been suggested to play a role in the development of chronic kidney disease. Here we describe step-by-step data analysis protocols for MRI monitoring of renal oxygenation in rodents via the deoxyhemoglobin concentration sensitive MR parameters T2* and T2—a contrast mechanism known as the blood oxygenation level dependent (BOLD) effect. This chapter describes how to use the analysis tools provided by vendors of animal and clinical MR systems, as well as how to develop an analysis software. Aspects covered are: data quality checks, data exclusion, model fitting, fitting algorithm, starting values, effects of multiecho imaging, and result validation. This chapter is based upon work from the PARENCHIMA COST Action, 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 experimental protocol chapter is complemented by two separate chapters describing the basic concept and data analysis.


Introduction
The parametric mapping of the transverse relaxation times T 2 * and T 2 (or relaxation rates R 2 * ¼ 1/T 2 * and R 2 ¼ 1/T 2 ) has the potential to yield inferences regarding renal oxygenation, since both parameters are sensitive to blood oxygenation [1]. The underlying mechanism rests on the inherent difference in the magnetic properties of oxygenated hemoglobin (diamagnetic) vs. deoxygenated hemoglobin (paramagnetic) [2]. The presence of deoxyhemoglobin in a voxel decreases both relaxation times, T 2 * and T 2, which are the time constants governing the exponential signal decays due to spin dephasing in gradient-echo (GRE) and spin-echo (SE) MR measurements, respectively.
Possible sources for this dephasing are magnetic field inhomogeneities ranging from the microscopic to the macroscopic scale, which can be classified with respect to the echo time into microscopic dynamic spin-spin interactions-T 2 , that is, typically ranging a time span between 1 and 100 ms (Fig. 1) and static mesoscopic interactions-T 2 0 . In fact, T 2 * includes the dynamic (irreversible) dephasing effects described by T 2 plus the additional effects that are due to static (reversible) dephasing effects described by T 2 0 . Blood oxygenation affects primarily T 2 * (often referred to as Blood Oxygenation Level Dependent, BOLD), but to a lesser extend also T 2 , via water diffusion effects in the proximity of blood vessels due to local magnetic field gradients [1].
Calculation of T 2 * and T 2 requires a series of MR images with different echo times. Repeated measurements with a single echo Fig. 1 Spin dephasing due to dynamic (irreversible; represented by the parameter T 2 ) and static (reversible; represented by the parameter T 2 0 ) sources can be quantified by calculation of T 2 * and T 2 from series of MRI images with different echo times using 1/T 2 * ¼ 1/T 2 + 1/T 2 0 . The possible sources for this dephasing range from the microscopic to the macroscopic scale and their level dictates the observed attenuation in the signal intensity. Blood oxygenation affects primarily T 2 * (often referred to as Blood Oxygenation Level Dependent, BOLD), but to a lesser extend also T 2 via diffusion effects in the proximity of blood vessels. USPIO: Ultra-small Superparamagnetic Particles of Iron Oxide, which can be used as intravascular contrast agents sequence using either a GRE or SE method and variable echo times (TE) would provide the most accurate T 2 * and T 2 values but would require very long acquisition times [2]. In order to address this shortcoming fast multiecho MRI methods are commonly used for in vivo studies: multi-gradient-echo (MGE) for T 2 * and multi-spinecho (also called "MSME") for T 2 . They provide relaxation times that slightly underestimate T 2 * and overestimate T 2 (see Note 1).
The core of T 2 ( * ) mapping consists of fitting the exponential model curve S(TE) ¼ S 0 exp.(ÀTE/T 2 ( * ) ) to the signal intensities of each image pixel at increasing TE. For the calculation of such parameter maps from a series of T 2 ( * ) -weighted images software tools and plugins are provided by the MR system vendors but they may lack some features instrumental for precise T 2 ( * ) mapping. For example, some important preprocessing steps may not be available, such as eliminating the first echo of T 2 data, removing echoes acquired at high TEs that may have too low signal-to-noise ratio (SNR) or applying Rician noise bias correction. Also, details about the used processing and curve fitting algorithms may not be available. It is also a limitation that the algorithms are fixed and cannot be modified for specific purposes. Therefore, depending on application, it may be advantageous to develop a dedicated analysis software program in-house. The analysis software is developed assuming 2D data was acquired, for 3D data (multiple slices) data set are separated into individual slices.
The pixel-wise signal fitting with a monoexponential model is known to introduce T 2 bias, and so for studies that require an extremely rigorous T 2 value (error < 1%) a dictionary-based method has been suggested. In dictionary-matching methodologies, T 2 maps are calculated by matching the T 2 signal decay to precomputed Echo Modulation Curve (EMC), therefore accounting for all echo pathways. The use of dictionary-based methods has been suggested to improve T 2 accuracy, accounting for the stimulated echoes and effective B 1 + field [3]. The use of dictionary-based methods is still fairly novel and its implementation and use is not covered in this chapter. This chapter will describe both how to use the vendor's analysis tools and how to develop your own analysis software conceptually. An example of a MATLAB script described in this chapter together with a data sample (T2analysis.m and data.mat) can be downloaded from https://github.com/JoaoPeriquito/T2-s-Mapping.
This data analysis protocol is complemented by two separate chapters describing the basic concept and experimental, 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.

Optional Tools
An image processing software for checking data quality, such as ImageJ (free Java-based image processing program developed at the National Institutes of Health and the University of Wisconsin; https://imagej.net/). We recommend and describe here the use of Fiji (https://fiji.sc/ [4]), which is ImageJ packaged with a wide range of plugins included.

Source Data: Format Requirements and Quality Check
Before the analysis it is highly recommended to check the image quality. This check should include SNR measurements, particularly for the images with longest echo times, and the assessment of geometric image distortions, motion artifacts, or susceptibility artifacts.
The steps in this section can be performed either on the scanner console using the MR vendors system viewing software, or offline using a software such as Fiji (recommended as a practical tool).

Input Requirements
As already mentioned, in order to be able to calculate T 2 ( * ) -maps, multiple images acquired with different echo times are needed. As an example, please refer to Fig. 2 and Fig. 3, where a series of T 2 *and T 2 -weighted images suitable for mapping is presented (see Note 2). Access to specific acquisition parameters is also necessary, namely the TE of each image and in some cases the image intensity scaling parameters (offset, slope) used when storing the image data.

Open/Import Images
When using the scanner console: 1. Open the T 2 ( * ) -weighted image series in the image viewer.
When using Fiji for DICOM images: 2. Import the DICOM image series using the DICOM import plugin.
When using Fiji for images in Bruker format: 4. Open the "method" file of the scan in a text editor or by dragging it to the Fiji window.
6. Load the image series using the RAW import function (Import > RAW) and providing the above parameters. The MR data file is called 2dseq and located in the subfolder "pdata/nr" (number of reconstruction) of the data directory.

Motion Artifacts Check
Abdominal imaging comes with complex bulk physiological motion due to the interplay of respiratory and bowel movement. If respiratory triggering of the data acquisition was used (recommended) there should only be minor motion artifacts. Check images for artifacts from motion (nontriggered acquisition) or possible residual motion (triggered acquisition). In the image viewer or Fiji go to the central slice and scroll rapidly from the first to the last echo image.

Check for motion between the echoes. If motion is noticeable
then information close to structure boundaries may not be reliable and image registration should be used for motion correction (registration is a separate topic that is not explained in this chapter).

Susceptibility Artifacts Check on T 2 *-Weighted Images
Kidney regions adjacent to bowels or in close proximity to skin/ fat/muscle boundaries or air cavities are particularly challenging to T 2 * imaging and prone to loss of anatomical integrity due to geometric distortions and signal loss created by susceptibility artifacts induced by the air-filled bowels, cavities and tissue interfaces surrounding the kidneys.
1. Scroll to the images acquired at the later echoes, where susceptibility artifacts are more severe.
2. Check areas affected by susceptibility artifacts, that is, areas in the kidney with unusual and severe hypointensities (typically in the form of spherical dark "shadows"; Fig. 4).
3. Make a record of areas affected by such artifacts (notes, screenshots). They must be excluded from the analysis or at least be considered during interpretation of the result.

Signal-to-Noise Ratio Check
Before the analysis it is highly recommended to define an SNR acceptance threshold (see Note 3).
1. Draw a region-of-interest (ROI) over the inner medulla in the first echo image (in Fiji use the "Freehand Selection" tool) (see Note 4). 4. Measure the standard deviation of the signal within this background ROI.
5. Divide this value by 0.655 to obtain the noise standard deviation. This is necessary because the background noise in MR magnitude images differs from noise in signal regions (see Note 6 if you use a multichannel RF array).
6. Calculate the SNR by dividing the mean signal of the inner medulla ROI by the noise standard deviation.
7. If the SNR is significantly below your predefined acceptance threshold, this indicates that there might be a technical problem or that the animal is not positioned correctly relative to the receiver RF coils. For more information on adequate SNR thresholds refer to Note 7.

Methods
Tools for calculating T 2 ( * ) -maps are provided by most MRI vendors, but, as mentioned in the introduction, we recommend creating your own analysis tool to have the maximum level of control over the processing steps and freedom to adapt them to your specific needs. Hence, this section will describe both how to use vendor's analysis tools and how to develop your own analysis software.
In both cases, mapping of T 2 and T 2 *, it may be required to exclude some data from the analysis, such as data points with insufficient SNR in the tail of the exponential decay. In the case of T 2 -mapping from data acquired with a MSME method, it is also recommended to discard the first echo (see Note 8).

Data Exclusion and Model Fitting (Siemens Syngo)
When using the Siemens scanner console this is currently not possible, but one can manually select the echo images that should be included when computing the fit. Using this feature, one can exclude entire images at high TEs. We recommend using a cut-off for echo images with low SNR. Please exclude all echo images with SNR < 2.6 in the ROI of interest (see Note 7). When performing T 2 -mapping, the first echo should also be excluded from the data analysis (see Note 8).
1. Load image series in the Viewing Task Card.
2. Using the ROI tool to draw an ROI on the image with the largest TE over the region of interest, and register the mean signal.
3. On the same image, use the ROI tool draw an ROI on the background, and register the standard deviation of the signal. Divide this value by 0.655 to obtain the noise standard deviation (see Note 6 if you use a multichannel RF array).
4. Calculate the SNR by dividing the mean signal ROI by the noise standard deviation. If the SNR is below 2.6 this image should be excluded on step "8" (see Note 7).

Select a 3 Â 3 or 4 Â 4 layout in the
View tab of the Control Area (right panel) in order to see all acquired echo images. 6. Select all echo images for T 2 *-mapping or all echo images except the first for T 2 -mapping (Fig. 5).
7. In the Control Area choose Eval > T 2 > OK. Two new images will appear in the viewer, the T 2 ( * ) -map and the S 0 -map.
8. Exclude all echo images that have a SNR < 2.6 (see Note 7).   8. In the measurements panel, select the ROI icon and draw the desired region of interest. An annotation with the area, mean value and standard deviation of the ROI will appear.
3. By default, T 2 /T 2 * reconstruction includes the first echo and performs a MLE fitting of the signal curve. Research customers can modify these parameters under the Scan Control Parameters panel select Reconstruction and choose the option Include first echo for (SE) T2,R2 calc (yes/no), and reconstruction method (MLE [5,6] or Least squares (RLSQ)) for T 2 ( * ) , R 2 ( * ) . If available we recommend to use Include first echo ¼ no and the default MLE calculation method (see Note 9)-Warning: changing these options affects every scan until reset of the system (Fig. 8). After the acquisition, drag and drop the images in the viewer.

Data Preparation and Model Fitting (Custom Program)
The following custom program for T 2 ( * ) relaxometry contains code examples following MATLAB syntax. These were tested to be compatible with Octave and should be easily reproducible in Python or an equivalent SDE.

Data Import and Exclusion of First TE (Only for T 2 )
1. Import your scan with the images obtained for different TEs.
Bruker data: in MATLAB use the import function provided by Bruker's pvtools-to obtain Bruker's pvtools send an e-mail to Bruker mri-software-support@bruker.com.
DICOM data: in MATLAB use the dicomread functionhttps://mathworks.com/help/images/ref/dicomread.html; in Python use dcmread from the pydicom package-https:// pydicom.github.io/). As input you need to provide a string containing the full path of the DICOM file.

Rician Noise Bias Correction
MR signal intensities are overestimated in low SNR magnitude images. This bias is due to the Rician distribution of noisy MR magnitude data and can be corrected in post-processing. We can understand the processing step as multiplication with a signal-level dependent correction factor, which approaches 1 as higher SNRs are reached. To incorporate Rician noise bias correction into our custom program, we first need to determine the noise level and then utilize a MATLAB/Octave tool provided for this book under https://github.com/LudgerS/MRInoiseBiasCorrection: 1. Add the downloaded noise correction tool to the search path.

Compute unbiased data by executing the following:
imgData_corrected¼correctNoiseBias(imgData_forAnalysis, sigma, 1) If you work with data from a multichannel RF array, change the last argument of the function to the number of receive elements.

Fitting Model
The most common method is to iteratively fit the model of the T 2 ( * ) decay to the signal intensity (SI) data of each pixel using the following equation: where S 0 is a scaling factor that includes many parameters such as the proton density together with the signal gain of the system.

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 S 0 use the first and largest SI (startVal_S0 ¼ SI_vector (1)).
For estimating a starting value for T 2 ( * ) , one can exploit the fact that at TE ¼ T 2 ( * ) the SI has fallen to 37% of its value at TE ¼ 0. Estimate the starting value by determining the TE whose SI is closest to 37% of the largest SI value: 4. Normalize the SI vector to a maximum value of 1 (divide by the largest SI value).
6. Get the absolute values of the SI vector (resulting in only positive values).

Find the index of the smallest value in the vector (in MATLAB:
[min_Value,min_Index] ¼ min(SI_norm_abs_37pct)).
8. Get the respective echo from your echoes using the index from the last step (startVal_T2¼TE_forFitting(min_Index)).

Fitting Algorithm
Least squares algorithms, such as the Levenberg-Marquardt [7] and Trust-region methods [8], are the most commonly used curve fitting algorithms for T 2 ( * ) -mapping. They work by minimizing a cost function, which describes the deviation of the fitted curve from the corresponding 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 an erroneous 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 2 ( * )mapping, where such limits can easily be defined, the Trust-region method is well suited (see Note 10) 1. Create for loops for pixel coordinates x and y to step through all pixels of the image.

Visual Display
1. Display the parameter map, which is a matrix with floating point numbers, as an image (in MATLAB: imagesc(mapT2);).
2. Remove axis labels and ensure that the axes are scaled such that the aspect ratio of the image is identical to that of the acquired FOV (For the case of square pixels in MATLAB: axis off; axis equal;).

Quantification
Quantitative values can be obtained from the parameter maps using manually drawn ROI for the different morphological regions (cortex, outer medulla, and inner medulla). Because manual ROI drawing can introduce unwanted additional variability or bias, we recommend the use of semiautomated methods, such as the concentric objects technique [9] or the morphology-based ROI placement technique [10]. For further details and step-by-step protocols for these two techniques please refer to the chapter by Riazy L et al. A simple approach to validate the data analysis procedure is to generate synthetic image data and then evaluate how close the results produced by the analysis are to the known true value. A series of artificial images could be created (using an in-house software program) that mimic a multiecho MR experiment. Each image could represent a chosen TE for which the signal intensity in a virtual phantom (could simply be a circle in the center of the image) has been calculated using the known exponential equation that describes the T 2 ( * ) relaxation. Alternatively, more sophisticated simulation approaches as previously described in [3] could be used, accounting for the timings and characteristics of the used RF pulses so as to evaluate the effect of using the simpler monoexponential decay model. Having generated the images, Rician noise is added. This synthetic data could then be analyzed, evaluating both accuracy (how close the resulting T 2 ( * ) is to the true T 2 ( * ) ) and precision (evaluating dispersion over multiple instances with the same level of SNR). The obtained values can be compared with those reported for healthy animals in the literature. Table 1 provides ranges of T 2 * and T 2 for different magnetic field strengths based on current literature.

Notes
1. With increasing echo number additional diffusion weighting is added due to the repeated magnetic field gradients used for spatial encoding (NB: this effect is much more pronounced in small animal systems, which use much stronger gradients than clinical systems). A greater bias is introduced by the superposition of stimulated echoes in multi-spin-echo imaging, which leads to significant overestimation of T 2 . It also demands the exclusion of the first echo (pure spin echo) from the analysis, because its magnitude is usually much smaller than that of the following echoes, which are combinations of spin-echoes and stimulated echoes. Although these biases depend on the acquisition parameters, they are fixed and reproducible for each protocol and hence acceptable in studies focusing on relative differences/changes, where precision is far more important than accuracy. If needed, more accurate T 2 values can be obtained from multiecho data using sophisticated postprocessing, such as dictionary-matching methodologies [3].
2. Make sure your data is ordered from the lowest to the highest TE as in Figs. 2 and 3.
3. When establishing the MR technique, define an SNR acceptance threshold. The aim is to have at least three (preferably five or more) number of echoes with an SNR > 5. This threshold will depend on the expected T 2 ( * ) values, which in turn depend on parameters like the magnetic field strength, shim quality, and tissue properties (pathology). Example: for a rat imaged at 9.4 T using a 4-element rat heart array receiver surface RF coil together with a volume resonator for excitation in combination with interventions leading to strong hypoxia an SNR > 60 was needed for the image with the lowest TE. Always draw the ROI over the same region of the kidney, that is, cortex, outer medulla, or inner medulla (we use the inner medulla, which is the brightest region of the kidney). If reaching a sufficiently high SNR is a problem, consider creating ROIs and performing the T 2 ( * ) -fitting on the average values of these ROIs in order to boost the SNR (SNR scales with the square root of the number of pixels). Another option is to modify the acquisition protocol, reducing the attained spatial resolution to gain SNR. 4. As a good practice you should always save your ROIs for later reference. After creating an ROI click on Analyze > Tools > ROI Manager > Add. Then rename your ROI to a meaningful name (using the button "Rename on the ROI Manager).

5.
When you create an ROI over a region with no signal make sure no artifacts are present. Click on Image > Adjust > Brightness and Contrast and drag the scroll bar "Minimum" to the left.
6. Sum-of-squares reconstruction of multichannel data changes the distribution of noisy MR magnitude data. Substitute the factor 0.655 by 0.682 or 0.695 for two channel and four channel data, respectively [13].
7. The Rician distribution of MR magnitude images introduces a signal level-dependent positive bias. The relative error introduced by this bias increases with decreasing SNR. For single receive element coils, an error of 10% is reached at SNR ¼ 2.6. For two-channel or four-channel coils, the corresponding threshold is at SNR ¼ 4.2 or 6.4, respectively, assuming sumof-squares reconstruction. Under pathophysiological conditions T 2 ( * ) can become significantly shorter and images acquired at high TEs may have a rather low SNR or may even reach noise level. This could result in a poor model curve fit and overestimate the calculated T 2 ( * ) . Therefore, data points at high TEs below a predefined SNR threshold should be excluded from the data analysis. Data exclusion can be done either during the acquisition (by not acquiring some echoes) or performed prior to data analysis. The latter is advantageous because the SNR cutoff to exclude data from the fitting can also be adjusted for each pixel individually. 8. When a multiecho sequence is used for T 2 -mapping, the first echo must be excluded from the data analysis. The reason is that the first echo is a pure spin echo, unlike all following echoes, which are superpositions of spin-echoes and stimulated echoes. As a result, the first echo has significantly smaller amplitude than the second echo and fitting the model curve to all data points would produce a poor curve fit and an inaccurate T 2 value [10].
9. MLE option refers to "Maximum-likelihood estimator" that takes into account the Rician noise. The noise level is measured as part of the scan. This approach should be robust against noise floor estimation bias and therefore there is no need to exclude later echoes from the analysis [5,6]. 11. Pseudocolor representations can be extremely useful for analyzing T 2 ( * ) -maps, since they generally enhance 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 emphasised 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 color-scale. It is recommended to always use the same color map and scale to improve comparability. PARENCHIMA (renalmri.org) is a community-driven Action in the COST program of the European Union, which unites more than 200 experts in renal MRI from 30 countries with the aim to improve the reproducibility and standardization of renal MRI biomarkers.