Pulse sequence method for MRI

Abstract

The present invention describes a method for magnetic resonance (MR) and/or MR imaging, comprising acquisition of signals and MR images originating from a RF and gradient sequence causing isotropic diffusion weighting of signal attenuation, wherein the isotropic diffusion weighting is achieved by one time-dependent dephasing vector q(t) having an orientation, wherein the isotropic diffusion weighting is proportional to the trace of a diffusion tensor D, and wherein the orientation of the time-dependent dephasing vector q(t) is either varied discretely in more than three directions in total, or changed continuously, or changed in a combination of discretely and continuously during the gradient pulse sequence, 0≦t≦echo time, where t represents the time. The method may be performed during a single shot (single MR excitation).

Claims

1. A method for magnetic resonance (MR) and/or MR imaging, the method comprising: performing a diffusion encoding scheme comprising a radio frequency (RF) and gradient sequence causing an isotropic diffusion weighting; and acquiring signals resulting from the diffusion encoding scheme, wherein the isotropic diffusion weighting is proportional to a trace of a diffusion tensor D, wherein the isotropic diffusion weighting is achieved by one time-dependent dephasing vector q(t) having an orientation, wherein the orientation of the time-dependent dephasing vector q(t) is varied discretely in more than three directions in total, changed continuously, or changed in a combination of discretely and continuously during a gradient pulse sequence, 0≦t≦t.sub.E, where t represents time and t.sub.E represents echo time, wherein the isotropic diffusion weighting is invariant under rotation of the diffusion tensor D and wherein an equation ${\int }_0^{t_E}{F\left(t\right)}^2{\hat{q}}^T\left(t\right).Math.D.Math.\hat{q}\left(t\right)dt=t_d\stackrel{\_}{D}$ is fulfilled, F(t) being a time-dependent normalized magnitude of the dephasing vector, {circumflex over (q)}(t) being a time-dependent unit direction vector, t.sub.d being an effective diffusion time given by $t_d={\int }_0^{t_E}{F\left(t\right)}^{2}dt$ and D being an isotropic mean diffusivity.

2. The method according to claim 1, wherein a time-dependent normalized magnitude F(t) of the dephasing vector is |F(t)|≦1 during an echo time t.sub.E from t=0 to t=t.sub.E and wherein the orientation of the dephasing vector at time 0 is identical to the orientation of the dephasing vector at time t.sub.E.

3. The method according to claim 1, wherein orientation of the time-dependent dephasing vector q(t) is changed with discrete steps in azimuth angle ψ, providing q(t) vector steps through at least four orientations with unique values of e.sup.iψ, such that ψ modulus 2π are equally spaced values.

4. The method according to claim 1, wherein the isotropic diffusion weighting is achieved by a continuous sweep of the time-dependent dephasing vector q(t) where an azimuth angle ψ(t) and a magnitude thereof is a continuous function of time so that the time-dependent dephasing vector q(t) spans an entire range of orientations parallel to a right circular conical surface and so that the orientation of the time-dependent dephasing vector q(t) at time 0 is identical to the orientation at time t.sub.E.

5. The method according to claim 4, wherein a particular inclination ζ is chosen to be a constant, time-independent value.

6. The method according to claim 5, wherein a particular inclination ζ is chosen so that ζ=ζ.sub.m=a cos(1/√{square root over (3)}).

7. The method according to claim 3, wherein the time-dependent normalized magnitude of the dephasing vector, F(t), is chosen as a harmonic function of time.

8. The method according to claim 1, wherein the method is performed in a single shot.

9. The method according to claim 1, wherein the orientation of the time-dependent dephasing vector q(t) is changed continuously.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIGS. 1A-D to 6A-D show examples of different gradient modulation schemes for isotropic diffusion weighting according to the present invention. Insets A depict components of the normalized dephasing vector, q.sub.x/|q| (dashed line), q.sub.y/|q| (dotted line) and q.sub.z/|q| (dash dotted line) and the normalized magnitude of the dephasing vector, F(t) (solid line). Insets B depict components of the normalized effective gradient vector, g.sub.x/|g| (dashed line), g.sub.y/|g| (dotted line) and g.sub.x/|g| (dash dotted line). Insets C depict time dependence of the azimuth angle. Insets D depict the evolution of the anisotropic diffusion weighting terms (16) as a function of time; the first term in Eq. (16) is shown as a dotted line, the second term is shown as a dashed dotted line, the third term as a solid line and the fourth term is shown as a dashed line.

(2) FIG. 7A-C show schematic representations of signal decays vs. b for isotropic (dashed line) and non-isotropic (solid line) diffusion weighting for different types of materials. The inset A depicts signal attenuation curves in case of anisotropic materials with 1D or 2D curvilinear diffusion. The attenuation curves are multi-exponential for non-isotropic diffusion weighting, while they are mono-exponential for isotropic diffusion weighting. The deviation between the attenuation curves for isotropic and non-isotropic diffusion weighting provides a measure of anisotropy. The inset B depicts an example of isotropic material with several apparent diffusion contributions resulting in identical and multi-exponential signal attenuation curves for isotropic and non-isotropic diffusion weighting. The inset C depicts an example of material with a mixture of isotropic and anisotropic components resulting in multi-exponential signal decays for both isotropic and non-isotropic diffusion weighting, while the deviation between the attenuation curves for isotropic and non-isotropic diffusion weighting provides a measure of anisotropy.

(3) FIG. 8A-C show experimental results with analysis for different types of materials. Experimental results for isotropic (circles) and for non-isotropic (crosses) diffusion weighting are shown in all the insets. Experimental results and analysis are shown for a sample with free isotropic diffusion (inset A), for a sample with restricted isotropic diffusion (inset B) and for a sample with high degree of anisotropy (inset C).

(4) FIGS. 9A and 9B show a Monte-Carlo error analysis for the investigation of systematic deviations and precision as a function of the range of diffusion weighting b for estimating the degree of micro-anisotropy with the disclosed analytical method.

DESCRIPTION; BACKGROUND AND SOME SPECIFIC EMBODIMENTS OF THE PRESENT INVENTION

(5) Assuming that spin diffusion in a microscopically anisotropic system can locally be considered a Gaussian process and therefore fully described by the diffusion tensor D(r), the evolution of the complex transverse magnetization m(r,t) during a diffusion encoding experiment is given by the Bloch-Torrey equation. Note that the Bloch-Torrey equation applies for arbitrary diffusion encoding schemes, e.g. pulse gradient spin-echo (PGSE), pulse gradient stimulated echo (PGSTE) and other modulated gradient spin-echo (MGSE) schemes. Assuming uniform spin density and neglecting relaxation, the magnetization evolution is given by

(6) $\begin{array}{cc}\frac{\partial m\left(r,t\right)}{\partial t}=-i\gamma g\left(t\right).Math.\mathrm{rm}\left(r,t\right)+▽.Math.\left[D\left(r\right).Math.▽m\left(r,t\right)\right],& \left(1\right)\end{array}$
where γ is the gyromagnetic ratio and g(t) is the time dependent effective magnetic field gradient vector. The NMR signal is proportional to the macroscopic transverse magnetization

(7) $\begin{array}{cc}M\left(t\right)=\frac{1}{V}{\int }_V^{}m\left(r,t\right)dr.& \left(2\right)\end{array}$

(8) If during the experiment each spin is confined to a domain characterized by a unique diffusion tensor D, the macroscopic magnetization is a superposition of contributions from all the domains with different D. Evolution of each macroscopic magnetization contribution can thus be obtained by solving Eqs. (1, 2) with a constant and uniform D. The signal magnitude contribution at the echo time t.sub.E is given by

(9) $\begin{array}{cc}\begin{array}{c}I\left(t_E\right)=I_0\exp \left(-{\int }_0^{t_E}{q}^T\left(t\right).Math.D.Math.q\left(t\right)dt\right)\\ =I_0\exp \left(-q^2{\int }_0^{t_E}{F\left(t\right)}^2{\hat{q}}^T\left(t\right).Math.D.Math.\hat{q}\left(t\right)dt\right),\end{array}& \left(3\right)\end{array}$
where I.sub.0 is the signal without diffusion encoding, g=0, and q(t) is the time-dependent dephasing vector

(10) $\begin{array}{cc}q\left(t\right)=\gamma {\int }_0^{t}g\left(t^{\prime }\right)d{t}^{\prime }=\mathrm{qF}\left(t\right)\hat{q}\left(t\right),& \left(4\right)\end{array}$
defined for the interval 0<t<t.sub.E. The dephasing vector in Eqs. (3) and (4) is expressed in terms of its maximum magnitude q, the time-dependent normalized magnitude |F(t)|≦1 and a time-dependent unit direction vector {circumflex over (q)}(t). Note that in spin-echo experiments, the effective gradient g(t) comprises the effect of gradient magnitude reversal after each odd 180° radio frequency (RF) pulse in the sequence. Eq. (3) assumes that the condition for the echo formation q(t.sub.E)=0 is fulfilled, which implies F(t.sub.E)=0. In general there might be several echoes during an NMR pulse sequence.

(11) The echo magnitude (3) can be rewritten in terms of the diffusion weighting matrix,

(12) $\begin{array}{cc}b=q^2{\int }_0^{t_E}{F\left(t\right)}^2\hat{q}\left(t\right).Math.{\hat{q}}^T\left(t\right)dt,\mathrm{as}& \left(5\right)\\ I\left(t_E\right)=I_0\exp \left(-b:D\right)=I_0\exp \left(-\underset{i}{.Math.}\underset{j}{.Math.}{b}_{\mathrm{ij}}{D}_{\mathrm{ij}}\right).& \left(6\right)\end{array}$
Integral of the time-dependent waveform F(t).sup.2 defines the effective diffusion time, t.sub.d, for an arbitrary diffusion encoding scheme in a spin-echo experiment

(13) $\begin{array}{cc}{t}_d={\int }_0^{t_E}{F\left(t\right)}^{2}dt.& \left(7\right)\end{array}$

(14) In the following we will demonstrate that even for a single echo sequence, gradient modulations g(t) can be designed to yield isotropic diffusion weighting, invariant under rotation of D, i.e. the echo attenuation is proportional to the isotropic mean diffusivity,
D=tr(D)/3=(D.sub.xx+D.sub.yy+D.sub.zz)/3.  (8)

(15) In view of what is disclosed above, according to one specific embodiment of the present invention, the isotropic diffusion weighting is invariant under rotation of the diffusion tensor D.

(16) According to the present invention, one is looking for such forms of dephasing vectors F(t){circumflex over (q)}(t), for which

(17) $\begin{array}{cc}{\int }_0^{t_E}{F\left(t\right)}^2{\hat{q}}^T\left(t\right).Math.D.Math.\hat{q}\left(t\right)dt=t_d\stackrel{\_}{D}& \left(9\right)\end{array}$
is invariant under rotation of D. If diffusion tenor D is expressed as a sum of its isotropic contribution, DI, where I is the identity matrix, and the anisotropic contribution, i.e. the deviatoric tensor D.sup.A, so that D=DI+D.sup.A, the isotropic diffusion weighing is achieved when the condition

(18) $\begin{array}{cc}{\int }_0^{t_E}{F\left(t\right)}^2{\hat{q}}^T\left(t\right).Math.D^A.Math.\hat{q}\left(t\right)dt=0& \left(10\right)\end{array}$
is fulfilled.

(19) In spherical coordinates, the unit vector {circumflex over (q)}(t) is expressed by the inclination ζ and azimuth angle ψ as
{circumflex over (q)}.sup.T(t)={{circumflex over (q)}.sub.x(t),{circumflex over (q)}.sub.y(t),{circumflex over (q)}.sub.z(t)}={ sin ζ(t)cos ψ(t),sin ζ(t)sin ψ(t),cos ζ(t)}.  (11)
The symmetry of the diffusion tensor, D=D.sup.T, gives
{circumflex over (q)}.sup.T.Math.D.Math.{circumflex over (q)}={circumflex over (q)}.sub.x.sup.2D.sub.xx+{circumflex over (q)}.sub.y.sup.2D.sub.yy+{circumflex over (q)}.sub.z.sup.2D.sub.zz+2{circumflex over (q)}.sub.x{circumflex over (q)}.sub.yD.sub.xy+2{circumflex over (q)}.sub.x{circumflex over (q)}.sub.zD.sub.xz+2{circumflex over (q)}.sub.y{circumflex over (q)}.sub.zD.sub.yz  (12)
or expressed in spherical coordinates as
{circumflex over (q)}.sup.T.Math.D.Math.{circumflex over (q)}=sin.sup.2 ζ cos.sup.2 ψD.sub.xx+sin ζ.sup.2 sin ψ.sup.2D.sub.yy+cos.sup.2 ζD.sub.zz+2 sin ζ cos ψ sin ζ sin ψD.sub.xy+2 sin ζ cos ψ cos ζD.sub.xz+2 sin ζ sin ψ cos ζD.sub.yz.  (13)
Equation (13) can be rearranged to

(20) $\begin{array}{cc}\hat{q}.Math.D.Math.\hat{q}=\stackrel{\_}{D}+\frac{3{\cos }^2\zeta -1}{2}\left(D_{\mathrm{zz}}-\stackrel{\_}{D}\right)+{\sin }^2\zeta \left[\frac{D_{\mathrm{xx}}-D_{\mathrm{yy}}}{2}\cos \left(2\psi \right)+D_{\mathrm{xy}}\sin \left(2\psi \right)\right]+\sin \left(2\zeta \right)\left(D_{\mathrm{xz}}\cos \psi +D_{\mathrm{yz}}\sin \psi \right).& \left(14\right)\end{array}$
The first term in Eq. (14) is the mean diffusivity, while the remaining terms are time-dependent through the angles at) and ζ(t) and ψ(t) which define the direction of the dephasing vector (4). Furthermore, the second term in Eq. (14) is independent of ψ, while the third and the forth terms are harmonic functions of ψ and 2ψ, respectively (compare with Eq. (4) in [13]). To obtain isotropic diffusion weighting, expressed by condition in Eq. (9), the corresponding integrals of the second, third and fourth terms in Eq. (14) must vanish. The condition for the second term of Eq. (14) to vanish upon integration leads to one possible solution for the angle ζ(t), i.e. the time-independent “magic angle”
ζ.sub.m=a cos(1/√{square root over (3)}).  (15)

(21) By taking into account constant ζ.sub.m, the condition for the third and the fourth term in Eq. (14) to vanish upon integration is given by

(22) 0$\begin{array}{cc}{\int }_0^{t_E}{F\left(t\right)}^2\cos \left[2\psi \left(t\right)\right]dt=0{\int }_0^{t_E}{F\left(t\right)}^2\sin \left[2\psi \left(t\right)\right]dt=0{\int }_0^{t_E}{F\left(t\right)}^2\cos \left[\psi \left(t\right)\right]dt=0{\int }_0^{t_E}{F\left(t\right)}^2\sin \left[\psi \left(t\right)\right]dt=0.& \left(16\right)\end{array}$
Conditions (16) can be rewritten in a more compact complex form as

(23) $\begin{array}{cc}{\int }_0^{t_E}{F\left(t\right)}^2\exp \left[ik\psi \left(t\right)\right]dt=0,& \left(17\right)\end{array}$
which must be satisfied for k=1, 2. By introducing the rate {dot over (τ)}(t)=F(t).sup.2, the integral (17) can be expressed with the new variable τ as

(24) $\begin{array}{cc}{\int }_0^{t_d}\exp \left[ik\psi \left(\tau \right)\right]d\tau =0.& \left(18\right)\end{array}$
Note that the upper integration boundary moved from t.sub.E to t.sub.d. The condition (18) is satisfied when the period of the exponential is t.sub.d, thus a solution for the azimuth angle is

(25) $\begin{array}{cc}\psi \left(\tau \right)=\psi \left(0\right)+\frac{2\pi }{t_d}n\tau ,& \left(19\right)\end{array}$
for any integer n other than 0. The time dependence of the azimuth angle is finally given by

(26) $\begin{array}{cc}\psi \left(t\right)=\psi \left(0\right)+\frac{2\pi n}{t_d}{\int }_0^t{F\left(t^{\prime }\right)}^{2}d{t}^{\prime }.& \left(20\right)\end{array}$
The isotropic diffusion weighting scheme is thus determined by the dephasing vector q(t) with a normalized magnitude F(t) and a continuous orientation sweep through the angles ζ.sub.m (15) and ψ(t) (20). Note that since the isotropic weighting is invariant upon rotation of D, orientation of the vector q(t) and thus also orientation of the effective gradient g(t) can be arbitrarily offset relative to the laboratory frame in order to best suit the particular experimental conditions.

(27) As understood from above, according to yet another specific embodiment, the isotropic diffusion weighting is achieved by a continuous sweep of the time-dependent dephasing vector q(t) where the azimuth angle ψ(t) and the magnitude thereof is a continuous function of time so that the time-dependent dephasing vector q(t) spans an entire range of orientations parallel to a right circular conical surface and so that the orientation of the time-dependent dephasing vector q(t) at time 0 is identical to the orientation at time t.sub.E. Furthermore, according to yet another embodiment, the inclination ζ is chosen to be a constant, time-independent value, i.e. the so called magic angle, such that ζ=ζ.sub.m=a cos(1/√{square root over (3)}). It should be noted that the method according to the present invention may also be performed so that ζ is chosen to be time-dependent, as far as condition (10) is fulfilled, however, this is not a preferred implementation.

(28) What is disclosed above implies that according to one specific embodiment of the present invention, the orientation of the dephasing vector, in the Cartesian coordinate system during the diffusion weighting sequence, spans the entire range of orientations parallel to the right circular conical surface with the aperture of the cone of 2*ζ.sub.m (double magic angle) and the orientation of the dephasing vector at time 0 is identical to the orientation of the dephasing vector at time t.sub.E, i.e. ψ(t.sub.E)−ψ(0)=2*π*n, where n is an integer (positive or negative, however not 0) and the absolute magnitude of the dephasing vector, q*F(t), is zero at time 0 and at time t.sub.E. Therefore, according to one specific embodiment, the time-dependent normalized magnitude F(t) of the dephasing vector is |F(t)|≦1 during an echo time t.sub.E from t=0 to t=t.sub.E and the orientation of the dephasing vector at time 0 is identical to the orientation of the dephasing vector at time t.sub.E.

(29) With reference to what is disclosed above it should be said that the concept of the magic angle is used in other types of methods in MR today. For instance in US2008116889 there is disclosed a method for magnetic resonance analysis or in fact MRI spectroscopy suggesting a magic angle technique. The turning around the magic angle as disclosed in US2008116889 is required to achieve higher spectroscopic resolution (reduce anisotropic susceptibility broadening). The method does no relate to diffusion weighting. According to the present invention the dephasing vector may be turned around the magic angle to achieve isotropic diffusion weighting, and is hence not related to turning the magnetic field or sample around the magic angle as described in US2008116889.

(30) According to the present invention, the isotropic weighting can also be achieved by q-modulations with discrete steps in azimuth angle ψ, providing q(t) vector steps through at least four orientations with unique values of e.sup.iψ, such that ψ modulus 2π are equally spaced values. Choice of the consecutive order and duration of the time intervals during which ψ is constant is arbitrary, provided that the magnitude F(t) is adjusted to meet the condition for isotropic weighing (10, 16).

Specific Implementations and Embodiments of the Present Invention

(31) The pulsed gradient spin-echo (PGSE) sequence with short pulses offers a simplest implementation of the isotropic weighting scheme according to the present invention. In PGSE, the short gradient pulses at times approximately 0 and t.sub.E cause the magnitude of the dephasing vector to instantaneously acquire its maximum value approximately at time 0 and vanish at time t.sub.E. The normalized magnitude is in this case given simply by F(t)=1 in the interval 0<t<t.sub.E and 0 otherwise, providing t.sub.d=t.sub.E. A simplest choice for the azimuth angle (20) is the one with n=1 and ψ(0)=0, thus

(32) $\begin{array}{cc}\psi \left(t\right)=\frac{2\pi t}{t_E}.& \left(21\right)\end{array}$
The dephasing vector is given by

(33) $\begin{array}{cc}{q}^T\left(t\right)=q\left\{\sqrt{\frac{2}{3}}\cos \left(\frac{2\pi t}{t_E}\right),\sqrt{\frac{2}{3}}\sin \left(\frac{2\pi t}{t_E}\right),\sqrt{\frac{1}{3}}\right\}.& \left(22\right)\end{array}$
The corresponding effective gradient, calculated from

(34) $\begin{array}{cc}{g}^T\left(t\right)=\frac{1}{\gamma }\frac{d}{dt}{q}^T\left(t\right)\mathrm{is}& \left(23\right)\\ g^T\left(t\right)=\frac{q}{\gamma }\left[\delta \left(0\right)-\delta \left(t_E\right)\right]\left\{\sqrt{\frac{2}{3}},0,\sqrt{\frac{1}{3}}\right\}+\sqrt{\frac{2}{3}}\frac{2\pi }{t_E}\frac{q}{\gamma }\left\{-\sin \left(\frac{2\pi t}{t_E}\right),\cos \left(\frac{2\pi t}{t_E}\right),0\right\}.& \left(24\right)\end{array}$
Here δ(t) is the Dirac delta function. Rotation around the y-axis by a tan(√2) yields

(35) $\begin{array}{cc}{g}^T\left(t\right)=\frac{q}{\gamma }\left[\delta \left(0\right)-\delta \left(t_E\right)\right]\left\{0,0,1\right\}+\sqrt{\frac{2}{3}}\frac{2\pi }{t_E}\frac{q}{\gamma }\left\{-\sqrt{\frac{1}{3}}\sin \left(\frac{2\pi t}{t_E}\right),\cos \left(\frac{2\pi t}{t_E}\right),-\sqrt{\frac{2}{3}}\sin \left(\frac{2\pi t}{t_E}\right)\right\}.& \left(25\right)\end{array}$
The effective gradient in Eqs. (24, 25) can conceptually be separated as the sum of two terms,
g(t)=g.sub.PGSE(t)+g.sub.iso(t).  (26)
The first term, g.sub.PGSE, represents the effective gradient in a regular PGSE two pulse sequence, while the second term, g.sub.iso, might be called the “iso-pulse” since it is the effective gradient modulation which can be added to achieve isotropic weighting.

(36) As may be seen from above, according to one specific embodiment of the present invention, the method is performed in a single shot, in which the latter should be understood to imply a single MR excitation.

(37) Analysis

(38) Below there will be disclosed a suggested analysis method which may be performed subsequent to the method disclosed above.

(39) Fractional anisotropy (FA) is a well-established measure of anisotropy in diffusion MRI. FA is expressed as an invariant of the diffusion tensor with eigenvalues λ.sub.1, λ.sub.2, and λ.sub.3,

(40) $\begin{array}{cc}FA=\sqrt{\frac{{\left({\lambda }_1-{\lambda }_2\right)}^2+{\left({\lambda }_1-{\lambda }_3\right)}^2+{\left({\lambda }_2-{\lambda }_3\right)}^2}{2\left({\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2\right)}}.& \left(27\right)\end{array}$
In typical diffusion MRI experiments, only a voxel average anisotropy can be detected. The sub-voxel microscopic anisotropy is often averaged out by a random distribution of main diffusion axis. Here we introduce a novel parameter for quantifying microscopic anisotropy and show how it can be determined by diffusion NMR.

(41) Information about the degree of micro-anisotropy can be obtained from comparison of the echo-attenuation curves, E(b)=I(b)/I.sub.0, with and without the isotropic weighting. Multi-exponential echo attenuation is commonly observed in diffusion experiments. The multi exponential attenuation might be due to isotropic diffusion contributions, e.g. restricted diffusion with non-Gaussian diffusion, as well as due to the presence of multiple anisotropic domains with varying orientation of main diffusion axis. The inverse Laplace transform of E(b) provides a distribution of apparent diffusion coefficients P(D), with possibly overlapping isotropic and anisotropic contributions. However, in isotropically weighed diffusion experiments, the deviation from mono-exponential attenuation is expected to originate mainly from isotropic contributions.

(42) In practice, the diffusion weighting b is often limited to a low-b regime, where only an initial deviation from mono-exponential attenuation may be observed. Such behaviour may be quantified in terms of the kurtosis coefficient K (Jensen, J. H., and Helpern, J. A. (2010). MRI quantification of non-Gaussian water diffusion by kurtosis analysis. NMR Biomed 23, 698-710.),

(43) 0$\begin{array}{cc}\ln E=-\stackrel{\_}{D}b+\frac{{\stackrel{\_}{D}}^{2}K}{6}{b}^2-\mathrm{.Math.}& \left(28\right)\end{array}$
The second term in Eq. (28) can be expressed by the second central moment of the distribution P(D).

(44) Provided that P(D) is normalized,

(45) $\begin{array}{cc}{\int }_0^{\infty }P\left(D\right)dD=1,& \left(29\right)\end{array}$
the normalized echo signal is given by the Laplace transform

(46) $\begin{array}{cc}E\left(b\right)={\int }_0^{\infty }P\left(D\right)e^{-\mathrm{bD}}dD.& \left(30\right)\end{array}$
The distribution P(D) is completely determined by the mean value

(47) $\begin{array}{cc}\stackrel{\_}{D}={\int }_0^{\infty }\mathrm{DP}\left(D\right)dD& \left(31\right)\end{array}$
and by the central moments

(48) $\begin{array}{cc}{\mu }_m={\int }_0^{\infty }{\left(D-\stackrel{\_}{D}\right)}^{m}P\left(D\right)dD.& \left(32\right)\end{array}$

(49) The second central moment gives the variance, μ.sub.2=σ.sup.2, while the third central moment, μ.sub.3, gives the skewness or asymmetry of the distribution P(D). On the other hand, the echo intensity can be expressed as a cumulant expansion (Frisken, B. (2001). Revisiting the method of cumulants for the analysis of dynamic light-scattering data. Appl Optics 40) by

(50) $\begin{array}{cc}\ln E=-\stackrel{\_}{D}b+\frac{{\mu }_2}{2}{b}^2-\mathrm{.Math.}& \left(33\right)\end{array}$
The first-order deviation from the mono-exponential decay is thus given by the variance of P(D).

(51) Assuming diffusion tensors with axial symmetry, i.e. λ.sub.1=D.sub.∥ and λ.sub.2=λ.sub.3=D.sub.⊥, and an isotropic distribution of orientation of the tensor's main diffusion axis, the echo-signal E(b) and the corresponding distribution P(D) can be written in a simple form. In case of the single PGSE experiment, using a single diffusion encoding direction, the distribution is given by

(52) $\begin{array}{cc}P\left(D\right)=\frac{1}{2\sqrt{\left(D-D_{\perp }\right)\left(D_{\parallel }-D_{\perp }\right)}},& \left(34\right)\end{array}$
with the mean and variance,

(53) $\begin{array}{cc}\stackrel{\_}{D}=\frac{D_{\parallel }+2D_{\perp }}{3}\mathrm{and}{\mu }_2=\frac{4}{45}{\left(D_{\parallel }-D_{\perp }\right)}^2.& \left(35\right)\end{array}$
The echo-attenuation for the single PGSE is given by

(54) $\begin{array}{cc}E\left(b\right)=\frac{\sqrt{\pi }}{2}\frac{e^{-{\mathrm{bD}}_{\perp }}}{\sqrt{b\left(D_{\parallel }+D_{\perp }\right)}}\mathrm{erf}\left(\sqrt{b\left(D_{\parallel }+D_{\perp }\right)}\right).& \left(36\right)\end{array}$

(55) For a double PGSE with orthogonal encoding gradients, the distribution P(D) is given by

(56) $\begin{array}{cc}P\left(D\right)=\frac{1}{\sqrt{\left(D_{\parallel }+D_{\perp }-2D\right)\left(D_{\parallel }-D_{\perp }\right)}},& \left(37\right)\end{array}$
with the same mean value as for the single PGSE but with a reduced variance,

(57) 0$\begin{array}{cc}{\mu }_2=\frac{1}{45}{\left(D_{\perp }-D_{\parallel }\right)}^2.& \left(38\right)\end{array}$
As in the single PGSE, also in double PGSE the echo-attenuation exhibits multi-component decay,

(58) $\begin{array}{cc}E\left(b\right)=\frac{\sqrt{\pi }}{2}\frac{e^{-b\frac{D_{\perp }+D_{\parallel }}{2}}}{\sqrt{b\frac{D_{\perp }-D_{\parallel }}{2}}}\mathrm{erf}\left(\sqrt{b\frac{D_{\perp }-D_{\parallel }}{2}}\right).& \left(39\right)\end{array}$
For randomly oriented anisotropic domains, the non-isotropic diffusion weighting results in a relatively broad distribution of diffusion coefficients, although reduced four-fold when measured with a double PGSE compared to the single PGSE. On the other hand the isotropic weighting results in

(59) $\begin{array}{cc}P\left(D\right)=\delta \left(D-\frac{D_{\parallel }+2D_{\perp }}{3}\right),\mathrm{with}& \left(40\right)\\ {\mu }_2=0& \left(41\right)\end{array}$
and a mono-exponential signal decay
E(b)=e.sup.−bD.  (42)

(60) The variance μ.sub.2 could be estimated by applying a function of the form (33) to fitting the echo attenuation data. However, in case of randomly oriented anisotropic domains, the convergence of the cumulant expansion of (36) is slow, thus several cumulants may be needed to adequately describe the echo attenuation (36). Alternatively, the distribution (34) may be approximated with the Gamma distribution

(61) $\begin{array}{cc}P\left(D\right)=D^{\alpha -1}\frac{e^{-D/\beta }}{\Gamma \left(\alpha \right){\beta }^{\alpha }},& \left(43\right)\end{array}$
where α is known as the shape parameter and β is known as the scale parameter. For the Gamma distribution, the mean diffusivity is given by D=α.Math.β, while the variance is given by μ.sub.2=α.Math.β.sup.2. The Gamma distribution is an efficient fitting function. With the two parameters it can capture a wide range of diffusion distributions, with both isotropic as well as anisotropic contributions. Conveniently, the Laplace transform of the Gamma function takes a simple analytical form,

(62) $\begin{array}{cc}\begin{array}{c}E\left(b\right)={\left(1+b\beta \right)}^{-\alpha }\\ ={\left(1+b\frac{{\mu }_2}{\stackrel{\_}{D}}\right)}^{-\frac{{\stackrel{\_}{D}}^2}{{\mu }_2}}.\end{array}& \left(44\right)\end{array}$

(63) The variance, μ.sub.2.sup.iso, obtained by fitting the function (44) to the isotropic diffusion weighted echo-decay is related to the isotropic diffusion contributions, since the variance is expected to vanish with isotropic weighting in a pure microscopically anisotropic system (see Eq. 41). The same fitting procedure on non-isotropically weighted data will yield the variance μ.sub.2 due to both isotropic and anisotropic contributions. The difference μ.sub.2−μ.sub.2.sup.iso vanishes when all diffusion contributions are isotropic and therefore provides a measure of micro-anisotropy. The mean diffusivity D, on the other hand, is expected to be identical for both isotropically and non-isotropically weighted data. The difference μ.sub.2−μ.sub.2.sup.iso is thus obtained by using the μ.sub.2.sup.iso and μ.sub.2 as free fit parameters when Eq. (44) is fitted to isotropically and non-isotropically weighted data sets, respectively, while a common parameter D is used to fit both data sets.

(64) The difference μ.sub.2−μ.sub.2.sup.iso along with D provide a novel measure for the microscopic fractional anisotropy (μFA) as

(65) $\begin{array}{cc}\mu \mathrm{FA}=\sqrt{\frac{3}{2}}\sqrt{\frac{{\mu }_2-{\mu }_2^{\mathrm{iso}}}{{\mu }_2-{\mu }_2^{\mathrm{iso}}+2{\stackrel{\_}{D}}^2/5}}.& \left(45\right)\end{array}$
The μFA is defined so that the μFA values correspond to the values of the well-established FA when diffusion is locally purely anisotropic and determined by randomly oriented axially symmetric diffusion tensors with identical eigenvalues. Eq. (45) is obtained by setting μFA=FA (27), assuming μ.sub.2−μ.sub.2.sup.iso=μ.sub.2 and expressing the eigenvalues D.sub.∥ and D.sub.⊥ in terms of D and μ.sub.2 (see Eq. 35). In the case of a one-dimensional curvilinear diffusion, when D.sub.∥>>D.sub.⊥, μFA=FA=1 and in the case of two-dimensional curvilinear diffusion, when D.sub.∥<<D.sub.⊥, μFA=FA=1/√{square root over (2)}.

(66) The difference μ.sub.2−μ.sub.2.sup.iso in Eq. (45) ensures that the micro-anisotropy can be quantified even when isotropic diffusion components are present. Isotropic restrictions, e.g. spherical cells, characterised by non-Gaussian restricted diffusion, are expected to cause a relative increase of both μ.sub.2 and μ.sub.2−μ.sub.2.sup.iso by the same amount, thus providing the difference μ.sub.2−μ.sub.2.sup.iso independent of the amount of isotropic contributions. The amount of non-Gaussian contributions could be quantified for example as the ratio √{square root over (μ.sub.2.sup.iso)}/D

(67) For anisotropic diffusion with finite orientation dispersion, i.e. when local diffusion tensors are not completely randomly oriented, the D and μ.sub.2−μ.sub.2.sup.iso are expected to depend on the gradient orientation in the non-isotropic diffusion weighting experiment. Furthermore, variation of the apparent diffusion coefficient (ADC), i.e. volume weighted average diffusivity, dependent on the gradient orientation and given by the initial echo decay of the non-isotropic diffusion weighting experiment, may indicate a finite orientation dispersion. Non-isotropic weighting experiment performed in several directions, similar to the diffusion tensor and diffusion kurtosis tensor measurements, performed with a range of b values to detect possibly multi-exponential decays, combined with the isotropic weighting experiment, is thus expected to yield additional information about micro-anisotropy as well as information related to the orientation dispersion of anisotropic domains.

(68) Eq. (44) could be expanded by additional terms in cases where this is appropriate. For example, the effects of a distinct signal contribution by the cerebrospinal fluid (CSF) in brain could be described by adding a mono-exponential term with the isotropic CSF diffusivity D.sub.1 to Eq. (44),

(69) $\begin{array}{cc}E\left(b\right)=f{e}^{-{\mathrm{bD}}_1}+\left(1-f\right){\left(1+b\frac{{\mu }_2}{\stackrel{\_}{D}}\right)}^{-\frac{{\stackrel{\_}{D}}^2}{{\mu }_2}},& \left(46\right)\end{array}$
where f is the fraction of the additional signal contribution. Eq. (46) could be used instead of Eq. (44) to fit the experimental data.

(70) When an extended fitting model described in Eq. (46) is applied, then the mean diffusivity, D, the additional diffusion contribution (f) and the diffusivity of the additional contribution (D.sub.1) are constrained to be equal for the isotropic and the non-isotropic diffusion weighted data.

(71) The method may involve the use of additional terms in Eq. (44), such as Eq. (46), applied to the analysis described in the above paragraphs. Eq. (46) comprises two additional parameters, i.e. fraction of the additional diffusion contribution (f) and diffusivity of the additional contribution (D.sub.1). One such example may be the analysis of data from the human brain, where the additional term in Eq. (46) could be assigned to the signal from the cerebrospinal fluid (CSF). The parameter D in Eq. (46) would in this case be assigned to the mean diffusivity in tissue while the parameter D.sub.1 would be assigned to the diffusivity of the CSF. The isotropic diffusion weighting could thus be used to obtain the mean diffusivity in the brain tissue without the contribution of the CSF.

(72) In addition, valuable information about anisotropy may be obtained from the ratio of the non-isotropically and the isotropically weighted signal or their logarithms. For example, the ratio of the non-isotropically and the isotropically weighted signals at intermediate b-values, might be used to estimate the difference between the radial (D.sub.⊥) and the axial (D.sub.∥) diffusivity in the human brain tissue due to the diffusion restriction effect by the axons. Extracting the information about microscopic anisotropy from the ratio of the signals might be advantageous, because the isotropic components with high diffusivity, e.g. due to the CSF, are suppressed at higher b-values. Such a signal ratio or any parameters derived from it might be used for generating parameter maps in MRI or for generating MR image contrast.

DETAILED DESCRIPTION OF THE DRAWINGS

(73) FIGS. 1 to 6 show examples of different gradient modulation schemes for isotropic diffusion weighting according to the present invention. In all of the FIGS. 1-6 the following is valid: A) Normalized dephasing magnitude F(t) (solid line), components of the normalized dephasing vector, q.sub.x/|q| (dashed line), q.sub.y/|q| (dotted line) and q.sub.z/|q| (dash dotted line). B) Azimuth angle ψ(t). C) Components of the normalized effective gradient vector, g.sub.x/|g| (dashed line), g.sub.y/|g| (dotted line) and g.sub.z/|g| (dash dotted line). Note that if a 180° RF pulse is used at t=t.sub.E/2, the actual hardware generated gradients are inverted compared to the ones shown in C) for times t>t.sub.E/2. D) The anisotropic weighting contributions from Eq. (16) as a function of time; the first term in Eq. (16) is shown as a dotted line, the second term is shown as a dashed dotted line, the third terms as a solid line and the fourth term is shown as a dashed line. The different presented gradient modulation schemes were constructed by first choosing the dephasing magnitude modulation, F(t), then calculating the corresponding time-dependent azimuth angle, ψ(t), followed by the calculation of the different components of the dephasing and gradient vectors. Note that, in this particular example, due to the choice of the time independent magic angle ζ.sub.m and the orientation of the laboratory axis, the z-components of the effective gradient vector and the dephasing vector are proportional to |g(t)| and F(t), respectively. This suggests that equivalent diffusion weighting values, b, can be achieved in an isotropic diffusion weighting experiment, utilizing gradients in all the three directions, and in a non-isotropic diffusion weighting experiment, utilizing only gradients in z-direction, if the z-gradient for non-isotropic diffusion weighting is larger than the z-gradient for isotropic diffusion weighting by factor √{square root over (3)}.

(74) The first example depicts the PGSE sequence with approximately constant F(t)=1, i.e. short z-gradient pulses (g.sub.z/|g|) at the beginning and at the end of the diffusion encoding interval. The gradient sequence is augmented by a sinusoidal gradient modulation in x-direction and with a cosine modulation in y-direction to achieve isotropic diffusion weighting. Note that, as in typical PGSE diffusion experiments, the non-isotropic diffusion weighting is achieved when x and y gradients are not active. In this example, the gradient modulations are identical in the intervals 0<t<t.sub.E/2 and t.sub.E/2<t<t.sub.E, when a 180° refocusing RF pulse is used, which is a preferred implementation for many applications, e.g. to achieve spectroscopic resolution. This may be advantageous due to possible asymmetries in gradient generating equipment. However, the use of short gradient pulses as well as the need to quickly increase the cosine gradient component to its maximum value following excitation and following the possible application of a 180° RF pulse, as well as quickly decrease its value to zero before a possible application of a 180° RF pulse, may be a disadvantageous implementation for some applications.

(75) The second example may be viewed as a PGSE with long gradient pulses in z-direction or a spin-echo experiment in a constant z-gradient (which may be provided by a stray field of the magnet) augmented with the gradient modulation in x and y directions for isotropic diffusion weighting. Similarly as in the first example, the possible need for fast rising and vanishing of some of the gradient components may be disadvantageous also in this case. Furthermore, unlike in the first example, modulations of some gradient components are not identical in the intervals 0<t<t.sub.E/2 and t.sub.E/2<t<t.sub.E.

(76) In relation to the description above and below it should be mentioned that also multi-echo variants of course are possible according to the present invention. Such may in some cases be beneficial for flow/motion compensation and for compensation of possible asymmetry in gradient generating equipment.

(77) In examples 3-6, we make use of harmonic gradient modulations for all gradient and dephasing components. These examples may be advantageous compared to the first two examples by using a more gradual variation in the dephasing magnitude. However, these examples do suffer from non-identical modulations of some gradient components in the intervals 0<t<t.sub.E/2 and t.sub.E/2<t<t.sub.E. While in examples 3-5 there may be the need for fast rising and vanishing of some of the gradient components immediately after and before the application of the RF pulses, the situation is more favourable in the sixth example, since all the gradient components conveniently vanish at times 0, t.sub.E/2 and t.sub.E. As my be understood from above, according to one specific embodiment of the present invention, the time-dependent normalized magnitude F(t) is chosen as a harmonic function of time. It should, however, be noted that this is not a must, as may be seen in FIGS. 1 and 2, where this is not the case.

(78) In FIG. 7A-C there is shown a schematic representation of signal decays vs. b for isotropic and non-isotropic diffusion weighting for different types of materials. In FIG. 7 the following is valid: A) Solid lines represent decays in a non-isotropic diffusion weighting experiment for 1D and 2D curvilinear diffusion (e.g. diffusion in reversed hexagonal phase H2 (tubes) and in lamellarphase Lα (planes), respectively). Dashed lines are the corresponding decays with isotropic diffusion weighting. The initial decay (D) is identical for the isotropic weighting as for the non-isotropic diffusion weighting. B) The decay for a system with 70% free isotropic diffusion and 30% restricted isotropic diffusion. In this case the isotropic and non-isotropic diffusion weighting result in identical signal decays in the entire b range. C) Decays for a system with 70% anisotropic diffusion (2D) and 30% restricted isotropic diffusion. Solid line corresponds to the non-isotropic diffusion weighting while the dashed line corresponds to the isotropic diffusion weighting. The initial decays are identical for the isotropic and for the non-isotropic diffusion weighting, while the deviation between the decays at higher b values reveals the degree of anisotropy.

(79) In FIG. 8A-C are shown experimental results with analysis of micro-anisotropy for different types of materials. Shown are normalized signal decays vs. bD for isotropic (circles) and non-isotropic (crosses) diffusion weighting. Solid lines represent optimal fits of Eq. (44) to the experimental data, with constraint of equal initial decays (shown as dashed lines) for isotropic and non-isotropic weighted data. All experiments were performed at 25° C. In all experiments, signal intensities were obtained by integration of the water peak. A) free water; data from the isotropic and non-isotropic diffusion weighting overlap and give rise to mono-exponential signal decays. The analysis gives D=2.2×10.sup.−9 m.sup.2/s and μLFA=0. B) Suspension of yeast cells from baker's yeast (Jästbolaget AB, Sweden) in tap water with restricted water diffusion; data from the isotropic and non-isotropic diffusion weighting overlap and give rise to multi-exponential signal decays. The analysis gives D=1.4×10.sup.−9 m.sup.2/s and μFA=0. C) Diffusion of water in a liquid crystal material composed by the Pluronic surfactant E5P68E6 with very high microscopic anisotropy, corresponding to a reverse hexagonal phase; data from the isotropic and non-isotropic diffusion weighting diverge at higher b-values and give rise to multi-exponential signal decay in case of the non-isotropic diffusion weighting and mono-exponential signal decay in case of the isotropic diffusion weighting. The analysis gives D=4.8×10.sup.−10 m.sup.2/s and μFA=1.0.

(80) In FIGS. 9A and 9B, the results of the Monte-Carlo error analysis show systematic deviations and precision of the D (A) and μFA (B) parameters estimated for the 1D (dots) and 2D (circles) curvilinear diffusion according to what has been disclosed above. The ratio of the estimated mean diffusivity to the exact values D, labelled as D/D (A) with the corresponding standard deviation values and the estimated μFA values (B) with the corresponding standard deviations are shown as dots/circles and error bars, respectively, as a function of the maximum attenuation factor b.sub.maxD for signal to noise level of 30.

(81) For μFA estimation, the optimal choice of the b-values is important. To investigate the optimal range of b-values, a Monte-Carlo error analysis depicted in FIGS. 9A and 9B has been performed. The echo-signal was generated as a function of 16 equally spaced b-values between 0 and b.sub.max for the cases of 1D and 2D curvilinear diffusion with randomly oriented domains. The upper limit, b.sub.max, was varied and the attenuation factors bD were chosen to be identical for the 1D and 2D case. The signal was subjected to the Rician noise with a constant signal to noise, SNR=30, determined relative to the non-weighted signal. Isotropic and non-isotropic weighed attenuation data were analysed with the protocol described herein to obtain D and μFA parameters. This analysis was repeated in 1000 iterations by adding different simulated noise signals with the given SNR. The procedure yields the mean and the standard deviation of the estimated D and μFA, shown as dots/circles and error bars respectively in FIG. 9B.

(82) The optimal range of the diffusion weighting b is given by a compromise between accuracy and precision of the μFA analysis and it depends on the mean diffusivity. If the maximum b value used is lower than 1/D, the μFA tends to be underestimated, while for maximum b values larger than 1/D the μFA tends to be overestimated. On the other hand the accuracy of μFA is compromised particularly at too low values of the maximum b, due to increased sensitivity to noise. See FIG. 9B.