Houssen Haddar

Marwa Kchaou

Maher Moakher

The Derivation of Homogenized Diffusion Kurtosis Models for Diffusion MRI

We use homogenization theory to establish a new macroscopic model for the complex transverse water proton magnetization in a voxel due to diffusion-encoding magnetic field gradient pulses in the case of biological tissue with impermeable membranes. In this model, new higher-order diffusion tensors emerge and offer more information about the structure of the biological tissues. We explicitly solve the macroscopic model to obtain an ordinary differential equation for the diffusion MRI signal that has similar structure as diffusional kurtosis imaging models. We finally present some validating numerical results on synthetic examples showing the accuracy of the model with respect to signals obtained by solving the Bloch-Torrey equation.

Introduction

Diffusion Magnetic Resonance Imaging (dMRI) is a non-invasive imaging technique that is applied in material sciences to investigate structural and transport properties of porous media as well as in medicine and neuroscience to study anatomical, physical, and functional properties of biological tissues and organs such as brain, muscles and bones [START_REF] Chen | Regional ventricular wall thickening reflects changes in cardiac fiber and sheet structure during contraction: quantification with diffusion tensor MRI[END_REF][START_REF] Grebenkov | From the microstructure to diffusion MRI, and back[END_REF][START_REF] Le Bihan | Diffusion MRI at 25: Exploring brain tissue structure and function[END_REF][START_REF] Moseley | Diffusion-weighted MR imaging of acute stroke: correlation with T2-weighted and magnetic susceptibility-enhanced MR imaging in cats[END_REF][START_REF] Rohmer | Reconstruction and visualization of fiber and sheet structure with regularized tensor diffusion MRI in the human heart[END_REF][START_REF] Warach | Fast magnetic resonance diffusionweighted imaging of acute human stroke[END_REF]. The classical dMRI experiment consists in applying two pulsed-gradient magnetic fields with a 180 degree spin reversal between the two pulses in order to encode the displacement of water molecules between them [START_REF] Callaghan | Diffraction-like effects in nmr diffusion studies of fluids in porous solids[END_REF][START_REF] Callaghan | Frequency-domain analysis of spin motion using modulatedgradient NMR[END_REF][START_REF] Mori | Introduction to Diffusion Tensor Imaging[END_REF][START_REF] Price | Pulsed-field gradient nuclear magnetic resonance as a tool for studying translational diffusion: Part 1. basic theory[END_REF].

The aim of this paper is to propose a new justification of the use of higher order diffusion tensor in the modeling of Diffusion Tensor Imaging (DTI) (See for instance [START_REF] Li | Numerical study of a macroscopic finite pulse model of the diffusion MRI signal[END_REF][START_REF] Liu | Characterizing non-Gaussian diffusion by using generalized diffusion tensors[END_REF][START_REF] Peerlings | Higher-order asymptotic homogenization of periodic linear elastic composite materials at low scale separation[END_REF] and the references therein), through homogenization theory [START_REF] Bensoussan | Asymptotic Analysis for Periodic Structures[END_REF]. Following the approach proposed in [START_REF] Coatléven | A new macroscopic model including membrane exchange for diffusion MRI[END_REF][START_REF] Haddar | A macroscopic model for the diffusion MRI signal accounting for time-dependent diffusivity[END_REF], we consider the periodic Bloch-Torrey partial differential equation as a microscopic model for the (complex) transverse water proton magnetization diffusion due to magnetic field gradient pulses. The goal is to determine higher-order asymptotic expansion of the solution in terms of the ratio cell-size/voxel-size. The latter is assumed to be small compared to the size of a voxel in a dMRI image. The obtained model can be seen as an improvement of the one derived in [START_REF] Coatléven | A new macroscopic model including membrane exchange for diffusion MRI[END_REF][START_REF] Li | Numerical study of a macroscopic finite pulse model of the diffusion MRI signal[END_REF]. We remark that the homogenization model obtained in the latter is time-independent and therefore is physically relevant only for the long-time limit of the dMRI signal. Going to higher-order terms in the homogenization process, we end up with models that depend on diffusion time and therefore would be valid for larger time intervals.

We mention here that there are other approaches to homogenization than the periodic homogenization theory used in this study. For example, in the case of random coefficients, effective medium theory can be used to obtain averaged quantities based on a statistical approach (see for instance [START_REF] Novikov | Effective medium theory of a diffusion-weighted signal[END_REF] and the references therein). The latter is more general as it does not need the periodicity assumption. However, it is based on physically oriented modeling ansatz and leads to computable expressions only in specific regimes for time profile.

Higher-order homogenization expansions has been employed in many other contexts in the literature and we refer the reader to [START_REF] Allaire | A comparison between two-scale asymptotic expansions and Bloch wave expansions for the homogenization of periodic structures[END_REF][START_REF] Allaire | Optimization of dispersive coefficients in the homogenization of the wave equation in periodic structures[END_REF][START_REF] Conca | On Burnett coefficients in periodic media[END_REF][START_REF] Moskow | First-order corrections to the homogenized eigenvalues of a periodic composite medium. aconvergence proof[END_REF][START_REF] Santosa | A dispersive effective medium for wave propagation in periodic composites[END_REF] and the references therein for an account of the vast literature on the subject. Our main contribution is to exploit this concept to justify the so-called diffusion kurtosis imaging (DKI) proposed in the literature of dMRI (See for instance [START_REF] Andrianov | Higher order asymptotic homogenization and wave propagation in periodic composite materials[END_REF][START_REF] Hui | A high order homogenization model for transient dynamics of heterogeneous media including micro-inertia effects[END_REF][START_REF] Jensen | Diffusional kurtosis imaging: the quantification of non-Gaussian water diffusion by means of magnetic resonance imaging[END_REF][START_REF] Peerlings | Higher-order asymptotic homogenization of periodic linear elastic composite materials at low scale separation[END_REF]. This technique is used to directly measure the non-Gaussian property of water diffusion characterized by a fourth-order tensor referred to as diffusion kurtosis tensor. In three-dimensional space, the kurtosis tensor has 81 components. But, owing to the full symmetry with respect to an interchange of indices, only 15 elements are independent. The DKI method provides a complementary information to DTI that can be useful in diagnosis [START_REF] Jensen | Diffusional kurtosis imaging: the quantification of non-Gaussian water diffusion by means of magnetic resonance imaging[END_REF][START_REF] Lu | Three dimensional characterization of non Gaussian water diffusion in humans using diffusion kurtosis imaging[END_REF]. It is based on the same type of pulse sequences employed for conventional diffusion-weighted imaging (DWI), but the required b values are somewhat larger than those usually used to measure diffusion coefficients. Generalized diffusion tensors (GDTI) have been also employed in [START_REF] Ozarslan | Generalized diffusion tensor imaging and analytical relationships between diffusion tensor imaging and high angular resolution diffusion imaging[END_REF] in order to assess high angular resolution diffusion imaging (HARDI) data. Among the differences between GDTI and DKI is that the generalized diffusion tensors defined in [START_REF] Ozarslan | Generalized diffusion tensor imaging and analytical relationships between diffusion tensor imaging and high angular resolution diffusion imaging[END_REF] depend on the choice of b-value while the one in [START_REF] Peerlings | Higher-order asymptotic homogenization of periodic linear elastic composite materials at low scale separation[END_REF][START_REF] Lu | Three dimensional characterization of non Gaussian water diffusion in humans using diffusion kurtosis imaging[END_REF] does not.

We shall see that our modeling lead to a DKI that depends on the time profile characteristics and therefore makes the bridge between the two mentioned above models. Moreover, we obtain computable expressions of DKI in terms of the cell geometry and intrinsic diffusion coefficients. In addition to homogenization technique, we rely on the interpretation of the macroscopic model in the Fourier domain to derive an ordinary differential equation (ODE) model that directly implies the dMRI signal over a given voxel (as for instance in [START_REF] Coatléven | A new macroscopic model including membrane exchange for diffusion MRI[END_REF][START_REF] Haddar | Adapting the Kärger model to account for finite diffusion-encoding pulses in diffusion MRI[END_REF]). In these first investigations, we restrict ourselves to the cases where the exchange across cell membranes is negligible, which is a good approximation for white matter or for relatively small b-values. The DKI expression also applies to heterogeneous intrinsic diffusion coefficients. We validate the accuracy of our models through some two-dimensional numerical simulations. We show in particular that the new model provides much better accuracy for the signal produced by simulating the full Bloch Torrey model.

Our paper is organized as follows. We present the classical Bloch-Torry model and fix notation in Section 2. We outline the periodic homogenization technique in Section 3 and provide the key steps to compute the different terms of the two-scale asymptotic expansions. We provide explicit expressions up to order 2. We present in Section 4 the ODE model that can be derived from the homogenized equations where we show that the signal from ODE model approximates the signal of the microscopic model at long diffusion times. We finally give in Section 5 some validating numerical results for toy models where we show that the signal from the homogenized model approximates the signal of the microscopic model at long diffusion times. For the reader convenience, we complement our paper with an appendix containing some technical details associated with Section 3.

Bloch-Torrey equations

For a volume Ω ⊂ R d (d = 2 or 3) occupied by a biological tissue, we denote by Γ the union of the boundaries of all the biological cells included in Ω. The domain Ω ext then represents the union of two subdomains the extra-cellular, Ω e , and intra-cellular, Ω c , domains

Ω ext = Ω \ Γ = Ω e ∪ Ω c .
The complex transverse water proton magnetization M = M x + ιM y , is modeled by the Bloch-Torrey equation [START_REF] Grebenkov | From the microstructure to diffusion MRI, and back[END_REF][START_REF] Torrey | Bloch equations with diffusion terms[END_REF]:

∂ ∂t M (x, t) + ιf (t)q • x M (x, t) -div(D(x)∇ M (x, t))+ 1 T 2 M (x, t) = 0, in Ω ext ×]0, T [, (1)
where T 2 is the spin-spin relaxation time, T > 0 is the final time, ι in the imaginary unit and D is the intrinsic diffusion function:

D = D e , in Ω e , D c , in Ω c . (2)
The constant vector, q ∈ R d , contains the amplitude and direction of the applied diffusion-encoding magnetic field gradient multiplied by the gyro-magnetic ratio of water protons, γ = 2.67513 × 10 8 rads -1 T -1 .

If we assume that the relaxation time T 2 is a constant independent of spatial variable x, then, by representing the transverse magnetisation M by a product of a purely time-dependent function exp(-1 T2 t) and a function of time and space M (x, t), i.e., M (x, t) = exp(-1

T 2 t)M (x, t),
the Bloch-Torrey equation (1) becomes

∂ ∂t M (x, t) + ιf (t)q • xM (x, t) -div(D(x)∇M (x, t)) = 0, in Ω ext ×]0, T [. (3)
Although that, in a very recent study, Veraart et al. [START_REF] Veraart | TE dependent diffusion imaging (TEdDI) distinguishes between compartmental T2 relaxation times[END_REF] showed that the relaxation time T 2 depends on the spatial position and hence it might not be valid to assume that T 2 is spatially uniform, in this article we make this assymption. The time profile, f (t) is chosen to satisfy the refocusing (or rephasing) condition

T E 0 f (t)dt = 0,
where T E is the time at which the signal is measured and is called the echo time. In this work, we consider the case where there is no water exchange between compartments, i.e., we consider the boundary condition

D(x)∇M (x, t) • ν| Γ = 0, on Γ, (4)
where ν is the outward unit normal to the domain Ω c . As a matter of fact, the diffusion inside and outside the cells are decoupled and one can only consider the diffusion in one domain (the other one is obtained in the same way). The PDE (3) and the boundary condition (4) must be supplemented with the initial condition:

M (x, 0) = M init , in Ω ext . (5)
The dMRI signal at echo time T E is the integral of magnetization M (x, T E ), normalized by the mass:

S := 1 |Ω| Ω M (x, T E ) dx,
where M is the solution of equation ( 3) and |Ω| denotes the volume of Ω. Numerically, the signal is usually plotted against a quantity called the b-value b(|q|

) := |q| 2 T E 0 t 0 f (s) ds 2 dt.
The homogenization method is a mathematical theory of averaging. It is an approach that allows to infer variations at macroscopic scales (pixel-level) from microscopic scale (cellular-level) description of the PDE model. For simplicity, we shall employ periodic homogenization theory [START_REF] Allaire | Optimization of dispersive coefficients in the homogenization of the wave equation in periodic structures[END_REF][START_REF] Bensoussan | Asymptotic Analysis for Periodic Structures[END_REF]. We assume that the volume to be modeled, Ω, can be described at the microscopic scale as a periodic domain. The periodicity scale ε, related to the size of cell versus the size of the voxel, is assumed to be small. We define the normalized periodicity box

Y =]0, 1[ d such that Y = Y e ∪ Y c
and Γ m = ∂Y c = ∂Y e \ ∂Y . As indicated earlier, due to the choice of the boundary conditions at the interface, applying the homogenization process to the exterior or the interior domain is similar. However, if the intra-cellular domain does not communicate with the extra-cellular domain and has a geometry Y e strictly included inside Y , then the homogenized tensor associated with this compartment vanishes. Therefore, in what follow we concentrate on applying the homogenization process to one subdomain that is the exterior domain. Taking into account the periodicity at microscale we define

Ξ ε = {ξ ∈ Z d |ε(ξ + Y ) ∩ Ω = ∅}, Ω ε e = ∪ ξ∈Ξ {ε(ξ + Y e ) ∩ Ω}.
Therefore, the initial boundary-value problem (3), ( 4) and ( 5) can be rewritten as:

       ∂ ∂t M ε (x, t) + ιf (t)Q(x)M ε (x, t) -div(D ε (x)∇M ε (x, t)) = 0, in Ω ε e ×]0, T [, D ε (x)∇M ε (x, t) • ν = 0, on Ξ ε , M ε (•, 0) = M init , in Ω ε e , (6)
where

D ε (x) := D( x ε ) and Q(x) = q • x.
We assume that the initial data M init defined on Ω is independent of ε.

For the choice of the volume to be modeled we shall take Ω = R d to avoid complications associated with imposing boundary conditions on Ω. This is justified by the fact that the voxel of interest is itself surrounded by other voxels and diffusion time is such that diffusion length is smaller than the voxel dimensions. This assumption avoids the modeling of boundary layer effects and justifies the use of higher-order asymptotics in the periodic homogenization ansatz.

A macroscopic model through two-scale asymptotic expansion

The present section is dedicated to an outline of the periodic homogenization machinery for deriving macroscopic models. We are interested in models that involve higher-order diffusion tensors, which requires going at least up to the second order in the asymptotic expansion. This type of calculation is by now very classical in the mathematical literature and its justification is well understood in the setting where no macroscopic boundary conditions are involved. The specificity of the considered model with respect to the literature is the presence of the term ιf (t)Q(x)M ε which induces new homogenized tensors for the second-order model as compared to the classical ones (see for instance [START_REF] Coatléven | A new macroscopic model including membrane exchange for diffusion MRI[END_REF]).

Following the periodic homogenization [START_REF] Bensoussan | Asymptotic Analysis for Periodic Structures[END_REF][START_REF] Santosa | A dispersive effective medium for wave propagation in periodic composites[END_REF], we postulate a two-scale asymptotic expansion for M ε in Ω ε e as

M ε (x, t) = +∞ n=0 ε n M n x, x ε , t in Ω ε e , (7)
where the function M n (x, y, t) is defined on Ω × Y e ×]0, T [ and is Y -periodic in y (the variable x indicates macroscopic variations and the variable y = x ε indicates the microscopic ones). Our goal is to characterize term by term the functions M n up to n = 2 and then write a model for the truncated sum (up to n = 2) of the average with respect to the y variable. This would give us (a formally) second-order approximation of M ε with respect to ε. For second-order PDE, in order to obtain the macroscopic equations for n = 2 we need to write the microscopic equations that involve the terms up to the order n + 2.

In the following we shall use the notation

M n (x, t) = M n (x, •, t) Ye = 1 |Y e | Ye M n (x, y, t)dy,
Also, 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 C and P

(C : P) i,l = j,k C iljk P iljk , (C ∴ P) i = l,j,k
C iljk P iljk and (C ::

P) = i,l,j,k C iljk P iljk .
The tensor product of two vectors u and v is defined by

(u ⊗ v) ij = u i v j .
3.1 Macroscopic model for M 0 , M 1 and M 2

We here summarize the equations satisfied by M 0 , M 1 and M 2 which are obtained after substituting ( 7) into ( 6) and equating the same powers of ε. The formal procedure is detailed in Sections A.1-A.3 of the Appendix. These equations necessitate the calculation of some macroscopic quantities through the solutions of time-independent problems posed on the unit box Y . We first introduce the solutions to the following classical unit-box problems

     -div y (D∇ y ω i + σe i ) = 0, in Y e , D∇ y ω i • ν + De i • ν = 0, on Γ m , ω i is Y -periodic, ω i Ye = 0 ( 8
)
where e i , i = 1, . . . , d are the vectors of the canonical basis of R d . These solutions allow us to define the homogenized second-order tensor D by

D ij := D∇ω j • e i + Dδ ij Ye , i, j = 1, . . . , d.
We observe that D is symmetric as its components can be written as D ij = D(∇ω j + e j ) • (∇ω i + e i ) Ye . With the aid of this tensor we can write the homogenized equations for M 0 as

   ∂ ∂t M 0 + ιQf (t) M 0 -div x (D∇ x M 0 ) = 0, in Ω×]0, T [, M 0 (•, 0) = M init , in Ω. (9)
It turns out that, likewise, the macroscopic model for M 1 can also be written as

   ∂ ∂t M 1 + ιQf (t) M 1 -div x (D∇ x M 1 ) = 0, in Ω×]0, T [, M 1 (•, 0) = 0, in Ω. (10)
The macroscopic model for M 2 involve two additional homogenized tensors and can be written as

       ∂ ∂t M 2 + ιQf (t) M 2 -div x (D∇ x M 2 ) = div x (T ∴ ∇ 3 x M 0 ) -div x D∇ x (S : ∇ 2 x M 0 ) + ιf (t)S : q ⊗ ∇ x M 0 , in Ω×]0, T [, M 2 (•, 0) = 0, in Ω, (11)
where the fourth-order tensor T is defined by

T ijkl = D ∂ ∂y l (Q ijk + R ijk ) Ye , i, j, k, l = 1, . . . , d,
and the symmetric second-order tensor S is defined by

S ij = D ∂ ∂y j θ i Ye = -ω i • ω j Ye .
The functions θ i , Q ijk and R ijk for i, j, k = 1, . . . , d are respectively defined as solutions to the following problems

     -div y (D∇ y θ i ) = ω i , in Y e , D∇ y θ i • ν = 0, on Γ m , θ i is Y -periodic, θ i Ye = 0. ( 12
)      -div y (D∇ y Q ijk ) = Dω j δ ki + D ∂ ∂y k (χ ij + ϕ ij ) -H ijk , in Y e , in Y e , D∇ y Q ijk • ν = 0, on Γ m Q ijk is Y -periodic, Q ijk Ye = 0 (13)      -div y (D∇ y R ijk ) = ∂ ∂y k (D(χ ij + ϕ ij )) , in Y e , D∇ y R ijk • ν + D(χ ij + Dϕ ij )ν k = 0, on Γ m , R ijk is Y -periodic, R ijk Ye = 0. ( 14
)
These problems make use of the functions χ ij and ϕ ij i, j = 1, . . . , d that are defined as solutions to the problems

       -div y (D∇ y χ ij ) = ∂ ∂y j (Dω i ), in Y e , D∇ y χ ij • ν + Dω i ν j = 0, on Γ m , χ ij is Y -periodic, χ ij Ye = 0, (15)      -div y (D∇ y ϕ ij ) = D ∂ ∂yi ω j + Dδ ij -D ij , in Y e , D∇ y ϕ ij • ν = 0, on Γ m , ϕ ij is Y -periodic, ϕ ij Ye = 0. ( 16
)
We refer to Section A.4 in the appendix for the expression of the third order tensor H.

The second-order macroscopic model

Using the three macroscopic models for M 0 , M 1 and M 2 given in ( 9), ( 10) and ( 11), we obtain the following second-order model that is satisfied, up to O(ε 3 )

       ∂ ∂t M ε,2 -div x D∇ x (1 -ε 2 S : ∇ 2 x )M ε,2 -ε 2 div x (T ∴ ∇ 3 x M ε,2 ) + ιf (t) QM ε,2 -ε 2 S : q ⊗ ∇ x M ε,2 = 0, in Ω×]0, T [, M ε,2 (•, 0) = M init , in Ω. (17)
We then formally expect that [START_REF] Li | Numerical study of a macroscopic finite pulse model of the diffusion MRI signal[END_REF] is in fact unstable for computing M ε,2 since the fourth-order operator is non negative [START_REF] Allaire | A comparison between two-scale asymptotic expansions and Bloch wave expansions for the homogenization of periodic structures[END_REF][START_REF] Conca | On Burnett coefficients in periodic media[END_REF]. However, since we are only interested in the mean value of M ε,2 at echo time, one is able to deduce from this model an ODE for the mean value that is stable for sufficiently small ε (see next section).

M ε = M ε,2 + O(ε 3 ) in Ω ε e . The model
Remark 1. We observe that the homogenized model ( 17) is independent from the choice of the measuring unit (m, mm, µm, etc.) of the unit box Y . For instance, if one replaces Y by Y = ηY where η is a scaling factor, then one should change by ε/η (since the physical periodicity is εY ). If we denote by D η , S η and T η the homogenized tensors computed by replacing Y with Y , then one can easily check that

D η = D, S η = η 2 S, and T η = η 2 T.
This shows that the coefficients in [START_REF] Li | Numerical study of a macroscopic finite pulse model of the diffusion MRI signal[END_REF] remain the same under the scaling η.

Asymptotic dMRI signal model for DKI

In this section, we explain how to use the macroscopic model [START_REF] Li | Numerical study of a macroscopic finite pulse model of the diffusion MRI signal[END_REF] to obtain an ODE model that allows us to compute the Kurtosis tensor. This approach has also been used in [START_REF] Coatléven | A new macroscopic model including membrane exchange for diffusion MRI[END_REF] for DTI and in [START_REF] Haddar | Adapting the Kärger model to account for finite diffusion-encoding pulses in diffusion MRI[END_REF] to obtain kurtosis approximations for small diffusion times. In practice, the reference normalized signal is measured at the echo time T E and it is given by

S BT ε (q) := R d M ε (x, T E ) dx R d M init dx
, where M ε is the solution of the reference Bloch-Torrey equation ( 6). Then, it is natural to define the approximate normalized signal by

S ε (q) := R d M ε,2 (x, T E ) dx R d M init dx
.

We start by transforming the third-order equation ( 17) by introducing the new unknown variable M ε,2 defined almost everywhere on R d ×]0, T [ by

M ε,2 (x, t) = M ε,2 (x, t)e ιqx•nF (t) , (18)
where F (t) := t 0 f (s)ds and n is a unit vector such that q = qn. Since F (T E ) = 0, one observes that M ε,2 (•, T E ) = M ε,2 (•, T E ) and therefore

S ε (q) := R d M ε,2 (x, T E ) dx R d M init dx . Setting m ε (t) := R d M ε,2 (x, T E ) dx,
one obtains after inserting the change of variable ( 18) in ( 17) and integrating over

R d    d dt m ε (t) + F (t) 2 q 2 (D : (n ⊗ n))m ε (t) + ε 2 F (t) 4 q 4 (B :: (n ⊗ n ⊗ n ⊗ n))m ε (t) -ε 2 f (t)F (t)q 2 (S : (n ⊗ n))m ε (t) = 0, (19)
where the fourth-order tensor B, called the Burnett tensor [START_REF] Conca | On Burnett coefficients in periodic media[END_REF], is defined by

B ijkl = D ij S kl -T ijkl . (20)
We see clearly that the ODE model ( 19) can be solved explicitly as

m ε (t) = m ε (0) exp - t 0 L ε (q, s)ds , (21)
where

L ε (q, t) =F (t) 2 q 2 (D : (n ⊗ n)) + ε 2 F (t) 4 q 4 (B :: (n ⊗ n ⊗ n ⊗ n)) -ε 2 f (t)F (t)q 2 (S : (n ⊗ n)).
The approximate normalized signal can then be rewritten as

S ε (q) = m ε (T E ) m ε (0) = exp - T E 0 L ε (q, s)ds .
We observe that

T E 0 f (t)F (t)dt = (F (T E ) 2 -F (0)
2 )/2 = 0 and therefore the expression simplifies to

S ε (q) = exp (-K ε (q, T E )) , (22)
where

K ε (q, T E ) := b(q)(D : (n ⊗ n)) + ε 2 q 4 (B :: (n ⊗ n ⊗ n ⊗ n)) T E 0 F (t) 4 dt. (23)
The expression for -ln ( S ε (q)) resembles the one proposed in [START_REF] Jensen | Diffusional kurtosis imaging: the quantification of non-Gaussian water diffusion by means of magnetic resonance imaging[END_REF] for DKI. The fourth-order tensor W referred to as diffusion kurtosis tensor is defined in [START_REF] Jensen | Diffusional kurtosis imaging: the quantification of non-Gaussian water diffusion by means of magnetic resonance imaging[END_REF] such that

ln[ S ε (b, n)] = -b(q)(D : (n ⊗ n)) + 1 6 b(q) 2 ( 1 3 TrD) 2 (W :: (n ⊗ n ⊗ n ⊗ n)). (24)
Then, according to our model, the fourth-order diffusion kurtosis tensor W and the fourth-order Burnett tensor B are multiple of one another and are related by

W = -6ε 2 T E 0 F (t) 4 dt ( 1 3 TrD) 2 ( T E 0 F (t) 2 dt) 2 B. (25)
In addition to having explicit dependence on the microstructure geometry, we observe that this expression depends also on the time profile f through ∆ and δ. The latter is a notable difference with [START_REF] Jensen | Diffusional kurtosis imaging: the quantification of non-Gaussian water diffusion by means of magnetic resonance imaging[END_REF][START_REF] Lu | Three dimensional characterization of non Gaussian water diffusion in humans using diffusion kurtosis imaging[END_REF][START_REF] Veraart | More accurate estimation of diffusion tensor parameters using diffusion kurtosis imaging[END_REF]. We observe that for fixed q, K ε (q, T E ) is indeed an approximation of b(q)Dn • n as ε → 0. However, this approximation is not uniform with respect to q → ∞. For instance, one can prove that (B ::

(n ⊗ n ⊗ n ⊗ n)) ≤ 0,
exactly in the same way as it is proved in [START_REF] Conca | On Burnett coefficients in periodic media[END_REF] for the case Γ m = ∅. Therefore, K ε (q, T E ) can become negative for large q and fixed ε. A domain of validity can be roughly fixed by the values of q and ε such that K ε (q, T E ) > 0. One can also replace K ε (q, T E ) by the Padé approximation for small ε, namely

Kε (q, T E ) = (b(q)(D : (n ⊗ n))) 2 b(q)(D : (n ⊗ n)) -ε 2 q 4 T E 0 F (t) 4 dt(B :: (n ⊗ n ⊗ n ⊗ n))
.

We observe that K ε (q, T E ) = Kε (q, T E ) + O(ε 4 ). Therefore, while this substitution does not affect the accuracy of the model in terms of ε it leads to a stable signal for all values of ε and q.

Numerical results

Our objective in this numerical section is to measure the accuracy of the second-order model and the improvement it yields over using the zeroth-order model corresponding to ε = 0. Our reference solution will be provided by numerical solution of the full Bloch-Torrey model [START_REF] Callaghan | Frequency-domain analysis of spin motion using modulatedgradient NMR[END_REF]. Since our concern is to measure the accuracy of the model, using K ε or Kε is the same as long as we restrict ourselves with K ε > 0. We shall use the classical Pulsed Gradient Spin Echo (PGSE) sequence, with two rectangular pulses of duration δ, separated by a time interval ∆ -δ, see e.g. [START_REF] Mori | Introduction to Diffusion Tensor Imaging[END_REF]:

f (t) =      1, if 0 < t ≤ δ, -1, if ∆ < t ≤ ∆ + δ; 0, elsewhere, (26)
We assume that the RF pulse is produced at t = δ+∆ 2 . Consequently, we take the echo time

T E = δ + ∆.
For this time profile

b(q) = q 2 T E 0 F (t) 2 dt = q 2 δ 2 (∆ - δ 3 ),
and

T E 0 F (t) 4 dt = δ 4 (∆ - 3 5 δ).
We then deduce that

K ε (q, T E ) = b(q)(D : (n ⊗ n)) + ε 2 b(q) 2 (B :: (n ⊗ n ⊗ n ⊗ n)) (∆ -3 5 δ) (∆ -δ 3 ) 2 .

and the coefficients of the diffusion kurtosis tensor are given by

W ijkl = -6ε 2 (∆ -3 5 δ) (∆ -δ 3 ) 2 ( 1 3 TrD) 2 B ijkl . (27)
Therefore, to obtain the approximate normalized signal, S ε , we need to compute the tensors D and B by solving the box-problems ( 8), ( 13) and [START_REF] Hui | A high order homogenization model for transient dynamics of heterogeneous media including micro-inertia effects[END_REF]. For the numerical tests presented in the following section, the solutions to these problem is numerically evaluated using the finite element solver FreeFem++ [13].

In order to solve the Bloch-Torrey model ( 6), we use the same change of unknown as [START_REF] Liu | Characterizing non-Gaussian diffusion by using generalized diffusion tensors[END_REF], namely

M ε (x, t) = M ε (x, t)e ιqx•nF (t) .
We impose constant initial magnetization M init . Since the coefficients in the equation satisfied by M ε are all periodic, the solution is also periodic and therefore the computational domain can be restricted to one period Y ε = εY (with periodicity conditions on the boundary). The normalized signal is then computed as

S BT ε (q) = Yε M ε (x, T E ) dx |Y ε |M init . ( 28
)
For the numerical simulations, once again we use the finite-element solver FreeFem++ [13] with a Crank-Nicolson scheme for time stepping to evaluate the solution to the full Bloch-Torrey model and the associated normalized signal S BT ε (q). All our numerical examples are restricted to 2D configurations (d = 2).

Examples with impermeable cells

We consider two examples of geometries where we vary the number of impermeable components inside a period. Figure [START_REF] Allaire | A comparison between two-scale asymptotic expansions and Bloch wave expansions for the homogenization of periodic structures[END_REF] shows the configurations that we choose and the finite-element mesh that has been used to compute the solution for the Bloch-Torrey model ( 6) as well as the solutions for the different problems associated with the homogenized tensors. ε , the approximate signal for zero order, S 0 and the approximate signal for order two, S ε , where the values of the model parameters are indicated. The curve -ln( S 0 ) is the horizontal line at b(q)D 2,2 that corresponds to the limiting problem as ε goes to 0. One clearly observes the improvement of the second-order model over this zeroth-order approximation for the three configurations. This indeed validates the expressions derived for the second-order model and attest the relevance of using DKI for improving dMRI imaging.

Examples mimicking membrane effects

Our model also applies to the heterogeneous case where the diffusion tensor is piece-wise constant and there is no barrier. Formally, it corresponds to replacing Y e by Y in equations ( 8) and ( 12)-( 16), removing the boundary conditions on Γ m and considering a piece-wise constant diffusion tensor D defined on the whole box Y . We use this type of configuration to model a three-compartment configuration with extracellular domain, intra-cellular domain and a membrane domain. The membrane domain corresponds to a thin layer around the intra-cellular domain. As test examples, we consider two different configurations that are depicted in Figure 4 as well as associated finite-element mesh that was used to compute the homogenized tensors. The thickness of the layer in the unit box Y (1 mm 2 ) is L 0 = 0.01mm. Outside the circles, the intrinsic diffusion is D e = 3 × 10 -3 mm 2 /s, inside the circles D c = 1.7 × 10 -3 mm 2 /s and inside the membrane layer D layer = 1 × 10 -5 mm 2 /s. for different values of the parameter ∆ = 5 and 15ms. This test shows the same type of results as in the case of impermeable cells. It confirms that the second-order model provides better fit between the "exact signal" and the approximate model.

Conclusion

In this work we have used homogenization theory to derive a higher-order diffusion tensors for dMRI signals. We obtain numerically computable expressions for these tensors in terms of the microscopic geometry and the signal time profile. Our model has similar structure as the one used in diffusion kurtosis imaging. The numerical results on toy configurations with realistic values for the physical parameters show the superiority of the model over the zeroth-order homogenized problem. In a future work, we shall extend this approach to the case of biological cells with membrane barriers where the homogenized model is much more complex and would couple "double exponential models" with diffusion kurtosis imaging.

A Appendix

We here give some details of the calculations that lead to the macroscopic equations for M 0 , M 1 and M 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

∇M n (x, x ε , t) = ∇ x M n (x, x ε , t) + ε -1 ∇ y M n (x, x ε , t), and
div(D∇M n ) = div x (D∇ x n ) + ε -1 div x (D∇ y M n ) + ε -1 div y (D∇ x M n ) + ε -2 div y D∇ y M n ).
Substituting M ε by the asymptotic expansion [START_REF] Chen | Regional ventricular wall thickening reflects changes in cardiac fiber and sheet structure during contraction: quantification with diffusion tensor MRI[END_REF] in problem ( 6) and multiplying the first two equations by ε 2 , we then obtain By equating the zeroth-order terms in ε we get

           +∞ n=0 ε n+2 ∂ ∂t M n + ιQf (t) +∞ n=0 ε n+2 M n - +∞ n=0 ε n+2 div x (D∇ x M n ) - +∞ n=0 ε n+1 (div x (D∇ y M n ) + div y (D∇ x M n )) - +∞ n=0 ε n div y (D∇ y M n ) = 0, that holds for (x, y) in Ω × Y e and      ε 0 D∇ y M 0 • ν + ε 1 D ∇ y M 1 • ν + ∇ x M 0 • ν + ε 2 D ∇ y M 2 • ν + ∇ x M 1 • ν + +∞ n=0 ε n+1 D ∇ y M n+1 • ν + ∇ x M n • ν = 0,
     -div y (D∇ y M 0 ) = 0, in Y e , D∇ y M 0 • ν = 0, on Γ m , M 0 is Y -periodic. ( 31
)
We then get that M 0 does not depend on y and therefore

M 0 = M 0 .
By equating the first-order terms in ε we obtain

     -div y (D∇ y M 1 + D∇ x M 0 ) = 0, in Y e , D (∇ y M 1 + ∇ x M 0 ) • ν = 0, on Γ m , M 1 is Y -periodic. ( 32
)
The compatibility condition [START_REF] Veraart | TE dependent diffusion imaging (TEdDI) distinguishes between compartmental T2 relaxation times[END_REF] is automatically verified which guarantees a unique (up to an additive constant) solution M 1 that depends on ∇ x M 0 . One can express the solution to [START_REF] Warach | Fast magnetic resonance diffusionweighted imaging of acute human stroke[END_REF] as

M 1 = ω • ∇ x M 0 + M 1 , in Y e . (33)
Now, by equating the second-order terms in ε we get

         -div y (D∇ y M 2 + D∇ x M 1 ) = div x (D∇ y M 1 + D∇ x M 0 ) -∂ ∂t M 0 -ιQf (t) M 0 , in Y e , D∇ y M 2 • ν + D∇ x M 1 • ν = 0, on Γ m , M 2 is Y -periodic. ( 34
)
The compatibility condition [START_REF] Veraart | TE dependent diffusion imaging (TEdDI) distinguishes between compartmental T2 relaxation times[END_REF] for this problem implies

Ye div x (D∇ y M 1 + D∇ x M 0 ) - ∂ ∂t M 0 -ιQf (t) M 0 dy = 0.
Substituting M 1 by its expression (33) we obtain the macroscopic problem (9) satisfied by M 0 .

A.2 Macroscopic model for M 1

To derive the macroscopic model for M 1 we need to have an expression of M 2 in terms of M 0 and M 1 by solving (34). To simplify the calculations, problem (34) is divided into well-posed sub-systems such that

M 2 = M 1 2 + M 2 2 + M 3 2 + M 2 ,
where the different terms are defined as follows.

-M 1 2 is the solution of:

     -div y (D∇ y M 1 2 ) = div y (D∇ x M 1 ), in Y e , D∇ y M 1 2 • ν + D∇ x M 1 • ν = 0, on Γ m , M 1
2 is Y -periodic.

Thus, using for i = 1, . . . , d the cell problems [START_REF] Coatléven | A new macroscopic model including membrane exchange for diffusion MRI[END_REF], which are well posed according to [START_REF] Veraart | TE dependent diffusion imaging (TEdDI) distinguishes between compartmental T2 relaxation times[END_REF], one can express the solution as

M 1 2 = ω • ∇ x M 1 + M 1 2 , in Y e .
-M 2 2 is the solution of:

     -div y (D∇ y M 2 2 ) = div y D∇ x (ω • ∇ x M 0 ) , in Y e , D∇ y M 2 2 • ν + D∇ x (ω • ∇ x M 0 ) • ν = 0, on Γ m , M 2
2 is Y -periodic.

For i, j = 1, . . . , d, the cell problems ( 15) are well posed according to the compatibility condition [START_REF] Veraart | TE dependent diffusion imaging (TEdDI) distinguishes between compartmental T2 relaxation times[END_REF]. Therefore, one can express the solution as

M 2 2 = (∇ 2 x M 0 ) ij χ ij + M 2 2 , in Y e .
Observe that, since χ ij Ye = 0, we have

M 2 2 = M 2 2 Ye .
-M 3 2 is the solution of:

     -div y (D∇ y M 3 2 ) = div x (D∇ y M 1 + D∇ x M 0 ) -∂ ∂t M 0 -ιQf (t) M 0 , in Y e , D∇ y M 3 2 • ν = 0, on Γ m , M 3 2 is Y -periodic.
Let us define

F := div x (D∇ y M 1 + D∇ x M 0 ) - ∂ ∂t M 0 -ιQf (t) M 0 .
Then, according to the macroscopic model [START_REF] Conca | On Burnett coefficients in periodic media[END_REF] for M 0 , we remark that

F Ye := div x (D∇ x M 0 ) - ∂ ∂t M 0 -ιQf (t) M 0 = 0.
Then, we have

F = F -F Ye = (∇ 2 x M 0 ) ij D ∂ ∂y i ω j + Dδ ij -D ij .
We consider, for i, j = 1, . . . , d, the cell problems [START_REF] Jensen | Diffusional kurtosis imaging: the quantification of non-Gaussian water diffusion by means of magnetic resonance imaging[END_REF]. The compatibility condition [START_REF] Veraart | TE dependent diffusion imaging (TEdDI) distinguishes between compartmental T2 relaxation times[END_REF] guarantees a unique (up to an additive constant) solution M 3 2 that depends on

∇ x ∇ x M 0 through M 3 2 = (∇ 2 x M 0 ) ij ϕ ij + M 3 2 in Y e .
By regrouping the expressions of M 1 2 , M 2 2 and M 3 2 we obtain the expression for M 2

M 2 = ω • ∇ x M 1 + (∇ 2 x M 0 ) ij (χ ij + ϕ ij ) + M 2 in Y e . (35)
By equating the third-order terms in ε we get

         -div y (D∇ y M 3 + D∇ x M 2 ) = div x (D∇ y M 2 + D∇ x M 1 ) -∂ ∂t M 1 -ιQf (t)M 1 , in Y e , D∇ y M 3 • ν + D∇ x M 2 • ν = 0, on Γ m , M 3 is Y -periodic. ( 36
)
The use of the compatibility condition on M 3 for this problem implies

Ye div x D(∇ y M 2 + ∇ x M 1 ) - ∂ ∂t M 1 -ιQf (t)M 1 dy = 0.
By substituting the expressions for M 1 and M 2 and making some recollections we obtain the following macroscopic model for M 1 .

   ∂ ∂t M 1 + ιQf (t) M 1 -div x (D∇ x M 1 ) = div x (H : ∇ 2 x M 0 ), in Ω×]0, T [, M 1 (•, 0) = 0, in Ω. ( 37
)
where we have defined the homogenized third-order tensor H through

H ijk := D ∂ ∂y k (χ ij + ϕ ij ) Ye , i, j, k = 1, . . . , d.
Since H is antisymmetric with respect to its last two arguments, i.e., H ijk = -H ikj for i, j, k ∈ {1, . . . , d} (see [START_REF] Moskow | First-order corrections to the homogenized eigenvalues of a periodic composite medium. aconvergence proof[END_REF] and also Section A.4), one can easily check that div x (H : ∇ 2 x M 0 ) = 0. Therefore, the macroscopic equations (37) simplify to [START_REF] Grebenkov | From the microstructure to diffusion MRI, and back[END_REF].

A.3 Macroscopic model for M 2

With the same reasoning as above, to derive the macroscopic model for M 2 , we first need to have an expression of M 3 in terms of M 0 , M 1 and M 2 . We consider problem (36) and decompose its solution M 3 as

M 3 = M 1 3 + M 2 3 + M 3 3 + M 4 3 + M 3 ,
where the different terms are defined as follows.

-M 1 3 is the solution of:

     -div y (D∇ y M 1 3 ) = div y (D∇ x M 2 ), in Y c , D∇ y M 1 3 • ν + D∇ x M 2 • ν = 0, on Γ m , M 1
3 is Y -periodic.

Using for i = 1, . . . , d the cell problems (8) which are well posed according [START_REF] Veraart | TE dependent diffusion imaging (TEdDI) distinguishes between compartmental T2 relaxation times[END_REF], one can express the solution as

M 1 3 = ω • ∇ x M 2 + M 1 3 in Y e .
-M 2 3 is the solution of:

     -div y (D∇ y M 2 3 ) = div y (D∇ x M 1 2 ), in Y e , D∇ y M 2 3 • ν + D∇ x M 1 2 • ν = 0, on Γ m , M 2
3 is Y -periodic.

Since we have div

y (D∇ x M 1 2 ) = (∇ 2 x M 1 ) ij ∂ ∂y i (Dω j ),
the solution can be expressed as

M 2 3 = (∇ 2 x M 1 ) ij χ ij + M 2 3 , in Y e .
-M 3 3 is the solution of:

     -div y (D∇ y M 3 3 ) = div y (D∇ x (M 2 2 + M 3 2 )), in Y e , D∇ y M 3 3 • ν + D∇ x (M 2 2 + M 3 2 ) • ν = 0, on Γ m , M 3
3 is Y -periodic.

For i, j, k = 1, . . . , d, we consider the cell problems ( 14) that allow us to express the solution as

M 3 3 = (∇ 3 x M 0 ) ijk R ijk + M 3 3 , in Y e .
-M 4 3 is the solution of:

     -div y (D∇ y M 4 3 ) = div x (D∇ y M 2 + D∇ x M 1 ) -∂ ∂t M 1 -ιQf (t)M 1 , in Y e , D∇ y M 4 3 • ν = 0, on Γ m , M 4
3 is Y -periodic.

Let us define

G := div x (D∇ y M 2 + D∇ x M 1 ) - ∂ ∂t M 1 -ιQf (t)M 1 .
According to the macroscopic model [START_REF] Grebenkov | From the microstructure to diffusion MRI, and back[END_REF] for M 2 , we remark that

G Ye = - ∂ ∂t M 1 -ιQf (t) M 1 + div x (D∇ x M 1 ) = 0.
Then, we have

G =G -G Ye = div x (D∇ y M 2 + D∇ x M 1 ) -div x ( D∇ y M 2 + D∇ x M 1 Ye ) -ω i ∂ ∂t (∇ x M 0 ) i -ιQf (t)ω • ∇ x M 0 .
Using [START_REF] Conca | On Burnett coefficients in periodic media[END_REF] and that ∂ ∂xi (QM 0 ) = Q ∂ ∂xi M 0 + q i M 0 , for each i we get ∂ ∂t (∇ x M 0 ) i + ιQf (t)(∇ x M 0 ) i = (∇ x div x (D∇ x M 0 )) i -ιq i f (t) M 0 .

Substituting the expressions for M 1 and M 2 and making some recollection, we obtain

G =(∇ 2 x M 1 ) ij D ∂ ∂y i ω j + Dδ ij -D ij + (∇ 3 x M 0 ) ijk Dω k δ ij + D ∂ ∂y i (χ jk + ϕ jk ) -H ijk -ω • ∇ x div x (D∇ x M 0 ) + ιω • qf (t) M 0 .
We finally can express the solution as

M 4 3 = (∇ 2 x M 1 ) : ϕ + (∇ 3 x M 0 ) ∴ Q -θ • ∇ x div x (D∇ x M 0 ) + ιθ • qf (t) M 0 + M 4 3 , in Y e ,
where M 4 3 = M 4 3 Ye and, for i, j, k = 1, . . . , d, θ i and Q ijk are solutions of cell problems ( 12) and (13), respectively. Collecting all expressions obtained for M i 3 one gets the expression of M 3

M 3 =ω • ∇ x M 2 + (∇ 2 x M 1 ) : (χ + ϕ) + (∇ 3 x M 0 ) ∴ (R + Q) -θ • ∇ x (div x (D∇ x M 0 )) + ιθ • qf (t) M 0 + M 3 , in Y e . (38)
The next step is to equate the fourth-order terms in ε which gives. Substituting the expressions for M 2 and M 3 and after simplification with the use of ω i Ye = χ ij Ye = ϕ ij Ye = 0, we obtain the macroscopic model ( 11) for M 2 .

A.4 Symmetry properties of H By using ( 8), ( 15), ( 16) and the fact that ω i Ye = 0, the components of the homogenized thirdorder tensor H can be rewritten as Since ∂ ∂yi (σω j ) = D∇ω j • e i = div(Dω j e i ), using the divergence theorem for the last term we finally get H ijk = -D∇ω j • e i ω k + Dω j e i • ∇ω k Ye .

H ijk = 1 | Y e | Ye D∇(χ ij + ϕ ij ) • e k dy = - 1 | Y e | Ye D∇ω k • ∇(χ ij + ϕ ij ) dy = 1 | Y e | Ye div(D∇(χ ij + ϕ ij ))ω k dy - 1 | Y e | Γm D∇(χ ij + ϕ ij ) • νω k ds = - 2 | Y e | Ye
This expression clearly shows the anti-symmetry property H ijk = -H ikj which when combined with the symmetry of ∇ 2 M 0 yields div x (H : ∇ 2 x M 0 ) = 0.

A.5 Alternative expression for T

The expression for the fourth-order tensor can be further simplified. Using (8), ( 13) and ( 14) we have Applying the divergence theorem for the surface integral we end up with T ijkl = -D∇(χ jk + ϕ jk )e i ω l + D(χ jk + ϕ jk )e i ∇ω l -Dω k ω l δ ij Ye .

T ijkl = 1 | Y e | Ye D∇(Q ijk + R ijk ) • e l dy = - 1 | Y e | Ye D∇ω l • ∇(Q ijk + R ijk ) dy = - 2 | Y e | Ye

Fig. 1 .

 1 Fig. 1. Periodicity box Y of 1 mm 2 and associated finite-element mesh for the two different configurations of impermeable cells.

Figure 2

 2 Figure2displays, for each geometry in Figure1, -ln(S) according to ε for three signals the reference signal, S BT ε , the approximate signal for zero order, S 0 and the approximate signal for order two, S ε , where the values of the model parameters are indicated. The curve -ln( S 0 ) is the horizontal line at b(q)D 2,2 that corresponds to the limiting problem as ε goes to 0. One clearly observes the improvement of the second-order model over this zeroth-order approximation for the three configurations. This indeed validates the expressions derived for the second-order model and attest the relevance of using DKI for improving dMRI imaging.Figure3displays the signal in terms of the b-value for fixed value of ε and for the geometry configuration indicated in Figure1(b). We show the results obtained for different values of the parameter ∆ = 2.5, 5ms. These examples show that our approximate kurtosis model has the correct dependence on the parameter ∆ since we obtain a better approximation of the signal for each case.

Figure 3

 3 Figure2displays, for each geometry in Figure1, -ln(S) according to ε for three signals the reference signal, S BT ε , the approximate signal for zero order, S 0 and the approximate signal for order two, S ε , where the values of the model parameters are indicated. The curve -ln( S 0 ) is the horizontal line at b(q)D 2,2 that corresponds to the limiting problem as ε goes to 0. One clearly observes the improvement of the second-order model over this zeroth-order approximation for the three configurations. This indeed validates the expressions derived for the second-order model and attest the relevance of using DKI for improving dMRI imaging.Figure3displays the signal in terms of the b-value for fixed value of ε and for the geometry configuration indicated in Figure1(b). We show the results obtained for different values of the parameter ∆ = 2.5, 5ms. These examples show that our approximate kurtosis model has the correct dependence on the parameter ∆ since we obtain a better approximation of the signal for each case.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Logarithm of normalized signals for two different geometries (see Figure 1) versus ε. The experiment parameters are De = 3 × 10 -3 mm 2 /s, b = 2500s/mm 2 , n = ey, ∆ = 5ms and δ = 2.5ms.

Fig. 4 .

 4 Fig. 4. Periodicity unit box Y (1 mm 2 ) and associated finite-element mesh for the two different configurations of permeable cells.

Figure 5

 5 Figure 5 displays, for each geometry in Figure 4, the signal according to b-value for the three signals S 0 , S ε and S BT ε

Fig. 5 .

 5 Fig. 5. Normalized signals versus the b-value for two different geometries correspond with Figure 4. Two different choices of the parameter ∆ = 5ms and 15ms. The experiment other parameters are De = 3 × 10 -3 mm 2 /s, Dc = 1.7 × 10 -3 mm 2 /s, D layer = 1 × 10 -5 mm 2 /s, L0 = 0.01mm, ε = 6 × 10 -3 , n = ey and δ = 2.5ms.

- 30 )A. 1

 301 that holds for (x, y) in Ω × Γ m and finally the initial conditionsM 0 (x, y, 0) = M init (x) and M n (x, y, 0) = 0 for n > 0.The derivation of averaged equations relies on writing the compatibility conditions for problems of div y (D(y)∇ y u(y)) = g(y), in Y e , D(y)∇ y u(y)• ν = h(y), on Γ m , u(y) is Y -periodic. (29)for some regular data g and h. If problem (29) admits a solution (up to an additive constant) then Ye g(y) dy + Γm h(y) ds(y) = 0. (Macroscopic model for M 0

-

  div y (D∇ y M 4 + D∇ x M 3 ) = div x (D∇ y M 3 + D∇ x M 2 ) -∂ ∂t M 2 -ιQf (t)M 2 , in Y e , D∇ y M 4 • ν + D∇ x M 3 • ν = 0, on Γ m , M 4 is Y -periodic. (39)Writing the compatibility condition[START_REF] Veraart | TE dependent diffusion imaging (TEdDI) distinguishes between compartmental T2 relaxation times[END_REF] for this problemYe div x D(∇ y M 3 + ∇ x M 2 ) -∂ ∂t M 2 -ιQf (t)M 2 dy = 0.

  ∂ ∂y i (Dω j )ω k dy + Γm Dω j • ν i ω k ds.

1 |- 1 | 1 |

 111 χ jk + ϕ jk ))ω l dy -1 | Y e | Ye Dω k δ ij ω l dy + Y e | Γm D(χ jk + ϕ jk )ν i ω l ds = Y e | Ye D∇(χ jk + ϕ jk ) • e i ω l dy -1 | Y e | Ye div(D(χ jk + ϕ jk )e i )ω l dy -Y e | Ye Dω k ω l δ ij dy + 1 | Y e | Γm D(χ jk + ϕ jk )ν i ω l ds.


S Laqrichi

F Marmier

D Gourc

Software cost and duration estimation based on distributed project data: A general framework

Keywords: neural network, regression, duration estimation, cost estimation, comparison

Effort estimation is one of the most challenging tasks in the process of software project management. Enhancing the accuracy of effort estimation remains a serious problem for software professionals. Accurate estimation is difficult to achieve. The main difficulty is to collect distributed knowledge as data and information are often dispersed over different services, departments or organisations. Other main difficulty is to propose a model representative enough of this multi-partner behaviour. The objective of this study is to propose a general framework of the estimation starting from the analysis of the available projects database, the choice and establishment of estimation model, up to the use of this model to make estimation for new projects. In this paper, a comparative study between regression models and neural network models is performed. The proposed study is applied on a dataset of an automotive company.

Introduction

Effort estimation is an important activity in software project management. Estimation in software projects consists in predict likely amount of effort, time and staffing level that are required to build a software system. It is used in the whole development life cycle of the software project from the bidding until the maintenance of the software. Both project managers and clients use effort estimation to predict the effort, the duration and the cost required to develop their software projects in order to establish contracts. Under estimating the effort and the cost required to develop a project results in budget overruns, while over estimating can lead to miss of biddings. Accurate estimation is then very important for companies' benefit and success.

Estimation is a complex activity that requires a high level of interoperability in both steps of model establishment and estimation for new projects. Indeed, the modeling step needs various data about previous projects that are dispersed over different services and entities (financial data, technical project data…). The estimation step requires various distributed data about new projects to estimate as well. Thus, entities have to communicate and collaborate to provide required information.

Various effort and duration estimation approaches were been developed. Traditional and well-known ones include expert judgment, Delphi, COCOMO and Putnam's SLIM. However, estimation methods did not produce sufficiently accurate results, this is why approximately 44% of software projects according to the Standish Group International fail on meeting the commitment on quality, time and cost.

A set of factors can influence the estimation accuracy and lead to estimations far from the reality. It includes, among others: lack of information about completed previous projects; use of new technologies; lack of experience with similar projects; choice of estimation approach and more [START_REF] Heemstra | Software cost estimation[END_REF].

The challenge of improving estimates accuracy has led to the development of several new methods and techniques for effort, duration and cost estimation. These methods are based on artificial intelligence such as NN (Neural Network) models.

Our work first focuses in formalizing the general estimation framework. This framework enables to compare different models. In this paper, we study the case of regression and NN estimation models applied on a big and diversified case study. This case study does not contain size project that is usually considered to be an important cost driver in estimation model establishment.

The present paper is organized into three sections: the first section presents literature review on effort estimation process, regression and NN models and their comparisons work carried by researchers. In the second section, a general framework for estimation is proposed. Finally, in the third section, the proposed framework is applied to a case study from the automotive industry.

Literature review 1.2.1 Estimation in software projects

The estimation process is based on two principal activities that are: (i) project size measure and (ii) effort, cost and duration estimation.

Project size (i) expresses the size of the software that is derived from the quantification of functional requirements specified by users [ISO/IEC14143]. Project size can be calculated by several methods and techniques of functional measurement such as FPA (Function Point Analysis) and COSMIC FP (COSMIC Function Point), thus it can be expressed in different units such as function points (FP) and source lines of code (SLOC).

The development effort (ii) is a function of the project size; it is expressed in man-hours, man-days or man-months. Duration estimation is either a function of project size or can be derived from the development effort. Effort and duration estimation, once estimated enable to calculate the project cost and staffing.

Various effort estimation methods can be used in the estimation process. They can be grouped in three main categories: (1) experience based methods, which is based on the expert intuition and experience drawn from previous executed project, such as expert judgment and analogy, (2) algorithmic model based methods, which are mainly based on equations expressing the effort as a function of discriminant parameters influencing the effort called effort drivers. Parametric models are established using historical data from complete projects, some of commonly used models are regression based models and Bayesian analysis based models [START_REF] Laqrichi | Toward an effort estimation model for information system project integrating risk[END_REF]. (3) Non-algorithmic model based methods, which model the relationship between the estimated variable and cost drivers using artificial intelligence techniques like NN and fuzzy logic. The relationship is not assumed to be well known or modelizable to specific shapes or equations [START_REF] Idri | Un modèle intelligent d'estimation des coûts de développement de logiciels[END_REF].

Regression is a widely used modeling technique and NN is a recent and evolutionary modeling technique. In this study, our attention was drawn to these two modeling techniques for the estimation activity of the global process (ii) because they seem to provide good estimates.

Regression models versus Neural Network models

Regression models are still the most popular models in literature; they include COCOMO [START_REF] Boehm | Software engineering economics[END_REF] and Putnam [START_REF] Basha | Analysis of Empirical Software Effort Estimation Models[END_REF]. Regression aims to model the relationship between inputs and outputs. In software estimation, the inputs are the discriminant parameters influencing the estimated variable called effort or cost drivers, the mean cost driver is the software size that is usually expressed in the Source Lines of Code. The output is the estimated variable that can be effort or duration or cost.

There are various types of regression that have been used in effort estimation models namely linear or multi linear regression [START_REF] Kok | The MERMAID Approach to software cost estimation[END_REF], non-linear regression that was [START_REF] Boehm | Software engineering economics[END_REF], and ordinal regression [START_REF] Sentas | Software productivity and effort prediction with ordinal regression[END_REF].

NN is a massively parallel adaptive network of simple nonlinear computing elements called Neurons, which are intended to abstract and model some of the functionality of the human nervous system in an attempt to partially capture some of its computational strengths [START_REF] Haykin | Neural Networks: A Comprehensive Foundation[END_REF].

NN is used in effort estimation due to its ability to learn from previous data. It is also able to model complex relationships between the dependent (effort or duration or cost) and independent variables (cost drivers). In addition, it has the ability to learn from the training data set thus enabling it to produce acceptable results for unseen data [START_REF] Boehm | Software engineering economics[END_REF]. But NN has one short coming that its estimation reason or relation between inputs and outputs cannot be justified. Different NN based model are proposed to predict and estimate effort and duration as COCOMO based NN model [START_REF] Barcelos Tronto | Comparison of artificial neural network and regression models in software effort estimation[END_REF][10] and radial basis function NN model [START_REF] Idri | Design of radial basis function neural networks for software effort estimation[END_REF].

Several research works have compared NN models with regression models. Finnie [START_REF] Finnie | A comparison of software effort estimation techniques: Using function points with neural networks, case-based reasoning and regression models[END_REF] compared regression model with two AI based estimation models that are case based reasoning and NN for software development effort. Authors found that AI based models perform better for complex software projects and outliers in training dataset than regression model. However, they gave no justification or explanation for the obtained results.

Heiat [START_REF] Heiat | Comparison of artificial neural network and regression models for estimating software development effort[END_REF] experimented FFNN (Feed-Forward Neural Network) with function point and RBNN (Radial Basis Neural Network) with SLOC for a dataset of 67 projects. Author concluded that NN approach is competitive with regression for same cases and significantly more accurate for others. This study presents some limitations: the size of data sets is small, the data sets used varied only in size and complexity varied in terms of language platform.

Framework of estimation of effort, duration and cost

The study prensented in this paper relies on the framework shown below (Fig 1 .1.)

Data preparation (a)

Organization's projects database is built over years by projects teams in order to capitalize the experience and information related to completed projects. It contains information about previous achieved projects such as project duration, project cost, project type, and platform development.

Steps for database preparation can be summarized as follows:

• Cleaning database: datasets related to irrelevant project parameters, such as parameters concerning information capitalization are discarded. Also, duplicate projects, that are projects with the same parameters but different estimates, are reduced. • Performing statistical tests: The projects database is explored to determine cost drivers. For this purpose, the statistical test of Pearson correlation and one-way ANOVA can be used [START_REF] Huang | A comparative evaluation on the accuracies of software effort estimates from clustered data[END_REF], they enables to examine the significance between the projects parameters and the variables to estimate in order to select the parameters with significant influence on these variables. The Pearson's correlation test is used for parameters with the ratio scale [START_REF] Neto | Pearson's Correlation Coefficient for Discarding Redundant Information in Real Time Autonomous Navigation System[END_REF] whereas One-Way Analysis of Variance (ANOVA) is used for parameters with the nominal scale. After determining cost drivers, the other parameters are discarded from the database. Projects with missing values in cost driver fields are discarded. Then, the variables to estimate are adjusted to be normal by discarding projects identified as outliers. A step of data normalization is required for the NN model establishment. • Dividing database: the projects data should be divided into two segments, one used to establish and train the effort estimation model and the other used to test and validate it. The holdout method or the k fold crossvalidation [START_REF] Refaeilzadeh | Cross-Validation[END_REF] approach can be used for this purpose. In this stage, the projects database is prepared in order to establish effort estimation models.

Establishment of estimation models (b1) Regression model establishment

The regression model establishment consists on modeling the relationship between dependant variables Y and independent variables X i in the form of a linear equation as:

Y=a 1 X 1 +a 2 X 2 +…+a n X n (1.2)
For our study, dependent variables are variables to estimate and independent variables are cost drivers. In order to establish the multi linear regression between these variables based on database, many statistical tools can be used such as XLSTAT. inputs, the number of output, the number of hidden layers and nodes, activation function). The number of inputs is the number of the projects characteristics, the number of outputs is the number of variable to estimate, the numbers of hidden layer and hidden layer nodes are less than or equal to twice the number of inputs [START_REF] Boetticher | An assessment of metric contribution in the construction of a neural network-based effort estimator[END_REF]. Activation function is used to transform and squash the amplitude of the output signal of a neuron to some finite value. Some of the most commonly used activation functions are sigmoid, Tanh, and Gaussian [START_REF] Karlik | Performance analysis of various activation functions in generalized MLP architectures of neural networks[END_REF]. There are a multitude of NN architecture and structure; the most used one is called Multilayer Perceptron (MLP) that is a feed forward artificial NN, i.e. the network is structured in a hierarchical way. It consists of different layers where the information flows only from one layer to the next layer. In this NN, Each node in one layer connects with a certain weight w ij to every node in the following layer. Input nodes distribute the signals to the nodes of the first hidden layer without processing it while nodes of hidden layers are neurons (or processing elements) with a nonlinear activation function [START_REF] Trenn | Multilayer perceptrons: approximation order and necessary number of hidden units[END_REF][START_REF] Rosenblatt | Principles of neurodynamics: perceptrons and the theory of brain mechanisms[END_REF] (view Fig 1 .2.).

Evaluation and comparison of models (c)

Model evaluation aims to measure how much model fits the context of the study. This context is defined by the variables to estimate, the database and the cost drivers used in estimation process.

The evaluating of estimation models makes it possible to compare them in order to choose the most adequate one. For this purpose, different accuracy indicators can be used for this study such as the Mean Magnitude of Relative Error (MMRE), the Pred(0.25) [START_REF] Conte | Software Engineering Metrics and Models[END_REF] [22] [START_REF] Moores | Developing a software size model for rule-based systems: a case study[END_REF].

Estimation of effort, duration or cost (d)

This last step consists on the implementation of the established models in order to estimate effort, duration or cost for new projects. For estimation using regression model, parameters of the model must be determined for the project to estimate.

Then the variable to estimate is simply calculated. For estimation using NN model, parameters must be determined then normalized in the same manner as projects database was normalized during database preparation (a). After that, the variable (s) to estimate is (are) then calculated then converted back to their real scale.

Case Study

The experiment described in this paper is carried on the data provided by an industrial company operating in the automotive sector. The main mission of the company is the design, production and sale of vehicles and mechanical components. The company is also involved in financing vehicle sales and dealership inventories. The database used in this case study consists on 6078 projects that are either carried or under way. These projects concern several domains of software industry (Cars, finance, commerce…) and their informations are organized into fourteen categories that involve 364 attributes. Due to lack of information about development effort, our focus will be put on the estimation of duration and cost.

Implementation of the methodology

The analysis of the database and the statistical tests enables to obtain a database of 214 projects and 4 attributes that are: project type, project BU, project difficulty, and domain. All cost drivers are qualitative, hence they should be transformed into dummy or binary variables [START_REF] Garavaglia | A smart guide to dummy variables: four applications and a macro[END_REF].

For regression model establishment (b1), software of data and statistic analysis is used in this work; it makes it possible to automatically establish the regression model using the database. As cost drivers are qualitative in this case study, a special case of linear regression called Analysis of variance (ANOVA) [START_REF] Gelman | Analysis of variance: Why it is more important than ever. Discussions. Author's reply[END_REF] is used to both transform cost drivers into binary variables and establish the estimation model.

For NN model establishment (b2), after transforming cost drivers into binary variables, outputs, that are duration and cost, are normalized to values between zero and one. The resulting database for NN model establishment consists of 27 binary cost drivers and 2 outputs. We use these models to estimate duration and cost of a project for which real duration and cost has been measured at the end of the project by the project management service. The proposed approach gives the results in Table 1.2.

Discussion

In this study case, effort estimation model cannot be established because information about effort is not provided in the database. Table 1.1 shows that, compared to regression model, NN model provides more accurate estimation for both duration and cost. This can be explained by the capacity of NN to model complex relationship between cost drivers and variables to estimate. A second explanation has to do with the complexity of the case study. This complexity manifests through the lack of relevant parameters due to the lack of information in database and the use of only qualitative cost drivers.

Our approach was applied on an example of a project with known achievement characteristics to concretely observe the estimation results (table1.2.), for this example NN model shows better results than regression as well. The differences between the estimated variables and the real ones may be due to the uncertainty in the model parameters and components. For neural network there is a significantly small uncertainty because NN has the ability to deal with the lack of data and cost drivers, in fact, it adjusts the model's weights so that it covers this lack. But for regression models, a complete database and a complete list of cost drivers is necessary in order to achieve good results. Otherwise, there will be a bigger uncertainty in regression model's coefficients.

Conclusion

The more accurate estimation is, the better the software project complies with the contractual commitments in terms of budget and duration. The model used for estimation is a crucial factor that affects estimation accuracy.

We presented a framework for estimation starting from the analysis of the available database up to the selection of the estimation model and its use on new projects. This framework is sufficiently flexible to provide estimation of different variables such as effort, duration and cost, depending on the available database about previous completed projects.

The attention was drawn to two models that are regression and NN models. The proposed framework was then applied on an industrial study case that consists of multisite IT projects. This study has shown that NN model is more accurate than regression model even with an important lack of information about previous projects. This lack of information can explain the uncertainty in estimations. Thus, it will be important to be able to measure this uncertainty in order to take it into account in the estimation process.

Future research studies can focus on the need of more realistic estimations by providing not a single value but an interval of estimation and a degree of trust associated to this interval.

Fig 1 . 1 .

 11 Fig 1.1. Framework of estimation for a software project This methodology consists on three major steps: • (a) Preparation of data: based on raw database containing information collected from previous achieved software projects • (b) Establishment of estimation models: consists on the establishment of estimation models to compare. For this study, the estimation models processed and compared within this framework are (b1) regression model and (b2) NN model. • (c) Evaluation and comparison of models using evaluation criteria and performance indicators. • (d) Estimates of effort, duration or cost using the selected estimation model These four steps are detailed below.

(Fig 1 . 2 .

 12 Fig 1.2. Algorithm for NN model establishment

Fig 1 . 3 .

 13 Fig 1.3. Architecture of a Multilayer Perceptron (MLP) 2. Training: The two main techniques employed by neural networks are known as supervised learning and unsupervised learning. In unsupervised learning, the NN requires no initial information regarding the correct classification of the data it is presented with. Supervised training works in much the same way as a human learns new skills, by showing the network a series of examples. The most used supervised training algorithm is back propagation algorithm [13]. The objective of training is to determine the weights of the NN to bring the outputs of the whole network closer to the desired outputs. 3. Test: After training NN, the test is performed on the database reserved for testing. It makes it possible to measure the potential of success of the trained NN using evaluation criteria. As the figure show, if the results of test are not satisfying, the architecture of the NN model is modified until reaching good error.

Table 1 . 1 .

 11 Comparison of results

	Evaluation criteria	Duration estimation models	Cost estimation models
		Using	Using neural	Using	Using neural
		regression	network	Regression	network
	MMRE %	69	23	658	14
	Pred(0.25) %	20	72	7	76

Table 1 . 2 .

 12 Estimation for a project

	Variables	to	Real achievement Estimation using	Estimation using
	estimate			neural network	regression
	Duration (day)		73	62,8	36
	Cost (K€)		14,02	17,18	4556,5

Acknowledgment

This work has been funded by the Fund Unique Interministerial (FUI) through the project Projestimate. We wish to acknowledge our gratitude and appreciation to all project partners for their contribution.

For this case study, MLP architecture is used with 27 inputs and 2 outputs. In order to determine the hidden layers numbers, training is repeated many times with the variation of the hidden layers number. The best number of hidden layers is that which provide the best performance in test phase. For training, as many experiences have shown that most NN are trained enough in less than 1000 epochs [START_REF] Boetticher | A neural network paradigm for characterizing reusable software[END_REF], the number of epochs in this study is set to 1000 epochs.

Experimental results

The established regression and NN model for estimating duration and cost are applied on the case study, and then the evaluation criteria are calculated. Table 1.1. presents accuracy indicators calculated for duration and cost estimation using regression and NN.