Technical field

The present inventive concept relates to a method for quantifying isotropic diffusion and/or anisotropic diffusion in a sample.

Background

A wide range of porous materials, from lyotropic liquid crystals [1 ] to brain tissue [2], contain anisotropic pores with varying sizes, shapes, and degrees of alignment on mesoscopic length scales. A complete characterization of the material requires estimation of all these parameters, but unfortunately their effects on the detected MRI (Magnetic Resonance Imaging) signal are hopelessly entangled when using conventional diffusion MRI methods based on the Stejskal-Tanner sequence [3] with two magnetic field gradient pulses. This sequence may in the following be referred to as the single pulsed field gradient (sPFG) sequence or experiment.

In diffusion MRI (dMRI), each voxel (which typically may be of a millimeter-size) of the image contains information on the micrometer-scale translational displacements of the water [15]. sPFG is used in diffusion tensor imaging (DTI), enabling quantification of mean diffusion (MD, also apparent diffusion coefficient, ADC) and diffusion anisotropy (Fractional Anisotropy, FA). Although sPFG-based DTI measures are very sensitive to changes in the cellular architecture, sPFG generally provides robust estimations only in highly organized white matter bundles. In less ordered tissue, it may provide little insight into the nature of that change, leading to common

misinterpretations. For example, changes in FA are thought to represent white matter integrity, however, many factors (cell death, edema, gliosis, inflammation, change in myelination, increase in connectivity of crossing fibers, increase in extracellular or intracellular water, etc) may cause changes in FA. The limited specificity of measures such as FA and MD hinders our ability to relate the measurements to neuropathologies or to local anatomical changes such as differences in connectivity [24, 25, 26, 27]. In contrast to sPFG, non-conventional dMRI sequences can begin to bridge between the macro and micro levels of scale in the brain by providing information about distributions of cellular shapes, sizes and membrane properties within a voxel.

Building on the formal analogy between the chemical shift and diffusion anisotropy tensors, it has been shown that solid-state NMR (Nuclear

Magnetic Resonance) techniques, such as "magic-angle spinning", can be adapted to diffusion MRI [4]. In its simplest form, magic-angle spinning of the g-vector allows for estimation of the distribution of isotropic diffusivities free from the confounding influence of anisotropy.

WO 2013/165312 discloses how isotropic diffusion weighting of a diffusion weighted echo signal attenuation may be achieved by a continuous or discrete modulation of the dephasing vector q(t) such that an anisotropic contribution to the echo signal is minimized, for example by employing magic- angle spinning. WO 2013/165313 discloses a method for quantifying microscopic diffusion anisotropy and/or mean diffusivity by analysis of echo attenuation curves acquired with two different gradient modulation schemes, wherein one gradient modulation scheme is based on isotropic diffusion weighting and the other gradient modulation scheme is based on non- isotropic diffusion weighting. WO 2013/165313 discloses that non-isotropic diffusion weighting may be achieved for example using single-pulse gradient spin echo (PGSE).

Although these prior art methods enable separation of isotropic and anisotropic contributions to the echo signal attenuation and quantification of inter alia microscopic fractional anisotropy, it would in some cases be desirable to have a greater freedom in terms of the gradient modulation schemes used for causing the diffusion weighting and still be able to analyze and quantify microstructure properties such as microscopic diffusion anisotropy and/or mean diffusivity e.g. for the purpose of tissue

characterization using diffusion spectroscopy. For example isotropic diffusion encoding may in some cases impose high requirements on the hardware with respect to slew rate and maximum magnitude which are difficult to meet with older and less expensive equipment.

Summary of the inventive concept

An objective of the present inventive concept is to provide a method for quantifying isotropic diffusion and/or anisotropic diffusion in a sample which does not require the use of diffusion encoding magnetic gradient pulse sequences causing isotropic diffusion encoding. Further objectives may be understood from the following summary of the inventive concept.

According to an aspect of the present inventive concept there is provided a method for quantifying isotropic diffusion and/or anisotropic diffusion in a sample, the method comprising:

performing diffusion weighted magnetic resonance measurements on the sample using diffusion encoding magnetic gradient pulse sequences

wherein each magnetic gradient pulse sequence Gj is generated such that a diffusion encoding tensor bj for the magnetic gradient pulse sequence Gj has one to three non-zero eigenvalues, where b, = i q, (t)q (t)dt , q, (t) is

The method further comprises collecting data representing magnetic resonance echo signals resulting from said measurements on the sample, wherein at least a subset of said data represents echo signals being acquired with a set of magnetic gradient pulse sequences causing anisotropic diffusion weighting, and wherein the diffusion encoding tensor for each gradient pulse sequence of said set of magnetic gradient pulse sequences has three nonzero eigenvalues, at least one of the three eigenvalues being different from the other two eigenvalues. The method further comprises calculating a degree of isotropic diffusion and/or a degree of anisotropic diffusion using said data.

The inventive method is, among others, based on the insight that performing diffusion encoding using a magnetic gradient pulse sequence (or shorter "pulse sequence") such that a diffusion encoding tensor has three non-zero eigenvalues, at least one different from the others, and causing anisotropic diffusion weighting makes it possible to control the effect of diffusion anisotropy in a sample material on the echo signal. As will be described in greater detail below, this enables accurate characterization of microscopic diffusion properties (such as isotropic and anisotropic diffusion), in particular of microscopic diffusion properties of microscopic compartments within the sample which are smaller than the spatial resolution of the magnetic resonance measurements. Moreover, such characterization is enabled while not being reliant on the use of isotropic diffusion weighting as in the prior art (which may be achieved by diffusion weighting tensors having three non-zero eigenvalues which all are equal). This may enable accurate diffusion measurements to be performed on a greater range of equipment.

According to the inventive method, each magnetic gradient pulse sequence Gj is generated such that a diffusion encoding tensor bj for the magnetic gradient pulse sequence Gj has 1 to 3 non-zero eigenvalues. In other words each magnetic gradient pulse sequence Gj is generated such that there is a diffusion encoding tensor representation bj of the pulse sequence Gj which has 1 to 3 non-zero eigenvalues. Analogously, for each magnetic gradient pulse sequence of the above-mentioned set of magnetic gradient pulse sequences causing anisotropic diffusion weighting, there is a diffusion encoding tensor representation which has 3 non-zero eigenvalues, at least one being different from the other two eigenvalues. The set of magnetic gradient pulse sequence causing anisotropic diffusion weighting may form at least a subset of the diffusion encoding magnetic gradient pulse sequences The at least one eigenvalue of the diffusion encoding tensor for each of said gradient pulse sequence causing anisotropic diffusion weighting may advantageously differ by at least 5%, and even more preferably by at least 10%, from any of the other two eigenvalues. This may ensure a sufficient degree of anisotropic diffusion weighting in the sample, facilitate subsequent calculations and reduce the hardware requirements.

Said subset of data may represent echo signals acquired from a same portion of the sample, the portion including a plurality of partial volumes presenting different degrees of isotropic diffusion or different degrees and/or orientations of anisotropic diffusion, wherein the calculation of a degree of isotropic diffusion and/or a degree of anisotropic diffusion may include calculation of an estimate of a degree of isotropic diffusion and/or an estimate of a degree of anisotropic diffusion for at least one of said partial volumes.

Especially, said portion may have a spatial extension which matches a spatial resolution of the diffusion weighted magnetic resonance

measurements. Hence, each one of the partial volumes may have an extension which is less than the spatial resolution. Such partial volumes may in the following be referred to as "microscopic partial volumes". Thus, in the above and in the following, the degree of isotropic diffusion and/or the degree of anisotropic diffusion calculated using said data may be referred to as a degree of isotropic diffusion and/or a degree of anisotropic diffusion for a sub- resolution or "microscopic" partial volume of the sample. A diffusion of each partial volume may have a diffusion tensor representation D. In other words each partial volume may have a diffusion which is definable by a respective diffusion tensor D. Thus, within the portion may be represented by a distribution (for example Gaussian distribution) of diffusion tensors D.

Preferably, a plurality of diffusion weighted magnetic resonance measurements may be performed on the sample. At least two, preferably a plurality, of the diffusion encoding magnetic gradient pulse sequences have tensor representations bj having one to three non-zero eigenvalues. At least two, preferably a plurality, of the diffusion encoding magnetic gradient pulse sequences are different from each other.

According to one embodiment said set of magnetic gradient pulse sequences causing anisotropic diffusion weighting forms a first set of magnetic gradient pulse sequences and said subset of data forms a first subset of data representing a first echo attenuation curve acquired with the first set of magnetic gradient pulse sequences, and wherein said data further includes at least a second subset of data representing a second echo attenuation curve acquired with a second set of magnetic gradient pulse sequences causing isotropic or anisotropic diffusion weighting. Thereby, isotropic and/or anisotropic diffusion may be quantified based on echo signals representing two different echo attenuation curves. It should be noted that the determinations "first" and "second" merely is to be construed as labels for the respective sets of magnetic gradient pulse sequences and subsets of data. They do not necessarily imply any particular ordering, i.e. that the first set of pulse sequences are applied to the sample prior to the second set of pulse sequences. Indeed, they may be applied in the reverse order or even in an arbitrary interleaved manner.

Each pulse sequence of the first set may be generated such that a first eigenvalue and a second eigenvalue of the diffusion encoding tensor for said pulse sequence are equal to each other. The third eigenvalue of a diffusion encoding tensor of the first set is different from the first and second

eigenvalues (thereby causing anisotropic diffusion weighting). Moreover, each pulse sequence of the second set may be such that a first and a second eigenvalue of the diffusion encoding tensor for said pulse sequence are equal to each other. Using different magnitudes of the first and the second eigenvalues of a diffusion encoding tensor translates to varying the effect diffusion anisotropy in the sample material will have on a resulting echo signal. Hence, generating magnetic gradient pulse sequences in this manner makes it possible to probe the diffusion properties of the sample with respect to isotropy and anisotropy.

The pulse sequences of the first set and the pulse sequences of the second set of magnetic gradient pulse may have varying maximum gradient magnitudes. With respect to the diffusion weighting tensor, this may be expressed as a trace of the diffusion encoding tensor for a pulse sequence of the first (or second) set varies throughout the first (or second) set. Thereby, the strength of the diffusion weighting may be varied.

According to one embodiment, for each pulse sequence of the first set there is a first diffusion encoding tensor invariant Δ definable by:

PAS

where b _xx represents the first eigenvalue of the diffusion encoding tensor for said pulse sequence, b _yy ^PAS represents the second eigenvalue of the diffusion encoding tensor for said pulse sequence and b _zz ^PAS represents the third eigenvalue of the diffusion encoding tensor for said pulse sequence, and

wherein the first set of pulse sequences is generated such that the first diffusion encoding tensor invariant A _{b 1} of the pulse sequences of the first set are equal to each other. By controlling the generation of the diffusion encoding magnetic gradient pulse sequences of the first set in this manner, the first echo attenuation curve represented by the first subset of data may represent an echo attenuation curve acquired using diffusion encoding tensors having a same degree of anisotropy A _{b l} .

Analogously, for each pulse sequence of the second set there may be a second diffusion encoding tensor invariant A _{b 2} definable by:

where b _xx ^PAS represents a first eigenvalue of the diffusion encoding tensor for said pulse sequence, b _yy ^PAS represents the second eigenvalue of the diffusion encoding tensor for said pulse sequence and b _zz ^PAS represents a third eigenvalue of the diffusion encoding tensor for said pulse sequence, and wherein the second set of pulse sequences is such that the second diffusion encoding tensor invariant A _{b 2} of the pulse sequences of the second set are equal to each other and A _{b 2} is different from A _{b 1} .

According to one embodiment calculating a degree of isotropic diffusion and/or a degree of anisotropic diffusion includes:

calculating a degree of isotropic diffusion and/or a degree of anisotropic diffusion by analyzing a change, a variation or a difference between a first echo signal acquired with a pulse sequence of the first set and a second echo signal acquired with a pulse sequence of the second set. Since the pulse sequences of the first set and the pulse sequences of the second set may be generated to present different degrees of anisotropy (i.e. Δ and A _{b 2} ), a change, a variation or a difference (e.g. with respect to amplitude) between the first and the second echo signals enables estimation of the degree of isotropic and/or degree of anisotropic diffusion. To simplify the estimation, the first and the second echo signals may be acquired with gradient pulses of equal maximum gradient magnitude (in other words equal values of the diffusion weighting magnitude b).

Within the scope of the inventive method, said data may, in addition to the first subset of said data and the second subset of said data, include at least a third subset of data representing a third echo attenuation curve acquired with a third set of magnetic gradient pulse sequences causing anisotropic diffusion weighting,

wherein the diffusion encoding tensor for each gradient pulse sequence of the third set has 3 non-zero eigenvalues of which a first eigenvalue and a second eigenvalue are equal to each other and different from a third eigenvalue, and

wherein, for each pulse sequence of the third set there is a third diffusion encoding tensor invariant Δ _{Α 3} definable by:

where b _xx ^PAS represents a first eigenvalue of the diffusion encoding tensor for said pulse sequence, b _yy ^PAS represents the second eigenvalue of the diffusion encoding tensor for said pulse sequence and b _zz ^PAS represents a third eigenvalue of the diffusion encoding tensor for said pulse sequence, and wherein the third set of pulse sequences is such that the third diffusion encoding tensor invariant Δ _{Α 3} of the pulse sequences of the third set are equal to each other and Δ _{Α 3} is different from A _{b 2} and A _{A l} . Acquiring data representing further echo attenuation curves along further "lines of constant encoding tensor anisotropy (e.g. Δ _{Α 3} ) enables an extended probing of the diffusion properties of the sample.

According to one embodiment each pulse sequence of the first set is such that A _{b 1} > 0 , each pulse sequence of the second set is such that

A _{b 2} = 0 , thereby causing isotropic diffusion weighting, and each pulse sequence of the third set is such that A _{b 3} < 0. This enables the "shape" of the diffusion characteristics of the sample to be estimated and analyzed. In particular, it becomes possible to estimate whether the diffusion of is mainly isotropic (i.e. spherical), mainly unidirectional (i.e. oblate) or mainly planar (i.e. prolate).

According to one embodiment the method further comprises

calculating, based on the data (which according to the above may include at least the first, the second and the third subset of data) representing said echo signals, a probability distribution indicating a probability of each one of said echo signals being associated with each one of a plurality of different values of a model isotropic diffusion parameter D, _so and/or a model anisotropic diffusion parameter A _D . This enables, among others, analysis of

measurement results obtained during measurements on a sample including domains presenting different degrees of isotropic and/or anisotropic diffusion. From the probability distribution the number of such domains (i.e.

components) may be identified.

The probability distribution may be calculated by determining a numerical solution to a system of equations relating the echo signals represented by said data to a product of a kernel function and said probability distribution. The system of equations may especially be a linear system of equations.

The probability distribution may be a joint probability distribution p and the kernel function may be a matrix K including at least M x N elements, each of said elements being based on an integration of

where A = 3bD _1S0 A _b A _D , for a combination of values of a diffusion weighting magnitude b, a diffusion encoding tensor invariant A _b , the model isotropic diffusion parameter D, _so and the model anisotropic diffusion parameter A _D . The elements of the matrix may be calculated for different combinations of values of the diffusion weighting magnitude b, the diffusion encoding tensor invariant A _h , the model isotropic diffusion parameter D, _so and the model anisotropic diffusion parameter A _D .

According to one embodiment the method may further comprise:

applying each pulse sequence of said first set of magnetic gradient pulse sequences a plurality of times to the sample, with different orientations of the gradient pulse with respect to a fixed laboratory frame, and forming said first subset of data by averaging echo signal measurements acquired for the different orientations. This may be referred to as "powder averaging" whereby, in cases where there is some preferential alignment of domain orientations, it is possible to mimic the effects of random domain orientations. Such "powder averaging" may be performed also for the second set of magnetic gradient pulse sequences. Namely, by applying each pulse sequence of said second set of magnetic gradient pulse sequences a plurality of times to the sample, with different orientations of the gradient pulse with respect to a fixed laboratory frame, and forming said second subset of data by averaging echo signal measurements acquired for the different

orientations.

According to one embodiment each one of said diffusion encoding magnetic gradient pulse sequences Gj forms part of a (separate) triple stimulated echo sequence. This may be particularly advantageous when performing measurements on a sample including material domains with comparably short transverse relaxation time T ₂ .

According to one embodiment the method further comprises:

forming a system of equations based on an expansion of a function relating an echo signal E to a diffusion encoding tensor b and a diffusion tensor D,

calculating an average diffusion tensor < D > and a diffusion tensor covariance tensor S by determining a solution to the system of equations using echo signal measurements represented by said data and

representations of at least a subset of the diffusion encoding tensors bj,

calculating an invariant bulk component Sbuik of the covariance tensor S by projecting S_onto a bulk basis E _bu ik,

calculating an invariant shear component S _S hear of the covariance tensor S by projecting S_onto a shear basis E _S hear, and

calculating a degree of isotropic diffusion and/or a degree of anisotropic diffusion using the invariant bulk component Sbuik and/or the invariant shear component S _S hear-

This embodiment enables quantifying differences in microstructure diffusion properties, for example in terms of microscopic diffusion anisotropy, without requiring axially symmetric diffusion encoding tensors. Especially, the invariant bulk component Sbuik may form an estimate of a variation of degrees of isotropic diffusion between the above-mentioned different partial volumes of the portion of the sample. The invariant shear component S _S hear may form an estimate of a variation of directions of anisotropic diffusion between the above-mentioned different partial volumes of the portion of the sample.

The calculated average diffusion tensor < D > may be calculated for the above-mentioned portion including a plurality of partial volumes. < D > may represent an estimate of an average diffusion tensor for the portion. Similarly the diffusion tensor covariance tensor S may represent an estimate of a covariance of the distribution of diffusion tensors for the portion.

The system of equations may be a linear system of equations, wherein the echo signal measurements represented by the data may be used to form constants of the linear system of equations and said at least a subset of the diffusion encoding tensors may be used to form parameters of the linear system of equations. Especially, the system of equations may be equivalent to a cumulant expansion of the function i?(b) = (exp( - < b,D>)) .

The degree of anisotropic diffusion may be calculated based on a sum of the invariant shear component S _S hear and a projection of the square of the average diffusion tensor < D > onto the shear basis E _S hear- Especially, the degree of anisotropic diffusion may be calculated based on a variance of said sum.

The degree of calculated anisotropic diffusion may further be based on a ratio between a projection of the square of the average diffusion tensor < D > onto the bulk basis E _bu ik and said sum. Especially, the degree of anisotropic diffusion may be calculated as an estimate of a microscopic fractional anisotropy ^FA based on said ratio.

The projection of S_onto the bulk basis E _bu ik, may be calculated by calculating an inner product between a matrix representation of the covariance tensor S_and a matrix representation of the bulk basis E _bu ik- The projection of S_onto the shear basis E _S hear, may be calculated by calculating an inner product between the matrix representation of the covariance tensor S and a matrix representation of the shear basis E _S hear-

The projection of the square of < D > onto the bulk basis E _bu iik, may be calculated by calculating an inner product between a matrix representation of the square of < D > and the matrix representation of the bulk basis E _bu ik-

The microscopic fractional anisotropy ^FA may in particular be calculated as

where (D ^¾2 ) = S+(D)'

The inventive method and the above-disclosed embodiments thereof may be methods for diffusion MRI, wherein the calculated degree of isotropic diffusion and/or degree of anisotropic diffusion may be used as a contrast parameter for a voxel of diffusion MRI data. Brief description of the drawings

The above, as well as additional objects, features and advantages of the present inventive concept, will be better understood through the following illustrative and non-limiting detailed description of preferred embodiments of the present inventive concept, with reference to the appended drawings. In the drawings like reference numerals will be used for like elements unless stated otherwise.

Fig. 1 illustrates some representative examples of gradient waveforms. Fig. 2 illustrates a theoretical echo attenuation signal S{b, A _b ) vs. JbD _iS o. Figs 3-4 illustrate a number of experimental results. Detailed description of preferred embodiments

To facilitate understanding of the present inventive concept, a discussion of some theoretical concepts will now be provided with reference to the drawings. Anisotropic Gaussian diffusion

The directionality of Gaussian diffusion processes is captured in the diffusion tensor D [11]. In its principal axis system (PAS), the tensor is diagonal with the elements D _XX ^PAS , D _yy ^PAS , and D _ZZ ^PAS . In the context of this disclosure, it is convenient to characterize the diffusion tensor with its isotropic value D,.

anisotropy AD, and asymmetry η ₀ ^' .

D l.so = (D x ^r x ^AS + D y ^r y ^AS + D z ^r z ^AS ) jj/3

1 D + D PAS

Δ PAS

D. yy

zz

^D 3D. (1)

Q ^FAS — D ^FAS

γ. _ yy xx_

2D l.so Δ D ₀

It may be noted that this formalism is reminiscent of the one used for the chemical shift tensor in solid-state NMR (e.g. as in [12] and [13]).

The elements are ordered according to the convention | D _zz ^PAS -Di _S0 | >

| Dyy ^PAS -Djso| > | Dxx ^PAS -Diso| . The numerical factors in Eq. (1 ) are selected to get parameters in the ranges -1 /2 < A _D ≤ 1 and 0 < η ₀ < 1■ The tensor is axially symmetric when η ₀ = 0. Positive and negative values of A _D correspond to prolate and oblate tensor shapes, respectively. The convention for ordering the elements assures that the z-axis is the main axis of symmetry for both prolate and oblate tensors. It should however be noted that other conventions (e.g. using the x- or y-axis as the symmetry axis) may be used without departing from the scope of the inventive concept.

A general rotation of the PAS from the lab-frame using the Euler rotation matrices R _z (a), R _y (P), and R _z (/) gives the following zz-element of the lab-frame diffusion tensor:

£) ^PAS _|_ £) ^PAS D ^PAS - D ^PAS

D _zz (a, p) = D ^? cos ² β + ^{w M} sin ² β ^» ** sin ² β cos2a (2)

Using the relations in Eq. (1 ), this expression can be rearranged to

D _2Z (a, p) = D _1S0 [l + A _D (2P ₂ (cosp) -¾ sin ² p cos2a)] (3) where P ₂ (x) = (3x ² -1 )/2 is the second Legendre polynomial. Eq. (3) reduces to

¾β) = ₈₀ [ΐ + 2Δ ₀ Ρ ₂ (οο _δ β)] (4) for axially symmetric tensors, η ₀ = 0. As a consistency check, it may be noted that insertion of Eq. (1 ) into Eq. (4) gives D _zz (0) = D _ZZ ^PAS and

Ο _ζζ (π/2) = (D _XX ^PAS + D _yy ^PAS )/2 as expected. Diffusion-weighting tensor b

The NMR signal is encoded with information about translational motion using a time-varying magnetic field gradient G ⁷ (f) = {G _x {t), G _y {t), G _z {t)}. The

instantaneous dephasing vector q(f) is given by the time integral (5) where γ is the magnetogyric ratio of the studied nucleus. During a

measurement, the echo signal may be recorded at the echo time r when the spin magnetization is rephased, i.e. q(z) = 0, in other words at the time instant t = r at which the spin magnetization is rephased. Assuming Gaussian diffusion, the signal amplitude S can be written as (6) where So is the signal intensity at zero gradient amplitude (i.e. the non- diffusion weighted echo signal) and b:D denotes a generalized scalar product defined as

b : D = y∑ .D..

(7)

The diffusion-weighting matrix b is given by b = } q(t) q ^T (t)dt

0

In analogy with Eq. (1 ), the Jb-matrix can be characterized with the total diffusion weighting b, anisotropy A _b , and asymmetry η^.

b = b b i

"-¾ ^rA _

( ^S +,

"y ^r y ^AS + , b i

^z ^r z

bA _b

It may be noted that the definition of the Jb-matrix through Eqs. (5) to (8) can be found in standard text books on diffusion NMR and MRI, see, e.g., chapter 4.4.1 in Price [14] or chapter 9.7.2 in Callaghan [15]. However,

characterization of the Jb-matrix using solid-state NMR terminology is novel and, as shown below, both simplifies the notation and provides a framework for designing measurement protocols and analysis methods. The elements of the Jb-matrix transform under rotation according to the rules for rank-2 tensors [16]. In the following, the Jb-matrix (i.e. b) may hence be referred to as the b- tensor. Variable-angle spinning of the q-vector

In spherical coordinates, the g-vector may be defined by its inclination ζ(ΐ), azimuth y/(f), and magnitude qF(t), where q is the maximum magnitude and F(t) the time-dependent magnitude normalized to the interval 0 < F(t)≤ 1 . The Cartesian components may be obtained from the relation

= qF (t) {sinC (t) cos (t), sinC (t) sin (t), cosC (t)} Insertion of Eq. (10) into Eq. (8) and application of standard trigonometric relations gives the following expressions for the b-tensor elements:

| 2 (,) 1 - cos ² ζ it) + sin ² ζ (t)cos[2y/ (t)}

q ²

o ^

K ₌ f ) _F > (,) 1 - co ¹ C ft) - sin ' C (Qcos[2y ^/ ¾

0 ^

As shown in Eriksson et al. [4], all terms containing ψ(ί) vanish if ^ is constant and the trajectory of the g-vector has at least three-fold symmetry, i.e.

, ητ 2πη

ψ(ί) =ψ (0) +ψ t+—

V N J N

(12)

F(t) = F

N J

Λ/ > 3, n = \,2,..., N

Another way of nulling the terms with y/(t) in Eq. (1 1 ) is if the trajectory fulfills the condition (see reference [4])

where n is an integer other than zero and f _d is an effective diffusion time g ^' by

Geometrically, this modulation of ψ(ί) and F(t) corresponds to a g-vector spinning about the z-axis at an angular velocity dy/(f)/df proportional to F(t) while following a path on the surface of a cone with aperture 2ζ.

Explicit evaluation of the b-tensor elements in Eq. (1 1 ) gives

2 l - cos ²

b xx = b yy = q ± ² t d. ^ —

K = q cos ζ

h , = h y„x = b„ = b which upon insertion into Eq. (9) yields

A _b = P ₂ (cos^ ) (16)

Assuming that the trajectory of the g-vector obeys Eqs. (12) or (13), the b- tensor is axially symmetric, with the z-axis as the main axis of symmetry, and has an anisotropy that can be tailored by varying the angle ζ from

conventional single-directional diffusion-weighting at ζ = 0, via isotropic diffusion-weighting at the magic-angle ζ = acos(1/3 ^{1 2} ) (i.e. using the terminology from the reference [4]), to so-called circular encoding at ζ = π/2 (i.e. using the terminology from the reference [17]). As long as the time- modulations of y/(f) and F(t) remain the same, adjustment of ζ only affects Δ*, without influencing the values of q, f _d , or b.

Numerical optimization of the gradient modulation functions

For a given modulation F(t) of the q-vector, its time-dependent orientation is given by the angle y/ (t), obtained via the integral in Eq. (13), and the selected constant value of ζ. The Cartesian components of the g-vector may be calculated with Eq. (10) and the gradient modulation functions are given by the derivative

A conceptually simple modulation function is F(t) = 1 in the interval 0 < t≤ τ and F(t) = 0 otherwise, corresponding to infinitely short and strong gradient pulses at t = 0 and τ. For practical implementation on MRI hardware with limited gradient capabilities it is necessary to have less abrupt transitions between F(t) = 0 and 1 . One possible procedure for finding optimal gradient waveforms for g-MAS diffusion weighting on clinical MRI scanners is described in [10]. In brief, F(t) may be expanded as

where the coefficients a _m are optimized iteratively to get maximum diffusion weighting within the constraints of a given waveform duration τ and maximum gradient amplitude G on each of the three gradient channels. The details of the optimization routine are described in the reference [10]. A final result of a joint optimization for the cone angles ζ = 0 and ζ = π/2, yielding one axial and two radial gradient modulation functions that can be superposed to give g- vector modulations at arbitrary cone angles and orientations in the lab frame are listed in Table 1 .

The coefficients a _m in Table 1 may be referred to as coefficients for a numerically optimized g-VAS gradient waveform. Using Table 1 , explicit gradient modulation functions may be obtained by selecting values of the echo time τ as well as the g-vector inclination and magnitude g, and subsequently calculating the normalized g-magnitude modulation F(t) with Eq. (18). The effective diffusion time f _d may be calculated with Eq. (14). The g- vector azimuthal angle y/(f) may be calculated with Eq. (13). The Cartesian components of the g-vector may be calculated with Eq. (10). Finally the

Cartesian components of the gradient vector may be calculated with Eq. (17).

Some representative examples of gradient waveforms are shown in Fig. 1 , including the corresponding g-vector trajectories and Mensor elements. As the azimuthal angle ψ{ί) of the g-vector is varied as a function of time the modulations exemplified in Fig. 1 may be referred to as variable- angle spinning of the g-vector (g-VAS). Rows 1 -4 of Fig. 1 correspond to the angles ζ = 0°, 35.3°, 54.7°, and 90°, as indicated in column 1 of Fig. 1 .

Column 2 of Fig. 1 shows gradient modulation functions G _x (t), G _y (t), and G _z {t) (dotted line, dashed line, and full line, respectively) obtained from the coefficients listed in Table 1. Column 3 in Fig. 1 shows corresponding q- vector modulation functions q _x (t) (dotted line), q _y (t) (dashed line), q _z (t) (full line), and q(t) (dash-dotted line). Column 4 of Fig. 1 shows a 3D plot of the q- vector trajectory (black line) in relation to the positive x-, y-, and z-axes. As may be seen the q-vector trajectories in the second and third row lie on the surface of a cone with aperture 2ζ. The g-vector trajectory in the first row is aligned with the z-axis. The q-vector trajectory in the fourth row lies in the xy- plane. Column 5 of Fig. 1 shows the corresponding ensor anisotropies A _b = 1 , 0.5, 0, and -0.5. m a _m l 10 ^"3

1 345.5666

2 276.0000

3 186.3600

4 102.7100

5 28.8030

6 13.0610

7 3.8476

8 5.0764

9 7.8710

10 10.9790

11 7.2733

12 4.6788

13 2.0443

14 0.2220

15 1.8141

16 0.0001

17 1.0098

18 1.2173

19 0.8931

20 0.7950

21 -0.2223 Table 1

Effective diffusion coefficient D _zz ^eff

Insertion of Eq. (15) into Eq. (7) gives

which can be rewritten as

D +A D - D ) (20) using the relations in Eq. (16) and the rotational invariance of the trace of the diffusion tensor D, _so = (D _xx + D _yy + D _zz )/3. Assuming axial symmetry of the diffusion tensor, insertion of Eq. (4) yields

The factors following the Jb-value can be interpreted as an effective diffusion coefficient D _zz ^eff (/¾ that depends on the anisotropy of the Mensor as well as the orientation of the diffusion tensor in the lab frame through the angle according to

Comparing Eqs. (4) and (22) shows that the effect of the inherent diffusion anisotropy is scaled by the value of A _b , which in turn depends on the angle ζ between the z-axis and the spinning g-vector. The values of D _zz ^eff are in the range between D _zz ^eff (0) = = Ζ _ΒΟ (1-Δ/ ο).

Powder-averaged signal attenuation and effective diffusivity distribution Consider a macroscopic sample consisting of an ensemble of randomly oriented microscopic anisotropic domains having the same values of D _IS0 and AD. In cases where there is some preferential alignment of domain

orientations, it is possible to mimic the effects random domain orientations by "powder-averaging" the data, i.e. record data for a series of directions of the symmetry axis of the g-trajectory and subsequently average the results over the various directions. In a measurement on the macroscopic sample encoded by a magnetic gradient pulse sequence represented by a g-vector having the z-axis as symmetry axis (i.e. 0) each domain gives rise to a signal that may be calculated by inserting Eq. (21 ) into Eq. (6)

Integrating the contributions from all the domains gives = 03)exp{- bD _lS0 [l + 2A _b A _D P ₂ (cos 0)])dj3

where Ρβ(β) is the angular distribution function, normalized in the interval 0 < β≤ π/2. A random distribution of domain orientations corresponds to Ρρ(β) = sin/?, which yield

= 3bD- upon evaluation of the integral in Eq. (24). In Eq. (25), T(S,X) is the lower incomplete gamma function. This function may for example be conveniently evaluated numerically with the "gammainc" function in Matlab. It should be noted that although the gamma and square-root factors in Eq. (25) are imaginary when the argument is negative, their ratio remains real and positive. The y(s,x) factor could also be written in terms of the error function "erf", utilizing the fact that y(1 /2,x) oc erf(x ^{1 2} ).

Fig. 2 illustrates the theoretical signal S{b, A _b ) vs. JbD _iS o according to Eq. (25), where b is the diffusion-weighting magnitude and D _IS0 the isotropic diffusivity, and A _b A _D , where Ab and A _D are the anisotropies of the diffusion- weighting and diffusion tensors b and D, respectively. The surface is calculated with Eq. (25) and the parameters b, D _IS0 , A _b and A _D are defined from the respective tensor eigenvalues in Eqs. (1 ) and (9).

The expression in Eq. (25)provides a basis for analysis of experimental data. When Ab = 0, Eq. (25) reduces to a single-exponential decay ^ * _Ρ (- _6¾0 ) ₍₂₆₎

Mi

thus providing a simple way of extracting the isotropic part of the diffusion tensor without the confounding effects of anisotropy. When b is held constant at some finite value on the order 1 /Dj _S0 , the value of A _D can be determined from the characteristic variation of S as a function of Ab.

The multi-exponential signal decay in Eq. (25) can be interpreted as being the Laplace transformation of a distribution of effective diffusivities P(D _zz ^eff ) according to

where

2j3D _lsa A _b A _D [D ? - D l - A _b A _D )\ in the range of D _zz ^eTT from -Δ/,Δο) to D _iso (1 +2Δ _ό Δο) and P(D _zz ^eTT ) = 0 otherwise. The distribution has a singularity at D _ZZ ^EFF = Ο _Ϊ50 (1 -Δ*,ΔΟ),

corresponding to the domain orientation β = π/2. Eq. (28) is analogous to the "powder-pattern" NMR spectrum obtained for an axially symmetric chemical shift anisotropy tensor. The mean value of the distribution is D _IS0 , while the 2 ^r and 3 ^rd central moments, μ ₂ and are

and

respectively. Signal intensity for a multi-domain material

A material consisting of a collection of domains with different values of D _IS0 and A _D gives rise to a powder-averaged signal that can be expressed as the integral transform P(D _lso ,A _D )dD _iso dA _D (31) where Ρ(Οί ₅₀ ,Δο) is the 2D joint probability distribution of D _IS0 and AD. The kernel K(b,A _b ,Di _S0 ,AD) is given by the right-hand side of Eq. (25) and maps the 2D "analysis-space" (Οί ₅₀ ,Δο) onto the 2D "acquisition-space" {b,A _B ).

Estimating Ρ(Οί ₅₀ ,Δο) from an experimental set of data l{b,A _B ) may be considered as an ill-posed problem and may benefit from special procedures to ensure numerically stability. Based on the approaches used for tackling similar problems in NMR diffusion and relaxation correlation methods [18-20], with some additional inspiration from compressed sensing [8, 9], the procedure outlined below may be used. Eq. (31 ) can be discretized and written in matrix form as

S= Kp (32) where s is a column vector of signal amplitudes measured for a M

combinations of (b,Ab), p is the sought-for column vector of probabilities for N discrete pairs of (Ο _Ϊ50 ,ΔΟ), and K is the matrix version of the kernel

K(b,Ab,Di _S0 ,AD) calculated for an MxN grid of {b,Ab) and (Οί ₅₀ ,Δο) pairs. Eq. (32) is a set of linear equations and can in principle be solved by straightforward matrix inversion if the problem is overdetermined, i.e M > N.

Unfortunately, the "smoothness" of the kernel renders this direct approach extremely sensitive to experimental noise, with wild fluctuations of the solution vector p for minor changes in the input data vector s. Assuming that the material consists of a few discrete components, it may be advantageous to look for a solution that is sparse, meaning that most of the elements of p are zero.

Based on the considerations above, a sparse solution p that is consistent with experimental data s may be estimated by minimizing the function

N M M

^ (P) =∑∑| / - ^ + ^λ ∑| ; I (33) where the 1 ^st term is the least-squares misfit and the 2 ^nd term is the /-i-norm weighted by the regularization parameter λ. Combined with non-negativity constraint on the elements of p, Eq. (33) can be formulated as a quadratic programming problem

/ ^: (p) = p ⁷ K ⁷ Kp- 2s ⁷ Kp+s ⁷ s+ 1 ⁷ p (34) which is readily solved by, e.g. , the "quadprog" function in the optimization toolbox of Matlab.

Example experiments

In the present and the subsequent section, a number of examples of proof-of- principle experiments will be described as well as results thereof. According to these examples, experiments were carried out on lyotropic liquid crystals with the detergent Aerosol-OT and water consisting of an equimolar mixture of H ₂ 0 and D ₂ 0. Based on the equilibrium phase diagram (see e.g. reference [21 ]), the detergent concentration was chosen to give three different liquid crystalline phases: lamellar (25 and 75 wt%), bicontinuous cubic (80 wt%), and reverse hexagonal (85 wt%). The samples were initially weighed into 10 ml_ vials, allowing for thorough mixing, and subsequently 400 μΙ_ was transferred to 5 mm disposable NMR tubes. The phase symmetry was independently verified by recording small-angle x-ray scattering patterns and ² H NMR spectra. A sample with two distinct diffusion tensor components was prepared by inserting the 5 mm NMR tube with the 25 wt% Aerosol-OT sample into a 10 mm NMR tube with decanol.

Diffusion magnetic resonance experiments were performed on a Bruker Avance-ll 500 MHz spectrometer with an 1 1 .7 T magnet equipped with a Bruker MIC-5 microimaging probe giving maximum magnetic field gradients of 3 T/m in three orthogonal directions. The g-VAS gradient modulation was implemented by including the waveforms shown in Fig. 1 on both sides of the 180° RF pulse in a standard ¹ H spin echo pulse sequence. The signal recorded during the second half of the spin echo yields a high-resolution spectrum after Fourier transformation, thus permitting the separation of the water ¹ H signal from the ones originating from the detergent. For studies of systems with short transverse relaxation time T ₂ , it is advantageous to implement the variation of the diffusion-weighting anisotropy A _b by adjusting the directions of the pulsed field gradients in a triple-stimulated echo sequence. The three gradient directions have the azimuth angle ψ = 0°, 120°, and 240° and the inclination ζ, which gives Δ*, through Eq. (16).

The 75, 80, and 85 wt% Aerosol-OT/water samples were investigated with the triple-stimulated echo implementation using three pairs of gradient pulses with duration δ = 1 ms and leading-edge separation Δ = 100 ms. A rectangular grid of the acquisition space (b,Ab) was sampled by varying the amplitude of the gradient pulses and the angle ζ. The maximum gradient amplitude was on the order of 1 T/m and adjusted for the different samples to get approximately the same maximum signal attenuation. Both the cubic and reverse hexagonal samples give sufficiently narrow Aerosol-OT resonance lines to be detectable with the triple-stimulated sequence. Empirically, it was found that the water and Aerosol-OT resonance lines overlap at 25 °C, but that the overlap is rendered insignificant by increasing the temperature to 80 °C. For consistency, the 75, 80, and 85 wt% samples were all studied at 80 °C. The samples remain in the same liquid crystalline phase as at 25 °C according to the equilibrium phase diagram and confirmed by ² H

spectroscopy measurements. The 25 wt% Aerosol-OT/water/decanol sample was studied at 25 °C with the spin echo version at a gradient modulation duration τ = 140 ms and maximum gradient amplitude 0.090 T/m. The (jb,A _/ ,)-space was sampled in a zig-zag pattern by varying the maximum gradient amplitude and ζ. Although the resonance lines could be separated by Fourier transformation, they were recorded jointly to give a signal containing multiple components with different diffusion behavior. In practice, only water and the terminal methyl group of the decanol have sufficiently long T ₂ to survive the lengthy spin echo sequence.

In order to assure that the data corresponds to a random distribution of domain orientations, as required by Eq. (25), the acquisition was repeated and averaged for 31 different "cone orientations", i.e. orientations of the main symmetry axis of the g-vector trajectory. These directions were chosen according to the electrostatic repulsion scheme (see reference [22] and [23]). Results of the example experiments

Fig. 3 shows experimental data for lamellar, bicontinuous cubic, and reverse hexagonal liquid crystalline phases. The data was fitted with Eq. (25), yielding values of the diffusion anisotropy A _D consistent with the known

microstructure. It is noteworthy that the sign of A _D can be extracted from the characteristic variation of the signal as a function of Δ*, as long as JbD _iS o is on the order of unity and above.

In more detail, Fig. 3 shows data representing the measured signal attenuations S{b, A _b ) vs. the diffusion-weighting magnitude b and anisotropy A _b for AOT/water liquid crystals of the (a) lamellar, (b) bicontinuous cubic, and (c) reverse hexagonal types. The top row shows schematic illustrations of these types. These geometries characterize the respective water

compartment geometries on the length scale of tens of nanometers. The filled circles represent experimental data points sampled on a rectangular grid in the (jb,A _/ ,)-space. The grid illustrates a fit of Eq. (25) to the experimental data using the initial signal intensity So, the isotropic diffusivity D _IS0 , and the diffusion anisotropy A _D as adjustable parameters. The fit yields D _iS o 10 ^"9 m ² s ^{" 1} = 3.53 (lamellar), 2.37 (cubic), and 1 .22 (reverse hexagonal), as well as A _D = -0.38 (lamellar), 0.00 (cubic), and 0.80 (reverse hexagonal).

The results for the liquid crystal/decanol sample are indicated in Fig. 4. The zig-zag sampling of the (jb,A _/ ,)-space facilitates display of the data in a 2D plot with signal vs. b. In such a plot, the presence of multiple components with different D _IS0 can be discerned as curvature in the lower envelope of the data, while non-zeros values of A _D result in oscillations at a frequency given by the sampling pattern. The amplitude of the oscillations is related to the magnitude of AD, while the ratio between the local maxima corresponding to Δ*, = 1 and - 0.5 gives the sign of AD.

The S{b,A _B ) data is converted to a probability distribution Ρ(Οί ₅₀ ,Δο) through the ^-regularized model-free approach as described above. The resulting distribution contains components at (D _IS0 = 10 ^"10 m ² /s, A _D = 0) and (Djso = 0 ^"9 m ² /s, A _D = -0.5), corresponding to decanol and water, respectively. Once it is known that the distribution contains two components, the

coordinates and amplitudes of these are more accurately determined by a two-component model fit where the signal from each of the components is described by Eq. (25). The result of such a fit is also shown in Fig. 4b. The obtained results are consistent with the known isotropic diffusion of decanol and the lamellar symmetry of the liquid crystalline phase. In more detail Fig. 4a illustrates the experimental water and decanol signal S{b, A _B ) vs. the diffusion-weighting magnitude b for a tube-in-tube sample with an AOT/water lamellar liquid crystal (inner tube) and decanol (outer tube). The zig-zag sampling pattern of the (jb,A _/ ,)-space is shown on top. Fig. 4b shows the probability density Ρ(Οί ₅₀ ,Δο) consistent with the S{b, A _B ) data. The contour lines show the result of a /i-regularized model-free estimate, while the crosses indicate the results of a two-component model fit giving

Diso 0 ^"9 m ² s ^"1 = 0.083 (decanol) and 1 .33 (water), as well as Δ ₀ = 0.00 (decanol) and -0.496 (water). The points in Fig. 4a indicate the experimental data while the black line represents the two-component model fit.

Description of embodiments

In accordance with the present inventive concept there is provided a method for quantifying isotropic diffusion and/or anisotropic diffusion in a sample. With reference to the preceding description, the isotropic diffusion may for example be quantified by the isotropic value D, _so for a diffusion tensor D (as defined in connection with Eq. 1 above). The anisotropic diffusion may for example be quantified by the anisotropy A _D for the diffusion tensor D.

The various calculation performed in the method may for example be implemented using a set of software instructions which may be stored on or embodied on a non-transitory computer storage medium.

The method may be performed using a state-of-the-art NMR

spectrometer or MRI device. As is well-known in the art, such devices may include one or more processors for controlling the operation of the device, inter alia the generation of the magnetic gradient pulse sequences, the acquisition of echo signals as well as sampling and digitizing the measured signals for forming data representing the acquired echo signals. The generation of the diffusion encoding magnetic gradient pulse sequences may be implemented using software instructions which may be stored on a computer readable media (e.g. on a non-transitory computer readable storage medium) and be executed by the one or more processors of the device. The software instructions may for example be stored in a program/control section of a memory of the device, to which the one or more processors of the device has access. Collected data representing the measurements may be stored in a data memory of the device, or of a computer or the like which may be connected to the device. The calculations of the method may be implemented software instructions which may be stored on a computer readable media and be executed by the one or more processors of the device. However it is equally possible to carry out the calculations on a device which is separate from the NMR spectrometer or MRI device, for example on a computer. The device and the computer may for example be arranged to communicate via a communication network such as a LAN/WLAN or via some other serial or parallel communication interface. It should further be noted that, instead of using software instructions, the operation of the method may be implemented in dedicated circuitry of the device such in one or more integrated circuits, in one or more application-specific integrated circuits (ASICs) or field- programmable gate arrays (FPGAs), to name a few examples.

The method comprises performing diffusion weighted magnetic resonance measurements on the sample using diffusion encoding magnetic gradient pulse sequences which may be denoted ■ With reference to the preceding description, each magnetic gradient pulse sequence Gj is definable by, or has a representation in the form of, a diffusion encoding

τ

tensor bj having three non-zero eigenvalues, where b, = J q, (t)q (t)dt

0

(Equation (5)) and q, (t) is the dephasing vector which is proportional to {t')dt' (Equation (8)).

For the purpose of acquiring echo signals, each diffusion encoding magnetic gradient pulse sequences may be supplemented with one or imaging magnetic gradients and optionally magnetic gradient correction gradients, as is well-known in the art. Hence each diffusion weighted magnetic resonance measurement may be performed using a magnetic pulse gradient sequence including the diffusion encoding magnetic gradient pulse sequence Gi, an imaging magnetic gradient sequence and optionally a correction magnetic gradient sequence. In some cases these sequences may overlap in time. However, even in such a case at least a part of the sequence may be described or represented by a diffusion encoding tensor bj having the above-mentioned properties.

The method further comprises collecting data representing the echo signal measurements. The data collection may include sampling and digitizing an echo signal received from the portion of interest of the sample. In line with the above, each echo attenuation measurement may be in the form of a echo signal amplitude S, which, assuming Gaussian diffusion, has a dependence on the diffusion encoding tensor bj and D as given by Eq. (6). The

measurements may be spin echo measurements or stimulated echo measurements, for example triple stimulated echo measurements.

The portion may as discussed above include a plurality of

"microscopic" partial volumes presenting different degrees of isotropic diffusion (e.g. Diso as defined above) or different degrees and/or orientations of anisotropic diffusion (e.g. A _D as defined above). As is well-known in the art the spatial resolution of an NMR spectrometer or MRI device is limited by inter alia the strength of the magnetic field, the magnitude of the gradient pulse sequence applied to the sample and slew rate. The data analysis disclosed in the below enables estimation of diffusion properties of the microscopic partial volumes within the portion, i.e. beyond the traditional resolution limitations of the measurements. To identify the echo signal component corresponding to the portion, the measurement signals from the sample may be subjected to a Fast Fourier Transform as is well-known in the art, thereby transforming the spectral components of each echo signal from the sample into a plurality of spatial regions of the sample.

According to the inventive method, at least a subset of the data representing the echo signal measurements is acquired with a set of magnetic gradient pulse sequences causing anisotropic diffusion weighting, wherein the diffusion encoding tensor for each gradient pulse sequence of the set of magnetic gradient pulse sequences has 3 non-zero eigenvalues, at least one of the 3 eigenvalues being different from the other eigenvalues. With reference to Eq. (9) the 3 non-zero eigenvalues of each pulse sequence may be denoted b _xx ^PAS , b _yy ^PAS and b _zz ^PAS . Analogously, each pulse sequence of said set is definable by a total diffusion weighting b, an anisotropy Δ*,, and an asymmetry ¾ following the definitions in Eq. (9). These parameters form invariant parameters for each pulse sequence. The parameter b represents the diffusion weighting magnitude for the gradient sequence.

An embodiment will now be described wherein echo signals are acquired using magnetic gradient pulse sequences generated such that the asymmetry parameter η _ύ equals 0. This implies that at least two of the eigenvalues of the corresponding diffusion weighting encoding tensor are equal (i.e. b _xx ^PAS = b _yy ^PAS ). Data representing the echo signal measurements may be collected as described above. The data may include a number of distinct subsets of data, each acquired using a different set of pulse sequences. For example, the data may include a first subset of data, a second subset of data and a third subset of data. The first, second and third subset of data may represent echo signal measurements acquired using a first (denoted and a third set (denoted of magnetic gradient pulse sequences, respectively. Although reference is made to these three subsets of data and sets of pulse sequences, the method is applicable also to both greater and smaller numbers of such subsets of data and sets of pulse sequences.

The first set of pulse sequences may each be generated to have an equal degree of anisotropy A _{b 1} but different maximum amplitudes b (defined in accordance with Eq. (9) and (16)). Analogously, the second set of pulse sequences may each be generated to have an equal degree of anisotropy A _{b 2} but different maximum amplitudes b. Likewise, the third set of pulse sequences may each be generated to have an equal degree of anisotropy Δ _{Α 3} but different maximum amplitudes b. For example, the pulse sequences may be generated using the approach described in connection with Table 1 . With reference to Fig. 2, it may be understood that the first, the second and the third subsets of data each represent respective echo attenuation curves acquired along respective lines of constant anisotropy, namely A _{b 1} , A _{A 2} and

Δ _{Α 3} . Using the thusly collected data, D, _so for the diffusion tensor D and/or A _D for the diffusion tensor D may be calculated.

For a sample including an ensemble of randomly oriented microscopic partial volumes of generally similar degrees of anisotropic and isotropic diffusion, (estimates of) D _\so and A _D for a randomly oriented diffusion tensor D characterizing diffusion in each partial volume, may be calculated by fitting the data using Eq. (25). If for example the second subset of data is acquired with A ² _ft = 0 the isotropic value D, _so for the diffusion tensor D may be directly calculated using Eq. 26 and the second subset of data.

In case there is some preferential alignment of the orientations of diffusion in the partial volumes in the sample, powder averaging may be applied, wherein each pulse sequence of each one of the first, second and third sets of magnetic gradient pulse sequences may be applied a plurality of times to the sample, with different orientations of the gradient pulse with respect to a fixed laboratory frame. The first, second and third subsets of data may thereafter be formed by averaging echo signal measurements acquired for the different orientations. Equations (25) and (26) may thereafter be used in a same manner as for a sample including randomly oriented domains.

For a sample including a collection of microscopic partial volumes exhibiting different degrees of anisotropic and/or isotropic diffusion (i.e. two- or multi-component material), the data (with or without powder averaging depending on whether there is a preferential alignment of orientations or not) may be used to estimate D, _so and/or the Δ for each component. In particular, the data analysis approach described in connection with Eqs. (32-34) may be used. The kernel matrix may be calculated using pairs of the model parameters {b,A _b ) and the model parameters (Οί ₅₀ ,Δο), where the values of {b,A _b ) correspond to the values for each gradient pulse sequence used during the measurements, and the values of (Οί ₅₀ ,Δο) are selected to cover the region of interest for the sample, for example based on a priori knowledge of the range of possible values of (Οί ₅₀ ,Δο) for the particular sample. If used in connection with diffusion MRI, D, _so and/or the A _D may be used to generate contrast for a voxel representing the portion of the sample. Similar

calculations may be performed to generate contrast for voxels representing other portions of the sample.

Although in the above, calculation of D, _so and A _D is described employing e.g. Eqs. (25-26) it should be noted that the parameters may be calculated or estimated also in other ways. The measurement data could also be fitted using a different model function relating diffusion weighted signals to the relevant diffusion metrics, e.g. D, _so or A _D . As one example of an alternative to Eq. (25), an expansion of Eq. (25) in terms of moments μ ₂ , μ ₃ , etc. of the distribution of diffusivities could be used instead. Another example would be to approximate the distribution of diffusivities by the gamma distribution function.

In the above, various embodiments have been disclosed wherein diffusion parameters are calculated based on echo signal measurements acquired using diffusion encoding tensors being axially symmetric. In the following, embodiments of the inventive method will be disclosed which do not require axially symmetric diffusion encoding tensors. One objective is to provide methods that can quantify differences in microstructure, for example in terms of microscopic diffusion anisotropy on the basis of diffusion

properties of "microscopic" partial volumes of a portion of the sample. To facilitate understanding of the following embodiments, a discussion of some theoretical concepts will now be provided. Theory

Consider a portion of a sample (on which diffusion NMR/MRI is to be performed) which includes a collection of partial volumes (e.g. "microscopic" partial volumes), where in each partial volume the diffusion is Gaussian and described by the diffusion tensor D. Diffusion properties in these

microenvironments within the portion may be modeled with a Gaussian distribution over tensors. The tensor D may thus be referred to as a stochastic variable with expectation (D), where ( ^■ ) represents integration over the distribution in the portion. The covariance of D may then be given by a 4th- order tensor S defined using a standard definition of covariance according to where D = D ® D is a fourth-order tensor obtained from the outer product of D with itself. The diffusion encoded MR-signal E from a portion including multiple such microscopic partial volumes, each having Gaussian diffusion, may be estimated by

E(b) = (exp ( - < b, D >)), (36) where < ^■ , ^■ > is the inner product. To facilitate understanding it may be noted that equation (36) is based on equation (6) however differs in that E denotes the normalized echo signal intensity (i.e. S / So) and includes the signal contribution from each of the microscopic environments of the portion. Furthermore, to simplify notation in the following the inner product will be used instead of the generalized scalar product used in Eq. (7). From a cumulant expansion of equation (36) follows that, where B = b ⁰² and S is the covariance of tensors within the portion (which in the case of diffusion MRI is represented by a voxel). To facilitate

understanding, a detailed derivation of the cumulant expansion is provided below.

The approximation in Eq. 37 is a cumulant expansion where e{b) = log E(b) is expanded around b = 0 (where b corresponds to the total diffusion weighting in analo with Eq. 9) according to where

For b = 0, these functions evaluates to where D=D® ² . Hence f"(0)=<N,(D® ² )-(D)® ² >=<N,S>. In case of the diffusion tensor D, two common invariant representations are the mean diffusivity (MD) and the fractional anisotropy (FA). MD may be calculated as a projection onto an isotropic base Ε. _∞ , according to

MD(D)= < D,E _1S0 > (44) where E _iso = ± \ , i.e., a third of the identity tensor. This is equivalent to the parameter D, _so defined in connection with Equation (1 ). Analogously the fourth order covariance tensor S may be projected onto an isotropic base to obtain a rotationally invariant parameter. The isotropic 4th-order tensor S, however, has two isotropic components, which in analogy with the field of mechanics may be interpreted as bulk and shear modulus of the 4th-order stress tensor (see e.g. reference [28]). These two bases may be defined by

and

shear (46) where I is the identity tensor. Note that Ei = Ebuik and E2 = E _S hear are orthogonal, i.e.,

< Ε,, Ε _; > = δ„. (47) where 6 _/y is unity if / = j, and zero otherwise.

Similarly to estimating the mean diffusivity MD as in Eq. (44), the 4th- order covariance tensor may be projected onto its two isotropic basis elements. Projecting onto Ebuik yields the variance in mean diffusivities VMD (which also may be denoted Sbuik) according to,

< S,E _bulk >=< MD(D) ² ^-MD((D)) ² = V _MO . (48)

This follows from the following equations Eqs. (49-53):

< S, E _bulk >=<(D-)-(D) ^e2 , E- >= (49)

< (D-), E- > - <(D) ⁰² , E- >= (50)

=(MD ² )-(MD) ² = V _N

Projecting S onto E _S hear yields another invariant parameter related to the variance of tensor eigenvalues,

<§Ε _5ΐΚΜ >=(\/ _λ (ϋ))-Κ((ϋ)) = Δ\/ _λ ,

(which also may be denoted S _S hear) where VK{ ) yields the variance of diffusion tensor eigenvalues (λ, , i = 1 , 2, 3),

V,(D)=<D ⁸² , E _shear >=

This follows from following equations Eqs. (57-64): Considering the projection not of S but of D on E _shear ,

<D ² ,E. >-<D,E > = (59) where the following relationship was utilized:

< D ² , E. >= ! o Tr( VDD) / = - o < D,' D >=< D ⁰² ' Io -| > (61)

Based on the above, the projection of S follows by:

< _¾

= (< D ^e \ E _shear >)- < (D) ⁰² , E _shear >= (63)

As realized by the inventors, V _M D can be interpreted as the bulk variation of diffusion tensors (i.e. variation in size) and Δ¼ as the shear of them (i.e. variation between directions).

In conventional diffusion tensor imaging (DTI), an often considered invariant parameter is the fractional anisotropy (FA). It is defined by the normalized variance of the eigenvalues of the average tensor (see e.g.

reference [29])

By utilizing Eqs. (44-48) and (54-56), the FA may be expressed as projections of the conventional average diffusion tensor, raised to a higher order by taking the outer product with itself, on the bulk and shear bases according to

Hence, by replacing the variance of the average diffusion tensor, V _\ ((D)), with the average variance of the eigenvalues of the microscopic compartment tensors (i.e. of the diffusion tensors of the partial volumes of the portion of the sample), (\A(D)), the microscopic FA (^FA) may be obtained according to where^D ^®2 ^ = S+(D) ^®2 . If all microscopic tensors share the same set of eigenvalues they share the same value of \A(D), and in that case ^FA will yield the exact FA of the microscopic tensors.

By calculating the FA using Equation (67), the parameter becomes insensitive to orientation dispersion since the outer product in the ^FA calculation acts on the local diffusion tensor, not the globally averaged tensor. This also allows for simple implementation and estimation.

Voigt notation of tensors

For implementation purposes, the tensors D and S may advantageously be represented in Voigt notation, which allows D to be represented as a column vector of size 6x1

The fourth order 4th-order tensor S can now be represented by a 6 χ 6 variance-covariance matrix, defined in terms of d according to since D ⁰² = dd ^T . In full, S is given by (70)

A tensor of rank four, such as B, may either be represented by a 6 χ 6 matrix (analogous to Equation (70)), B = b® ² = b _v b _v ^T , or in Voigt notation by a 21 x 1 column vector b according to:

Using Voigt notation the bulk and shear bases may be represented by:

and

As indicated above, these two isotropic bulk and shear bases are by design orthogonal, and as can be noted from Eqs. (45) and (46), adding these basis functions gives the simple structure of the diagonal identity matrix .

Inner and outer products

One advantage of using matrix and vector representations of the tensors in implementations is that inner and outer products become straightforward to implement in software. The outer product of a tensor (e.g. D) may be calculated according to

D ⁸² = dd ^T ( ⁷⁴ )

One may thus define, for example, N = n® ² or in Voigt notation n = n® ² .

Inner products are represented by < ^■ , ^■ >, and may be defined as an element-wise multiplication followed by summation according to

< D, N >= d ^T n ( ⁷⁵ ) or

< S, B> s b (76) The inner product of two matrices may also be defined according to

Description of embodiments

Taking the above into account, according to a preferred embodiment, a plurality of diffusion weighted echo attenuation measurements are performed on the sample using a plurality of different diffusion encoding magnetic gradient pulse sequences, wherein each magnetic gradient pulse sequence Gj is generated such that a diffusion encoding tensor bj for the magnetic gradient pulse sequence Gj. The plurality of pulse sequences may include a combination of pulse sequences with diffusion encoding tensors having 1 to 3, and preferably 2 to 3, non-zero eigenvalues. The above discussion concerning implementation of the method on an NMR spectroscope or MRI device applies also to the present embodiment.

In a general case the q-vector may be built up by a time-dependent gradient to traverse an arbitrary path in q-space. The rank (i.e. the number of non-zero eigenvalues) of the diffusion encoding tensor depends on the path, and becomes 1 in the case of sPFG, 2 for dPFG when the first and the second gradient pulse are applied along orthogonal directions, and 3 in the isotropic encoding case such as the triple-PFG [4] or q-MAS [30]. For example, a planar diffusion encoding tensor, i.e. an encoding that is rotationally symmetric in the plane can be achieved by a set of time varying gradients that produce a planar q-space trajectory. Constant angular b-value encoding can be ensured by varying the speed of the traversal in q-space, by using slower speed at low q-values, since the b-value is a function of both time and q-value. At a low q, a long diffusion time can build up the same encoding power (b-value), as a higher q-value with a shorter diffusion time.

A set of data {Ei, ... E _m } representing the echo attenuation

measurements acquired using the pulse sequences Gj may be collected and arranged in a column vector representation, b and (D) may (using Voigt notation) be arranged in column vector representations of size 6 x 1 , denoted by b _v and (d), and similarly, the fourth-order tensors B and S may be arranged as column vector representations of size 21 x 1 , denoted by b and s. The elements of each bj may be obtained using the definitions in Eqs. (5) and (8). The inner products in Eq. (37) can now be expressed by simple matrix operations according to

< b,(D) >= b ^T ,(d) (78) and

(79)

< B,S>= b s

Since Eq. (37) is a linear model, (d) and s may be estimated using pseudoinversion to solve the following equation system logE,

logE,

1

Equation (80) forms a linear system of equations based on the cumulant expansion of Equation (36). The data {Ei, ... Em} forms constants of the linear system. The diffusion encoding tensor representations b _v i... _m of the pulse sequences Gj and the fourth-order tensor representations bi... _m thereof form parameters of the linear system.

In total, the model has 1 + 6 + 21 = 28 free parameters in EQ, (d), and s. To enable estimation of a solution to Equation (80) using pseudoinversion, the data should be acquired with measurement tensors of varying shapes, such that the correlation between bulk component (b ¾ _u]k ) and shear component b^e _shear ) become less than unity. This explains why separation of the two isotropic components of S would not be possible with conventional sPFG. Assuming that the bulk and shear components of the encoding tensors are not fully correlated, the pseudoinversion may be performed if the number of measurements exceeds 28. However, fewer measurements can be used if only projections of <d>, and s are sought for.

In the above the inventive concept has mainly been described with reference to a limited number of examples. However, as is readily

appreciated by a person skilled in the art, other examples than the ones disclosed above are equally possible within the scope of the inventive concept, as defined by the appended claims. List of references

In the above disclosure, one or more numbers between a pair of brackets "[ ]" refer to a correspondingly numbered reference document in the following list of references:

[1 ] P.T. Callaghan, 0. Soderman, Examination of the lamellar phase of

Aerosol OT/water using pulsed field gradient nuclear magnetic resonance, J. Phys. Chem., 87 (1983) 1737-1744.

[2] M.E. Moseley, J. Kucharczyk, H.S. Asgari, D. Norman, Anisotropy in diffusion-weighted MRI, Magn. Reson. Med., 19 (1991 ) 321 -326.

[3] E.O. Stejskal, J.E. Tanner, Spin diffusion measurements: Spin echoes in the presence of a time-dependent field gradient, J. Chem. Phys., 42 (1965) 288-292.

[4] S. Eriksson, S. Lasic, D. Topgaard, Isotropic diffusion weighting in PGSE NMR by magic-angle spinning of the q-vector, J. Magn. Reson., 226 (2013) 13-18.

[5] L. Frydman, G.C. Chingas, Y.K. Lee, P.J. Grandinetti, M.A.

Eastman, G.A. Barrall, A. Pines, Variable-angle correlation spectroscopy in solid-state nuclear magnetic resonance, J. Chem. Phys., 97 (1992) 4800- 4808.

[6] A. Bax, N.M. Szeverenyi, G.E. Maciel, Chemical shift anisotropy in powdered solids studied by 2D FT NMR with flipping of the spinning axis, J. Magn. Reson., 55 (1983) 494-497.

[7] Y.-Q. Song, L. Venkataramanan, M.D. Hurlimann, M. Flaum, P. Frulla, C. Straley, T1 -T2 correlation spectra obtained using a fast two- dimensional Laplace inversion, J. Magn. Reson., 154 (2002) 261 -268.

[8] M. Lustig, D. Donoho, J.M. Pauly, Sparse MRI: The application of compressed sensing for rapid MR imaging, Magn. Reson. Med., 58 (2007) 1 182-1 195.

[9] M. Mobli, M.W. Maciejewski, A.D. Schuyler, A.S. Stern, J.C. Hoch, Sparse sampling methods in multidimensional NMR, Phys. Chem. Chem. Phys., 14 (2012) 10835-10843.

[10] D. Topgaard, Isotropic diffusion weighting in PGSE NMR:

Numerical optimization of the q-MAS PGSE sequence, Microporous

Mesoporous Mater., 178 (2013) 60-63.

[1 1 ] P.J. Basser, J. Mattiello, D. Le Bihan, MR diffusion tensor spectroscopy and imaging, Biophys. J., 66 (1994) 259-267. [12] K. Schmidt-Rohr, H.W. Spiess, Multidimensional solid-state NMR and polymers, Academic Press, San Diego, (1994).

[13] M. Duer, Introduction to solid-state NMR spectroscopy, Blackwell Publishing Ltd, Oxford, (2004).

[14] W.S. Price, NMR studies of translational motion, Cambridge

University Press, Cambridge, (2009).

[15] P.T. Callaghan, Translational dynamics & magnetic resonance, Oxford University Press, Oxford, (201 1 ).

[16] G.B. Arfken, H.-J. Weber, Mathematical methods for physicists, 4th ed., Academic Press, San Diego, (1995).

[17] V.J. Wedeen, D.L. Rosene, R. Wang, G. Dai, F. Mortazavi, P. Hagmann, J.H. Kaas, W.-Y. l. Tseng, The geometric structure of the brain fiber pathways, Science, 335 (2012) 1628-1634.

[18] K.P. Whittal, Analysis of large one-dimensional and two- dimensional relaxation data sets, J. Magn. Reson. A, 1 10 (1994) 214-218.

[19] D. Bernin, D. Topgaard, NMR diffusion and relaxation correlation methods: New insights in heterogeneous materials, Curr. Opin. Colloid Interface Sci., 18 (2013) 166-172.

[20] J. Mitchell, T.C. Chandrasekera, L.F. Gladden, Numerical estimation of relaxation and diffusion distributions in two dimensions, Prog. Nucl. Magn. Reson. Spectrosc, 62 (2012) 34-50.

[21 ] B. Jonsson, B. Lindman, K. Holmberg, B. Kronberg, Surfactants and polymers in aqueous solution, John Wiley & Sons Ltd, (1998).

[22] M. Bak, N.C. Nielsen, REPULSION, A novel approach to efficient powder averaging in solid-state NMR, J. Magn. Reson., 125 (1997) 132-139.

[23] D.K. Jones, M.A. Horsfield, A. Simmons, Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging, Magn. Reson. Med., 42 (1999) 515-525.

[24] Le Bihan, D., J.B.H.: Diffusion MRI at 25: exploring brain tissue structure and function (2012)

[25] Assaf, Y., P.O.: Diffusion MRI at 25: exploring brain tissue structure and function (2008)

[26] Jones, D. K., K.T.R.T.R.: White matter integrity, fiber count, and other fallacies: the dos and don't's of diffusion mri (2013)

[27] Vos, S. B., J.D.K.J.B.V.M.A.L.A.: The influence of complex white matter architecture on the mean diffusivity in diffusion tensor mri of the human brain (2012) [28] Moakher, M.: Fourth-order cartesian tensors: old and new facts, notions and applications. The Quart. Jour, of Mechanics and Applied Math. 61 (2) (2008) 181-203

[29] Basser, P.J., Pierpaoli, C: Microstructural and physiological features of tissues elucidated by quantitative-diffusion-tensor mri. Journal of Magnetic Resonance, Series B 1 1 1 (3) (1996) 209-219

[30] Valette, J., Giraudeau, C, Marchadour, C, Djemai, B., Geffroy, F., Ghaly, M.A., Le Bihan, D., Hantraye, P., Lebon, V., Lethimonnier, F.: A new sequence for single- shot diffusion-weighted nmr spectroscopy by the trace of the diffusion tensor. MRM 68(6) (2012) 1705-1712