Asymptotic Behavior of the Direction-Averaged Diffusion-WeightedMRI Signal using Different B-Tensor Encoding Schemes

1Cardiff University Brain Research Imaging Centre (CUBRIC), School of Psychology, Cardiff University, Cardiff, United Kingdom 2Laboratorio de Procesado de Imagen, ETSI Telecomunicación Edificio de las Nuevas Tecnologías, CampusMiguel Delibes s/n, Universidad de Valladolid, 47011 Valladolid, Spain 3MaryMacKillop Institute for Health Research, Faculty of Health Sciences, Australian Catholic University, Melbourne, Victoria, 3065, Australia


| INTRODUCTION
Pathological disorders happen at the cellular level, therefore, it is important to obtain information on the micrometer scale. Diffusion MRI provides a tool to study brain microstructure based on the Brownian motion of water molecules [1] and it is therefore sensitive to the changes in the microstructure of the tissue [2, 3, 4]. Different mathematical representations are proposed to describe this relationship between the diffusion signal and the changes in the microstructure [5, 6,7] the most prominent are the bi-exponential [8,9,10,11,12], the stretched exponential [13] and the power law [14,15,16,17] models. The mathematical forms of these approaches are quite different. In the bi-exponential approach, the large b-value behavior is assumed to be dominated by the intracellular compartment. For stretched exponentials, the signal relationship with the b-value is exp [−(k b) a ], where k is a constant and a < 1 is the stretching parameter. In the statistical model developed by Yablonskiy et al. [14], the signal decays as 1/b for large b, while the other studies [15,16,17], have reported that the signal at high b-value, decays as 1/ √ b.
The aforementioned studies all used single diffusion encoding (SDE). Since the development of the Pulsed Gradient Spin Echo (PGSE) sequence [18], there have been many works aimed at maximizing the information that can be obtained from a dMRI experiment by exploring different acquisition protocols [19,20]. One such modification is the addition of multiple gradient pairs. We can use two pairs of pulsed-field gradients (PFGs) to obtain a Double Diffusion Encoding (DDE) [21,22]. It has been shown that DDE, as well as other multiple encoding schemes such as triple diffusion encoding (TDE) [23], provide information that is not accessible with SDE [24].
This approach has been utilized by several groups for extracting microstructural information [25,26,27,28]. A framework was recently proposed [29] to probe tissue using different q-space trajectory encodings which can be described by a second-order b-tensor. SDE, DDE, and TDE can be characterized by b-tensors, with one, two and three non-zero eigenvalues, respectively. In this framework, SDE is also called Linear tensor encoding (LTE), DDE with perpendicular directions is called planar tensor encoding (PTE) and TDE with equal eigenvalues is called spherical tensor encoding (STE).
In this study, we investigate the effect of different b-tensor encodings on the diffusion signal at high b-values. To remove the effect of fiber orientation distribution [30], the acquired signal is averaged over all diffusion directions for each shell. This so-called powder-averaged signal [31,32] has less complexity than the direction-dependent signal.
Powder averaging yields a signal whose orientation-invariant aspects of diffusion are preserved but with an orientation distribution that mimics complete dispersion of anisotropic structures.
For the powder-averaged signal, the diffusion attenuation is a function of the orientation-invariant aspects of the diffusion and the encoding. The compartment diffusion attenuation is (Eq. (2) in [33]): where S is the normalized diffusion signal and D || and D ⊥ are the the parallel and perpendicular diffusivities respectively.
We use the subscript " e " and " a " to denote the parameters of the extra-and intra-axonal compartments respectively.
Here, we study the effect of b-tensor shape on the diffusion-weighted signal at high b-values. We begin with linear, planar and spherical before considering general b-tensor shapes.

| LTE
In linear tensor encoding, b ∆ = 1 and assuming stick-like geometry, D ⊥ = 0 in Eq. (1), therefore S i c has the following form: Because bD a || 1 for large b-values (b > 7000 s/mm 2 ), erf( bD a || ) = 1 and the power-law scaling [17] for the intraaxonal compartment becomes: The sensitivity of MR to axon radius would alter the b −1/2 scaling [34] because there will be a perpendicular diffusivity and the exponential term in Eq.
(1) will not be zero. For large b-values, the extra-axonal signal (see Appendix A) decays exponentially faster than the intra-axonal compartment (exp (−bD e ⊥ ) 1) and therefore, can be neglected.

| PTE
In planar tensor encoding, b ∆ = −1/2 and S i c has the following form: For large b-values, bD a || 1, therefore the diffusion signal can be approximated by the following equation (see Appendix A): where !! denotes the double factorial and N depends on the bD a || value ( Fig. 1 and Table 1).
For large b-values, the extra-axonal signal decays exponentially faster than the intra-axonal compartment, exp (−bD e ⊥ ) 1, (see Appendix A) and can be neglected.

| STE
In spherical tensor encoding, b ∆ = 0 and S i c has the following form: For large b-values, both intra-and extra-axonal signals decay exponentially fast, exp 1 and both of them are negligible (see Appendix A). Therefore, the spherical tensor encoding does not provide a considerable signal for large b-values in a two-compartment model.

| General Case
Here, we consider the general case of b ∆ 0, to cover all b-tensor shapes between b ∆ = −0.5 (PTE) to b ∆ = 1 (LTE).

2.2.1
As noted above, in this range the error function in Eq.
which is only physically plausible (i.e. the ratio of diffusion coefficients has to be ≥ 0) for b ∆ − 1 ≥ 0, but the maximum value that b ∆ can take is one, and therefore D ⊥ has to be zero to satisfy Eq. (7) i.e. the geometry has to be that of a stick, and the b-tensor has to be a pure LTE to have a power-law relationship.

2.2.2
Conversely, in the range −0.5 ≤ b ∆ < 0, as in Eq. (4), the error function becomes imaginary. Similar to the first scenario, to have a power-law relationship the exponential term has to be one. By replacing the first term of the approximation in Eq. (16) into Eq.
(1), we have: To have the exponential equal to one: where the right side of the equation is negative for −0.5 < b ∆ < 0 which is not physically plausible for the left side of the equation (i.e. ratio of diffusivities). Therefore, the only possible case is to have D ⊥ = 0 which again means stick-like geometry and b ∆ = −0.5 which is pure PTE. Clearly the exponential term will become zero if and only if b ∆ = −0.5, and thus the 1/b signal form will occur if and only if the b-tensor shape has just 2 non-zero eigenvalues, i.e. pure PTE. Thus, for stick-like geometries, there are only two b-tensor shapes for which a power-law exists: pure linear and pure planar.
As the above equations show, we do not observe a power-law for non-stick-like geometries.

| Simulations
Synthetic data were generated with 30 diffusion encoding gradient orientations uniformly distributed on the unit sphere [35,36] and 21 b-values spaced in the interval [0, 10000 s/mm 2 ] with a step-size of 500 s/mm 2 . The noise is considered Rician with SNR = 250 for the b0 image, which is practically feasible using the Connectom scanner [37] with an echo time of 88ms [38]. A two-compartment model with a Watson orientation distribution function is used: where W (n) is the Watson ODF and S c y l is the signal attenuation of the impermeable cylinders [39] in the presence of b-tensor encoding [40] (Appendix B). The ground truth parameter values defined by a set of parameters [f = 0.65, D a || = 2 µm 2 /ms, D || e = 2 µm 2 /ms, D ⊥ e = 0.25, 0.5, 0.75 µm 2 /ms and κ = 11] and axon radius r i , come from the bins of the histograms in [41]. We average the signal over all r i s weighted by r 2 i . In histology, there is a possibility of tissue shrinkage. To account for this change, the axon radius values are multiplied with three shrinkage factors η = 0, 1, 1.5 [41,42]. The η = 0 case simulates the effect of zero-radius axons. The noisy diffusion signal is modeled according to the following: where S n and S are the noisy and noise-free signal respectively and N r and N i are the normal distributed noise in the real and imaginary images respectively with a standard deviation of σ [43,44].

| In vivo Data
Two healthy participants who showed no evidence of a clinical neurologic condition were scanned in this study. Diffusionweighted images were acquired with 30 gradient directions for planar tensor encoding (PTE) on a 3T Connectom MR imaging system (Siemens Healthineers, Erlangen, Germany). Twenty axial slices with a voxel size of 4mm isotropic (given the strong signal attenuations investigated here, a low resolution of 4 mm isotropic was used) and a 64×64 matrix size, TE = 88 ms, TR = 4300 ms, were obtained for each individual.
Diffusion data were acquired for 10 b-value shells from 1000 to 10000 s/mm 2 with a step size of 1000 and each shell had the same 30 diffusion encoding gradient orientations uniformly distributed on the unit sphere. One b0 image was acquired between each b-value shell as a reference.
The data were corrected for Gibbs ringing [45], eddy current distortions and subject motion [46]. We normalized the direction-averaged signal based on the b0 signal in each voxel.

| RESULTS
The asymptotic expansion of erfi(x ) in Eq. 16 (see Appendix A) is valid when x → ∞, but large values of bD a || would suppress the signal to immeasurable levels, and therefore there are practical bounds on the value of bD a || that can be achieved. Therefore, we compared the original signal in Eq. 4 and the approximated signal using Eq. 5 for different values of N and bD a || ( Fig. 1 and Table 1). We use a normalized error to compare the original (Eq. 4) and the approximated signal (Eq. 5): where S is the original signal obtained from Eq. 4 andŜ is the approximated signal from Eq. 5. Fig. 1 showsŜ /S for 3 < bD a || < 20 and 4 < N < 21. The selected range of bD a || is compatible with the range of b-values that we can obtain from the Connectom scanner and also the range of D a || that exist in the brain [47]. Based on Fig. 1 the number of terms in Eq. 10 should be smaller than or equal to the bD a || (N ≤ bD a || where ... denotes the floor function) to have the minimum error (Ŝ /S is close to one). As the number of terms goes beyond the bD a || , the error increases.    ≥ 14). One of the challenges with the high b-values is the noise, as the signal amplitude can be close to the noise floor. Therefore, here we find the maximum value of bD a || that we can use before hitting this rectified noise floor (see Appendix C).
The noise in complex MR data is normally distributed, whereas the noise in magnitude images is Rician distributed [43,44]. Here, we select a minimum SNR value equal to 2 (see Appendix C). By setting the diffusion-weighted intensity in Eqs.
(2), (4), and (6) to the mean background signal, we obtain the b-value that makes the signal equal to the noise floor. Fig. 2 shows the maximum bD a || as a function of SNR for different encoding schemes and different noise floors. The maximum value of bD a || that can be used while staying above the noise floor increases when SNR increases, but the rate of this change is different for different encoding schemes. The maximum bD a || value (bD a || max ) is proportional to the square of SNR, (bD a || max ∼ S N R 2 ) for LTE, where this relationship is linear for PTE (bD a || ∼ S N R ) and it is logarithmic for STE (bD a || ∼ ln (S N R )). Based on this plot, if SNR = 50 the values of bD a || max for linear, planar and spherical tensor encoding schemes are around 312, 21 and 9 respectively. The SNR in our data is around 250 [38] therefore the measured signal values in our experiment are higher than the noise level. For this SNR, the bD a || max for linear, planar and spherical tensor encoding schemes are around 15625, 100 and 16 respectively.   F I G U R E 3 Simulated direction-averaged PTE signal for 7000 < b < 10000s/mm 2 and the results of the power-law fit. [48]. In the WM, the α value is close to one, supporting the theory. In grey matter and CSF, the exponent is larger (1.16 and 1.26, respectively). According to the theory outlined above, this would be consistent with a lack of pure 'stick-like' geometry in these tissue components. The spatial resolution of the data must be recognized, i.e., at 4 mm isotropic voxels, obtaining a 'pure' GM signal and 'pure' CSF signal is challenging. It is likely that the intermediate exponent in the GM between that of the WM and CSF is partly attributable to a partial volume effect, and partly attributable to the inadequacy of the model for grey matter architecture. Further investigation of this phenomenon in grey matter is beyond the scope of this work.

| DISCUSSION
The main finding of this paper is a theoretical derivation, and confirmation in silico and in vivo of a power-law relationship between the direction-averaged DWI signal and the b-value using planar tensor encoding, as given by Eq. 5, for b-values ranging from 7000 to 10000 s/mm 2 . In white matter, the average value of the estimated exponent is α = 1.02 ± 0.02.
For smaller b-values, this behavior must break down as the DWI signal of PTE cannot be approximated by Eq. 5 (Fig. 1) and also we cannot neglect the contribution of the extracellular compartment. It could also fail for very large b-values, if there were immobile protons that contributed a constant offset to the overall signal or if there is any sensitivity to the axon diameter [34].
The exponent of approximately one for white matter using PTE is consistent with the large b-value limit predicted for a model of water confined to sticks Eq. 5, which is used to describe the diffusion dynamics of intra-axonal water. Our results confirm this relationship between the diffusion signal and the b-value (Fig. 6).
The b −1/2 -scaling has previously been suggested by [17,16] for linear tensor encoding. Two other proposed models predict power law signal decay, for large b-values using a linear tensor encoding. One of these is the statistical model [14], where the signal decays as 1/b for large b. Some other models [49,50,51], assume a gamma distribution for the diffusion coefficients and the signal has the following form: which scales as (b c /b) for b b c . However, in this case, the exponent does not have a universal value, it depends on the distribution.
This work interprets the diffusion-weighted MRI signal decay at high b-values in the form of S ∼ b −1 for planar tensor encoding. An important application of this finding is using the combination of linear and planar tensor encodings to characterize the intra-axonal diffusivity and the signal fraction as it is proposed by [52] using triple diffusion encoding. encoding. The different exponents of these encodings could be used to provide independent validation of the presence of stick-like geometries in vivo where, due to the slower decay, LTE is the most SNR efficient method. Any deviation from the power-law could indicate deviation from stick-like geometry with both LTE and PTE encoding, as exploited by [34] for estimating the effective radius of axons. Again, for such applications of the power-law, the LTE approach is to be favoured over the PTE approach, on account of the higher SNR per unit time.

A.1 | LTE
In linear tensor encoding, b ∆ = 1 and S ec has the following form:
bD a || 1 for large b, therefore we have: where N depends on the bD a || value ( Fig. 1 and table 1).
The extra-axonal signal decays exponentially faster than the intra-axonal compartment for large b, e −bDe ⊥ 1, and can be neglected.

A.3 | STE
In spherical tensor encoding, b ∆ = 0 and S ec has the following form:

B | SIGNAL ATTENUATION IN A CYLINDRICAL PORE USING PTE
The signal attenuation of the impermeable cylinders [39] using PTE is generated using the following equation [40]: ln (S ⊥ c y l ) = − where α n is the root of the derivatives of the first order Bessel function J 1 (α n ) = 0 and and

C | SNR AND ERROR
Let us assume that a real signal S follows a Rician distribution with parameters A and σ with PDF [43] p(x |A, σ) = x σ 2 e − x 2 +A 2 where A is the (absolute value) of the original signal (without noise) and σ 2 is the variance of the complex Gaussian noise. It can be seen as The question in MRI of how low can we go with the signal (i.e., when do we reach the noise floor) will always depend on the application and on the estimator we are using. However, we can always consider a lower bound to the SNR related to the error of the measured signal.
To calculate an SNR threshold independent of the particular application, we can use two different definitions of error: (1) the Mean Square Error (MSE) or (2) the mean error (ME). We define the MSE as where S is the measured signal and A is the original signal. We use the mean value to assure that this error is a statistical property and not an isolated measure. The ME is alternatively define as: For the SNR calculation, we will consider that the error committed is a percentage of the original signal (to make it signal dependent), i.e.
Alternatively, for the ME: Assuming a Rician distribution of parameters A and σ, the errors become: The MSE: The ME: For the sake of simplicity, in this paper, we will consider ME as an error measure, since MSE is more restrictive. The relation between ME and SNR for different errors can be seen in Table 2.

A C K N O W L E D G E M E N T S
The data were acquired at the UK National Markus Nilsson for providing the pulse sequences for b-tensor encoding.

R E F E R E N C E S
[1] Tanner J. Self diffusion of water in frog muscle. Biophysical journal 1979;28(1):107.