Analytical performance bounds for multi-tensor diffusion-MRI

    Analytical Performance Bounds for Multi-tensor Diffusion-MRI Farid Ahmed Sid, Karim Abed-Meraim, Rachid Harba, Fatima OulebsirBoumghar PII: DOI: Reference:

S0730-725X(16)30161-8 doi: 10.1016/j.mri.2016.10.014 MRI 8636

To appear in:

Magnetic Resonance Imaging

Received date: Revised date: Accepted date:

15 June 2016 9 September 2016 5 October 2016

Please cite this article as: Ahmed Sid Farid, Abed-Meraim Karim, Harba Rachid, Oulebsir-Boumghar Fatima, Analytical Performance Bounds for Multi-tensor DiffusionMRI, Magnetic Resonance Imaging (2016), doi: 10.1016/j.mri.2016.10.014

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Analytical Performance Bounds for Multi-tensor Diffusion-MRI

RI P

T

Farid AHMED SID 1∗ , Karim ABED-MERAIM 2 , Rachid HARBA 2 and Fatima OULEBSIR-BOUMGHAR 1

1

NU

∗

ParIMéd/LRPE,FEI,USTHB, BP 32 El Alia, Bab Ezzouar, 16111 Algiers, Algeria. PRISME Laboratory, University of Orléans, 12 Rue de Blois, 45067 Orléans, France.

SC

2

Corresponding author at:

AC

CE

PT

ED

MA

ParIMéd/LRPE,FEI,USTHB, BP 32 El Alia, Bab Ezzouar, 16111 Algiers, Algeria. E-mail address : [email protected]

ACCEPTED MANUSCRIPT

ED

MA

NU

SC

RI P

T

Purpose: To examine the effects of MR acquisition parameters on brain white matter fiber orientation estimation and parameter of clinical interest in crossing fiber areas based on the Multi-Tensor Model (MTM). Material and Methods: We compute the Cramér-Rao Bound (CRB) for the MTM and the parameter of clinical interest such as the Fractional Anisotropy (FA) and the dominant fiber orientations, assuming that the diffusion MRI data are recorded by a multi-coil, multi-shell acquisition system. Considering the sum-of-squares method for the reconstructed magnitude image, we introduce an approximate closed-form formula for Fisher Information Matrix that has the simplicity and easy interpretation advantages. In addition, we propose to generalize the FA and the mean diffusivity to the multi-tensor model. Results: We show the application of the CRB to reduce the scan time while preserving a good estimation precision. We provide results showing how the increase of the number of acquisition coils compensates the decrease of the number of diffusion gradient directions. We analyze the impact of the b-value and the Signal-to-Noise Ratio (SNR). The analysis shows that the estimation error variance decreases with a quadratic rate with the SNR, and that the optimum b-values are not unique but depends on the target parameter, the context, and eventually the target cost function. Conclusion: In this study we highlight the importance of choosing the appropriate acquisition parameters especially when dealing with crossing fiber areas. We also provide a methodology for the optimal tuning of these parameters using the CRB. Cramér-Rao Bound; Non-central Chi distribution; Multiple-receiver coils; Multi-Tensor

AC

CE

PT

Model.

ACCEPTED MANUSCRIPT

1 Introduction

AC

CE

PT

ED

MA

NU

SC

RI P

T

It is a real challenge to find in vivo and non invasively how the brain is connected. Unfortunately, with conventional Magnetic Resonance (MR) imaging, the white matter of the human brain has a homogeneous appearance due to the low spatial resolution which is at a scale of millimeters. Diffusion Magnetic Resonance Imaging (dMRI) is the only MRI modality that is able to provide information about the brain’s white matter structure connectivity, usually at a scale of micrometers. This modality is based on the observation of the average displacement of water molecules in every voxel. The water molecules motion alongside the fiber-bundles is fast compared to the perpendicular motion, giving rise to diffusion anisotropy. This property can be exploited to reconstruct the white matter bundles using the Diffusion Tensor (DT) MRI [1] so that the principal direction of the tensor may be related to the fiber orientation afterwards. This model provides useful anisotropy indices such as Fractional Anisotropy (FA) and Mean Diffusivity (MD) that are widely used to diagnose neurodegenerative diseases like Multiple Sclerosis [2], Alzheimer disease [3], or tumors [4]. However, the diffusion tensor model is a relatively good approximation for voxels containing a single fiber-bundle, but is inaccurate in the case of fiber crossing. To overcome this inherent limitation of the DT-MRI model, various approaches based on High Angular Resolution Diffusion Imaging (HARDI) techniques have been proposed. Some of them are model-based and the others are model-free (more details can be found in [5], [6], [7]). However, these approaches have not yet been adopted in routine clinical use due to their long scan times; in addition, most of them are limited to the assessment of connectivity and focus on the characterization of the diffusion profile at each voxel. To deal with crossing fibers, we chose to work with the Multi-Tensor Model (MTM) proposed in [8] which is an extension of the DT model. The key idea behind this model is the fact that voxels containing a single fiber-bundle are, relatively, well represented by the DT model. Then voxels containing more than one fiber-bundles are expected to be well represented by a sum of several tensors. This enables to extract the orientation information and the diffusion properties of each fiber-bundle independently. In contrast to many representation of diffusion signal based on HARDI acquisition, the MTM allows the computation of anisotropy factors such as the FA and MD for each fiber-bundle. So, the MTM seems to be a good candidate for clinical routine use, due to its relatively short acquisition time. But, it has been reported that the related model parameter estimation is quite complex and challenging compared to the DT model [8], [9], [10]. The reason, as shown in [11], is that with a single-shell HARDI acquisition, the model parameters are collinear leading to an infinite number of solutions. Scherrer [11] developed a new optimal acquisition scheme CUSP (CUbe and SPhere) to accurately estimate the full MTM parameters, based on multi-shell acquisitions in a clinically acceptable scan time. To study the precision of the parameter estimation for the multi-tensor model and to determine and analyze the influence of the tunable parameters of MRI scanner as well as the noise level on the estimation accuracy, we use the Cramèr-Rao Bound tool which provides a theoretical lower bound on the variance of any unbiased estimator. In all of our derivations, we considered that the diffusion signal is acquired with an MRI system equipped with multiple coils and using multi-shell (i.e. multiple b-values). Previous works that compute the CRB for diffusion MRI models, particularly [12], [13], [14], [15], provide expressions for the Fisher Information Matrix (FIM) involving integral formulae that require heavy evaluations, and that are prone to numerical errors. The lack of a

ACCEPTED MANUSCRIPT

ED

MA

NU

SC

RI P

T

closed-form expression for the integral formula is the major difficulty that hinders providing analytical analysis. The present paper proposes a simplified approximation for the integral formula and consequently the FIM expression. The presented approximation is not only a more elegant way to solve the problem at hand and speed up the computations but also allows the study of the influence of reconstruction methods and tunable system parameters on the estimation precision. Our contributions can be summarized as follows: • We have first generalized the CRB derivation to the case of multi-tensor, multi-shell and multi-coil acquisition system, using the sum-of-squares (SoS) reconstruction technique. And also, besides the general CRB formula of the MTM model, we paid a special attention to the parameters of clinical interest and developed dedicated CRB expressions for them. • As mentioned above, we introduce an original approximate analytical expression for the FIM which has the advantage of being simple and provides easy computations. • The latter analytical expression is exploited to conduct a thorough analysis of the performance bound. This study leads to many interesting observations related to the optimal selection of the b-value parameter and to the impact of the number of coils and the number of diffusion gradients on the estimation performance. • One of the most important measures for the diagnosis of neurodegenerative diseases is the Fractional Anisotropy (FA), which has been defined only for the DT model [16]. In this article, we propose to generalize this measure to the multi-tensor model before being used for the performance analysis we conducted based on CRB derivation.

2 Material and methods

PT

This section provides a brief review of the noise model in the case of the SoS reconstruction technique. It also contains a description of the diffusion signal model and the parameters of clinical interest. Finally, the Cramér-Rao Bound (CRB) theory is briefly reviewed.

CE

2.1 Noisy Data Model

AC

MRI systems currently provided by the major vendors are equipped with multiple acquisition coils, where the coils are the “antenna" of the MRI system used to transmit signal to the patient and/or to receive the return signal. Each coil, after demodulation and filtering gives two signals treated as the real and imaginary parts of a complex raw data. Afterwards, the two dimensional inverse discrete Fourier transform, of the raw data, results in 𝐿 complex images in the image space, 𝐿 being the number of coils. The complex signal for the 𝑙 𝑡ℎ coil, in the image domain, at a given voxel can be expressed as 𝑆𝑙 = 𝑐𝑙 𝐴 + 𝑛𝑙 , 𝑙 = 1, … , 𝐿

(1)

where 𝑐𝑙 is the complex sensitivity factor of the 𝑙 𝑡ℎ coil at the voxel under consideration, 𝐴 is the image voxel intensity in the absence of noise and 𝑛𝑙 is an additive complex white Gaussian noise process with zero-mean and covariance matrix 𝚺. In diffusion MRI, it is a common practice to work with magnitude images rather than the complex ones, while the phase component is ignored due to the high sensitivity of the phase to motion and eddy currents. If the Fourier space, known as k-space is fully sampled, the composite magnitude signal (CMS) may be obtained by combining the images from all channels using the root-sum-of-squares (SoS) method as described in [17] 2 𝑀 = √∑𝐿𝑙=1 𝑆𝑙∗ 𝑆𝑙 = √∑𝐿𝑙=1 (𝑆𝑅𝑙 + 𝑆𝐼𝑙2 ),

(2)

ACCEPTED MANUSCRIPT

𝒫(𝑀) = 𝜎2 𝐴𝐿−1 𝑒

−

𝑀2 +𝐴2 𝑇 2𝜎2

𝑇

𝑀𝐴𝑇

𝐼𝐿−1 (

𝜎2

1

)=𝜎

𝜂2

𝜂1−2𝐿 𝑒 − 2

𝑥2

− 𝑥 𝐿 𝑒 2𝜂2 𝐼𝐿−1 (𝑥),

RI P

𝑀𝐿

T

where the asterisk denotes complex conjugation and 𝑆𝑅 and 𝑆𝐼 are the real and imaginary components of the complex image 𝑆𝑙 . If the variance of the noise at each coil is the same and there are no correlations (𝚺 = 𝜎 2 𝐈), then the magnitude 𝑀 follows the Nc-Chi distribution [18] with a Probability Density Function (PDF) given by [19], [20], [18] (3)

all the coils and is given by 𝐴𝑇 = (∑𝐿𝑙=1 𝑐𝑙 𝐴2 ) 𝐴

1

SC

where 𝐼𝐿−1 (⋅) stands for the modified Bessel function of the first kind of order (𝐿 − 1), 𝐴𝑇 is the total image voxel intensity in the absence of noise taking into account the contribution of 2=𝐴∥𝑐∥2

where 𝑐 = [𝑐1 , ⋯ , 𝑐𝐿 ]ú , 𝑥 =

𝐴𝑇 𝑆𝑇 𝜎2

and

PT

Tr(𝚺) 𝐿eff

,

(4)

(5)

AC

2 𝜎eff =

𝐴2𝑇 Tr(𝚺)+(Tr(𝚺))2 , 2 ̅ú 𝚺𝐀 ̅ +𝚺𝐹 𝐀

CE

𝐿eff =

ED

MA

NU

𝜂 = 𝜎𝑇 is the SNR. The PDF given by Eq. ((3)) is known also as the generalized Rice distribution, since the Rice distribution corresponds to the case 𝐿 = 1. Interestingly, the Nc-Chi distribution applies to a larger class of dMRI signals including the correlated noise case and parallel MRI image reconstruction cases, as explained in the following paragraphs. Correlated noise case: We have assumed previously that the noise terms are of equal variance and are uncorrelated. However, in multi-coil systems noise correlations may exist as highlighted in [21], [17], [22]. Taking noise correlations into account, Aja-Fernández et al. [23], [24], [25] showed that the Nc-Chi distribution is still an accurate approximation of the obtained signal PDF, but with the following “effective" parameters:

̅ = [𝐴1 , … , 𝐴𝐿 ]. In with ⋅𝐹 is the Frobenius norm and Tr(⋅) stands for the trace operator, 𝐀 this case, the PDF given by Eq. ((3)) may be written as: 𝒫(𝑀) = with 𝑥 =

𝑥2

2

𝜂 1−2𝐿eff − eff 𝐿eff −2𝜂2 2 𝑥 𝜂 𝑒 𝑒 eff 𝐼𝐿eff −1 (𝑥), 𝜎eff eff

1

𝐴𝑇 𝑀 2 𝜎eff

and 𝜂eff =

𝐴𝑇 𝜎eff

(6)

is the effective SNR. As shown by Eq. ((4)) and Eq. ((5)),

the noise correlation reduces the effective number of coils 𝐿eff ≤ 𝐿 (where 𝐿 is the actual number of coils) and increases the effective variance of the noise, 𝜎eff ≥ 𝜎 with equality if the noise is uncorrelated (spatially white) with equal variance. Parallel MRI (pMRI) case: The current trend in multi-coil systems is the undersampled acquisition of the k-space known as parallel MRI. This results in a shortened scan time due to the reduction of the number of gradients. Since the beginning of clinical use of pMRI, various different algorithms have been developed by researchers to solve the problem of aliased reconstructed images resulting from undersampled k-space. The most famous are the SENSitivity Encoding (SENSE) [26] and the GeneRalized Autocalibrating Partially Parallel Acquisition

ACCEPTED MANUSCRIPT

2.2 Diffusion Signal Model

NU

SC

RI P

T

(GRAPPA) [27]. However, other techniques have also been proposed (for a review see [28], [29] ). For the GRAPPA algorithm, missing k-space lines are first estimated (interpolated) for each receiver channel to fill the whole k-space, resulting in particular in correlated noise terms. In the second step, the CMS is calculated using the SoS of the resulting complex signal in each coil (Eq. ((2))) resulting in Nc-Chi distribution for the CMS [23] as explained in the previous paragraph. The SENSE algorithm directly reconstructs one complex final image without going through intermediate coil images like GRAPPA. The intermediary reconstructed signal will follow a complex Gaussian distribution. Thus, the final magnitude image has a Rician distribution which is a particular case of Nc-Chi distribution for 𝐿 = 1. In brief, the Nc-Chi distribution is a valid model for a large range of situations. Consequently, eventhough presented for the SoS reconstruction of multi-channel image data, the proposed study has a wider impact and can be adapted and used in many other contexts with appropriate values for 𝐿eff and 𝜎eff . For convenience, in the remainder of this paper we drop the subscript (eff) from our notations.

PT

ED

MA

Different representations are considered in the literature to model the diffusion signal in the human brain white matter (see [5], [6]). In this work, we use the multi-tensor model to deal with crossing fibers. We assume that the diffusion in a single fiber bundle is mono-exponential and can be described by an anisotropic Gaussian profile. The signal intensity in the noiseless case, taking into account the contribution of all scanner coils, is then a finite mixture of signals, described by the Stejskal-Tanner equation [30] ú

(7)

CE

−𝑏𝑗 𝑔𝑘𝑗 𝐷𝑖 𝑔𝑘𝑗 𝐴𝑘𝑗 = 𝐶𝐿 𝐴0 ∑𝑚 , 𝑖=1 𝑓𝑖 𝑒

AC

where 𝑚 is the number of fiber-bundles, 𝐶𝐿 =∥ 𝑐 ∥2 , 𝐴0 is the signal intensity without diffusion weighting, i.e. when 𝑏𝑗 = 0 , 𝐠 𝑘𝑗 = [𝑔𝑥 , 𝑔𝑦 , 𝑔𝑧 ]ú is a unit vector representing the coordinates of the 𝑘 𝑡ℎ diffusion gradient direction applied using the b-value 𝑏𝑗 , which includes the main parameters of the diffusion sequence. In Eq. ((7)), 𝑗 refers to the shell index and 𝑘 refers to the diffusion gradient index (see Fig. 1). 𝑓𝑖 is the volume fractions of occupancy of the 𝑖 𝑡ℎ fiber-bundle with ∑𝑚 𝑖=1 𝑓𝑖 = 1. 𝐷𝑖 are 3 × 3 positive definite symmetric matrices representing the diffusion tensors associated with the 𝑖 𝑡ℎ fiber-bundles.

2.3 Parameters of Clinical Interest Most tractography algorithms employ the directional information contained in the eigenvectors of the diffusion tensor, to approximate the white matter fiber orientation. Therefore, the first parameter of clinical interest is the principal eigenvectors of 𝐷𝑖 , 𝑖 = 1, ⋯ , 𝑚. The second parameter of clinical interest is the anisotropy index. Indeed, in the literature, many factors were defined to measure the anisotropy of the brain’s tissues including the Fractional Anisotropy (FA), the Relative Anisotropy (RA), the Mean Diffusivity (MD), the Volume Ratio (VR) first described in [16], the Ellipsoidal Area Ratio (EAR) [31] and the MODe of anisotropy (MOD) [32]. They all quantify the degree to which the observed diffusion is direction-dependent. It is of great interest to evaluate and analyze the impact of the acquisition parameters (b-value, SNR, number of diffusion gradients and number of coils) on the estimation precision of these factors, especially when they

ACCEPTED MANUSCRIPT

RI P

T

are used for clinical diagnosis [33]. This is one of the objectives of our work. Note that the previous anisotropy indices have been defined for the DT model and their extensions to the multi-tensor model need to be considered. In the DT model, voxels with two or more fiber-bundles will have a lower FA than voxels with one single fiber-bundle (Fig. 2-(a)) and therefore, the loss of some fibers will increase the FA, which might be confusing. This has been highlighted by [34], where it is shown that high FA does not necessarily point out to the existence of a single fiber-bundle, as it might reflect the degeneration of one of the fiber-bundles in crossing areas. To overcome this limitation we suggest to extend the FA definition to the MTM as (8)

SC

MT − FA = ∑𝑚 𝑖=1 𝑓𝑖 FA𝑖 ,

where, MT-FA stands for Multi-Tensor FA. Indeed, this definition preserves the properties of the original fractional anisotropy index FA1 .

NU

since:

• It includes the DT model case corresponding to 𝑓1 = 1 for which MT − FA = FA =

CE

PT

ED

MA

• If the diffusion tensors are highly anisotropic, i.e. FA1 ≈ FA2 , ⋯ , ≈ FA𝑚 ≈ 1, then the global index will reflect this situation, i.e. MT − FA ≈ 1. • If the diffusion tensors are highly isotropic, i.e. FA1 ≈ FA2 , ⋯ , ≈ FA𝑚 ≈ 0, then the global index will satisfy MT − FA ≈ 0. Moreover, this definition has the advantage of providing the appropriate anisotropy index value (i.e. MT-FA close to one) in crossing fiber regions (90% of the white matter voxels contain two or more fiber orientations [35]) when the standard definition based on the DT model may lead to wrong values. Also, the FA in crossing fiber voxels depends on the crossing angles (see Fig. 2-(a)), while our proposed factor does not (see Fig. 2-(b)). In the same way, we define the Multi-Tensor Mean of Diffusivity (MT-MD) as MT − MD = ∑𝑚 𝑖=1 𝑓𝑖 MD𝑖 .

(9)

AC

From Eqs. ((8)) and ((9)), loss of fibers in crossing area, will decrease the MT-FA and MT-MD. This property can be used as a biomarker to identify changes to white matter that are specific to certain neurodegenerative diseases like the Alzheimer disease [3].

2.4 Cramér-Rao Lower Bound Our main objective in this study is to derive theoretical bounds for the estimation error variance of the desired model parameters we need to obtain from the set of observed data, and analyze the impact of certain system parameters on the overall estimation performance. To achieve this objective, we rely on the Cramér-Rao Bound (CRB) theory, briefly reviewed here. The CRB provides fundamental limits to the precision with which a parameter can be determined from any experimental data using unbiased estimation methods. Indeed, for a given statistical model depending on a parameter vector 𝚯, the mean squares error of any unbiased estimator of 𝚯 is lower bounded by the CRB computed as the inverse of the FIM. The likelihood function is denoted by 𝐿(𝑆, 𝚯), where 𝑆 is the observation vector. The (𝑖, 𝑗)𝑡ℎ entry of the FIM matrix is given by [36]:

ACCEPTED MANUSCRIPT 𝐹𝑖,𝑗 (Θ) = −𝔼[(

∂2 ln(𝐿(𝑆,Θ)) ∂𝜃𝑖 ∂𝜃𝑗

)]

(10)

𝐶𝑜𝑣(𝜉) ≥ CRB(𝜉) = ∇𝜉(Θ)CRB(Θ)(∇𝜉(Θ))ú ,

RI P

T

where 𝔼[⋅] stands for the statistical expectation operator. If 𝜉̂ represents an unbiased estimator of 𝜉(Θ) a function of Θ, then its error covariance matrix is greater than or equal to the CRB, i.e.

SC

with CRB(Θ) = 𝐹 −1 (Θ) and ∇𝜉(Θ) is the gradient matrix which (𝑖, 𝑗)𝑡ℎ entry is

(11) ∂𝜉𝑖 ∂𝜃𝑗

. In

NU

particular, this means that the error variance in the estimation of parameter 𝜉𝑖 is greater than the 𝑖 𝑡ℎ diagonal element of the CRB matrix in Eq. ((11)).

2.5 FIM Expression for the Multi-Tensor Model

𝐾

∂𝜂𝑘𝑗 ∂𝜂𝑘𝑗

CE

𝑗 F𝑚,𝑛 = ∑𝑁 𝑗=1 ∑𝑘=1

PT

ED

MA

In this subsection, we derive the expression of the FIM and CRB matrices for the model in Eqs. ((6)) and ((7)). For that, we assume that the dMRI measurements are statistically independent and, for simplicity, we assume the noise power 𝜎 2 known. We consider 𝑁 sets of diffusion gradient directions {𝐠 𝑘𝑗 , 𝑗 ∈ [1, ⋯ , 𝑁],1 ≤ 𝑘 ≤ 𝐾𝑗 } . 𝐾 = ∑𝑁 1 𝐾𝑗 is then the total number of diffusion gradients used for the acquisition of the diffusion MRI data, 𝑆 = [𝑆1ú , ⋯ , 𝑆𝑁ú ]ú where 𝑆𝑗ú denotes a vector of 𝐾𝑗 measurements acquired at b-value 𝑏𝑗 , 𝐾𝑗 being the total number of diffusion gradients of the 𝑗 𝑡ℎ shell (see Fig. 1). The general expression of the FIM in this context (multi-shell and multi-coil system when using the SoS) is given by

AC

with

∂𝜃𝑚 ∂𝜃𝑛

𝑥 2 𝐼𝐿2 (𝑥𝑘𝑗 )

ℳ𝑘𝑗 ,

(12)

2 ] − 𝜂𝑘𝑗 ,

ℳ𝑘𝑗 = 𝔼 [𝜂2𝑘𝑗𝐼2

(13)

𝑘𝑗 𝐿−1 (𝑥𝑘𝑗 )

𝐴

𝜂𝑘𝑗 = 𝜎𝑘𝑗 where 𝐴𝑘𝑗 is the intensity value at the 𝑘 𝑡ℎ measurement of the 𝑗 𝑡ℎ shell given by Eq. ((7)). Using Eq. (6), ℳ𝑘𝑗 can be expressed as an integral formula −2(𝐿+1)

ℳ𝑘𝑗 = 𝜂𝑘𝑗 The

𝑒

parameter

𝜂2 𝑘𝑗 − 2

+∞ ∫0

𝑥2 2𝜂2 𝑘𝑗

𝐼𝐿2 (𝑥) 𝐼𝐿−1 (𝑥)

2 𝑑𝑥 − 𝜂𝑘𝑗 .

(14)

vector in the case of the multi-tensor model is (𝑖) (𝑖) where ( 𝑑𝑖 = [𝑑1 , … , 𝑑6 ]ú is the six-element vector

ú [𝑑1ú , ⋯ , 𝑑𝑚 , 𝑓1 , ⋯ , 𝑓𝑚−1 ]ú

𝚯= rearrangement of

𝑥 𝐿+2 𝑒

−

ACCEPTED MANUSCRIPT

𝐷𝑖 = (𝑑4(𝑖) (𝑖) 𝑑6

(𝑖)

𝑑4

(𝑖)

𝑑2

(𝑖) 𝑑5

(𝑖)

𝑑6

(𝑖) 𝑑5 ) ; 𝑖 = 1, ⋯ , 𝑚,

(15)

(𝑖) 𝑑3

to be estimated from the observed data vector 𝑆. We define ú

RI P

−𝑏𝑗 𝐠̃ 𝑘𝑗 𝐝𝑖 𝜂𝑘𝑗 = 𝐶𝐿 𝜂0 ∑𝑚 , 𝑖=1 𝑓𝑖 𝑒

T

(𝑖)

𝑑1

𝐴0

(16)

is the SNR for no diffusion weighting (𝑏 = 0). 𝐠̃ is obtained from the rewriting of the quadratic form 𝐠 ú 𝐃𝐠 as a dot product between two vectors as follows 𝐠 ú 𝐃𝐠 = 𝐠̃ ú 𝐝, i.e. 𝐠̃ = [𝑔𝑥2 , 𝑔𝑦2 , 𝑔𝑧2 , 2𝑔𝑥 𝑔𝑦 , 2𝑔𝑦 𝑔𝑧 , 2𝑔𝑥 𝑔𝑧 ]ú. For notational simplicity1, we consider in the sequel the bi-tensor case, for which the parameter vector reduces to 𝚯 = [𝑑1ú , 𝑑2ú , 𝑓]ú , where 𝑓1 = 𝑓 and 𝑓2 = (1 − 𝑓). Using Eq. ((12)), the FIM can be expressed in matrix form as: 𝐹 = ∑𝑁 𝑗=1 𝐹𝑗 , where 𝐹𝑗 is a symmetric matrix with an upper bloc-triangular part given by 𝜎

2 ú 𝑏𝑗2 𝐆𝑗 𝚼𝑗 𝚼1𝑗 𝐆𝑗

𝑏𝑗2 𝐆𝑗 𝚼𝑗 𝚼1𝑗 𝚼2𝑗 𝐆𝑗ú

𝑏𝑗 𝐆𝑗 𝚼𝑗 𝚼1𝑗 𝐯𝑗

2 ú 𝑏𝑗2 𝐆𝑗 𝚼𝑗 𝚼2𝑗 𝐆𝑗

𝑏𝑗 𝐆𝑗 𝚼𝑗 𝚼2𝑗 𝐯𝑗 ).

ED

𝐹𝑗𝑈 = 𝐶𝐿2 𝜂02 (

MA

NU

SC

where 𝜂0 =

(17)

𝐯𝑗ú 𝚼𝑗 𝐯𝑗

PT

𝐆𝑗 is a 6 × 𝐾𝑗 matrix derived from the 𝐾𝑗 diffusion gradient components associated with 𝑏𝑗 given by 𝐆𝑗 = [𝐠̃1𝑗 , ⋯ , 𝐠̃ 𝐾𝑗𝑗 ]. 𝚼𝑗 , 𝚼1𝑗 , 𝚼2𝑗 are 𝐾𝑗 × 𝐾𝑗 diagonal matrices, given by = diag(ℳ1𝑗 , ⋯ , ℳ𝐾𝑗𝑗 ).

𝚼𝑖𝑗

= 𝑓𝑖 diag(𝑒 −𝑏𝑗𝐠̃1𝑗𝐝𝑖 , ⋯ , 𝑒

ú

= [𝑒

−𝑏𝑗 𝐠̃ ú1𝑗 𝐝2

AC

𝐯𝑗

CE

𝚼𝑗

−𝑒

−𝑏𝑗 𝐠̃ ú1𝑗 𝐝1

−𝑏𝑗 𝐠̃ ú𝐾 𝑗 𝐝𝑖 𝑗

,⋯,𝑒

), 𝑖 = 1,2.

−𝑏𝑗 𝐠̃ ú𝐾 𝑗 𝐝2 𝑗

−𝑒

−𝑏𝑗 𝐠̃ ú𝐾 𝑗 𝐝1 ú 𝑗

(18)

] .

Remark: In the case of the diffusion tensor model (DT), the signal model is given by ú

𝐴𝑘𝑗 = 𝐶𝐿 𝐴0 𝑒 −𝑏𝑔𝑘𝑗 𝐷𝑔𝑘𝑗 .

(19)

The parameter vector in this case is 𝚯 = 𝐝 = [𝑑1 , … , 𝑑6 ]ú and the FIM reduces to 2 ̃2 ú 𝐹 = 𝐶𝐿2 𝜂02 ∑𝑁 𝑗=1 𝑏𝑗 𝐆𝑗 𝚼𝑗 𝚼𝑗 𝐆𝑗 ,

(20) ú

̃𝑗 = diag(𝑒 −𝑏𝑗𝐠̃1𝑗𝐝 , ⋯ , 𝑒 −𝑏𝑗𝐠̃𝐾𝑗𝑗𝐝 ). where 𝚼 ú

2.6 CRB Expression for Parameters of Clinical Interest 1

Except for a more complex and heavier notational system, there is no theoretical limitation or difficulty to deal with 𝑚 > 2.

ACCEPTED MANUSCRIPT Based on the CRB theory, the unbiased estimation of the diffusion anisotropy index for a given tensor 𝐷 cannot be obtained with an error variance lower than: CRB(FA) = ∇(FA) CRB(𝜆) ∇(FA)ú ,

T

(21) ∂(FA) ∂(FA) ∂(FA)

RI P

where ∇(FA) is the gradient vector of the anisotropy factor given by [ with

∂𝜆1

,

∂𝜆2

,

∂𝜆3

],

SC

1 2 ∂(FA) 𝜆𝑖 (1 − FA ) − 2 ∑𝑞≠𝑖 𝜆𝑞 = , FA ∑𝑞 𝜆2𝑞 ∂𝜆𝑖

and

NU

CRB(𝜆) = ∇(𝜆) CRB(𝐝) ∇(𝜆)ú ,

(22)

MA

where 𝜆 = [𝜆1 , 𝜆2 , 𝜆3 ]ú and CRB(𝐝) denotes the CRB of the tensor under consideration given by its corresponding 6 × 6 bloc of the inverse of the FIM computed in Eq. ((12)). ∇(𝜆) denotes the 3 × 6 gradient matrix computed in appendix A, and given by (23)

ED

∇(𝜆) = (𝐄e𝐄)ú 𝐽.

PT

The symbol e being the Khatri-Rao product and 𝐽 is a 9 × 6 selection matrix defined by 𝑣𝑒𝑐(𝐷) = 𝐽𝑑. 𝐄 is the 3x3 eigenvector matrix of tensor 𝐷. For the 𝑖 𝑡ℎ eigenvector of 𝐷, its unbiased estimation error covariance is lower bounded by:

CE

CRB(𝐞𝑖 ) = ∇(𝐞𝑖 ) CRB(𝑑) ∇(𝐞𝑖 )ú .

(24)

AC

The first-order perturbation theory allows us to compute the eigenvector derivative as shown in appendix A CRB(𝐞𝑖 ) = (𝐞ú𝑖 ⊗ 𝚪𝑖 ) 𝐉 CRB(𝑑) 𝐉 ú (𝐞𝑖 ⊗ 𝚪𝑖 ),

(25)

𝚪𝑖 = (𝐃 − 𝜆𝑖 𝐈)# , where (#) denotes the pseudo inverse and ⊗ refers to the Kronecker product. Using the Frobenius norm, the mean square error achievable by any unbiased estimator 𝐞̂𝑖 of the eigenvector 𝐞𝑖 satisfies: MSE(𝐞̂𝑖 ) = 𝔼{∥ 𝐞̂𝑖 − 𝐞𝑖 ∥2𝐹 }Tr(CRB(𝐞𝑖 )),

(26)

2.7 Approximate Analytical FIM Expression Unfortunately, the integral in Eq. ((14)) does not admit a closed-form expression. For this reason, the authors in [12], [13], [14], [15] rely on numerical evaluation of these integral terms. Using MATLAB, the Bessel function 𝐼𝐿 (𝑥) is available under the command 𝑏𝑒𝑠𝑠𝑒𝑙𝑖(𝐿, 𝑥) but its use is prone to numerical errors due to the large dynamical range of Bessel function values. Fortunately, the scaled version of the Bessel function in MATLAB under the command

ACCEPTED MANUSCRIPT 𝑏𝑒𝑠𝑠𝑒𝑙𝑖(𝐿, 𝑥, 1) equal to 𝐼̆𝐿 (𝑥) = 𝐼𝐿 (𝑥)𝑒 −𝑥 helps to mitigate this numerical issue according to 2

ℳ𝑘𝑗 =

𝜂 −2(𝐿+1) − 𝑘𝑗 ∞ 𝜂𝑘𝑗 𝑒 2 ∫0

𝑥

𝐿+2

𝑒

𝑥−

𝑥2 2𝜂2 𝑘𝑗

𝐼̆𝐿2 (𝑥)

𝐼̆𝐿−1 (𝑥)

2 𝑑𝑥 − 𝜂𝑘𝑗 .

(27)

RI P

T

To avoid this integral calculus (cumbersome and prone to numerical errors), and to get a more tractable FIM formula, we propose a very simple approximate analytical expression for the expectation in Eq. ((13)), where the main steps are given in this subsection. Thus, our approach is based on the decomposition of the square ratio of the Bessel function of the first kind (𝑓(𝑥, 𝐿) = 𝐼𝐿 (𝑥)

appendix B for details) 𝑓𝐿𝑜𝑤 (𝑥, 𝐿) ≈

𝑥2 4𝐿2

𝑥4

(5𝐿+6) 𝑥 6

SC

𝐿−1

2

) ) into two analytical series expressions for low and high SNRs, respectively (see (𝑥) (7𝐿+12) 𝑥 8

(21𝐿3 +118𝐿2 +214𝐿+120) 𝑥 10

− 8𝐿3 (𝐿+1) + 64𝐿4 (𝐿+1)2 (𝐿+2) − 128𝐿5 (𝐿+1)2 (𝐿+2)(𝐿+3) + 512𝐿6 (𝐿+1)3 (𝐿+2)2 (𝐿+3)(𝐿+4) + ⋯ (2𝐿−1)

𝑓𝐻𝑖𝑔ℎ (𝑥, 𝐿) ≈ 1 −

+

NU

(𝐼

(2𝐿−1)(2𝐿−3) 𝑥2

−

(2𝐿−1)(2𝐿−3)2

MA

𝑥 (2𝐿−7)(2𝐿−5)(2𝐿−3)(2𝐿−1)(2𝐿+3) 128 𝑥 5

(28)

−

+

(𝐿−2)(2𝐿−3)(2𝐿−1)

+

4 𝑥4 (𝐿−3)(2𝐿−3)(2𝐿−1)(−9+4𝐿(𝐿−2)) 8 𝑥3

16 𝑥 6

+⋯

𝜂2 𝐿

ℳHigh ≈ 1 −

(2𝐿+3) 𝜂 6

−

(𝐿+3) 𝜂 8

+⋯

2𝐿3 (𝐿+1) 𝐿4 (𝐿+1) (2𝐿−1) (2𝐿−3)(2𝐿−1) (𝐿−3)(2𝐿−3)(2𝐿−1)

+

2𝜂 2

4𝜂 4

−

4𝜂 6

+

(2𝐿−3)(2𝐿−1)(4𝐿(𝐿−8)+57) 16𝜂 8

+⋯

CE

(29)

𝜂4

− 𝐿2 +

PT

ℳLow ≈

ED

We inject the above series expansions into Eq. ((13)) and by applying the general formula of the moments of the Nc-Chi distribution [37], we obtain the following expressions:

AC

where ℳLow , ℳHigh refer to the approximate expressions of ℳ at low and high SNR, respectively. For large values of 𝐿, the previous expressions can be further simplified to: ℳLow ≈ − ∑𝑞𝑖=1 (−

𝜂2

𝐿

𝑖

𝑖

) + 𝑂(𝜂2 )𝑞+1 , 𝐿 1 𝑞

ℳHigh ≈ ∑𝑞−1 𝑖=0 (− 𝜂 2 ) + 𝑂 (𝜂 2 ) ,

(30) (31)

𝑞 denotes here the series expansion order. When 𝑞 goes to infinity, both expressions lead to the same limit value 𝐿

−1

ℳ ≈ (1 + 𝜂2 ) .

(32)

In Fig. 3 we show a graphical comparison between the expected value ℳ based on the integral formula given by Eq. ((14)) and our expression given by Eq. ((32)). The two curves are almost identical. However, for moderate values of 𝐿 (i.e. 𝐿 < 5). The resulting error in the bending area is not null but can be neglected (more details can be found in appendix B). We draw the reader’s attention that modern clinical scanners are usually equipped with more than 4 acquisition coils. Thus the evaluation of the integral given by Eq. ((14)), is practically reduced to the evaluation of the expression given by Eq. ((32)).

ACCEPTED MANUSCRIPT By plugging this expression into Eq. ((12)) we get the approximate analytical FIM formula 𝐾

𝑗 F𝑚,𝑛 ≈ ∑𝑁 𝑗=1 ∑𝑘=1

∂𝜂𝑘𝑗 ∂𝜂𝑘𝑗 ∂𝜃𝑚 ∂𝜃𝑛

𝐿

−1

(1 + 𝜂2 ) .

(33)

𝑘𝑗

T

As shown in appendix B, this FIM approximation provides an excellent fit to the exact value given by Eq. ((12)).

RI P

2.8 Performance analysis

−1

NU

𝐿

SC

Based on the CRB expressions, we propose here a theoretical investigation of the performance impact of the coils number, the SNR and the b-value, respectively. 2.8.1 Effect of the Number of Coils For multi-channel systems, the image reconstruction technique combines information from different acquisition coils. From Eq. ((32)), ℳ can be written as ℳ ≈ (1 + 𝐶 2 𝜂̃2 ) , 𝐿

(34) 𝐴

𝐶𝐿2

𝐿

≈ 1, which might occur for example if 𝑐𝑙 ≈ 1, ∀𝑙, then ℳ will be 1

−1

ED

we take

𝐿

MA

where 𝜂̃ is the SNR of the diffusion-weighted MRI data normalized by 𝐶𝐿 i.e. 𝜂̃ = 𝜎𝐶 . If

ℳ ≈ (1 + 𝜂̃2 ) .

(35)

CE

PT

ℳ no longer depends on 𝐿, and the FIM given by Eq. ((17)) for the bi-tensor model and Eq. ((20)) for the DT model will be of the form: FIM ≈ 𝐶𝐿2 FIM0 ≈ 𝐿 FIM0 and hence the CRB CRB = 𝐿 0 . This means that the estimation error variance will be reduced by a factor 𝐿 (as compared to the single coil system) when considering an acquisition system using 𝐿 coils.

AC

2.8.2 Effect of the SNR Diffusion-Weighted Images (DWI) have an inherently low SNR, which makes DWI analysis complicated and prone to errors in clinical interpretation. Therefore in [38] the authors suggest that the SNR should never be below 3 in any of the DWI. Descoteaux et al. [39] maintain the SNR greater than 4 along the most attenuated direction to limit the Rician noise bias in the data. Since the SNR cannot be set by the user, (it depends on several parameters: b-value, the direction of the diffusion gradient, the medium ... ), in the following we study the influence of the intrinsic SNR (denoted 𝜂0 ) (i.e. without diffusion weighting) on parameters estimation using the ú

CRB. From Eq. ((32)) we can write: ℳ𝑘 = (1 + ú

(1 +

ú

(𝑓𝑒 −𝑏𝐠̃𝑘 𝐝1 +(1−𝑓)𝑒 −𝑏𝐠̃ 𝑘 𝐝2 )−2 𝜂02

𝑒 2𝑏𝐠̃𝑘 𝐝 𝜂02

−1

)

for the DT model and ℳ𝑘 =

−1

)

for the bi-tensor model (with the assumption that

𝐿 𝐶𝐿2

≈ 1). We

notice that for large 𝜂0 , ℳ𝑘 will go to one, and 𝚼𝑗 will be close to the identity matrix. In that case, using Eq. ((17)), the CRB can be written in the form CRB≈

CRB0 𝜂02

. This means that the

estimation error variance decreases with a quadratic rate with the SNR value (for example, if the SNR value is doubled the estimation MSE will be reduced by a factor 4).

ACCEPTED MANUSCRIPT

MA

NU

SC

RI P

T

2.8.3 Effect and Optimal Selection of b-value It is well known that higher b-values require long echo time (TE) and the study of the effect of the b-value on minimum TE was previously performed (see [40], [11]) which is beyond the scope of this work. In this subsection, we study only the impact of the b-value on the desired parameter estimation precision. We refer to the b-value that minimizes the CRB expression as the “optimal b-value". More precisely, our purpose is to seek for the appropriate b-value and investigate the effect of the SNR on the choice of this value when considering first an isotropic then an anisotropic tensor. Unfortunately, for the bi-tensor model, we do not have a closed-form expression for the optimal 𝑏-values and consequently we will assess its effect in the next section through numerical simulation. In fact, as shown in the following, the b-values selection depends on several parameters and might be optimized only for specific applications (or diagnoses purposes) as done in [41]. Thus, we consider here the single tensor case. It has been shown that the optimal diffusion weighting for diffusivity measurements is approximately 1.11/ADC [42], (ADC denotes the Apparent Diffusion Coefficient) and 1.1 × 3/Tr(𝐷) for the anisotropic tensor [43]. To check this, we will use the CRB to assess the optimal b-values (𝑏𝑜𝑝𝑡 ) that minimize the CRB. In the sequel, we analyze the CRB for the mono-shell case (for simplicity we drop the index 𝑗). When considering an isotropic tensor 𝐃 = 𝜆𝐈, with 𝐈 being the identity matrix. The expression of the FIM given by Eq. ((20)) will be

ED

𝐿

−1

FIM = 𝐿𝜂02 𝑏 2 𝑒 −2𝑏𝜆 (1 + 𝐶 2 𝜂2 𝑒 2𝑏𝜆 ) 𝐿 0

𝐆𝐆ú ,

(36)

PT

and the CRB expression in this case is 1

𝐿

0

CE

CRB = 𝐿𝜂2 𝑏2 𝑒 2𝑏𝜆 (1 + 𝐶 2 𝜂2 𝑒 2𝑏𝜆 ) 𝐇, 𝐿 0

(37)

AC

where 𝐇 = (𝐆𝐆ú )−1. To find the optimal value of 𝑏 that minimizes Eq. ((37)), the first 1 𝐿 derivative of 𝑏2 𝑒 2𝑏𝜆 (1 + 2 2 𝑒 2𝑏𝜆 ) with respect to 𝑏 is set equal to zero leading to 1+𝛼

𝐶𝐿 𝜂0

1

𝐿

𝑏𝑜𝑝𝑡 = (1+2𝛼) 𝜆 , with 𝛼 = 𝐶 2 𝜂2 𝑒 2𝑏𝑜𝑝𝑡𝜆 . 𝐿 0

(38)

For high values of the SNR 𝜂0 , 𝛼 will be close to zero and the optimal b-value will be 1 approximately 𝑏𝑜𝑝𝑡 ≈ 𝜆 . For very low SNR 𝜂0 , 𝛼 will be greater than one and the optimal 1

b-value is close to 𝑏𝑜𝑝𝑡 ≈ 2𝜆. So, the optimal b-value for an isotropic tensor does not depend on the number of diffusion gradients but depends strongly on the isotropic diffusion value and the 1 1 SNR. Which means that the optimal b-value will be in the range 2𝜆 ≤ 𝑏𝑜𝑝𝑡 ≤ 𝜆. In the case of an anisotropic tensor, it is difficult to get an exact analytical expression for 𝑏𝑜𝑝𝑡 , due to its dependency on the SNR, the intrinsic diffusivities, the number of diffusion gradients, the tensor orientation and the parameters to be optimized (i.e. CRB of 𝜆𝑖 , 𝐞𝑖 , FA, RA, ...). To approach the optimal b-value, we will replace the anisotropic tensor by its closest isotropic tensor 𝐃𝑖𝑠𝑜 = 𝜆𝐈, which allows us to use the formula given in Eq. ((38)). For that, we will decompose the tensor 𝐃 into its isotropic and anisotropic parts 𝐃 =

ACCEPTED MANUSCRIPT 𝐃𝑖𝑠𝑜 + 𝐷𝑎 , where 𝐷𝑎 = (𝐃 − 𝜆𝐈) is an anisotropic tensor. 𝜆 is chosen so that to minimize the Tr(𝐃) distance between 𝐃 and 𝐃𝑖𝑠𝑜 . Thus, 𝜆 = 3 . By plugging this result in Eq. ((38)), the optimal b-value for an anisotropic tensor will be in the range 3

T

≤ 𝑏𝑜𝑝𝑡 ≤ Tr(𝐷).

(39)

RI P

3 2Tr(𝐷)

3 Results

ED

MA

NU

SC

In this section, we exploit the developed CRB to analyze numerically the impact of the controllable system parameters on the estimation of the bi-tensor parameters. We consider the relative error measure of the desired parameter 𝜃 (𝜃) STDmin 𝑒𝜃 = × 100, 𝜃 (⋅) STDmin is the minimum standard deviation in the estimation of the parameter under consideration obtained from its CRB. Unless otherwise stated, we simulated two orthogonally crossing fiber bundles using the model given in Eq. ((7)) with, [𝜆1𝑖 , 𝜆𝑖2 , 𝜆𝑖3 ] = [1708,303,114] × 10−6 𝑚𝑚2 . 𝑠 −1 , 𝑖 = 1,2 diffusivities typical values taken from [44]. The number of acquisition coils is assumed constant 𝐿 = 8, 𝑓 = 0.5 and the SNR 𝜂0 is set equal to 30.

3.1 Effect of the Number of Diffusion Gradients per Shell

𝐾

CE

PT

In [45], the authors proposed a design method for the multi-shell sampling protocol, with a uniform coverage on each shell and a global uniform angular coverage. Here, we use this protocol to study the impact of the number of diffusion gradients per shell on the estimation 𝐾 𝐾 precision. Three configurations from [45] are considered corresponding to (𝐾1 ≈ 2 ), (𝐾1 ≈ 3 ),

AC

and (𝐾1 ≈ 5 ), where 𝐾 denotes the total number of diffusion gradients. The plots from Fig. 4 for a fixed total number of diffusion gradients equal to 𝐾 = 120 and 𝑏1 = 1000 𝑠. 𝑚𝑚−2 show that the diffusion gradients distribution scheme has a non negligible impact on the estimation precision. Moreover, the optimal configuration depends on the tensors shape and hence is not unique . The optimal choice of the diffusion gradient distribution is an open problem which might be achieved only for specific objectives and application context. We adopt in the sequel the third 𝐾 configuration corresponding to (𝐾1 ≈ 5 ).

3.2 Influence of the Number of Diffusion Gradients From Fig. 5, we observe a large gain when varying the number of diffusion gradients until close to 50 after which the gain becomes less and less significant. In a clinical routine, the total scan time must be as short as possible, due to some pathologies that may prevent the patient from staying for long time periods. Hence, to reduce the scan time while preserving a good estimation precision, it is preferable to increase the number of acquisition coils rather than the number of diffusion gradients. For example, from Fig. 5 we can see that using 16 coils and about 50 diffusion gradient directions is as good as working with 4 coils and 185 diffusion gradients. Note that we have observed such behaviour (i.e. fast initial decrease of the error function followed by a

ACCEPTED MANUSCRIPT saturation effect) in all simulation contexts we have conducted.

3.3 Optimal b-values

MA

NU

SC

RI P

T

For the bi-tensor case, it was shown in [11] that at least two b-values are necessary for the estimation of all tensor parameters. In this subsection we aim to find the b-value pair (𝑏1 , 𝑏2 ) with 𝑏1 < 𝑏2 that minimizes the CRB. We keep 𝑏1 and 𝑏2 ≤ 3000 𝑠. 𝑚𝑚−2 to preserve the Gaussian diffusion assumption [14], [11]. In this first experiment, the results in Figs. 6-(a) and 6-(b) represent the relative error measure for the MT-FA index defined in Eq. ((8)) and the principal tensor direction (𝑒1 ) of 𝐷1 , respectively, with respect to the values of 𝑏1 and 𝑏2 . We can see that the minimum error for MT-FA is reached for 𝑏1 ∈ [300,1500] 𝑠. 𝑚𝑚−2 and 𝑏2 ∈ [2100,3000] 𝑠. 𝑚𝑚−2 while the minimum error for 𝑒1 is obtained for 𝑏1 ∈ [300,3000] 𝑠. 𝑚𝑚−2 and 𝑏2 ∈ [2500,3000], 𝑠. 𝑚𝑚−2 . In Fig. 7, we consider a “quasi isotropic" second tensor given by [𝜆12 , 𝜆22 , 𝜆23 ] = [1685,1528,1410] × 10−6 𝑚𝑚2 . 𝑠 −1 . We observe performance degradation and the range of optimal b-values is different from the previous situation. This illustrates the fact that the optimum pair is not unique but depends on the target parameter, the context, and eventually the target cost function. However, in most cases we have observed that the pair (𝑏1 , 𝑏2 ) ≈ (800,2000) 𝑠. 𝑚𝑚−2 represents a good ad-hoc choice for the model parameter estimation.

ED

3.4 Effect of the Crossing Angle

CE

4 Discussion

PT

In the following, the estimation performance is investigated with respect to the crossing angle. We keep 𝑏1 = 1000 𝑠. 𝑚𝑚−2 and we vary 𝑏2 in the range [1100,3000] 𝑠. 𝑚𝑚−2 . Fig. 8 shows that the crossing angle has a strong impact on the estimation of the bi-tensor parameters and the best performance is obtained for orthogonally crossing fibers.

AC

It is well known that the DT parameters and its derived factors (FA, MD, …), which are usually used in clinical studies, are highly dependent on the acquisition system parameters. Hence, motivated by: • The significant number of tractography-based applications like neurosurgical planning and post-surgery evaluation [46], white-matter Atlas [47] and humain brain connectivity [48], • The use of many DT derived factors to diagnose neurodegenerative diseases like multiple sclerosis [2], Alzheimer [3], or tumors [4], • The limitation of the DT model in bundle-crossing areas which represent about 90% of the white matter voxels [35], • The complexity and challenge of the MTM parameter estimation problem, we propose to use the Cramér-Rao Bound as a tool to analyze the influence of the tunable system parameters on the precision estimation of MTM parameters and its derived factors regardless of the method used for the estimation. We have mainly focused on the SoS reconstruction method throughout the computation of the CRB for the MTM and its derived indices such as the MT-FA and the dominant orientations used in tractography. The same approach can also be used with the other reconstruction methods and the other multi-compartment models.

ACCEPTED MANUSCRIPT

ED

MA

NU

SC

RI P

T

One of the major contribution of this study is the derivation of an analytical FIM formula (section 2.7). This is a major improvement over existing works using the CRB for dMRI applications, particularly [12], [13], [14], [15]. Despite that our analytical expression is an elegant way to solve the computational problem and to speed up the computations, it allows to easily study the performance bounds. In this article, we propose to extend the FA and MD based on the DT model to the MTM as given by Eqs. ((8)) and ((9)). The standard definition based on the DT model may lead to wrong values. Voxels with two or more fiber-bundles will have a lower FA than the single fiber-bundle, which depends on the crossing angles. The proposed fractional anisotropy appears to be a useful tool in crossing-fibers region since it preserves the high anisotropic values. The presented framework allows to study the influence of the number of diffusion gradients and the number of diffusion gradients per shell. We observed fast initial decrease of the error function followed by a saturation effect, which means that there is a limit in the precision when increasing the number of diffusion gradients. To cross this boundary, we must look to another parameter. Based on the CRB, we found that the estimation error variance decreases with a quadratic rate with the SNR value, and is reduced by a factor 𝐶𝐿2 (as compared to the single coil system) when considering an acquisition system using 𝐿 coils. So, in order to reduce the scan time while preserving a good estimation precision, it is preferable to increase the number of acquisition coils rather than the number of diffusion gradient directions. Another point analyzed in this work is the impact of the b-value: the optimal b-value to be used outside the fiber crossing areas, can be approximated by Eq. ((38)). The solution will be in the 3 3 3 3 range 2𝑇𝑟(𝐷) ≤ 𝑏𝑜𝑝𝑡 ≤ 𝑇𝑟(𝐷). For high SNR, 𝑏𝑜𝑝𝑡 ≈ 𝑇𝑟(𝐷) and for low SNR 𝑏𝑜𝑝𝑡 ≈ 2𝑇𝑟(𝐷). In a

CE

PT

crossing-fibers region, multi-shell acquisitions are required [11] where a minimum of two b-values (𝑏1 , 𝑏2 ) are needed. From the study carried out in subsection 3.3, the optimum pair of (𝑏1 , 𝑏2 ) is not unique but depends on the target parameter, the context, and eventually the target cost function.

5 Conclusion

AC

This paper derives the Cramèr-Rao Bound (CRB) for the multi-tensor model associated with multi-coil, multi-shell dMRI data. Besides the CRB of the model parameter vector, we derived the CRBs of other related parameters that are of clinical interest. After providing the exact theoretical expression, we derived an approximate closed-form formula for the FIM. The simplicity of the latter allowed us to conduct a thorough analysis of the impact of different parameters on the estimation performance. Based on this, several key observations have been made on the optimal selection of some controllable system parameters. In particular, the use of these results should lead to significant reduction of the dMRI data acquisition time without affecting the estimation performance.

ACCEPTED MANUSCRIPT

Appendices

T

6 Derivatives of Eigenvalues and Eigenvectors

NU

6.1 Eigenvectors derivatives

SC

RI P

In this appendix, we present mathematical details for the computation of the eigenvectors and the eigenvalues derivatives using the first-order perturbation theory. We consider the eigendecomposition of matrix 𝐷 = 𝐄𝚲𝐄ú where 𝐸 = [𝑒1 , 𝑒2 , 𝑒3 ] is an orthonormal matrix containing the eigenvectors of 𝐷 and Λ is a diagonal matrix whose entries are the ordered eigenvalues (in descending order) of 𝐷, i.e. 𝚲 = diag(𝜆1 , 𝜆2 , 𝜆3 ). We define as in [49]: 𝚷𝑖 = 𝐞𝑖 𝐞ú𝑖 , the orthogonal projector onto the eigensubspace associated with the eigenvalue 𝜆𝑖 and 𝚪𝑖 = (𝐃 − 𝜆𝑖 𝐈)# , where (#) denotes the pseudo-inverse operator.

CE

PT

ED

MA

The first-order perturbation of 𝚷𝑖 is (for details see appendix I of [49]): 𝛿𝚷𝑖 = −𝚷𝑖 𝛿𝐃𝚪𝑖 − 𝚪𝑖 𝛿𝐃𝚷𝑖 , (40) where 𝛿(⋅) refers to the first order perturbation. By right multiplying the terms in Eq. ((40)) by 𝑒𝑖 and using the fact that, at the first order, 𝛿𝑒𝑖ú 𝑒𝑖 = 0 (this is due to the unit norm constraint of 𝑒𝑖 ), we obtain: 𝛿𝑒𝑖 = −𝚪𝑖 𝛿𝐃 𝐞𝑖 . (41) ú Using the relation 𝑣𝑒𝑐(𝐴𝐵𝐶) = (𝐶 ⊗ 𝐴)𝑣𝑒𝑐(𝐵), where 𝑣𝑒𝑐(𝑥) denotes the vectorization of the matrix 𝑥. Then, Eq. ((41)) leads to: 𝛿𝑒𝑖 = −(𝐞ú𝑖 ⊗ 𝚪𝑖 )𝑣𝑒𝑐(𝛿𝐃), (42) which can be written as follows 𝛿𝑒𝑖 = −(𝐞ú𝑖 ⊗ 𝚪𝑖 )𝐽𝛿𝑑. (43) The symbol ⊗ being the Kronecker product and 𝐽 is a 9 × 6 selection matrix defined by 𝑣𝑒𝑐(𝐷) = 𝐽𝑑.

6.2 Eigenvalues derivatives

AC

For the 𝑖 𝑡ℎ eigenvalue, we can write 𝜆𝑖 = Tr(𝐷 𝚷𝑖 ). (44) and hence: 𝛿𝜆𝑖 = Tr(𝛿𝐃 𝚷𝑖 ) + Tr(𝐃 𝛿𝚷𝑖 ). (45) Using the fact that Π𝑖 Γ𝑖 = 0 and consequently, at the first order, Tr(𝐃 𝛿𝚷𝑖 ) = 0. Then, Eq. ((45)) leads to: 𝛿𝜆𝑖 = Tr(𝛿𝐃 𝚷𝑖 ). (46) ú Using the property 𝑣𝑒𝑐(𝐴𝐵𝐶) = (𝐶 ⊗ 𝐴)𝑣𝑒𝑐(𝐵), we obtain (47) Tr(𝛿𝐃 𝚷𝑖 ) = (𝐞ú𝑖 ⊗ 𝐞ú𝑖 )𝑣𝑒𝑐(𝛿𝐃). Finally: 𝛿𝜆𝑖 = (𝐞ú𝑖 ⊗ 𝐞ú𝑖 )𝐽𝛿𝐝. (48)

7 Analytical Approximation Our aim in this appendix is to propose an analytical approximation for the expectation given by Eq. ((13)) to avoid its evaluation by the integral formula given by Eq. ((14)) that is cumbersome and prone to numerical errors. For convenience, we drop the indexes 𝑘 and 𝑗 in

ACCEPTED MANUSCRIPT ((13)), so ℳ = 𝔼[

𝑥 2 𝐼𝐿2 (𝑥) ] − 𝜂2 . 2 𝜂2 𝐼𝐿−1 (𝑥)

T

For small x values:

In that case we use the infinite series of the modified Bessel functions given in [50] 𝑠!(𝑠+𝐿)!

RI P

𝑥 2𝑠+𝐿

1

𝐼𝐿 (𝑥) = ∑∞ 𝑠=0

(2 )

(49)

which the first terms correspond to 𝑥2

𝑥4

(5𝐿+6) 𝑥 6

SC

which allows us to find the MacLaurin expansion of the fraction 𝑓(𝑥, 𝐿) = (𝐼 (7𝐿+12) 𝑥 8

𝐼𝐿 (𝑥) 𝐿−1

2

) for (𝑥)

(21𝐿3 +118𝐿2 +214𝐿+120) 𝑥 10

NU

𝑓𝐿𝑜𝑤 (𝑥, 𝐿) ≈ 4𝐿2 − 8𝐿3 (𝐿+1) + 64𝐿4 (𝐿+1)2 (𝐿+2) − 128𝐿5 (𝐿+1)2 (𝐿+2)(𝐿+3) + 512𝐿6 (𝐿+1)3 (𝐿+2)2 (𝐿+3)(𝐿+4) + ⋯ (50)

MA

Now, combining Eq. ((50)) with Eq. ((13)) and by applying the general expression of the moments of the Nc-Chi distribution [37] 𝑚

1 𝐹1 (. )

where

2 1

2

2

Γ(𝐿)

ED

𝔼(𝑥𝐿𝑚 ) =

𝑚 𝑚 𝜂2 2 2 𝜂 𝑚 Γ(𝐿+ ) 𝐹1 (− ,𝐿,− )

(51)

represents the confluent hypergeometric function of the first kind [51], we obtain 𝜂2

ℳLow =

𝐿

𝜂4

(2𝐿+3) 𝜂 6

− 𝐿2 +

2𝐿3 (𝐿+1)

−

(𝐿+3) 𝜂 8 𝐿4 (𝐿+1)

+⋯

(52)

PT

We note here that the moments turn out to be simple polynomials. We observe that as 𝐿 grows one can simplify Eq. ((52)) to 𝜂2

𝑖

𝜂2 𝐿

).

AC

term (−

CE

ℳLow ≈ − ∑𝑞𝑖=1 (− 𝐿 ) + 𝑂(𝜂2 )𝑞+1 (53) which represents the sum of the first 𝑞 terms of a truncated geometric series with as common

ℳLow ≈

𝑞 𝜂2 ) 𝐿 𝐿 1+ 2 𝜂

1−(−

(54)

The latter is a convergent series at low SNR, i.e. for 𝐿

−1

𝜂2 𝐿

< 1 ⇔ 𝜂 < √𝐿 and leads to the limit

value ℳ ≈ (1 + 𝜂2 ) .

For large x values:

1

We use here the following asymptotic expansion in powers of (𝑥) of the modified Bessel function of the first kind [51] when 𝑥 → ∞ and 𝐿 held fixed: 𝐼𝐿 (𝑥)~

𝑒𝑥

√2𝜋𝑥

(1 −

(4𝐿2 −1) 8𝑥

+

(4𝐿2 −1)(4𝐿2 −9) 2!(8𝑥)2

−

(4𝐿2 −1)(4𝐿2 −9)(4𝐿2 −25) 3!(8𝑥)3

+. . . )

This expansion allows us to calculate an asymptotic expression of 𝑓(𝑥, 𝐿) = (𝐼 the form

(55) 𝐼𝐿 (𝑥) 𝐿−1

2

) in (𝑥)

ACCEPTED MANUSCRIPT 𝑓𝐻𝑖𝑔ℎ (𝑥, 𝐿)

≈1−

(2𝐿−1)

+

(2𝐿−1)(2𝐿−3) 𝑥2

−

(2𝐿−1)(2𝐿−3)2

+

(𝐿−2)(2𝐿−3)(2𝐿−1)

8 𝑥3 4 𝑥4 (𝐿−3)(2𝐿−3)(2𝐿−1)(−9+4𝐿(𝐿−2))

𝑥 (2𝐿−7)(2𝐿−5)(2𝐿−3)(2𝐿−1)(2𝐿+3)

(56) + − + ⋯ 128 𝑥 5 16 𝑥 6 Plugging Eq. ((56)) in Eq. ((13)) and using the general expectation expression Eq. ((51)) leads to ≈ 2𝐿 − −

𝜂2

Γ(𝐿) 𝜂

𝜂2 (𝐿−2)(2𝐿−3)(2𝐿−1) 1 𝐹1 (1;𝐿;− ) 2 8(𝐿−1) 𝜂 4

+

(2𝐿−1)(𝐿−1) 𝜂2

−

+⋯

1 2

𝜂2 ) 2

8√2 Γ(𝐿) 𝜂 3

SC

(57)

1 2

(2𝐿−1)(2𝐿−3)2 Γ(𝐿− ) 1 𝐹1 ( ;𝐿;−

T

ℳHigh

1

RI P

1

√2(2𝐿−1) Γ(𝐿+2) 1 𝐹1 (−2;𝐿;− 2 )

NU

By applying the asymptotic expansion of the confluent hypergeometric function of the first kind for large SNR to Eq. ((57)), one obtains (2𝐿−1) (2𝐿−3)(2𝐿−1) (𝐿−3)(2𝐿−3)(2𝐿−1) (2𝐿−3)(2𝐿−1)(4𝐿(𝐿−8)+57) ℳHigh ≈ 1 − 2𝜂2 + − + + ⋯ (58) 4 6 4𝜂 4𝜂 16𝜂 8 The latter for larger values of 𝐿 reduces to the following simplified expression 𝐿

𝑖

1 𝑞

MA

ℳHigh ≈ ∑𝑞−1 𝑖=0 (− 𝜂 2 ) + 𝑂 (𝜂 2 )

(59) 𝐿

which represents the sum of the first 𝑞 terms of a geometric series of common terms (− 𝜂2 ) 𝐿

ED

ℳHigh ≈

𝑞

1−(− 2 ) 𝜂 𝐿

1+ 2 𝜂

(60)

𝐿

ℳ ≈ (1 + 𝜂2 )

−1

PT

This sum converges for high SNR values (i.e. for 𝜂 > √𝐿) and leads to the limit value (The same limit as for Eq. ((54))).

AC

CE

In Fig. 9, we compare the approximate integral values (i.e. ℳHigh and ℳLow ) to the exact one, evaluated numerically. These plots illustrate the effectiveness of the approximations for different values of 𝐿. In order to quantify the approximation error between the integral formula given by Eq. ((14)) and our analytical one given by Eq. ((32)), we plot in Fig. 10 the relative error as a function of the SNR (𝜂) for different 𝐿 values. We note from this figure that the relative error is relatively small and decreases significantly when 𝐿 or the SNR increases.

References [1] Basser PJ, Mattiello J, LeBihan D, MR diffusion tensor spectroscopy and imaging, Biophys J 66 (1) (1994) 259–67.

[2] Emilia Sbardella, Francesca Tona, Nikolaos Petsas, and Patrizia Pantano, DTI Measurements in Multiple Sclerosis: Evaluation of Brain Damage and Clinical Implications , Multiple Sclerosis International 2013:671730.

ACCEPTED MANUSCRIPT

T

[3] Acosta-Cabronero J, Alley S, Williams GB, Pengas G, Nestor PJ, Diffusion tensor metrics as biomarkers in Alzheimer’s disease, PLoS One (2012) 7(11) e49072:1–14.

RI P

[4] Schmainda KM, Diffusion-weighted MRI as a biomarker for treatment response in glioma, CNS oncology (2012) 1(2):169–180.

NU

SC

[5] Assemlal HE, Tschumperl̩ D, Brun L, Siddiqi K, Recent advances in diffusion MRI modeling: Angular and radial reconstruction, Medical Image Analysis 15 (04) (2011) 369Р96.

MA

[6] Panagiotaki E, Schneider T, Siow B, Hall MG, Lythgoe MF, Alexander DC, Compartment models of the diffusion MR signal in brain white matter: a taxonomy and comparison, Neuroimage 59 (3) (2012) 2241–54.

PT

ED

[7] Megherbi T, Kachouane M, Oulebsir-Boumghar F, and Deriche R, Crossing Fibers Detection with an Analytical High Order Tensor Decomposition, Computational and Mathematical Methods in Medicine 2014 (476837) (2014) 18.

CE

[8] Tuch DS, Diffusion MRI of Complex Tissue Structure, Ph.D. thesis, MIT (2002).

AC

[9] Kreher BW, Schneider JF, Mader I, Martin E, Hennig J, Il’yasov KA, Multitensor approach for analysis and tracking of complex fiber configurations, Magn Reson Med (2005) 54(5):1216–25.

[10] Hosey T, Williams G, Ansorge R, Inference of multiple fiber orientations in high angular resolution diffusion imaging, Magn Reson Med (2005) 54(6):1480–9.

[11] Scherrer B, Warfield SK, Parametric representation of multiple white matter fascicles from cube and sphere diffusion MRI, PLoS One (2012) 7(11):e48232.

[12] Brihuega-Moreno O, Heese FP, Hall LD, Optimization of diffusion measurements using Cramer-Rao lower bound theory and its application to articular cartilage, Magn Reson Med (2003) 50(5):1069–76.

[13] Alexander DC, A general framework for experiment design in diffusion MRI and its

ACCEPTED MANUSCRIPT application in measuring direct tissue-microstructure features, Magn Reson Med (2008) 60(2):439–48.

RI P

T

[14] Caan MW, Khedoe HG, Poot DH, den Dekker AJ, Olabarriaga SD, Grimbergen KA, van Vliet LJ, Vos FM, Estimation of diffusion properties in crossing fiber bundles, IEEE Trans Med Imaging (2010) 29(8):1504–15.

SC

[15] Beltrachini L, von Ellenrieder N, Muravchik CH, Error bounds in diffusion tensor estimation using multiple-coil acquisition systems, Magn Reson Imaging (2013) 31(8):1372–83.

MA

NU

[16] Basser PJ, Pierpaoli C, Microstructural and physiological features of tissues elucidated by quantitative-diffusion-tensor MRI, J Magn Reson B (1996) 111(3):209–19.

ED

[17] Roemer PB, Edelstein WA, Hayes CE, Souza SP, Mueller OM, The NMR phased array, Magn Reson Med 16 (2) (1990) 192–225.

PT

[18] Andersen AH, Kirsch JE, Analysis of noise in phase contrast MR imaging, Medical Physics (1996) 23(6):857–69.

CE

[19] Constantinides CD, Atalar E, McVeigh ER, Signal-to-noise measurements in magnitude images from NMR phased arrays, Magn Reson Med (1997) 38(5):852–7.

AC

[20] Wilkinson J. H, Probability distributions involving Gaussian random variables, Kluwer: Norwell, MA, 2002.

[21] Hayes CE, Roemer PB, Noise correlations in data simultaneously acquired from multiple surface coil arrays, Magn Reson Med (1990) 16(2):181–91.

[22] Brown R, Wang Y, Spincemaille P, Lee RF, On the noise correlation matrix for multiple radio frequency coils, Magn Reson Med (2007) 58(2):218–24.

[23] Aja-Fernández S, Brion V, Tristán-Vega A, Statistical noise analysis in GRAPPA using a parametrized noncentral Chi approximation model, Magn Reson Med (2011) 65(4):1195– 206.

ACCEPTED MANUSCRIPT

T

[24] Aja-Fernández S, Brion V, Tristán-Vega A, Influence of noise correlation in multiple-coil statistical models with sum of squares reconstruction, Magn Reson Med (2012) 67(2):580–5.

RI P

[25] Aja-Fernández S, Vegas-Sánchez-Ferrero G, Tristán-Vega A, Noise estimation in parallel MRI: GRAPPA and SENSE, Magn Reson Imaging (2014) 32(3):281–90.

SC

[26] Pruessmann KP, Weiger M, Scheidegger MB, Boesiger P., SENSE: sensitivity encoding for fast MRI., Magn Reson Med. (1990) 42(5):952–62.

MA

NU

[27] Griswold MA, Jakob PM, Heidemann RM, Nittka M, Jellus V, Wang J, Kiefer B, Haase A, Generalized autocalibrating partially parallel acquisitions (GRAPPA), Magn Reson Med (2002) 47(6):1202–10.

ED

[28] Tanuj Kumar Jhamb, Vinith Rejathalal, V.K. Govindan, A Review on Image Reconstruction through MRI k-Space Data, IJIGSP 7 (7) (2015) 42–59.

PT

[29] Hollingsworth KG, Reducing acquisition time in clinical MRI by data undersampling and compressed sensing reconstruction, Phys Med Biol (2015) 60(21):R297–322.

AC

CE

[30] Stejskal, E. O. and Tanner, J. E., Spin Diffusion Measurements: Spin Echoes in the Presence of a Time―Dependent Field Gradient, The Journal of Chemical Physics 42 (1) (1965) 288–292.

[31] Xu D, Cui J, Bansal R, Hao X, Liu J, Chen W, Peterson BS, The ellipsoidal area ratio: an alternative anisotropy index for DTI, Magn Reson Imaging (2009) 27(3):311–23.

[32] Ennis DB, Kindlmann G, Orthogonal tensor invariants and the analysis of diffusion tensor magnetic resonance images, Magn Reson Med (2006) 55(1):136–46. [33] Goveas J, O’Dwyer L, Mascalchi M, Cosottini M, Diciotti S, De Santis S, Passamonti L, Tessa C, Toschi N, Giannelli M, Diffusion-MRI in neurodegenerative disorders, Magn Reson Imaging 33 (7) (2015) 853–76.

[34] Douaud G, Jbabdi S, Behrens TE, Menke RA, Gass A, Monsch AU, Rao A, Whitcher B, Kindlmann G, Matthews PM, Smith S, DTI measures in crossing-fibre areas: increased

ACCEPTED MANUSCRIPT diffusion anisotropy reveals early white matter alteration in MCI and mild Alzheimer’s disease, Neuroimage (2011) 55(3):880–90.

RI P

T

[35] Jeurissen B, Leemans A, Tournier JD, Jones DK, Sijbers J, Investigating the prevalence of complex fiber configurations in white matter tissue with diffusion magnetic resonance imaging, Hum Brain Mapp 34 (11) (2013) 2747–66.

SC

[36] Steven M. Kay, Fundamentals of statistical signal processing: estimation theory, Prentice-Hall: USA, 1993.

MA

NU

[37] Sijbers J, Signal and noise estimation from magnetic resonance images, Ph.D. thesis, Universitaire Instelling Antwerpen (Belgium) (1998).

ED

[38] Jones DK, Basser PJ, Squashing peanuts and smashing pumpkins: How noise distorts diffusion-weighted MR data, Magn Reson Med (2004) 52(5):979–93.

PT

[39] Descoteaux M, Deriche R, Le Bihan D, Mangin JF, Poupon C, Multiple q-shell diffusion propagator imaging, Med Image Anal (2011) 15(4):603–21.

CE

[40] Alexander DC, Barker GJ, Optimal imaging parameters for fiber-orientation estimation in diffusion MRI, Neuroimage 27 (2) (2005) 357–67.

AC

[41] Kingsley PB, Monahan WG, Selection of the optimum b factor for diffusion-weighted magnetic resonance imaging assessment of ischemic stroke, Magn Reson Med 51 (5) (2004) 996– 1001.

[42] Bito Y, Hirata S, Yamamoto E., Optimal gradient factors for ADC measurements, in: 3rd Annual Meeting of International Society for Magnetic in Medicine, 1995, p. 913.

[43] Jones DK, Horsfield MA, Simmons A, Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging, Magn Reson Med 42 (1999) 515–525.

[44] Pierpaoli C, Jezzard P, Basser PJ, Barnett A, Di Chiro G, Diffusion tensor MR imaging of the human brain, Radiology (1996) 201(3):637–48.

ACCEPTED MANUSCRIPT [45] Caruyer E, Lenglet Ch, Sapiro G, Deriche R, Design of multishell sampling schemes with uniform coverage in diffusion MRI, Magn Reson Med 69 (6) (2013) 1534–1540.

RI P

T

[46] Golby AJ, Kindlmann G, Norton I, Yarmarkovich A, Pieper S, Kikinis R, Interactive Diffusion Tensor Tractography Visualization for Neurosurgical Planning, Neurosurgery (2011) 68(2):496–505.

SC

[47] Catani M, Thiebaut de Schotten M, A diffusion tensor imaging tractography atlas for virtual in vivo dissections, Cortex (2008) 44(8):1105–32.

MA

NU

[48] Lazar Mariana, Mapping brain anatomical connectivity using white matter tractography, NMR Biomed (2010) 23(7):821–35.

ED

[49] Abed-Meraim K, Loubaton P, Moulines E, A subspace algorithm for certain blind identification problems, IEEE Transactions on Information Theory 43 (2) (1997) 499–511.

PT

[50] George B. Arfken, Hans J. Weber and Frank E. Harris, Chapter 14 - Bessel Functions , in: Mathematical Methods for Physicists (Seventh Edition) , Academic Press, Boston, 2013, pp. 643 – 713.

AC

CE

[51] Abramowitz M, Stegun IA, Handbook of mathematical functions, New York: Dover Publications, 1970.

ACCEPTED MANUSCRIPT

Figures

𝐿

NU

SC

RI P

T

Figure 1: Mono-shell (left) versus Multi-shell (right). The dots represent the encoded diffusion gradient directions: On the left, single-shell, j=1, dMRI data are acquired with one b-value and one set of diffusion gradients (e.g. 1000𝑠. 𝑚𝑚−2 × 120𝑑𝑖𝑟𝑒𝑐𝑡𝑖𝑜𝑛𝑠). On the right, double shell, j=1,2, two b-values associated with two set of diffusion gradient directions (e.g. 1000𝑠. 𝑚𝑚−2 × 24𝑑𝑖𝑟𝑒𝑐𝑡𝑖𝑜𝑛𝑠 and 3000𝑠. 𝑚𝑚−2 × 96𝑑𝑖𝑟𝑒𝑐𝑡𝑖𝑜𝑛𝑠) [45]. Figure 2: Fractional anisotropy in crossing fibers. (a) Diffusion tensor FA: high in voxels with single fiber-bundle, the tensor has a prolate form. In regions of crossing fibers, FA is reduced. With two fiber-bundles, the tensor is more oblate in crossing at 90∘ than at 60∘ . Thus, FA varies w.r.t. the crossing angle. (b) MT-FA: In crossing area, the MT-FA is high and is a weighted average of the crossing bundles FA. Figure 3: Graphical comparison of the obtained approximation and the integral formula. We −1

plot the final results ℳ ≈ (1 + 𝜂2 )

in dashed blue line and the integral formula in a continued

AC

CE

PT

ED

MA

green line, as a function of the SNR (𝜂) for different L values. Figure 4: The impact of the number of diffusion gradients per shell: The relative error measure of MT-FA when (a) The two tensors are highly anisotropic, (b) The second tensor is quasi isotropic. In (c) we show the relative error measure of e 1 when the two tensors are highly anisotropic. Figure 5: The relative error measure of (a) MT-FA, and (b) e 1 vs the number of diffusion gradients for different acquisition coil numbers. Figure 6: The relative error measure of (a) MT-FA, and (b) e 1 w.r.t. (𝑏1 , 𝑏2 ). The two tensors are highly anisotropic. Figure 7: The relative error measure of (a) MT-FA, and (b) e 1 w.r.t. (𝑏1 , 𝑏2 ). The first tensor is highly anisotropic while the second is “quasi-isotropic". Figure 8: The impact of the crossing angle on the estimation of the MT-FA and e 1 . Figure 9: The graph of the integral formula, evaluated numerically, together with our two proposed expressions as a function of the SNR (𝜂) for different L values. Figure 10: Graph of the relative approximation error between the integral formula and our analytical one.

MA

NU

SC

RI P

T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

Figure 1

NU

SC

RI P

T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

MA

Figure 2

ED

MA

NU

SC

RI P

T

ACCEPTED MANUSCRIPT

AC

CE

PT

Figure 3

AC

Figure 4

CE

PT

ED

MA

NU

SC

RI P

T

ACCEPTED MANUSCRIPT

NU

SC

RI P

T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

MA

Figure 5

NU

SC

RI P

T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

MA

Figure 6

NU

SC

RI P

T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

MA

Figure 7

SC

RI P

T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

MA

NU

Figure 8

ED

MA

NU

SC

RI P

T

ACCEPTED MANUSCRIPT

AC

CE

PT

Figure 9

ED

MA

NU

SC

RI P

T

ACCEPTED MANUSCRIPT

AC

CE

PT

Figure 10