Julien Coatl Éven

Houssem Haddar

Jing-Rebecca Li

A MACROSCOPIC MODEL INCLUDING MEMBRANE EXCHANGE FOR DIFFUSION MRI

Diffusion Magnetic Resonance Imaging is a promising tool to obtain useful information on the microscopic structure and has been extensively applied to biological tissues. We establish a new macroscopic model from homogenization theory for the complex transverse water proton magnetization in a voxel due to diffusion-encoding magnetic field gradient pulses in the case of intermediate water exchange across biological cellular membranes. Based on a particular scaling of the permeability condition modeling cellular membranes, this macroscopic model reproduces the memory effects often observed in experiments. Explicit formulae given by homogenization for the coefficients of this model emphasize their link to the relevant physiological quantities. In addition, we explicitly solve the macroscopic model to obtain an ODE model for the dMRI signal. This ODE model is numerically easy to invert, and the inverse problem of retrieving model coefficients from synthetic dMRI signal data is considered.

1. Introduction. The image contrast in diffusion Magnetic Resonance Imaging (dMRI) comes from the differing water diffusion characteristics in the imaged tissue at different spatial positions. DMRI can be used to detect, for example, cerebral ischemia [START_REF] Warach | Fast magnetic resonance diffusion-weighted imaging of acute human stroke[END_REF], demyelinating disorders [START_REF] Horsfield | Applications of diffusion-weighted and diffusion tensor mri to white matter diseases, a review[END_REF], and tumors [START_REF] Schnapauff | Diffusion-weighted echo-planar magnetic resonance imaging for the assessment of tumor cellularity in patients with soft-tissue sarcomas[END_REF][START_REF] Sugahara | Usefulness of diffusion-weighted mri with echo-planar technique in the evaluation of cellularity in gliomas[END_REF][START_REF] Tsushima | Magnetic resonance (mr) differential diagnosis of breast tumors using apparent diffusion coefficient (adc) on 1.5-t[END_REF]. For a survey, see [START_REF] Le | Diffusion mri at 25: Exploring brain tissue structure and function[END_REF].

The signal measured by the MRI scanner is a mean-value measurement in a physical volume, called a voxel, whose size is much larger than the scale of the microscopic variations of the cellular structure. The resolution of dMRI is on the order of 1 mm 3 , meaning the dMRI signal combines the diffusion characteristics of a tissue volume (voxel) of 1 mm 3 . This is very large compared to cell features, which vary from sub-µm (diameter of neurites) to tens of µm (diameter of neuronal bodies and glial cells) in the brain (see for example, [START_REF] Hellwig | A quantitative analysis of the local connectivity between pyramidal neurons in layers 2/3 of the rat visual cortex[END_REF]). In other words, dMRI is used to show the averaged characteristics of the microscopic structure on a macroscopic scale. Another very important spatial scale to consider for this work is the diffusion distance, which gives an indication of the displacement of water molecules during the measured diffusion time. At physically realistic dMRI diffusion times of 10-100 ms, the average diffusion distance is, assuming an average diffusivity of 10 -3 mm 2 /s, between 8-25 µm. Thus, from the point of view of diffusing water molecules, the cellular structure beyond tens of µm from their starting positions do not contribute to their displacements.

Most of the macroscopic models for the dMRI signal were obtained in the limit of long diffusion time. For example, if the diffusion time is long enough that all the water molecules experience the same Gaussian diffusion environment, the macroscopic model is simply described by a single effective diffusion tensor. In the presence of permeable cell memebranes, various analytical formulae for the effective diffusion tensor for special geometries such as cubes and spheres can be found in [START_REF] Hasselman | Effective thermal conductivity of composites with interfacial thermal barrier resistance[END_REF][START_REF] Torquato | Effect of the interface on the properties of composite media[END_REF][START_REF] Latour | Time-dependent diffusion of water in a biological model system[END_REF]. For general cell shapes, while assuming no contiguous extra-cellular space, an analytical formula for the effective diffusion tensor can be found in [START_REF] Dmitry | Random walks with barriers[END_REF], where the time dependence of the average diffusion distance is considered. In [START_REF] Cheng | Effective conductivity of periodic arrays of spheres with interfacial resistance[END_REF] the effective diffusion tensor for general cell shapes and allowing a contiguous extra-cellular space was formulated based on the solution of a set of Laplace PDEs.

There is ample experimental evidence that the diffusion in brain tissue at diffusion times relevant to dMRI is not Gaussian: the dMRI signal has been fitted as multiple compartmental Gaussian [START_REF] Niendorf | Biexponential diffusion attenuation in various states of brain tissue: Implications for diffusionweighted imaging[END_REF][START_REF] Mulkern | Multi-component apparent diffusion coefficients in human brain[END_REF][START_REF] Clark | Water diffusion compartmentation and anisotropy at high b values in the human brain[END_REF][START_REF] Maier | Normal brain and brain tumor: Multicomponent apparent diffusion coefficient line scan imaging1[END_REF], or by including a Kurtosis term to quantify the non-Gaussianness [START_REF] Chabert | Diffusion tensor imaging of the human optic nerve using a non-cpmg fast spin echo sequence[END_REF][START_REF] Jensen | Diffusional kurtosis imaging: The quantification of non-gaussian water diffusion by means of magnetic resonance imaging[END_REF], or by fitting with fractional order diffusion [START_REF] Magin | Anomalous diffusion expressed through fractional order differential operators in the bloch-torrey equation[END_REF].

There have been many works on formulating macroscopic models based on multiple compartment diffusion. For example, the no-exchange models of [START_REF] Pabitra | A model for diffusion in white matter in the brain[END_REF][START_REF] Assaf | Axcaliber: A method for measuring axon diameter distribution from diffusion mri[END_REF][START_REF] Sune | Modeling dendrite density from magnetic resonance diffusion measurements[END_REF] separate the cylindrical-shaped neurites from the space outside them to make different diffusion compartments. When there is no water exchange between the compartments, the total dMRI signal is just the sum of the signals from each of the compartments. We note that in these no-exchange models, one can in fact easily allow non-Gaussian diffusion in each compartment, and this flexibility can potentially result in the applicability of the models at shorter diffusion times. However, the weakness of the no-exchange models is that water is assumed to stay (mostly) in a single compartment during the diffusion time. Even though no-exchange models may be appropriate for the brain white matter, where the myelin around axons prevents significant water exchange between the axons and the extra-cellular space, it is not clear that in the gray matter, where the cell membranes are permeable, one can neglect the water exchange between cells and the extra-cellular space.

In this paper, we develop a macroscopic model based on multiple diffusion compartments, where the diffusion within each compartment is considered Gaussian. To our knowledge, up to the present, there has been only one multiple compartment macroscopic model for dMRI in heterogeneous media that allows exchange between the compartments. That model is the Karger model [START_REF] Karger | Principles and application of self-diffusion measurements by nuclear magnetic resonance[END_REF], which was originally developed for micro-porous crystallites. The Karger model was obtained on the basis of phenomenological modeling of the experimentally obtained signal curves and only works under a strong and often unrealistic assumption on the applied diffusion-encoding magnetic field gradient sequence. This assumption is that the duration of the gradient pulse is very short compared to the diffusion time and is called the narrow pulse assumption. Nevertheless, the Karger model has been used in biological tissue imaging [START_REF] Waldeck | Nmr diffusion measurements to characterise membrane transport and solute binding[END_REF][START_REF] Stanisz | An analytical model of restricted diffusion in bovine optic nerve[END_REF][START_REF] Pfeuffer | Restricted diffusion and exchange of intracellular water: theoretical modelling and diffusion time dependence of 1h nmr measurements on perfused glial cells[END_REF][START_REF] Lee | Effects of equilibrium exchange on diffusion-weighted nmr signals: The diffusigraphic shutter-speed?[END_REF][START_REF] Quirk | Equilibrium water exchange between the intra-and extracellular spaces of mammalian brain[END_REF][START_REF] Meier | Diffusion in compartmental systems. i. a comparison of an analytical model with simulations[END_REF][START_REF] Roth | Quantification of water compartmentation in cell suspensions by diffusion-weighted and t2-weighted mri[END_REF][START_REF] Aslund | Filter-exchange pgse nmr determination of cell membrane permeability[END_REF] to invert for model parameters.

In this paper, starting from a microscopic model of the complex transverse water proton magnetization due to diffusion-encoding magnetic field gradient pulses in a spatial volume on the order of a voxel, which mathematically takes the form of a Bloch-Torrey partial differential equation (PDE) in a heterogeneous medium with barriers, varying at the scale of the biological cell features, we will provide a macroscopic limit model, through homogenization, that is applicable to general diffusion-encoding magnetic field gradients, i.e., not subject to the narrow pulse assumption.

The main difficulty in the modeling of diffusion at the microscopic level is the choice of the membrane law. We adopted here the so-called two compartment model where the membrane is infinitesimally thin and induces a discontinuity of the field proportional to the flux. The nature of the homogenized model then would depend on the scaling for the proportionality coefficient with respect to the periodicity size parameter. A comprehensive presentation as well as a fully rigorous mathematical justification of all the possible limit models is the subject of a forthcoming work [START_REF] Coatléven | Mathematical justification of macroscopic models for diffusion mri through the periodic unfolding method[END_REF].

We will concentrate in this article on the scaling that provides the limit model that is closest to the Karger model. We shall restrict ourselves to only a formal derivation of the homogenized model and put more emphasis on obtaining the simplest possible expression of the dMRI signal, which is the spatial integral of the complex transverse magnetization over the voxel at one time point called the echo time. The signal associated to the new macroscopic (and homogeneous) model that we derive in this paper, being the limit, when the size of biological cells tends to zero, of the signal of the microscopic model, is thus expected to be an approximation of the non-homogenized one at long diffusion times. The equations governing the magnetization are of course much simpler than those involving the microscopic description of the medium. Their coefficients can be explicitly linked, in the context of periodic homogenization, to relevant properties of the underlying biological tissue.

In addition, if we enforce periodic boundary conditions on the voxel to be modeled, we can explicitly solve the macroscopic (PDE) model to obtain an ODE model for the dMRI signal. The Karger model can be seen as an approximate version of our ODE model in the narrow pulse limit. The main interest of this ODE model is indeed to solve the inverse problem of identifying the model coefficients that make our approximate signal fit best the measured signal. Some preliminary inversion tests are performed on synthetic examples.

We plan a thorough numerical study of the macroscopic model we obtained here in the context of biological tissue dMRI, using more complex geometries than those shown in this paper, and will conduct a study of the convergence properties of our model as well as a comparison with the Karger model. Thus, these aspects are not included in this paper.

The paper is organized as follows. In Section 2, we review the Bloch-Torrey PDE with barriers that describes the complex transverse water proton magnetization due to diffusion-encoding magnetic field gradient pulses at the microscopic scale. For simplicity, we then make the hypothesis that the domain to be modeled is periodic, to allow us to apply periodic homogenization theory. In Section 3, we make the formal homogenization of our model problem in the periodic context, using a particular scaling for the permeability condition on the cellular membranes. In Section 4, we provide an ODE description of the dMRI signal. We then show with a few numerical experiments that the signal from the homogenized model accurately approximates the signal of the microscopic model at long diffusion times. In Section 5 we solve numerically the inverse problem of finding the coefficients of our ODE model from the signal obtained without homogenization, and show that through these coefficients we can recover some relevant biological properties of the tissue. In Section 6 we discuss some issues regarding the use of periodic homogenization and the periodic extension of the voxel to R d , where d is the spatial dimension, in the context of the spatial and time scales of dMRI. Section 7 contains our conclusions.

2. Model Problem. Biological tissues contain abundant water and are composed of cells of various sizes and shapes surrounded by the extra-cellular space. The cells can be represented, at the scale of interest for our problem, as a bounded domain surrounded by a very thin membrane (see Figure 2.1 left). The diffusion coefficients for water in the interior part of the cells, in the membrane, and in the extra-cellular space may be different from each other.

Γ I Ω e Ω c
While the diffusion coefficients in the interior part of the cells and in the extra-cellular space remain of the same order of magnitude, it is much smaller inside the membranes. It is in fact very difficult to obtain direct measurements of this coefficient inside the cellular membranes, as the membranes are extremely thin. This is why the membranes are most of the time spatially neglected and replaced by a permeability condition between the extra-cellular and intra-cellular parts [START_REF] Price | Pulsed-field gradient nuclear magnetic resonance as a tool for studying translational diffusion: Part 1. basic theory[END_REF]. Consequently, we replace the realistic, "three-compartment" biological cell by the idealized "twocompartment" cell displayed in Figure 2.1 (right).

For a volume Ω of biological tissue, we denote by Γ I the union of the boundaries of all the "two-compartment" cells included in Ω. Γ I thus delimits two subdomains: the extracellular domain Ω e (e standing for extra-cellular) and the intra-cellular domain Ω c (c standing for cellular) . The domain Ω ext then represents the union of the open extra-cellular and open intra-cellular domains, i.e.

Ω ext = Ω \ Γ I = Ω e ∪ Ω c .
A classic dMRI experiment consists of applying two pulsed (meaning short duration in time) gradient (meaning linearly varying in space) magnetic fields with a 180 degree spin reversal in between the two in order to mark the positions of the water molecules between the two pulses [START_REF] Stejskal | Spin diffusion measurements: Spin echoes in the presence of a timedependent field gradient[END_REF].

The complex transverse water proton magnetization M is modeled by the following Bloch-Torrey PDE [START_REF] Hc | Bloch equations with diffusion terms[END_REF] with jump:

∂M (x, t) ∂t + ıq • xf (t)M (x, t) -div (σ(x)∇M (x, t)) = 0 in Ω ext ×]0, T [, σ∇M • ν |Γ I = κ [M ] |Γ I on Γ I , [σ∇M • ν] |Γ I = 0 on Γ I , M (•, 0) = M init in Ω ext ,
where ν is the normal exterior to the intra-cellular domains, [.] |Γ I is the jump (extra-cellular minus intra-cellular) on Γ I for a quantity defined on both parts of the domain, κ is the permeability coefficient, ı = √ -1, M init is the initial magnetization. The constant vector in R d , q, contains the amplitude and direction information of the applied diffusion-encoding magnetic field gradient multiplied by the gyro-magnetic ratio of the water proton, and f , where max t f (t) = 1, is the normalized time profile of the diffusion-encoding magnetic field gradient sequence. The coefficient σ is the intrinsic diffusion coefficient. The time profile of the classic Pulsed Gradient Spin Echo (PGSE) [START_REF] Stejskal | Spin diffusion measurements: Spin echoes in the presence of a timedependent field gradient[END_REF] sequence (simplified to include only the parameters relevant to diffusion, i.e., the imaging gradients are ignored) is the following:

     f (t) = 1, 0 < t ≤ δ, f (t) = -1, ∆ < t ≤ ∆ + δ, f (t) = 0, elsewhere, (2.1)
where we made f (t) negative in the second pulse to include the effect of the 180 degree spin reversal between the pulses. The time at which the signal is measured is called the echo time T E > δ +∆.

For simplicity and in order to be able to apply a well known theoretical framework, we will assume that the volume to be modeled, Ω, can be described as a periodic domain. More precisely, we will assume that there exists a period ε, which represents the average size of a "representative" volume of Ω, and which is small compared to the size of Ω. We will discuss in Section 6 this choice of periodic homogenization and the applicability of our obtained model to non-periodic geometries.

We now define the unit periodicity cell

Y =]0, 1[ d such that Y = Y e ∪ Y c
, where Y c is the intracellular domain and is an open set which can be made of several connected parts (i.e it is the union of the interiors of the selection of biological cells that are included in the periodicity cell).

Ξ ε = ξ ∈ Z d | ε(ξ + Y ) ∩ Ω = ∅ ,
so that Ω will contain almost ♯(Ξ ε ) periodicity cells. Finally, we denote

Ω ε e = ξ∈Ξε ε(ξ + Y e ) ∩ Ω, Ω ε c = ξ∈Ξε ε(ξ + Y c ) ∩ Ω, Ω ε ext = Ω ε e ∪ Ω ε c . (2.2)
Notice that

∂Ω ε c ∩ ∂Ω ε e = ∂Ω ε e \ ∂Ω = ξ∈Ξε Γ ε,ξ m where Γ ε,ξ m = ε(ξ + Γ m ) ∩ Ω.
Of course, the diffusion coefficient will be assumed to be periodic as well, i.e., there exists σ

∈ L ∞ (Y ) such that σ(x) = σ x ε , with σ = σ c in Y c , σ e in Y e . (2.3)
The most common practical choice for σ e and σ c is to take them constant, so that σ is piecewise constant. For both physical and technical reasons, we will further assume that there exists two constants 0 < σ -< σ + < +∞ such that, for almost every y ∈ Y :

σ -≤ σ c ≤ σ + and σ -≤ σ e ≤ σ + .
With these more precise description of the domain, our model problem can be rewritten as:

∂M ε (x, t) ∂t + ıq • xf (t)M ε (x, t) -div (σ ε (x)∇M ε (x, t)) = 0 in Ω ε ext ×]0, T [, σε ∇M ε • ν | Γ ε,ξ m = κ ε [M ε ] | Γ ε,ξ m ∀ ξ ∈ Ξ ε , [σ ε ∇M ε • ν] | Γ ε,ξ m = 0 ∀ ξ ∈ Ξ ε , M ε (•, 0) = M init in Ω ε ext , (2.4)
where σε = σ x ε . Finally, we will assume that the time profile f belongs to L ∞ (]0, T [) and that the initial data M init is defined on Ω independently of ε and belongs to L 2 (Ω) ∩ L 1 (Ω). For notation convenience we set Q(x) := q • x.

It only remains to choose the volume to be modeled, Ω. Obviously, it should be on the order of a voxel. However, in the derivation of the macroscopic model that follows, to avoid complications associated with imposing boundary conditions on Ω, we make the simplification that Ω = R d , where d is the space dimension. This is a reasonable simplification in context of the mathematical derivation that follows because we develop the macroscopic model in the limit of a period ε that is small compared to the size of Ω. Thus, from the point of view of the "representative" volume, the error associated with replacing the true Ω with R d is small. This choice of Ω also allows the formulation of a system of ODEs for the dMRI signal, which is the spatial integral of the magnetization. If we had chosen a finite Ω, with associated boundary conditions, the form of the macroscopic PDE model for the magnetization would not change, only the associated boundary conditions would have to be added to the macroscopic PDE. However, in that case, there would not be a simple ODE model for the dMRI signal. A discussion of the boundary condition on Ω in the context of dMRI is included in Section 6.

Finally, we conclude this section by explaining our choice for the permeability coefficient κ ε . As usual in homogenization, the choice of a scaling for this coefficient, driven by the fact that it is experimentally a very small parameter, will be determinant with regard to the limit model. More precisely, the dependence of κ ε on ε → 0 will influence the nature of the limit macroscopic model. A comprehensive presentation as well as a fully rigorous mathematical justification of all the possible limit models corresponding to a scaling κ ε = σ m ε m , for all m ∈ Z, is the subject of a forthcoming work [START_REF] Coatléven | Mathematical justification of macroscopic models for diffusion mri through the periodic unfolding method[END_REF]. We will concentrate in this article on the relevant case for dMRI, i.e., the scaling that provides the limit model that shows a "memory" effect, which corresponds with m = 1. We shall restrict ourselves to only a formal derivation of the homogenized model and put more emphasis on obtaining results on the integral of the magnetization, which approximates the dMRI signal.

3. A macroscopic model through two scale asymptotic expansions. In this section, we derive our macroscopic model through formal homogenization from the two compartment model. We note here that M ε does not satisfy the Bloch-Torrey equation in all Ω ε ext , but separately in Ω ε c

and Ω ε e , with jump conditions on the Γ ε,ξ m 's. This remark will be important to our formal derivation of the macroscopic models.

3.1. Formal derivation of the macroscopic models through two-scale asymptotic expansions. Using classical periodic homogenization techniques [START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF], we are going to develop M ε using two-scale asymptotic expansions. According to our previous remark, we are not going to introduce a unique expansion in all Ω ε ext , but two expansions, one for each phase Ω ε c and Ω ε e . Consequently, we write

M ε (x, t) = M e ε (x, t) = +∞ i=0 ε i M i,e x, x ε , t in Ω ε e , M ε (x, t) = M c ε (x, t) = +∞ i=0 ε i M i,c x, x ε , t in Ω ε c , (3.1)
where the functions M i,e (x, y, t) and M i,c (x, y, t) are defined respectively on Ω × Y e ×]0, T [ and Ω × Y c ×]0, T [, and the M i,e (x, y, t) are assumed Y -periodic in y. The aim of such an ansatz is to separate the "macroscopic" variations (the x variable) form the "microscopic" ones (the y variable), and then try to obtain a new problem involving only the macroscopic one.

3.1.1. Periodic cell equations and jump conditions. To get the equations for each of the M i,e 's and the M i,c 's, we start by noticing that, α ∈ {c, e},

∇M i,α (x, x ε , t) = ∇ x M i,α (x, x ε , t) + ε -1 ∇ y M i,α (x, x ε , t),
and therefore

div σ α x ε ∇M i,α (x, x ε , t) = div x σ α x ε ∇ x M i,α (x, x ε , t) + ε -2 div y σ α x ε ∇ y M i,α (x, x ε , t) +ε -1 div x σ α x ε ∇ y M i,α (x, x ε , t) + ε -1 div y σ α x ε ∇ x M i,α (x, x ε , t) .
Substituting these formulae into the two compartment model, using the asymptotic expansions (3.1) then matching the terms in front of the same power of ε, one gets for the three first orders,

-div y (σ α ∇ y M 0,α ) = 0 in Ω × Y α ×]0, T [, (3.2)
-div y (σ α ∇ y M 1,α ) = div y (σ α ∇ x M 0,α ) + div x (σ α ∇ y M 0,α ) in Ω × Y α ×]0, T [, (3.3)
-div y (σ α ∇ y M 2,α ) = div y (σ α ∇ x M 1,α ) + div x (σ α (∇ y M 1,α + ∇ x M 0,α )) - ∂M 0,α ∂t -ıQf (t)M 0,α in Ω × Y α ×]0, T [. (3.4)
Notice that the equations have been formally expanded to all Ω for the macroscopic variable x.

Next, we make the following ansatz for the jumps of M ε and its fluxes, for

x ∈ Γ ε m , [M ε ] | Γ ε m (x, t) = ∞ i=0 ε i M i,e (x, x ε , t) -M i,c (x, x ε , t) , (3.5)
[σ ε ∇M ε • ν] | Γ ε m = ∞ i=0 ε i σ e ∇M i,e (x, x ε , t) • ν -σ c ∇M i,c (x, x ε , t) • ν . (3.6)
Then, the no-jump relation for the fluxes becomes

ε -1 (σ e ∇ y M 0,e • ν -σ c ∇ y M 0,c • ν) + ∞ i=0 ε i (σ e ∇ y M i+1,e • ν + σ e ∇ x M i,e • ν -σ c ∇ y M i+1,c • ν -σ c ∇ x M i,c • ν) = 0,
which gives, after identifying each power of ε and expanding to all (x, y)

∈ Ω × Γ m σ e ∇ y M 0,e • ν = σ c ∇ y M 0,c • ν in Ω × Γ m ×]0, T [, (3.7)
and for i ≥ 1,

σ e ∇ y M i,e • ν + σ e ∇ x M i-1,e • ν = σ c ∇ y M i,c • ν + σ c ∇ x M i-1,c • ν in Ω × Γ m ×]0, T [. (3.8)
To write the equivalent condition for the traces, recall that we have assumed that κ ε = εσ m , where σ m > 0 is a constant independent on ε. Now we write the jump relations for traces, using (3.5) and (3.6) as

+∞ i=0 ε i+1 σ m (M i,e -M i,c ) = ε -1 σ α ∇ y M 0,α • ν + ε 0 σ α (∇ y M 1,α • ν + ∇ x M 0,α • ν) +ε 1 σ α (∇ y M 2,α • ν + ∇ x M 1,α • ν) + +∞ i=2 ε i σ α (∇ y M i+1,α • ν + ∇ x M i,α • ν).
3.1.2. Formal limit problem. As usual when performing formal homogenization through two-scale asymptotic expansions, we only need (as we shall see) to obtain the equations for the first three terms, M 0,α , M 1,α and M 2,α . The problem for M 0,α is then given by

-div y (σ α ∇ y M 0,α ) = 0 in Y α , σ α ∇ y M 0,α • ν = 0 on Γ m , M 0,e is Y -periodic.
(3.9)

Let us recall the following classical lemma (see for instance [START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF][START_REF] Brezis | Analyse fonctionnelle, théorie et applications[END_REF]) on solutions to this type of problems that will also be useful to derive the macroscopic model.

Lemma 1. Let f c ∈ L 2 (Y c ), f e ∈ L 2 (Y e ), ψ c and ψ e ∈ H -1/2 (Γ m ). Then, there exists a unique solution u c ∈ H 1 (Y c ), up to a constant, to -div y (σ c ∇ y u c ) = f c in Y c , σ c ∇ y u c • ν = ψ c on Γ m , if and only if, for any i ∈ 1, N c , Yc,i f c + < ψ c , 1 > H -1/2 (Γm∩∂Yc,i),H 1/2 (Γm∂Yc,i) = 0,
f e -< ψ e , 1 > H -1/2 (Γm∩∂Ye,i),H 1/2 (Γm∩∂Ye,i) = 0, where H 1 ♯ (Y e )
is the space of H 1 and Y -periodic functions on Y e . The above lemma tells us that neither M 0,e nor M 0,c depend on y. The problem for M 1,α is consequently

-div y (σ α (∇ y M 1,α + ∇ x M 0,α )) = 0 in Y α , σ α ∇ y M 1,α • ν + σ α ∇ x M 0,α • ν = 0 on Γ m , M 1
,e is Yperiodic in the y variable.

(3.10) Thus, introducing for i = 1, d the cell problems

-div y (σ α (∇ y w i,α + e i )) = 0 in Y α σ α ∇ y w i,α • ν + σ α e i • ν = 0 on Γ m w i,e is Y -periodic, (3.11)
which are well posed (in appropriate H 1 spaces) according to Lemma 1 since, using the divergence theorem, one can express the solutions to (3.10) as

M 1,α = d i=1 w i,α ∂M 0,α ∂x i in Y α . (3.13)
Now, we consider the equations satisfied by M 2,α where the variables x and t are treated as parameters.

-div

y (σ α (∇ y M 2,α + ∇ x M 1,α )) = div x (σ α (∇ y M 1,α + ∇ x M 0,α )) - ∂M 0,α ∂t -ıQf (t)M 0,α in Y α , σ e ∇ y M 2,e • ν + σ e ∇ x M 1,e • ν = σ m (M 0,e -M 0,c ) on Γ m , σ c ∇ y M 2,c • ν + σ c ∇ x M 1,c • ν = σ m (M 0,e -M 0,c ) on Γ m , M 2
,e is Yperiodic in the y variable.

(3.14) Using lemma 1, it is easy to see that the compatibility condition for M 2,e implies

Ye div x (σ e (∇ y M 1,e + ∇ x M 0,e )) - ∂M 0,e ∂t -ıQf (t)M 0,e - Γm σ m (M 0,e -M 0,c ) = 0, (3.15)
while the compatibility condition for M 2,c implies

Yc div x (σ c (∇ y M 1,c + ∇ x M 0,c )) - ∂M 0,c ∂t -ıQf (t)M 0,c + Γm σ m (M 0,e -M 0,c ) = 0. (3.16)
We define the homogenized tensors D α through

D α,ij := 1 |Y α | Yα σ α ∇w j,α • e i + σ α e j • e i ,
i, j = 1, d which can be also expressed in a symmetric form

D α,ij = 1 |Y α | Yα σ α (∇w j,α + e j ) • (∇w i,α + e i ), (3.17)
since using (3.11),

Yα σ α (∇w j,α + e j ) • ∇w i,α = - Yα div y (σ α (∇w j,α + e j ))w i,α = 0.
We also define the two coefficients

η c := σ m |Γ m | |Y c | and η e := σ m |Γ m | |Y e | . ( 3.18)
Then, using the fact that M 0,e and M 0,c does not depend on y, we obtain from (3.15) and (3.16) the macroscopic model

∂M 0,e ∂t + ıQf (t)M 0,e -div x (D e ∇ x M 0,e ) + η e (M 0,e -M 0,c ) = 0 in Ω×]0, T [, M 0,e (•, 0) = M init in Ω, ∂M 0,c ∂t + ıQf (t)M 0,c -div x (D c ∇ x M 0,c ) + η c (M 0,c -M 0,e ) = 0, in Ω×]0, T [, M 0,c (•, 0) = M init in Ω, (3.19)
which is a coupled system of modified Bloch-Torrey equations with homogeneous diffusion tensors.

3.2.

A simplified macroscopic model. For our setting of the problem where we assumed that Y c does not touch the boundary of Y the homogenized system (3.19) can be simplified since in that case D c = 0. This can be easily seen by checking that w i,c = -y i and therefore ∇ y w i,c + e i = 0. Consequently our model simplifies to

∂M 0,e ∂t + ıQf (t)M 0,e -div x (D e ∇ x M 0,e ) + η e (M 0,e -M 0,c ) = 0 in Ω×]0, T [, M 0,e (•, 0) = M init in Ω, ∂M 0,c ∂t + ıQf (t)M 0,c + η c (M 0,c -M 0,e ) = 0 in Ω×]0, T [, M 0,c (•, 0) = M init in Ω. (3.20)
Remark that the equation for M 0,c can be solved in terms of M 0,e

M 0,c (x, t) = M init G c (x, t, 0) + t 0 G c (x, t, s)M 0,e (x, s)ds, (3.21)
where we have set

G c (x, t, s) := exp - t s (ıq • xf (r) + η c )dr . (3.22)
Thus, we can decouple the system (3.20) into

∂M 0,e ∂t + (ıQf (t) + η e )M 0,e -div x (D e ∇ x M 0,e ) -η e t 0 G c (t, s)M 0,e (s)ds = η e M init G c (t, 0) in Ω×]0, T [, M 0,e (•, 0) = M init in Ω, M 0,c = M init G c (t, 0) + t 0 G c (t, s)M 0,e (s)ds in Ω×]0, T [. (3.23)
The first equation of (3.23) emphasizes the fact that this macroscopic model will behave quite differently from a classical Bloch-Torrey equation. In particular, the presence of the integral with respect to time gives birth to memory effects for M 0,e . Let us also remark that since the boundary of Y is contained in the boundary of Y e , then the homogenized tensor D e is positive definite as soon as σ c is also positive definite [START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF]. Then we can expect some well posedness of the homogenized problem. We recall that for simplicity we take Ω = R d . With this choice, problem (3.20) involves an unbounded coefficient, namely Q(x), and therefore its analysis is not classical. In particular, it is not clear in which sense problem (3.20) is to be understood. To make it more precise, we set q = qn, where n is a unitary vector of R d , and define M α almost everywhere on R d ×]0, T [, for α ∈ {e, c} by

M 0,α (x, t) = M 0,α (x, t)e iqx•n t 0 f (s)ds , (3.24)
denoting q = qn where n is a unitary vector. Then, at least formally, and with

F (t) := t 0 f (s)ds, ∂ M 0,α ∂t = ∂M 0,α ∂t + ıqn • xf (t)M 0,α e ıqn•xF (t) , ∇ M 0,α = (∇M 0,α + ıqnM 0,α F (t)) e ıqn•xF (t) , div(D α ∇ M 0,α ) = (div(D α ∇M 0,α ) + ıqD α ∇M 0,α • nF (t)) e ıqn•xF (t) + div ıqD α M 0,α • nF (t) .
Then, since

ıqD α ∇M 0,α • nF (t)e ıqn•xF (t) = ıqD α ∇ M 0,α • nF (t) + q 2 D α M 0,α F (t) 2 ,
we get

div(D α ∇ M 0,α ) = div(D α ∇M 0,α )e ıqn•xF (t) +ıqD α ∇ M 0,α •nF (t)+q 2 D α F (t) 2 M 0,α +div ıqD α M 0,α • nF (t) ,
and finally obtain

∂ M 0,e ∂t -div D e ∇ M 0,e -ıqnF (t) M 0,e + ıqD e ∇ M 0,e • nF (t) + q 2 D e F (t) 2 M 0,e +η e ( M 0,e -M 0,c ) = 0 in R d ×]0, T [, ∂ M 0,c ∂t + η c ( M 0,c -M 0,e ) = 0 in R d ×]0, T [, M 0,e (•, 0) = M init and M 0,c (•, 0) = M init in R d . (3.25)
While it may seem a more complicated problem, notice that it does not involve any unbounded coefficients, and its analysis enters a well-known theoretical framework.

Theorem 1. Let f ∈ L ∞ (0, T ) and M init ∈ L 2 (Ω). We denote H = L 2 (Ω) × L 2 (Ω), X = H 1 (Ω) × L 2 (Ω), X ′ = H -1 (Ω) × L 2 (Ω) and W (0, T, X) = U ∈ L 2 (0, T, X) | ∂ t U ∈ L 2 (0, T, X ′ )
There exists a unique solution M = ( M e , M c ) ∈ W (0, T, X) ∩ C 0 (0, T, L 2 (R d )) 2 to problem (3.25) which moreover satisfies for some C, C ′ > 0,

|| M || L ∞ (0,T,H) ≤ √ 2e CT ||M init || H , (3.26)
and

|| M || L 2 (0,T,X) ≤ C ′ (1 + e CT )||M init || H . (3.27)
Proof. We denote < •, • > H -1 ,H 1 the duality product between H -1 (R d ) and H 1 (R d ) and (•, •) L 2 the scalar product on L 2 (R d ). If there exists a solution M ∈ W (0, T, X), we can use any function V ∈ X as a test function for problem (3.25), and we get after formally integrating by parts

∂ t M 0,e , V e H -1 ,H 1 + D e ∇ M 0,e , ∇V e L 2 + η e ( M 0,e -M 0,c ), V e L 2 -ıqF (t)D e M 0,e n, ∇V e L 2 + ıqF (t)D e ∇ M 0,e • n, V e L 2 + q 2 F (t) 2 D e n • n M 0,e , V e L 2 = 0 in D ′ (]0, T [), ∂ t M 0,c , V c L 2 + η c ( M 0,c -M 0,e ), V c L 2 = 0 in D ′ (]0, T [),
We introduce for almost every t ∈]0, T [, the bilinear form a(•; •, •) defined, for every (U, V ) ∈ X 2 by

a(t; U, V ) = (D e ∇U e , ∇V e ) L 2 + (η e (U e -U c ), V e ) L 2 -(ıqF (t)D e U e n, ∇V e ) L 2 + (ıqF (t)D e ∇U e • n, V e ) L 2 + q 2 F (t) 2 D e n • nU e , V e L 2 + (η c (U c -U e , V c ) L 2 . As f ∈ L ∞ (]0, T [) (and consequently F ∈ L ∞ (]0, T [)), for every (U, V ) ∈ X 2 , t -→ a(t; U, V ) is measurable.
Moreover, using the Cauchy-Schwarz inequality,

|a(t; U, V )| ≤ ||D e || ||∇U e || L 2 ||∇V e || L 2 + ||D e ||q||F || L ∞ (]0,T [) ||U e || L 2 ||∇V e || L 2 +||D e ||q||F || L ∞ (]0,T [) ||∇U e || L 2 ||V e || L 2 ||D e ||q 2 ||F || 2 L ∞ (]0,T [) ||U e || L 2 ||V e || L 2 +η e (||U e || L 2 + ||U c || L 2 ) ||V e || L 2 + η c (||U e || L 2 + ||U c || L 2 ) ||V c || L 2 ,
which globally gives

|a(t; U, V )| ≤ C(η e , η c , ||D e ||, q, F )||U || X ||V || X ,
for some constant C(η e , η c , ||D e ||, q, F ) > 0 where X is endowed with its natural Hilbert space norm

||U || 2 X := ||U e || 2 H 1 + ||U c || 2 L 2 . Next, we have a(t; U, U ) = (D e ∇U e , ∇U e ) L 2 -ıqF (t) (U e D e n, ∇U e ) L 2 + ıqF (t) (D e ∇U e • n, U e ) L 2 +q 2 F (t) 2 D e n • n (U, U ) L 2 + η e (U e -U c , U e ) L 2 + η c (U c -U e , U c ) L 2 .
Then,

ℜ(a(t; U, U )) ≥ σ -||∇U e || 2 L 2 + q 2 F (t) 2 σ -||U e || 2 L 2 -2q|F (t)|||D e ||||∇U || L 2 ||U || L 2 +η e (U e -U c , U e ) L 2 + η c (U c -U e , U c ) L 2 ,
where σ -is such that D e ξ • ξ ≥ σ -||ξ|| 2 for any ξ ∈ R d , which exists as D e is positive definite. Using Young's inequality, we get

ℜ(a(t; U, U )) ≥ σ - 2 ||∇U e || 2 L 2 + q 2 F (t) 2 σ --2||De|| 2 σ - ||U e || 2 L 2 + η e -ηe+ηc 2 ||U e || 2 L 2 + η c -ηe+ηc 2 ||U c || 2 L 2 .
Then, for any λ > 0, such that for some fixed δ > 0 and almost every t ∈]0, T [,

λ + q 2 F (t) 2 σ -- 2||D e || 2 σ - - η e + η c 2 ≥ δ > 0, (such a λ exists as f ∈ L ∞ (]0, T [), ℜ(a(t; U, U )) + λ||U || 2 L 2 ≥ min σ - 2 , δ ||U || 2 X .
Then, Lion's theorem (see for instance [14, Chapter XVIII]) ensures the existence and uniqueness of the solution of the variational problem

d dt < M (t), V > X ′ ,X +a(t; M (t), V ) = 0 in D ′ (]0, T [), ∀ V ∈ X, M (0) = M init .
(3.28)

Now, in order to get estimates (3.26) and (3.27), we use V = M (t) as a test function, which is possible since M ∈ W (0, T, X). We integrate over time, and get

1 2 || M 0,e (t)|| 2 L 2 + 1 2 || M 0,c (t)|| 2 L 2 -||M init || 2 L 2 + t 0 a(s; M (s), M (s))ds = 0.
Proceeding as the proof of coercivity and using Young's inequality,

1 2 || M 0,e (t)|| 2 L 2 + 1 2 || M 0,c (t)|| 2 L 2 + σ - 2 t 0 ||∇ M 0,e (s)|| 2 L 2 ds + q 2 σ - t 0 F (s) 2 || M 0,e (s)|| 2 L 2 ds +η e t 0 || M 0,e (s)|| 2 L 2 ds + η c t 0 || M 0,c (s)|| 2 L 2 ds ≤ ||M init || 2 L 2 + 2||D e || 2 q 2 σ - t 0 F (s) 2 || M 0,e (s)|| 2 L 2 ds + η e + η c 2 t 0 || M 0,e (s)|| 2 L 2 ds + η e + η c 2 t 0 || M 0,c (s)|| 2 L 2 ds. (3.29) In particular, || M (t)|| 2 H ≤ 2||M init || 2 L 2 + t 0 4||D e || 2 q 2 σ - F (s) 2 + (η e + η c ) || M (s)|| 2 H ds.
Then, using Gronwall's lemma,

|| M (t)|| H ≤ √ 2||M init || L 2 exp 2||D e || 2 q 2 σ - t 0 F (s) 2 ds + (η e + η c )t ≤ 2||M init || L 2 exp 2||D e || 2 q 2 ||F || 2 L ∞ (]0,T [) σ - + η e + η c t,
which implies (3.26). To prove (3.27), we notice that (3.29) yields

t 0 ||∇ M (s)|| 2 H ds ≤ 2 σ -||M init || 2 L 2 + t 0 4||D e || 2 q 2 σ - F (s) 2 + (η e + η c ) || M (s)|| 2 H ds.
Using (3.26), we deduce

t 0 ||∇ M (s)|| 2 H ds ≤ 2 σ -||M init || 2 L 2 1 + exp 4||D e || 2 q 2 ||F || 2 L ∞ (]0,T [) σ - + η e + η c t ,
which corresponds with(3.27). Finally, using the well-known injection W (0, T, H 1 ) ֒→ C 0 (0, T, L 2 ), it is now a classical exercise to prove rigorously the equivalence between problems (3.28) and (3.25). Now, the multiplication by e ıqn•xF (t) and its inverse e -ıqn•xF (t) clearly defines a bijection between W (0, T, X) with its natural norm

||U || 2 W (0,T,X) := ||U || 2 L 2 (0,T,X) + ||∂ t U || 2 L 2 (0,T,X ′ )
and the space

W q,f (0, T, X) = U ∈ L 2 (0, T, X) | ∂ q,f t U ∈ L 2 (0, T, X ′ ) ,
where ∂ q,f t denotes the operator ∂ q,f t = ∂ t + ıqn • xf (t), if we endow W q,f (0, T, X) with the norm

||U || 2 W q,f (0,T,X) = ||U || 2 L 2 (0,T,X) + ||∂ q,f t U || 2 L 2 (0,T,X ′ ) .
Thus, the previous result dictates the following notion of solution for problem (3.20):

A variational solution of problem (3.20) is an element M 0 = (M 0,e , M 0,c ) of W q,f (0, T, X) ∩ C 0 (0, T, L 2 (R d )) 2 such that, for any V ∈ X, we have

∂ q,f t M 0,e , V e H -1 ,H 1 + (D e ∇M 0,e , V e ) L 2 (R d ) + (η e (M 0,e -M 0,c ), V ) L 2 (R d ) = 0 in D ′ (]0, T [), ∂ q,f t M 0,c , V c L 2 + (η c (M 0,c -M 0,e ), V ) L 2 (R d ) = 0 in D ′ (]0, T [), M 0,e (0) = M init and M 0,c (0) = M init in R d . (3.30) We immediately deduce from Theorem 1 . Theorem 2. Let f ∈ L ∞ (0, T ) and M init ∈ L 2 (Ω).
There exists a unique variational solution M = (M 0,e , M 0,c ) ∈ W q,f (0, T, X) ∩ C 0 (0, T, L 2 (R d )) 2 to problem (3.20).

If Ω is bounded, then the previous result also holds true if we replace H 1 (Ω) by H 1 0 (Ω) and we supplement equations (3.20) by homogeneous Dirichlet boundary conditions for M 0,e on ∂Ω. In that case, the proof can be be conducted directly on M 0 , as the problem on a bounded domain does not involve unbounded coefficients anymore.

Remark 1. In the cases where D c is also positive definite matrix (which corresponds to the presence of elongated cells or axons), one has to rather consider system (3.19) and Theorem 2 holds true in this case with

H = L 2 (Ω) × L 2 (Ω) and X = H 1 (Ω) × H 1 (Ω).
From the solution to (3.20) one can built an approximation of the magnetization as

M ε ≈ |Y e | |Y | M 0,e + |Y c | |Y | M 0,c
This formally account for

M ε ≃ M 0,e in Ω ε e and M ε ≃ M 0,c in Ω ε c .
More precisely one can prove that under the hypothesis of Theorem 2,

M ε ⇀ |Y e | |Y | M 0,e + |Y c | |Y | M 0,c weakly in L 2 (0, T, L 2 (R d )). (3.31)
The proof is given in [START_REF] Coatléven | Mathematical justification of macroscopic models for diffusion mri through the periodic unfolding method[END_REF] and uses the periodic unfolding method [START_REF] Cioranescu | The periodic unfolding method in homogenization[END_REF], extended to the time dependent cases.

4. An ODE model for the dMRI signal. In this section, we give more details about the measured signal in a dMRI experiment and we explain how to use the model introduced in the previous section to obtain an approximation to the signal in a simple manner. In practice, the diffusion-encoding magnetic field gradient has a time profile that satisfies, at the echo time T E ,

T E 0 f (t)dt = 0, ( 4.1)
and the signal is measured at t = T E . The measured signal is given by

S ε (q, n) = R d M ε (x, T E ) dx.
Then, it is natural to define the approximate signal by:

S(q, n) = R d M e (x, T E ) + M c (x, T E ) dx,
where we have set

m 0 e = R d M e (x, 0) dx, m 0 c = R d M c (x, 0) dx and m 0 = R d M e (x, 0) + M c (x, 0) dx.
These initial moments will naturally appear in the following and are of great interest for the inverse problem since

θ e = m 0 e m 0 = |Y e | |Y | and θ c = m 0 c m 0 = |Y c | |Y |
represents respectively the extra-cellular volume fraction and the intra-cellular volume fraction. We now establish that the signal obtained from solving (3.20) can be equivalently obtained by solving a simpler system of coupled ODEs. Theorem 3. Under the hypothesis of Theorem 2, if we assume in addition that f satisfies (4.1) and assume that the Fourier transform of M init is continuous in the neighborhood of the origin, then the signal S(q, n) is well defined and can computed as

S(q, n) = m(T E ),
where σ e := D e n • n.

Proof. Let (M 0,e , M 0,c ) be the unique solution of (3.20) in C 0 (0, T E , H) ∩ W q,f (0, T E , X). We set again Then, as we have seen before, we obtain the system of equations for M e and M c ,

∂ M e ∂t -div(D e ∇ M e ) + 2ıqD e ∇ M e • n t 0 f (s)ds +q 2 D e n • n t 0 f (s)ds 2 M e + η e M e -η c M c = 0 in Ω×]0, T E [, ∂ M c ∂t + η c M c -η e M e = 0 in Ω×]0, T E [,
(note that we used here the identity η e θ e = η c θ c ) with the initial conditions

M e (•, 0) = θ e M init , M c (•, 0) = θ c M init in Ω.
We denote by M e and M c the Fourier transform with respect to the space variable of respectively M e and M c and denote the dual variable ξ. Applying the Fourier transform to the system of equations satisfied by M e and M c yields

∂M e ∂t + D e ξ • ξM e -2qD e ξ • nM e t 0 f (s)ds +q 2 D e n • n t 0 f (s)ds 2 M e + η e M e -η c M c = 0 in Ω×]0, T E [, ∂M c ∂t + η c M c -η e M e = 0 in Ω×]0, T E [.
From our hypothesis, M e (•, 0) and M c (•, 0) are continuous functions. We therefore deduce from classical standard results on parametric first order linear differential systems that M e and M c belong to C 1 ([0, T E ], C 0 (Ω)). We set m e (t) := M e (0, t) and m c := M c (0, t), which can be formally written as

m e (t) = R d M e (x, t) dx and m c (t) = R d M c (x, t) dx (4.4)
where the integrals has to be understood as duality pairing in the distributional sense. We observe that, for ξ = 0, the above system of equations directly implies

dm e dt + q 2 σ e t 0 f (s)ds 2 m e + η e m e -η c m c = 0, dm c dt + η c m c -η e m e = 0,
with σ e = D e n • n. Finally, for the signal we have for α = e, c, using

T E 0 f (s)ds = 0, R d M α (x, T E )dx = R d M α (x, T E )e -ıqn•x T E 0 f (s)ds dx = R d M α (x, T E ) = m α (T E )
which concludes the proof. Our ODE model is reminiscent of the so-called Karger model [START_REF] Karger | Principles and application of self-diffusion measurements by nuclear magnetic resonance[END_REF]. In fact, the latter can be seen as a special case of our model when the pulse duration δ is very small compared to the interval between the two pulses, in other words, δ ≪ ∆. However, this assumption is not always true in dMRI experiments due to safety concerns and/or hardware limitations. In addition, the Karger model is obtained heuristically, whereas we derived the ODE model of (4.3) more rigorously from the microscopic Bloch-Torrey PDE.

Remark 2. Indeed, for the general case where D c = 0, the system (4.3) has to be replaced with This allows us to obtain an expansion of the measured signal with respect to z as

dm e dt + q 2 σ e t 0 f (
S( √ z, n) = m 0 - +∞ n=1 σ e z n T E 0 t 0 f 2 m e,n-1 (t)dt. (4.8)
It is difficult to use this expression to get explicit links between measured signals and tissue properties in the general case. For instance, the first term of this expansion leads to, using

S(0, n) = m 0 , log S(q, n) S(0, n) = - σ e q 2 m 0 T E 0 t 0 f 2 m e,0 (t)dt + O(q 4
).

If we use the relation η e m 0 eη c m 0 c = 0 (verified by the homogenized coefficients due to (3.18)), then we simply have log S(q, n)

S(0, n) = -σ e θ e q 2 T E 0 t 0 f 2 + O(q 4 ).
The product σ e θ e can then be estimated directly by computing the slope of the log of the measured signal. The higher order terms give much more involved relations and we do not pursue this.

4.3. Numerical validation of the ODE model. Now we will compare our macroscopic model to the two compartments microscopic model, using very simple geometries and biologically reasonable parameters for the intrinsic diffusion coefficients and the membrane permeability. Simulations on more complex cellular geometries will be the subject of a future paper.

We will take the commonly used PGSE diffusion-encoding sequence in (2.1). By construction, we have:

T E 0 f (t)dt = 0,
and for simplicity we will take T E = δ + ∆.

The dMRI signal can be measured for several values of q, δ, ∆ and directions n. We will consequently denote by S(q, δ, ∆, n) the signal to emphasize this dependency and also normalize this signal by dividing it by R d M init dx. It is common in the dMRI community not to display S(q, δ, ∆, n) as a function of q, but as a function of the so called b-value:

b(q) := q 2 ∆+δ 0 t 0 f (s)ds 2 dt = q 2 δ 2 ∆ - δ 3 .
To understand why this quantity is used, notice that in the case of a single Bloch-Torrey equation in R d with a constant diffusion tensor D,

∂M ∂t + iqn • xf (t)M -div x (D∇ x M ) = 0 in R d ×]0, T E [, M (•, 0) = M init in R d ,
the Fourier transform M is the solution of:

∂ M ∂t -f (t)qn • ∇ ξ M + ξ T Dξ M = 0 in R d ×]0, T E [, M (•, 0) = Minit in R d .
This problem can easily be solved using the method of characteristics, and we obtain:

M (ξ, t) = Minit ξ + qn t 0 f (s)ds exp - t 0 D ξ + qn t s f (µ)dµ • ξ + qn t s f (µ)dµ ds .
Then, since R d M (x, ∆ + δ)dx = M (0, ∆ + δ), we can get, in this particular case, an expression for the corresponding normalized signal:

S(q, δ, ∆, n) := ( R d M (x, ∆ + δ)dx)/( R d M init (x)dx) = M (0, ∆ + δ)/ Minit (0),
that simplifies to, using

∆+δ 0 f (s)ds = 0, S(q, δ, ∆, n) = exp -q 2 Dn • n ∆+δ 0 s 0 f (µ)dµ 2 ds = exp (-(Dn • n)b(q)) .
We then observe that log( S(q, δ, ∆, n)) is a linear function of the b-value and the slope of this linear function is -(Dn • n), which gives access to the coefficients of the tensor D by varying n. We now provide some numerical results using biologically reasonable values for the intrinsic diffusion coefficients and the membrane permeability on some simple geometries. We compare the signal obtained using the system of ODEs (4.3) to the signal obtained by solving the two-compartment periodic PDE model (2.4).

Since (2.4) is a problem on R d , we cannot numerically solve it for arbitrary initial condition. For the numerical simulations we will rather consider the case of constant M init . It can be easily shown that a solution of (2.4) is quasi-periodic in space, with period ε and quasi-periodic coefficient exp(-ıεq i t 0 f (s)ds) in each direction e i . Thus, if we extend the definition of a normalized signal as

S ε = lim N →∞ Ω N M ε (x, T E )dx Ω N M init (x)dx ,
where Ω N denotes the union of N cells, one gets, for constant M init ,

S ε = 1 M init |Y ε | Yε M ε (x, T E )dx,
where Y ε denotes a given periodicity cell. Then one can reduce the computations to only one periodicity cell with quasiperiodidicty condition. For the homogenized model we indeed equivalently extend the definition of the normalized signal S similarly to S ε by replacing M ε with the homogenized field M . Even if the case of constant initial conditions does not enter into the previous theoretical framework (initial conditions is not L 2 (R d )) one can easily verify that if we denote (m e , m c ) the solution of (4.3) associated with the initial conditions: 2), with (M 0,e , M 0,c ), solution of (3.20). One then easily check that the normalized signal is still given by S(q, δ,

∆, n) = m e (T E ) + m c (T E ).
We test two situations in dimension 2, for simplicity. One is a single circular biological cell of radius R m placed at the center of the periodicity cell, for several values of ∆, in the gradient direction e x . The second one is the case of several circular biological cells of variable radii, in two gradient directions, e x and e y . We solve the ODE model (4.3), along with the periodic reference model (2.4) on a very refined mesh. Both computations are performed using FreeFem++ [START_REF] Pironneau | freefem++[END_REF].

A comparison of the signals obtained from the exact and approximate models is displayed on Figures 4.1 and 4.2, where the values of the model parameters are indicated. These values are chosen close to the values used often in the literature [START_REF] Xu | Numerical study of water diffusion in biological tissues using an improved finite difference method[END_REF][START_REF] Harkins | Assessment of the effects of cellular tissue properties on adc measurements by numerical simulation of water diffusion[END_REF] for dMRI numerical simulations. We can see on these figures that the ODE model (4.3) provides an excellent approximation when modeling the dMRI signal at all four diffusion times: ∆ = 5, 15, 25, 35 ms. In particular, this approximate model accurately reproduces the 'curvature' of the signal: the obtained signals with the two models are indeed not a single exponential (the log is not a straight line), as it would be in the case of single Bloch-Torrey equation without membranes. This phenomenon, well-known in the medical imaging community, can be reinterpreted from the ODE model: it is the fact that we have two coupled equations that induces this memory-like effect (see the end of Section 3).

5.

The inverse problem : retrieving macroscopic properties of tissues from dMRI measurements. Now that we have checked that our ODE model accurately reproduces the dMRI signal in the example geometries, the natural idea is to check whether we can retrieve the coefficients of our ODE model from the measured signals, i.e. if we can solve the inverse problem of finding the unknown coefficients β = (η e , η c , D e , m e,0 , m c,0 ). These coefficients are of great practical importance. Indeed, from initial moments, we can recover the cellular volume fractions, thus giving information on the concentration of cells, or their potential swelling. The coupling coefficients are also a way to obtain information on the permeability of cellular membranes, which is a very difficult quantity to measure in practice, while the homogenized diffusion tensor may give some information on the orientation distribution of elongated cells, for instance.

Let (S i ) 1≤i≤Nexp be a set of measures. For instance, we can take

S i = S(q i , δ i , ∆ i , n i )
where S is the signal from the homogenized model, or

S i = S ε (q i , δ i , ∆ i , n i )
where we use measures obtained from the two compartment microscopic model.

A natural way to formulate the inverse problem is to solve the least squares problem: find

G(β) = Nexp i=1 |S(q i , δ i , ∆ i , n i ) -S i | 2 .
This problem can be tackled numerically using one of the many methods in the literature for solving inverse problems using least squares. As our purpose here is to illustrate the feasibility of this inverse problem, we have chosen to use the well tested interior-reflective Newton method described in [START_REF]An interior, trust region approach for nonlinear minimization subject to bounds[END_REF] and [START_REF] Coleman | On the convergence of reflective newton methods for large-scale nonlinear minimization subject to bounds[END_REF] and implemented in Matlab under the name lsqnonlin. To generate the measures, we have used the following values for the physical parameters, σ e = 3 × 10 -3 mm 2 /s, σ c = 1.6 × 10 -3 mm 2 /s, κ = 5 × 10 -5 m/s, and computed the signals from the models (4.3) and (2.4) in the directions e x and e y , for q i = 0.000015×i×10 6 (mm×ms) -1 , i = 0, 1, • • • , 25, and for (δ, ∆) = (3.5, 5), (3.25, 7.5), [START_REF] Aslund | Filter-exchange pgse nmr determination of cell membrane permeability[END_REF][START_REF] Clark | Water diffusion compartmentation and anisotropy at high b values in the human brain[END_REF], (2.75, 12.5)ms. We now provide results for the two examples described in the previous section.

First, we consider the case of a single circular biological cell of radius R m = 0.49µm. To solve our minimization problem, we use a random starting set of parameters β 0 that are at most 50% different in the l ∞ norm from the exact set of parameters β * we are looking for. We provide results for several values of β 0 on Table 5.1 (we display (η e , η c , θ c , D e (0, 0), D e (0, 1), D e (1, 0), D e (1, 1)) in each entry of the table, as it is equivalent in our case to knowing θ c and m 0 = S(0, δ i , ∆ i , n i ) or m e,0 and m c,0 ). We see from these results that, as can be expected, the estimated parameters are closer to the real, homogenized values, when we use the signals from model (4.3) than when we use the nonhomogenized signal coming from model (2.4). Indeed, the modeling error, corresponding to the fact that (4.3) is an approximation of (2.4) and thus the signals coming from these two models, while being close, are not equal, can be reinterpreted as (2.4) providing "noisy" data when we try to solve the inverse problem for (4.3).

T ype of data

It can also be observed on Table 5.1 that the quality of the estimated coefficients seems to depend on our initial guess. To limit the effect of this arbitrary choice, we have solved the inverse problem for several initial guesses, and took the values of the estimated coefficients that correspond to the smallest residual. Results are displayed on Table 5.2. We see on Table 5.2 that the results are in this case more robust, and that we obtained a quite good approximation of our homogenized coefficients. To further improve the results, we notice that: Thus, we can eliminate m c,0 in our inversion process. We reproduce the above study after this elimination. We display on Table 5.3 the inversion from one randomly chosen set of initial guesses, while on Table 5.4 we again display the values of the estimated parameters that correspond to the smallest residual from several randomly chosen initial guesses. We see on these two Tables that the estimated parameters are much closer to the exact homogenized values when we use this constraint, even when using a single randomly chosen initial guess (the results being of course still better when Numerical inversion for the case of a single circular biological cell, for two randomly chosen inital sets β 0 , using the constraint θc = ηe ηe+ηc .

η e η c = θ c θ e = θ c 1 -θ c and m c,0 = θ c m 0 , T ype
we use several initial guesses). Thus, it seems that adding this constraint can improve parameter estimation when solving the inverse problem. Finally, we display on Table 5.5 the results obtained for the periodicity cell of Figure 4.2, where there are 4 circular biological cells of different sizes, with the parameters : σ e = 3 × 10 -3 mm 2 /s, σ c = 1.6 × 10 -3 mm 2 /s, κ = 5 × 10 -5 m/s, and obtained the dMRI signal for the same input parameters, q, n, δ, ∆, as previously. We use the aforementioned constraint and several initial guesses, taking again the results corresponding to the smallest residual. We see that even in this more complicated situation, the inverse problem can still be solved to give good estimates of the model parameters.

6. Discussion of two assumptions made in the derivation of the macroscopic model. Now we discuss two assumptions we made on the volume Ω for which we formulated our macroscopic model in the context of biological tissue dMRI. One assumption is that Ω is periodic with a period ε that is small compared to the size of Ω. The second is that Ω can be taken to be R d . The first assumption was made so that we could apply periodic homogenization theory, which is technically 0.0017 Table 5.5 Numerical inversion for the case of the periodicity cell of Figure 4.2, using several intial sets β 0 , under the constraint θc = ηe ηe+ηc simpler than the generic case. The second assumption was made purely for the simplicity of presentation, and for obtaining an ODE system governing the dMRI signal of the macroscopic PDE model. We discuss them in turn.

In this paper, we applied periodic homogenization theory, which technically simpler than the generic (non-periodic) case. Clearly, biological tissue is not periodic with any period, however, it is known in some similar contexts (see for instance [START_REF] Arbogast | Gravitational forces in dual-porosity systems: I. model derivation by homogenization[END_REF] and [START_REF]Gravitational forces in dual-porosity systems: Ii. computational validation of the homogenized model[END_REF] for the case of porous media) that even for media that are not truly periodic, the homogenized limit obtained in the periodic case remains formally relevant for describing generic (non periodic) cases, i.e. the obtained macroscopic equations for both configurations have the same analytical form (only the exact values of the coefficients would change). We conjecture that this is also true for our model but showing this is beyond the scope of this paper. We just note that the pheonomenologically derived macroscopic Karger model has the same form (but different coefficients) as the macroscopic signal model we derive here, and has been widely used by MR physicists already in tissue dMRI (references given in Section 1). In fact, a validation of the Karger model by numerical simulations on a computational domain that is not periodic was performed in [START_REF] Fieremans | Monte carlo study of a two-compartment exchange model of diffusion[END_REF]. Now, concerning taking Ω to be R d , what is meant by this, in the case when Ω is not itself periodic, is actually the periodic extension of Ω to R d . To obtain the dMRI signal in a specific voxel V one should sum the water proton magnetization inside V at the echo time. Thus it suffices to consider the volume Ω to contain V and a buffer region around V . The buffer region contains parts of the voxels close to V . The width of the buffer region should be a small multiple of the diffusion distance. In this way, one would have accounted for essentially all the water molecules that enter, exit, or stay in V during the diffusion time. The spatial integral for the dMRI signal should be taken over V only. Then the boundary condition on ∂Ω: Dirichlet, Neumann, periodic, etc., has insignificant effect on the dMRI signal because we created a large enough buffer around V . At physically realistic dMRI diffusion times of 10-100 ms, the average diffusion distance is, assuming an average diffusivity of 10 -3 mm 2 /s [START_REF] Niendorf | Biexponential diffusion attenuation in various states of brain tissue: Implications for diffusionweighted imaging[END_REF], between 8-25 µm. Thus, the width of the buffer region is much smaller than the side lengths of a voxel (around 1 mm), and there is very little influence on the dMRI signal in V from adjacent voxels, all the more so if there is not a large tissue variation between adjacent voxels. One can then simply place copies of V next to it rather than the true neighboring voxels to get the signal in V , meaning: Ω can be taken to be exactly V , and periodic conditions can be put on V . The error associated with a periodic extension of the voxel (instead of using the true neighboring voxels) will be small. We illustrate this point below with a numerical example. Very importantly, the choice of periodic boundary conditions on Ω = V has the advantage of enabling the formulation of a system of ODEs for dMRI signal, which is not the case if we imposed Dirichlet or Neumann boundary conditions on Ω ⊃ V .

Let V = [-500µm, 500µm] 2 be a voxel, and let B = [-750µm, 750µm] 2 contain V and a buffer region around V . We placed 20000 circular biological cells in B whose centers are uniformly randomly distributed in B. The radii of the circles were chosen randomly with uniform distribution between 1 and 5 µm. Thus, B itself is not a periodic domain. See Figure 6.1(a). We computed the dMRI signal in V for the PGSE sequence, δ = 80 ms and ∆ = 80 ms, in two different ways. One is solving the microscopic Bloch-Torrey PDE in B, with periodic conditions on ∂B. The second is solving the microscopic PDE in V , with periodic conditions on ∂V . The spatial integral of the PDE solutions was taken over V in both cases to produce the dMRI signal in V . The intrinsic diffusion coefficient is 3 × 10 -3 mm 2 /s inside and outside the biological cells. The membrane permeability is κ = 10 -5 m/s. In Figure 6.1(b) we see that the two dMRI signals are indistinguishable at the diffusion time simulated. The diffusion distance is bounded by 4 × 3 × 10 -3 mm 2 /s × 160ms ≈ 44µm, which is small compared to 1 mm, the side length of V . This gives numerical support to our claim that even though it is more strictly correct to use the domain B and sum the magnetization over V , the numerical difference compared with using the domain V with periodic conditions on ∂V is insignificant, due to the fact that the diffusion distance is much smaller than the size of the voxel in tissue dMRI.

Conclusion.

We have formulated a new macroscopic model for the complex transverse water proton magnetization in a voxel due to diffusion-encoding magnetic field gradient pulses for diffusion MRI in biological tissue. This new model, derived through periodic homogenization of the two-compartment Bloch-Torrey equation, reproduces the memory effects commonly observed experimentally, and explains it as the influence of cellular membranes on the diffusion of water. In addition, by choosing periodic boundary conditions on the boundary of the voxel, we showed that the dMRI signal can be obtained after solving a system of ODEs. The advantage of formulating the ODE model for the dMRI signal is that quantitative macroscopic information on the probed tissue can be easily obtained through solving the inverse problem using the ODE model. We provided a preliminary numerical investigation of the inverse problem and showed the feasibility and relative accuracy of this procedure for synthetically generated dMRI signal data.

Fig. 2 . 1 .

 21 Fig. 2.1. Schematic of a biological cell, with its membrane (left) and an idealized two-compartment cell (right).

Fig. 2 . 2 .

 22 Fig. 2.2. A periodicity cell Y , containing simplified biological cells with an infinitely thin interface for membraneNow, we denote

  Yc div y σ c e i dy -Γm σ c e i • νds = Yc div y σ e e i dy + Γm σ e e i • νds = 0, (3.12)

M

  e := |Y e | |Y | M 0,e and M c := |Y c | |Y | M 0,c . (4.2) 4.1. A coupled ODE model for dMRI's signals. Let us denote

e dt + q 2 σ e t 0 f (s)ds 2 m

 02 where m = m e + m c , and (m e , m c ) ∈ C 1 (0, T E ) 2 is the unique solution of dm e + η e m eη c m c = 0, dm c dt + η c m cη e m e = 0, m e (0) = m 0 e and m c (0) = m 0 c ,

M

  e (x, t) = M e (x, t)e ıqn•x t 0 f (s)ds and M c (x, t) = M c (x, t)e ıqn•x t 0 f (s)ds .

s)ds 2 mq 2 σ c t 0 f (s)ds 2 m( 4 . 5 ) 4 . 2 . 0 f (s)ds 2 m( 4 . 6 ) 2 ηe η e η e + η c t 0 s 0 f 2 1

 2024542024622 e + η e m eη c m c = 0, dm c dt + c + η c m cη e m e = 0, m e (0) = m 0 e and m c (0) = m 0 c , where σ c = D c n • n.Analytic expansion of the ODE model solution. We extend the definition of our ODE model for any complex number z as dm e (z, •) dt (t) + zσ e t e (z, t) + η e m e (z, t) -η c m c (z, t) = 0, dm c (z, •) dt (t) + η c m c (z, t)η e m e (z, t) = 0,m e (z, 0) = m 0 e and m c (z, 0) = m 0 c . The solutions (m e (z), m c (z)) belong to C 1 (0, T E ) and thus to L ∞ (0, T E ). It is quite easy to see that (m e (•, t), m c (•, t)) can be differentiated with respect to z for any z ∈ C. Then m e (•, t) and m c (•, t) are holomorphic on all C and admit an expansion of the form m e (z, t) = +∞ n=0 z n m e,n (t) and m c (z, t) = +∞ n=0 z n m c,n (t). Inserting these expansions in (4.3), straightforward calculations lead to m e,0 (t) = η c (m 0 e + m 0 c ) η e + η c + η e m 0 eη c m 0 c η e + η c e -(ηe+ηc)t , m c,0 (t) = η e (m 0 e + m 0 c ) η e + η c -η e m 0 eη c m 0 c η e + η c e -(ηe+ηc)t , m e,n (t) = -σ e η e + η c t e e -(ηe+ηc)(s+t) + η c m e,n-1 (s)ds, m c,n (t) = -σe -(ηe+ηc)(s+t) m e,n-1 (s)ds. (4.7)

m-ıq•x t 0 f-ıq•x t 0 f

 00 e (0) = θ e and m c (0) = θ c , then M e (x, t) = M init m e (t)e (s)ds and M c (x, t) = M init m c (t)e (s)ds correspond, through (4.

Fig. 4 . 1 .

 41 Fig. 4.1. Normalized signals for a single circular biological cell of radius Rm = 0.49µm in a periodicity cell of 1µm 3 , with σe = 3 × 10 -3 mm 2 /s, σc = 1.6 × 10 -3 mm 2 /s, κ = 5 × 10 -5 m/s, δ = 3.5ms, n = ex, for various ∆.

Fig. 4 . 2 .

 42 Fig. 4.2. Normalized signals and the periodicity cell, with σe = 3 × 10 -3 mm 2 /s, σc = 1.6 × 10 -3 mm 2 /s, κ = 5 × 10 -5 m/s, δ = 3.5ms, ∆ = 5ms, n = ex and ey.

  and consequently:m c,0 = η e m 0 η e + η c .

  Fig. 6.1. (a) The voxel V = [-500µm, 500µm] 2 and a larger domain B = [-750µm, 750µm] 2 containing V and a buffer region around V . There are 20000 circular biological cells in B whose centers are uniformly randomly distributed in B. The radii of the circles were chosen randomly with uniform distribution between 1 and 5 µm. (b) The dMRI signals in V obtained by solving the microscopic PDE in B and in V subject to periodic boundary conditions are indisguishable. The signals are computed by integrating the PDE solutions over V . The intrinsic diffusion coefficient is 3 × 10 -3 mm 2 /s both inside and outside the biological cells. The membrane permeability is κ = 10 -5 m/s. The gradient direction is ex.

Table 5 . 1

 51 Numerical inversion for the case of a single circular biological cell, for two randomly chosen inital sets β 0

Table 5 . 2

 52 Numerical inversion for the case of a single circular biological cell, using several initial sets β 0 .

	of data	Initial guesses Homogenization	Inversion
			0.00062	0.00085
			0.00020	0.00023
			0.75	0.60
	N on-homogenized signal S ε	10	0.00090 7.39 10 -9	0.00117 7.28 10 -9
			7.39 10 -9	10.27 10 -9
			0.00090	0.00117
			0.00062	0.00064
			0.00020	0.00022
			0.75	0.69
	N on-homogenized signal S ε	15	0.00090 7.39 10 -9	0.00097 7.81 10 -9
			7.39 10 -9	6.01 10 -9
			0.00090	0.00097
			0.00062	0.00063
			0.00020	0.00022
			0.75	0.70
	N on-homogenized signal S ε	20	0.00090 7.39 10 -9	0.00096 6.88 10 -9
			7.39 10 -9	11.06 10 -9
			0.00090	0.00096364
			0.00062	0.00081
			0.00020	0.00022
			0.75	0.95
	Homogenized signal S	10	0.00090 7.39 10 -9	0.00106 7.37 10 -9
			7.39 10 -9	10.35 10 -9
			0.00090	0.00106
			0.00062	0.00063
			0.00020	0.00020
			0.75	0.74
	Homogenized signal S	15	0.00090 7.39 10 -9	0.00091 9.36 10 -9
			7.39 10 -9	9.96 10 -9
			0.00090	0.00091

Table 5 . 3

 53

	T ype of data	Initial guess Homogenization	Inversion
		0.00060	0.00062	0.00056
		0.00030	0.00020	0.00022
			0.75	0.72
	N on-homogenized signal S ε	0.00059 4.52 10 -9	0.00090 3.34 10 -9	0.000897 4.52 10 -9
		3.82 10 -9	3.34 10 -9	3.82 10 -9
		0.00079	0.00090	0.000897
		0.00043	0.00062	0.00052
		0.00019	0.00020	0.00021
			0.75	0.71
	N on-homogenized signal S ε	0.000887 2.07 10 -9	0.00090 3.34 10 -9	0.00087 2.07 10 -9
		3.64 10 -9	3.34 10 -9	3.64 10 -9
		0.00066	0.00090	0.00087
		0.00074	0.00062	0.00070
		0.00020	0.00020	0.00021
			0.75	0.77
	Homogenized signal S	0.00087 3.40 10 -9	0.00090 3.34 10 -9	0.00096 3.99 10 -9
		2.5 10 -9	3.34 10 -9	2.5031 10 -9
		0.00071	0.00090	0.00096
		0.00032	0.00062	0.00055
		0.00030	0.00020	0.00020
			0.75	0.74
	Homogenized signal S	0.00060 2.02 10 -9	0.00090 3.34 10 -9	0.00084 2.02 10 -9
		2.91 10 -9	3.34 10 -9	2.91 10 -9
		0.00063	0.00090	0.00084

Table 5 . 4

 54 Numerical inversion for the case of a single circular biological cell, using several intial sets β 0 , under the constraint θc = ηe ηe+ηc .

	T ype of data	Initial guesses Homogenization	Inversion
			0.00062	0.00055
			0.00020	0.00022
			0.75	0.72
	N on-homogenized signal S ε	10	0.00090 7.39 10 -9	0.00089 8.13 10 -9
			7.39 10 -9	4.97 10 -9
			0.00090	0.00089
			0.00062	0.00059
			0.00020	0.00020
			0.75	0.745
	Homogenized signal S	10	0.00090 7.39 10 -9	0.00087 10.35 10 -9
			7.39 10 -9	10.63 10 -9
			0.00090	0.00087
	T ype of data	Initial guesses Homogenization Inversion
			0.00045	0.00032
			0.00046	0.00038
			0.51	0.45
	N on-homogenized signal S ε	10	0.0019 7.78 10 -6	0.0018 8.18 10 -6
			7.78 10 -6	4.15 10 -6
			0.00175	0.00159
			0.00048	0.00046
			0.00046	0.00045
			0.51	0.50
	Homogenized signal S	10	0.0019 7.78 10 -6	0.0019 9.55 10 -6
			7.78 10 -6	6.64 10 -6
			0.00175

We only considered the case where the biological cells are of a small size (which corresponds to D c = 0 in our model). Numerical investigations of this approach for more complex cell geometries (including elongated cells) and three dimensional configurations will be the subject of a future work. While the model we develop in this paper supposes Gaussian diffusion in each compartment, one possible way of generalizing this model is to make the diffusion in each compartment non-Gaussian.

On the mathematical side, we only provided a formal justification of the periodic model. A complete mathematical justification of the homogenization process is the object of a separate publication [11]. For the inverse problem, studying the uniqueness and robustness of determining the coefficients of the ODE model from noisy data is another future direction of research.