A SECOND ORDER ASYMPTOTIC MODEL FOR DIFFUSION MRI IN PERMEABLE MEDIA

. Starting from a reference partial differential equation model of the complex transverse water proton magnetization in a voxel due to diffusion-encoding magnetic field gradient pulses, one can use periodic homogenization theory to establish macroscopic models. A previous work introduced an asymptotic model that accounted for permeable interfaces in the imaging medium. In this paper we formulate a higher order asymptotic model to treat higher values of permeability. We explicitly solved this new asymptotic model to obtain a system of ordinary differential equations that can model the diffusion MRI signal and we present numerical results showing the improved accuracy of the new model in the regime of higher permeability

reference fine scale model for the complex-valued transverse water proton magnetization due to diffusionencoding magnetic field gradient pulses. The goal of this paper is to determine higher-order asymptotic expansions of the solution in terms of the ratio of the cell-size to the voxel-size. This ratio is assumed to be small in the dMRI application. The new model of this paper can be seen as a higher order extension of the lower ordered model derived in [8,19]. Similarly to [8,11], we then rely on the interpretation of the macroscopic model in the Fourier domain to derive a system of ordinary differential equations (ODEs) whose solution directly gives the dMRI signal. We then validated the accuracy of the new higher ordered model through some two-dimensional numerical simulations. We show in particular that the new model is more accurate in the regime of higher interface permeability.
The paper is organized as follows. We present the reference fine-scale Bloch-Torrey PDE model and periodic homogenization in Section 2. In Section 3 we provide the key steps to compute the different terms of the two-scale asymptotic expansions and give explicit expressions up to order 2. We present in Section 4 the system of ODEs that results from the order 2 asymptotic model. We give in Section 5 numerical validation using two dimension geometries and in Section 6 we show improvement for higher permeability. For the readers' convenience, we complement our paper with an Appendix A containing some technical details associated with Section 3.

Mathematical context
Let Ω be a domain (open set) containing biological tissue, excluding the boundaries of the biological cells, and let Γ be the union of the cell boundaries. The open set Ω is the union of the extra-cellular domain Ω and the intra-cellular domain Ω : Suppose is the intrinsic diffusion coefficient: The effective time profile, ( ), represents the sequence of applied diffusion-encoding gradient pulses, and obeys the following property: ∫︁ 0 ( ) d = 0, (2.1) where is the time at which the signal is measured called the echo time. A reference model for the complex-valued transverse water proton magnetization = + subject to diffusion-encoding is the following Bloch-Torrey equation [29]: in Ω, (2.2) where is the unit normal vector pointing exterior to the intra-cellular domain, the brackets [·] denote the jump across the interfaces Γ, is the membrane permeability, in the imaginary unit, > 0 is the echo time, init is the initial spin density. The constant vector in R , = , contains the amplitude , which is the product of the diffusion-encoding gradient amplitude multiplied by the gyro-magnetic ratio of the water proton, the unit vector points in the direction of the applied diffusion-encoding magnetic field gradient.
The dMRI signal at echo time is the integral of magnetization ( , ), normalized by the mass: 3) Figure 1. Periodicity box containing 3 biological cells. It is composed of , the exra-cellular domain (red), the intra-cellular domain (green), Γ the boundary of (dark green), and the sides of the box (dark red).
Numerically, the signal is usually plotted against a quantity called the -value

Periodic homogenization
In this work, we will apply periodic homogenization to derive asymptotic solutions of the Bloch-Torrey equation, assuming that the size of the period, , is small. We define the unit periodicity box with the analogous intra-cellular and extra-cellular subsets. See illustration in Figure 1.
The geometric assumption of periodic homogenization is that the domain Ω contains the union of periodic copies of . If the domain Ω is finite, the notation is more complicated than if Ω = R dim . See illustration in Figure 2.
First, we give the notation when Ω is finite. The intra-cellular and extra-cellular domains are defined as periodic extensions of and : The interior boundaries are defined as the periodic extension of Γ : However, to simplify the presentation in this paper by avoiding complications coming from imposing external boundary conditions on Ω, we will choose to take Ω = R , is the space dimension. This assumption is reasonable if the diffusion distance is small compared to the size of the imaging voxel. Finally, on Ω defined as the union of infinite periodic extensions, the PDE becomes: where ( ) := ( ). We assume that the initial data init defined on Ω are independent of and belong to 2 (Ω) ∩ 1 (Ω).

Derivation of asymptotic model
Following the classical periodic homogenization technique [4,17,28], we use the two-scale asymptotic expansions for the solution in Ω and Ω : where the component functions ( , , ) are defined on Ω × ×]0, [, ∈ { , } and are -periodic in . The aim is to separate the macroscopic variations and the microscopic ones to obtain a new problem involving only macroscopic quantities.

Preliminaries
We substitute the asymptotic expansions (3.1) into the periodic model (2.7), and multiply the first equation in (2.7) by 2 to obtain for ( , ) in Ω × . In order to obtain the condition for the traces on the interface Γ , we observe that and By identifying each power of , the no-jump relation for the flux becomes (3.5) Using (3.3), the jump relation for the fields becomes We define the following volume fractions and surface to volume ratios: := |Γ | | | · (3.10)

Notation
In order to shorten the formulas, we respectively denote by :, ∴ and :: the contraction products with respect to the last 2, 3 and 4 indices of the given two tensors. For instance, for any fourth-order tensors and , .
The tensor product of two vectors and is defined by 3.2. Second order asymptotic model under scaling = 0 From homogenization theory, we know that the choice of a scaling for the permeability coefficient has an influence on the nature of the asymptotic model. We take a special form of the scaling of the permeability coefficient corresponding to = 0 , (3.11) that was used previously in [8]. This case corresponds to moderately permeable membranes and leads to a two-compartment macroscopic model [8]. In this paper, we derive higher order asymptotic for this choice of . First we need to define averages and jumps of functions that are solutions of PDEs that will appear later. Now we summarize the equations satisfied bŷ︁ 0 ,̂︁ 1 and̂︁ 2 which are obtained after substituting (3.1) into (2.7) and equating the same powers of . The formal procedure is detailed in Appendix A.
We obtain the homogenized equations for̂︁ 0 : in Ω. (3.16) with the homogenized second-order tensor D defined by where to obtain D we need to compute solutions of time-independent problems posed on the unit box : is -periodic, ⟨ ⟩ = 0, (3.18) with , = 1, . . . , are the vectors of the canonical basis of R .
Likewise, the equations for̂︁ 1 can be written as in Ω. (3.19) No new homogenized tensors are needed.
Finally, the equations for̂︁ 2 involve new homogenized tensors and can be written as (3.20) The needed tensors are B , S , Z , , I and they require the solutions , , and for , , = 1, . . . , . The functions are defined as solutions to the problem (3.18) and the rest of the required functions are solutions to the following box problems: is -periodic, ⟨︀ ⟩︀ = 0, is -periodic, ⟨ ⟩ = 0. (3.23) The symmetric second-order tensor S is defined by The symmetric fourth-order tensor B is the Burnett tensor that is defined by Allaire et al.
where the fourth-order tensor T is defined by The coefficients of the second-order tensor Z , and I are defined by and respectively. We will use ± sign where + corresponds to = and − corresponds to = and ∓ inversely. These simplified expressions (3.26), (3.24), (3.27) and (3.28) are justified in the Appendix A. By combining the equations for̂︁ 0, ,̂︁ 1 and̂︁ 2 given in (3.16), (3.19) and (3.20), we get the following equations for the second order approximation of (2.7), namely, in Ω.
(3.30) We note in the case that the intra-cellular domain is finite, then the tensor D defined by (3.17) is identically zero [8], hence B = −T .

Coupled ODE model for dMRI signal
In this section, we derive a system of ODEs whose solution gives the dMRI signal of the new asymptotic model. The reference normalized signal is measured at the echo time and it is given by is the solution of the reference Bloch-Torrey equation (2.7). For the new second order asymptotic model, we havẽ︀ where ,2 and ,2 are the solutions of (3.30). We now prove that the second order signal obtained from solving (3.30) can be equivalently obtained by solving a simpler system of coupled ODEs.
Theorem 4.1. If we assume that satisfies (2.1), the signal̃︀ 2 is well defined and can computed as︀ where ( , ) is the solution of the second-order ODE model To give an idea of the proof, we first transform the macroscopic model (3.30) by introducing new unknown variables̃︁ ,2 and̃︁ ,2 defined almost everywhere on After the calculations detailed in Section A.4 of the appendix, we obtain the system of equations for̃︁ ,2 and︁ ,2 . After that, we denote by ℳ ,2 and ℳ ,2 the Fourier transform with respect to the space variable of respectivelỹ︁ ,2 and̃︁ ,2 , and denote the dual variable by : We apply the Fourier transform to the system of equations satisfied bỹ︁ ,2 and̃︁ ,2 , and set which can be formally written as We observe that for = 0, the system of equations for ℳ ,2 and ℳ ,2 directly implies the ODE model (4.3).
Finally, for the signal, we have the result that for = , , using ( ) = 0, We refer to Section A.4 in the appendix for more details on the proof of Theorem 4.1.
Remark 4.2. We recall that the Bloch-Torrey signal is in theory complex-valued. The real part of represents the solution due to diffusion, the imaginary part is not relevant to diffusion. We note that the function G 1 , ∈ { , } is imaginary.
For this reason, the second-order ODE model (4.3) can be simplified by removing G 1 and G 1 to get the ODE system (4.11) Suppose one only wants the first order approximation, the signal̃︀ 1 can be computed as︀ where ( , ) is the solution of the first-order ODE model. (4.13) Suppose one wants the zeroth order approximation, the signal̃︀ 0 that can be computed as︀ where ( 0 , 0 ) is the solution of the zeroth-order ODE model (4.15) The above is precisely the asymptotic model derived previously in [8].
Remark 4.3. We recall that the Bloch-Torrey signal is in theory complex-valued. The real part of represents the solution due to diffusion, the imaginary part is not relevant to diffusion. For this reason, the real part of the dMRI signal can be obtained just by solving the zeroth-order ODE model (4.15), there is no need to compute the first-order ODE model (4.13).

Numerical results
Our objective is to validate the accuracy of the new second-order asymptotic model and measure its improvement over the previously derived zeroth-order asymptotic model of [8]. Our reference solution will be provided by the numerical solution of the reference Bloch-Torrey model (2.7). The three models we compare will be denoted by: (1) "Ref": the reference solution from solving the Bloch-Torrey PDE.
(3) " 2": the newly derived second-order asymptotic model given in (4.11). We use the classical Pulsed Gradient Spin Echo (PGSE) sequence, with two rectangular pulses of duration , separated by a time interval ∆ − : We take the echo time = + ∆. We remind that for this time profile the -value is written as We choose the gradient direction to be = . We fix 0 , ∆, and and vary the non-dimensional parameter in order to stay within physically reasonable parameters. The values of the intrinsic diffusivities = 3 × 10 −3 mm 2 /s, = 1.7 × 10 −3 mm 2 /s, are chosen close to the values often used in the literature [13,31]. We validate the convergence in simple two-dimensional geometries, see Figure 3.

Finite element solution
To compute numerical solutions, we used the finite element method ( 1 elements on triangular elements in two dimensions), implemented in the open software FreeFem++ [14]. For example for the geometry on the left in Figure 3 we use 15 380 triangles in and 15664 triangles in .

Reference signal
In order to solve the Bloch-Torrey problem (2.7), we will consider the case of constant and recall that we take Ω = R . Following the approach of [24], we transform the Bloch-Torrey equation by defining a new unknowñ︁ almost everywhere on R ×]0, [ bỹ︁ Multiplying the equation of the system (2.7) by · ( ) and using the definition of̃︁ ( , ), we obtain the following transformed PDE on the box: The normalized signal is

Second-order reduced model
To obtain the second order reduced model, we observe that the ODE model (4.11) is Numerically, we solve these four box problems by using the finite element method ( 1 on triangular elements in two dimensions). We solved the variational problems and obtained the tensors. To solve the ODE model we used the -method.

Zeroth-order Reduced model
To obtain the zeroth order reduced model, the ODE model (4.15) only requires )︂ .

Convergence
In Figure 4 we display the logarithm of the computed signals as a function of for the reference signal Ref, the order zero approximation 0 and the order two approximation 2. The zeroth-order signal is a horizontal line because the model (4.15) does not depend on . One clearly observes the improvement of the second-order approximation over the zeroth-order approximation.
In Figure 5 we show the convergence of the approximation errors | Ref − 0 | and | Ref − 2 | as functions of for three different choices of (∆, ). The convergence order is numerically computed to be 1.8 for the zerothorder model and it is numerically computed to be around 3.5 for the new second-order asymptotic model. The convergence orders are consistent in that R2 converges 2 orders faster than R0. The apparence of the half order is typical for diffusion problems.

Improvement for higher permeability
For application to dMRI, we observe that the ODE models (4.15) and (4.11) are independent from the choice of the measuring unit of the unit box . For instance, if one replaces by ′ = where is a scaling factor and by ′ = , then one should change by / . If we denote by , , , and the solution of the box problems computed by replacing with ′ , one can easily check that We will choose a value of = 1 × 10 −2 which allows us to use the periodicity box = of 10 m 2 . This means that we assume that the periodicity is on the order of the sizes of brain cells ( ( m)). Let us emphasize that we did not check that this periodicity length assumption is valid with experimental data. Our comparison will be done with the solution to the Bloch-Torrey equation. Physically, diffusion time in the brain in the range of 1 ms-100 ms can be measured which correspond to average diffusion distance of 2.5 m-25 m. In addition, the permeability coefficient varies among the different tissues from sub 10 −6 m/s to more than 10 −3 m/s. For instance, for values of between 10 −6 m/s and 10 −4 m/s, the membranes are semi-permeable while bigger values of correspond to very permeable membranes. If is below 10 −6 m/s, the membranes can be considered as impermeable.
In our numerical tests, we will vary the diffusion time ∆ and the cell permeability . We see clearly in Figures 6a and 6b that the signal of 2 is closer to the reference signal when we use a long diffusion time ∆. In Figures 6c and 6d, we consider the effect of the membrane permeability. We can see with varing that at higher permeability, 2 is closer to Ref.
Two interesting quantities in dMRI are the ADC (the apparent diffusion coefficient) and the AK (the apparent kurtosis). According to [16], the ADC and the AK can be fitted from the signal using In Figure 7, we see that the ADC and AK of the 2 model are closer to the reference values than the 0 model. We show in Figure 8 that in a configuration of 4 biological cells, the 2 model is more accurate than the 0 model for a wide range of -values and especially at higher permeability. Therefore, to have stability, both eigenvalues of ( , ) must be non-negative. In Figure 9, we plot the two eigenvalues as a function of time for two values of and different values of and ∆. According to Figure 9, we can see that both eigenvalues are always non-negative which shows numerical stability.

Conclusions
In this paper, we formulated a second order asymptotic model to treat higher values of permeability. We explicitly solved this new asymptotic model to obtain a system of ordinary differential equations that can model the diffusion MRI signal and we present numerical results showing the improved accuracy of the new model in the regime of higher permeability.
Higher-order homogenization expansions have been employed in many other contexts in the literature and we made ample use of higher dimensional tensors in the development of our asymptotic model. We hope our work will offer insight into the role higher dimensional tensors play in multiple scale diffusion phenomena.

Appendix A.
We here give some details of the calculations that lead to the macroscopic equations for̂︁ 0 ,̂︁ 1 and̂︁ 2 as well as some symmetry properties for the homogenized tensors. The convention of summation on repeated indices is adopted in the following expressions. Using the differentiation rules we observe that Substituting by the asymptotic expansion (3.1) in problem (2.7) and multiplying the first two equations by 2 , we then obtain (3.5)-(3.7).
To derive the averaged equations, we need the following theorem, which can be proved using the Green's identities. By equating the zeroth-order terms in we get It is clear that 0 does not depend on and therefore By equating the first-order terms in we get The compatibility condition (A.1) is automatically verified which guarantees a unique (up to an additive constant) solution 1 that depends on ∇̂︁ 0 . Then, by linearity we obtain wherê︁ 1 = ⟨ 1 ⟩ and for = 1, . . . , are solutions of box problems (3.18). Next, by equating the second-order terms in we get We use the notatioñ︀̂︁ It is easy to see that the compatibility condition (A.1) for problem (A.5) implies Replacing 1 by its expression (A.4) we obtain the macroscopic model (3.16) satisfied bŷ︁ 0 .

A.2. Macroscopic model for̂︁ 1
We here would like to determine the correction given by first-order terms. To derive the macroscopic model for̂︁ 1 , we need to get the expression of 2 . By inserting the expression of 1 (A.4) into (A.5), we can look for a solution 2 such that So, we can express the solution of (A.13) as is the solution of By equating the third-order terms in we get Let us define Then, the compatibility condition for (A.18) implies (A.20) By replacing 1 and 2 by their expressions (A.4) and (A.17) and after recollection we obtain Then,̃︀̂︁ Since H is antisymmetric in its last two indices, it follows that div (H : ∇ 2̂︁ 0 ) = 0. We introduce the vector through

A.3. Macroscopic model for̂︁ 2
To derive the macroscopic model for̂︁ 2 in this case of scaling, we need to have an expression for 3 in terms of̂︁ 0 ,̂︁ 1 and̂︁ 2 . So we consider the problem (A.18) and replace 1 and 2 by their expressions (A.4) and (A.17). We look for a solution 3 such that which can be expressed as where is the solution of the cell problem (3.18). Since we have ⟨ ⟩ = 0,̂︁ ,1  is the solution of is -periodic, (A. 28) which can be expressed as The third-order tensor ℛ is defined as the solution to: where and are solutions to (3.21) and (3.22), respectively.
Therefore, the expression of 3 is written as (A.42). To derive the macroscopic model for̂︁ 2 , we will look for the expression of 3 . Therefore, we can express the solution to (A.18) as (A.42) By equating the fourth-order terms in we get The compatibility condition for this problem implies and After some simplifications, we find div We denote by ℳ ,2 and ℳ ,2 the Fourier transform with respect to the space variable of respectivelỹ︁ ,2 and̃︁ ,2 and denote the dual variable by . We observe that for = 0, the system of equations for ℳ ,2 and ℳ ,2 directly implies the ODE model (4.3). Finally, for the signal we have for = , , using ( ) = 0,