Data-driven classification of tissue water populations by massively multidimensional diffusion-relaxation correlation MRI.

1
Introduction
Diffusion magnetic resonance imaging (dMRI) provides a non-invasive means to access microstructural information about brain tissue by measuring the displacement of water molecules. The interaction of these water molecules with the complex cellular architecture and composition gives rise to multiple water populations associated with various cellular structures, including cell somas, neurites, and diverse lipidic and chemical compositions. The MRI signal from a single voxel reflects a distinct degree of heterogeneity, resulting from a sum of the contributions from these water populations. In pathological conditions, alterations in cellular structures—both in morphology and composition—can further increase the complexity of the dMRI signal by affecting the underlying water populations.
Pierpaoli and Pierpaoli et al. (1996) initially demonstrated the separation of meaningful water populations using only diffusion MRI data through diffusion tensor imaging parameters, trace and anisotropy. The use of per-voxel single trace and anisotropy values allowed to produce the segmentation of WM, GM, and FW regions. Expanding on this basis, the diffusion MRI community has made significant efforts in addressing intravoxel heterogeneity by developing both biophysical models (Novikov et al., 2018a; Jelescu et al., 2020) and advanced dMRI sequences (Lundell et al., 2019; Aggarwal, 2020; Reymbaut, 2020; Henriques et al., 2021; Jiang et al., 2023). Biophysical modeling separates water populations by defining cellular compartments. In white matter (WM), axons are often represented as sticks, extra-axonal space as an axially symmetric tensor, and free water (FW) as a spherical tensor (Zhang et al., 2012; Novikov et al., 2018b; Coelho et al., 2022). In gray matter (GM), additional parameters account for the soma compartment (Palombo et al., 2020). Each model has specific assumptions, requiring compatible acquisition protocols for optimal application. Consequently, these models are limited to the water populations they are designed to analyze and may not be applicable to tissues outside the brain. In practice, studies focusing on WM require specific parameter estimators to capture different diffusion regimes, enhancing the fit between the models and the tissue signal. In addition to biophysical modeling, valuable information regarding water populations can also be obtained by adjusting the scanning acquisition parameters.
Recently, there has been an increased application of multidimensional acquisitions that combine b-tensor encoding diffusion and relaxation (De Almeida Martins and Topgaard, 2018; De Almeida Martins et al., 2021; Martin et al., 2021). One approach is massively multidimensional diffusion-correlation MRI (Narvaez et al., 2022), which facilitates the acquisition of per-voxel non-parametric distributions of diffusivities, anisotropies, restrictions, and chemical compositions (D(ω)–R1-R2-distributions). In this approach, the signal is obtained using free gradient waveforms by modifying the trace (b), normalized anisotropy (bΔ), orientation (Θ, Φ), centroid frequency (ωcent) with variable repetition (τR), and echo time (τE). Once the distributions are obtained, an arbitrary set of “bins” limits are selected on the isotropic diffusivity (Diso) and squared normalized anisotropy () space to separate and visualize water populations into meaningful tissue fractions (i.e., WM, and GM) and free water (FW). Similarly, earlier binning approaches included 1D R2-distributions for separating myelin water fractions (Kim et al., 2017) or 2D “spectral regions of interest” for differentiating tissue types from 2D correlation experiments (Benjamini and Basser, 2017; Slator et al., 2019; Pas et al., 2020). In non-parametric D(ω)–R1-R2-distributions, Diso and contains homologous microstructural information to mean diffusivity (MD) (Basser et al., 1994). and microscopic fractional anisotropy (μFA) (Lasič et al., 2014; Shemesh et al., 2016a; Topgaard, 2016). Even though the “bin” division of the per-voxel distributions is practical and yields adequate separation of WM-like, GM-like and FW-like water populations, the selection of the “bin” limits needs to be manually modified in order to provide meaningful tissue fractions under different conditions (i.e., in vivo, ex vivo, type of sample analyzed); and it misses the opportunity to exploit the multiple dimensions available. Moreover, the “bin” selection falls short in disentangling the complexity inherent within the voxel. In addition to well defined water population in tissue, there are potential image artifacts and not physical plausible distributions due to under-sampling of bΔ at low b-values (De Almeida Martins and Topgaard, 2018) that could be mixed with true water population using manual binning approaches. This limitation highlights the need for unsupervised methods capable of capturing biological relevant features from multidimensional diffusion-relaxation correlation experiments.
Previous studies have utilized unsupervised approaches incorporating full tensor distribution data to further characterize tissue heterogeneity in diffuse gliomas (Song et al., 2022). Full tensor distributions account for both the symmetric and asymmetric components of the tensor, in contrast to diffusion tensor distribution (Topgaard, 2019) (DTD) approaches, which consider only the axially symmetric components. Most recently, the simultaneous acquisition of diffusion data while modifying repetition and echo times has demonstrated the ability to differentiate water populations and provide greater specificity in the context of astrogliosis (Benjamini et al., 2023). Kundu et al. (2023) have introduced an unsupervised learning approach on the optimal transpose distances (LOT) from the multidimensional T1-T2, T1-MD, and T2-MD spaces to distinguish cortical layers that could not be distinguishable using only 1D relaxation or diffusion parameters.
In this study, we used Gaussian mixture model (GMM) to perform data-driven clustering of non-parametric multidimensional diffusion-relaxation distributions, extending beyond the traditional fractions of WM, GM, and FW. We explore the classification of WM, GM, and FW tissue fractions while also identifying additional detailed tissue clusters within both WM and GM. These additional clusters can represent shared region-specific water diffusion populations or possible artifacts within a voxel. The data-driven clustering was applied to a comprehensive acquisition protocol for ex vivo rat brain data, as well as shorter protocols for both in vivo rat and human brains. This work represents a significant advancement in understanding the efficacy of this approach compared to previous methodologies.

2
Materials and methods
2.1
Animals and data acquisition
2.1.1
Ex vivo healthy rat brain at 11.7 T
An adult male Sprague-Dawley rat was used for the ex vivo study (10 weeks old, 350 g, Harlan Netherlands B.V.). The rat was housed individually in a cage with climate-control room under a 12:12 h light-dark cycle with ad libitum diet. The animal procedures were approved by the Animal Ethics Committee of the Province Government of Southern Finland and carried out according to the guidelines set by the European Community Council Directives 2010/63/EEC.
First, the rat was anesthetized with 5% isoflurane in a mixture of 70% nitrogen/30 % oxygen and transcardially perfused with 0.9% NaCl for 5 min (30 ml/min) followed by 4% paraformaldehyde (PFA) in 0.1 M phosphate buffer (PB) pH 7.4 for 25 min (30 ml/min). After the perfusion, the brain was removed from the skull and postfixed in 4% PFA for 4 h. Before the ex vivo imaging, the brain was immersed in 0.1 M phosphate-buffered saline (PBS) pH 7.4 to remove the excess of PFA for 24 h. Only before the scan, the brain was divided into two hemispheres. Then, the left hemisphere was placed inside a 10-mm diameter tube filled with perfluoropolyether (Galden, Solvay, Italy) to suppress background signal in the images.
The dataset was acquired on a Bruker Avance-III HD 11.7 T spectrometer with a MIC-5 probe giving 3 T/m maximum gradient amplitude on-axis using a 10-mm diameter coil. The images were obtained using Bruker's multi-slice multi-echo (MSME) sequence customized to use double-rotation gradient waveforms (Jiang et al., 2023) running on Paravision 6.01. The multidimensional acquisition of the rat brain consisted of 1,095 images varying b-value (0.038–7.13·109 sm−2), centroid frequency ωcent/2π within 10–90% (29–110 Hz), normalized anisotropy bΔ (−0.5, 0, 0.5, and 1), orientation (Θ,Φ), recovery time τR (200–7,000 ms), and echo time τE (6–180 ms) in a total scanning time of 87 h. The whole acquisition protocol is shown in Supplementary Figure 1. The images were acquired at 18 ± 1 °C. The acquisition was a single slice in the coronal plane with a resolution of 70 × 70 × 250 μm3 (matrix size 122 × 128 × 1) located at approximately −3.5 mm from bregma.
The acquired images were preprocessed as follows: The data was reconstructed using Paravision 6.01 into the real and imaginary domains and denoised, using MPPCA approach included in DESIGNER toolbox (Chen et al., 2024). After denoising, the magnitude image was calculated and Gibbs-ringing artifact removal was applied with the sub-voxel shifts (Kellner et al., 2016) using MRtrix3 toolbox (Tournier et al., 2019).

2.1.2
In vivo rat brain at 7 T
One male Sprague-Dawley rat was used for the in vivo experiment (10 weeks old, 350 g, Harlan Netherlands B.V). It was housed individually in a cage with climate-control room under a 12 h/12 h light/dark cycle with ad libitum diet. The animal procedures were approved by the Animal Ethics Committee of the Province Government of Southern Finland and carried out according to the guidelines set by the European Community Council Directives 2010/63/EEC.
The rat was anesthetized under ~5% isoflurane in a mixture of 70% nitrogen/30% oxygen for 5 min, after which it was placed on the animal bed equipped with heated water circulation. The anesthesia then was maintained at ~2% isoflurane with the same nitrogen/oxygen mixture. The body temperature (36–37 °C) and respiration rate (35–60 breaths per min) were monitored during the scans.
The dataset was acquired on a 7–T/16-cm horizontal Bruker PharmaScan 7 Tesla preclinical system with a maximum gradient strength of 760 mT/m running on Paravision 6.01. An actively decoupled standard volume coil was used for RF-transmission and a quadrature surface coil for receiving. The images were obtained using Bruker's Spin Echo-Echo Planar Imaging (SE-EPI) sequence customized to use double-rotation gradient waveforms (Jiang et al., 2023; Yon et al., 2025). The images were acquired in a horizontal plane with a resolution of 250 × 250 × 700 μm3 with a slice gap of 300 μm and matrix size of 80 × 129 × 5. Double sampling and two segments were used. The images were acquired with two repetitions to include the reversed phase-encode images for later preprocessing steps. The multidimensional acquisition of the rat brain consisted of 389 pair of images varying b-value (0.031–2.65·109 sm−2), centroid frequency ωcent/2π within 10–90% (18–92 Hz), normalized anisotropy bΔ (−0.5, 0, 0.5, and 1), orientation (Θ,Φ), repetition time τR (1,000–3,500 ms), and echo time τE (21–77 ms) in a total scanning time of 1 h 9 min. The whole acquisition protocol is shown in Supplementary Figure 2.
The preprocessing of the in vivo rat images was as it follows: (1) The images were reconstructed in Paravision 6.01. (2) The images were separated into their respective phase-encoding directions using Matlab (Natick, Massachusetts: The MathWorks Inc.); (3) each set of phase-encoding directions were denoised separately using the MPPCA approach (Chen et al., 2024); (4) Gibbs rings artifact removal (Kellner et al., 2016; Tournier et al., 2019); and (5) correction of susceptibility-induced off-resonance field using FSL top up tool as described in (Yon et al., 2025).

2.1.3
In vivo human brain at 3 T
Six healthy volunteers were scanned (2 females, 4 males: age range 23–51 years) in compliance with the Finnish ethical permit 134/2023, 23.05.2023 and Finnish Medicines Agency (FIMEA) permit, FIMEA/2023/004206 in a 3-T MAGNETOM Vida Siemens scanner, equipped with a 70 cm bore size and 60/200 XT gradients. T1w images with voxel size of 1 × 1 × 1 mm3, repetition time τR = 2.3 s, and echo time τE = 1.9 ms were evaluated by an experienced radiologist, and no abnormalities were detected. The multidimensional images were obtained with a SE-EPI sequence customized to use diffusion general gradient waveforms (Martin et al., 2020). The data was acquired with a voxel size of 2 × 2 × 2 mm3 (matrix size 114 × 114 × 55 mm3). The multidimensional acquisition of the human brain consisted of 139 images varying b-value (0–3·109 sm−2), centroid frequency ωcent/2π within 10–90% (4–9 Hz), normalized anisotropy bΔ (−0.5, 0, and 1) using numerically optimized waveforms (Sjölund et al., 2015), orientation (Θ, Φ), repetition time τR (620–7,000 ms), and echo time τE (40–150 ms) in a total scanning time of 40 min. The whole acquisition protocol is shown in Supplementary Figure 3.
The preprocessing of the in vivo human images was as it follows: (1) The images were denoised using the MPPCA approach (Chen et al., 2024); (2) Gibbs rings artifact removal (Kellner et al., 2016; Tournier et al., 2019); and (3) correction of susceptibility-induced off-resonance field using FSL top up tool using a couple of non-diffusion weighted images with two different phase-encoding directions (Andersson et al., 2003).

2.2
Estimation of non-parametric D(ω)-R1-R2-distributions and visualization
Non-parametric D(ω)-R1-R2-distributions are derived through the approximation of the b(ω)-τR-τE encoded signal S according to the following equation:
The inversion of Equation 1 allows the acquisition of the ω-dependent self-diffusion tensor D(ω) by approximating Di(ω) as an axisymmetric tensor,
where R(θi, ϕi) is a rotation matrix and ω-dependent parallel and perpendicular eigenvalues, D||, i(ω) and D⊥, i(ω), are described by Lorentzian transitions (Narvaez et al., 2024) according to:
Using the Lorentzian approximation, each component in the D(ω)-R1-R2-distributions is described with its weight w and parameter set [D||, D⊥, θ, φ, D0, Γ||, Γ⊥, R1, R2]. We use the Monte Carlo inversion (Prange and Song, 2009) described in detail in Narvaez et al. (2022) and Narvaez et al. (2024) to retrieve the ensemble of solutions using Matlab (Natick, Massachusetts: The MathWorks Inc.) with the md-dmri toolbox (Nilsson et al., 2018).
The inversion limits used for the ex vivo rat brain at 11.7 T, in vivo rat at 7 T, and in vivo human brain at 3 T are presented in Table 1.
The Monte Carlo inversion settings for the three experiments were set to 20 steps of proliferation, 20 steps of mutation/extinction, 200 input components per step of proliferation and mutation/extinction, 10 output components, and bootstrapping by 100 repetitions using random sampling with replacement. As a result, the Monte Carlo inversion outputs an ensemble of 100 independent per-voxel solutions. Each of the solutions comprised the weights w and coordinates [D||, D⊥, θ, φ, D0, Γ||, Γ⊥, R1, R2] for 10 components within the pseudo-randomly sampled primary analysis space. We then evaluated D(ω) distributions at chosen values within 10 and 90 percentiles of the ω range probed by the gradient waveforms (Narvaez et al., 2022, 2024; Yon et al., 2025), giving ω-dependent distributions [D||(ω), D⊥(ω), θ, φ, R1, R2]. To follow well stablished metrics in the field (Pierpaoli et al., 1996), we obtained isotropic diffusivity Diso(ω), normalized anisotropy DΔ(ω), and squared normalized anisotropy (ω) through:
and
respectively.
In line with previous results from oscillating gradient encoding (Aggarwal et al., 2012, 2020; Xu et al., 2016; Aggarwal, 2020; Arbabi et al., 2020) and frequency-dependent multidimensional studies (Narvaez et al., 2022, 2024; Johnson et al., 2024; Manninen et al., 2024; Yon et al., 2025), restriction effect in ω-dependent Diso, and is quantified using the difference within the investigated frequency window from low to high ω.
The distributions can be visualized using contour plots in the 2D Diso-, Diso-R1, and Diso-R2 space, as well as bin-resolved signal fractions fbinn, means Ebinn[X], variances Vbinn[X], and covariances Cbinn[X,Y] appropriate to create per-voxel parameter maps (Topgaard, 2019). In this study, we used only per-voxel fbinn encoded as brightness and Ebinn[X] encoded as color for visualization, these were calculated as follows:
and

2.3
Unsupervised clustering of D(ω)-R1-R2-distributions
Derived from 100 bootstraps with 10 components each, we obtain a total of 1,000 per-voxel components. To consider the weights (w), each parameter value (Diso, DΔ or , R1, R2, and Δω/2πDiso) was replicated by its associated weight normalized and discretized between 0 and 100. Then, the resulted weighted parameters Diso, DΔ or , R1, R2, and Δω/2πDiso were organized into a 2D matrix and clustered using GMM approach, which is a probabilistic model that assumes that the data is a result of a mixture of normally distributed processes.
Using multidimensional normal distribution,
where μ is the mean of the distribution, Σ is the standard deviation of the distribution, and x is data, the GMM can be written as
where xn is a data sample, K is the number of clusters, N is the number of datapoints, and . The GMM was solved by an expectation maximization algorithm that maximizes the likelihood of a model given the number of mixture components k. The different dimensions in the bootstrapped data were normalized to the same range by dividing each variable value with the standard deviation of the variable.
The data was clustered in two different ways. In the first one, we only include the variables Diso and into the clustering using k = 3 to compare the results with the manual “3-bin” resolved fractions as used in previous studies (De Almeida Martins et al., 2020; Yon et al., 2020). For the second approach, we used Diso, DΔ, R1, R2 and Δω/2πDiso variables for the ex vivo and in vivo rat brain, and Diso, DΔ, R1, R2 for the in vivo human brain using k = 2–10. We selected DΔ for clustering to clearly differentiate the oblate distribution, facilitating the possible identification of inversion artifacts or biologically relevant features. As it is well known, unsupervised clustering methods rely in the selection of the k components. We calculated Schwarz's Bayesian Information Criterion (BIC) to aid the selection of the ideal k components using k = 2–15 (Supplementary Figure 4). The spatial information from the resulting clusters, clusn, can be visualized using the per-voxel normalized frequency, Fclusn. Using a similar formulation as bin-resolved approach, we obtained the cluster-resolved signal fractions fclus which can be visualized using contour plots in the 2D Diso-, Diso-R1, and Diso-R2 space and per-voxel means Eclusn[X] parametric maps calculated as follow:
and
In the in vivo human data, we computed the clusters for every subject independently to compare the cluster centroids and covariance matrices across solutions.

2.4
Tissue processing and histology
The hemisphere scanned ex vivo was washed 2 times for 25 min in 0.1 M PB pH 7.4 and then cryoprotected for 48 h in a solution containing 20% glycerol in 0.02 M potassium phosphate-buffered saline (KPBS) pH 7.4. After cryoprotection, the hemisphere was frozen in dry ice and stored at −70 °C until further processing. The hemisphere was coronally sectioned into 30-μm thick section and 1:2 series using sliding microtome. The first series of sections were stored in 10% formalin at room temperature and used for Nissl staining using thionin. The second series of sections were collected into tissue collection solution (30% ethylene glycerol, 25% glycerol in 0.05 M sodium phosphate buffer) and stored at −20 °C until staining with gold chloride for myelin (Laitinen et al., 2010).

2.5
MRI region-of-interest (ROI) and histology analyses
ROI analysis was performed on the ex vivo rat brain using the cluster fractions resulting from the full space clustering with k = 7. ROIs were selected in the auditory cortex, primary somatosensory cortex, granule cell layer of the dentate gyrus, corpus callosum and cingulum. Diso, DΔ, R1, and R2-distributions of fclusA, fclusB, fclusC, fclusD, and fclusG were obtained for each ROI. The projections were normalized by the maximum intensity within the cluster fractions. This normalization allowed to compare the cluster fraction contribution to each ROI independently of the ROI size.
We used a Nissl-stained section and a myelin-stained section from the same level in MRI to perform a comparison of cellular features with ROIs. Cells segmentation was performed on the Nissl-stained section, in which the per-cell area was calculated using QuPath v.0.5.0 software (Bankhead et al., 2017). Structure tensor (ST) analysis was calculated on the myelin staining (Budde and Frank, 2012; Molina et al., 2020), then, the photomicrograph was divided in patches of 70 × 70 μm2 to match with the in-plane voxel resolution as in ex vivo MRI. Angular central Gaussian (ACG) probability method was used to parametrize the results derived from the ST (Salo et al., 2021) and the coherency map ( ) was calculated (Püspöki et al., 2016). In addition, the optical density (OD) on the myelin staining was computed by manually selecting a background threshold and normalizing.

3
Results
Figures 1A, B shows the division of D(ω)–R1-R2-distributions into the manual three “bin” limits on the Diso- space, as in previous studies (De Almeida Martins and Topgaard, 2018; De Almeida Martins et al., 2020, 2021; Yon et al., 2020; Martin et al., 2021; Reymbaut et al., 2021; Narvaez et al., 2022, 2024). These bins can be considered as WM-like (bin1), GM-like (bin2), and FW-like (bin3). In the present study, the term WM to refer to such as the corpus callosum, cingulum, internal capsule, and optic tract (see the myelin-stained section in Figures 6A–C). Bin1 also included highly myelinated regions in the deep cortical and thalamic areas. Bin2 represented the remaining GM regions of the brain, while bin3 included the ventricles and liquid on the periphery of the brain.
Our first step was to compare the bin-resolved mean parametric maps, and the maps obtained with our unsupervised classification approach using k = 3. Figures 1A, B shows the per-voxel mean parametric maps of bin1 with characteristic low E[Diso], high E[], and both high E[R1] and E[R2] values in WM. bin2 had similar E[Diso] as in bin1, and lower E[], E[R1] and E[R2] than in bin1, as expect in GM. bin3 showed higher E[Diso] than in bin1 and bin2, while E[], E[R1], and E[R2] values were lower than in bin1 and bin2, corresponding to FW. Figures 1C, D display the cluster-resolved mean parametric maps using k = 3 in the Diso- space. The parametric maps of three cluster-fractions resemble the bin-resolved counterparts. The major difference between these two approaches was found in the FW fraction (Figure 1D). bin3 was only present in the ventricles and brain periphery, while clus3 showed the ventricles, brain periphery but also GM. In clus3, the highest E[Diso] values were found in the ventricles and the periphery of the brain. In GM, the E[Diso] values were lower than in clus3, but not as low as in clus2. The frequency-dependent parameter, Δω/2π E[Diso], also showed differences between bin- and cluster-resolved maps. While the highest Δω/2πE[Diso] values were mainly found in the granule cell layer of the dentate gyrus in the bin-resolved bin1 and bin2, these high values were found in the cluster-resolved clus2 and clus3 (white arrowhead in Figure 1).
Figure 2 displays the resulted normalized cluster frequency by setting k from 2 to 10 when applied the full Diso, DΔ, R1, R2 and Δω/2πDiso space. As showed in Figure 1D, using the Diso- space generated biological relevant fractions (i.e., a clear distinction of WM, GM, and FW regions) using k = 3; however, when using the full space required at least to k = 4 to get such biologically relevant fractions.
Using k = 4, FclusA (included 43% of the data) contains mainly GM, while major WM regions were in FclusB (29%) including the corpus callosum, external and internal capsule or optic tract. FclusC (16%) comprised FW regions in the ventricles and periphery of the brain, and the granule and pyramidal cell layers. FclusD (12%) contained WM areas, especially in the cingulum and internal capsule, as well as highly myelinated thalamic regions. The difference between k = 4 and k = 5 clusters was that FclusC in k = 4 split into 2 fractions in k = 5: FclusC displaying GM regions and FclusE mainly ventricles and periphery of the brain. clusD in k = 5 exclusively showed WM regions. With k = 6, an additional cluster, FclusF (2%), was identified throughout the entire brain. Using k = 7, we found the separation of FclusA (30%) into FclusG (19%) displaying GM regions. The clusters obtained with k = 8, 9, and 10 remained constant as with k = 7, with the addition of FclusH, FclusI, and FclusJ; however, these clusters contain small percentages of data of approximately 1–7%.
From k = 7 onwards, the clusters FclusA to FclusG remained similar across regions and percentages. For this reason, we selected the k = 7 cluster to visualize the parameter distributions that define each of the fractions, and obtained 2D contour plots between Diso-DΔ, Diso-R1, and Diso-R2 (Figure 3). For plotting the frequency dependance, we used the frequencies at 10, 50, and 90 % (ω/2π = 29, 69.5, and 110 Hz). The cluster fractions were ordered by amount of data classified into each fraction from largest to smallest.
fclusA, fclusG, and fclusC contained GM regions; however, those fractions were defined by different parameters (Figure 3A). fclusA mainly consisted of DΔ values close to zero; fclusG contained negative DΔ with the mode (Mo) = −0.46, and fclusC was characterized by variable DΔ at different frequencies ω/2π Hz. fclusB and fclusD were non-zero in WM regions, the main distinction between these two fractions was that R1 (Mo = 0.65 s−1) and R2 (Mo = 24.38 s−1) values in fclusB were lower than in fclusD. In fE, non-zero values were present in FW regions with low R1 and R2 values than previous cluster fractions. Diso values across cluster fractions showed that fclusA, fclusB, and fclusG had a similar Diso distributions (Mo = 0.27 × 10−9 m2s−1) (Figure 3B). fclusC showed higher Diso mode values at 0.62 × 10−9 m2s−1 at ω/2π = 110 Hz. FW cluster (fE) presented values over 1 × 10−9 m2s−1, while fclusD showed the lowest Diso values (Mo = 0.12 × 10−9 m2s−1) than the rest of the fractions.
Figure 4 shows all the parameter maps using the per-voxel cluster-resolved means in the k = 7 cluster fractions using the full space. We obtained the color-coded WM, GM, and FW-like fractions (Figure 4A) using the fclusB, fclusA and fE, respectively. This visualization is a representation of the already described distributions in Figure 3. Hence, we can observe in Figure 4B WM-like regions (clusB, and clusD) and GM-like regions (clusA, clusG, and clusC) in color-coded according to the per-voxel means. It is worth emphasizing that previously, in Figures 1B–D the high Δω/2πE[Diso] were located mainly in the granule cell layer of the dentate gyrus within other tissue or FW fractions. Here (Figure 4B), clusC is the cluster that contains high Δω/2πE[Diso] not only in the granule cell layer but also across other GM-like regions.
The cluster-resolved parameter maps E[Diso], E[], and E[R2] from the WM-like fractions clusB and clusD were compared with myelin staining from the same animal (Figures 5A–C). Although all are major WM fiber bundles, we observed differences in organization and density across the corpus callosum, cingulum, optic tract, and internal capsule. The optic tract demonstrated a more compact and dense organization than the corpus callosum. In contrast, the internal capsule exhibited complex crossing bundles, whereas to other regions like the corpus callosum were characterized by more uniformly parallel fiber organization (Figure 5C). The cingulum exhibited a rostro-caudal orientation, compared to the corpus callosum, which displayed an in-plane orientation in coronal sections (Figure 5C). Notably, although less dense than WM regions, certain GM regions—such as the thalamus, retrosplenial cortex, and deep cortical layers—exhibit high myelin content, as shown in Figure 5A. Regarding the parametric maps that can be comparable with the myelin staining by their presence in WM-like regions (Figure 5B), E[] in clustB showed similar values across the corpus callosum, optic tract and internal capsule, while the E[R2] parametric map in clusD showed higher values in the optic tract and internal capsule than the corpus callosum (Figure 5C).
The representative Nissl-stained section revealed regions with high cell density, such as the granule cell layer, CA3 and medial habenular nucleus (Figures 5D–F), which corresponded to high Δω/2πE[Diso] values in clusC (Figure 5F). In contrast, clusA exhibited limited differentiation between GM regions in the Δω/2πE[Diso] and E[Diso] parametric maps.
Figure 6 presents the results of ROI analyses for selected brain areas using both MRI and histology (Figure 6A), revealing distinct distribution patterns between WM and GM regions (Figure 6B). Here, the projections of the distributions were normalized by the maximum intensity within the cluster fractions. This normalization allowed to compare the proportion of data included in the cluster-fractions between ROIs independently of the number of voxels. fclusA distributions for all four parameters (Diso, DΔ, R1, and R2) displayed consistent parameter values and similar data proportions across the auditory cortex, primary somatosensory cortex, and granule cell layer. The most notable difference emerged in fclusC, in which the proportion of data was larger in the granule cell layer, followed by the auditory cortex and then, by the primary somatosensory cortex. Comparison of GM regions with Nissl staining (Figure 6C) revealed the highest cell density in the granule cell layer. Although segmentation of individual cell bodies in the granule cell layer was not successful, the auditory cortex was found to have a larger cellular area compared to the primary somatosensory cortex (Figure 6D), consistent with observations in fclusC. Small differences were noted between GM regions in fclusB and fclusG, where the auditory cortex showed a higher data proportion than primary somatosensory cortex and granule cell layer. As mentioned above, fclusB is associated with high myelin content. Within GM regions, the auditory cortex demonstrated higher density of myelinated axons, as showed by the myelin staining and optical density map (Figures 6E, F). The coherency map further revealed increased fiber orientation coherence toward the superficial cortical layers in both the auditory cortex and primary somatosensory cortex (Figure 6G).
The proportion of data within the corpus callosum and cingulum in fclusA, fclusG, and fclusC was relatively small as compared to GM areas (Figure 6B). However, among WM regions, the cingulum displayed a slightly higher proportion of data within those cluster-fractions than the corpus callosum. Both the corpus callosum and cingulum demonstrated lower cell density and smaller cellular area than GM regions (Figures 6C, D). In fclusB, the corpus callosum and cingulum showed consistent distributions across the DΔ, R1, and R2 parameters, with only small differences in Diso. Notably, the cingulum displayed a broader Diso distribution compared to the corpus callosum. The most prominent difference was observed in fclusD, where the cingulum had a higher proportion of data and broader R1 and R2 distributions than the corpus callosum. Both, the corpus callosum and cingulum contained similar high densities of myelinated axons, as confirmed by the myelin staining and optical density map (Figures 6E, F). However, the coherency map revealed lower coherence values in the cingulum compared to the corpus callosum (Figure 6G).
We obtained the per-voxel cluster frequencies from the in vivo rat brain using k = 2 to k =10 (Supplementary Figure 5). Compared to the ex vivo rat brain data, the in vivo data were characterized by lower signal-to-noise ratio (SNR), larger voxel size, and a different acquisition protocol. Consequently, the optimal k value for biologically relevant clustering differed between in vivo and ex vivo datasets. Despite these differences, the resulting cluster at k = 6 (Figure 7) remained similar in terms of the overall parameter distributions.
In the in vivo dataset (Figure 7A), components belonging to clusA, clusB, and clusC were primarily located in GM with smaller proportions in WM regions (33, 19, and 18%, respectively). clusD was observed in the ventricles and periphery of the brain (14%) while clusE was predominantly located in WM regions (13%). clusF was distributed across the entire brain but accounted for only a small portion of the clustered data (3%), consistent with findings in the ex vivo clusF results (Figure 3A). ClusD from ex vivo rat was not present in the in vivo rat clusters.
For plotting the frequency dependance, the centroid frequencies values at 10, 50 and 90% (ω/2π = 18, 52 and 92 Hz) were used. Contour plots for each fraction (Figure 7A) showed that fclusA was mainly characterized by the mode on DΔ = 0. The cluster fraction fclusB, had higher Diso and mode of DΔ = −0.5 compared to fclusA. A clear frequency dependance in both Diso and DΔ was displayed in fclusC. The cluster fraction considered as the FW-like fraction, fclusD, displayed the highest values of Diso > ~1.5 × 10−9 m2s−1 together with low R1, and R2. Regarding the WM-like fraction fE, it contained the highest DΔ > ~0.5 and high R2. Cluster-resolved parameter maps (Figure 7D) of fclusC and fE allowed to visualize that the highest values of Δω/2πE[Diso] in fclusC were in the cerebellum whereas fE exhibited high values of E[DΔ] and E[R2] in WM regions.
In the in vivo human brain, the parameters included into the clustering were Diso, DΔ, R1 and R2. The parameter Δω/2πDiso was excluded due to the limited frequency range available in the dataset (4–9 Hz). Per-voxel cluster frequency was obtained using k = 2 to k =10 (Supplementary Figure 6). The results from k = 6 were used to display the cluster fractions contour plots and cluster-resolved maps from a representative subject (Figure 8). Although the in vivo human data resulted in similar number of clusters as in vivo rat brain data, there is a clear absence of a frequency-dependent fraction. However, the clustering effectively delineated spherical and linear shapes as observed in the DΔ parameter. fclusA (22%) was primarily in GM regions such as cortex and deep GM characterized by low Diso and the mode on spherical shapes DΔ = 0. fclusB (21%) in WM regions had low Diso, DΔ > 0.5, and higher R2 than fclusB. fclusC (21%) was across WM and GM regions characterized by mode on DΔ = −0.5. fclusD was mainly located in ventricles and periphery of the brain (18%) with mode in Diso ≅ 3 × 10−9 m2s−1, while fE had a larger range of Diso between 1 × 10−9 m2s−1 and 4 × 10−9 m2s−1. Finally, fclusF (4%) was found in the periphery of the brain.
We computed the clusters for all the 6 human subjects to test repeatability (Figure 9). For all the subjects the clusters were calculated independently and with k = 2–10. We found that at k = 6, the per-voxel cluster frequencies were similar across subjects. Hence, we plot the clusters centroids and covariance matrices in the normalized Diso-DΔ space for two out of the six clusters (Figure 9A). The selected clusters to display were those in WM and GM regions. For the WM clusters, normalized DΔ across subjects were approximately on 1.6 with some variation of Diso around 0.7 and 0.9. The GM clusters displayed more consistent centroid on the Diso dimension around 0.65. We observed that the WM and GM cluster fractions (Figure 9B) have a similar mode across subjects. Additionally, the Diso and cluster-resolved maps of the WM and GM clusters are displayed from subject 1 (Figure 9C).

4
Discussion
In this study, we identified distinct water populations from non-parametric D(ω)-R1-R2-distributions using unsupervised clustering. Our method was evaluated in both in vivo and ex vivo datasets. The number of discernible clusters was influenced by the acquisition protocol, with broader relaxation and frequency-dependent parameters enabling more detailed clustering outcomes. Despite protocol differences, at least four biological relevant fractions—WM-like, two GM-like, and FW-like—were consistently observed across healthy human and rat in vivo brain data, as well as ex vivo rat brain. This flexible, data-driven approach offers broad potential for application to other tissues, including brain tumors, heart, prostate, and breast.
Rather than depending on single per-voxel trace and anisotropy values to separate water populations, bin division of non-parametric multidimensional diffusion-relaxation MRI uses the per-voxel Diso- distribution space to differentiate WM, GM and FW contribution to a voxel. Additionally, it has been observed the advantages of adding dimensions such as in (De Almeida Martins et al., 2020), in which the manual binning of the Diso--R2 space allowed the distinction of cortical GM and deep GM regions. This technique facilitates the differentiation of various water populations in a single voxel, as well as the assessment of their relative contributions. Nonetheless, bin limits are still selected manually and required adjustment depending on whether the study involves in vivo or ex vivo tissue, as well as if the data type varies from brain tissue.
Clustering methods address the limitations of manual bin selection by providing an unbiased, data-driven approach. Unsupervised methods using tensor distribution data have improved the characterization of diffuse gliomas (Song et al., 2022). Simultaneous diffusion and relaxation measurements enhance water population differentiation, aiding astrogliosis analysis (Benjamini et al., 2023) and cortical layers differentiation (Kundu et al., 2023). In addition to relaxation and diffusion parameters, our study incorporates restriction properties characterized by the frequency-dependent domain in the non-parametric multidimensional distributions. Considering variable diffusion time regimes enables a more thorough assessment of tissue microstructure. Previous work has demonstrated that in vivo rat GM can contain multiple water populations, each exhibiting distinct restriction properties, reflecting the heterogeneity of the tissue environment (Yon et al., 2024).
The incorporation of relaxation and frequency dependence significantly enhances the ability to resolve tissue heterogeneity within a voxel. Particularly on the relaxation side, variable TE allows a more efficient sampling of the multidimensional space and provides detail on the meso- and microscale tissue structures, while variable TR increases the sensitivity to the local chemical composition and component quantification (De Almeida Martins and Topgaard, 2018; Weiskopf et al., 2021). Furthermore, to increase the frequency dependence range is particularly important to include variable diffusion encoding waveform durations, and this is better achieved by variable TE. In this context, applying the unsupervised approach to the non-parametric D(ω)-R1-R2-distributions enables the extraction of per-voxel distributions of multiple water populations with greater biological relevance, as evidenced by their correspondence to histological counterparts. In contrast, relying on manual Diso and bin-limits in high-dimensional data can lead to the mixing of tissue fractions that are unlikely to biologically coexist within the same spatial domain. For example, in the manual bin approach, WM-like fraction (bin1) exhibited high Δω/2πDiso values within the granule cell layer. This ambiguity was resolved through unsupervised clustering, which correctly assigned the granule cell layer to a separate tissue fraction. A similar challenge arises in GM, where variations in cellular density, cell size, axonal content, and chemical composition differ across distinct brain regions. The unsupervised clustering approach provided characteristic cluster fraction distributions, offering a more accurate and biologically meaningful delineation of tissue heterogeneity.
We applied the unsupervised clustering approach to both in vivo and ex vivo datasets, revealing differences in the number of clusters identified between ex vivo and in vivo rat and in vivo human brains. These differences were primarily attributed to variations in acquisition protocols. The ex vivo settings allowed for higher spatial resolution and the collection of more extensive data due to the longer acquisition times, whereas in vivo acquisitions were constrained by shorter imaging sessions, MRI gradient strength, and motion-related limitations. Specially in clinical setups, tensor-valued diffusion encoding prohibits access to very short TE, decreasing the sensitivity of T2 sampling. Nevertheless, brain region differences in both R1 and R2 dimensions were still observed in our results showing the added value of those dimensions, even in degraded sampling conditions. Despite protocol-dependent differences, we consistently identified similar biologically relevant fractions across both ex vivo and in vivo datasets. A recent work demonstrated the reliability and reproducibility of non-parametric D(ω)-R1-R2-distributions in the human brain (Manninen et al., 2024). Additionally, our repeatability experiments showed that unsupervised clustering within the same protocol and acquisition settings results in similar centroids and covariance matrices.
Notably, certain components were relatively straightforward to interpret, such as the high DΔ distributions, which were predominantly associated with WM regions. In contrast, the zero and negative anisotropic distributions without frequency-dependence, were more challenging to interpret biologically. However, these components were notably more prominent in GM, suggesting underlying variations in cellular and microstructural properties. When frequency dependence data were available, as in the ex vivo and in vivo rat datasets, the clustering approach enabled the extraction of a fraction resembling the apparent intra-soma compartment. This fraction was characterized by increased Diso at higher frequencies, consistent with previous findings using oscillating gradient spin echo acquisitions (OGSE) (Does et al., 2003; Aggarwal et al., 2012; Baron and Beaulieu, 2014; Oshiro et al., 2024). Due to the limited frequency range in the in vivo human data, this intra-soma-like feature could not be resolved, this could be addressed by using clinical MRI scanner with stronger gradients (Tan et al., 2020; Hennel et al., 2021; Xu, 2021; Dai et al., 2023).
The clusters obtained in this study should not be interpreted as isolated features corresponding to a single microstructural property. For instance, a cluster characterized by high Δω/2πDiso may be attributed to the apparent soma fraction, analogously to the time-dependent microscopic kurtosis (μK) (Novello et al., 2022). However, it is essential to acknowledge that the pronounced restriction properties revealed by the frequency dependence of diffusion can also reflect underlying exchange effects. This phenomenon has been demonstrated by Chakwizira et al. (2023) who showed that the sensitivity to restriction or exchange can be modulated by varying the design of gradient waveforms. In the current study, we employed a broad spectrum of gradient waveforms, which inherently introduces the potential for confounding effects between restriction and exchange. To disentangle these contributions more effectively, a targeted approach could be considered, such as selecting a subset of the data with waveforms specifically optimized to enhance sensitivity to either restriction or exchange.
The inclusion of frequency dependence is crucial for probing the microstructure of GM microstructure. Specifically, for Δω/2πDΔ, we observed that increasing the frequency shifts the anisotropy distribution closer to 0. A tendency for FA to slightly decrease at lower diffusion times has been previously noted in time and frequency-dependent diffusion studies (Tan et al., 2020; Tétreault et al., 2020; Ba et al., 2024). However, in contrast with FA and μFA, DΔ distributions provides a more detailed representation of tensor shapes, distinguishing between prolate and oblate (Lawrenz et al., 2010; Shemesh et al., 2016b). In our observations, at low frequencies, the Δω/2πDΔ values in GM regions are predominantly closer to 1 (stick-like structures) and −0.5 (planar-like features). As the frequency increases, a distinct shift occurs in the distribution toward DΔ = 0 (spherical structures). This observation highlights the value of frequency dependence in identifying and characterizing specific microstructural compartments within GM. Investigating the soma more precisely—whether through biophysical modeling or signal representation—requires selecting an optimal frequency or time-dependent range. Intermediate ranges may incorporate contributions from other microstructural features, thereby complicating the interpretation and potentially leading to suboptimal signal fitting. There are only few studies about the time-dependent effect in the microscopic anisotropy (Ianuş et al., 2018; Narvaez et al., 2022, 2024; Johnson et al., 2024; Yon et al., 2025).
The use of Monte Carlo inversion algorithms to resolve multidimensional parameter spaces is a well-established method (Prange and Song, 2009), supported by prior studies and simulations (De Almeida Martins and Topgaard, 2018; Narvaez et al., 2024). This approach explores the entire solution space through multiple component decomposition, which, while comprehensive, can introduce challenges such as data overfitting and the generation of highly complex distributions. To address potential overfitting, the resulting solutions are typically weighted and summarized into statistical descriptors, such as medians, variances, and covariances, to enable the calculation of parametric maps (Yon et al., 2020; De Almeida Martins et al., 2021; Martin et al., 2021; Narvaez et al., 2022, 2024). In the current study, we leveraged the weighted solutions directly as input for clustering, acknowledging that not all data points contribute equally to the overall biological relevance of the solution space. By doing so, we aimed to retain the granularity of the multidimensional data while reducing the influence of less significant components. The Gaussian Mixture Model (GMM) first iteratively determines the centroids (μ) of the clusters, with each centroid representing a meaningful data point. In contrast, the covariance (σ) reflects the spread and shape of the clusters, which plays a key role in preventing data overfitting. When applied to non-parametric D(ω)-R1-R2-distributions, these GMM parameters correspond to the mean and variance of the distributions, both of which require careful interpretation. As the number of clusters increases within a given tissue, additional cluster fractions with less than approximately 8% of the data start to emerge. To avoid overfitting, we limited the number of clusters selected for each experiment. The validity of the identified cluster fractions can only be confirmed through spatial correlation with well-established biological structures, which necessitates histological comparisons for accurate biological interpretation.
In this study, we utilized two histological stainings, Nissl staining for cell bodies and myelin staining to identify myelinated fibers. While these stainings provide valuable insight into cellular density and myelin distribution, they do not capture the full complexity of the tissue microstructure. Brain tissue consists of diverse components, including various cell populations, dendrites, unmyelinated axons, vasculature, and extracellular matrix, which were not accounted for in this study, thus limiting the interpretation of the cluster fractions. Furthermore, while light microscopy is a well-established tool, its resolution is limited to visualize certain microstructural features, such as subcellular components and neuropil. These features are critical for understanding the full extent of tissue heterogeneity, particularly in GM regions with diverse cellular, axonal and synaptic characteristics. To strengthen the associations between MRI and histology, more comprehensive analyses incorporating additional stainings (e.g., for neuronal and glial markers, vasculature, or extracellular components) and advanced imaging techniques are needed. Electron microscopy (EM) could play a crucial role in this process, providing much higher resolution images that allow for visualization of ultrastructural features, such as synaptic junctions, fine axonal structures, and the organization of myelin at a sub-micrometer scale. Such insights would significantly enhance our understanding of the tissue microstructure and improve the biological interpretation of MRI data. Moreover, the integration of advanced image analysis techniques, including 3D reconstructions, machine learning-based segmentation, and quantitative morphometric analyses, could offer more detailed and precise evaluations of tissue organization (Abdollahzadeh et al., 2021; Larsen et al., 2021; Salo et al., 2021). Future work should thus incorporate these complementary methods to provide a more comprehensive view of tissue organization and its relationship to multidimensional MRI-based measures.

5
Conclusions and outlook
In conclusion, this study highlights the potential of unsupervised classification of multidimensional diffusion-relaxation correlation data to enhance tissue characterization at the voxel level. The method effectively disentangles tissue components and provides valuable insights into diffusivity, anisotropy, and chemical properties, which are crucial for identifying tissue changes in pathologies such as tumors, demyelination, Wallerian degeneration, and cell swelling. This approach offers significant advantages over traditional methods, enabling a more comprehensive understanding of tissue microstructure without relying on assumptions. While promising, further validation is needed in clinical settings, particularly with complementary histological data. Expanding this method to other tissues and diseases could broaden its applications, making it a powerful tool for both research and clinical diagnostics, especially in neuroimaging and other medical fields.