Tensor-valued diffusion encoding for diffusional variance decomposition (DIVIDE): Technical feasibility in clinical MRI systems

Microstructure imaging techniques based on tensor-valued diffusion encoding have gained popularity within the MRI research community. Unlike conventional diffusion encoding—applied along a single direction in each shot—tensor-valued encoding employs diffusion encoding along multiple directions within a single preparation of the signal. The benefit is that such encoding may probe tissue features that are not accessible by conventional encoding. For example, diffusional variance decomposition (DIVIDE) takes advantage of tensor-valued encoding to probe microscopic diffusion anisotropy independent of orientation coherence. The drawback is that tensor-valued encoding generally requires gradient waveforms that are more demanding on hardware; it has therefore been used primarily in MRI systems with relatively high performance. The purpose of this work was to explore tensor-valued diffusion encoding on clinical MRI systems with varying performance to test its technical feasibility within the context of DIVIDE. We performed whole-brain imaging with linear and spherical b-tensor encoding at field strengths between 1.5 and 7 T, and at maximal gradient amplitudes between 45 and 80 mT/m. Asymmetric gradient waveforms were optimized numerically to yield b-values up to 2 ms/μm2. Technical feasibility was assessed in terms of the repeatability, SNR, and quality of DIVIDE parameter maps. Variable system performance resulted in echo times between 83 to 115 ms and total acquisition times of 6 to 9 minutes when using 80 signal samples and resolution 2×2×4 mm3. As expected, the repeatability, signal-to-noise ratio and parameter map quality depended on hardware performance. We conclude that tensor-valued encoding is feasible for a wide range of MRI systems—even at 1.5 T with maximal gradient waveform amplitudes of 33 mT/m—and baseline experimental design and quality parameters for all included configurations. This demonstrates that tissue features, beyond those accessible by conventional diffusion encoding, can be explored on a wide range of MRI systems.


Introduction
Diffusion magnetic resonance imaging (dMRI) enables non-invasive imaging of tissue microstructure. The vast majority of dMRI techniques rely on the conventional Stejskal-Tanner experiment [1] that employs a pair of pulsed field gradients to yield diffusion encoding along a single direction for each preparation of the signal, so called 'linear' encoding. Methods such as diffusion tensor imaging (DTI) [2] yield voxel-scale average parameters which are sensitive to alterations of the tissue microstructure in both healthy development and disease [3]. Diffusional kurtosis imaging (DKI) [4] is an extension to DTI that can probe non-gaussian processes and sub-voxel heterogeneity-a feature that has the potential to differentiate, for example, high and low-grade gliomas [5,6]. Information on the microstructure can also be inferred from biophysical models, which may provide more specific information based on assumptions about the tissue composition [7]. However, techniques that rely on linear diffusion encoding alone all have a fundamental drawback: they convolve the effects of microscopic diffusion anisotropy, orientation dispersion, and heterogeneous isotropic diffusivity [8][9][10][11]. Consequently, markedly different tissue archetypes, such as elongated cell structures that are randomly oriented and isotropic tissues with varying cell density, may be indistinguishable regardless of the modelling approach [8,9].
This limitation can be mitigated by performing diffusion encoding along multiple directions in a single acquisition. This approach was introduced by Cory et al. [12], who proposed parallel and orthogonal double diffusion encoding (DDE) as a probe of local pore geometry. This concept has since been elaborated on by many contributors [13]. Mori and van Zijl [14] and Wong et al. [15] introduced isotropic diffusion encoding by pulsed field gradients in order to sensitize the signal to the direction-average diffusivity without the need for multiple encoding directions, which may be useful in rapid diffusion-weighted imaging [16,17]. Eriksson et al. [18] later demonstrated that a combination of linear and isotropic encoding is sensitive to microscopic anisotropy [19] and that isotropic encoding can be obtained by continuously varying the gradient waveform, a premise that facilitates waveform optimization based on arbitrary gradient trajectories [20].
To emphasize that the encoding is no longer described by a vector (a single gradient direction), we refer to it as tensor-valued diffusion encoding and describe it by a b-tensor. The tensor can be characterized by its trace (b-value) and its symmetry axis (encoding direction). It can also be assigned a shape defined from the b-tensor eigenvalues, similar to shape parameters related to the diffusion tensor [21,22]. The b-tensor formalism was proposed by Westin et al. [21] and generalizes to arbitrary gradient trajectories. Variation of the b-tensor shape is central to probing microscopic anisotropy, as demonstrated for angular DDE that yields linear and planar tensor encoding (LTE and PTE) [12,13,[23][24][25][26][27][28], and for combinations of LTE and spherical tensor encoding (STE) [9,19].
The gradient waveforms necessary to yield, for example, planar or spherical b-tensors, are generally more demanding on the gradient hardware compared to encoding along a single direction. As such, these techniques have primarily been used with high-performance gradient systems that facilitate sufficient data quality and acceptable acquisition times. For example, an early clinical study where microscopic anisotropy was investigated based on a combination of linear and spherical encoding was performed at a 3 T scanner with 80 mT/m gradients-currently considered a high-performance system. Even so, it required an echo time of 160 ms to achieve a b-value of 2.8 ms/μm 2 which limited the spatial resolution and coverage [32]. This warrants development of efficient gradient waveform design [20], parsimonious signal sampling protocols, and investigations of the technical feasibility in MRI systems with different performance.
In this study we aim to survey the technical feasibility of tensor-valued diffusion encoding for DIVIDE analysis, specifically linear and spherical b-tensors, adapted for a whole-brain acquisitions over a wide range of MRI systems. We consider four configurations with field strengths between 1.5 and 7 T, and maximum gradient amplitudes between 45 and 80 mT/m. The most challenging configurations were expected to be those where the SNR was low due to strong T2 relaxation, i.e. configurations where a low maximal gradient amplitudes led to long echo times or with ultra-high field strengths that led to shorter transversal relaxation times [33][34][35].

Methods
We investigated the technical feasibility of tensor-valued diffusion encoding on three MRI scanners with field strengths between 1.5 and 7 T, and maximal gradient amplitudes between 45 and 80 mT/m. The three scanners were used in four configurations (A-D) as described in Table 1. Notably, configurations B and C use the same MR scanner, but configuration B used the protocol and waveforms optimized for configuration A. Furthermore, we emphasize that configuration D employed experimental hardware in a clinical setting.
All configurations were investigated using the clinically relevant requirements of whole-brain coverage (120 mm contiguous coverage in feet-head direction) and acquisition times below 10 minutes. Technical feasibility was quantified in terms of the SNR across the brain parenchyma, Table 1. Configuration of hardware, imaging protocols and gradient waveforms. Configurations A-C used 20-channel receive head/neck coil arrays. Configuration D used a 32-channel transmit/receive head coil array. The maximal gradient waveform amplitude and slew rate allowed by the hardware (G max and S max ) and the maximum that was used for STE and LTE are stated along with imaging and waveform timing parameters.

Protocol design and data acquisition
The protocols were designed using the following three steps: (i) set the minimal and maximal b-values and number of b-values, (ii) determine the minimal echo and repetition times (TE and TR) for each configuration, (iii) distribute the available signal samples across b-values, encoding directions and b-tensor shapes. Details for each step are given in the Supporting Information. Protocols and waveforms were separately optimized for each configuration and are summarized in Table 1. The following imaging parameters were kept constant for all configurations: acquisition matrix 112×112 in 30 contiguous slices, partial-Fourier 0.75, in-plane acceleration factor 2 (A-C: GRAPPA, D: SENSE), bandwidth 1800 Hz/pixel, and fat suppression was engaged (A-C: 'strong', D: 'strong' with slice gradient reversal). Notably, the default spatial resolution was 2×2×4 mm 3 ; the effect of anisotropic voxels is commented on in the discussion. The strength of the diffusion encoding was b = 0.1, 0.7, 1.4 and 2 ms/μm 2 for both LTE and STE. For these b-values, the number of directions for LTE was 6, 6, 12, and 16, respectively. The STE was not rotated since it is assumed to be rotation invariant, but was repeated 6, 6, 12, and 16 times, respectively. This resulted in a total of 80 signal samples per voxel. The sampling order was volume-interleaved to distribute energy consumption and heating more evenly over time [36] and to reduce the potential bias caused by system drift [37]. The total acquisition time (T tot ) was calculated as the number of samples multiplied by TR to remove dependency on preparation phases that may be different across configurations. All experiments were performed using an in-house developed prototype pulse sequence based on the diffusion-weighted spin-echo sequence with echo-planar imaging readout provided by Philips Healthcare (Best, the Netherlands) and Siemens Healthcare (Erlangen, Germany). Asymmetric gradient waveforms for STE were optimized for each configuration to minimize TE using the software described by Sjölund et al. [20], available at https://github. com/jsjol/NOW. The optimization of gradient waveforms used the max norm, heat dissipation factor 0.6, and a slew rate limit of 100 T/m/s, to comply with both duty cycle and peripheral nerve stimulation limits [38,39]. The actual slew rate was determined by the temporal resolution of the waveform optimization as well as the amplitude and timing of the waveform. Although this imposed relatively low slew rates compared to the limit used in optimization, it is expected to have marginal effect on encoding efficiency. As noted in Table 1, the STE gradient amplitude was not maximized because such waveforms could exceed the duty cycle of the scanner. Instead, several amplitudes were tested empirically, and the most efficient waveform that was robustly executed was selected. The waveform optimization included concomitant field effect compensation, and all waveforms had Maxwell indices below 200 (mT/m) 2 ms which ensures negligible signal error [40]. The LTE used symmetric bipolar trapezoid waveforms to approximately match the gradient amplitude and diffusion time of the STE experiment [41]. An example of a spin-echo sequence with both types of waveforms is shown in Fig  1. Comprehensive definitions of waveforms, timing settings, signal sampling schemes and model fit settings are available at https://github.com/filip-szczepankiewicz/Szczepankiewicz_ PONE_2019.
All data were acquired twice within the same session for each configuration, in a single healthy volunteer (male, age 27 y) with previous experience as an MRI volunteer. There was no compensation for volunteering. All experiments were approved by the Regional Ethical Review Board (Lund, Sweden, EPN 2014-735 for systems I and II, EPN 2012-428 for system III), and informed consent was obtained prior to participation. Diffusion-weighted data were corrected for motion and eddy currents in Elastix [42] using an extrapolation-based reference volume [43]. To retain the impact of signal noise on the images, no smoothing was applied at any stage of the processing.
To investigate the interplay between SNR and spatial resolution, we created two alternative configurations. Based on configuration A, we mitigated the low SNR by reducing the spatial resolution. A low-resolution configuration (A � ) was achieved by setting the matrix size to 88×88 resulting in 2.5×2.5×4 mm 3 voxels and the bandwidth was reduced to 1180 Hz/pixel to preserve the readout time and total acquisition time. Based on configuration C, we sacrificed SNR in favor of an increased spatial resolution. A high-resolution configuration (C � ) was achieved by reducing the slice thickness to 2 mm. To preserve coverage, the number of slices was doubled and the repetition time was set to 8900 ms. This yielded a spatial resolution of 2×2×2 mm 3 and an acquisition time of 12 minutes.

DIVIDE parameter estimation
DIVIDE employs a joint analysis of data acquired with b-tensors (B) with multiple b-values (b = Trace(B)) and shapes (b Δ ) [22]. Assuming that the observed diffusion process can be approximated by a distribution of diffusion tensors (multi-Gaussian) [4,44], the fitting parameters can be interpreted in terms of the signal at b = 0 (S 0 ), mean diffusivity (MD) and the anisotropic and isotropic diffusional variance (V A and V I ) [19,32]. Although several approaches for estimation of these parameters exist [21,24,29,45], we estimate the parameters by fitting the Laplace transform of the gamma distribution [46,47] to the powder-averaged LTE and STE signal ( � S) simultaneously [19,32], according to For LTE, the powder average is the arithmetic average across directions for each b-value (see Supporting Information), whereas STE is natively rotation invariant and requires no such averaging. Note that Eq 1 is only valid for axisymmetric b-tensors, and that b Δ is 1 and 0 for LTE and STE, respectively [48]. The average and variance of the gamma distribution produced by the fit are estimates of the average and variance of the true distribution of diffusion coefficients. If the gamma and true distributions have different functional forms, the estimates may be biased. The fitting was initialized at random values, accounted for heteroscedasticity (weighted by the square root of the number of signal samples that make up each average), and was constrained to yield real-but not necessarily positive-values for V I or V A . Negative Note that the STE waveform is asymmetric around the refocusing pulse [20] and that it must therefore be designed to compensate for errors caused by concomitant fields [40]. The LTE waveform is bipolar to match the effective diffusion times for STE and LTE [41]. Note that crusher gradients (dotted lines) are not engaged when the diffusion encoding acts as a crusher. https://doi.org/10.1371/journal.pone.0214238.g001 variances are not physical but can be caused by noise and were allowed to avoid a positive bias for values close to zero. The diffusional variance was reported in terms of a normalized metric similar to the mean kurtosis from DKI [4], defined according to [19,49] where 'x' denotes the isotropic (I) or anisotropic (A) diffusional variance. The normalized signal at b = 2 ms/μm 2 averaged over directions or repetitions ( � S=S 0 ) was also considered in the analysis of repeatability, and denoted A(L) and A(S) when using LTE and STE, respectively. Finally, we calculated the microscopic fractional anisotropy (μFA), according to which can be interpreted as the FA that would be observed if all structures in the sample were aligned in parallel [9,21]. Thus, μFA captures the microscopic diffusional anisotropy even in samples that are isotropic on the voxel scale [24,32,50]. We use the definition by Westin et al. [21] (including the contribution from V I ) rather than the one by Lasič et al. [19]. The post-processing and analysis software [51] was written in Matlab (The MathWorks, Natick, MA, USA) and is available at https://github.com/markus-nilsson/md-dmri.

Analysis of repeatability
For each configuration, data from two identical acquisitions were compared to investigate repeatability. Spatial correspondence across acquisitions was achieved by not repositioning the subject between acquisitions. To gauge the signal and DIVIDE parameter repeatability, we calculated the differences between the first and second acquisition (ΔX = X 1 -X 2 ) [52], such that ΔX is a distribution of voxel-wise differences. The repeatability was visualized in maps of ΔX as well as Bland-Altman plots. The repeatability was summarized by the mean and standard deviation of ΔX to capture the overall parameter bias and precision. To avoid inflated variability due to a misaligned brain periphery and partial-volume effects with CSF, only voxels where μFA > 0.7 and MD < 1.5 μm 2 /ms were considered.

Analysis of SNR
At low SNR, the MR signal is positively biased due to the rectified noise floor, with detrimental effects on parameter accuracy [53][54][55]. Since STE was repeated several times for each b-value, we could estimate the SNR at each b-value by computing the ratio between the mean of the STE signal and its standard deviation. By assuming that signal is approximately Rice distributed, we can use a threshold of SNR < 3 to mark regions where signal bias is likely to influence signal accuracy [56]. Doing so may identify regions of the brain that require action to improve SNR. To avoid overestimation of the SNR, it was calculated based on image data that had not been motion corrected. As a summary parameter of data quality, we report the fraction of the brain parenchyma where SNR was above 3 and 6 at the highest b-value (b = 2 ms/μm 2 ), denoted Q 3 and Q 6 . Only voxels where MD < 1.5 μm 2 /ms were included in the analysis to avoid regions dominated by cerebrospinal fluid.

Simulation of parameter accuracy and precision
Simulations were used to investigate the impact of signal noise on the bias, variability and negative values of DIVIDE parameters. To this end, we used toy models of tissue described by two Noise sampled from a Rice distribution was added to yield SNR between 10 and 100 in S 0 . Each simulation used the same signal sampling scheme and fitting procedure as for the in vivo case. Simulations were repeated 10 4 times using independent realizations of noise for each SNR value. The mean and standard deviation of model parameters across realizations was calculated to quantify their accuracy and precision.

Results
Tensor-valued diffusion encoding was successfully performed using all four configurations. Fig 2 shows the calculated DIVIDE parameter maps. There is an appreciable difference in image quality depending on hardware performance. As expected, higher field strengths and stronger gradients generally render lower echo times and higher image quality. The influence of noise was most prominent for configuration A. This is likely due to the long echo time (TE = 115 ms) incurred by the relatively low gradient amplitude. By contrast, SNR was relatively high for configuration D, but its parameter maps exhibited contrast that differed from other configurations. For example, MK I was negative in large regions throughout the parenchyma, and MK A was higher than at other systems, most prominently in coherent white matter pathways. Image artifacts were also noticeable, such as geometric distortions, likely associated to issues commonly observed for dMRI at ultra-high fields [57]. Fig 3 shows maps of the SNR at the highest b-value (b = 2 ms/μm 2 ) in three axial slices for all configurations. Configuration A exhibited SNR < 3 close to the lateral ventricles and in the inferior parts of the brain. Configurations B-D exhibited SNR > 3 in more than 96% of the brain parenchyma. A complete list of the quality parameters is reported in Table 2. Fig 4 shows the parameters estimated from simulations of noisy signal. The results show an association between SNR and parameter accuracy as well as precision. Overall, signal noise is a plausible cause for the parameter variations seen in the parameter maps in vivo, including occurrences of negative diffusional variance, especially at low SNR and for true values close to zero. The precision in MK I was consistently lower than for MK A , regardless of their true values, in accordance with in vivo results (Table 2). Moreover, the simulations suggested that regions where SNR is low may render a positive parameter bias, especially in MK I . None of the tested cases showed that noise can induce a negative bias in MK I . Simulations of tissue where the true MK I = 0 indicated that noise is unlikely to account for the positive MK I values found throughout the brain. This indicates that the true isotropic diffusional variance is detectable and positive in brain parenchyma.
Resulting parameter and SNR maps at lower and higher resolutions on configurations A � and C � , are shown in Fig 5. As expected, the larger voxels used at configuration A � improved the SNR. For the lower resolution (2.5×2.5×4 mm 3 voxels), we found Q 3 = 99%, Q 6 = 61%, compared to Q 3 = 85%, Q 6 = 16% with the baseline resolution (2×2×4 mm 3 voxels). For configuration C � , the smaller voxels (2×2×2 mm 3 ) reduced the overall SNR, but it remained relatively high at Q 3 = 94% and Q 6 = 45%. These results show that protocols with image resolutions tailored to the system performance are also technically feasible.

Fig 2. DIVIDE parameter maps in transversal and coronal slices for configurations A-D.
The image quality generally improves with higher field strength and gradient amplitude. Most notably there is a discernably higher level of noise for configuration A. Configuration D generally exhibited higher MK A , a pronounced geometrical distortion at the anterior part of the brain, and more pronounced ghosting artifacts, as compared with the other configurations. Furthermore, large regions of negative MK I were observed only for configuration D.

SNR maps at b = 2 ms/μm 2 in three transversal slices for configurations A-D and SNR distributions in the brain parenchyma at b = 0.1 and 2 ms/μm 2 (histograms).
White outlines show the outer perimeter within which SNR > 3, and the red outlines show regions where SNR < 3. At high b-values, images from all systems exhibited low SNR in the ventricles due to the high diffusivity of CSF. Configuration A shows low SNR in the central and inferior parts of the brain. Configuration C exhibited the highest SNR across all bvalues, but also a posterior rim of low SNR caused by poor fat suppression (red arrows). Configuration D generally exhibited high SNR in the peripheral regions; in the lower parts of the brain, it exhibited a heterogeneous SNR. It also exhibited an irregular perimeter in the inferior part of the brain, which was likely caused by local field heterogeneity in proximity to the ear canal and poor RF homogeneity, both commonly associated to dMRI at ultra-high field strengths [58].
https://doi.org/10.1371/journal.pone.0214238.g003 Table 2. Quality and repeatability parameters for configurations A-D. Q 3 and Q 6 are the fraction of tissue where SNR is above 3 and 6 at b = 2 ms/μm 2 . The voxel-wise parameter difference between acquisition 1 and 2 (ΔX) is reported as the average ± one standard deviation. The lowest and highest SNR was observed for configurations A and C, respectively, and parameter precision follows the same pattern. Neither the signal attenuation nor DIVIDE parameters had a relevant bias compared to their standard deviation. Furthermore, MK I consistently exhibits a larger standard deviation compared to MK A . agreement with the simulations (Fig 4). Notably, the precision of MD and μFA are both visibly dependent on their absolute values, which warrants a careful interpretation of the global standard deviation. The mean and standard deviation of voxel-wise differences are reported in Table 2.

Discussion
We have demonstrated that tensor-valued diffusion encoding with high b-values is technically feasible across a wide range of different MRI scanners and configurations, and we have shown the range of results that can be expected in terms of DIVIDE parameter maps. These results can serve as a baseline in the design of future experiments. In addition to presenting imaging protocols with high efficiency, this work also demonstrates that scanners with lower gradient strengths-not commonly associated with 'advanced diffusion imaging'-have sufficient performance and stability for tensor-valued diffusion encoding by asymmetric gradient waveforms at reasonable SNR and resolution. Even at 1.5 T and a maximum gradient waveform amplitude of 33 mT/m (configuration A), a minor reduction in spatial resolution was enough to facilitate relatively high SNR. Likewise, the high performance of configuration C could be harnessed to access a higher resolution, albeit at a prolonged scan time related to an increase in number of slices (Fig 5). We acknowledge that the present implementation at configuration D (7 T) displays a prominent inconsistency, namely repeatable regions where MK I is negative, which are not physical and cannot be explained by low SNR (Fig 4). The current results cannot be contrasted to literature since, to the best of our knowledge, b-tensor encoding at 7 T has not been investigated on a group level in vivo. We speculate that this may be caused by, subject motion or poor eddycurrent compensation. It may also be related to field-dependent relaxation times in different tissue compartments [59][60][61][62] which may cause large values of MK A and an underestimation of MK I , predicated on a longer decay time in the most anisotropic compartment. Although further investigations are warranted to understand the source and impact of this bias, data that are estimated for a noise-free signal; deviation from the line indicates parameter bias caused by noise. As expected, precision and accuracy both improve with increasing SNR. The estimation of MD and μFA appear to be the most accurate and precise, although the μFA shows a deterioration of accuracy and precision when its true value is low (panel II) [19]. Both MK A and MK I suffer a positive bias when SNR in the b = 0 image is approximately 20 or less. Interestingly, MK A is always more precise and accurate than MK I , indicating that it is generally less sensitive to noise. The simulations show that signal noise can cause negative values for both MK I and MK A , which is especially likely when the true values are close to zero. Finally, in the case where MK I = 0 (panel III), signal noise did not cause a strong positive bias in MK I for SNR levels that match configuration C (where the majority of voxels have SNR > 20 at b = 0.1 ms/μm 2 , Fig 3). This suggests that signal noise alone is not likely to explain the positive MK I that is observed throughout the brain parenchyma. https://doi.org/10.1371/journal.pone.0214238.g004 Tensor-valued diffusion encoding for diffusional variance decomposition suggests that tensor-valued diffusion encoding is feasible also at ultra-high fields if the bias can be neglected or addressed.
Simulations based on realistic signal characteristics and SNR levels established a relation between SNR and the accuracy and precision of DIVIDE parameters. As shown in the parameter maps and simulations, low SNR introduces a positive parameter bias, most prominently seen for MK I . This is especially relevant in the central parts of the brain where the SNR tends to be lower compared to peripheral tissue. The simulations also showed that noise had a markedly smaller impact on MK A compared with MK I , a feature that can be appreciated also in the parameter maps.
The technical feasibility, especially in a clinical setting, hinges mainly on the use of parsimonious sampling protocols and efficient gradient waveforms. In this study we achieved efficient tensor-valued encoding by using optimized asymmetric gradient waveforms [20]. The benefit of using the currently proposed waveforms and protocols can be appreciated by comparing them to our initial implementation [32], where a 10-minute protocol at a 3T scanner with 80 mT/m gradients covered only 5 slices at a resolution of 3×3×3 mm 3 . For reference, using the previous waveform design [18] at configurations A, C, and D would render echo times of For configuration C � , the resolution was increased, while maintaining high SNR in the superior and peripheral parts of the brain, although inferior and central parts showed regions where SNR was below 3 and may therefore suffer from non-negligible signal bias. https://doi.org/10.1371/journal.pone.0214238.g005 Tensor-valued diffusion encoding for diffusional variance decomposition approximately 230, 140 and 170 ms, respectively. This translates to a loss in SNR of 60%, 50% and 70%, assuming white matter relaxation times presented by Cox and Gowland [33]. Although asymmetric designs excel in efficiency, they are susceptible to errors caused by concomitant fields [63,64]. Such errors can be difficult to detect by the naked eye, especially in vivo, but lead to quantification errors for several popular asymmetric gradient designs. To avoid this, we used so called Maxwell-compensated waveforms [40].
We emphasize that this study surveys technical feasibility, whereas investigation of the clinical feasibility, pertaining to specific applications, was outside the scope of this work. Preliminary work suggests that tensor-valued diffusion encoding may add novel information to investigations of schizophrenia [21], brain tumors [9], multiple sclerosis [28,65], cortical malformations [66], prostate tumors [67], microstructure imaging of the healthy brain [32,68], and kidneys [69]. The present protocol design was intended for investigations of healthy brain tissue, and therefore accounts for the approximate diffusivity and anisotropy of brain tissue, as described in the Supporting Information. Naturally, investigations of tissues with vastly different characteristics may require a different protocol. For example, a protocol designed for tumors where MD � 1.8 μm 2 /ms and FA � 0.1 [9] would generally entail a lower maximal bvalue (b max < 1.3 ms/μm 2 ) and fewer diffusion-encoding directions (n min = 3 at b max )-a markedly different premise for optimization compared to normal brain tissue (Supporting Information). A comprehensive protocol design should also consider the relaxation characteristics of the tissue. For example, transversal relaxation rates and SNR may depend on the tissue iron content [70,71], which is associated with ageing. Specific applications can also benefit substantially from adapting the image quality and scan time to the expected level of variance in the studied population [72]. Tensor-valued diffusion encoding can be accelerated by simultaneous multi-slice imaging [73] and interleaving techniques [36] that reduce acquisition times and improve duty cycle constraints, as well as improvement in terms of reducing eddy-current effects [74,75].
In this study we used DIVIDE to investigate repeatability of parameters related to tissue microstructure which have a straightforward interpretation if effects of diffusion time and intra-compartment kurtosis are both negligible. Although studies of the brain performed at relatively long diffusion times indicate that such effects are subtle [76][77][78], we emphasize that the multi-Gaussian framework may be incomplete [79], as demonstrated by specialized diffusion encoding schemes that have detected time-dependent diffusion coefficients in biological tissues [80][81][82][83][84][85]. Future studies will aim to design waveforms to probe the interplay between diffusion-time dependence, compartment kurtosis and exchange [41,79,86,87].
We acknowledge several limitations of the present study. First, we used asymmetric voxels (2×2×4 mm 3 ) as an effective means to increase the SNR and coverage of the protocols. Asymmetric voxels are known to introduce confounding effects in measures of voxel-scale diffusion anisotropy, for example in DTI, due to the interaction between voxel and structure geometry [88]. A similar limitation is not true for DIVIDE since it renders parameters that are independent of the orientation dispersion of the underlying tissue [19]. However, other issues pertaining to asymmetric voxels remain, such as increased partial-volume effects in the through-plane The voxel-wise difference (Diff.) between the first and second acquisition (Acq. 1 and 2) is color coded in red-green. The normalized powder-averaged signal at b = 2 ms/μm 2 (A(L) and A(S)) is given in percent, and the MD is given in units of μm 2 /ms. The difference maps show the largest differences in tissue interfaces where small misregistration between the first and second acquisition causes large parameter discrepancy. Bland-Altman plots show the distributions of voxel-wise differences in tissue where μFA < 0.7 and MD < 1.5 μm 2 /ms. Solid and dashed lines show the average and two standard deviations of the distributions. All configurations showed negligible bias in repeatability of signal and DIVIDE parameters, and a configuration-dependent parameter precision. Note that the estimated precision pertains to the per-voxel parameter uncertainty; analyzing the average over multiple voxels is expected to markedly improve the precision. https://doi.org/10.1371/journal.pone.0214238.g006 Tensor-valued diffusion encoding for diffusional variance decomposition direction. This is especially relevant if the data is also intended to support tractography, where isotropic spatial resolution is preferable [89,90]. Second, this study aimed to represent a wide range of MRI scanners, but several relevant hardware features were omitted in this study. For example, we did not account for the effects of coil design, RF and gradient systems, field heterogeneity, shimming, image acquisition, or reconstruction technique [34,57,58,91]. Furthermore, the settings for the waveform optimization are not likely to be optimal because optimality depends on experimental settings, hardware configuration, and duty cycle, all of which may be site specific. Nevertheless, the waveforms suggested herein may serve as a starting point for future optimization. Finally, we emphasize that although the included MRI systems span a wide range of capabilities, current results are based on a single healthy volunteer, a single analysis method, post-processing pipeline, and a limited set of scanners, and may therefore lack generalizability. However, two critical aspects of technical feasibility and robustness [92], namely repeatability and sufficiently high SNR, could be established. The test-retest data in Fig 6 show that the tensor-valued encoding yields repeatable results for all configurations, except for the normalized signal for configuration D. The generalizability can also be strengthened by a separate analysis that we performed on a cohort of ten subjects, previously acquired with hardware specifications similar to configuration B and used in a proof-of-concept implementation [20]. Although that study had limited coverage (44 mm feet-head) and longer echo time (TE = 130 ms), the experiments yielded homogeneous data quality on the group level, where Q 3 = 95 ± 2% and Q 6 = 21 ± 4%; consistent with the current quality estimated at configuration B when considering the longer TE.

Conclusions
Tensor-valued diffusion encoding can probe microstructural features beyond those available with conventional methods, but it requires non-conventional diffusion-encoding waveforms. The preferred platform for in vivo experiments based on such experiments has been 3 T scanners with high-performance gradients. In this study, we demonstrated that tensor-valued encoding that supports DIVIDE analysis is technically feasible over a wide range of MRI scanners, with main magnetic fields between 1.5 and 7 T, and gradient waveform amplitudes as low as 33 mT/m. We have also reported baseline repeatability values for whole-brain DIVIDE protocols with acquisition times between 5 and 9 minutes. The implementation was facilitated by efficient asymmetric gradient waveforms and parsimonious signal sampling protocols. By enabling tensor-valued diffusion encoding and DIVIDE at a wide range of scanners at clinically acceptable acquisition times, we expect that such methods may be more broadly used to facilitate new and exciting venues for dMRI research.

S1 Fig. Minimal number of directions (n min ) required to yield a rotation invariant powderaveraged signal for linear and planar tensor encoding (LTE and PTE)
. n min can be read from the figure for combinations of the tissue fractional anisotropy (FA) and attenuation factor (b�MD), where the numbers in the circles show n min for a color-coded interval. As expected, higher anisotropy and attenuation both demand a larger number of diffusion-encoding directions. Spherical tensor encoding always requires n min = 1, since it is inherently invariant to rotation. Note that we consider the signal to be rotation invariant when CV < 1%. The PTE plot is included for completeness, even if it was not used in the data acquisition of this study. (DOCX)