ASSESSMENT OF THE EFFECT OF TISSUE MOTION IN DIFFUSION MRI: DERIVATION OF NEW APPARENT DIFFUSION COEFFICIENT FORMULA

. We investigate in this paper the diﬀusion magnetic resonance imaging (MRI) in deformable organs such as the living heart. The diﬃculty comes from the hight sensitivity of diﬀusion measurement to tissue motion. Com-monly in literature, the diﬀusion MRI signal is given by the complex mag- netization of water molecules described by the Bloch-Torrey equation. When dealing with deformable organs, the Bloch-Torrey equation is no longer valid. Our main contribution is then to introduce a new mathematical description of the Bloch-Torrey equation in deforming media. In particular, some numer- ical simulations are presented to quantify the inﬂuence of cardiac motion on the estimation of diﬀusion. Moreover, based on a scaling argument and on an asymptotic model for the complex magnetization, we derive a new apparent diﬀusion coeﬃcient formula. Finally, some numerical experiments illustrate the potential of this new version which gives a better reconstruction of the diﬀusion than using the classical one.


1.
Introduction. Diffusion magnetic resonance imaging (diffusion MRI) is an imaging technique that is capable of providing MRI images with a contrast sensitive to the diffusion motion of water molecules [8]. Diffusion MRI plays a very important role in studying the microscopic structure of biological tissues through the measure of the diffusion coefficient of water molecules in the tissues. This technique was successfully applied to static organs such as in brain imaging ( [6], [30], [7]). However, its implementation on moving organs like the beating heart is very difficult because of the cardiac motion during acquisition. The heart has a pseudoperiodic and fast motion that can be tracked with dynamic series of images ( [23], [31]). The tissue motion induces an attenuation of the signal measured in diffusion MRI which can considerably degrade the quality of estimation of the diffusion coefficient. Because of the sensitivity to motion, it is difficult to assess to what extent the diffusion characteristics obtained from diffusion MRI reflect the real properties of the cardiac tissues. Some experimental studies have investigated the influence of cardiac motion on the diffusion measurements ( [19]- [18]).
The signal measured in diffusion MRI is described by the magnetization density of the water molecules which can be modeled by the complex valued Bloch-Torrey partial differential equation (PDE) [27]. This equation provides a mathematical description of the diffusion phenomenon of water molecules in the field of magnetic resonance imaging. In static organs, such as the human brain, several approaches have been proposed in the literature to model the diffusion MRI signal from the Bloch-Torrey equation. These different approaches offer the possibility of understanding the link between the diffusion of water molecules and the geometric structure of biological tissues. For example, analytical ( [24]- [16]) and numerical models ( [32]- [13]) of the solution of the Bloch-Torrey equation have been introduced under simplified assumptions on the geometry of spatial domains in which the diffusion process takes place. Recently, a multi-compartment model for the Bloch-Torrey equation has been developed in [14] within a more complex geometry and under more realistic experimental conditions. However, very few models have been proposed to study the influence of physiological motion on the signal measured in diffusion MRI. For example, in [21] the Bloch-Torrey equation was expressed in generalized curvilinear coordinates to describe the behavior of magnetization in the heart during its deformation over the cardiac cycle and a change of basis formulas were used in order to take into account the effect of motion on diffusion.
The focus of this work is firstly presenting a mathematical model based on the Bloch-Torrey equation, which is capable of taking into account the physiological motion of the heart, and allowing the quantification of tissue deformation on the complex magnetization. This model is obtained by writing the Bloch-Torrey PDE in a domain that deforms over time according to the laws of continuum mechanics. Then from the resulting model, the diffusion MRI images are simulated at different moments in the cardiac cycle by using an analytical motion field mimicking a realistic deformation of the heart, and this to validate the model by comparing the results of numerical simulations with the existing experimental results. Secondly, we focus on the mitigation issue of the effect of cardiac motion on the diffusion MRI images. For this purpose, we propose a method for correcting the diffusion MRI images which were reconstructed in the presence of cardiac motion. This method is derived from a singular perturbative model based on an asymptotic analysis of the Bloch-Torrey equation with motion. The resulting model provides an asymptotic model of ordinary differential equation (ODE) that allows to derive a new formula of apparent diffusion coefficient for the diffusion correction and reconstruction of diffusion MRI images free form motion effect. This paper is organized as follows, in Section 2 the Bloch-Torrey equation in static medium is recalled and diffusion measurement technique used in practice is described. The Bloch-Torrey equation established in the presence of physiological motion is introduced in Section 3. Section 4 is dedicated to numerical simulation of the diffusion MRI images during the cardiac cycle. These images are reconstructed by means of the complex magnetization obtained by solving the Bloch-Torrey equation with motion by a finite element method. In section 5, we present a singular perturbation method for correcting the diffusion MRI images that are affected by heart motion during the cardiac cycle. This method provides an asymptotic ODE model that allows to derive a relationship between the complex magnetization influenced by motion and the exact diffusion. This relationship is then used in numerical simulations to correct the diffusion images reconstructed in the presence of motion. C of water molecules including the effect of molecular diffusion. This equation is given in [27] for all x ∈ Ω by: where x = (x, y, z) t is the spatial position, T 2 is a physical parameter called transverse relaxation time, and D(x) represents a second order diffusion tensor which is a 3x3 symmetric positive definite matrix: where D xx is the diffusion coefficient of water molecules in the x-direction and D xy is the correlation between the diffusion coefficients in the x and y-directions, etc. In diffusion MRI, the water molecules are subjected to a time dependent magnetic field gradient (called diffusion encoding gradient): in order to encode the diffusion motion of water molecules in the different directions of space. The Bloch-Torrey equation reads then as: where γ is the gyromagnetic ratio of the hydrogen atom such that γ 2π = 42.58 MHz/Tesla. To measure the diffusion motion of water molecules with magnetic resonance imaging, a special shape of diffusion encoding gradients is performed. This shape is usually the spin echo diffusion gradient sequence introduced by Stejskal and Tanner in [24]. The general idea of this sequence is to excite the water molecules with a 90 o radio-frequency (RF) pulse [15], dephase them and encode their positions in a specific spatial direction with a time-constant magnetic field gradient of intensity G and duration δ. After that, a 180 o RF pulse [15] is applied to invert the phase of magnetization of water molecules. A second diffusion gradient with equal intensity and duration to those of the first gradient is applied at a time ∆ after the first gradient to rephase the water molecules (implying an effect of −G with the 180 o RF pulse). The signal is finally acquired at a time called echo time T E (Fig.  1). Noting that to produce a signal at the time t = T E, we have the condition that If the water molecules move due to diffusion during the time ∆, then the effect of the second gradient will not be exactly opposite to the first gradient. This leads to a phase incoherence of water molecules resulting in a reduction of the magnitude of the complex magnetization M xy , and consequently, the signal will be smaller relative to the case in which there was no diffusion encoding gradients. The decrease in the acquired signal will reflect the amount of diffusion which occurred during the application of the diffusion encoding gradients. Hence, measuring the signal decay makes possible the measure of the diffusion coefficient in the direction of application of the diffusion gradient sequence. In order to estimate the diffusion coefficient, two MRI images are acquired, one with diffusion encoding gradient and the another without diffusion encoding gradient. These two images correspond to the MRI signals calculated over each voxel V at the echo time T E with diffusion gradient as well as with no diffusion gradient: Then, the diffusion coefficient can be estimated for each point of the MRI image through the so-called Apparent Diffusion Coefficient (ADC) ( [24], [9]) as: with the diffusion weighting factor (or b-value): In the case of anisotropic diffusion, the apparent diffusion tensor D A can be estimated according to the relation [11]: where the diffusion weighting B-matrix is given by: The apparent diffusion tensor D A is a 3x3 symmetric positive definite matrix. This implies that there are six apparent diffusion coefficients to be estimated. To do so, a set of diffusion MRI signals is acquired by the application of the diffusion encoding gradients in at least six non-coplanar spatial directions to get a linear system of at least six equations that can be resolved to find the coefficients of the tensor D A .
3.1. Mathematical model. In this section we consider an evolving set Ω(t) ⊂ R d (d=2 or 3) for all t ∈ [0, T ] with T > 0. Let us introduce the geometric transformation ϕ which is a regular time-space dependent function: and assume that for any point x, the curve t → (t, x) satisfies the system: for a given velocity field v : R d → R d . We assume that the evolving set {Ω(t)} t∈[0,T ] is defined from an initial domain Ω(0) ⊂ R d by: The following notations will be used in this section: If for some time T > 0, a function g is defined in Ω(t) as: then it can be defined in Ω(0) as: We now give a description of the Bloch-Torrey equation in the moving set Ω(t). The idea consists to analyze the magnetic flux of M in a elementary volume V (t) ⊂ R d which evolves by the geometric transformation ϕ from an elementary volume V (0) at time t 0 = 0. We denote by M (X, t) the magnetization at X ∈ V (t) at time t. In order to establish the conservation law for the magnetization density we have to specify the change of magnetization in the volume V (t) which is given by d dt V (t) M (X, t)dX. We write the equation of conservation taking into account losses J due to the flux through the boundary of V (t), which is assumed to be: where D(X) is the diffusion tensor at X. We get: where n X is the normal unit vector to the surface dS at point X. By using the divergence theorem we get: Proposition 1. If the magnetization M satisfies Equation (4) in the domain V (t), then it satisfies the following equivalent equations that describe M in the deformed and initial configurations: where F(t, x) := Dϕ(t, x) is the Jacobian matrix of ϕ and F −t the transpose of F −1 .
Proof. By Reynolds transport theorem we have: By identification with Eq. (4) we find: Hence equation (5) is obtained for every volume V (t). Now to get Eq. (6), Eq. (4) is expressed in the domain V (0). We use the change of variable formulas for the gradient and divergence [4]: and the Liouville-Jacobi formula for the derivative of the determinant: which gives: Then the second member of Eq. (4) can be written as: Now, by applying the change of variable formula for the integrals and the Liouville-Jacobi formula, the first member of Eq. (4) gives: By combining Eqs. (7) and (8), we get: Which is true for every V (0): Equation (9) allows to take into account the effect of motion in the Bloch-Torrey equation (2). The frequency term γx · G(t) in Eq. (2) is transformed into γϕ(t, x) · G(t) according to the motion. Then the Bloch-Torrey equation with motion effect is written as follows: When the frequency term γϕ · G is high, the magnetization M xy becomes highly oscillating and this limits the time step that can be used in the numerical simulations. To address this limitation, Eq. (10) is demodulated by writing the magnetization M xy as follows: Substituting (11) into Eq. (10) and multiplying the resulting equation by exp iγ t 0 ϕ(x, t ) · G(t )dt , we get: where we noted the new diffusion tensor by: and the phase term by: We have: where D∇Φ is the Jacobian matrix of ∇Φ(x, t), and K : (D∇Φ) is the double dot product of K and D∇Φ. By replacing Eqs. (13) and (14) in Eq. (12), we get the following equation for the magnetization m: and Noting that to get the solution of Eq. (10), we solve Eq. (15) to obtain m and by using the expression (11) we deduce the solution of Eq. (10).

3.2.
Existence and uniqueness. Assume that Ω(0) be a bounded Lipschitz subset of R d , and Q = Ω(0) × [0, T ] for T > 0. We associate to the Bloch-Torrey equation (10) an initial condition and homogenous Dirichlet boundary condition. In this case the magnetizations M xy and m have the same initial and boundary values. Then the same initial and boundary conditions can be associated to Eq. : For each t ∈ [0, T ], we denote by L the second order differential operator given by: We recall that the tensor K is given by where F is the Jacobian of the deformation ϕ that can be written in function of the displacement u. We get then F := I + Du.
In order to guarantee the definite positivity of K we suppose that Du(x, t) is small compared to 1 for each x ∈ Ω(0) and t ∈ [0, T ]. This allows to satisfy the following ellipticity hypothesis for the operator L: The coefficients c 1i for i = 1, ..., d, and c 2 are supposed in L ∞ (Q, C). In order to de formulate problem (16) in a weak sense, let v be a test function in V = H 1 0 (Ω(0), C) equipped with the scalar product and the induced norm: By multiplying (16) by v, by integrating by parts on Ω(0), and by using the Green's formula, a weak formulation of problem (16) is obtained: where m 0 ∈ H = L 2 (Ω(0), C) and the sesquilinear form a is given by: The existence and unicity of this problem can be obtained as a consequence of the following theorem.
and satisfies the properties: there exists a constants α > 0 and σ ∈ R such that: Then, for all u 0 ∈ H, the following problem 3.3. Finite element approximation. We note that to calculate the magnetization M xy , we solve numerically Eq. (16) for the magnetization m and by using the relation (11) we deduce M xy . We use two-dimensional Lagrangian finite element method for the resolution of Eq. (16). A weak formulation for this problem is written as follows: where v ∈ V is a real test function. We define a mesh for the domain Ω(0), that enables to introduce an approximation space V h which is a finite dimension subspace of V . To discretize Ω(0), we choose a triangular mesh T h (Fig. 2), where h is the maximum mesh size. The space of approximation V h is the set of continuous functions on Ω(0) and which are polynomials on each triangle K of the mesh. By choosing an approximation by polynomials of degree less than or equal to 1, we have: This space admits a base functions (ϕ j ) j that take the value 1 at node x j and 0 at the other nodes x k of the mesh. We get then the following approximate problem: where m 0h ∈ V h is an approximation of the initial condition m 0 . The unknown function m h is discretized as follows: where m j is the value of the solution m h at the node x j , and N is the total number of nodes. We set v h = ϕ i , the weak formulation in (18) gives: T ], we introduce the matrices: for i, j = 1, . . . , N . The weak formulation (19) can be written in matrix form: , and the initial condition m(0) is given. For the temporal discretization of ODE (20) an Euler implicit scheme is used. The interval [0, T ] is subdivided into n T ∈ N subintervals of equal length ∆t := T /n T . We therefore have the sequence (m n ) 0≤n≤n T satisfying: with A n = A(n∆t) and m 0 = m(0) is given.
4. Application to cardiac diffusion MRI.

4.1.
Analytical 2D short-axis cardiac motion model. To study the effect of cardiac motion on diffusion MRI images, we implement a 2D cardiac image simulator introduced in [2]. This simulator is based on a spatial-temporal analytical cardiac motion model mimicking a realistic deformation of the heart. It can be applied to a real short-axis image of the left ventricle to generate from it a synthetic sequence of MRI images of the beating heart. The resulting images resemble to a real MRI sequence and can be used to control the dynamics of the heart (Fig. 2). The initial form of a short-axis slice of the left ventricle can be modeled by a ring delimited by two circles of internal R int and external R ext radii, respectively (Fig. 2). The heart deformation is defined analytically during the cardiac cycle and controlled by the temporal function S given by: where the time T s is the cardiac contraction duration that is called systolic phase and T d is the cardiac expansion duration known as diastolic phase. The duration of one cardiac cycle in human heart is about 1000ms. The systolic duration covers one third of the cardiac cycle and the diastolic duration covers two thirds (Fig.3). We denote by (R 0 , θ 0 ) the polar coordinates at an initial time t 0 of a material point c 0 in the ring Ω(0) = [R int , R ext ] × [0, 2π], and by (R(t), θ(t)) its coordinates during deformation. The heart deformation ϕ is defined in polar coordinates by [2]: The variation of thickening of the myocardial wall over time is given by λ(t) = 1 − 0.2S(t). The radial contraction is: with g(θ 0 ) = 1 + 0.1(cos(θ 0 + 3π/4) + 1), R min = R int − 5(R ext − R int )/6 is the minimal limit value of the internal radius R int during systole, and Then, the heart deformation field in Cartesian coordinates is:

Diffusion measurements in cardiac diffusion MRI.
There are several forms of diffusion encoding sequences proposed in literature for diffusion measurement in the beating heart. In this paper we choose a diffusion encoding sequence called STEAM (see [20], [29]) for the numerical simulations of the diffusion MRI images. The idea of this sequence is to divide the spin echo diffusion encoding sequence (Fig. 1) into two parts and replace the 180 o RF pulse by two RF pulses of 90 o which will be applied in two successive cardiac cycles (Fig. 4). To use this diffusion encoding sequence, the heartbeat should be regular and the motion is periodic during two consecutive heart beats. The diffusion encoding gradients must be applied exactly at the same time in two successive cardiac cycles in order to fix the cardiac position and compensate the two phase terms that are equal in amplitude and opposite in sign. A duration equal to T E/2 separates between the first and second RF pulse 90 o , and between the third RF pulse and the acquisition of diffusion MRI signal. The time ∆ is equal to the duration of one cardiac cycle, and the two parts of the diffusion encoding sequence are separated by a period called mixing time T M := ∆ − T E/2. In cardiac diffusion imaging, the diffusion encoding sequences are synchronized with the electrocardiogram (ECG). This serves to trigger acquisitions at different times of the cardiac cycle by varying the trigger delay (TD) that corresponds to the time between the beginning of the cardiac cycle and the beginning of the sequence of diffusion encoding magnetic field gradients. Synchronization makes it possible to define an optimal acquisition windows during which the influence of cardiac motion on the diffusion MRI images is reduced.

Numerical simulations.
In this section we study numerically the influence of cardiac motion on diffusion MRI images at different moments of the cardiac cycle. A STEAM diffusion encoding gradients (Fig. 4) is applied and for the cardiac motion, the analytical cardiac deformation field introduced in Section 4.1 is used.  The diffusion MRI images are reconstructed from the acquired diffusion MRI signal at the echo at time T = T E + T M , according to the ADC formula (3), for different values of the trigger delay TD. The resulting images are shown on Figure  6 as well as the exact diffusion coefficient. Noting that these images are presented in the deforming configuration in order to show the results according to the cardiac motion.
The relative error in diffusion is calculated on the whole domain Ω(0), over the cardiac cycle, and presented on Figure 7(a). This result confirms that the diffusion is unaffected by the cardiac deformation at two points in the cardiac cycle called sweet spots that correspond to times when the cardiac deformation approximates its time average during the cardiac cycle. These time points are identified at mid-systole and mid-diastole ( Fig. 7(b)). These numerical simulations were already published in [12] in a comparative framework for evaluating different forms of diffusion encoding sequences used in cardiac diffusion MRI. These results are in agreement with the previous experimental studies [29], [3], [25], in which it has been shown that there exist two time points in the cardiac cycle, sweet spots, relative to which the effect of myocardial deformation on the diffusion measurement is approximately zero. However, it would be necessary to correct the diffusion if one wanted to reconstruct it at times other than sweet spots, and this is what we will address in the next section.

5.1.
Derivation of the asymptotic model. As we have seen in the previous section, the classical ADC formula does not give a sufficiently good reconstruction of the diffusion in the case of moving organ. To address this issue, we propose a new motion correction method of the diffusion MRI images which were reconstructed in the presence of motion. This method is established from an asymptotic model based on a particular scale of coefficients of the Bloch-Torrey equation with motion for the demodulated magnetization m (Eq. (16)).
The derivation of the asymptotic model for the complex magnetization M xy is based on an analysis of order of magnitude of the diffusion D and the gradient of phase ∇Φ which are present in the coefficients of problem (16) for the magnetization m. In the biological tissues, the diffusion D is usually in the order of 10 −8 cm 2 /ms [17]. The phase term Φ in the case of tissue deformation is given by: where G is the diffusion encoding gradient and ϕ is the cardiac deformation. The constant γ is the gyromagnetic ratio of the hydrogen atom that satisfies γ 2π = 42.58MHz/Tesla, and thus γ = 2.67513 × 10 5 rad ms −1 Tesla −1 . We calculate the evolution of the squared norm of ∇Φ(x, t) at different times of the cardiac cycle. The results are shown in Figure 8.
As we can see, the quantity ∇Φ(x, t) 2 2 is in the order of 10 6 . Thus we can choose a strictly positive small parameter ε that is in the order of 10 −4 , and consider D = ε 2D and ∇Φ = 1 ε ∇Φ. Then problem (16) can be transformed into the following problem, where its solution will be denoted m ε : with K = ε 2K and the coefficientsc 1 ,C 2 , and C 3 are given by: Now we study the asymptotic behavior of the solution m ε of problem (21) when ε tends to 0. We write this solution as an asymptotic series with ε: We substitute this approximation in (21), we find that m ε tends formally to m 0 when ε → 0 where m 0 satisfies:

5.2.
Existence and unicity of the solution m . For each time t ∈ [0, T ], we denote by L ε the differential operator given by: where the components ofc 1 , and the coefficientsK ij ,C 2 , C 3 are supposed in L ∞ (Q). Under the assumption of small deformations, as in Section 3.2 the operator L ε , ε > 0 is elliptic in the sense: for t ∈ [0, T ]. A weak formulation for problem (21) can be written as follows: where a(t; ε; m ε , v) : V × V → C, ε > 0, is a family of sesquilinear forms given by: Then as previously, the unique existence of solution u of the problem (24) can be obtained as a consequence of theorem (3.1). notice that the coercivity of a can then be obtained by the following lemma.
Proof. We have: We take v = m , we get: According to Lemma 5.1, the form a is coercive. Then: This implies: According to (25) we have: ∃ε 0 such that ∀ε < ε 0 , σ(ε) ≤ 2 C 3 L ∞ (Q) . By integrating (27) between 0 and T , we get the estimation: with C a constant independant of ε. We deduce that the function m ε est bounded in L 2 ((0, T ), H) independently of ε. By substituting the inequality (27) in (26) and integrating between 0 and T , we find: with C a constant independent of ε. Then α(ε)m ε is bounded in L 2 ([0, T ]; V ), and likewise for the function εm ε since we have ε 2 α < α(ε).
We show now that ∂ t m ε is bounded in L 2 ([0, T ], V ). We have: ) Since 0 < ε ≤ 1, then: By integrating between 0 and T , and using the fact that εm ε L 2 ([0,T ],V ) and m ε L 2 ([0,T ],H) are bounded, we have: We can then extract a subsequence of m ε that converges weakly in L 2 ([0, T ], H) to a function u 1 ∈ L 2 ([0, T ], H), and a subsequence of α(ε)m ε weakly converging to a function u 2 in L 2 ([0, T ], V ). By using the result of compactness of the canonical injection of V in H (Theorem Rellich-Kondrashov), the convergence of α(ε)m ε to u 2 is in L 2 ([0, T ], H). As the convergence of m ε to u 1 is in L 2 ([0, T ], H), and α(ε) tends to 0 with ε, we have that α(ε)m ε tends to 0 in L 2 ([0, T ], H). By uniqueness of the limit, it comes u 2 = 0. Similarly, we have the weak convergence of ∂ t m ε in L 2 ([0, T ], V ) to a function u 3 . By the continuity of derivation in the space distributions D (Q), ∂ t m ε converges weakly to ∂ t u 1 in L 2 ([0, T ], V ). We use these weak convergences to go to the limit as ε tends to 0 in the weak formulation (24): and to show that the second term tends to zero we use the fact that ε∇m ε is bounded in L 2 ([0, T ], H). For the third term we integrate by parts which is possible becausec 1 ∈ C 1 (Q, C) d , then we use the fact that m ε is bounded in L 2 ([0, T ], H) to show that it tends to zero. Finally, we get the problem: where u 1 (·, 0) has a sense because we have u 1 , ∂ t u 1 ∈ L 2 ([0, T ], V ), then u 1 ∈ C 0 ([0, T ], V ). By comparing this problem to Eq. (23), we get u 1 = m 0 . Hence the weak convergences of m ε to m 0 in L 2 ([0, T ], H), and of ∂ t m ε to ∂ t m 0 in L 2 ([0, T ], V ).

5.4.
Analysis of the ADC formula and motion correction. We recall that in the case static domain, the ADC formula reads as: where the diffusion weighting B-matrix is given by: We have seen in the previous section that this formula is not valid in the case of moving media. By using the explicit expression of m 0 and the convergence result of m to m 0 , we can derive the following generalized ADC formula: where T = T E + T M , and B ϕ represents a B-matrix modified according to the deformation ϕ: Indeed, remark that using the asymptotic ODE model (23), the magnetization m 0 satisfies: and then we have the magnetization: From a practical point of view the deformation field ϕ is supposed to be known at each time interval during which the diffusion encoding sequence is applied, and this allows to calculate the modified diffusion weighting matrix B ϕ for the estimation of the diffusion tensor D.

Numerical simulations.
5.5.1. Correction in the case of homogenous diffusion. In Section 4.3 we have seen that the diffusion images (see Fig. 6) are influenced by cardiac motion during the cardiac cycle. This results show the need to correct the diffusion images by comparing them with the image of exact diffusion. In this section the diffusion MRI images will be corrected according to Eq. (29) to reconstruct new diffusion images without motion effect. On figures 9 and 10 we present the diffusion images reconstructed at different times in the cardiac cycle before and after motion correction. These results are shown in the deforming configuration because the image reconstruction is done when the deformation takes place. The images of the error between the exact diffusion and the diffusion reconstructed without motion effect are also shown. The error value is ranging form 0.25×10 −5 to 0.4×10 −5 mm 2 /s. We observe that the corrected diffusion is close to the exact diffusion shown on Figure 6. 5.5.2. Correction in the case of non-homogeneous diffusion. A case of non homogeneous isotropic diffusion in also tested. On Figure 11 the exact diffusion coefficient is shown. Figures 12 and 13 present the reconstruction result of the diffusion MRI images before and after motion correction as well as the induced error, at times TD=250ms and TD=850ms, respectively. The error value ranges form 0.5×10 −5 to 2×10 −5 mm 2 /s. As we can see on the figures, the corrected diffusion obtained in systole and diastole is close to the exact diffusion coefficient.

5.5.3.
Influence of noise and cardiac motion variability on diffusion correction. In this part, we present a numerical test to study the influence of a noisy motion on the reconstructed diffusion images. The observed diffusion image shown in Fig.  12(a) is corrected and presented in Figure 14 when the cardiac motion used in the diffusion correction formula (Eq. (29)) is impacted by a noise for two values of signal-to-noise ration (SNR), 40dB and 30dB. Images of error between the exact and reconstructed diffusion are also presented. The error value ranges form 0.5×0.7 −4 to 5×10 −4 mm 2 /s. On Figure 15 we show corrected diffusion images of the observed diffusion presented on Fig. 12(a) as well as error images when the parameters T s and T d of the cardiac motion are decreased by 10% and 20% with respect to their true values used perviously. Here, the error in diffusion ranges form 0.2×10 −4 to 1.2×10 −4 mm 2 /s.
The results presented in figures 14 and 15 prove that the correction method can improve the diffusion reconstruction compared to the diffusion before correction even in the presence of noise and error committed on the motion. 5.5.4. Numerical influence of the parameter ε. Figure 16 presents diffusion images reconstructed at TD=250ms, before correction, after correction, and error between the exact and the corrected diffusion with different values of the parameter ε of the asymptotic expansion (22) introduced in Section 5.1. These results show that when ε becomes much larger than 10 −4 , the diffusion correction method degrades in efficiency (linearly with respect to ε) because the model that we use is of order of 1 with respect to ε and hence it requires that ε is sufficiently small. Higher order asymptotic model should be used to deal with large values of ε.

5.5.5.
Diffusion correction on an irregular set Ω. Now, a numerical test is performed on an irregular ring in order to approximate a real shape of the left ventricle. Figure 17 represents the distribution of the diffusion coefficient in this domain. The reconstructed diffusion images are shown in Figure 18 as well as the error images in diffusion. The value of error is ranging from 0.5×10 −5 to 1.5×10 −5 mm 2 /s. These images show that even when the domain has an irregular shape, the motion correction method gives a good results, in terms of error values that range the same order of magnitude than for a regular shape domain. We note that to improve the image quality and properly represent the points located on the edge of the domain, the spatial resolution of the images is increased when the reconstruction is done.    6. Conclusion. In this paper, we introduced a modified model of the Bloch-Torrey equation that describes the behavior of the complex magnetization and the diffusion MRI signal in the presence of physiological motion. Then, we quantified, by using numerical simulations, the effect of cardiac motion on the diffusion MRI images reconstructed at different moments of the cardiac cycle. Two time points called sweet spots have been identified in the cardiac cycle, wherein the diffusion images are not affected by cardiac deformation. However, when the diffusion encoding sequence is performed at time points of the cardiac cycle other than sweet spots, a motion correction is necessary after the reconstruction of the diffusion MRI images. For this purpose, we have proposed a motion correction method by means of a new apparent diffusion coefficient formula which is based on an asymptotic ODE model for the complex magnetization. This model has been obtained by an asymptotic analysis of the Bloch-Torrey equation with motion. This motion correction method has been used to correct simulated diffusion MRI images in several situations. This method is valid at any time of the cardiac cycle, since it is based on a mathematical expression that shows explicitly the influence of cardiac deformation on the diffusion. The variability of the parameters of the model has been investigated and the proposed diffusion correction method is shown to be stable.