Multidimensional correlation of nuclear relaxation rates and diffusion tensors for model-free investigations of heterogeneous anisotropic porous materials

Despite their widespread use in non-invasive studies of porous materials, conventional MRI methods yield ambiguous results for microscopically heterogeneous materials such as brain tissue. While the forward link between microstructure and MRI observables is well understood, the inverse problem of separating the signal contributions from different microscopic pores is notoriously difficult. Here, we introduce an experimental protocol where heterogeneity is resolved by establishing 6D correlations between the individual values of isotropic diffusivity, diffusion anisotropy, orientation of the diffusion tensor, and relaxation rates of distinct populations. Such procedure renders the acquired signal highly specific to the sample’s microstructure, and allows characterization of the underlying pore space without prior assumptions on the number and nature of distinct microscopic environments. The experimental feasibility of the suggested method is demonstrated on a sample designed to mimic the properties of nerve tissue. If matched to the constraints of whole body scanners, this protocol could allow for the unconstrained determination of the different types of tissue that compose the living human brain.

Magnetic Resonance Imaging (MRI) has been firmly established as the method of choice for non-invasive investigations of the structure of the living human brain. In its simplest form, MRI yields maps of the 1 H spin density weighted by nuclear relaxation. The detected signal originates mainly from water because of its high concentration and favourable relaxation properties. The longitudinal relaxation rate R 1 of water is exquisitely sensitive to the local chemical composition of the tissue 1-6 while the transverse relaxation rate R 2 additionally carries information about micro-and mesoscale tissue structures [7][8][9] . The sensitivity of R 1 and R 2 to a plethora of tissue properties is a great advantage when aiming for contrast between tissue types and detection of subtle pathological changes, but it is also a challenge when attempting to interpret the observed contrast in terms of a specific property. Diffusion MRI, on the other hand, has a more direct relation between tissue microstructure and the measurable parameters, such as the diffusion coefficient D and, for anisotropic materials, the diffusion tensor D 10,11 . The method is particularly powerful for studies of white matter; being widely used to correlate microstructural changes with diseases 12,13 , normal brain development 14,15 , and learning 16 .
With a volume of approximately 1 mm3, most imaging voxels contain multiple types of tissue with varying diffusion and relaxation properties 17,18 , all of which contributing to the measured MRI signal. Because of this heterogeneity, the signal from each voxel is more accurately described with distributions rather than unique values of the observables. Encoding of the MRI signal for information about positions, displacements, and relaxation requires tens to hundreds of milliseconds, a time during which the observables are averaged by exchange processes and micrometre-scale displacements. The underlying microstructure is captured by a distribution of effective observables that is related to the signal data via a Laplace transformation. Estimating a distribution from discretely sampled and noisy data by inverse Laplace transformation is a highly challenging problem in numerical analysis and in general requires regularization to improve the stability of the solution at the expense of introducing inversion artefacts [19][20][21] . Despite these difficulties, continuous R 2 -distributions have proven useful for mapping the fraction of myelin water 4,5 . Multidimensional Laplace NMR [22][23][24][25][26][27] yields improved characterization of heterogeneous materials by estimating joint distributions of observables such as R 1 , R 2 , and D. In order to resolve a discrete component, it is sufficient that just one of the available observables has a value that is substantially different from the other components. Resolution in one of the dimensions allows for estimation of minor differences in the other observables that would be challenging to detect in the one-dimensional Laplace approach.
All multidimensional Laplace approaches to date have been limited to correlations between one or more relaxation rates and the scalar diffusion coefficient D rather than the diffusion tensor D. Diffusion encoding in a single direction convolves the effects of diffusion anisotropy and the orientation of the diffusion tensor eigenvectors, giving rise to D-distributions with complex shapes 28 that are impossible to reproduce with regularized Laplace inversion because of the well-known artefacts. Neglecting this fact when applying multidimensional Laplace methods to heterogeneous anisotropic materials, such as brain tissue, will invariably produce distributions fraught with artefacts that, although being numerically stable, are more related to the choice of regularization than the sought-for structure of the investigated material. Alternatively, the number of the components, as well as some of their properties, can be assumed using prior knowledge of the underlying structure, leaving a few parameters to be estimated by fitting the model to the experimental data 29 . Within the typical signal-to-noise ratio (SNR) in MRI, the same data can often be equally well described with fundamentally different models 30,31 , or even distinct parameter sets for the same model 32,33 . The validity of the analysis is difficult to ascertain since the results depend just as much on the assumptions as on the data.
We have recently used principles from solid-state nuclear magnetic resonance (NMR) spectroscopy 34 to design new diffusion MRI methods for investigating heterogeneous anisotropic materials [35][36][37][38][39] . The main insight from solid-state NMR is that the effects of diffusion tensor anisotropy and orientation can be disentangled by acquiring multidimensional data with correlations between isotropic and directional diffusion encoding 37,38 . Model-free and unconstrained inversion of the multidimensional data allows estimating a diffusion tensor distribution (DTD) [39][40][41][42][43] that is directly related to the distribution of microscopic environments in the material. Here, we join the DTD and Laplace frameworks into a single experimental method that combines the simple relation between diffusion tensors and microstructure with the sensitivity of R 1 and R 2 to the chemical composition of the tissue 1-6 . In this novel approach, a heterogeneous sample is characterized through correlations between the isotropic diffusivity D iso , normalized diffusion anisotropy D Δ , orientation of the diffusion tensor (θ, φ), and relaxation rates R 1 and R 2 .
The method at hand is tested with both simulations and experiments conducted on a phantom composed of an aligned liquid crystal and a yeast suspension, a schematic of which can be found elsewhere 37 . The used phantom mimics the diffusion properties of the nerve tissue model proposed by Stanisz et al. 29 ; the liquid crystalline and intracellular yeast components resemble water within axons and glial cells respectively, while the extracellular environment of the yeast suspension replicates the diffusion properties of nerve tissue's extracellular water. Spectroscopic proof-of-concept experiments were performed using Bruker microimaging hardware and a modified version of the NMR pulse sequence introduced in ref. 44 . While the experimental implementation presented here is too demanding for the hardware constraints of clinical MRI scanners, the proposed method can be redesigned using the smooth gradient waveforms from similar in vivo studies 40,45,46 . We consider this contribution as an important step towards the development of a protocol for the model-free quantification of the fractions of different structural components of the brain.

Theory
Diffusion MRI techniques are sensitive to molecular displacements at the 10-100 ms time-scale. Such sensitivity allows us to probe the underlying microstructure of a biological tissue by measuring its effects on the translational motion of water. The measured displacements can be approximated by an apparent diffusion tensor D that depends on both the geometry of the pore space 38,47 and the timing parameters of the used MRI sequence. Within the DTD model [39][40][41][42][43] , a heterogeneous material is pictured as a collection of independent microscopic domains whose structure is characterized by a corresponding microscopic D. Combining the DTD and Laplace NMR frameworks, the composition of a voxel can be characterized by the distribution P(D, R 1 , R 2 ) which is mapped into the signal amplitude S by a kernel K(τ R , τ e , b, R 1 , R 2 , D) according to where τ R and τ e are, respectively, the repetition and echo times, and b denotes the diffusion encoding tensor 40,48,49 , all of which being parameters under direct experimental control. S 0 symbolizes the signal obtained when no diffusion or relaxation encoding is applied. The "biophysical modelling" 29-33 approach of analysing MRI data can be seen as a post-acquisition procedure where a functional form of P(…) is assumed and fitted to the signal attenuation curve. Because the dimensionality and form of K(…) are controlled by the choice of pulse sequence and sampling scheme, the acquired data is intrinsically linked to the experimental protocol.
In typical relaxation and diffusion MRI protocols, an "encoding block", wherein R 1 , R 2 , and D modulate the magnetization, precedes a block where the signal is read out. The modulation details are contained within an appropriate kernel, whose functional form is imprinted into the signal via equation (1). At the micrometre length scale, where relevant cellular features appear, the effects of relaxation and diffusion on the detectable transverse magnetization are well described by the Bloch-Torrey equation 50,51 . While the general solution of the Bloch equations is highly complex, a simple exact expression can be attained if one neglects the presence of boundaries and exchange phenomena 52,53 . Being commonplace in studies of biological tissues and heterogeneous SCientifiC REPORTs | (2018) 8:2488 | DOI:10.1038/s41598-018-19826-9 materials 30-32,39-43 , such approximation is normally referred to as the Gaussian or free diffusion assumption. Its use yields the following kernel 22 : where b:D denotes a generalized scalar product defined as and −b D exp( : ) describe the weighting introduced by longitudinal recovery, transverse relaxation, and diffusion, respectively. The diffusion-encoding factor shows that D is measured relative to a b-tensor, which is implemented through a set of carefully chosen gradient waveforms. As demonstrated below via an explicit expansion of b:D, the parameterization of b can be manipulated to target specific diffusion features.
A symmetric positive definite second order tensor can be represented as a 3 × 3 matrix, and fully characterized by 6 independent elements [38][39][40] . Imposing an axial symmetry constraint, the number of independent values is reduced to four; two of them characterize the tensor orientation, and the other two define its size and shape parameters. In its principal axis system (PAS), an axial symmetric tensor Λ can be parameterized by its axial and radial eigenvalues λ and λ ⊥ , respectively. Alternatively, a complete description of Λ can be devised on the basis of its isotropic average λ iso and normalized anisotropy λ Δ 38,39,48 where λ is the trace of Λ. While the relative "size" of distinct tensors are represented by their corresponding λ iso values, λ Δ reports on the "shape" of Λ with λ Δ < 0, λ Δ = 0, and λ Δ > 0 describing oblate, spherical, and prolate tensors, respectively. In the diffusion MRI literature, both b and D are usually approximated as axisymmetric second order tensors 30-32,39-43 that can then be conveniently parameterized by the above metric. Following refs 35,38 : x y z with θ and φ denoting the polar and azimuthal angles that define the orientation of the D-tensor PAS in respect to the laboratory frame of reference. Using a similar parameterization for b, we can write 48,54 : where b corresponds to the trace of the b-tensor, and b Δ is its normalized anisotropy. The l i quantities are defined through equation (5), with the (θ, φ) variables replaced by (Θ, Φ), the polar and azimuthal angles describing a rotation of b through the laboratory frame. Combining equations (4) and (6), the double scalar product of equation (2) can be expanded as iso 2 where P 2 (x) = (3x 2 − 1)/2 is the 2 nd Legendre polynomial, and its argument can be rewritten as the cosine of the arc-angle β between the b-and diffusion tensors. Equation (7) highlights the effects of the b-tensor shape on the measured signal decay. A linear b-tensor (b Δ = 1) returns an effective diffusion coefficient D eff = bD iso [1 + 2D Δ P 2 (β)] defined by the size, anisotropy, and orientation of D, whereas spherical diffusion encoding (b Δ = 0) averages out the effects of anisotropy and orientation, thus allowing the direct measurement of D iso . Insertion of equations (7) and (2) into (1) yields ( , , , , , , , , , , , ) wherein size, anisotropy and orientation of the tensors are written as separate dimensions of the acquisition, (8) belongs to a class of integral transforms whose inversion is notoriously difficult, with an infinite number of distinct distributions P(R 1 , R 2 , D) being consistent with the experimental signal 55 . To perform an unconstrained estimation of the assumedly sparse distribution P(R 1 , R 2 , D iso , D Δ , θ, φ), enough information has to be encoded into the multidimensional signal in order to render S(τ R , τ e , b, b Δ , Θ, Φ) specific to the sample's microstructure. However, the high dimensionality of our acquisition space complicates the task of comprehensively sampling the (τ R , τ e , b, b Δ , Θ, Φ) set within a reasonable experimental time. Traditional Laplace NMR sampling schemes consist of rectangular grids wherein the amount of sampling is given by n 1 × n 2 × … × n i , with n i being the number of experimental points in the i-th dimension of the acquisition space. In our 6D space, the use of a grid scheme would then allow very low n i before the total number of samples becomes prohibitively high.
To circumvent the aforementioned problem we opt to follow the pseudo-random strategy shown in Fig. 1, which establishes 6D correlations across the (R 1 , R 2 , D iso , D Δ , θ, φ) space by randomly selecting combinations of (τ R , τ e , b, b Δ , Θ, Φ). Compared to the grid design, the pseudo-random approach allows a more thorough sampling of each element of (τ R , τ e , b, b Δ , Θ, Φ) for the same total number of acquired data points. The increase of sampling density per dimension, coupled to the introduction of a bias towards higher signal values, increases the sensitivity of S(τ R , τ e , b, b Δ , Θ, Φ) to the analysed material. The specificity of the acquired data, coupled with the assumption of a sparse analysis space, allows the model-free determination of P(R 1 , R 2 , D iso , D Δ , θ, φ) via an appropriate inversion algorithm 21,[56][57][58] . A comparison between the pseudo-random scheme and a grid-like sampling scheme is presented in the Supplementary Information. While the proposed acquisition scheme is reminiscent of the non-uniform sampling protocols from multidimensional NMR spectroscopy 59,60 or MRI fingerprinting 61 , the data analysis of this contribution follows a different methodology than the one of the cited works.
Before progressing any further, one should discuss the validity of our signal model. The derivation of Eq. (8) rests upon the two base assumptions of the DTD model: diffusion inside the different micro-domains is Gaussian and exchange of diffusing particles between the various domains is non-existent. When those assumptions are not met, the protocol measures an effective set of (R 1 , R 2 , D iso , D Δ , θ, φ)-values that are averaged by the water's residence time in the different domains and are highly dependent on the choice of experimental time parameters. Furthermore, neglecting the effects of restricted diffusion may result in the estimation of different apparent diffusivities for distinct values of b Δ (see Supplemental Material of ref. 37 ) or a bias in the estimation of domain sizes 62 . However, despite its limitations, the DTD model has shown to be a useful departure point for MRI studies of the human brain 40 .
A pulse sequence capable of encoding the NMR signal for both diffusion and nuclear relaxation is displayed in Fig. 2. The sequence is based on a previous triple-stimulated echo protocol 44 , to which an initial 90° pulse and subsequent recovery time were added in order to grant access to the full 6D acquisition space. Variation of τ R and τ e allows signal modulation by R 1 and R 2 . Diffusion encoding is performed by three consecutive pairs of bipolar gradient pulses, whose direction vectors n 1 , n 2 , and n 3 are equally distributed on the surface of a cone of aperture 2ζ. The angle ζ is directly related to the anisotropy of b through 44

Results
For the proof-of-concept experiments we used the NMR pulse sequence in Fig. 2 to acquire 1500 experimental points sampled according to the scheme of Fig. 1. The water resonance lines from the liquid crystal and the yeast suspension are narrow enough to be resolved in the chemical shift dimension (Fig. 3a). Hence, the signal decays from the sample's two constituents can either be separated or merged depending on which spectral regions are integrated. Figure 3c and d show that samples with distinct microstructure give rise to different signal patterns, thus highlighting that the signal pattern is specific to the underlying P(R 1 , R 2 , D iso , D Δ , θ, φ) function. Pulse sequence for encoding diffusion tensor D and nuclear relaxation effects. The thin vertical lines denote 90°x RF pulse while the thick vertical lines represent 180°y pulses. The sequence relies on the detection of a stimulated echo, whose magnitude is encoded for longitudinal recovery R 1 , transverse relaxation R 2 and D in three individual blocks. Relaxation encoding is done through the variation of the τ R and τ e delays, which weight the effects of R 1 and R 2 , respectively. The signal is further affected by R 2 and R 1 during the constant τ 1 and τ 2 delays, respectively. Diffusion is modulated by three sets of bipolar gradient pulses whose unit vectors n 1 , n 2 , and n 3 are shown in the right panel. The colour-code of the direction vectors is consistent with that of the pulse sequence brackets, showing that (n 1 , n 2 , n 3 ) have a three-fold azimuthal symmetry around an axis over which they hold a constant polar angle ζ. The b-tensor anisotropy is tuned by changing ζ, as described by Eq. (9).
(Adapted from ref. 44 . Copyright 2015 by Elsevier.) Capitalizing on the specificity of our acquisition protocol, the distribution P(R 1 , R 2 , D iso , D Δ , θ, φ) was retrieved from the measured data through a model-free inversion of Eq. (8). As elaborated upon in the Supplementary Information, the inversion procedure followed the same approach of our previous correlation works 37,39 . Briefly, Eq. (8) is inverted through a non-negative least squares fitting algorithm augmented with bootstrap 63 resampling. Our procedure allows the estimation of an ensemble of equally valid P(R 1 , R 2 , D iso , D Δ , θ, φ) solutions, the average of which is reported in Fig. 4 as a set of contour maps. Besides the non-negative constraint and a restriction on the range of allowed distributions, no further regularization constraints were used.
The probability spectrum displayed in Fig. 4 resolves the three microscopic components of the composite sample on a (R 1 , R 2 , D iso , D || /D ⊥ = (1 + 2D Δ )/(1−D Δ ), θ, φ) basis. The choice to parameterize anisotropy with D || /D ⊥ is justified by the fact that, unlike the D Δ parameter, D || /D ⊥ > 0 for any given "shape" and can thus be represented in a logarithmic scale. Yielding similar relaxation rates (log(R 1 /s −1 ) ≈ 0 and log(R 2 /s −1 ) ≈ 2), the intra-and extracellular yeast compartments are separated as two isotropic components with log(D iso /m 2 s −1 ) ≈ −9 and −11, respectively. The measurement of a relaxation time T 1 = 1/R 1 longer than the intracellular lifetime of yeast, τ ≈ 0.5 s 64 , indicates that the R 1 values obtained for both intra-and extracellular water are partially averaged by molecular exchange. The anisotropic population at [log(R 1 /s −1 ), log(R 2 /s −1 ), log(D iso /m 2 s −1 ), log(D || /D ⊥ )] ≈ (0, 1.5, −9, 2) originates from the liquid crystal. The orientation distribution function P(θ, φ) of the anisotropic component is naturally recovered and is represented as a 3D surface plot in Fig. 4c. Its shape and color-coding are consistent with an alignment along the z-direction of the laboratory frame, which coincides with the main magnetic field.
Besides the three expected components, the (D iso , D || /D ⊥ ) map also displays an oblate population with the same R 1 , R 2 , and D iso values as the fast isotropic component. Inconsistent with the phantom's design, the existence of this component can be attributed to an inversion artefact resulting from the under-sampling of b Δ at low b-values (see second plot of Fig. 1b). Since not enough "shape" information was encoded into the initial decay, the noise fluctuations accommodate an extra anisotropic component that bears no physical meaning. As demonstrated in Fig. 5, the presence of the oblate population can be corrected by increasing the b Δ range at lower b-values. While such modifications could be implemented without increasing the sampling density, we choose to keep such artefact for illustrative purposes. The spread of isotropic components over the D || /D ⊥ axis offers an additional example of artefacts resultant from the data noise 65,66 , the presence of which clearly shapes the correlation map without decreasing the fit quality (See Supplementary Figure S1).
The limits of our method were further assessed with data simulated for a discrete three-compartment voxel compatible with the used phantom. The artificial signal was generated using the pseudo-random acquisition protocol displayed in Fig. 1 and the (τ R , τ e , b, b Δ , Θ, Φ) set described in the Methods section. Gaussian distributed noise with an amplitude of 1/SNR was added to the data. Figure 6 shows the (D iso , D Δ ) maps estimated from the model free inversion of signals simulated using different sampling densities and SNR values. Although the entire P(R 1 , R 2 , D iso , D Δ , θ, φ) distribution was inverted, we only display the P(D iso , D Δ ) projection to facilitate the interpretation of the figure. The overall properties of the P(R 1 , R 2 , D iso , D Δ , θ, φ) recovered from experimental data can be see in simulations with the same sample size and an intermediate SNR. As expected, the inversion quality increases with increasing sample size and SNR as the amount of information encoded in the signal directly correlates with those two properties.
In the spirit of Prange and Song 66,67 , we used the variability between solutions to devise histograms showing the intrinsic variability of the resolved populations (see Fig. 7). The statistical analysis was carried out for four different sub-volumes of our analysis space, each of them containing one of the various components identified in Fig. 6. Focusing on the distribution of the weight of each sub-volume, P vol (R 1 , R 2 , D iso , D Δ , θ, φ), one notices that regions comprising higher diffusivities ( Fig. 7c and d) carry a higher uncertainty. Such behaviour can be interpreted as a symptom of the "shape" undersampling at low b and is at the origin of the oblate artefact, which disappears for higher sampling densities. While the solutions from the 500 points dataset completely fail to resolve the various populations, both 1500 and 5000 points are shown to return average values close to the underlying ground-truth. In fact, despite the presence of an artefact component, the average values estimated from the 1500 points dataset are close to or within a standard deviation from the "true" distribution. The exception lies on the calculated D Δ values, whose higher uncertainty lies in the non-trivial interdependence of this variable with the D iso -, θ-, and φdimensions.

Discussion
The link between the dimensionality of the correlation space and resolution power is well established, being discussed in various multidimensional NMR reviews [23][24][25] . Recently, we have also shown that the projection of D unto a sparse basis of (D iso , D Δ , θ, φ) allows us to resolve details that would otherwise be concealed by broad distributions of D eff wherein size, shape and orientation are intrinsically entangled [37][38][39] . Here, high dimensionality and sparsity are combined to endow the presented protocol with a much superior resolution power when compared to classical Laplace NMR methods, which build correlations in the non-sparse space of D eff rather than probing the separate components of the full diffusion tensor. Focusing on the correlation maps and their respective projections in Fig. 4b we notice that the three components cannot be resolved in the 1D D || /D ⊥ -, R 1 -and R 2dimensions. Since the various microenvironments of our colloidal phantom possess very similar R 1 and R 2 rates, these dimensions do not significantly contribute to the resolution of the different components. The components of this particular phantom are only resolved in the D iso , D || /D ⊥ , and (θ, φ) dimensions, which are separated in our novel acquisition scheme. Besides aiding in the separation of water components, the information gained by separating the various eigenvalues of D allows the detection of subtle differences in R 1 -and R 2 -values that would not be visible in standard 1D or 2D relaxation correlation experiments. For example, the 2D projection P(R 2 , D iso ) reports a slightly higher value of R 2 for the intracellular yeast domain compared to the extracellular one, a behaviour consistent with previous findings 25,68 .
The proposed method can be schematized into three separate steps: pulse sequence, acquisition protocol, and inversion algorithm, all of which can be independently matched to the experimental conditions at hand. For instance, the pulse sequence in Fig. 2, which is useful in the study of materials with high values of R 2 , can be substituted by a gradient modulation more compliant with the hardware constraints of clinical scanners. The family of continuous gradient waveform methods recently used for in vivo diffusion MRI studies 40,45,46,69,70 offers a suitable alternative. Similarly, the pseudo-random scheme could be replaced by a rectangular sampling grid, or P(R 1 , R 2 , D iso , D Δ , θ, φ) be retrieved through an inversion algorithm exploiting one of the many available regularization metrics 20,21,71,72 . The 6D correlation approach is also a natural template for the design of lower dimensional correlation experiments. For example, the experiments introduced in refs 37,39 can be respectively seen as projections onto the 4D (D iso , D Δ , θ, φ) and 2D (D iso , D Δ ) spaces. Just as well, any other projection space can be chosen, e.g. (R 1 , R 2 , D iso ), albeit at a cost of resolution.
Being an inherently ill-posed mathematical problem, the inversion of Eq. (8) accommodates an ensemble of solutions that are compatible with S(τ R , τ e , b, b Δ , Θ, Φ) within the measured SNR 63,66,67 . In order to circumvent the inversion's non-unique character, the solution space is typically constrained through the introduction of regularization. The most standard constraints are the consideration of a non-negative distribution function, the restriction of the analysis space to physically realistic ranges, and the introduction of a bias towards smooth distributions 58 . Although capable of mitigating the influence of noise or reducing over-fitting, common regularization strategies yield well-known artefacts that affect the shape of the retrieved distribution 22,25,58 . Hence, the solution's non-uniqueness and prevalence of artefacts should always be kept in mind when designing or using such experiments, a fact overlooked in contemporary MRI literature 73 . In an attempt to explicitly show the non-unique character, we use a bootstrapping procedure to estimate hundreds of possible solutions, and display a select number of them in the Supplementary Information. Examples of inversion artefacts were singled out from Fig. 4 in the Results section. The existence of those artefacts can motivate a more precise strategy for the determination and quantification of the different domains. For example, the signal could be fitted to a discrete model, constructed according to any prior knowledge about the studied material, and whose number of components and constraints should be compatible with the model-free results. In the present case, that means we do not expect oblate domains. Alternatively, one could introduce regularization constraints penalizing the appearance of additional peaks 20,21,71,72 or use a refined version of our inversion where all components are constrained to either isotropic or prolate geometries.
In summary, we have suggested and demonstrated an NMR protocol capable of characterizing a sample's heterogeneity on the basis of the observables defined by the Bloch-Torrey equations. Through an exhaustive sampling of a 6D acquisition space, the full P(R 1 , R 2 , D iso , D Δ , θ, φ) distribution is accessed in a model-free fashion. Since the typical MRI voxel comprises multiple microscopic domains with varying chemical and diffusion properties, the presented method shows great potential for in vivo studies of the human brain. In particular, through the imposition of physiologically reasonable constraints, we expect that our protocol can serve as a basis for experiments capable of determining the composition of a voxel in terms of tissue and cell types.

Methods
Phantom Preparation. The liquid crystal preparation followed the recipe described in a recent publication 74 , having the composition 41.94 wt% water (Milli-Q quality), 13.94 wt% of the hydrocarbon 2, 2, 4-trimethylpentane (Sigma-Aldrich, Sweden), and 44.12 wt% of the detergent sodium 1, 4-bis(2-ethylhexoxy)-1, 4-dioxobutane-2-sulfonate (trade name AOT from Sigma-Aldrich, Sweden). At room temperature, the liquid crystal is in a reverse 2D hexagonal phase wherein water diffuses along cylindrical channels with ~5 nm diameter which span over lengths of hundreds of micrometres, thus giving rise to highly anisotropic diffusion 74 . The sample was heated to 25 °C, and left in the NMR magnet over the weekend in order to align the crystallites parallel to the spectrometer's magnetic field 75,76 . A cube of fresh baker's yeast (Kronjäst AB) was acquired at a local supermarket, and mixed with tap water at a 1:1 volume ratio. The mixture was kept in its original container, and stored for 24 hours at room temperature in order to enable the sedimentation of yeast. The sediment was transferred to a 10 mm NMR tube by a syringe equipped with a 1 mm diameter needle. The resulting yeast suspension is composed of microscopic spherical cells with water populations in the intra-and extra-cellular domains. Both populations exhibit isotropic diffusion, but with a few orders of magnitude difference in apparent diffusivity 35,77 .
The flame sealed 5 mm tube containing the liquid crystal was inserted into the 10 mm tube containing the yeast cell suspension, with the larger tube being successively sealed with Parafilm tape. The composite phantom was centrifuged for 10 minutes at 223 g to create a denser packing of yeast cells at the bottom of the wider tube. Before starting the NMR measurements, the sample was placed inside the NMR spectrometer and left to equilibrate for 30 minutes at 18 °C. NMR Measurements. The experiments were conducted on a Bruker Avance-II 500 MHz spectrometer equipped with a MIC-5 probe fitted with a 10 mm RF insert. The temperature was kept at a constant value of 18 °C throughout the entire experiment time of 161 minutes. Measurements were performed with the pulse sequence of Fig. 2, using a delay τ 1 of 6.4 ms, a longitudinal storage time τ 2 equal to 64 ms, and a recycle delay of 0.1 s. Diffusion encoding was performed by symmetric trapezoidal gradients with a ramp time of 0.2 ms, a plateau time of 1.9 ms, and a peak amplitude of 0.6 T/m. The τ R time interval is preceded by a set of three saturation pulses designed to null the longitudinal magnetization. Coherence pathway selection was performed with a two-step RF pulse and receiver phase cycle coupled with extensive use of spoiler gradients 44 .
All data processing was carried out with in-house code written in Matlab (The Mathworks, Natick, MA). A more detailed account on the data inversion algorithm can be found in the Supplementary Information. Data availability. The dataset described in the current study are available from the corresponding author on reasonable request. MATLAB code for the data inversion is freely available at https://github.com/daniel-topgaard/ md-dmri.