Repeatability of Quantitative Sodium Magnetic Resonance Imaging for Estimating Pseudo-Intracellular Sodium Concentration and Pseudo-Extracellular Volume Fraction in Brain at 3 T

The purpose of this study is to assess the repeatability of the quantification of pseudo-intracellular sodium concentration (C1) and pseudo-extracellular volume fraction (α) estimated in brain in vivo using sodium magnetic resonance (MRI) at 3 T. Eleven healthy subjects were scanned twice, with two sodium MRI acquisitions (with and without fluid suppression by inversion recovery), and two double inversion recovery (DIR) proton MRI. DIR MRIs were used to create masks of gray and white matter (GM, WM), that were subsequently applied to the C1 and α maps calculated from sodium MRI and a tissue three-compartment model, in order to measure the distributions of these two parameters in GM, WM or full brain (GM+WM) separately. The mean, median, mode, standard deviation (std), skewness and kurtosis of the C1 and α distributions in whole GM, WM and full brain were calculated for each subject, averaged over all data, and used as parameters for the repeatability assessment. The coefficient of variation (CV) was calculated as a measure of reliability for the detection of intra-subject changes in C1 and αfor each parameter, while intraclass correlation (ICC) was used as a measure of repeatability. It was found that the CV of most of the parameters was around 10–20% (except for C1 kurtosis which is about 40%) for C1 and α measurements, and that ICC was moderate to very good (0.4 to 0.9) for C1 parameters and for some of the α parameters (mainly skewness and kurtosis). In conclusion, the proposed method could allow to reliably detect changes of 50% and above of the different measurement parameters of C1 and αin neuropathologies (multiple sclerosis, tumor, stroke, Alzheimer’s disease) compared to healthy subjects, and that skewness and kurtosis of the distributions of C1 and αseem to be the more sensitive parameters to these changes.


Introduction
This study was approved by the institutional review board (IRB) of New York University Langone Medical Center and all volunteers signed informed consent form (#11783, Technical Development of up to 7T Magnetic Field Level) prior to the scans.

MRI hardware
All scans were performed at 3 T on a Tim Trio system (Siemens, Erlangen, Germany) using a dual-tuned 1 H/ 23 Na birdcage radiofrequency (RF) coil tuned at 128/33 MHz (Stark Contrast, Erlangen, Germany).

Proton MRI acquisition
Two double inversion recovery (DIR) MRI acquisitions were performed. The first DIR image was acquired in order to suppress both CSF and WM using the DIR Turbo Spin Echo SPACE sequence [24,25] with the following parameters: TR = 7500 ms, TE = 300 ms, field-of-view (FOV) = 220×320×320 mm 3 isocenter, resolution = 2.5 mm isotropic, inversion times TI 1 = 2650 ms and TI 2 = 550 ms, time of acquisition (TA) = 4:00 min. The second DIR image was acquired in order to suppress both CSF and GM with the same parameters as the first DIR except TI 1 = 2800 ms and TI 2 = 800 ms.
All sodium images were reconstructed offline in Matlab (MathWorks, Natick, MA, USA) with standard 3D regridding [28] and density compensation [29] with a nominal resolution of 2.5 mm isotropic (128×128×128 voxels), matching the nominal resolution of the DIR proton images.

Data processing summary
The data processing (in Matlab) for calculating the final pseudo-intracellular sodium concentration (C 1 ) and pseudo-extracellular volume fraction (α) maps is described in details in Ref. [23]. In the present article, we will just summarize this data processing in 4 steps: 1. Calibration phantoms and linear regression: The signal from five calibration phantoms (Agar gel 3% with 10, 30, 50, 70 and 100 mM NaCl, 17 mm diameter and 100 mm length) placed within the FOV on the side of the head was measured and averaged over 4 consecutive slices (10 voxels/phantom/slice). Their relaxation times were also measured as T1 = 38 ms and T2 Ã = 7 ms at 3 T. The loss of signal of the sodium phantoms due to relaxation during RF pulses and delays was estimated by full density operator simulation of spin 3/2 dynamics [30][31][32] during the RF pulse sequence and was used to correct the phantom signals, by factors the 1.10 (sequence 1) and 1.60 (sequence 2). Linear regression of the corrected signals of the phantoms versus their sodium concentrations was then performed (Fig. 1A). Linear regression was considered as valid only when the coefficients of determination R 2 ⩾ 0.99 and adjusted R 2 adj ⩾0:98 [33], in order to improve the robustness of the method against noise and signal variations in the phantoms.
2. Intermediate sodium maps: The apparent total sodium concentration (aTSC) and apparent pseudo-intracellular sodium concentration (aISC) maps were calculated from sequences 1 and 2, respectively, using the coefficients of the linear regression of phantom signals (Fig. 1A). Using average sodium relaxation times in brain from the literature [2,14,20,21,27] (T1*35 ms, and T2 short *5 ms, T2 long *25 ms, in parenchyma), correction factors of the sodium maps (0.85 for aTSC and 0.5 for aISC) were calculated using full density operator simulation of the sodium spin dynamics [30][31][32] during the RF pulse sequences 1 and 2 respectively.
3. GM and WM masks: 3D GM and WM and full brain (WM+GM) masks were calculated from the 1 H DIR acquisitions using SPM8 [34] in Matlab. From the resulting GM (or WM) SPM probability map, only pixels with a probability ⩾ 0.75 of being in the GM (or WM) were kept in the GM (or WM) mask. The aTSC and aISC maps were then multiplied by the GM, WM and full brain masks (Fig. 1A). These masked aTSC and aISC maps were therefore used for the quantification of C 1 and α in order to measure the distributions of these 2 parameters separately in WM, GM and full brain.
4. C 1 and α maps: Pseudo-intracellular sodium concentration (C 1 ) and pseudo-extracellular volume fraction (α) quantification was based on a simple three-compartment model (Fig. 1B). In this model, the extracellular compartment (including interstitial volume, CSF, plasma and blood) has a constant average sodium concentration C 2 *140 mM [2,18,20,21,35]. We also considered that the water (fluid) volume fraction is constant and take averages values w WM = 0.7, w GM = 0.85 and w brain = 0.775 (mean value from WM and GM) [36][37][38]. We also assumed here that all extracellular sodium signals are completely suppressed by inversion recovery in sequence 2, or within noise level of the image (due to imperfection of inversion pulse and possible variations of T1 from CSF and other extracellular fluids). The value of each voxel of the aTSC map is by definition equal to the total sodium concentration within each voxel: aTSC = (C 1 ×V 1 +C 2 ×V 2 )/V t (with V t = total volume of the voxel). The value of each voxel of the aISC map is by definition equal to the intracellular sodium concentration only: aISC=(C 1 ×V 1 )/V t . From these assumptions and equations, we can calculate the unknown parameters C 1 and α of interest, using the relationships given in Fig. 1B and with w taking the values w WM , w GM and w brain depending on the masked aTSC and aISC maps used: α = (aTSC−aISC)/C 2 and C 1 = (C 2 ×aISC)/(w×C 2 −aTSC+aISC). This calculation is performed for each voxel. All voxels are then recombined in maps of C 1 and α in WM, GM and full brain, as shown in Fig. 1C. The distribution of values of C 1 and α over whole WM, GM or full brain can also be plotted as histograms (Fig. 1D), and the properties of each distribution can be characterized for each volunteer in terms of histogram summary statistics: mean, median, mode, standard deviation (std), skewness and kurtosis. The reliability of the method will be assessed in terms of repeatability of these six summary statistics.
Note that no RF (B 1 ) map was acquired in this preliminary study, as it was not deemed necessary from preliminary data acquired while testing this technique. We used a birdcage coil with homogeneous transmit and receive B 1 fields, no improvement in the quality of the images, nor in the sodium data quantification was observed when B 1 correction was applied. Moreover, B 1 correction for the fluid suppressed images (from inversion recovery) is not trivial, as the effect of B 1 inhomogeneities on the images depends not only on the pulses (inversion and excitation) but also on the T1 of the tissues and the relaxation of their magnetization during the inversion time, which varies in different compartments of the brain.

Lesion simulations
Simulations of 'fluid' (cystic-like) and 'solid' (tumor-like) lesions, either compact (10×10×10 voxels) or randomly distributed over the whole brain (1000 voxels over 86571, which correspond to about 1.15% of all voxels in brain) were also performed and are presented in Supporting Information (S1 File). These simulations are similar to the ones presented in Ref. [23]  . Three-compartment model: compartment 1 = intracellular space of volume V 1 and sodium concentration C 1 ; compartment 2 = extracellular space of volume V 2 , sodium concentration C 2 and extracellular volume fraction α; compartment 3 = solid compartment of volume V s and no sodium. Total volume is V t = V 1 + V 2 +V s . Water volume fraction is denoted w. C. C 1 and α maps in full brain, GM, and WM calculated from aTSC, aISC and the 3-compartment model. D. Histograms of C 1 and α values over the whole full brain, whole GM and whole WM data. GM = gray matter, WM = white matter, GM+WM = full brain. Reproduced from Figs. 1-4 from Ref. [23].
Supplementary Material. These simulated lesions include voxels with partial volume effect from the CSF. Random 1-voxel lesions are distributed over the whole brain (including voxels with partial volume effect from CSF). Details are given in the captions of the figures A to D in S1 File.

Statistical analysis
Restricted maximum likelihood estimation of the variance components in a random effects model was used to estimate the intra-subject variance (i.e., the variance between results from replicate scans of the same subject) and inter-subject variance (the variance between results from different subjects) of each measure (i.e., mean, median, mode, std, skewness and kurtosis) within each tissue. The estimated variance components were used to compute the coefficient of variance (CV, in %) as the square root of the intra-subject variance expressed as a percentage of the mean, and the intra-class correlation (ICC) as the inter-subject variance divided by the sum of the intra-and inter-subject variances. By expressing the intra-subject variance relative to the overall mean of a given measure, the CV is considered an indicator of the utility of a measure for detecting within-subject changes over time; temporal changes that are not substantively greater in magnitude than the CV would be difficult to distinguish from changes attributable to random noise. It is noted that the CV is not appropriate for measures (e.g., skewness) that can be both positive and negative since the mean of such a measure does not represent the magnitude of a typical observed value of the measure. The ICC is a measure of the repeatability (or test-retest reliability) of the method [39]. It ranges from 0 to 1, with values near 1 implying that the ratio of the intra-subject variance to the inter-subject variance is close to 0. As a result, the ICC should be high (close to 1) if a measure is to be potentially useful for the detection of a mean difference between subject groups. In a similar fashion as described in Ref. [40], repeatability was regarded as very good if ICC ⩾ 0.  Fig. 2 show an example C 1 and α maps and distributions from two scans of the same volunteer. Note the higher intensities, mainly in α maps, near the ventricles, due to residual partial volume effect from CSF. The C 1 and α distributions from scan and re-scan present very similar shapes between the two scans, for GM, WM and full brain (GM+WM). Fig. 3 and 4 show scatter plots of mean, median, mode, std, skewness and kurtosis of C 1 and α from the eleven volunteers scanned twice (1 st scan vs. 2 nd scan), in GM, WM and full brain. Data from the histogram statistics for all volunteers and all scans can be seen in Supporting Information (S1 Dataset). On Fig. 3, we observe that all mean, median and mode C 1 values from scan/re-scan are in the range 5-20 mM. The mean, median and mode α values are in the range 0.10-0.25, over whole GM, WM and full brain. Std C 1 values are all regrouped around 5 mM, while std α values are in the range 0.05-0.10. On Fig. 4, we observe that most of the values of skewness and kurtosis are close to the diagonal (good match between scan and re-scan measurements), with the exception of a few data points. Skewness values for C 1 are all regrouped around 0, while they are in the range 1-3 for α. Kurtosis values for C 1 are tightly regrouped in the range 3-5, while they are more spread out for α, in the range 5-15.

Results
Tables 1 and 2 present the mean all , std all , CV and ICC of the C 1 and α measurements, over all data (n = 22, from 11 volunteers scanned twice). Note that the mean and std over all data were labeled mean all and std all to avoid confusion with the mean and std of the measurements of C 1 and α for each volunteer and scan. Results can be summarized as below.
Mean all and std all C 1 -We can observe that mean C 1 over the 22 scans had a mean all ±std all of 11.4±1.5 mM in GM and slightly higher 13.6±2 mM in WM, std was 4.5±0.7 mM in both GM and WM, skewness was 0±0.3 in both GM and WM, kurtosis was 3.4±0.5 in GM and slightly higher 4.0±1.7 in WM (but with higher variability). Median and mode values were in the same range as the mean value, in GM and WM respectively. For full brain, all the measures were of the same order of magnitude than in GM and WM. respectively. For full brain, all the measures were of the same order of magnitude than in GM and WM. CV C 1 -We can observe that the CV in GM was *12% for mean, median, mode and std, and lower *6% for kurtosis. CV values were generally higher in WM: *12% for mean and median, *19% for mode, *16% for std and *40% for kurtosis. The CV values in full brain were of the same order of the values for GM and WM, except in the case of kurtosis where CV is low *6%.
α-We can observe that CV in GM was *10% for mean and median and std, *15% for mode, and much higher *27% for kurtosis. The same trend occurred in WM, with CV *12% for mean, median and std, *18% for mode and *22% for kurtosis. The CV values in full ICC C 1 -In both GM and WM, ICC was moderate *0.5 for mean, median and std, and good to very good *0.6-0.8 for mode and skewness. For kurtosis, ICC was high *0.8 in GM, but much lower *0.4 in WM, due to intra-subject variance in this tissue. For full brain, ICC was moderate *0.5 for mean, median, mode and std, but good to very good for skewness (*0.7) and kurtosis (*0.9), respectively. α-In both GM and WM, ICC was low *0.1-0.2 for mean, median and mode, good *0.7 for std, and poor to moderate for skewness and kurtosis (*0.4-0.5). In full brain however, ICC for skewness and kurtosis were both very good *0.9, while it was still good for std (*0.6), and very low for mean, median and mode (*0.1).

Discussion
Mean, median, mode and std values of C 1 were all in good agreement with values from the literature for intracellular sodium concentrations in healthy brain tissues, which are generally in the range 5-15 mM [20,41]. The mean, median and mode α values were also in good agreement with values from the literature for extracellular volume fraction in healthy brain tissues, generally in the range 0.15-0.25 [8,11]. From the mean skewness and kurtosis values, we found that the C 1 distribution over whole GM, WM or full brain was very close to the normal (Gaussian) distribution (skewness = 0, kurtosis = 3), while the α distributions in GM and WM deviate from the normal distribution. The positive skewness (range 1-2) for α indicates a higher range of values above the mean value, while higher kurtosis (> 3) indicates a more 'peaked' distribution of α, with fatter tails, compared to normal distribution. These disparities for α from the normal distribution occur mostly because of both methodological and physiological reasons. From the methodology point-of-view, due to the low resolution of the 23 Na and 1 H images, both C 1 and α maps exhibit residual partial volume effect that can influence their distributions, mainly in regions close the CSF which will therefore have high α > 0.2 (increase skewness of α). Moreover, incomplete inversion of the magnetization and noisy data can also generate voxels with both higher and lower values compared to normal distribution and therefore create fat tails. For the physiological point-of-view, the difference in kurtosis of α between GM and WM could be explained by the differences in cell packing characteristics and structure [42] between these two tissues. It is generally considered that WM has a more anisotropic structure than GM [42,43], due to the presence of axons and glial cells, associated with a wider range of extracellular volume fractions (therefore higher kurtosis that GM), and potentially more high values of α > 0.2 (higher skewness than GM). As intracellular sodium concentrations in all cells in the brain don't vary too much from the 5-15 mM range, the C 1 distributions don't differ significantly between GM and WM.
The skewness and kurtosis of both the C 1 and α distributions measured with this protocol may also be in part the results of the discrepancy between the edge definitions of the 1 H DIR images (2.5 mm resolution), the 23 Na aTSC images (5 mm resolution for the data acquisition) and the 23 Na aISC images (6.7 mm resolution). Because of all these limitations, it can be difficult with the present method to assign a relevant biological or clinical meaning to the skewness  and kurtosis of the distributions. As we are more interested in changes in these variables from healthy to non-healthy tissues, and because the exact same protocol is applied to all subjects, this lack of accurate biological interpretation should not impair the ability of the method to assess changes in C 1 and α due to pathologies. Except for C 1 kurtosis where the CV > 40%, all measures of C 1 and α can be considered good variables for detecting changes over time in individuals (intra-subject variations), provided that the changes are above 20% of their mean values. From the error propagation calculated in [23], the uncertainty on the measurements of both C 1 and α was calculated as around 40% for standard uncertainties on water fraction (w) and extracellular sodium concentration (C 2 ) approximations. Taking CV and uncertainties into account, it is therefore reasonable to estimate that in general the proposed method would be able to detect changes in C 1 and α of about 50% and above (to be on the safe side), which corresponds to mean C 1 > 20-25 mM and mean α > 0.3, both for local (regions-of-interest) or global (whole tissue) measurements. Due to all the model assumptions, the low SNR of sodium MRI, the error propagation due to uncertainties in the relaxation times, in the water fraction and in other post-processing parameters, we can estimate that changes below 50% will not be a reliable. If we consider only changes of 50% and above for mean C 1 and α, as expected in many neuropathologies from the literature  doi:10.1371/journal.pone.0118692.t002 [5][6][7][8][9][10][11], the combination of C 1 and α measurements can also provide good classification between different groups of subjects (moderate to good ICCs), mainly if measurements are performed over full brain. For example, in stroke, it is generally estimated that an apparent total tissue sodium concentration (aTSC) of 65-70 mM is a marker of irreversible tissue damage (cell death). This would correspond to an intracellular sodium concentration of about 40-50 mM (see for example Ref. [3,17,18]). We could therefore reasonably expect that 20mM < C 1 < 50 mM would correspond to loss of homeostasis that might be reversible, and C 1 > 50 mM would correspond to irreversible loss of cell viability. For cancer, tumor cells have altered physiology, such as evasion of apoptosis, limitless replicative potential, tissue invasion [44], and therefore also altered homeostasis. In that case, it can be expected that increases in C 1 and α above 50% and even more can be linked to tumor malignancy, and potentially grading and staging (See Ref. [2,[13][14][15] and references within).
As ICC for skewness and kurtosis in full brain are very good, these two variables seem to be the best features for detecting differences between subjects. From the simulations of different lesions (see S1 File), we can see that compact fluid lesions (cystic-like) can be detected on the α maps, compact solid lesions (tumor-like) can be detected on the C 1 maps, but random lesions cannot be detected on either α or C 1 maps. However, these lesions (compact and random) can be detected on the C 1 and α distributions (over whole brain), mainly through changes in skewness and kurtosis, while other measures remained almost unchanged.
As a first approximation in this study, we assume that all extracellular sodium is either completely suppressed by IR (with TI optimized to suppress fluids) or reduced within noise level (for 'bound' extracellular sodium, which can be expected to have intermediate T1 between fluid sodium T1 and bound sodium T1 from within the cells). Future studies will take this into account by using a more complex tissue model (more compartments) which will include the presence of sodium with restricted motion within the interstitial space, that is different from extracellular CSF or plasma space.
Due to time constraints for the subjects (we try to keep scans of the volunteers in about 45-60 min maximum) and for scanner availability, the sodium scans were optimized to last 11 min and 17 min by using shorter TRs (80 ms and 100 ms). Although these TRs are not optimized for full recovery of the longitudinal magnetization, we found that the parameters used were a good compromise between SNR and total acquisition time for this pilot study. Relaxation times in pathologies are unknown and might affect the quantification, but we expect these relaxation times within and outside pathological cells to be still quite similar to the relaxation times within and outside normal cells, respectively. For example, assuming an average T1 *30 ms in parenchyma [2], a variation of ±20% in T1 of the tissue will induce a variation of signal of only 3.6% in sodium images acquired without fluid suppression and TR = 80 ms, and a variation of about 5% in sodium images acquired with fluid suppression and TR = 100 ms, compared to acquisitions with TR = 150 ms (fully recovered magnetization with TR = 5×T1). These signal variations will propagate to the C 1 and α quantification and are included in the general error propagation calculated in Ref. [23] for aTSC and aISC quantification, which add up to about 40% when combined with water fraction and extracellular sodium concentrations uncertainties.
Partial volume effect is a real concern for this data and was already partially discussed in Ref. [23] (p. [4][5]. Images were acquired with different resolutions in this pilot study for practical reasons: we wanted to be able to scan people within a reasonable time without losing too much SNR, mainly for sodium with fluid suppression. We are now working on ways to improve the SNR of the images (new multichannel coil, iterative reconstruction such as compressed sensing, data denoising) in order to increase the resolution of the images, with and without fluid suppression (we expect to reach 3-4 mm isotropic for both, with about 10-12 min acquisition time for each acquisition).
No sodium signal from zero concentration in reference phantoms was measured for the linear regression, as it is hidden in the noise of the images. As sodium images have low SNR, we observed that adding measurements of [Na] = 0 mM (therefore noise) can induce large errors in the quantifications, which are very sensitive to the slopes of the linear regressions. Future increases of the SNR of the images, as already described above, will improve the robustness of the method to noise.
In conclusion, estimating the pseudo-intracellular sodium concentration (C 1 ) and pseudoextracellular volume fraction (α) in brain in vivo is feasible using a combination of DIR proton MRI and sodium MRI, with moderate to good repeatability (ICC) and CV. This preliminary method will be furthermore improved in future studies by implementation of compressed sensing reconstruction [45] and data denoising methods [46] for increasing the signal-to-noise ratio, accelerating sodium acquisitions, and/or increasing the resolution of the sodium images. Water fraction estimation with proton MRI [37,38], which was considered as fixed in the present work, will also be included in the future, as this parameter also changes with pathologies. Future studies with multichannel RF coils will also include the calculation of sodium B 1 maps and the development of B 1 inhomogeneities correction for the calculation of C 1 and α, in combination with increased SNR and resolution, to reduce CSF partial volume effect and inefficient CSF suppression effect on data quantification. Next steps will include estimating the efficiency of the method on patients for assessing different neuropathologies such as multiple sclerosis, brain tumors, traumatic brain injuries, or Alzheimer's disease, either locally on C 1 and/or α maps, or globally using the distributions of C 1 and/or α over whole GM, whole WM or whole brain.
Supporting Information S1 Dataset. Dataset of C 1 and α measurements (mean, median, mode, standard deviation, skewness, kurtosis) in GM, WM and full brain, for all scans of all volunteers. (XLSX) S1 File. Figure A. Pseudo-intracellular sodium concentration (C 1 ) maps of the brain of a healthy volunteer (1 axial slice) with artificial 'fluid' and 'solid' inclusions. Figure B. Distributions of all pseudo-intracellular sodium concentration (C 1 ) values in full brain (GM+WM, black), GM (blue), WM (red) from a volunteer, with artificial 'fluid' and 'solid' inclusions. Figure C. Pseudo-extracellular volume fraction (α) maps of the brain of a healthy volunteer (1 axial slice) with artificial 'fluid' and 'solid' inclusions.