Diffusion model

Given the concentration ${\displaystyle \rho }$ and flux ${\displaystyle J}$, Fick's first law gives a relationship between the flux and the concentration gradient:

${\displaystyle J(x,t)=-D\nabla \rho (x,t)}$

where D is the diffusion coefficient. Then, given conservation of mass, the continuity equation relates the time derivative of the concentration with the divergence of the flux:

${\displaystyle {\frac {\partial \rho (x,t)}{\partial t}}=-\nabla \cdot J(x,t)}$

Putting the two together, we get the diffusion equation:

${\displaystyle {\frac {\partial \rho (x,t)}{\partial t}}=D\nabla ^{2}\rho (x,t).}$

Magnetization dynamics

With no diffusion present, the change in nuclear magnetization over time is given by the classical Bloch equation

${\displaystyle {\frac {d{\vec {M}}}{dt}}=\gamma {\vec {M}}\times {\vec {B}}-{\frac {M_{x}{\vec {i}}+M_{y}{\vec {j}}}{T_{2}}}-{\frac {(M_{z}-M_{0}){\vec {k}}}{T_{1}}}}$

which has terms for precession, T2 relaxation, and T1 relaxation.

In 1956, H.C. Torrey mathematically showed how the Bloch equations for magnetization would change with the addition of diffusion.^[22] Torrey modified Bloch's original description of transverse magnetization to include diffusion terms and the application of a spatially varying gradient. Since the magnetization ${\displaystyle M}$ is a vector, there are 3 diffusion equations, one for each dimension. The Bloch-Torrey equation is:

${\displaystyle {\frac {d{\vec {M}}}{dt}}=\gamma {\vec {M}}\times {\vec {B}}-{\frac {M_{x}{\vec {i}}+M_{y}{\vec {j}}}{T_{2}}}-{\frac {(M_{z}-M_{0}){\vec {k}}}{T_{1}}}+\nabla \cdot {\vec {D}}\nabla {\vec {M}}}$

where ${\displaystyle {\vec {D}}}$ is now the diffusion tensor.

For the simplest case where the diffusion is isotropic the diffusion tensor is a multiple of the identity:

${\displaystyle {\vec {D}}=D\cdot {\vec {I}}=D\cdot {\begin{bmatrix}1&0&0\\0&1&0\\0&0&1\end{bmatrix}},}$

then the Bloch-Torrey equation will have the solution

${\displaystyle {M}={M}_{\text{bloch}}e^{-{\frac {1}{3}}\gamma ^{2}G^{2}t^{3}D}\sim e^{-bD_{0}}}$

The exponential term will be referred to as the attenuation ${\displaystyle A}$. Anisotropic diffusion will have a similar solution for the diffusion tensor, except that what will be measured is the apparent diffusion coefficient (ADC). In general, the attenuation is:

${\displaystyle A=e^{-\sum _{i,j}b_{ij}D_{ij}}}$

where the ${\displaystyle b_{ij}}$ terms incorporate the gradient fields ${\displaystyle G_{x}}$, ${\displaystyle G_{y}}$, and ${\displaystyle G_{z}}$.

Grayscale

The standard grayscale of DWI images is to represent increased diffusion restriction as brighter.^[23]

Mathematical foundation—tensors

Diffusion MRI relies on the mathematics and physical interpretations of the geometric quantities known as tensors. Only a special case of the general mathematical notion is relevant to imaging, which is based on the concept of a symmetric matrix.^{[notes 1]} Diffusion itself is tensorial, but in many cases the objective is not really about trying to study brain diffusion per se, but rather just trying to take advantage of diffusion anisotropy in white matter for the purpose of finding the orientation of the axons and the magnitude or degree of anisotropy. Tensors have a real, physical existence in a material or tissue so that they do not move when the coordinate system used to describe them is rotated. There are numerous different possible representations of a tensor (of rank 2), but among these, this discussion focuses on the ellipsoid because of its physical relevance to diffusion and because of its historical significance in the development of diffusion anisotropy imaging in MRI.

The following matrix displays the components of the diffusion tensor:

${\displaystyle {\bar {D}}={\begin{vmatrix}D_{\color {red}xx}&D_{xy}&D_{xz}\\D_{xy}&D_{\color {red}yy}&D_{yz}\\D_{xz}&D_{yz}&D_{\color {red}zz}\end{vmatrix}}}$

The same matrix of numbers can have a simultaneous second use to describe the shape and orientation of an ellipse and the same matrix of numbers can be used simultaneously in a third way for matrix mathematics to sort out eigenvectors and eigenvalues as explained below.

Physical tensors

The idea of a tensor in physical science evolved from attempts to describe the quantity of physical properties. The first properties they were applied to were those that can be described by a single number, such as temperature. Properties that can be described this way are called scalars; these can be considered tensors of rank 0, or 0th-order tensors. Tensors can also be used to describe quantities that have directionality, such as mechanical force. These quantities require specification of both magnitude and direction, and are often represented with a vector. A three-dimensional vector can be described with three components: its projection on the x, y, and z axes. Vectors of this sort can be considered tensors of rank 1, or 1st-order tensors.

A tensor is often a physical or biophysical property that determines the relationship between two vectors. When a force is applied to an object, movement can result. If the movement is in a single direction, the transformation can be described using a vector—a tensor of rank 1. However, in a tissue, diffusion leads to movement of water molecules along trajectories that proceed along multiple directions over time, leading to a complex projection onto the Cartesian axes. This pattern is reproducible if the same conditions and forces are applied to the same tissue in the same way. If there is an internal anisotropic organization of the tissue that constrains diffusion, then this fact will be reflected in the pattern of diffusion. The relationship between the properties of driving force that generate diffusion of the water molecules and the resulting pattern of their movement in the tissue can be described by a tensor. The collection of molecular displacements of this physical property can be described with nine components—each one associated with a pair of axes xx, yy, zz, xy, yx, xz, zx, yz, zy.^[30] These can be written as a matrix similar to the one at the start of this section.

Diffusion from a point source in the anisotropic medium of white matter behaves in a similar fashion. The first pulse of the Stejskal Tanner diffusion gradient effectively labels some water molecules and the second pulse effectively shows their displacement due to diffusion. Each gradient direction applied measures the movement along the direction of that gradient. Six or more gradients are summed to get all the measurements needed to fill in the matrix, assuming it is symmetric above and below the diagonal (red subscripts).

In 1848, Henri Hureau de Sénarmont^[31] applied a heated point to a polished crystal surface that had been coated with wax. In some materials that had "isotropic" structure, a ring of melt would spread across the surface in a circle. In anisotropic crystals the spread took the form of an ellipse. In three dimensions this spread is an ellipsoid. As Adolf Fick showed in the 1850s, diffusion exhibits many of the same patterns as those seen in the transfer of heat.

Mathematics of ellipsoids

At this point, it is helpful to consider the mathematics of ellipsoids. An ellipsoid can be described by the formula: ${\displaystyle ax^{2}+by^{2}+cz^{2}=1}$. This equation describes a quadric surface. The relative values of a, b, and c determine if the quadric describes an ellipsoid or a hyperboloid.

As it turns out, three more components can be added as follows: ${\displaystyle ax^{2}+by^{2}+cz^{2}+dyz+ezx+fxy=1}$. Many combinations of a, b, c, d, e, and f still describe ellipsoids, but the additional components (d, e, f) describe the rotation of the ellipsoid relative to the orthogonal axes of the Cartesian coordinate system. These six variables can be represented by a matrix similar to the tensor matrix defined at the start of this section (since diffusion is symmetric, then we only need six instead of nine components—the components below the diagonal elements of the matrix are the same as the components above the diagonal). This is what is meant when it is stated that the components of a matrix of a second order tensor can be represented by an ellipsoid—if the diffusion values of the six terms of the quadric ellipsoid are placed into the matrix, this generates an ellipsoid angled off the orthogonal grid. Its shape will be more elongated if the relative anisotropy is high.

When the ellipsoid/tensor is represented by a matrix, we can apply a useful technique from standard matrix mathematics and linear algebra—that is to "diagonalize" the matrix. This has two important meanings in imaging. The idea is that there are two equivalent ellipsoids—of identical shape but with different size and orientation. The first one is the measured diffusion ellipsoid sitting at an angle determined by the axons, and the second one is perfectly aligned with the three Cartesian axes. The term "diagonalize" refers to the three components of the matrix along a diagonal from upper left to lower right (the components with red subscripts in the matrix at the start of this section). The variables ${\displaystyle ax^{2}}$, ${\displaystyle by^{2}}$, and ${\displaystyle cz^{2}}$ are along the diagonal (red subscripts), but the variables d, e and f are "off diagonal". It then becomes possible to do a vector processing step in which we rewrite our matrix and replace it with a new matrix multiplied by three different vectors of unit length (length=1.0). The matrix is diagonalized because the off-diagonal components are all now zero. The rotation angles required to get to this equivalent position now appear in the three vectors and can be read out as the x, y, and z components of each of them. Those three vectors are called "eigenvectors" or characteristic vectors. They contain the orientation information of the original ellipsoid. The three axes of the ellipsoid are now directly along the main orthogonal axes of the coordinate system so we can easily infer their lengths. These lengths are the eigenvalues or characteristic values.

Diagonalization of a matrix is done by finding a second matrix that it can be multiplied with followed by multiplication by the inverse of the second matrix—wherein the result is a new matrix in which three diagonal (xx, yy, zz) components have numbers in them but the off-diagonal components (xy, yz, zx) are 0. The second matrix provides eigenvector information.

Measures of anisotropy and diffusivity

Visualization of DTI data with ellipsoids.

In present-day clinical neurology, various brain pathologies may be best detected by looking at particular measures of anisotropy and diffusivity. The underlying physical process of diffusion causes a group of water molecules to move out from a central point, and gradually reach the surface of an ellipsoid if the medium is anisotropic (it would be the surface of a sphere for an isotropic medium). The ellipsoid formalism functions also as a mathematical method of organizing tensor data. Measurement of an ellipsoid tensor further permits a retrospective analysis, to gather information about the process of diffusion in each voxel of the tissue.^[32]

In an isotropic medium such as cerebrospinal fluid, water molecules are moving due to diffusion and they move at equal rates in all directions. By knowing the detailed effects of diffusion gradients we can generate a formula that allows us to convert the signal attenuation of an MRI voxel into a numerical measure of diffusion—the diffusion coefficient D. When various barriers and restricting factors such as cell membranes and microtubules interfere with the free diffusion, we are measuring an "apparent diffusion coefficient", or ADC, because the measurement misses all the local effects and treats the attenuation as if all the movement rates were solely due to Brownian motion. The ADC in anisotropic tissue varies depending on the direction in which it is measured. Diffusion is fast along the length of (parallel to) an axon, and slower perpendicularly across it.

Once we have measured the voxel from six or more directions and corrected for attenuations due to T2 and T1 effects, we can use information from our calculated ellipsoid tensor to describe what is happening in the voxel. If you consider an ellipsoid sitting at an angle in a Cartesian grid then you can consider the projection of that ellipse onto the three axes. The three projections can give you the ADC along each of the three axes ADC_x, ADC_y, ADC_z. This leads to the idea of describing the average diffusivity in the voxel which will simply be

${\displaystyle (ADC_{x}+ADC_{y}+ADC_{z})/3=ADC_{i}}$

We use the i subscript to signify that this is what the isotropic diffusion coefficient would be with the effects of anisotropy averaged out.

The ellipsoid itself has a principal long axis and then two more small axes that describe its width and depth. All three of these are perpendicular to each other and cross at the center point of the ellipsoid. We call the axes in this setting eigenvectors and the measures of their lengths eigenvalues. The lengths are symbolized by the Greek letter λ. The long one pointing along the axon direction will be λ₁ and the two small axes will have lengths λ₂ and λ₃. In the setting of the DTI tensor ellipsoid, we can consider each of these as a measure of the diffusivity along each of the three primary axes of the ellipsoid. This is a little different from the ADC since that was a projection on the axis, while λ is an actual measurement of the ellipsoid we have calculated.

The diffusivity along the principal axis, λ₁ is also called the longitudinal diffusivity or the axial diffusivity or even the parallel diffusivity λ_∥. Historically, this is closest to what Richards originally measured with the vector length in 1991.^[33] The diffusivities in the two minor axes are often averaged to produce a measure of radial diffusivity

${\displaystyle \lambda _{\perp }=(\lambda _{2}+\lambda _{3})/2.}$

This quantity is an assessment of the degree of restriction due to membranes and other effects and proves to be a sensitive measure of degenerative pathology in some neurological conditions.^[34] It can also be called the perpendicular diffusivity (${\displaystyle \lambda _{\perp }}$).

Another commonly used measure that summarizes the total diffusivity is the Trace—which is the sum of the three eigenvalues,

${\displaystyle \mathrm {tr} (\Lambda )=\lambda _{1}+\lambda _{2}+\lambda _{3}}$

where ${\displaystyle \Lambda }$ is a diagonal matrix with eigenvalues ${\displaystyle \lambda _{1}}$, ${\displaystyle \lambda _{2}}$ and ${\displaystyle \lambda _{3}}$ on its diagonal.

If we divide this sum by three we have the mean diffusivity,

${\displaystyle \mathrm {MD} =(\lambda _{1}+\lambda _{2}+\lambda _{3})/3}$

which equals ADC_i since

${\displaystyle {\begin{aligned}\mathrm {tr} (\Lambda )/3&=\mathrm {tr} (V^{-1}V\Lambda )/3\\&=\mathrm {tr} (V\Lambda V^{-1})/3\\&=\mathrm {tr} (D)/3\\&=ADC_{i}\end{aligned}}}$

where ${\displaystyle V}$ is the matrix of eigenvectors and ${\displaystyle D}$ is the diffusion tensor. Aside from describing the amount of diffusion, it is often important to describe the relative degree of anisotropy in a voxel. At one extreme would be the sphere of isotropic diffusion and at the other extreme would be a cigar or pencil shaped very thin prolate spheroid. The simplest measure is obtained by dividing the longest axis of the ellipsoid by the shortest = (λ₁/λ₃). However, this proves to be very susceptible to measurement noise, so increasingly complex measures were developed to capture the measure while minimizing the noise. An important element of these calculations is the sum of squares of the diffusivity differences = (λ₁ − λ₂)² + (λ₁ − λ₃)² + (λ₂ − λ₃)². We use the square root of the sum of squares to obtain a sort of weighted average—dominated by the largest component. One objective is to keep the number near 0 if the voxel is spherical but near 1 if it is elongate. This leads to the fractional anisotropy or FA which is the square root of the sum of squares (SRSS) of the diffusivity differences, divided by the SRSS of the diffusivities. When the second and third axes are small relative to the principal axis, the number in the numerator is almost equal the number in the denominator. We also multiply by ${\displaystyle 1/{\sqrt {2}}}$ so that FA has a maximum value of 1. The whole formula for FA looks like this:

${\displaystyle \mathrm {FA} ={\frac {\sqrt {3((\lambda _{1}-\operatorname {E} [\lambda ])^{2}+(\lambda _{2}-\operatorname {E} [\lambda ])^{2}+(\lambda _{3}-\operatorname {E} [\lambda ])^{2})}}{\sqrt {2(\lambda _{1}^{2}+\lambda _{2}^{2}+\lambda _{3}^{2})}}}}$

The fractional anisotropy can also be separated into linear, planar, and spherical measures depending on the "shape" of the diffusion ellipsoid.^[35][36] For example, a "cigar" shaped prolate ellipsoid indicates a strongly linear anisotropy, a "flying saucer" or oblate spheroid represents diffusion in a plane, and a sphere is indicative of isotropic diffusion, equal in all directions.^[37] If the eigenvalues of the diffusion vector are sorted such that ${\displaystyle \lambda _{1}\geq \lambda _{2}\geq \lambda _{3}\geq 0}$, then the measures can be calculated as follows:

For the linear case, where ${\displaystyle \lambda _{1}\gg \lambda _{2}\simeq \lambda _{3}}$,

${\displaystyle C_{l}={\frac {\lambda _{1}-\lambda _{2}}{\lambda _{1}+\lambda _{2}+\lambda _{3}}}}$

For the planar case, where ${\displaystyle \lambda _{1}\simeq \lambda _{2}\gg \lambda _{3}}$,

${\displaystyle C_{p}={\frac {2(\lambda _{2}-\lambda _{3})}{\lambda _{1}+\lambda _{2}+\lambda _{3}}}}$

For the spherical case, where ${\displaystyle \lambda _{1}\simeq \lambda _{2}\simeq \lambda _{3}}$,

${\displaystyle C_{s}={\frac {3\lambda _{3}}{\lambda _{1}+\lambda _{2}+\lambda _{3}}}}$

Each measure lies between 0 and 1 and they sum to unity. An additional anisotropy measure can used to describe the deviation from the spherical case:

${\displaystyle C_{a}=C_{l}+C_{p}=1-C_{s}={\frac {\lambda _{1}+\lambda _{2}-2\lambda _{3}}{\lambda _{1}+\lambda _{2}+\lambda _{3}}}}$

There are other metrics of anisotropy used, including the relative anisotropy (RA):

${\displaystyle \mathrm {RA} ={\frac {\sqrt {(\lambda _{1}-\operatorname {E} [\lambda ])^{2}+(\lambda _{2}-\operatorname {E} [\lambda ])^{2}+(\lambda _{3}-\operatorname {E} [\lambda ])^{2}}}{{\sqrt {3}}\operatorname {E} [\lambda ]}}}$

and the volume ratio (VR):

${\displaystyle \mathrm {VR} ={\frac {\lambda _{1}\lambda _{2}\lambda _{3}}{\operatorname {E} [\lambda ]^{3}}}}$

Summary

For DTI, it is generally possible to use linear algebra, matrix mathematics and vector mathematics to process the analysis of the tensor data.

In some cases, the full set of tensor properties is of interest, but for tractography it is usually necessary to know only the magnitude and orientation of the primary axis or vector. This primary axis—the one with the greatest length—is the largest eigenvalue and its orientation is encoded in its matched eigenvector. Only one axis is needed as it is assumed the largest eigenvalue is aligned with the main axon direction to accomplish tractography.