A General Linear Relaxometry Model of R1 Using Imaging Data

Purpose The longitudinal relaxation rate (R1) measured in vivo depends on the local microstructural properties of the tissue, such as macromolecular, iron, and water content. Here, we use whole brain multiparametric in vivo data and a general linear relaxometry model to describe the dependence of R1 on these components. We explore a) the validity of having a single fixed set of model coefficients for the whole brain and b) the stability of the model coefficients in a large cohort. Methods Maps of magnetization transfer (MT) and effective transverse relaxation rate (R2*) were used as surrogates for macromolecular and iron content, respectively. Spatial variations in these parameters reflected variations in underlying tissue microstructure. A linear model was applied to the whole brain, including gray/white matter and deep brain structures, to determine the global model coefficients. Synthetic R1 values were then calculated using these coefficients and compared with the measured R1 maps. Results The model's validity was demonstrated by correspondence between the synthetic and measured R1 values and by high stability of the model coefficients across a large cohort. Conclusion A single set of global coefficients can be used to relate R1, MT, and R2* across the whole brain. Our population study demonstrates the robustness and stability of the model. Magn Reson Med, 2014. © 2014 The Authors. Magnetic Resonance in Medicine published by Wiley Periodicals, Inc. Magn Reson Med 73:1309–1314, 2015. © 2014 Wiley Periodicals, Inc.


INTRODUCTION
At a given field strength, the rate of longitudinal recovery (R 1 ¼ 1/T 1 ) is determined by the microstructural tissue environment (eg, the local mobility of water molecules), the macromolecular content, and the local concentration of paramagnetic ions such as iron or gadolinium-based contrast agents. The linear dependence of R 1 on these tissue properties has been modeled and verified with independent laboratory measures from excised tissue (1)(2)(3)(4)(5)(6)(7) or by reference to literature values in the case of iron content (8).
Rapid cross-relaxation with macromolecules facilitates energy exchange between excited protons and their local environment, increasing the observed R 1 recovery rate. When macromolecular protons are selectively saturated using off-resonance radiofrequency (RF) irradiation, the MR water signal is attenuated by magnetization transfer (MT) (9). In MRI of humans, MT contrast is invoked by application of off-resonant RF pulses prior to excitation. In voxels with a higher macromolecular content, the mobile water will experience a greater percentage loss of signal (MT saturation) (10) as a consequence of a given prepulse and the dynamic of the MT (11). Consequently, measures of MT provide information about the macromolecular content of the microstructural environment and provide us with surrogate markers or model estimates for the bound water fraction (12). In a gradient echo acquisition, the observed MR signal created after RF excitation decays over time by the effective transverse relaxation rate (R 2 * ¼ 1/T 2 *). The presence of paramagnetic metals, including iron, leads to local distortion of the static B 0 field causing a more rapid decay (ie, higher R 2 *) (13).
In this study, we used a quantitative multiparameter mapping (MPM) protocol to acquire whole brain quantitative maps of the longitudinal relaxation rate (R 1 ), effective transverse relaxation rate (R 2 *), and percent saturation due to MT. We used a linear model of R 1 that is motivated by the linear dependence of the measured R 1 on the relaxation rates of the individual spin pools under the conditions of fast exchange between pools. The measured MT saturation and R 2 * estimates served as surrogate concentrations for macromolecules and iron, respectively, to calculate the model coefficients. The validity of estimating a single set of model coefficients for the whole brain-including gray matter, white matter, and deep gray matter structures-was assessed by calculating the Pearson coefficient between the measured R 1 values and synthetic R 1 values generated using the model coefficients. A set of model coefficients was calculated individually for each of 138 volunteers to examine their stability across a population.

Data Acquisition
A whole brain quantitative MPM protocol (14) was used to acquire 1-mm isotropic data from 138 healthy volunteers (male, n ¼ 49; female, n ¼ 89; age range, 19-75 y [mean age, 46.6 y; standard deviation, 21 y]) on two 3T whole body MR systems (Magnetom TIM Trio, Siemens Healthcare, Erlangen, Germany, 69 volunteers per scanner) each equipped with a standard 32-channel head coil for receive and RF body coil for transmission. The protocol consisted of PD-, T 1 -, and MT-weighted multiecho FLASH acquisitions and additional B 1 field calibrations as described by Weiskopf et al. (14). The total scanning time of the MPM protocol was $25 min. Informed written consent was obtained prior to scanning with approval from the local ethics committee.
Quantitative maps were derived from the MPM protocol using bespoke MATLAB tools (Mathworks, Natick, Massachusetts, USA). In brief, regression of the log signal from the eight PD-weighted echoes was used to calculate a map of R 2 *. The first six echoes of each of the three acquired weightings were then averaged to increase the signal-tonoise ratio (15). The resulting PD-weighted, T1-weighted, and MT-weighted volumes were used to calculate maps of MT and R 1 as described previously (10,14,16). To maximize the accuracy of the R 1 maps, inhomogeneity in the flip angle was corrected by mapping the B 1 þ transmit field according to the procedure detailed in the study by Lutti et al. (17) and the intrinsically imperfect spoiling characteristics were corrected using the approach described by Preibisch and Deichmann (18). Use of the measured R 1 maps in in vivo histological studies illustrated the high level of accuracy of the technique (19)(20)(21).
The semiquantitative MT map depicts the specific percentage loss of magnetization caused by a Gaussian RF pulse (4 ms duration, 220 nominal flip angle) applied 2 kHz off-resonance prior to nonselective excitation. This differs from the commonly used MT ratio (MTR, percentage reduction in steady state signal) by explicitly accounting for spatially varying T 1 relaxation times and flip angles (10) and results in higher contrast in the brain than the MT ratio (22). Additional minor corrections for flip angle inhomogeneity in the MT maps were applied as described by Weiskopf et al. (14).

Linear Relaxometry Model
Under the conditions of fast exchange, the crossrelaxation time between different water components is assumed to be much shorter than the MR relaxation times. In this case, the measured longitudinal relaxation rate is a weighted sum of the relaxation rates of the various contributory components (7): where f i is the fraction of spins in pool i with relaxation rate R 1i . In the absence of any exogenous contrast agents, the measured R 1 will predominantly depend on the fraction of free water spins, as well as the fraction of bound water spins at macromolecular sites and, on a smaller contribution, from iron sites (8,23). In this case, equation 1 becomes: where R 1f is the relaxation rate of free water; f M is the fraction of spins bound to macromolecules; r 1M is the relaxivity at macromolecular sites (ie, R 1M À R 1f ) where R 1M is the relaxation rate at macromolecular sites, f FE is the fraction of spins at iron sites, and r 1FE is the relaxivity at iron sites; the index j sums over all potential unspecified contributions that remain.

Estimating Model Parameters
A model of R 1 based purely on imaging data can be constructed from Eq. 2 by replacing the known contributors to R 1 with surrogate markers. Using the MT and R 2 * maps from the MPM protocol as surrogate markers for the macromolecular and iron concentrations, respectively, the model can be expressed as a function of spatial position, r, as: Here the set of b parameters are global constants and e(r) is the spatially specific residual encompassing the unspecified contributions to R 1 and noise. For a given subject, a model matrix, M, is constructed with three col-umns: unity, MT, and R 2 * (Fig. 1). In this case, the b parameters are the model coefficients and least squares solution to the matrix equation R 1 ¼ Mb, and the residuals are the difference between the measured and synthetic R 1 values (ie, e ¼ R 1 À Mb).
The parameters of the linear relaxometry model were estimated for each volunteer separately by fitting the relationship across all brain voxels (excluding cerebrospinal fluid, which does not contain appreciable amounts of macromolecules or iron and is not wellcharacterized by the particular acquisition protocol used in this study). Brain voxels were identified by automated segmentation using SPM8 (24). All voxels with a gray or white matter probability >50% and a cerebrospinal fluid probability <50% were pooled to estimate the global, tissue-independent b parameters as well as a spatial map of the residual error on a subject-by-subject basis. The accuracy of the model was assessed by calculating a synthetic R 1 map from the b parameters (ie, Mb). The Pearson coefficient between the synthetic and measured R 1 values was calculated for each subject.

RESULTS
A single subject example data set is shown in Fig. 2 along with the model estimated R 1 and the residuals (ie, the difference between the estimated and measured R 1 values). The linear model fitted well with a mean Pearson coefficient across the entire cohort of 0.93 6 0.03. The resulting b parameters are summarized for the group in Table 1. The coefficient of variation across volunteers was <5.5% for both b 0 and b 1 but rises to 37.2% for b 2 .
The residuals were generally close to or within the acquisition noise level (14). Within the whole brain, they are centered at À0.13 6 0.34% across the cohort. Examining the residuals by tissue class revealed a bias of 1.62 6 0.38% toward lower model R 1 values in white matter and a bias of 0.97 6 0.52% toward higher model R 1 values in gray matter. The residual maps show that the greatest deviation of the model to the measured R 1 occurred within the basal ganglia, where it was underestimated by the model to varying degrees but in the order of 10%.
The b 0 parameter corresponded directly to the relaxation rate of free water. The estimated longitudinal relaxation time (T 1 ) of free water in brain tissue was 3.736 6 0.198s.

DISCUSSION
A single set of b parameters sufficed to model R 1 across the whole brain, incorporating voxels from both gray and white matter. The model parameters describing the influence of free water and MT on R 1 were highly stable across this large cohort, suggesting that underlying microstructural differences due to macromolecular content are well captured by the MT map across both gray and white matter. This may be due to the fact that the transfer times are similar in gray and white matter (25). In agreement with the well-established two-site-fastexchange model (2,7,23,26), the contribution from macromolecular components has been shown to be the dominant factor within the model. In addition, it has been shown that when gray and white matter are considered concurrently water content is a better predictor of R 1 variation than is iron content (8). The b 0 parameter appears to give a reliable estimate of the longitudinal relaxation rate of free water (R 1f ) that is in keeping with values reported in the literature [see Table 1 in Rooney et al. (23)]. The remaining components of the model contain inseparable components: the fraction of the spins at the site and the relaxivity of the site. For example, the b 1 MT term of the model corresponds to the product of the macromolecular bound water fraction and the relaxivity at these macromolecular sites. Either of these components may change depending on the local microstructure. To separate the contributions and validate this model further, one or the other factor must be known. The bound water fraction can be measured histologically using calorimetric approaches, as has been done for white matter (3). Using the histological estimate of water fraction and combining it with our results yields a relaxivity estimate of 3.687 6 0.198s À1 for WM, which is in good agreement with the literature (23).
Post mortem validation has shown high correspondence between R 2 * and iron content (13) and between MTR and myelin content (12). The fact that the b parameters largely showed good stability across this large cohort (age range, 19-75 y) validates our use of MT and R 2 * as surrogate markers for macromolecular and iron content. When the MTR measure was used in the model, the model coefficients were less stable and the model residuals were increased with anatomical and bias field structure present (Fig. 3). The stability of the model coefficients also suggests that they could be used to generate an estimated R 1 map directly from MT and R 2 * maps. Applying the model in such a way would facilitate a reduction in the total scan time for the quantitative MPM approach. It would also allow a synthetic quantitative map to still be calculated from suboptimal data (eg, if motion occurred during part of the protocol). The volunteers included in this analysis were all healthy volunteers with no evidence of cognitive impairment. Because we did not study patients, we cannot conclude whether a different set of model parameters may be required in pathological conditions. If they were, then the resulting parameters may become a diagnostic measure.
The model parameter describing the R 1 dependence on R 2 * (ie, b 2 ) had a markedly higher coefficient of variation across volunteers than the free water or macromolecular terms. For the majority of the brain, the model performs well without the inclusion of the iron term, and the Pearson coefficient remains high (mean, 0.92). However, inclusion of the iron term via the R 2 * measure reduces the residuals for high iron structures, such as the basal ganglia and dentate nucleus (Fig. 4). This is to be expected, since it is in these structures that iron will contribute most significantly to R 1 . Nonetheless, the residuals remain highest in these structures, suggesting that the iron effects are not described fully by the model. Additional or higher-order terms (eg, accounting for the different chemical forms in which iron is present in vivo (27)), may be required to more fully model the true R 1 variation in these structures. Additionally, the comparatively smaller contribution from iron sites may be poorly estimated because of the higher noise level in the R 2 * map estimated by our MPM protocol which uses a maximal echo time of 19.70 ms and because the multiecho three-dimensional FLASH volumes used to calculate the MT and R 1 maps are averaged resulting in an effective echo time of 8.45 ms and therefore residual R 2 * weighting. It has also been suggested that the measured R 2 * will depend on the orientation of white matter fibers with respect to the main field (28)(29)(30). These effects are currently not accounted for, though it may be possible to extend the model to include them.
The residuals of the model fit can have a variety of sources, typically falling into two categories: unknown or unspecified model terms, as discussed above, and noise. In some cases, the residuals are high because they capture a spatially coherent form of noise (ie, artifact) present in one or more of the maps used in the model. This is most typically caused by subject motion during the acquisition of one of the constituent three-dimensional FLASH volumes of the MPM protocol. Because the volumes combine in different ways to calculate each of the quantitative maps such motion can lead to inconsistencies across the model matrix components, which may in turn be captured by the residuals. This feature raises the interesting possibility of using the general linear relaxometry model for artifact characterization and correction. Many opportunities arise from the application of this model such as further understanding R 1 dependencies, potential disease biomarkers, and artifact correction.

CONCLUSIONS
We have shown that a general linear model of longitudinal relaxation can be applied voxel-wise across the whole brain by using MT and effective transverse relaxation maps as surrogate concentrations for macromolecules and iron, respectively. This model fits well and provides a single set of model parameters per individual that is remarkably stable across the cohort.