[

A Bloch-Torrey Equation for Diffusion in a Deforming Media

Damien Rohmer, Grant T. Gullberg
The authors are with the Division of Life Science, E.O. Lawrence National Berkeley Laboratory,
1 Cyclotron Road, Berkeley, CA 94720, USA.

]

Introduction

2DIFFUSION phenomenon is a macroscopic observation of the microscopic structure of the material in which it occurs. Knowledge of this diffusion in cases of various biological tissues such as cardiac muscle [1,2,3,4,5] , skeletal muscles [6,7,8,9], brain [10,11,12,13], tongue [14,15,16], spinal cord [17,18] or other [19] is therefore a probe into their intrinsic structures. Diffusion measurement using the Diffusion Tensor Magnetic Resonance Imaging (DTMRI) technique is a powerful tool for exploring the microscopic structure of tissues. This technique is well established for static cases and is widely used in brain imaging where the change in the diffusion coefficient is related to brain diseases [20,21,22]. Its utilization for other organs are currently in the research process and one of the most challenging applications is using the DTMRI technique on a living heart. The fast and complex motion of the heart [23] prevents accurate measurements of the diffusion parameters. Static acquisitions of diffusion measurements on an ex-vivo heart have already demonstrated the value of measuring diffusion in cardiac tissue for delineating the fibrous and laminar structures of cardiac muscle [1,24,25,26] and for the remodeling that occurs after heart diseases [27,28,29]. Some experiments have been done on the in-vivo heart using some approximations and configurations of the motion of the heart [30,31,32,33], however the problem of correlation between the results and real data is still an open problem, even for the case of static organs [22]. Knowing the mathematical behavior of the signal in a deforming media might help to interpret results obtained on the living heart. Simulations could then be performed to determine what parameters could be optimized to have a more accurate simulation of the DTMRI acquisition of the in-vivo cardiac system.
After presenting the mathematical context of the diffusion process, there will be a review of the method of diffusion acquisition in the static and dynamic case. Then the equation defining the behavior of the MRI signal in the context of a deforming media will be introduced, an application will be provided in the case of the prolate spheroidal coordinates [34,35], and the method of the numerical solution of this equation will be developed.

Diffusion Process

From the Diffusion Coefficient to the Tensor of Diffusion

When looking at the microscopic resolution, fluid molecules are in constant motion because of thermal agitation. This random motion was initially noticed by the Scottish naturalist Robert Brown in $1810$ when he recorded the oscillatory motion of pollen grains suspended in water. This observation remained a mystery until Albert Einstein's kinetic-molecular theory was established in $1905$ . It stated that each molecule is agitated due to thermodynamic effects and each follows a random path. At a macroscopic resolution, each small region of water diffuses around its initial position. This diffusion can occur from one substance into another but it can also occur from one region of water to another region of water. To relate the microscopic resolution to the macroscopic resolution, Einstein introduced the diffusion coefficient $D$ such that in the one dimensional case

$$6\tau\,D=\left[x\left(t+\tau\right)-x\left(t\right)\right]^2$$

where the molecule of water is positioned at the position $x(t)$ at the time $t$ and moves to the position $x(t+\tau)$ at the time $t+\tau$ due to diffusion. In the case of three dimensional displacement for many molecules, the relative displacement of a molecule at a position ${\bf x}$ (where bold letter refers to vectors) is designated by the vector ${\bf u}(\tau)={\bf x}(t+\tau)-{\bf x}(t)$ , which is a function of time $\tau$ . The diffusion coefficient is therefore expressed as the dot product:

$$6\tau\,D={\bf u}(\tau)\cdot{\bf u}(\tau)$$

The characteristic length $l$ is then introduced such that $l(\tau)=\sqrt{6 D \tau}$ . This length is characteristic of the expansion taken by any small section of water during a time $\tau$ when diffusion occurs. If some molecules of water initially positioned at the origin are followed during a time $\tau$ , they will spread homogeneously around a circle of radius $l(\tau)$ .

Thus far, the diffusive media has been considered to be isotropic, meaning that diffusion proprieties are identical in all directions, which is true for fluids. In biological tissues (such as fibrous tissues) [36,37], polymers [38] or nematic liquid crystals (where all the molecules point in the same direction but do not share the same row formation) [39], material microstructures are aligned along some directions. This alignment enables the molecule of water to move more easily in the direction of the alignment than in the perpendicular one [22]. This diffusion differs from the ``intrinsic'' diffusion in that it depends on the interaction of the molecule of water with the physical structure of the material. This type of diffusion should therefore be called ``apparent'' diffusion or ``restricted'' diffusion and measures only the mobility of the molecule of water in the given physical structure (cellular structure in the case of the biological sample). The given coefficient $\mathrm{D}$ should therefore be called ``Apparent'' Diffusion Coefficient (ADC). However, the two terms are usually mixed to simplify the expression. In anisotropic diffusion, this diffusion depends on the orientation and an initial region of water will not spread as a spherical shape. Taking the hypothesis of a simple model where only one direction is promoted at each position, the diffusion tensor can be introduced by using the notation:

$$6\tau \mathrm{D} = {\bf u}\,{\bf u}^{T} = {\bf u}\otimes {\bf u}$$ (1)

where $\otimes$ represents the tensorial product. This coefficient $\mathrm{D}$ is therefore a tensor of rank two and is represented by a symetric $3\times 3$ positive definite matrix by its definition in Eq. 1. This matrix depends unfortunately on the basis in which the acquisition has been performed¹. This matrix is build as a quadratic form and maps the unit sphere into an ellipsoid. This ellipsoid can always be separated in three main directions by eigenvalues decomposition. The first direction is the elongated part of the ellipsoid at a distance $\lambda_1$ of the origin and in the direction ${\bf d}_1$ where ${\bf d}_1$ is a unitary vector. Along this direction ${\bf d}_1$ , the displacement of the molecules of water are the largest, therefore it is denoted as the Principal Direction of the apparent Diffusion (PDD). Two other directions can be differentiated using the same notation, the medium and the smallest side of the ellipsoid are then defined by the direction ${\bf d}_2$ and ${\bf d}_3$ with the respective size $\lambda_2$ and $\lambda_3$ . The diffusion tensor can therefore be decomposed in an appropriate basis where the tensor is diagonal

$$\mathrm{D} = \sum_{i=1}^{3}{\lambda_i\,{\bf d}_i\otimes{\bf d}_i} = R\,\Lambda\,R^{T}$$

where $\Lambda$ is the diagonal matrix $\mathrm{diag}(\lambda_1,\lambda_2,\lambda_3)$ with $\lambda_1\geq\lambda_2\geq\lambda_3$ . $R=({\bf d}_1,{\bf d}_2,{\bf d}_3)$ is the rotation matrix which transforms the initial Cartesian basis into the local basis of the ellipsoid. This visualization is shown in Fig. 1

Fig. 1: Visualization of the diffusion tensor as an ellipsoid with the corresponding basis vector.

However, this formalism associated with the tensor of diffusion reaches some limits as soon as the material exhibits a more complex configuration. If more than one direction is promoted, as in the case of a fiber crossing in the biological tissue such as with the white matter in the brain [40], the diffusion tensor is unable to solve the problem accurately. In this case, the principal direction of diffusion will be aligned with the average of the two crossing directions leading to a large error of interpretation as illustrated in Fig. 2. This problem is especially limiting in the case of real studies because of the partial volume effect. If large voxels are taken in account during the acquisition, the diffusion tensor has a high probability for averaging the directions that follow different paths inside this voxel. In order to have a more accurate model, higher order tensors can be used (q-ball) [41,42,43]. However, these methods still require a longer acquisition time and a strong gradient system [44] for a good signal to noise ratio (SNR). Moreover, for the purpose of the muscular fibers of the heart, the bundles are supposed to follow the same orientation in a local neighboring and therefore the direction of diffusion varies smoothly. For this reason, the tensor notation will be kept for the purpose of the current study.

Fig. 2: Example of fiber crossing. The red tube represents an example of fiber. The blue circle shows the intersection of the fiber crossing. With the tensor approximation, the principal direction of diffusion is the average between the two directions of the fiber.

The Diffusion Equation

Performing a numerical or mathematical study of a macroscopic structure taking into account every microscopic molecule is not feasible. The studied element, such as water, or spin orientation in the case of the MR signal, is studied as a continuous function. The macroscopic behavior is therefore known by equations that are linked to the microscopic structure via the diffusion tensor.

A scalar density of the diffusing material is introduced using the notation $\phi$ , which represents the water density or the spin orientation in one direction. This density depends on the position ${\bf x}$ and on the time $t$ . Denoting ${\bf j}$ the density flux, Fick's first law states that this flux is given by the product of the diffusion coefficient or tensor with the spatial variation of the density (the diffusion occurs from the highest concentration to the lowest):

$${\bf j}=-\mathrm{D}\,\nabla \phi$$ (2)

Therefore ${\bf j}$ is a vector oriented through the local lowest concentration at each position. The diffusion equation can be derived from the continuity equation. If no sink or sources of diffusing material exist, the time rate of change of the density equals the opposite of the divergence of the density flux:

$$\frac{\partial \phi}{\partial t}=-\nabla \cdot {\bf j}$$

And using Eq. 2 to express the flux,

$$\frac{\partial \phi}{\partial t}=\nabla\cdot\left(\mathrm{D}\nabla\,\phi\right)$$ (3)

In the case of a constant isotropic diffusion $\mathrm{D}=D\,\mathrm{I}$ , where $\mathrm{I}$ is the $3\times 3$ identity matrix and $D$ is the scalar coefficient of diffusion, Equation 3 simplifies then to the Laplace equation:

$$\frac{\partial \phi}{\partial t}=D\,\nabla\cdot\left(\nabla\,\phi\right)$$

$$\Rightarrow \frac{\partial \phi}{\partial t}=D\,\triangle \phi$$ (4)

In the case where the boundary conditions equal zero at infinity, the solution can be expressed analytically by the convolution of the initial condition with a Gaussian function:

$$\phi({\bf x},t)=\frac{1}{\sqrt{D}\left(4\pi\,t\right)^{\frac{N}{2}}}e^{-\frac{{\bf x}^T {\bf x}}{4\,D\,t}}\;\ast\,\phi({\bf x},0)$$ (5)

where $N$ is the number of dimensions of the problem (in the space ${\bf x}=(x,y,z)$ , $N=3$ ). The exponential function is called the Gaussian kernel of the equation and solves the problem for any initial condition given by $\phi$ . In the case of an anisotropic diffusion independent of time, this fundamental solution can also be expressed by [45]:

$$\phi({\bf x},t)=\frac{1}{\sqrt{\vert\mathrm{D}\vert}\left(4\pi\,t... ...}{2}}}e^{-\frac{{\bf x}^T \mathrm{D}^{-1}{\bf x}}{4\,t}}\;\ast\,\phi({\bf x},0)$$

In the case of a vectorial diffusive density, each component of the quantity $\boldsymbol \phi$ is then diffused in the material. In this case the equation of diffusion is given by the three scalar equation of diffusion:

$$\phi^{i}=\nabla\left(\mathrm{D}\nabla\,\phi^i\right)$$

where $\boldsymbol{\phi}^i$ is the $i^{th}$ component of the vectorial quantity $\boldsymbol{\phi}$ , with usually $i\in[\![1,3]\!]$ . To simplify the notation, the vectorial equation will be noted abusively:

$$\frac{\partial \boldsymbol{\phi}}{\partial t}=\nabla\left(\mathrm{D}\nabla\,\boldsymbol{\phi}\right)$$

where the first operator $\nabla$ acts on the scalar quantity. To further simplify, the notation ${\bf u}_{,v}=\partial_{v}{\bf u}=\frac{\partial {\bf u}}{\partial v}$ will be used:

$$\partial_{t}\,{\boldsymbol{\phi}}=\nabla\cdot\left(\mathrm{D}\nabla\,{\boldsymbol{\phi}}\right)$$ (6)

The case of inhomogeneous isotropic diffusion is most often used in filtering in image processing. An example is given by taking the original noisy picture and the result given by a normal Gaussian filter comparing them to the inhomogeneous diffusion in Fig. 3. The inhomogeneous diffusion coefficient is given by $D({\bf x})=D_0\,\frac{\lambda^2}{\lambda^2+\Vert\nabla f({\bf x})\Vert^2}$ .

Fig. 3: Examples of two dimensional diffusion filtering. The original image is on the top left and Gaussian noise was added to this to obtain the image on the top right. The two images at the bottom show the results of the diffusion filtering. The image on the left was obtained after an isotropic homogeneous diffusion corresponding to a convolution by a constant Gaussian kernel. The image on the right is the result of the inhomogeneous diffusion. The results were obtained numerically using the technique explained at the end of the report.

An example of anisotropic diffusion is also given. A scalar function is given with an initial value of zero or one. Then an inhomogeneous anisotropic diffusion is performed and shown in Fig. 4. The principal directions of the diffusion are plotted on the second picture.

Fig. 4: The first image shows the initial scalar density function composed of values zero and one. The second image shows the result of the anisotropic diffusion. The principal direction of diffusion depends on the position and is shown in the third image (the other direction of diffusion can be neglected for being $15$ times smaller).

Link between Signal Attenuation and Diffusion

DTMRI and Bloch-Torrey Equation

Diffusion Tensor Magnetic Resonance Imaging (DTMRI) provides a way to access the tissue microstructure by having knowledge of the diffusion with a large number of samples. DTMRI is a link to the classical MRI [46]. Atoms such as ${}^1\mathrm{H}$ , which are especially abundant in water and therefore in biological tissue, possess a nuclear spin angular momentum. These can be visualized as spinning charged spheres that give rise to a small magnetic moment. When a magnetic field ${\bf B}$ is applied to these spins, the magnetic moment vector tends to align along the direction of the field. The nuclear spins are in resonance at the pulsation called the Lamor frequency:

$$\omega = \gamma\,\Vert{\bf B}\Vert$$ (7)

where $\omega$ is the pulsation of the spin and $\gamma$ is the gyromagnetic ratio such that $\frac{\gamma}{2\pi}=42.58\,\mathrm{MHz/T}$ . By convention, the $z$ axis is aligned with the largest main magnetic field ${\bf B}_0$ and a Cartesian basis $({\bf e}_{1},{\bf e}_{2},{\bf e}_3)$ is introduced such that any vector ${\bf x}=\sum_{i=1}^{3}x^i{\bf e}_i$ and the $z$ axis is aligned with the vector ${\bf e}_3$ . In order to separate the spins from the huge magnetic field ${\bf B}_0$ ($1-9$ T), the radiofrequency field pulse ${\bf B}_1$ is applied at the Lamor frequency in the normal plane of ${\bf B}_0$ to excite the spins out of equilibrium and make them more aligned in this normal $({\bf e}_1,{\bf e}_2)$ plane. In order to create a dephasing between the spins, and therefore to differentiate them, a magnetic gradient ${\bf G}$ is added in the spatial direction during the reading process. Then the acquisition is performed and the spin signal coming from the component situated in the normal plane of ${\bf B}_0$ is acquired. The magnetization vector of each spin ${\bf M}$ is characterized by the Bloch equation [47]:

$${\bf M}_{,t} = {\bf M}\times \gamma{\bf B} - \frac{ M^{1}\,{\bf e}_1 + M^{2}\,{\bf e}_2 } {T_2} - \frac{M^3-M_0^3}{T_1}\,{\bf e}_3$$

where ${\bf M}=M^1{\bf e}_1+M^2{\bf e}_2+M^3{\bf e}_3$ , ${\bf M}(\bf x,t=0)={\bf M}_0$ and ($T_1$ ,$T_2$ ) are the relaxation coefficients which depend on the tissue. In this equation the total magnetic field ${\bf B}$ can be written as:

$${\bf B}({\bf x},t)={\bf B}_0({\bf x},t)+{\bf G}({\bf x},t)\cdot{\bf x}$$

In order to take into account the diffusion process, the Bloch equation is modified into the Bloch-Torrey equation [48] by adding the diffusion contribution given by Eq. 6 where the vectorial density is given by the magnetization vector ${\bf M}$ :

$$\begin{array}{*{2}{l}} {\bf M}_{,t} = & \displaystyle {\bf ... ...t\left(\mathrm{D}\,\nabla {\bf M}\right) \mbox{.} \end{array}$$ (8)

Expression of the Attenuation

The DTMRI measurement is performed for the component of the magnetization transverse to the ${\bf B}_0$ field. Then only the $M^1$ and $M^2$ components are of interest for the signal acquisition. The complex notation is used where $\underline{M}({\bf x},t)=M^1({\bf x},t)+i\, M^2({\bf x},t)$ . Equation 8 is rewritten as:

$$\underline{M}_{,t}=-i\gamma{\bf G}\cdot{\bf x}\;\underline{M}-\frac{\underline{M}}{T_2}+\nabla\cdot\left(\mathrm{D}\nabla\underline{M}\right)$$ (9)

In order to find an analytical solution, it is supposed that for a given position of the molecule, the diffusion tensor $\mathrm{D}$ does not vary with the displacement of the molecule due to the diffusion. It is also supposed that this tensor does not depend on $t$ . Then, for every position ${\bf x}$ , the tensor $\mathrm{D}$ is constant. The attenuation expresses the loss of acquired signal caused by the diffusion. In some conditions, this attenuation can be analytically expressed. The complex magnetization is therefore separated as attenuation caused by the diffusion $A$ , exponential decreasing of intensity of characteristic time $\frac{1}{\alpha}$ because of the relaxation parameter $T_2$ , and a phase component given by the spin phase $\varphi$ . The general expression of the complex magnetization vector is therefore given by:

$$\underline{M}({\bf x},t)=A_{\bf x}(t)\,e^{-\alpha(t)}\,e^{i\varphi({\bf x},t)}$$ (10)

where the expression $A_{\bf x}$ indicates that this attenuation is given for a fixed position as an argument and not as a variable. Plugging this substitution into Eq. 9 gives the expression:

$$\begin{array}{c} \left( i\varphi_{,t}({\bf x},t) -\alpha'(... ...e^{i\varphi({\bf x},t)}\,e^{-\alpha(t)} \mbox{ .} \end{array}$$ (11)

The spin phase $\varphi$ is defined as:

$$\varphi({\bf x},t) = \gamma\,\int_{0}^{t}{\bf x}\cdot{\bf G}(t')\,\mathrm{d}t'$$

and the attenuation $\alpha$ is such that :

$$\alpha(t)=-\frac{t}{T_2}$$

Equation 11 can therefore be rewritten as:

$$A_{\bf x}'(t)=-A_{\bf x}(t)\,\left(\nabla \varphi\right)^T\mathrm{D}\left(\nabla\varphi\right)$$ (12)

where, in the case of the Cartesian coordinate system,

$$\nabla\varphi({\bf x},t)=\gamma\int_{0}^{t}{\bf G}(t')\;\mathrm{d}t'$$ (13)

The matrix $B$ can be defined as the outer product ${\bf B}=\left(\nabla\varphi\right)\left(\nabla\varphi\right)^T$ . Then Eq. 12 can be written using the trace operator $\mathrm{tr}$ (The following equation is true in every coordinate system and will be proved in the Cartesian and curvilinear case in the following parts.):

$$A'_{\bf x}(t)=-\mathrm{tr}\left(\mathrm{B}({\bf x},t)\,\mathrm{D}\right)\;A_{\bf x}(t)$$ (14)

This is a first order differential equation that can be solved analytically to give the attenuation at every position compared to the initial signal:

$$A_{\bf x}(t) = A(0)\,e^{-\int_{0}^{t}\mathrm{tr}\left(\mathrm{B}({\bf x},t')\mathrm{D}\right)\,\mathrm{d}t'}$$

The DTMRI technique can access the value of the attenuation $A_{\bf x}(t)$ ; and the value of $A_{\bf x}(0)$ is approximated by an acquisition performed by a normal MRI technique (where the signal is supposed to be unattenuated by the diffusion). Then the results are usually presented by the relation that links the attenuation to the diffusion tensor:

$$\ln\left(\frac{A_{\bf x}(t)}{A_{\bf x}(0)}\right) = -\int_{0}^{... ...)\right)^T \mathrm{D} \left(\nabla\varphi({\bf x},t)\right) \, \mathrm{d}t'$$ (15)

or with the matrix $B$ , using the linearity of the integral and the fact that $\mathrm{D}$ is considered as a constant through time

$$\ln\left(\frac{A_{\bf x}(t)}{A_{\bf x}(0)}\right) = - \mathrm{... ...int_{0}^{t} \mathrm{B}({\bf x},t') \mathrm{d}t' \right) \mathrm{D} \right]$$ (16)

Diffusion Measurement

Static Case

In order to measure the effect caused by the diffusion in a static organ such as the brain or ex-vivo heart, a special pulse sequence is performed. This pulse sequence is usually the spin echo diffusion sequence introduced by Stejskal and Tanner [49]. After the classical RF pulse, which flips the magnetization vector at   $90^{\circ}$ [46], a constant gradient of diffusion ${\bf G}_d$ is applied during a time $\delta_d$ . During this time, the water molecules are dephased depending on their positions. At time $t=\frac{TE}{2}$ after the RF pulse, a $180^{\circ}$ pulse is applied to invert the phase of the spins. Then a second gradient ${\bf G}_d$ is applied during the same period of time $\delta_d$ . This is the rephasing part. If the molecule is fixed, the spin will be exactly rephased at the end of this second gradient of diffusion. However, if the molecule has moved because of the diffusion in the direction of the gradient, then the phase will not be exactly recovered. At the time $t=TE$ , after the initial RF pulse, and because of the phase inversion, a spin echo occurs and the MRI acquisition is performed. The sequence of the pulse is drawn in Fig. 5.

Fig. 5: The classical spin-echo diffusion weighted sequence. The magnetization vector is flipped by the RF pulse, then the diffusion gradient ${\bf G}_d$ is applied during the time $\delta_d$ . After a time $\frac{TE}{2}$ , the magnetization vector is inverted and a second gradient diffusion is applied during the same amount of time. After a time $TE$ , the spin echo occurs and the signal is acquired using the reading gradients.

The DTMRI technique reconstructs the diffusion tensor by measuring the attenuation of the received signal, whereas, a classical MRI acquisition is done without the diffusion tensor (in this case the attenuation caused by the diffusion is supposed to be negligible). Then a reconstruction process is done to determine the value of the diffusion tensor utilizing knowledge of the attenuation factor at each position.

In the case of a fixed organ defined by the Cartesian coordinates, the relation given in Eq. 13 can be used and the matrix $\mathrm{B}$ is therefore expressed as:

$$\mathrm{B}(t)= \gamma^2\left(\int_0^t{\bf G}_d(t')\,\mathrm{d}t'\right)\otimes\left(\int_0^t{\bf G}_d(t')\,\mathrm{d}t'\right)$$ (17)

where the integral of the gradient can be analytically expressed with knowledge of the pulse sequence. The origin of time is taken when the first gradient of diffusion begins:

$$\begin{array}{l} \displaystyle \int_{0}^{t}{\bf G}_d(t)\;\ma... ...x{if} &\Delta+\delta_d<t \\ \end{array} \right. \end{array}$$ (18)

The signal acquisition is performed at the time $\tau$ , then the matrix $\mathrm{B}$ is expressed by:

$$\int_0^{\tau} \mathrm{B}(t)\,\mathrm{d}t=\gamma^2\,\delta_d^2\left(\Delta-\frac{\delta_d}{3}\right)\;{\bf G}_d\otimes{\bf G}_d$$ (19)

which is proven in appendix A. The measurement of the attenuation gives the relation:

$$\ln\left(\frac{A_{\bf x}(t)}{A_{\bf x}(0)}\right) = - \gamma^2... ...rm{tr} \left[ \left( {\bf G}_d\otimes{\bf G}_d \right) \mathrm{D} \right]$$ (20)

where it is proven in appendix B that $\mathrm{tr}\left(\mathrm{B}\;\mathrm{D}\right)=\gamma^2\delta^2\,{\bf G}_d^T\mathrm{D}{\bf G}_d$ in Cartesian coordinates. In the case where $\Delta\gg\delta_d$ , the equation is usually simplified as:

$$\ln\left(\frac{A_{\bf x}(t)}{A_{\bf x}(0)}\right) = - \gamma^2\,\delta_d^2\Delta \;\; {\bf G}^T_d \mathrm{D} {\bf G}_d$$ (21)

where the time $\Delta$ is called the effective diffusion time and the scalar ${\bf G}^T\,\mathrm{D}{\bf G}$ is called the effective diffusion coefficient that represents the projection of the diffusion tensor onto the direction of the gradient. The knowledge of this projection of the tensor can be visualized by taking the length of the ellipsoid of diffusion in the direction given by the gradient. This information enables reconstruction of the diffusion tensor by recording this projection value in different directions [11,50]. Finally, denoting ${\bf k}=\gamma\delta\,{\bf G}$ , the equation can be given in the form:

$$\ln\left(\frac{A_{\bf x}(t)}{A_{\bf x}(0)}\right) = - \Delta \;\; {\bf k}^T \mathrm{D} {\bf k}$$ (22)

Deforming Media

Curvilinear Coordinates and Tensor Notation

To simplify the following notation, the Einstein summation convention is used such that for a repeated index $a^i\,b_i=\sum_{i=2}^3a^i\,b_i$ . The curvilinear coordinates are introduced to follow the deformations of the material. The fixed Cartesian basis $\left({\bf e}_1,{\bf e}_2,{\bf e}_3\right)$ will serve as the initial frame of reference. Every position in the initial reference frame is given by the vector ${\bf x}=x^i\,{\bf e}_i$ . When the material is deformed, the new position is located using curvilinear coordinates $(\xi^1,\xi^2,\xi^3)$ [51,35] with the contravariant basis $({\bf g}_1,{\bf g}_2,{\bf g}_3)$ . The covariant basis will be denoted $\left({\bf g}^1,{\bf g}^2,{\bf g}^3\right)$ . Every position in the deformed material is given by the vector $\boldsymbol{\xi}=\xi^i\,{\bf g}_i=\xi_i\,{\bf g}^i$ , where $\xi_i$ are the contravariant components of the vector. The basis is not necessarily an orthonormal basis, and therefore ${\bf g}^i\neq{\bf g}_i$ . The curvilinear basis is given by its covariant vectors:

$${\bf g}_i=\frac{\partial x^j}{\partial \xi^i}{\bf e}_{j}$$

or its contravariant vectors:

$${\bf g}^i=\frac{\partial \xi^i}{\partial x^j}{\bf e}^{j}$$

The arc element in the deforming material is given by $\mathrm{d}s=\sqrt{\mathrm{g}_{ij}\,\mathrm{d}x^i\,\mathrm{d}x^j}$ , where the contravariant metric tensor is introduced such that $\mathrm{g}_{ij}={\bf g}_i\cdot{\bf g}_j$ . The covariant metric tensor is defined by $\mathrm{g}^{ij}={\bf g}^i\cdot{\bf g}^j$ and the two metric tensors are related by a matrix inversion such that $[\mathrm{g}_{ij}]=[\mathrm{g}^{ij}]^{-1}$ , where $[ \;\;]$ indicates a matrix.

The special tensor $\delta$ will be used such that $\tensor{\delta}{_i^j}$ and $\tensor{\delta}{^i_j}$ equal one if $i=j$ and zero otherwise.

The tensor $\mathrm{D}$ is expressed in the original coordinate system by $\mathrm{\overline{D}}=\tensor{\mathrm{\overline{D}}}{_i^j}\,{\bf e}^i\otimes{\bf e}_j$ . The tensor expressed in the new coordinate system will be denoted $\mathrm{D}=\tensor{\mathrm{D}}{_i^j}{\bf g}^i\otimes{\bf g}_j$ . The diffusion tensor follows the rules of the transformation of a tensor of rank two:

$$\tensor{\mathrm{D}}{_i^j}({\bf x},t)=\frac{\partial x^j}{\partial... ...ac{\partial \xi^k}{\partial x^i}\,\tensor{\mathrm{D}}{_k^l}(\boldsymbol{\xi},t)$$ (23)

It is shown in appendix C that the transformation of the tensor is given by:

$$\tensor{\mathrm{D}}{_i^j}={\bf g}_i\,\overline{\mathrm{D}}\,{\bf g}^{j}=\mathrm{tr}\left({\bf g}^j\otimes {\bf g}_{i}\,\overline{\mathrm{D}}\right)$$ (24)

The $\nabla$ operator is given by:

$$\nabla={\bf g}^i\frac{\partial \;\;}{\partial \xi^i}={\bf e}^i\frac{\partial \;\;}{\partial x^i}$$ (25)

The Christoffel symbols of second kind $\Gamma$ [52,53] are also defined such that:

$${\bf g}_{i,j}=\tensor{\Gamma}{^k_{ij}}\,{\bf g}_k$$ (26)

and are related to the coordinate by the following relation:

$$\tensor{\Gamma}{^i_{jk}}=\frac{\partial^2 x^l}{\partial \xi^j\partial \xi^k}\,\frac{\partial \xi^i}{\partial x^l}$$ (27)

The derivative of the covariant basis vector will also be used and it is proven in appendix D that:

$${\bf g}^{i}_{,j}=-\tensor{\Gamma}{^i_{kj}}{\bf g}^k$$

The Christoffel symbols can also be calculated by the relation with the metric tensor

$$\tensor{\Gamma}{^k_{ij}} = \frac{1}{2}\,\mathrm{g}^{kl} \left( \mathrm{g}_{li,j} + \mathrm{g}_{lj,i} - \mathrm{g}_{ij,l} \right)$$ (28)

Some relations about tensors reminded here will be used later. ${\bf a}$ , ${\bf b}$ , ${\bf c}$ and ${\bf d}$ designate any vectors.

$$\left({\bf a}\otimes{\bf b}\right)\cdot{\bf c}=\left({\bf b}\cdot{\bf c}\right)\,{\bf a}$$ (29)

$$\left({\bf a}\otimes{\bf b}\right)\,\left({\bf c}\otimes{\bf d}\right)=\left({\bf b}\cdot{\bf c}\right)\,\left({\bf a}\otimes{\bf d}\right)$$ (30)

Mechanical Deformation

The media in which the diffusion occurs is now supposed to be deformable through time. In the following section, the implicit Einstein summation convention will be used such that $a^{i}\,b_{i}=\sum_{i}a^{i}\,b_{i}$ . To study the mechanical deformation, the referential (or Lagrangian) approach is used [54,51,55].

The fixed original referential is given by the basis $({\bf e}_{1},{\bf e}_{2},{\bf e}_{3})$ , and every initial position for this basis is given by ${\bf X}=X^{i}\,{\bf e}_{i}$ . The material can then be deformed; the curvilinear basis associated with the deformation was introduced in the previous section. The curvilinear basis is related to the initial coordinate by the vectorial function $\boldsymbol{\xi}=\boldsymbol{\xi}({\bf X},t)$ .

Fig. 6: Example of deformation, where the image on the left represents the undeformed initial referential and the image on the right is the deformed referential where the curvilinear basis depends on the position.

The deformation is expressed with the tensorial deformation gradient such that $\mathrm{d}\boldsymbol{\xi}=\mathrm{F}\,\mathrm{d}{\bf X}$ , where

$$\mathrm{F}=\frac{\partial \boldsymbol{\xi}}{\partial {\bf X}}=\frac{\partial \xi^i}{\partial X^j}{\bf e}^{j}\otimes{\bf g}_i$$

DTMRI Measurement

The measurement of the diffusion tensor in the case of the deforming heart has been performed by Reese et al. [30,31,32,33]. For this purpose, a special pulse sequence is performed. After the RF flip angle at the time $t=t_0$ , the initial diffusion gradient is turned on during the time $\delta_d$ . After a time $\frac{TE}{2}$ , a second RF pulse is performed, but only at $90^{\circ}$ . The cardiac cycle is trigged, and the sequence waits for the next R-wave that occurs at the time $\Delta$ . After a time $t=\Delta+t_0$ , a third RF pulse at $90^{\circ}$ is performed and the second diffusion gradient is applied. This diffusion gradient ${\bf G}_d$ is identical to the first one and occurs during the same configuration of the heart. After the time $\Delta+t_0+\frac{TE}{2}=\tau$ , the stimulated echo occurs and the acquisition is performed. The sequence is represented in Fig. 7

Fig. 7: The Pulse sequence used for the dynamic acquisition. The cardiac cycle is trigged on the Electrocardiogram (ECG). After the first R-wave, an RF pulse is used to flip the magnetization to $90^{\circ}$ . The diffusion gradient ${\bf G}_d$ is then applied during the time $\delta_d$ and a second RF flip is applied after the first one. The second R-wave is trigged and a third magnetization flip is applied. A second gradient of diffusion is performed (identical to the first one) and the signal acquisition is performed when the spin echo occurs at $t=\tau$ .

Equations 15 and 16 are still true (the proof in the curvilinear coordinate of Eq. 16 is exactly the same as in the Cartesian coordinates) but the spin phase is now defined in the deformed material $\varphi(\boldsymbol{\xi},t)=\gamma\,\int_{0}^{t}{\bf X}(\boldsymbol{\xi},t')\cdot{\bf G}(t')\,\mathrm{d}t'$ . Considering a smooth deformation, the differential relation for the spin phase is therefore:

$$\nabla^{T}\varphi({\boldsymbol{\xi}},t)\;\mathrm{d}{\boldsymbol{\xi}}=\nabla^{T}\varphi({\bf X},t)\;\mathrm{d}{\bf X}$$

Using the fact that the deformation gradient is defined by $\mathrm{d}\boldsymbol{\xi}=\mathrm{F}\;\mathrm{d}{\bf X}$ , the following relation links the gradient to the spin phase:

$$\nabla\varphi(\boldsymbol{\xi},t)=\mathrm{F}^{-T}({\bf X},t)\;\nabla\varphi({\bf X},t)$$ (31)

where $\nabla\varphi({\bf X},t)=\gamma\int_{0}^{t}{\bf G}(t')\,\mathrm{d}t'$ . The tensor of diffusion is also expressed in the deformed basis, denoting $\overline{\mathrm{D}}$ the tensor expressed in the initial referential[52]:

$$\mathrm{D}^{i}_{j}=\frac{\partial \xi^i}{\partial X^k}\;\frac{\pa... ...}^{k}_{l}=F^{i}_{k}\,\overline{\mathrm{D}}^{k}_{l}\,\left(F^{-1}\right)^{l}_{j}$$

that can be expressed in the matrix form:

$$\mathrm{D}=\mathrm{F}\,\overline{\mathrm{D}}\,\mathrm{F}^{-1}$$ (32)

Denoting $\nabla\varphi({\bf X},t)=\nabla \varphi$ , for the sake of simplicity, the complete measured attenuation is therefore given by applying Eq. 32 and 31 to Eq. 15:

$$\ln\left(\frac{A(t)}{A(0)}\right) = -\int_{0}^{t}\left(\nabla \va... ...T\;\overline{D}\;\mathrm{F}^{-1}\,\mathrm{F}^{-T}\;\nabla \varphi\;\mathrm{d}t'$$

The deformation tensor can always be decomposed between stretch and rotation such that $\mathrm{F}=\mathrm{R}\mathrm{U}$ [54], where $\mathrm{R}$ is a rotation matrix such that $\mathrm{R}^{-1}=\mathrm{R}^{T}$ and $\mathrm{U}$ is the right stretch tensor, which is a symmetric matrix. Then the following relation holds: $\mathrm{U}=(\mathrm{F}^{T}\,\mathrm{F})^{\frac{1}{2}}$ . The attenuation can therefore be decomposed such that:

$$\ln\left(\frac{A(t)}{A(0)}\right) = -\int_{0}^{t}\left(\nabla \va... ...{T}\;\overline{D}\;\left(\mathrm{U}^{-1}\right)^2\;\nabla \varphi\;\mathrm{d}t'$$ (33)

If it is assumed that only a linear deformation affects the diffusion tensor $\overline{D}$ from the referential system via rotation during the time of the acquisition[56], then the cardiac deformation is assumed to be treated locally as a rigid body. The attenuation given by Eq. 15 is:

$$\ln\left(\frac{A(t)}{A(0)}\right) = -\int_{0}^{t}\left(\nabla \va... ...m{R}\;\overline{D}\;\mathrm{R}^{T}\mathrm{F}^{-T}\;\nabla \varphi\;\mathrm{d}t'$$ (34)

and using $\mathrm{R}^{-1}=\mathrm{R}^{T}$ and $\mathrm{U}=\mathrm{U}^{T}$

$$\ln\left(\frac{A(t)}{A(0)}\right) = -\int_{0}^{t}\left(\nabla \va... ...}\;\mathrm{U}^{-1}\;\overline{D}\;\mathrm{U}^{-1}\;\nabla \varphi\;\mathrm{d}t'$$ (35)

The integral can be simplified with the knowledge of the pulse sequence used:

$$\begin{array}{*{1}{l}} \ln\left(\frac{A(\tau)}{A(0)}\right)=... ...}^{-1}\,{\bf G}\;\delta^2\,\mathrm{d}t' \mbox{ .} \end{array}$$

It is assumed that $\delta\ll\Delta$ where $\Delta$ covers a complete cardiac cycle close to $1\;s$ for a human and $\delta\simeq 6-8\;ms$ . Then the first part of the integral is negligible compared to the second part. Denoting ${\bf k}=\gamma\,\delta\;{\bf G}$ , the attenuation can therefore be approximated as:

$$\ln\left(\frac{A(\tau)}{A(0)}\right)=- {\bf k}^{T} \left( \int... ...,\mathrm{U}^{-1}\;\overline{D}\;\mathrm{U}^{-1}\,\mathrm{d}t' \right) {\bf k}$$ (36)

This equation can be rewritten in the classical form as in the static case by comparison with Eq. 22:

$$\ln\left(\frac{A(\tau)}{A(0)}\right)=- \Delta \;\; {\bf k}^{T} \mathrm{D}_{\mbox{\scriptsize {obs}}} {\bf k} \mbox{ ,}$$ (37)

where the observed diffusion $\mathrm{D}_{\mbox{\scriptsize obs}}$ is changed due to the deformation such that:

$$\mathrm{D}_{\mbox{\scriptsize {obs}}}({\bf X}) = \frac{1}{\Delta}... ...erline{\mathrm{D}}({\bf X})\;\mathrm{U}^{-1}({\bf X},t)\;\mathrm{d}t \mbox{ .}$$ (38)

Finally, the strain tensor can be introduced as defined by:

$$S=\frac{1}{2}\left(\mathrm{F}^T\mathrm{F}-\mathrm{I}\right)$$

where $\mathrm{F}^T\mathrm{F}=\mathrm{U}^2$ . If the strain is small, a first order linear approximation can be used:

$$\mathrm{U}\simeq \mathrm{I}-\mathrm{S}$$

Then, the first order approximation of the observed diffusion is given by:

$$\mathrm{D}_{\mbox{\scriptsize {obs}}}({\bf X}) = \frac{1}{\Delta}... ...overline{\mathrm{D}}\,\mathrm{S}({\bf X},t)\; \right) \mathrm{d}t \mbox{ .}$$ (39)

Expression of the MR Signal in Curvilinear Coordinates

In order to enable the simulation of the MR signal in a deforming media, the Bloch-Torrey equation is expressed in curvilinear coordinates that have been introduced in the previous section.

Cartesian Coordinates

In Cartesian coordinates, using the notation of Eq. 8, the three dimensional set of equations is given by:

$$\left\{ \begin{array}{l} \begin{array}{ll} M_{,t}^{x}=&\gamma\... ...nabla\left(\mathrm{D}\nabla M^{z}\right) \\ \end{array} \end{array} \right.$$ (40)

where $M^x$ , $M^y$ , and $M^z$ depend on the Cartesian components ${\bf x}$ . The part that depends on the coordinate system is the diffusive part. A general scalar $\phi$ is introduced to express any component of the magnetization vector. The expression $\nabla \cdot \left(\mathrm{D}\,\nabla\,\phi\right)$ is then expressed in the new coordinate system. In the case of the Cartesian coordinates, it is shown in appendix E, that the diffusion part can be expressed as:

$$\nabla\cdot\left(\mathrm{D}\nabla \phi\right) = \tensor{\mathrm{D}}{_i^j_{,i}}\,\phi_{,j}+\tensor{\mathrm{D}}{_i^j}\phi_{,ij}$$ (41)

This equation can be inserted in Eq. 40 where $\phi$ is respectively replaced by $M^x$ , $M^{y}$ , and $M^z$ in the three scalar equations.

General Curvilinear Coordinates

The curvilinear coordinates designated by the vector $\boldsymbol{\xi}=\xi^i{\bf g}_i$ are used. The tensor $\mathrm{D}$ is expressed in the new coordinate system by the relation given in Eq. 24. The density flux ${\bf j}=\mathrm{D}\nabla\phi$ can be expressed as in the case of the Cartesian coordinates (see beginning of appendix F).

$$\mathrm{D}\nabla\phi=\tensor{\mathrm{D}}{_i^j}\phi_{,j}{\bf g}^i$$ (42)

where $\phi_{,j}$ denotes $\frac{\partial \phi}{\partial \xi^j}$ . Finally, the diffusion is given by the divergence of this flux and gives the expression proven in appendix F:

$$\nabla\cdot\left(\mathrm{D}\nabla\phi\right) = \left[ \left(\t... ...{\Gamma}{^l_{ji}}\tensor{\mathrm{D}}{_l^k}\,\phi_{,k} \right] \mathrm{g}^{ij}$$ (43)

The complete Bloch-Torrey equation where the magnetization components are expressed on the general curvilinear coordinates is therefore given by the vectorial equation:

$$\begin{array}{*{2}{l}} \displaystyle {\bf M}_{,t}=& \displ... ...\bf M}_{,k} \right] \,\mathrm{g}^{ij} \mbox{ ,} \end{array}$$ (44)

where ${\bf M}$ is expressed as ${\bf M}(\xi^1,\xi^2,\xi^3)=M^{1}(\boldsymbol{\xi}){\bf e}_1+M^{2}(\boldsymbol{\xi}){\bf e}_2+M^{3}(\boldsymbol{\xi}){\bf e}_3$ .

Application in Prolate Spheroidal Coordinates

The prolate spheroidal coordinates [35] are defined by the change of variable:

$$\left\{ \begin{array}{*{2}{l}} X^1&=C\,\sinh(\xi^1)\sin(\xi^2)\... ...xi^2)\sin(\xi^3) \\ X^3&=C\,\cosh(\xi^1)\cos(\xi^2) \\ \end{array} \right.$$ (45)

The coordinates can easily model the geometry of the ventricle, as seen in fig. 8

Fig. 8: Figure of the geometry of the prolate spheroidal coordinates.

The contravariant basis vectors are defined as

$${\bf g}_1= C \left( \begin{array}{c} \cosh(\xi^1)\sin(\xi^2)\... ...1)\sin(\xi^2)\sin(\xi^3) \\ \sinh(\xi^1)\cos(\xi^2) \\ \end{array} \right)$$

$${\bf g}_2= C \left( \begin{array}{c} \sinh(\xi^1)\cos(\xi^2)\... ...)\cos(\xi^2)\sin(\xi^3) \\ -\cosh(\xi^1)\sin(\xi^2) \\ \end{array} \right)$$

$${\bf g}_3= C \left( \begin{array}{c} \sinh(\xi^1)\sin(\xi^2)\... ...\xi^3) \\ -\sinh(\xi^1)\sin(\xi^2)\cos(\xi^3) \\ 0\\ \end{array} \right)$$

Then, the contravariant metric tensor is given by: $g_{ij}=g_i \cdot g_j$ . It can be noticed that the prolate spheroidal basis is orthogonal, meaning that:

$$\forall (i,j)\in[\![1;3]\!]\;\;,\;i\neq j\Rightarrow {\bf g}_i\cdot {\bf g}_j=0$$

Therefore, only the diagonal entries of the metric tensor are not zero:

$$\left\{ \begin{array}{rl} \mathrm{g}_{11}=\mathrm{g}_{22}=&C^2\... ...\xi^2)\right)\\ g_{33}=&C^2\sinh^2(\xi^1)\sinh^2(\xi^2) \end{array} \right.$$

The covariant metric tensor is obtained by inverting the matrix of the contravariant metric tensor. Since it is diagonal, the inversion is obvious:

$$\left\{ \begin{array}{rl} \mathrm{g}^{11}=\mathrm{g}^{22}=&\fra... ...ht)}\\ g^{33}=&\frac{1}{C^2\sinh^2(\xi^1)\sin^2(\xi^2)} \end{array} \right.$$

In this case of orthogonal coordinates, the scalar factors $h$ are defined such that $h_i=\sqrt{g^{ii}}$ (no summation convention). Then the change of volume in the new coordinate system is given by:

$$\mathrm{d}\xi^1\mathrm{d}\xi^2\mathrm{d}\xi^3 = \left(\prod_i h_i\right) \mathrm{d}X^1\mathrm{d}X^2\mathrm{d}X^3$$

$$\begin{array}{*{2}{l}} \Rightarrow& \mathrm{d}X^1\mathrm{d}... ...\ & \mathrm{d}\xi^i\mathrm{d}\xi^2\mathrm{d}\xi^3 \end{array}$$

The covariant basis vectors given by $g^i=g^{ij}g_j$ are easily calculated in the orthogonal case by $g^i=g^{ii}g_{i}$ (no summation convention). Finally, the Christoffel symbols are expressed using Eq. 28. The non zero components are given by:

$$\left\{ \begin{array}{*{7}{l}} \tensor{\Gamma}{^1_{11}}=&& &n&\... ...\tensor{\Gamma}{^3_{32}}=&& & &\mathrm{cotan}(\xi^2)& \\ \end{array} \right.$$ (46)

where

$$n=\frac{1}{\sinh^2(\xi^1)+\sin^2(\xi^2)}$$

Finally, the equation can be rewritten taking into account that the metric tensor is diagonal. Writing the summation symbols, the equation for an orthogonal curvilinear system of Eq. 44 is therefore:

$$\begin{array}{*{2}{l}} {\bf M}_{,t}= & \displaystyle {\bf... ...{,k} \right) \Bigr) \biggr] \Biggr) \mbox{ .} \end{array}$$ (47)

Numerical Solution

Importance of the Implicit Method

The one-dimensional scalar equation with constant diffusion coefficient in the Cartesian case is given in Eq. 4. In order to use the finite differences methods, the numerical differentiation operator $\mathscr{D}$ is introduced [57,58]:

$$\mathscr{D}^1 \phi=\phi(x+\Delta\,x)-\phi(x-\Delta\,x)$$

where $x$ is the one dimensional variable. This operator is the central difference operator and is second order accurate provided that

$$\exists A \in\mathbb{R}^{+}\;,\; \left\vert\frac{\mathscr{D}^1 \p... ...elta x}-\frac{\partial \phi}{\partial x}\right\vert\leq A\,\vert\Delta x\vert^2$$

The second derivative in space is also a second order centered scheme:

$$\mathscr{D}^{11} \phi=\phi(x+\Delta\,x)-2\,\phi(x)+\phi(x-\Delta\,x)$$

where $\phi_{,11}\simeq\frac{\mathscr{D}^{11}\phi}{\left(\Delta x\right)^2}$ . The temporal derivative is approximated by the first order forward scheme [59]:

$$\mathscr{D}^t \phi = \phi(t+\Delta t)-\phi(t)$$ (48)

The function is discretized on a spatial grid of intervals of size $\Delta x$ . Every position $x=k\,\Delta x$ is denoted to simplify by $k$ . The grid is divided in $N$ intervals where the index starts at zero. The function is stored as a vector ${\bf u}$ such that $\phi(k\Delta x)=u[k]$ . The matrix operator that acts on the vector directly is introduced by using the same notation. The matrix $\mathscr{D}^1$ is for instance described by:

$$[\mathscr{D}^1]_{i,j}= \left\{ \begin{array}{rl} 1 & \mbox{if ... ... \\ -1 & \mbox{if }\;\; i=j+1 \\ 0 & \mbox{otherwise} \end{array} \right.$$

and

$$[\mathscr{D}^{11}]_{i,j}= \left\{ \begin{array}{rl} 1 & \mbox{... ...j \\ 1 & \mbox{if }\;\; i=j+1 \\ 0 & \mbox{otherwise} \end{array} \right.$$

Both can be easily computed by using the convolution mask such that $\mathscr{D}^1=\mathrm{I}_{N}\ast[-1,0,1]$ , $\mathscr{D}^{11}=\mathrm{I}_N\ast[1,-2,1]$ where $\mathrm{I}_N$ is the identity matrix of size $N$ , and the result of the convolution is cropped to the size $N\times N$ . Therefore, the matrix product $\frac{1}{2\,\Delta x}\mathscr{D}^1\,{\bf u}$ is an approximation of $\phi_{,1}$ , respectively for $\frac{1}{\left(\Delta x\right)^2}\mathscr{D}^{11}{\bf u}$ and $\phi_{,11}$ .

The one dimensional scalar diffusion equation given by:

$$\frac{\partial \phi}{\partial t}=D\,\triangle \phi$$

is rewritten in matrix form as:

$$\frac{1}{\Delta t}\mathscr{D}^t{\bf u}=\frac{D}{\left(\Delta x^2\right)}\,\mathscr{D}^{11}{\bf u}$$ (49)

To differentiate the unknown vector at the time $t+\Delta t$ and the previously known value at the time $t$ , the vector ${\bf u}(t+\Delta t)$ is denoted ${\bf u}^{+}$ . Then the first scheme possible to solve this equation is the first explicit Euler method, which states that:

$$\frac{{\bf u}^{+}-{\bf u}}{\Delta t} = \frac{D}{\left(\Delta x\right)^2}\mathscr{D}^{11}{\bf u}$$

and gives the numerical implementable recursive algorithm:

$${\bf u}^{+} = \left( \mathrm{I}_N + D \frac{\Delta t}{\left(\Delta x\right)^2} \mathscr{D}^{11} \right) {\bf u}$$

where knowledge of the initial condition ${\bf u}(t=0)={\bf u}_0$ enables the calculation of the function through time. This explicit method is easily implemented, however it can be shown [58,60,61] that it is only conditionally stable since the numerical solution tends to grow faster than the real solution, and can become unstable. It is shown that the ratio $\mu=2\,D\frac{\Delta t}{\left(\Delta x\right)^2}$ has to be such that $\mu\leq 1$ to have a stable scheme. In this condition, the time step that can be used for the evolution decreases as the square of the number of points on the grid. In the case of a variable diffusion coefficient in time and space, the equation of diffusion is then expressed as:

$$\phi_{,t}=D({\bf x},t)\,\phi_{,11}+D_{,1}({\bf x},t)\phi_{,1}$$

In this case, a stability condition also occurs on the derivative of the diffusion coefficient [58]. Calling $\mu=\frac{\Delta t}{\left(\Delta x\right)^2}D$ and $\nu=\frac{\Delta t}{\Delta x}D_{,1}({\bf x},t)$ . The condition of stability is then given by:

$$\left\{ \begin{array}{*{2}{l}} \vert\nu\vert&\leq\,2\mu \\ 2\,\mu&\leq1 \end{array} \right.$$

which implies:

$$\left\{ \begin{array}{*{2}{l}} \Delta x&\leq 2\frac{D}{\vert D_... ...} \\ \Delta t&\leq \frac{\left(\Delta x\right)^2}{2\,D} \end{array} \right.$$

If the derivative of the diffusion is large compared to the gradient itself (noise in the measurement of the diffusion tensor), the spatial resolution will increase significantly and the spatial time resolution will become too small to be computed. If this explicit method can still be applied in the one dimensional case, the two dimensional or the three dimensional cases become too complex to solve because the time step is too small and therefore requires too much computation (the condition of stability acts also on $x$ , $y$ , and $z$ ). A numerical example is shown in Fig 9.

Fig. 9: Example of comparison between explicit and implicit method for the diffusion in one dimension. The picture of the left represents an initial condition of $\phi$ represented by a step function (slice selection of the magnetization). For the study, $60$ numerical sampled are taken in account for $\Delta x=0.0167$ and $\Delta t=0.00149$ . The diffusion coefficient is constant set to $D=0.1$ . The picture of the right shows the results after 10 iterations. The red dots represents the numerical values obtained by the implicit method, the blue cross are obtained by the explicit method. The black dotted line represents the true solution obtained by Eq. 5.

Another method of solving this equation is given by the first order implicit Euler method [59], coming back to Eq. 49; this method is given by:

$$\frac{1}{\Delta t}\mathscr{D}^t{\bf u}=\frac{D}{\left(\Delta x^2\right)}\,\mathscr{D}^{11}{\bf u}^{+}$$

and give the following scheme:

$$\left( \mathrm{I}_N - D \frac{\Delta t}{\left(\Delta x\right)^2} \mathscr{D}^{11} \right) {\bf u}^{+} = {\bf u}$$

which can be seen as a matrix equation $\mathcal{A}{\bf u}^{+}={\bf b}$ , where $\mathcal{A}=\left(\mathrm{I}_N -D\frac{\Delta t}{\left(\Delta x\right)^2}\mathscr{D}^{11}\right)$ and ${\bf b}={\bf u}$ . Then the solution at the following time step is given by the solution of the linear system ${\bf u}^+=\mathcal{A}^{-1}{\bf b}$ . This complete implicit method is unconditionally stable and gives good results for the diffusion effect. Contrary to the explicit case, a slight decrease of the numerical solution occurs during the iterative process, then the numerical solution of a constant function decreases to zero. This is characteristic of a numerical diffusive effect added by this numerical method. This numerical diffusion doesn't cause a problem when the solution is supposed to be diffused; however, for the case of the Bloch-Torrey equation, the rotational term given by vectorial product needs to be solved without artificial decrease. Without diffusion, the magnitude of the solution should remain constant. In order to get a more accurate method, the Crank-Nichols scheme was used. This method states that [58,60]:

$$\frac{1}{\Delta t}\mathscr{D}^t{\bf u}=\frac{D}{\left(\Delta x^2\right)}\,\mathscr{D}^{11}\frac{{\bf u}^{+}+{\bf u}}{2}$$

giving the linear system:

$$\left( \mathrm{I}_N - D \frac{\Delta t}{2\,\left(\Delta x\rig... ...athrm{I}_N + D \frac{\Delta t}{2\,\left(\Delta x\right)^2} \right) {\bf u}$$

which can still be written as a matrix system: $\mathcal{A}{\bf u}^+={\bf b}$ , where $\mathcal{A}=\left(\mathrm{I}_N -D\frac{\Delta t}{\left(2\,\Delta x\right)^2}\mathscr{D}^{11}\right)$ , and ${\bf b}=\left(\mathrm{I}_N +D\frac{\Delta t}{\left(2\,\Delta x\right)^2}\mathscr{D}^{11}\right)$ u. Due to its symmetry, this method enables the problem to be solved without artificial decrease or increase of magnitude when the real solution has a constant magnitude.

Application to the Bloch-Torrey Equation

In order to apply the method to the Bloch-Torrey equation, the matrix operators of finite differences have to be defined to act on the three dimensional vector. Now the first order center difference is given by:

$$\mathscr{D}^i {\bf M}={\bf M}(\xi^i+\Delta\,\xi^i)-{\bf M}(\xi^i-\Delta\,\xi^i)$$

where ${\bf M}$ depends on $\boldsymbol{\xi}$ , but only the parameter on which the operator acts is reminded. The order is then given by the relation:

$$\exists A \in\mathbb{R}^{+}\;,\; \left\Vert\frac{\mathscr{D}^i {\... ...ac{\partial {\bf M}}{\partial \xi^i}\right\Vert\leq A\,\vert\Delta \xi^i\vert^2$$

The double second derivative is still given by:

$$\mathscr{D}^{ii} {\bf M}={\bf M}(\xi^i+\Delta\,\xi^i)-2\,{\bf M}(\xi^i)+{\bf M}(\xi^i-\Delta\,\xi^i)$$

but the mixed derivative also needs to be defined:

$$\begin{array}{*{2}{l}} \mathscr{D}^{ij} {\bf M} =& \Bigl( ... ...-\Delta\xi^i,\xi^j-\Delta\xi^j) \Bigr) \mbox{ ,} \end{array}$$

where $i\neq j$ if the boundary condition is given such that ${\bf u}$ is zero outside of the area covered by the limits of the discretization. The function ${\bf M}$ is also discretized on a grid of size $N_1\times N_2\times N_3$ where each voxel has a volume given by $\Delta\xi^1\times\Delta\xi^2\times\Delta\xi^3$ . The function ${\bf M}$ is stored on a large vector ${\bf u}$ of size $3 N_1 N_2 N_3$ such that:

$$\begin{array}{l} M^i(k_1\Delta\xi^1,k_2\Delta\xi^2,k_3\Delta... ...1+N_1\left(k_2+N_2\,k_3\right)\right)] \mbox{ .} \end{array}$$

The mixed derivative matrix operator is a large square matrix of size $\left(3 N_1 N_2 N_3\right)^2$ . Each line $I$ and column $J$ of this matrix refers to the parameters $(i^l,k^l_1,k^l_2,k^l_3)$ and $(i^c,k^c_1,k^c_2,k^c_3)$ for the lines and columns respectively. The matrices are sparse and $\mathscr{D}^{12}$ is, for instance, given by:

$$\begin{array}{l} [\mathscr{D}^{12}]_{I,J}= \\ \left\{ \b... ...\\ 0 & \mbox{ otherwise} \end{array} \right. \end{array}$$ (50)

In order to express the system with matrix operators, the matrix $\mathrm{G}$ is constructed such that $\mathrm{G}{\bf u}$ corresponds to ${\bf M}\times\gamma{\bf B}$ :

$$\begin{array}{l} [\mathrm{G}]_{I,J}= \\ \left\{ \begin{a... ...\\ 0 & & \mbox{ otherwise} \end{array} \right. \end{array}$$ (51)

The projection matrices on the ${\bf e}_i$ axis are also defined to take into account the parameters $T_1$ and $T_2$ . $\mathrm{P}^1$ , $\mathrm{P}^2$ , and $\mathrm{P}^3$ define the projection matrices onto ${\bf e}_1$ , ${\bf e}_2$ , and ${\bf e}_3$ , respectively. The matrices are defined as:

$$\begin{array}{l} [\mathrm{P}^m]_{I,J}= \\ \left\{ \begin... ...\\ 0 & & \mbox{ otherwise} \end{array} \right. \end{array}$$ (52)

The complete system can then be expressed as:

$$\begin{array}{*{2}{l}} \mathscr{D}^t\,{\bf u} = & \Bigl( ... ...igr) {\bf u} -\frac{{\bf u}^z_0}{T_1} \mbox{ ,} \end{array}$$ (53)

where ${\bf u}_0^z$ is the vector that corresponds to $M^z(\boldsymbol{\xi},t=0)$ and the matrix identity $\mathrm{I}$ is introduced to respect the dimensions of the matrices for the addition and substraction. The equation is written in matrix form as:

$$\mathscr{D}^t{\bf u}= \mathcal{S}{\bf u}+{\bf s}$$ (54)

where:

$$\begin{array}{*{2}{l}} \mathcal{S}=& \mathrm{G}-\frac{\math... ...athrm{D}}{_l^k}\mathscr{D}^{k} \right] \mbox{ ,} \end{array}$$

and ${\bf s}=\frac{{\bf u}_0^z}{T_1}$ . The equation is parabolic, therefore, the semi implicit Crank-Nicholson method introduced in the previous section is used such that:

$$\frac{{\bf u}^{+}-{\bf u}}{\Delta t} = \mathcal{S}\frac{{\bf u}^{+}+{\bf u}}{2}+{\bf s}$$ (55)

Reorganizing the terms, the expression is given by the linear system:

$$\left(\mathrm{I}-\Delta t\,\frac{\mathcal{S}}{2}\right){\bf u}^{+} = \left(\mathrm{I}+\Delta t\,\frac{\mathcal{S}}{2}\right){\bf u}+{\bf s}$$ (56)

which is a matrix system of type $\mathcal{A}\,{\bf u}^+={\bf b}$ , where the matrix has dimension $(3\,N_x\,N_y\,N_z)^2$ with:

$$\mathcal{A}=\left(\mathrm{I}-\Delta t\,\frac{\mathcal{S}}{2}\right)$$

and

$${\bf b}=\left(\mathrm{I}+\Delta t\,\frac{\mathcal{S}}{2}\right){\bf u}+{\bf s}$$

Conclusions and Future Works

An overview of the diffusion process and the current methods of measurement by DTMRI have been presented concerning the static and dynamic case. The dynamic method is still being researched and the results need to be correlated and repeated by other measurements.

In order to understand the behavior of the MR signal on a deforming media such as the heart, the equation has been transformed into general curvilinear coordinate system that can fit to the geometry and deformation of the cardiac system. An example of geometrical fitting was given with the application of the prolate spheroidal system. Motion can be introduced first by expanding these coordinate by evolving the value of $C$ as a parameter of the coordinate system through time. The equation expanded in the new coordinate system is linear and the diffusion part makes this equation parabolic.

The numerical implementation of the solution of a diffusion equation has been presented and the importance of the implicit solution, especially for the multidimensional solution, has been iterated. The system of equations were discretized and transformed in a linear matrix system. The solution is given by iterations following the Crank-Nicholson scheme. The rapidity of the results are however largely dependent on the size of the matrix. These matrices are very large and sparse. Different methods to solve the system have to be tested as the direct matrix inversion requires $O(N^3)$ calculations [61] where $N$ is the size of the matrix, which is given by $3\,N_1 N_2 N_3$ and where each $N_i$ gives the size in each dimension $\xi^1$ , $\xi^2$ , and $\xi^3$ . Least squares and conjugate gradient methods enable the process to be sped up; however, for large systems where $N_i \geq 50$ , the iteration process becomes very slow. More research on implementing an Alternative Direction Implicit (ADI) scheme should be performed. Then each iteration should follow three passes but the matrix on the left hand side of the equation becomes tridiagonal and could be rapidly inverted using the Thomas algorithm [58,60].

Acknowledgement

This work was supported in part by the Director, Office of Science, Office of Biological and Environmental Research, Medical Science Division of the U.S. Department of Energy (DEO) under Contract No. DE-AC02-05CH11231 and in part by Public Health Service (NIH) grant number R01 EB000121 awarded by the National Institute of Biomedical Imaging and Bioengineering, Department of Health and Human Services.

Appendix

Proof of Eq. 19

This proof shows that with the expression of the gradient given in Eq. 18, the integral $\int_{0}^{\tau} \mathrm{B}(t) \mathrm{d}t$ at the time of the acquisition, is given by Eq 19 where the tensor $\mathrm{B}$ is given by:

$$B(t)=\gamma^2\left(\int_{0}^{t}{\bf G}_d(t')\;\mathrm{d}t'\right)\otimes\left(\int_{0}^{t}{\bf G}_d(t')\;\mathrm{d}t'\right)$$

Using the linearity of the integral:

$$\begin{array}{*{2}{l}} \displaystyle \int_{0}^{\tau} \mathr... ...ta_d}^{\tau} \mathrm{B}(t)\,\mathrm{d}t \mbox{ .} \end{array}$$

Introducing Eq. 18 in the expression gives:

$$\begin{array}{l} \int_{0}^{\tau} \mathrm{B}(t)\,\mathrm{d}t=... ...\,\mathrm{d}t\right)\\ {\bf G}_d\otimes{\bf G}_d \end{array}$$

$$\Rightarrow \int_{0}^{\tau} \mathrm{B}(t)\,\mathrm{d}t=\\ \gam... ...t(\Delta-\delta_d\right)+\frac{\delta_d^3}{3}\right) {\bf G}_d\otimes{\bf G}_d$$

This giving the expected result: $\int_{0}^{\tau} \mathrm{B}(t)\,\mathrm{d}t=\gamma^2\delta_d^2\left(\Delta-\frac{\delta_d^2}{3}\right)\,{\bf G}_d\otimes{\bf G}_d$ .

Proof of Eq. 14

The proof of Eq. 14 in the case of fixed Cartesian coordinates show that:

$$\left(\nabla \varphi\right)^T\,\mathrm{D}\left(\nabla\varphi\right) = \mathrm{tr}\left(\mathrm{B}\,\mathrm{D}\right)$$

with

$$\mathrm{B}(t) = \left(\nabla \varphi\right)\left(\nabla \varphi\right)^T$$

It is also known that in this case of fixed organs, $\left(\nabla \varphi\right)({\bf x},t)=\gamma\int_{0}^{t}{\bf G}_d(t')\mathrm{d}t'$ . To simplify the notation the vector ${\bf v}$ is introduced such that ${\bf v}=\nabla\varphi$ . Then $\mathrm{B}\,\mathrm{D}=\sum_{i=1}^3\tensor{\mathrm{B}}{_i^k}\tensor{\mathrm{D}}{_k^j}$ and $\tensor{\mathrm{B}}{_i^k}=v_i\,v^k$ . The matrix product gives

$$\tensor{[\mathrm{B}\mathrm{D}]}{_i^j}=\sum_{k=1}^{3}v_i\,v^k\,\tensor{\mathrm{D}}{_k^j}$$

$$\Rightarrow \mathrm{tr}\left(\mathrm{B}\mathrm{D}\right) = \sum_{i=1}^{3} \sum_{k=1}^{3} v^k\,\tensor{\mathrm{D}}{_k^i}\,v_i$$

And on the other hand, the matrix product:

$${\bf v}^T\mathrm{D}{\bf v} = \sum_{j=1}^{3} \left[ \left( \sum_{k=1}^{3} v^k\,\tensor{\mathrm{D}}{_k^j}\right)v_j \right]$$

$$\Rightarrow {\bf v}^T\mathrm{D}{\bf v} = \sum_{i=1}^{3} \sum_{k=1}^{3} v^k\,\tensor{\mathrm{D}}{_k^i}v_i$$

The two expressions are then equal, and using the fact that: $\nabla\varphi=\gamma\int_{0}^t{\bf G}_d(t')\mathrm{d}t'$ , it follows that

$$\mathrm{tr}\left(\mathrm{B}\mathrm{D}\right)=\gamma^2\delta^2\,\mathrm{tr}\left({\bf G}_d{\bf G}_d^T\right)$$

$$\mathrm{tr}\left(\mathrm{B}\mathrm{D}\right)=\gamma^2\delta^2\,{\bf G}_d^T\mathrm{D}{\bf G}_d$$

Proof of Eq. 24

It is shown that the transformation of the tensor of diffusion follows the rule given by Eq. 24:

$$\tensor{\mathrm{D}}{_i^j}={\bf g}_i\,\overline{\mathrm{D}}\,{\bf g}^{j}=\mathrm{tr}\left({\bf g}^j\otimes {\bf g}_{i}\,\overline{\mathrm{D}}\right)$$

The first equality is developped to give:

$${\bf g}_i\overline{\mathrm{D}}{\bf g}^j = \frac{\partial x^k}{\... ...f e}^l\otimes{\bf e}_m\; \right) \frac{\partial \xi^j}{\partial x^n}{\bf e}^n$$

Using the relation from Eq. 29,

$${\bf g}_i\overline{\mathrm{D}}{\bf g}^j = \frac{\partial x^k}{\... ...i^j}{\partial x^n} \left({\bf e}_k\cdot{\bf e}^l\right)\,\tensor{\delta}{_m^n}$$

$$\Rightarrow {\bf g}_i\overline{\mathrm{D}}{\bf g}^j = \frac{\p... ...rac{\partial \xi^j}{\partial x^n} \tensor{\delta}{_k^l}\,\tensor{\delta}{_m^n}$$

$$\Rightarrow {\bf g}_i\overline{\mathrm{D}}{\bf g}^j = \frac{\p... ...i}\; \tensor{\overline{\mathrm{D}}}{_k^m} \;\frac{\partial \xi^j}{\partial x^m}$$

$$\Rightarrow {\bf g}_i\overline{\mathrm{D}}{\bf g}^j = \frac{\p... ...} \frac{\partial x^k}{\partial \xi^i}\; \tensor{\overline{\mathrm{D}}}{_k^l}$$

which corresponds to the relation given in Eq. 23. The second equality is then studied:

$${\bf g}^j\otimes{\bf g}_i\,\overline{\mathrm{D}} = \left[ \lef... ...tensor{\overline{\mathrm{D}}}{_m^n} \left( {\bf e}^m\otimes{\bf e}_n \right)$$

$$\Rightarrow {\bf g}^j\otimes{\bf g}_i\,\overline{\mathrm{D}} = ... ... \otimes {\bf e}_l \right) \cdot \left( {\bf e}^m\otimes{\bf e}_n \right)$$

Using the relation on the tensorial product given in Eq. 30, the equation is simplified to:

$${\bf g}^j\otimes{\bf g}_i\,\overline{\mathrm{D}} = \frac{\parti... ...m{D}}}{_m^n} \tensor{\delta}{_l^m} \left( {\bf e}^k\otimes{\bf e}_n \right)$$

$$\Rightarrow {\bf g}^j\otimes{\bf g}_i\,\overline{\mathrm{D}} = ... ...tensor{\overline{\mathrm{D}}}{_l^n} \left( {\bf e}^k\otimes{\bf e}_n \right)$$

The following property of the trace operator is also used $\mathrm{tr}\left({\bf a}\otimes{\bf b}\right)={\bf a}\cdot{\bf b}$ to have finally the expression of the tensor:

$$\mathrm{tr}\left({\bf g}^j\otimes{\bf g}_i\,\overline{\mathrm{D}}... ...l}{\partial \xi^i} \tensor{\overline{\mathrm{D}}}{_l^n} \tensor{\delta}{^k_n}$$

$$\Rightarrow \mathrm{tr}\left({\bf g}^j\otimes{\bf g}_i\,\overlin... ...x^k} \frac{\partial x^l}{\partial \xi^i} \tensor{\overline{\mathrm{D}}}{_l^k}$$

$$\Rightarrow \mathrm{tr}\left({\bf g}^j\otimes{\bf g}_i\,\overlin... ...x^l} \frac{\partial x^k}{\partial \xi^i} \tensor{\overline{\mathrm{D}}}{_k^l}$$

which is also the relation of the tensor given in Eq. 23

Proof of Eq. 26

It is proven in this appendix that ${\bf g}^{i}_{,j}=-\tensor{\Gamma}{^i_{kj}}{\bf g}^k$ . Using the relation ${\bf g}^i\cdot{\bf g}_j=\tensor{\delta}{^i_j}$ given by derivation

$$\forall (i,j,k)\in[\![0,3]\!]\;\;,\; \left({\bf g}^i\cdot{\bf g}_j\right)_{,k}=0$$

$$\Rightarrow {\bf g}^{i}_{,k}\cdot{\bf g}_j=-{\bf g}^i\cdot{\bf g}_{j,k}$$

$$\Rightarrow {\bf g}^{i}_{,k}=-\left({\bf g}^i\cdot{\bf g}_{j,k}\right)\,{\bf g}^j$$

The Christoffel symbol is introduced with the relation of Eq 26.

$${\bf g}^{i}_{,k}=-\left({\bf g}^i\cdot\left(\tensor{\Gamma}{^l_{jk}}{\bf g}_{l}\right)\right)\,{\bf g}^j$$

$$\Rightarrow {\bf g}^{i}_{,k}=-\tensor{\Gamma}{^l_{jk}}{\bf g}^j\tensor{\delta}{^i_l}$$

$$\Rightarrow {\bf g}^{i}_{,j}=-\tensor{\Gamma}{^i_{kj}}{\bf g}^k$$

which ends the proof.

Proof of Eq. 41

It is shown that when the differentiation is performed on a Cartesian grid:

$$\nabla\cdot\left(\mathrm{D}\nabla \phi\right) = \tensor{\mathrm{D}}{_i^j_{,i}}\,\phi_{,j}+\tensor{\mathrm{D}}{_i^j}\phi_{,ij}$$

where $u_{,i}=\frac{\partial u}{\partial X^i}$ and $u$ can represent $\mathrm{D}$ or $\phi$ .

Since the Cartesian coordinate system of basis $({\bf e}_1,{\bf e}_2,{\bf e}_3)$ is orthonormal, the contravariant and covariant vectors are equal meaning $({\bf e}_1,{\bf e}_2,{\bf e}_3)=({\bf e}^1,{\bf e}^2,{\bf e}^3)$ . Moreover, this basis doesn't depend on the position ${\bf X}$ , therefore:

$$\forall (i,j)\in[\![0;3]\!]\;\;,\; {\bf e}_{i,j}={\bf e}^{i}_{,j}={\bf0}$$

The first part $\mathrm{D}\,\nabla\phi$ can be expressed as:

$$\mathrm{D}\,\nabla\phi = \tensor{\mathrm{D}}{_i^j}\left({\bf e}^{i}\otimes {\bf e}_{j}\right)\;{\bf e}^{k}\frac{\partial \phi}{\partial X^k}$$

The relation given by Eq. 29 can then be recognized. Therefore $\left({\bf e}^{i}\otimes {\bf e}_{j}\right)\;{\bf e}^{k}={\bf e}^{i}\left({\bf e}_j\,{\bf e}^k\right)={\bf e}^{i}\,\tensor{\delta}{_j^k}$ , where $\tensor{\delta}{_j^k}=1$ only if $j=k$ and equals 0 otherwise. Then the first part of the diffusion is given by:

$$\mathrm{D}\,\nabla\phi = \tensor{\mathrm{D}}{_i^j}\phi_{,k}\tensor{\delta}{_j^k}\,{\bf e}^i$$

$$\Rightarrow \mathrm{D}\,\nabla\phi = \tensor{\mathrm{D}}{_i^j}\phi_{,j}\,{\bf e}^i$$

The complete diffusion equation is obtained by taking the divergence of this expression:

$$\nabla\cdot\left(\mathrm{D}\,\nabla\phi\right)= {\bf e}^{i}\frac... ...{i}} \cdot \left( \tensor{\mathrm{D}}{_j^k}\,\phi_{,k}\;{\bf e}^{j} \right)$$

$\mathrm{D}$ and $\phi$ depend on the coordinates but ${\bf e}$ doesn't.

$$\nabla\cdot\left(\mathrm{D}\,\nabla\phi\right)= \left( \tensor{... ...{_j^k}\,\phi_{,ki}\; \right) \left( {\bf e}^{i} \cdot {\bf e}^{j} \right)$$

Assuming the scalar $\phi$ , and therefore that each component of the magnetization vector ${\bf M}$ are differentiable twice to enable the Schwarz's theorem on the mixed derivative to hold [62] such that $\phi_{,ki}=\phi_{,ik}$ , the diffusion process is given by:

$$\nabla\cdot\left(\mathrm{D}\,\nabla\phi\right)= \left( \tensor{... ...hi_{,k}+\tensor{\mathrm{D}}{_j^k}\,\phi_{,ik}\; \right) \tensor{\delta}{^i^j}$$

where, in Cartesian coordinates, $\tensor{\delta}{^i^j}$ equals one only if $i=j$ and 0 otherwise. The diffusion is therefore given by the relation of Eq. 41:

$$\nabla\cdot\left(\mathrm{D}\,\nabla\phi\right)= \tensor{\mathrm{D}}{_i^k_{,i}}\,\phi_{,k}+\tensor{\mathrm{D}}{_i^k}\,\phi_{,ik}\;$$

Proof of Eq. 43

It is proven that in curvilinear coordinates $(\xi^1,\xi^2,\xi^3)$ with covariant basis vectors ${\bf g}^1,{\bf g}^2,{\bf g}^3$ , contravariant basis vectors ${\bf g}_1,{\bf g}_2,{\bf g}_3$ , and contravariant metric tensor $\mathrm{g}^{ij}$ , the expression of the diffusion is given by:

$$\nabla\cdot\left(\mathrm{D}\,\nabla\phi\right) = \left[ \left(... ...{\Gamma}{^l_{ji}}\tensor{\mathrm{D}}{_l^k}\,\phi_{,k} \right] \mathrm{g}^{ij}$$

where $\tensor{\Gamma}{^k_{ij}}$ is the Christoffel symbol of second kind. The first part of the equality $\mathrm{D}\nabla\phi$ follows the same proof as for the Cartesian case:

$$\mathrm{D}\nabla\phi = \tensor{\mathrm{D}}{_i^j}\left({\bf g}^i\otimes{\bf g}_j\right)\cdot{\bf g}^k\frac{\partial \phi}{\partial \xi^k}$$

$$\Rightarrow \mathrm{D}\nabla\phi = \tensor{\mathrm{D}}{_i^j}\phi_{,k}\;{\bf g}^i\tensor{\delta}{_j^k}$$

$$\Rightarrow \mathrm{D}\nabla\phi = \tensor{\mathrm{D}}{_i^j}\phi_{,j}\;{\bf g}^i$$

which corresponds to the expression given in Eq. 42. The next part of the derivation takes the divergence of this vector such that:

$$\nabla\cdot\left(\mathrm{D}\nabla\phi\right) = {\bf g}^i\frac{\... ...l \xi^i} \cdot \left( \tensor{\mathrm{D}}{_j^k}\phi_{,k}\;{\bf g}^j \right)$$

The derivative is expanded on the product where the vector ${\bf g}^j$ depends on the variable ${\bf x}$ :

$$\nabla\cdot\left(\mathrm{D}\nabla\phi\right) = {\bf g}^i \left... ...{,i}\;{\bf g}^j + \tensor{\mathrm{D}}{_j^k}\phi_{,k}\;{\bf g}^j_{,i} \right]$$

The derivative of the covariant basis vector given in Eq. 4.2.1 is inserted into the relation:

$$\nabla\cdot\left(\mathrm{D}\nabla\phi\right) = {\bf g}^i \left... ...tensor{\mathrm{D}}{_j^k}\,\phi_{,k}\;\tensor{\Gamma}{^j_{li}}{\bf g}^l \right]$$

$$\Rightarrow \nabla\cdot\left(\mathrm{D}\nabla\phi\right) = \le... ... \tensor{\mathrm{D}}{_j^k}\,\phi_{,k}\;\tensor{\Gamma}{^j_{li}}\mathrm{g}^{il}$$

where the indexes $l$ and $j$ can be inverted in the last part of the sum:

$$\Rightarrow \nabla\cdot\left(\mathrm{D}\nabla\phi\right) = \le... ...{\Gamma}{^l_{ji}}\tensor{\mathrm{D}}{_l^k}\,\phi_{,k} \right] \mathrm{g}^{ij}$$

which ends the proof. Finally, the expression can also be seen in its expanded form:

$$\Rightarrow \nabla\cdot\left(\mathrm{D}\nabla\phi\right) = \le... ...{\Gamma}{^l_{ji}}\tensor{\mathrm{D}}{_l^k}\,\phi_{,k} \right) \mathrm{g}^{ij}$$

Bibliography

1

Patrick Helm, Mirza Faisal Beg, Michael I.Miller, and Raimond L.Winslow.
Measuring and mapping cardiac fiber and laminar architecture using diffusion tensor MR imaging.
Annals of the New York Academy of Sciences, 1047:296-307, June 2005.

2

E. W. Hsu, A. L. Muzikant, S. A. Matulevicius, R. C. Penland, and C. S. Henriquez.
Magnetic resonance myocardial fiber-orientation mapping with direct histological correlation.
American Journal of Physiology. Heart and Circulatory Physiology, 274(5):1627-1634, May 1998.

3

Alex Holmes.
Development of Diffusion Tensor Magnetic Resonance Imaging Techniques for Rapid Determination of Myocardial Fiber Orientation.
PhD thesis, Johns Hopkins University, Baltimore, Maryland, 1999.

4

Yi Jiang, Kumar Pandya, Oliver Smithies, and Edward W. Hsu.
Three-dimensional diffusion tensor microscopy of fixed mouse hearts.
Magnetic Resonance in Medicine, 52(3):453-460, September 2004.

5

David Scollan.
Reconstructing the Heart: Development and Application of Biophysically-Based Electrical Models of Propagation in Ventricular Myocardium Reconstructed from Diffusion Tensor MRI.
PhD thesis, Johns Hopkins University, Baltimore, Maryland, 2002.

6

Vitaly J. Napadow, Qun Chen, Vu Mai, Peter T. C. So, and Richard J. Gilbert.
Quantitative analysis of three-dimensional-resolved fiber architecture in heterogeneous skeletal muscle tissue using NMR and optical imaging methods.
Biophysical Journal, 80(6):2968-2975, June 2001.

7

J. E. Tanner.
Self diffusion of water in frog muscle.
Biophysical Journal, 28(1):107-116, October 1979.

8

G. G. Cleveland, D. C. Chang, C. F. Hazlewood, and H. E. Rorschach.
Nuclear magnetic resonance measurement of skeletal muscle. anisotropy of the diffusion coefficient of the intracellular water.
Biophysical Journal, 16(9):1043-1053, September 1976.

9

C. C. Van Donkelaar, L. J. Kretzers, P. H. M. Bovendeerd, L. M. A. Lataster, K. Nicolay, J. D. Janseen, and M. R. Drost.
Diffusion tensor imaging in biomechanical studies of skeletal muscle function.
Journal of ANATOMY, 194(1):79-88, January 1999.

10

Denis Le Bihan.
Looking into the functional architecture of the brain with diffusion MRI.
Nature Reviews. Neuroscience, 4(6):459-80, June 2003.

11

Denis Le Bihan, Jean-Francois Mangin, Cyril Poupon, Chris A. Clark, Sabina Pappata, Nicolas Molko, and Hughes Chabriat.
Diffusion tensor imaging: Concepts and applications.
Journal of Magnetic Resonance Imaging, 13(4):534-546, April 2001.

12

Oliver Natt and Jens Frahm.
In vivo magnetic resonance imaging: Insights into structure and function of the central nervous system.
Measurement Science and Technology, 16:R17-R36, February 2005.

13

Thomas L. Chenevert, James A. Brunberg, and James G. Pipe.
Anisotropic diffusion in human white matter: Demonstration with MR techniques in vivo.
Radiology, 177(2):401-5, November 1990.

14

Richard J. Gilbert, Lee H. Magnusson, Vitaly J. Napadow, Thomas Benner, Ruopeng Wang, and Van J. Wedeen.
Mapping complex myoarchitecture in the bovine tongue with diffusion-spectrum magnetic resonance imaging.
Biophysical Journal, 91:1014-1022, May 2006.

15

Richard J. Gilbert and Vitaly J. Napadow.
Three-dimensional muscular architecture of the human tongue determined in vivo with diffusion tensor magnetic resonance imaging.
Dysphagia, 20(1):1-7, 2005.

16

Van J. Wedeen, Timothy G. Reese, Vitaly J. Napadow, and Richard J. Gilbert.
Demonstration of primary and secondary muscle fiber architecture of the bovine tongue by diffusion tensor magnetic resonance imaging.
Biophysical Journal, 80(80):1024-1028, February 2001.

17

Chad A. Holder, Raja Muthupillai, Srinivasan Mukundan, James D. Eastwood, and Patricia A. Hudgins.
Diffusion-weighted MR imaging of the normal human spinal cord in vivo.
American Journal of Neuroradiology, 21(10):1799-1806, November 2000.

18

Joon W. Lee, Jae H. Kim, Heung S. Kang, Jong S. Lee, Ja-Young Choi, Jin-Sup Yeom, Hyun-Jib Kim, and Hye Chung.
Optimization of acquisition parameters of diffusion-tensor magnetic resonance imaging in the spinal cord.
Investigative Radiology, 41(7):553-559, July 2006.

19

Robert L. Cooper, David B. Chang, Allan C. Young, Carroll J. Martin, and Betsy Ancker Johnson.
Restricted diffusion in biophysical systems.
Biophysical Journal, 14(3):161-177, March 1974.

20

Sverre Rosenbaum.
Evaluation of Human Stroke by MR Imaging.
PhD thesis, Faculty of Medicine, University of Copenhagen, Denmark, August 2000.

21

Josef Pfeuffer, Wolfgang Dreher, Eva Sykova, and Dieter Leibfritz.
Water signal attenuation in diffusion-weighted 1h NMR experiments during cerebral ischemia: Influence of intracellular restrictions, extracellular tortuosity, and exchange.
Magnetic Resonance Imaging, 16(9):1023-1032, November 1998.

22

Christian Beaulieu.
The basis of anisotropic water diffusion in the nervous system - a technical review.
NMR in Biomedicine, 15(7-8):435-55, November 2002.

23

Burkhard Wunsche.
The visualization and measurement of left ventricular deformation.
In Proceedings of the First Asia-Pacific Bioinformatics Conference on Bioinformatics, pages 119-128. Australian Computer Society, Inc., 2003.

24

Leonid Zhukov and Alan H. Barr.
Heart-muscle fiber reconstruction from diffusion tensor MRI.
IEEE , Proceedings of the 14th IEEE Visualization, page 79, October 2003.

25

Damien Rohmer, Arkadiusz Sitek, and Grant T. Gullberg.
Reconstruction and visualization of fiber and sheet structure with regularized tensor diffusion MRI in the human heart.
Technical Report LBNL-60277, Lawrence Berkeley National Laboratory Report, 2006.

26

Damien Rohmer, Arkadiusz Sitek, and Grant T. Gullberg.
Fiber tracking in the entire heart.
Available at http://cfi.lbl.gov/ asitek/fibertrack /3_fiber_tracking_entire_heart.html, 2006.

27

Patrick A. Helm, Laurent Younes, Mirza F. Beg, Daniel B. Ennis, Christophe Leclercq, Owen P. Faris, Elliot McVeigh, David Kass, Michael I. Miller, and Raimond L. Winslow.
Evidence of structural remodeling in the dyssynchronous failing heart.
Circulation Research, 98(1):125-132, January 2006.

28

Joseph C. Walker, Julius M. Guccione, Yi Jiang, Peng Zhang, Arthur W. Wallace, Edward Hsu, and Mark B. Ratcliffe.
Helical myofiber orientation after myocardial infarction and left ventricular surgical restoration in sheep.
Journal of Thoracic and Cardiovascular Surgery, 129(2):382-390, February 2005.

29

Marc A. Pfeffer and Eugene Braunwald.
Ventricular remodeling after myocardial infarction. experimental observations and clinical implications.
Circulation, 81(4):1161-1172, April 1990.

30

Jiangang Dou, Timothy G. Reese, Wen-Yih I. Tseng, and Van J. Wedeen.
Cardiac diffusion MRI without motion effects.
Magnetic Resonance in Medicine, 48(1):105-114, July 2002.

31

Wen-Yih I. Tseng, Timothy G. Reese, Robert M. Weisskoff, Thomas J. Brady, and Van J. Wedeen.
Myocardial fiber shortening in humans: Initial results of MR imaging.
Radiology, 216(1):128-139, July 2000.

32

Wen-Yih I. Tseng, Timothy G. Reese, Robert M. Weisskoff, and Van J. Wedeen.
Cardiac diffusion tensor MRI in vivo without strain correction.
Magnetic Resonance in Medicine, 42(2):393-403, July 1999.

33

Timothy G. Reese, Robert M. Weisskoff, R. Neil Smith, Bruce R. Rosen, Robert E. Dinsmore, and Van J. Wedeen.
Imaging myocardial fiber architecture in vivo with magnetic resonance.
Magnetic Resonance in Medicine, 34(6):786-91, December 1995.

34

Bing Feng, Alexander I. Veress, and Grant T. Gullberg.
The prolate spheroidal transform for gated SPECT.
IEEE Transaction on Nuclear Science, 48:872-875, June 2001.

35

Wolfram Neutsch.
Coordinates.
Walter de Gruyter, 1996.

36

Pabitra N. Sen and Peter J. Basser.
A model for diffusion in white matter in the brain.
Biophysical Journal, 89(5):2927-38, November 2005.

37

Pabitra N. Sen.
Diffusion and tissue microstructure.
Journal of Physics: Condensed Matter, 16:S5213-20, November 2004.

38

Mihail S. Iotov.
Diffusion in Amorphous Media.
PhD thesis, California Institute of Technology, Pasadena California, September 1997.

39

S. V. Dvinskikh and I. Furo.
Anisotropic self-diffusion in the nematic phase of a thermotropic liquid crystal by 1H-spin-echo nuclear magnetic resonance.
Journal of Chemical Physics, 115(4):1946-50, July 2001.

40

Mette R. Wiegell, Henrik B. W. Larsson, and Van J. Wedeen.
Fiber crossing in human brain depicted with diffusion tensor MR imaging.
Radiology, 217(3):897-903, December 2000.

41

David S. Tuch, Timothy G. Reese, Mette R. Wiegell, and Van J. Wedeen.
Diffusion MRI of complex neural architecture.
Neuron, 40(5):885-895, December 2003.

42

Evren Ozarslan and Thomas H. Mareci.
Generalized diffusion tensor imaging and analytical relationships between diffusion tensor imaging and high angular resolution diffusion imaging.
Magnetic Resonance in Medicine, 50(5):955-965, November 2003.

43

David Solomon Tuch.
Diffusion MRI of Complex Tissue Structure.
PhD thesis, Massachusetts Institute of Technology, Massachusetts, USA, January 2002.

44

J. R. Basser and D. K. Jones.
Diffusion-tensor MRI: theory, experimental design and data analysis - a technical review.
NMR in Biomedicine, 15(7-8):456-467, November 2002.

45

Bernd Jahne.
Digital Image Processing.
Springer, sixth edition, 2005.

46

Dwight G. Nishimura.
Principles of Magnetic Resonance Imaging.
Stanford University, California, 1996.

47

F. Bloch.
Nuclear induction.
Pysical Review, 70(7):460-474, October 1946.

48

H. C. Torrey.
Bloch equation with diffusion terms.
Phys Rev, 104(3):563-565, November 1956.

49

E.O. Stejskal and J.E. Tanner.
Spin diffusion measurements: Spin echoes in the presence of a time-dependent field gradient.
The Journal of Chemical Physics, 42(1):288-292, January 1964.
Spin Echo sequence from Stejskal and Tanner for diffusion tensor MRI.

50

Cyril Poupon.
Détection des faisceaux de fibres de la substance blanche pour l'étude de la connectivité anatomique cérébrale.
PhD thesis, Ecole Nationale Superieure des Telecommunications de Paris, Paris, France, 1999.

51

T. J. Chung.
Applied Continuum Mechanics.
Cambridge, 1997.

52

A. I. Borisenko and I. E. Tarapov.
Vector and Tensor Analysis. With Applications.
Dover, 1979.

53

Erwin Kreysig.
Differential Geometry.
Dover, 1991.

54

Jay D. Humphrey.
Cardiovascular Solid Mechanics. Cells, Tissues, and Organs.
Springer, 2001.

55

S. P. Timoshenko and J. N. Goodier.
Theory of Elasticity.
McGraw Hill, third edition, 1987.

56

Timothy G. Reese, Van J. Wedeen, and Robert M. Weisskoff.
Measuring diffusion in the presence of material strain.
Journal of Magnetic Resonance. Series B., 112(3):253-258, September 1996.

57

Alfio Quarteroni and Alberto Valli.
Numerical Approximation of Partial Differential Equations.
Springer, 1997.

58

K. W. Morton and D. F. Mayers.
Numerical Solution of Partial Differential Equations.
Cambridge, 2005.

59

E. Haier, S. P. Norsett, and G. Wanner.
Solving Ordinary Differential Equations I, Nonstiff Problems.
Springer, second edition, 2000.

60

J. W. Thomas.
Numerical Partial Differential Equations. Finite Difference Methods.
Springer, 1995.

61

Alfio Quarteroni, Riccardo Sacco, and Fausto Saleri.
Numerical Mathematics.
Springer, 1991.

62

I. N. Bronshtein, K. A. Semendyayev, G. Musiol, and H. Muehlig.
Handbook of mathematics.
Springer, fourth edition, 2003.

About this document ...

A Bloch-Torrey Equation for Diffusion in a Deforming Media

This document was generated using the LaTeX2HTML translator Version 2002-2-1 (1.71)

Copyright © 1993, 1994, 1995, 1996, Nikos Drakos, Computer Based Learning Unit, University of Leeds.
Copyright © 1997, 1998, 1999, Ross Moore, Mathematics Department, Macquarie University, Sydney.

The command line arguments were:
latex2html MRI_diffusion.tex -split 0

The translation was initiated by damien on 2007-03-03

Footnotes

... performed¹

Trace and determinant are independent of the basis.