Nicolas Moutal
email: nicolas.moutal@polytechnique.edu

Ivan I Maximov
email: ivan.maximov@psykologi.uio.no

Denis S Grebenkov
email: denis.grebenkov@polytechnique.edu

Surface-to-volume ratio of an anisotropic medium by diffusion NMR with general gradient encoding

Keywords: Time-dependent diffusion coefficient, NMR, MRI, Mitra's formula, Surface-to-volume ratio, Spherical encoding, Anisotropy

Since the seminal paper by Mitra et al., diffusion MR has been widely used in order to estimate surfaceto-volume ratios. In the present work we generalize Mitra's formula for arbitrary diffusion encoding waveforms, including recently developed q-space trajectory encoding sequences. We show that surface-to-volume ratio can be significantly misestimated using the original Mitra's formula without taking into account the applied gradient profile. In order to obtain more accurate estimation in anisotropic samples we propose an efficient and robust optimization algorithm to design diffusion gradient waveforms with prescribed features.

I. INTRODUCTION

I N a seminal paper, Mitra et al. have derived the short- time asymptotic behavior of the time-dependent diffusion coefficient in restricted geometries [1]:

D(t) = D 1 - 1 d 4 3 √ π S V √ Dt + O(t) , (1)
where D is the intrinsic diffusion coefficient, d is the space dimensionality, S/V is the surface-to-volume ratio of the medium, and O(t) means that the next term is at most of order of t. Mitra's formula describes the decrease of the diffusion coefficient due to restriction of spin-carrying molecules by the boundaries of the medium at short diffusion times. Higherorder terms of Mitra's formula expansion were analyzed as well and provided additional information about the medium structure such as mean curvature, permeability and surface relaxation [2]- [8].

Using pulsed-gradient spin-echo (PGSE) experiments [9] the diffusion coefficient D(t) can be estimated from the diffusion signal attenuation where the time t is substituted by ∆, a time between two very short gradient pulses. This sort of protocol was often used to estimate the surface-to-volume ratio of porous media [5], [START_REF] Fordham | Proceedings of the Second International Meeting on Recent Advances in MR Applications to Porous Media[END_REF]- [14]. However, such sequences typically require high gradients and do not take advantage of the experimental variety of gradient encoding schemes.

Mitra's formula (1) was extended to constant field gradient [15] which received experimental validation in [5]. An extension to an arbitrary gradient waveform was later derived in [8]. The particular case of oscillating gradients was considered in [16]. It was recently experimentally demonstrated that such sequences make the estimation of S/V accessible to smallgradient hardware, such as clinical scanners [17].

In the article by Mitra et al., the factor 1/d was claimed to be valid for any medium of dimensionality d, by extrapolating results obtained with a sphere (d = 3), a cylinder (d = 2), and a slab (d = 1). It was pointed out in the review [8] that an anisotropic medium should yield different S/V ratios depending on the gradient orientation with respect to the medium. As the structure of the medium is probed by diffusion, the diffusion length (typically of the order of microns for water) naturally distinguishes three types of anisotropy:

• The microscopic anisotropy on much smaller scales than the diffusion length (e.g., intracellular structure with submicron-sized organelles); • The mesoscopic anisotropy on scales comparable to the the diffusion length (this is typically the size of pores, cells, or other confining domains); • The macroscopic anisotropy on much larger scales than the diffusion length that correspond to the size of the voxel. The microscopic anisotropy is usually modeled via a nonisotropic diffusion tensor D [18]- [21]. Mesoscopic anisotropy, on the other hand, manifests itself in the shape of individual compartments or pores whereas macroscopic anisotropy is related to orientation dispersion. For instance, diffusion tensor imaging typically characterizes macroscopic anisotropy via fractional anisotropy (FA) and microscopic anisotropy via micro-fractional anisotropy (µFA) [22]- [24].

Since short-time experiments deal with small diffusion length scales (a few microns for liquids), anisotropy tends to be relevant at the mesoscopic and macroscopic scales rather than at the microscopic one. For this reason, throughout this article we focus on mesoscopic and macroscopic anisotropy of the confining medium. We extend previously obtained results to arbitrary linear gradient encoding schemes and obtain a generalization of Mitra's formula to gradient profiles that can change their amplitude in all directions. This is particularly important for the analysis of diffusion signal acquired by using q-space trajectory encoding schemes [24], including, e.g., isotropic diffusion weighting [25]- [START_REF] Topgaard | Proceedings of the 11th Internatioanl Bologna Conference on Magnetic Resonance in Porous Media (MRPM11)[END_REF]. In this derivation, we employ the methods discussed in reviews [8], [START_REF] Axelrod | [END_REF], [32].

The paper is organized as follows: in Sec. II, we introduce some notations and present our generalization of Mitra's formula. The technical computations are detailed in Appendix I. The proposed formula differs from the classical one (1) by a numerical factor η which is shown to be dependent on the structure of the medium and on the applied gradient waveform. In Sec. III, we study the effect of structure, in particular, of the anisotropy of the confining domain. We first consider a single domain and then evaluate the influence of orientation dispersion on the scale of a voxel. Section IV is devoted to a design of gradient waveforms and their influence on the estimated parameters. We start with the simpler case of linear encoding, for which we recover and extend earlier results. In particular, we show that the diffusion encoding waveform significantly influences the factor η, that may lead to important errors on the estimated S/V ratio. The minimal and maximal achievable values of η are explained in Appendix II. After that, we turn to 3D gradient encoding schemes, with a focus on spherical encoding techniques. We show that typical spherical encoding sequences do not perfectly average out the anisotropy of the medium in the generalized Mitra's formula. In the end, we present a simple algorithm that allows one to design optimal gradient sequences with prescribed properties.

II. RESULTS

We assume that spin-carrying molecules diffuse with diffusivity D in a restricted domain Ω, in the presence of a magnetic field gradient g(t), with t ∈ [0, T ]. Here, t = 0 corresponds to the beginning of the gradient sequence after the 90 • radiofrequency (rf) pulse and t = T corresponds to the echo time at which the signal is acquired (see Fig. 1). We presume that there are no magnetic impurities near the domain boundaries, so that the gradient is uniform in the domain. We define

q(t) = γ t 0 g(t ′ ) dt ′ , ( 2
)
where γ is the gyromagnetic ratio of the spin-carrying molecules, and

b = T 0 |q(t)| 2 dt (3)
is the conventional b-value. The gradient profile is supposed to obey the refocusing condition

q(T ) = T 0 g(t ′ ) dt ′ = 0. (4)
From an experimental point of view, this means that g(t) is the "effective" gradient which takes into account the effect of rf-pulses on the spins (for example, the gradient is effectively reversed by a 180 • rf pulse) [33]. This convention allows us to treat spin echo, gradient echo, stimulated echo, and other techniques, with the same formalism. At small b-values (that is, bD ≪ 1), the MR signal attenuation E can be written as

E ≈ exp(-bD(T )) , (5)
where D(T ) is the effective (or apparent) diffusion coefficient probed by diffusion MR. We generalize the Mitra's formula (1) as

D(T ) = D 1 -η 4 3 √ π S V √ DT + O(T ) (6)
by introducing the numerical prefactor η that depends both on the structure of the medium and on the gradient waveform. We stress that dependence on the waveform implies that one cannot, strictly speaking, interpret D(T ) as a measure of mean-squared displacement of randomly diffusing molecules, except for the theoretical case of two infinitely narrow gradient pulses.

Using mathematical methods discussed in [8], [START_REF] Axelrod | [END_REF], [32], we derived in Appendix I that

η = Tr(WF (3) ) b √ T , ( 7
)
where Tr is the trace. Here we introduced the "temporal" matrix F (3) which is a particular case of the general F (m) matrices,

F (m) = - γ 2 2 T 0 T 0 g(t 1 )⊗g(t 2 ) |t 2 -t 1 | m/2 dt 1 dt 2 , (8)
and the "structural" matrix

W = 1 S ∂Ω n ⊗ n dS , (9)
where the integration is performed on the boundary ∂Ω of the domain Ω and n is the unit outward normal vector to the boundary. We used the tensor product notation ⊗: if a and b are vectors, then a ⊗ b is a matrix with components

(a ⊗ b) ij = a i b j . (10)
With these notations, F (2) is actually the conventional bmatrix [18]- [20] and the b-value can be simply expressed as (see Appendix I from Eq. (58) to Eq. (60) for a detailed computation) b = Tr(F (2) ) .

Thus, the correction factor η in Eq. ( 6) is the combination of the matrices W and F (3) which characterize the medium structure and the gradient sequence, respectively. Note that if one multiplies the gradient waveform g(t) by a factor A, then b and F (3) are both multiplied by A 2 and η is unchanged. In the same way, if one applies a dilatation of the time interval [0, T ] with a factor B, then b is multiplied by B 3 , F (3) is multiplied by B 7/2 , √ T is multiplied by B 1/2 but η remains unchanged again. Finally, η is invariant under dilatation of the domain Ω due to the normalization factor 1/S in the definition of W. The higher-order terms in the asymptotic expansion (6) involve increasing powers of √ DT associated with the temporal matrices F (m) with increasing integer m. These matrices are coupled with structural matrices (such as in Eq. ( 7)) that characterize the medium structure and properties such as curvature, permeability or surface relaxation. However, these properties do not affect the first-order term in (6), on which we focus in this paper.

Mitra's formula (1) was derived for PGSE experiments with (infinitely) short gradient pulses, where t = ∆ is the interpulse time. We emphasize that for general gradient profiles, ∆ is not defined anymore, and we use instead echo time T in our generalized formula (6). If we compare the two formulas by setting t = T (which corresponds to the profile shown on the first panel in Fig. 1), we see that Mitra's formula corresponds to the simple expression

η Mitra = 1/d . ( 12
)
Below we show that this relation is incorrect for general medium structures and gradient profiles.

III. DEPENDENCE ON THE STRUCTURE

A. Simple shapes

For any bounded domain Ω, the matrix W is symmetric, positive-definite, and one has Tr(W) = 1. For example, if Ω is a sphere, one can obtain W sphere = I/3, which is invariant under any spatial rotations of the medium (i.e., proportional to the 3 × 3 unit matrix I), as expected. Throughout the article, we call such matrices "isotropic". However, the same result holds if Ω is a cube, i.e. the cube is also qualified as isotropic by the W matrix. Clearly, the matrix W does not uniquely characterize the shape of Ω.

Let us now consider the example of a rectangular parallelepiped. We choose its sides along the axes (e x , e y , e z ) and denote their length by a, b, c. Then the normal vector n is either ±e x , ±e y , or ±e z depending on the facet of the parallelepiped, and by integrating over each facet we get

W = 1 bc + ca + ab   bc 0 0 0 ca 0 0 0 ab   . ( 13
)
One can see from this simple example that, by varying a, b, c, the matrix W can be any symmetric positive-definite matrix with unit trace. We emphasize that a diagonal form of the matrix is the result of a convenient choice of the parallelepiped's axes.

In the limit when one side of the parallelepiped (say, c) tends to infinity (or is much bigger than the other two), the rectangular parallelepiped transforms into a cylindrical domain with a rectangular cross-section and the W matrix becomes

W = 1 a + b   b 0 0 0 a 0 0 0 0   . (14)
Note that in the special case a = b (cylindrical domain with square cross-section), one obtains the same result as for a circular cylinder of axis e z : W cyl = (Ie z ⊗ e z )/2. In the opposite limit where a and b are much bigger than c, the parallelepiped transforms into a slab perpendicular to e z and the W matrix becomes W slab = e z ⊗ e z . One recognizes in the previous examples the factor 1/d of Mitra's formula (1): 1/3 for a sphere, 1/2 for a circular cylinder, and 1 for a slab. However, even in these basic cases, the factor η remains affected by the gradient waveform, as discussed in Sec. IV.

B. The effect of orientation dispersion

Now we consider a random medium consisting of infinite circular cylinders with random orientation and random radii. For the sake of simplicity, we assume the radius of each cylinder to be independent from its orientation. Equations ( 5) and ( 6) describe the signal on the mesoscopic scale (one cylinder). In order to estimate the macroscopic signal, an average over the voxel is needed. Within the scope of small b-values (bD ≪ 1), one has for the macroscopic signal:

E = exp(-bD(T )) ≈ exp(-b D(T ) ) , (15)
where the average is performed over all cylinders inside the voxel. Coming back to Eqs. ( 6) and ( 7), we see that the average of D(T ) is obtained through the average of the W matrices of the cylinders, that we now compute.

From the previous section, the W matrix of a cylinder oriented along any direction u (where u is a unit vector) is

W cyl (u) = 1 2 (I -u ⊗ u) . (16)
Moreover, for a cylinder of radius R, one has S/V = 2/R, thus the voxel-averaged effective diffusion coefficient reads

D(T ) = D 1 - 4 3 √ π 2 R Tr W cyl F (3) b √ T + O(T ) .
The averaged matrix W cyl depends on the angular distribution of the cylinder orientations. For example, a distribution with a rotation symmetry around the z-axis yields

W cyl = 1 6   2 + OP 0 0 0 2 + OP 0 0 0 2 -2OP   , ( 17
)
where OP is the orientation order parameter of the medium that is defined as

OP = 3 cos 2 θ -1 /2 = 3u 2 z -1 /2 , ( 18
)
where θ is the random angle between the axis of the cylinder u and the symmetry axis e z . The parameter OP can take any value from -1/2 (for θ = π/2, i.e. all the cylinders are in the x -y plane) to 1 (for θ = 0, i.e. all the cylinders are aligned with e z ). The special value OP = 0 corresponds to an isotropic matrix W = I/3 and can be obtained, for example, with a uniform distribution [22]- [24].

The orientation order parameter has direct analogies with other diffusion models describing the water diffusion in strongly anisotropic medium such as the human brain. For instance, if one considers randomly oriented neurites with a Watson distribution of parameter κ [34], then one can compute

OP = 3 2 √ πκ e -κ erfi( √ κ) - 3 4κ - 1 2 , ( 19
)
where erfi is the imaginary error function. In the limit κ going to -∞, we obtain OP = -1/2. If κ = 0, the distribution is uniform and OP = 0. If κ tends to +∞, then OP = 1.

An important consequence of the above computations is that experiments at short diffusion times and small-amplitude gradients are unable to distinguish the mesoscopic anisotropy (the anisotropy of each cylinder) inside a macroscopically isotropic medium (uniform distribution of the cylinders). Therefore, regimes with longer diffusion times or higher gradients are needed in order to have an option for extracting mesoscopic diffusion information [22], [23].

IV. DEPENDENCE ON THE GRADIENT WAVEFORM

In this section we investigate the dependence of the correction factor η (and of higher-order terms) on the gradient waveform captured via the F (m) matrices. We begin with the simpler case, the so-called linear gradient encoding, where the gradient g(t) has a fixed direction and each F (m) matrix is reduced to a single scalar. We show that significant deviations from the classical formula (1) arise depending on the chosen waveform.

Next, in Sec. IV-B, we study how the correction factor is affected in the most general case when both gradient amplitude and direction are time dependent. In particular, we show that recently invented spherical encoding sequences [29], [START_REF] Topgaard | Proceedings of the 11th Internatioanl Bologna Conference on Magnetic Resonance in Porous Media (MRPM11)[END_REF] do not provide the full mixing effect in the sense that η still depends on the orientation of the (anisotropic) medium. In order to resolve this problem we describe in Sec. IV-C a simple and robust algorithm to design diffusion gradient profiles with desired features and constraints.

A. Constant gradient direction

If we set g(t) = g(t)e, with a constant unit vector e, the F (m) matrices become

F (m) = F (m) e ⊗ e , (20)
with the scalar

F (m) = T 0 T 0 g(t 1 )g(t 2 ) |t 2 -t 1 | m/2 dt 1 dt 2 . ( 21
)
The numerical factor η becomes

η = F (3) b √ T e T We . (22)
Thus, by keeping the same profile g(t) and only changing the direction of the applied gradient e, the factor F (3) /(b √ T ) is unchanged and the factor (e T We) allows one to probe the whole W matrix. For this purpose, one can transpose standard diffusion tensor imaging reconstruction techniques [18] to our case: by choosing multiple non-coplanar directions for e, one obtains a system of linear equations on the coefficients of W that can be solved numerically. Bearing in mind that W is symmetric positive-definite matrix with trace one, one needs at least 6 diffusion directions to estimate 5 independent coefficients of the W matrix and the S/V ratio.

For a W matrix such as the one of a parallelepiped (13), the factor η takes different values depending on the gradient direction e. For instance, if one applies the gradient in a direction perpendicular to the smallest facets of parallelepiped, one probes the S/V ratio of these facets, not of the whole structure. Although this example is specific, the conclusion is general: the mesoscopic anisotropy of a confining domain, captured via the matrix W, can significantly bias the estimation of the surface-to-volume ratio. This circumstance was ignored in some former studies with application of the classical Mitra's formula, which is only valid for isotropic domains. While spherical encoding scheme aims to resolve this issue by mixing contributions from different directions, we will see in Sec. IV-B that this mixing is not perfect for formerly proposed spherical encoding sequences.

In the remaining part of this subsection, we consider the particular case of isotropic (e.g., spherical) domains with W = I/3 so that the structural aspect is fully decoupled from the temporal one. In this case, Eq. ( 7) yields

η = α 3 , with α = F (3) b √ T , (23)
and we can focus on the temporal aspect (gradient waveform) captured via the factor α.

Figure 1 shows some examples of temporal profiles and the corresponding values of α. The maximum attainable value of α is slightly over 1 (around 1.006), the minimum value being 0 (see Appendix II for more details). One can achieve very small values of α by using very fast oscillating gradients. In fact, for sinusoidal gradient waveforms of angular frequency ω one has α ∼ ω -1/2 , in the limit ωT ≫ 1 [16], [17].

This finding has an important practical consequence: if one ignores the factor α and uses the original Mitra's formula (for which α = 1), one can significantly underestimate the surfaceto-volume ratio (by a factor 1/α) and, thus, overestimate the typical size of compartments.

B. Isotropy and spherical encoding

Microscopic anisotropy is usually modeled via a nonisotropic diffusion tensor D, and the expression (5) for the diffusion signal becomes

E = exp -Tr F (2) D . ( 24
)
Typical spherical encoding sequences [25]- [START_REF] Topgaard | Proceedings of the 11th Internatioanl Bologna Conference on Magnetic Resonance in Porous Media (MRPM11)[END_REF] aim to average out the microscopic anisotropy of the medium by applying an encoding gradient with time-changing direction. Mathematically, the goal is to obtain an "isotropic" F (2) matrix, F (2) ∝ I, so that the signal in Eq. ( 24) depends only on the trace Tr(D) and thus yields the same result for any orientation of the medium. We recall that throughout the paper, we call a matrix isotropic if it is proportional to the unit matrix I (in other words, its eigenvalues are equal to each other). Mesoscopic anisotropy manifests itself in the W matrix of individual compartments, as we explained in Sec. III. Thus, from Eq. ( 6) we can deduce that mesoscopic anisotropy is averaged out by the gradient sequence only if F (3) is isotropic. In this case, the factor η does not depend on the orientation of the mesoscopically anisotropic medium, and one can estimate precisely the surface-to-volume ratio of the medium. Moreover, from Eq. (7) we see that in this case, η can be read directly from the expression of F (3) :

F (3) iso /(b √ T ) = η I . (25)
Similarly, it is expected that the isotropy condition for the matrices F (4) , F (5) , . . . would be needed if the higher-order terms were considered.

Hence, the natural question arises: "Do the former spherical encoding sequences that were designed to get an isotropic F (2) (or b) matrix [24], produce isotropic F (m) matrices (or at least F (3) )?". For instance, for the q-Space Magic-Angle-Spinning sequence [29], [START_REF] Topgaard | Proceedings of the 11th Internatioanl Bologna Conference on Magnetic Resonance in Porous Media (MRPM11)[END_REF] we obtained

F (3) /(b √ T ) =   0.83 0 0.29 0 0.42 0 0.29 0 0.51   . (26)
This matrix has eigenvalues [0.34, 0.42, 1.00] and is thus not isotropic. Similarly, for the FAMEDcos sequence [36], we get

F (3) / √ T =   0.39 0 0.04 0 0.32 0 0.04 0 0.29   , (27)
which is also not isotropic. All spherical encoding schemes that we could find in the literature produce anisotropic F (3) matrices.

C. How to obtain isotropic matrices?

The question in the subsection title can be restated in an algebraic language: how to find three functions g x (t), g y (t), g z (t) with zero mean (see Eq. ( 4)) that are mutually "orthogonal" and have the same "norm" for a given set of symmetric bilinear forms φ m , m = 2, 3, . . ., with

φ m (f 1 , f 2 ) = T 0 T 0 f 1 (t 1 )f 2 (t 2 )|t 1 -t 2 | m/2 dt 1 dt 2 .
(28) Since the space of functions with zero mean is infinitedimensional, we can be confident in finding such three functions. However, Eq. ( 28) involves a non-integer power of time that prevents us from getting analytical solution for this problem. Thus, we design a simple algorithm for generating the gradient sequences that satisfy these conditions.

The idea is to choose a family of functions (f 1 , f 2 , . . . , f k ) (for example, sines or polynomials, possibly with phase jumps at T /2) and to search for g x (t), g y (t), g z (t) as linear combinations of the basis functions. This is a generalization of the classical sine and cosine decomposition which was already used in the context of waveform optimization [START_REF] Topgaard | Proceedings of the 11th Internatioanl Bologna Conference on Magnetic Resonance in Porous Media (MRPM11)[END_REF]. Mathematically, this means that

  g x (t) g y (t) g z (t)   = X      f 1 (t) f 2 (t)
. . .

f k (t)      , ( 29
)
where X is a 3 × k matrix of coefficients to be found. Now we define the k × k matrices A (m) by

A (m) i,j = φ m (f i , f j ) , m = 2, 3, . . . (30)
In this way, one can compute directly the F (m) matrices according to

F (m) = XA (m) X T , (31)
where X T is the transpose of X. The problem is then reduced to an optimization problem for the matrix X, which can be easily done numerically. In other words, one searches for a matrix X that ensures simultaneously the diagonal form of the matrices F (2) , F (3) , etc. The optimization algorithm can include various additional constraints. On one hand, one has a freedom to choose an appropriate family (f 1 , f 2 , . . . , f k ), for example, to ensure smoothness of the resulting gradient profile. Similarly, the refocusing condition can be achieved by choosing zero-mean functions. On the other hand, it is also possible to add some constraints as a part of the optimization problem. This is especially easy if the constraints can be expressed as linear or bilinear forms of the gradient profile g(t). For instance, each F (m) matrix in ( 8) is a bilinear form of the gradient profile allowing one to express them as the simple matrix multiplication [START_REF] Axelrod | [END_REF]. Another example of additional conditions consists in imposing zeros to the designed gradient profiles. For practical reasons, it is often easier to manipulate with the gradients that satisfy to

g(0) = g(T /2) = g(T ) = 0 . ( 32
)
This is a linear condition on the gradient profile. If one denotes by V the k × 3 matrix

V =      f 1 (0) f 1 (T /2) f 1 (T ) f 2 (0) f 2 (T /2) f 2 (T ) . . . . . . . . . f k (0) f k (T /2) f k (T )      , (33)
then Eq. ( 32) becomes

XV =   0 0 0 0 0 0 0 0 0   . (34)
In the following, we impose the above condition to produce our gradient waveforms.

It is worth to note that one can also generate flowcompensated gradients, or more generally, motion artifacts suppression techniques. In order to design motion insensitive gradient waveforms one needs to impose linear conditions on the gradient profile

T 0 t p g(t) dt = 0 , p = 1, 2, . . . , N , (35)
where p = 1 corresponds to the flow compensation, and higher values of p account for acceleration, pulsatility, etc. [37], [38]. This condition can be rewritten in the matrix form XM = 0, where the k × N matrix M is defined by

M i,p = T 0 t p f i (t) dt . ( 36
)
Another type of optimizaton constraints can be based on hardware limitations such as a need to minimize heat generation during the sequence execution which amounts to minimizing the following quantity

g, g = T 0 |g(t)| 2 dt , (37)
which is a bilinear form of the gradient profile. Similar to representation [START_REF] Axelrod | [END_REF] for F (3) , one can define a matrix H i,j = f i , f j to write Eq. ( 37) as g, g = Tr XHX T , and then to include it into the optimization procedure.

The previous examples showed how linear and bilinear forms of the gradient profile can be simply expressed in terms of the weights matrix X, which allows one to perform very fast computations. The matrix corresponding to each condition (for example, A (n) , V, H) has to be computed only once, then optimization is reduced to matrix multiplications. The size of the matrices involved in the computations is defined by the size of the chosen set of functions (f 1 , . . . , f k ). Note that the set size is independent of the numerical sampling of the time interval [0, T ] that controls accuracy of the computations. Some properties of the designed gradients do not fall into the category of aforementioned linear or bilinear forms, e.g., "max" amplitude-function (i.e., one cannot impose the maximal gradient constraint in this way). They can naturally be included in the optimization, however one cannot apply the previous techniques in order to speed up their computation.

We have to emphasize that the "optimal" solution is not unique and it depends on the choice of the set f 1 , . . . , f k . Moreover, if the set is sufficiently large and the number of degrees of freedom is greater than the number of constraints, then the algorithm will likely yield different solutions depending on the initial condition for X. This property can be advantageous in practice, as one can design many optimal solutions. The described optimization algorithm was implemented in Matlab (The MathWorks, Natick, MA USA) and is available upon request. It concatenates all the chosen constraints in a single vector-valued function f (X) of the weight matrix X, in such a way that the constraints are expressed by the condition f (X) = 0. This equation is then solved numerically with the Levenberg-Marquardt algorithm.

Figure 2 shows two examples of gradient waveforms that produce isotropic F (2) and F (3) matrices. These profiles were obtained from a set with k = 9 sine functions and k = 9 polynomials, respectively. The first set was composed of cos(πjt/T ) with j = 1, . . . , 5; sin(πjt/T ) with j = 2, 4, 6; and ε(t) sin(4πt/T ) where ǫ(t) is a piecewise constant function that is equal to 1 on [0, T /2] and -1 on (T /2, T ]. The second set was composed of a mixture of monomials, symmetrized odd monomials and antisymmetrized even monomials, with zero mean:

(t -T /2), (t -T /2) 2 - T 2 /12, (t -T /2)|t -T /2|, (t -T /2) 3 , |t -T /2| 3 -T 3 /32, (t -T /2) 4 -T 4 /80, (t -T /2) 3 |t -T /2|, (t -T /2) 5 , |t - T /2| 5 -T 5 /192.
Although the combination of symmetric and antisymmetric functions helped us to increase the number of basis functions while keeping low degree monomials or slowly oscillating sines, one could alternatively use just monomials, polynommials, or other basis functions as well.

The condition of isotropy for both matrices F (2) and F (3) yields 5+5 equations on the components of matrix X. Besides of matrices F (2,3) , Eq. ( 32) adds another 9 equations on the components of X. Moreover, we imposed the b-value so that the algorithm satisfied 20 conditions with 3k = 27 degrees of freedom.

In both cases, the gradient waveform corresponds to η ≈ 0.1. We observed that the algorithm could not produce gradient waveforms with arbitrary values of η. There are bounds for η-values outside of which the optimization process did not converge. This behavior was expected, because even in the one-dimensional case, one already had some mathematical limitations for the parameter η (see Sec. IV-A and Appendix II). These bounds can be extended by by adding more basis functions (i.e., by increasing the size k of their set). Another way to extend the bounds is to reduce the number of constraints, for example, by dropping out the condition of isotropic F (2) matrix and only keeping the condition on F (3) .

Interestingly, the F (4) matrix represents a special case: integrating by parts in Eq. ( 8) one can show that

F (4) = T 0 q(t) dt ⊗ T 0 q(t) dt . (38)
This implies that the matrix has rank one, so it cannot be proportional to the unit matrix unless it is null, that occurs under the simple condition

T 0 q(t) dt = 0 . ( 39
) 0 T/2 T 0 T/2 T
x y z y x z Fig. 2. Two examples of gradient waveforms that produce isotropic F (2) and F (3) matrices and that satisfy Eq. (32). Note that the gradients are "effective" gradients in the sense that we reversed them after the 180 • rf pulse at T /2. The bottom figure shows the corresponding q(t). (left) the profiles are combination of 9 piecewise sine functions with frequencies up to 6/T ; (right) the profiles are piecewise polynomials of order 5, and in addition they satisfy F (4) = 0. This condition can be easily included in our optimization algorithm. This is the case for the designed profile shown on the right panel in Fig. 2. As a consequence, if there is no effects of permeation and surface relaxivity, the corresponding term (of the order of DT ) vanishes in the expansion (6).

This property is well-known for cosine-based waveforms with an integer number of periods [39], and, indeed, such functions automatically satisfy to Eq. (39). However, this property is not exclusive to sine functions (for example, the right panel of Fig. 2 was obtained with polynomial functions). It is also easy to show that Eq. ( 39) is equivalent to condition (35) for n = 1. In other words, flow-compensated gradient profiles automatically cancel the DT -order correction term in the generalized Mitra's expansion, as it was already pointed out earlier in [38].

V. CONCLUSION

We presented a generalization of Mitra's formula that is applicable to any gradient waveform and any geometrical structure. This generalized formula differs from the classical one by a numerical factor in front of S/V . In the case of linear encoding schemes, we showed that this factor can significantly affect the estimation of S/V and lead to overestimated size of domain compartments.

We also discussed in detail the effect of anisotropy of the medium and the use of spherical encoding schemes. In particular, we showed that in order to estimate the surface-tovolume ratio of an anisotropic medium, the gradient should satisfy the isotropy condition (F (3) ∝ I) that is different from the usual one (F (2) ∝ I). As a consequence, typical spherical encoding schemes do not satisfy this new condition. We presented a simple and flexible algorithm that allows fast optimization of gradient waveforms and is well-suited for design of diffusion weighted sequences with specific features such as isotropy of F (3) , flow compensation, heat limitation, and others.

The developed extension of Mitra's formula has a significant practical impact, in particular, in medical applications where fast non-invasive imaging techniques are crucial. The proposed approach characterizes the underlying microstructure via novel quantitative metrics such as W-tensor, orientation order parameter and more accurate surface-to-volume ratio. The quantitative scalar maps based on those metrics possess a high potential as a novel set of biomarkers and allow one to apply for both well-known diffusion tensor formalism and further improvement of diffusion models based on compartmentization such as NODDI [34]. The practical advantages of the developed approach for designing new gradient encoding schemes for in vivo brain imaging on clinical scanners will be demonstrated in a separate publication.

APPENDIX I THEORETICAL COMPUTATIONS

The signal is proportional to the expectation of the transverse magnetization which has a form of the characteristic function of the random dephasing φ acquired by diffusing spin-carrying molecules:

E = E{e iφ }, φ = γ T 0 B(r t , t) dt , ( 40
)
where T is the echo time, r t is the random trajectory of the nucleus, γ is the gyromagnetic ratio, and γB(r, t) is the Larmor frequency corresponding to the magnetic field. In this work, we consider the most general form of the linear gradient g(t):

B(r, t) = (g(t) • r) = g x (t)x + g y (t)y + g z (t)z. (41)
In particular, the dephasing can be decomposed as

φ = φ x + φ y + φ z , φ i = γ T 0 dt g i (t)(e i • r t ) (i = x, y, z), (42)
where e x , e y and e z are the units vectors in three directions, and (e i • r t ) is the projection of the molecule position at time t onto the direction e i .

The effective diffusion coefficient is related to the second moment of the dephasing, i.e., we need to evaluate

E{φ 2 } = i,j=x,y,z E{φ i φ j }. (43)
We emphasize that the three components φ x , φ y and φ z are independent only for free diffusion, whereas confinement would typically make them correlated. In other words, one cannot a priori ignore the cross terms such as E{φ x φ y }.

In order to compute these terms, we use the following representation [8]:

E{φ i φ j } = γ 2 T 0 dt 1 T t1 dt 2 Ω dr 0 Ω dr 1 Ω dr 2 Ω dr 3 × ρ(r 0 )P t1 (r 0 , r 1 )P t2-t1 (r 1 , r 2 )P T -t2 (r 2 , r 3 ) × [B i (r 1 , t 1 )B j (r 2 , t 2 ) + B j (r 1 , t 1 )B i (r 2 , t 2 )] (44)
where P t (r, r ′ ) is the propagator in the domain Ω, and ρ(r 0 ) is the initial density of particles (the initial magnetization after the 90 • rf pulse). If the boundary is fully reflecting and ρ(r 0 ) is uniform, then the integrals over r 0 and r 3 yield 1, so that

E{φ i φ j } = γ 2 V T 0 dt 1 T t1 dt 2 Ω dr 1 Ω dr 2 P t2-t1 (r 1 , r 2 ) × [B i (r 1 , t 1 )B j (r 2 , t 2 ) + B j (r 1 , t 1 )B i (r 2 , t 2 )] , ( 45
)
where V is the volume of the domain. We get thus

E{φ i φ j } = γ 2 T 0 dt 1 g i (t 1 ) T t1 dt 2 g j (t 2 )K ij (t 2 -t 1 ) + T 0 dt 1 g j (t 1 ) T t1 dt 2 g i (t 2 )K ji (t 2 -t 1 ), (46)
where

K ij (t) = Ω Ω k i (r 1 )P t (r 1 , r 2 )k j (r 2 ) dr 1 dr 2 , (47)
with k i (r) = (e i • r). Since K ij (t) = K ji (t) due to the symmetry of the propagator, we can rewrite the moment as

E{φ i φ j } = γ 2 T 0 g i (t 1 ) T 0 g j (t 2 )K ij (|t 2 -t 1 |) dt 1 dt 2 .
(48) We rely on the general short-time expansion for the heat kernels (see [START_REF] Davies | Heat Kernels and Spectral Theory[END_REF]- [START_REF] Desjardins | [END_REF] and references therein)

K ij (t) = m≥0 c m (k i , k j ) (Dt) m/2 , ( 49
) with c 0 (f, g) = 1 V Ω f (r)g(r) dr , (50a)
c 1 (f, g) = 0 , (50b)
c 2 (f, g) = - 1 V Ω (∇f (r) • ∇g(r)) dr , (50c)
c 3 (f, g) = 4 3 √ π 1 V ∂Ω ∂f (r) ∂n ∂g(r) ∂n dS , (50d)
where ∂/∂n = (n•∇) is the normal derivative at the boundary, and n is the unit normal vector at the boundary oriented outward the domain. We note that the expansion (49) is an asymptotic series which has to be truncated. In our case, we get

c 0 (k i , k j ) = 1 V Ω (e i • r)(e j • r) dr , (51a)
c 1 (k i , k j ) = 0 , (51b) c 2 (k i , k j ) = -δ ij , (51c)
c 3 (k i , k j ) = 4 3 √ π 1 V ∂Ω (e i • n)(e j • n) dS (51d)
(in the last integral, the normal vector n depends on the boundary point). Combining these results, we get

E{φ i φ j } = γ 2 T 0 dt 1 g i (t 1 ) T 0 dt 2 g j (t 2 ) × -δ ij (D|t 2 -t 1 |) + 4 3 √ π S V W ij (D|t 2 -t 1 |) 3/2 + • • • , ( 52
)
where S is the surface area, the matrix W is defined by

W = 1 S ∂Ω n ⊗ n dS , (53)
and the zeroth order term (with c 0 ) vanished due to the rephasing condition

T 0 g i (t) dt = 0 (i = x, y, z). (54)
We can write this result more compactly as

E{φ i φ j /2} = Dδ ij F (2) ij - 4D 3/2 3 √ π S V W ij F (3) ij + • • • , (55)
where we introduced the matrices

F (m) = - γ 2 2 T 0 T 0 g(t 1 ) ⊗ g(t 2 ) |t 2 -t 1 | m/2 dt 1 dt 2 . ( 56
)
As a consequence, we compute the second moment as

E{φ 2 /2} = D Tr(F (2) ) - 4D 3/2 3 √ π S V Tr(WF (3) ) + • • • . (57)
Let us introduce the auxiliary function

h(t 1 ) = T 0 g(t 2 )|t 2 -t 1 | dt 2 . (58)
We split the above integral and perform an integration by parts

h(t 1 ) = 1 γ t1 0 q(t 2 ) dt 2 - 1 γ T t1 q(t 2 ) dt 2 ,
where we used the conditions q(0) = 0 and q(T ) = 0. Now we note that

T 0 dt 1 t1 0 g(t 1 ) ⊗ q(t 2 ) dt 2 = T 0 dt 2 T t2 g(t 1 ) ⊗ q(t 2 ) dt 1 = - 1 γ T 0 q(t 2 ) ⊗ q(t 2 ) dt 2 ,
where we used again q(T ) = 0. In the same way one gets

T 0 dt 1 T t1 g(t 1 ) ⊗ q(t 2 ) dt 2 = 1 γ T 0 q(t 2 ) ⊗ q(t 2 ) dt 2 .
Putting all the pieces together, one finally obtains

F (2) = T 0 q(t) ⊗ q(t) dt , (59)
so that F (2) is actually the b-matrix [18]- [20] and one can compute

Tr(F (2) ) = T 0 |q(t)| 2 dt = b , (60)
thus we recover the signal for free diffusion E = e -E{φ 2 /2} = e -bD in the absence of confinement. In turn, the effective diffusion coefficient, which is experimentally determined from the dependence of -ln E on b at small b-value, is expressed through the second moment as

D(T ) = lim b→0 -ln E b = lim b→0 E{φ 2 /2} b (61)
from which we obtain Eq. ( 6).

APPENDIX II OPTIMAL VALUE OF α

In the case of linear gradient encoding in a spherical domain, we obtain that the Mitra formula is corrected by a factor α which is computed from the gradient profile as

α = F (3) b √ T . (62)
In this section, we investigate the maximum and the minimum values of α. Integrating by parts (following the same procedure as in Eqs. ( 58)-( 60)), one obtains

F (3) = 3 8 T 0 T 0 q(t 1 )q(t 2 )|t 1 -t 2 | -1/2 dt 1 dt 2 . ( 63
)
Note that despite its singularity at 0, the function 1/ |t| is integrable, hence the above integral is well-defined. Next, we apply a change of variables from t ∈ [0, T ] to t/T ∈ [0, 1] and q(t) to Q(t/T ), which gives

α = 3 8 Q 2 1 0 1 0 Q(t 1 )Q(t 2 )|t 1 -t 2 | -1/2 dt 1 dt 2 , (64)
with the usual L 2 norm. We can understand the above expression as a scalar product

α = 3 8 Q, AQ Q, Q , (65)
with an integral operator with the kernel |t 1 -t 2 | -1/2

(AQ)(t 1 ) =

1 0 Q(t 2 )|t 1 -t 2 | -1/2 dt 2 . ( 66
)
One can see that A is a convolution operator because the kernel can be expressed as K(t 1 -t 2 ) (with K(t) = 1/ |t|).

Denoting by Q(ω) the Fourier transform of Q(t) and by K(ω) the Fourier transform of K(t), one gets

Q, AQ = 1 2π ∞ -∞ | Q(ω)| 2 K(ω) dω , (67)
with K(ω) = 2π/|ω|. This shows that α is always positive (in other words, the operator A is positive-definite). This result is expected from a physical point of view: if α were negative, then the effective diffusion coefficient would increase with time, what is nonphysical behavior. The minimum value 0 can be asymptotically obtained, for example, with very fast oscillating gradients. It is, indeed, clear from (67) that if g(t) is a cosine function with angular frequency ω 0 such that the number of periods N 0 = ω 0 T /(2π) ≫ 1, then | Q(ω)| 2 is concentrated around ±ω 0 , and we obtain α ≈ 3/(8 √ N 0 ) [16] (see also Fig 1). Now we turn to the maximum value of α. The condition that Q(t) is null outside of [0, 1] is difficult to take into account in Fourier space and we could not extract further information from Eq. (67). In order to bound the maximum value of α, one can use the Cauchy inequality:

|(AQ)(t 1 )| = 1 0 K(t 1 -t 2 )Q(t 2 ) K(t 1 -t 2 ) dt 2 ≤ 1 0 K(t 1 -t 2 ) dt 2 1 2 1 0 Q 2 (t 2 )K(t 1 -t 2 ) dt 2 1 2
.

One can easily compute the function

1 0 K(t 1 -t 2 ) dt 2 = 2 √ t 1 + 2 √ 1 -t 1 , ( 68
)
whose maximum is 2 √ 2. Thus, one gets

|(AQ)(t 1 )| ≤ 2 √ 2 1 0 Q 2 (t 2 )K(t 1 -t 2 ) dt 2 1 2
. (69)

Using again the Cauchy inequality, one obtains

Q, AQ ≤ 2 3/4 Q 1 0 1 0 Q 2 (t 2 )K(t 1 -t 2 ) dt 2 dt 1 1 2
.

The same reasoning about the maximum value of the integral of K yields The problem can be considered from another point of view. Due to the symmetry of the operator A, it is wellknown that the function Q maximizing Q, AQ / Q 2 is the eigenfunction of A with the highest eigenvalue. As a consequence, if one searches for a good estimation of the maximum α as well as the corresponding "optimal" gradient profile, then one can use the following procedure: (i) to choose an initial profile Q 0 which is sufficiently general or sufficiently close to a guessed optimal profile; (ii) to apply iteratively the operator A and to renormalize the result; (iii) to stop when the sequence has converged.

Q, AQ ≤ 2 3/2 Q 2 , (70)
For example, the initial profile Q 0 (t/T ) = 1, which corresponds to two infinitely narrow gradient pulses at time 0 and T , yields α = 1, which is close to the optimal value. Thus, it is a good initial condition for the iterative process. The result of such a procedure is shown in Fig. 3. This yields Fig. 3. The result of the iterative procedure in order to obtain the optimal profile that maximizes the value of α.

an optimum value of α of about 1.006, thus very close to 1. It is worth to note, however, that the optimal profile Q(t/T ) differs clearly from Q 0 (t/T ) = 1.

Fig. 1 .

 1 Fig.1. Illustration of some gradient profiles for spin echo experiments. We stress that these gradient profiles are "effective" in the sense that we reversed the gradients after T /2 to show the effect of the 180 • rf pulse. The corresponding values of α are given for each profile (see Sec. IV-A). Note that α = 1 for the NPA-case (first panel), which corresponds to Mitra's formula (1) where t = ∆ = T .

  from the examples in Fig.1that α = 1 can be achieved, which implies that the maximum value of α is in the interval [1, 1.06].