Spectral principal axis system (SPAS) and tuning of tensor-valued encoding for microscopic anisotropy and time-dependent diffusion in the rat brain Open Access

Diffusion magnetic resonance imaging (dMRI) is commonly used as a non-invasive probe of brain microstructure (Tournier et al., 2011). The technique is very sensitive to changes in microstructure, but its specificity is limited because the dMRI signal reflects a mix of tissue characteristics, including cell shape, size, orientation, and permeability. This problem is exacerbated when the tissue is heterogeneous, that is, when the relatively large imaging voxel contains a mixture of tissues with different microstructural features (Jeurissen et al., 2013; Mulkern et al., 2000). One approach to increase specificity is to employ increasingly more sophisticated biophysical models (Novikov, Kiselev, et al., 2018; Olesen et al., 2022; Palombo et al., 2020), which often rely on fragile assumptions (Jelescu et al., 2020; Jones et al., 2013; Lampinen et al., 2023). Another approach is to use encoding strategies that aim to increase specificity at the acquisition stage (Topgaard, 2020). The goal is to disentangle different confounding factors, such as molecular exchange (Callaghan & Furó, 2004; Chakwizira, Westin, et al., 2023; Lasič et al., 2011), restricted and hindered diffusion (Balinov et al., 1993; Lasič et al., 2006, 2009; Neuman, 1974; Stepišnik, 1993), incoherent flow (Ahlgren et al., 2016), microscopic and macroscopic anisotropy (Jespersen et al., 2013; Lasič et al., 2014; Nilsson et al., 2020; Szczepankiewicz et al., 2016; Topgaard, 2017, 2019; Westin et al., 2014, 2016), and relaxation (de Almeida Martins et al., 2020; De Almeida Martins & Topgaard, 2018; Lampinen, Szczepankiewicz, et al., 2020; Veraart et al., 2018).

Traditional diffusion encoding employs pulse field gradients (PFG) (Stejskal, 1965; Stejskal & Tanner, 1965) to sensitize the signal to diffusion along a single direction. Such encoding is used in diffusion tensor imaging (DTI) (Basser et al., 1994) and diffusion kurtosis imaging (DKI) (Jensen et al., 2005). However, these approaches entangle information about tissue anisotropy on microscopic and macroscopic scales. In tensor-valued encoding (Jespersen et al., 2013, 2014; Lasič et al., 2014; Nilsson et al., 2020; Szczepankiewicz et al., 2016; Topgaard, 2017, 2019; Westin et al., 2016), gradients are played out simultaneously but asynchronously along multiple orthogonal directions to yield b-tensors of varying shapes, which enables microstructure assessment unconfounded by the orientation distribution of anisotropic microstructures and morphological heterogeneity within macroscopic imaging voxels. Tensor-valued encoding can probe microscopic fractional anisotropy (µFA), for example, by using a combination of linear tensor encoding (LTE) and planar tensor encoding (PTE) (Jespersen et al., 2013) or with LTE and spherical tensor encoding (STE) (Lasič et al., 2014; Szczepankiewicz et al., 2016). Adopting more general approaches to tensor-valued encoding, not limited to the PFG paradigm, opens up venues for optimizations in terms of b-tensor shapes used for different hardware and physiological constraints (Szczepankiewicz, Westin, et al., 2021) and encoding waveforms suitable to maximize encoding strength (Sjölund et al., 2015), minimize effects of concomitant fields (Szczepankiewicz, Westin, et al., 2019), background gradients (Szczepankiewicz & Sjölund, 2021), or enable varying degree of motion compensation (Lasič et al., 2020; Szczepankiewicz, Teh, et al., 2021).

Analysis of tensor-valued encoding data often assumes multi-Gaussian diffusion. This means that the measured signals depend only on the b-tensor size and shape but not on the temporal characteristics of the diffusion encoding gradient waveforms. However, restricted diffusion in tissue could cause time-dependent diffusion (TDD), where the signals and the associated apparent diffusivities depend on the encoding waveform time scales. Different b-tensors may utilize encoding waveforms with different time scales and thereby feature different sensitivities to TDD (Lasič et al., 2022; Lundell et al., 2018, 2020; Lundell, Nilsson, Dyrby, et al., 2019; Lundell & Lasič, 2020; Teh et al., 2021). Combining b-tensors with different TDD sensitivities could significantly bias estimation of microscopic anisotropy (µA) or the separation of isotropic and anisotropic sources of the mean kurtosis. Such conditions may be prominent in preclinical settings using shorter encoding times (Ianuş et al., 2018; Lundell, Nilsson, Dyrby, et al., 2019) or in tissues with larger cells like in the heart (Lundell et al., 2020; Teh et al., 2021). For unbiased assessment of µA, b-tensors need to be tuned for equal TDD sensitivity to yield equal mean diffusivities (MD) (Lundell et al., 2020; Lundell, Nilsson, Dyrby, et al., 2019; Lundell & Lasič, 2020).

This work considers TDD due to restricted diffusion in tensor-valued diffusion encoding with two main objectives: (1) to reduce bias in µA assessment and (2) to obtain two independent contrasts due to µA and TDD using a minimal set of diffusion encoding waveforms. We propose two approaches for tuning LTE to an arbitrary b-tensor shape.

The first approach uses a shorter protocol with a single LTE waveform, ideal when the primary goal is to assess µA unconfounded by TDD. The tuned LTE waveform is selected as the one-dimensional projection of an STE waveform that yields the TDD profile most similar to that of the STE. This approach provides tuning relative to the sample, thus requiring prior knowledge of the expected range of restriction sizes. The second approach involves a longer protocol with multiple LTE waveforms, using geometric averaging to obtain the tuned LTE data. This method achieves tuning that is independent of restricted diffusion effects in the tissue. In addition, we aim to maximize the experimental yield, by allowing acquisition of independent contrasts sensitive to µA and TDD within a single multidimensional protocol. This is achieved by defining the spectral principal axis system (SPAS) of encoding (Lasič et al., 2022) and combining it with the geometric averaging approach. The effectiveness of these methods is demonstrated through simulations and experiments on both perfusion-fixed and in vivo rat brains.

In this section, we present a framework for reducing the confounding effects of time-dependent diffusion (TDD) in tensor-valued encoding for microscopic anisotropy (µA) assessment. We also introduce the spectral principal axis system (SPAS) to probe both µA and TDD using the same set of encoding waveforms.

We begin by revisiting the principles of tensor-valued encoding, highlighting the need to account for TDD due to restricted diffusion. Using frequency-domain analysis under the Gaussian phase approximation (GPA), we identify two critical encoding attributes: spectral trace and spectral anisotropy. Spectral trace is the key parameter for tuning b-tensors, while spectral anisotropy reflects how sensitivity to TDD varies with encoding orientation. These concepts underpin the SPAS, which defines an orthogonal set of projections optimized for distinct TDD sensitivities.

We then outline two approaches for tuning: one based on matching apparent mean diffusivity (MD) and another using geometric signal averaging. Both strategies help mitigate the confounding effect of TDD in µA assessment, with geometric averaging offering more robust tuning and enabling simultaneous probing of µA and TDD within a single encoding protocol.

When imaging voxels contain multiple compartments with different diffusivities, anisotropies, and orientations, the signal attenuation is generally multiexponential. At low b-values, the signal attenuation can be approximated by a cumulant expansion

where MD is the apparent mean diffusivity and $VD$ is the apparent diffusion variance, related to the “mean kurtosis” (MK) by $VD=13MD2MK$ (Szczepankiewicz et al., 2016; Westin et al., 2016). Note that MD and $VD$ are rotationally invariant measures derived from the powder-averaged signal. In tensor-valued encoding, the anisotropy or shape of the b-tensors provides an independent encoding dimension, allowing us to separate $VD$ into isotopic (⁠$VI$⁠) and anisotropic (⁠$VA$⁠) contributions, $VD=VI+VA.$ For axially symmetric encoding, we have

where $bΔ2$ is 0 for spherical (STE) and 1 for linear (LTE) tensor encoding (Eriksson et al., 2015; Topgaard, 2017). By varying the b-tensor shape, we can probe $VA$⁠, assuming the encoding of MD and $VI$ remains independent of b-tensor shape, which holds when the diffusion process is approximately Gaussian.

Now, consider a combination of STE and LTE encodings, commonly used to assess microscopic anisotropy (Lasič et al., 2014; Nilsson et al., 2020; Szczepankiewicz et al., 2016; Topgaard, 2017, 2019; Westin et al., 2016). The signal difference between STE and LTE at a given b-value is approximately

where we define the microscopic anisotropy contrast as

For Gaussian diffusion, $ΔMD=0$ and $ΔVD$ is proportional to the microscopic fractional anisotropy (µFA) (Jespersen et al., 2013; Lasič et al., 2014), defined as

In contrast to FA, µFA removes the effect of orientation dispersion, which is related to the order parameter (Lasič et al., 2014). Note that the definition of µFA above differs slightly from that in Szczepankiewicz et al. (2016), where the isotropic diffusion variance was added to $MD2$⁠. Furthermore, µFA² maps, visualized in this work, scale similarly as the µA maps as a function of $ΔVD$⁠.

Restricted diffusion may lead to $ΔMD≠0$⁠, introducing bias in the assessment of µA and µFA. For example, if LTE measures diffusion at longer time scales than STE, the apparent MD from LTE may be lower due to restricted diffusion, leading to $ΔMD<0.$ This could increase the signal separation in Eq. (3), falsely suggesting higher µA and positively biasing µFA estimates.

An illustrative example of the potential bias introduced by combining “oscillating” STE gradients (Topgaard, 2013) with more efficient pulsed-like LTE gradients (Lasič et al., 2014) is provided by experiments in densely packed isotropic yeast cells (Lundell, Nilsson, Dyrby, et al., 2019). In these experiments, high signal differences could have been erroneously interpreted as high cell anisotropy if TDD effects were not considered. Such a combination of different gradient timings for LTE and STE is reminiscent of experiments employing oscillating gradients (Callaghan & Stepišnik, 1996; Gore et al., 2001; Stepišnik, 1993, 1998, 1999) or a combination of oscillating and pulsed gradients (X. Jiang et al., 2017; Reynaud et al., 2016a; Reynaud, 2017) designed to probe restriction sizes. Although combining waveforms with different sensitivities to TDD is useful for probing cell sizes, it is generally undesirable in tensor-valued encoding. This underscores the importance of carefully tuning b-tensors to minimize the confounding effects of TDD to ensure unbiased assessments of µA. The following sections will address this challenge by introducing a frequency-domain analysis that forms the basis for understanding and managing TDD effects in tensor-valued encoding.

The effects of TDD due to restricted diffusion can be analyzed in the frequency domain using the Gaussian phase approximation (GPA) (Ianuş et al., 2013; Stepišnik, 1993, 1998, 1999). This approach can also be extended to the asynchronous diffusion encoding waveforms characteristic of tensor-valued encoding (Lundell & Lasič, 2020). In the first-order approximation, TDD is examined by analyzing the encoding and diffusion power spectra. Unlike LTE, tensor-valued encoding employs asynchronous gradients $g(t)$ along orthogonal axes, which are essential for tensorial diffusion weighting. In frequency domain, the signal attenuation is given via coupling of the dephasing and velocity cross-power spectral densities. For brevity we call them encoding and the diffusion power spectra (Lundell & Lasič, 2020). With the dephasing vector

and its Fourier transform

the tensorial 3 x 3 (i x j) components of the encoding power spectrum are defined as

The b-tensor components are given by

and the b-value by

where

is the encoding spectral trace (see Fig. 1A). This is the key encoding attribute for tuning different b-tensors. Since $s(ω)$ is real valued and symmetric around 0, we can use 0 as the lower frequency bound in Eq. (10).

The directional variation of the encoding power distribution is represented by spectral projections along a unit vector u as

This directional dependence of the encoding power distribution defines spectral anisotropy (SA), which reflects how sensitivity to TDD varies along different axes in tensor-valued encoding (Fig. 1). While the b-value (b-tensor trace) and b-tensor anisotropy suffice to characterize Gaussian diffusion, tuning (spectral trace) and spectral anisotropy add crucial encoding information to address TDD effects (Lundell & Lasič, 2020). These concepts are central to the encoding strategies discussed in the following sections.

Consider a single sub-ensemble of spins (tissue compartment) in the sense of an ergodic stochastic process of positions characterized by a common velocity autocorrelation function (Lundell & Lasič, 2020; Van Kampen, 2001). Such compartments could for example reside within intra- or extra-cellular spaces (Lasič et al., 2009; Novikov et al., 2011; Stepišnik, 1993) or represent incoherent flow (Kennan et al., 1994). With GPA (Ianuş et al., 2013; Stepišnik, 1993, 1998, 1999), the attenuation factor $β=−ln(S/S0)$ is obtained by integrating the inner product of the encoding and diffusion spectra over the entire frequency range,

Note that the symmetry of the encoding spectra, $sij(ω)=s¯ji(ω)$⁠, could be used to rewrite Eq. (13) with 0 as the low integration limit (see equations 2.100–2.103 in Lundell & Lasič, 2020). The diffusion spectrum from a single compartment can be expressed with rotation matrices R as

where $λk(ω)$ are diffusion eigenspectra (Lasič et al., 2009, 2022; Stepišnik, 1993). This expression is useful to evaluate ensemble average signals (Lundell & Lasič, 2020) and was used in the simulations presented in Figure 4.

In this section we elaborate on what we mean by tuning and suggest two approaches for tuning considering STE as the tuning reference and LTE as a tensor being tuned to the reference. Although we demonstrate the approaches for the special cases of STE and LTE, these approaches could be adapted to account for arbitrary b-tensor shapes.

Different b-tensors can be considered tuned to each other if they yield approximately equal MD values also in the presence of TDD. Tuning helps eliminate the first order confounding effect of TDD by ensuring $ΔMD=0$ in Eq. (3). The MD is derived from the first cumulant of the ensemble and direction average signal,

where $λ¯iso(ω)$ is the average isotropic diffusion spectrum, defined as

representing the average apparent isotropic diffusivity$D¯iso$ (Lundell & Lasič, 2020). The overbar denotes averaging over compartments with different eigen-spectra (Lasič et al., 2022). Eq. (15) highlights a key point: matching spectral traces is essential for tuning. We will demonstrate that this can be achieved through two different approaches.

Different LTEs can be derived from STE by considering waveform projections

with encoding power spectra $su(ω)$ from Eq. (12). For a powder sample, such LTEs would yield the attenuation factor as

The tuned LTE can be obtained by finding LTE projections yielding apparent diffusivity $Du$ equal to the MD from STE. This requires assuming a specific diffusion spectrum, for example, for spherical restrictions (see Figs. 1B and 2). To better suite hardware constraints, alternative projections, which sacrifice some tuning for reduced gradient amplitude, slew rate or their product, $max{|gu|}⋅max{|g˙u|}$⁠, may be advantageous, as shown in Figure 2 for the optLTE (Lasič et al., 2023).

Such tuning approach does not aim to directly match the spectral traces (see Eqs. (11) and (15)), but rather to match the apparent mean diffusivities, which depends on both the encoding and diffusion power spectra. Although this does not represent a universal tuning strategy, the smooth nature of the encoding and diffusion spectra warrant relative robustness of the proposed approach, where tuning is approximately invariant in a wide range of restriction sizes (see Fig. 3).

As we see from Eq. (15), given equal b-values, two otherwise different b-tensors are considered tuned if their spectral traces $s(ω)$ are exactly matching (see panels Figs. 1A and 2B). Such requirement is more stringent compared with matching MDs, but it ensures tuning independent of diffusion spectra. While matching spectral traces might be challenging, an equivalent tuning can be attained via geometrical averaging of signals from any orthogonal set of tensor projections. Consider a powder sample and the geometric/arithmetic $(〈⋯〉g​/​〈⋯〉a)$ averages of signals from LTE projections along orthogonal directions. Since the geometric average signal is in the first approximation given by the arithmetic average attenuation factor,

the geometric average of signals from LTE projections corresponds to the signal encoded with the spectral trace and yields equal MD values. To match the encoding powers (b-values), the LTE gradient waveforms need to be appropriately scaled as $gLTE(t)=c1/2 gu(t)$⁠, where $c=bbu=∫0∞s(ω) dω∫0∞su(ω) dω$⁠, which is equal to 3 when STE is the tuning reference. The geometric average signal in Eq. (19) is expected to approximate well the tuned LTE signal up to the second cumulant (diffusion variance), assuming negligible effects of intra-compartmental kurtosis and exchange.

In addition to the tuning property, b-tensors generally exhibit varying sensitivities to TDD along different tensor axes (Lasič et al., 2021; Lundell et al., 2018; Lundell & Lasič, 2020; Teh et al., 2021), a property referred to as the spectral anisotropy (SA) (Lundell & Lasič, 2020). SA represents the directional variation in encoding power distribution, as defined in Eq. (12). This can be visualized using color coding (Fig. 1), where the relative weights for red, green, and blue (RGB) are computed for each projection as $∫ω∈Ωsu(ω) dω∫0∞su(ω) dω$⁠, where $Ω$ represents the low-, mid-, and high-frequency bands. In the presence of TDD, SA may cause spherical tensor encoding (STE) to lose its rotational invariance, meaning it no longer yields truly isotropic diffusion encoding but instead depends on the relative orientation of the encoding (q-trajectory) and the diffusion compartment (Jespersen et al., 2019; Lasič et al., 2021). As a result, when averaging over orientations (powder averaging), SA can introduce additional variance of apparent diffusivities (⁠$VD$ in Eq. (1)) in samples with anisotropic structures, potentially biasing the assessment of microscopic anisotropy (Jespersen et al., 2019; Lasič et al., 2021; Lundell & Lasič, 2020; Lundell, Nilsson, Dyrby, et al., 2019; Lundell et al., 2018). The residual diffusion variance due to SA in STE is shown for simple restricted diffusion compartments under GPA in Figure 4D. The analysis of TDD in tensor-valued encoding was presented by Lundell and Lasič (Lundell & Lasič, 2020), and a summary of results for the second cumulant, highlighting the effect of SA, is provided in the Supplementary Materials.

We define the spectral principal axis system (SPAS) as the orthogonal set of directions that exhibit maximally different sensitivities to TDD. Because TDD effects arise from the interaction between the encoding and diffusion spectrum, TDD sensitivity can only be defined in relation to the diffusion spectrum. This implies that SPAS does not represent a universal encoding property, as its definition is context dependent. Although SA can be fully described by considering all the spectral components of $su(ω)$ or by analyzing the b-tensor spectrum—that is, breaking down the b-tensor into contributions from different frequency channels (see figure 2.13B in Lundell & Lasič, 2020)—the SPAS offers an approximate, yet practical, approach for evaluating SA.

To gauge TDD differences for different spectral projections $su(ω)$⁠, we split the entire frequency range into three bands with equal encoding power. The crossover frequencies were determined based on the spectral trace $s(ω)$⁠, such that the total cumulative power,

reaches 1/3 and 2/3 of the total b-value (see Fig. 1B and Eq. (10), c.f. figure 2.15 in Lundell & Lasič, 2020). The spectral bands each contain one-third of the total encoding power, and their exact spectral content is not critical. The gradual transition of encoding power across directions (see Eq. (12)) ensures that when one band is prominent, the others contribute less. The SPAS can be defined ad hoc simply as the eigenvectors of the low-pass filtered b-tensor, that is, the tensor resulting from setting the integration limits in Eq. (9) as the crossover frequency of the low-frequency band. For waveforms with pronounced SA, this results in distinct spectral projections along the eigenvector axes of the low-pass-filtered b-tensor, capturing different TDD sensitivities, as shown in Figure 1B. Although RGB is used for both spectral anisotropy weighting and SPAS projections, this does not imply a one-to-one correspondence. SPAS projections are not band limited but span the full frequency range with overlapping spectra (Fig. 1C), each emphasizing different TDD sensitivities. In general, waveforms optimized for efficient b-weighting tend to concentrate encoding power at lower frequencies. The q-MAS serves as an illustrative example of an isotropic b-tensor, where $q(t)$ oscillates with high frequency orthogonal to a low-frequency (PGSE-like) axis (Lasič et al., 2014; Lundell, Nilsson, Dyrby, et al., 2019; Topgaard, 2013). The SPAS approach assumes spectral anisotropy in the low-frequency band for meaningful separation of TDD sensitivities. Experimentally feasible waveforms naturally exhibit sufficient spectral anisotropy to ensure distinct spectral projections. In principle, if a waveform exhibited a recessed mid-range along at least one projection—where both low- and high-frequency bands are strong—SPAS may misrepresent the directional variation of TDD sensitivity. The SPAS approach may not be relevant in cases with generally low SA, such as for multiple rotations of the q-vector (H. Jiang et al., 2023), where TDD sensitivity is more evenly distributed.

Alternatively, the SPAS could be defined in terms of the encoding spectral moments (see chapter 2.5.9 in Lundell & Lasič, 2020),

which for p = 1,2… corresponds to Eq. (13) in the low-frequency approximation of the diffusion spectrum, applicable when encoding comprises power at relatively low frequencies (Lundell & Lasič, 2020). The case of p = 1 could, for example, apply to extracellular hindered diffusion (Novikov et al., 2011, 2014; Novikov, Fieremans, et al., 2018), where $D(ω)∝|ω|$⁠, while the case of p = 2 could be suited for restricted diffusion, where $D(ω)∝ω2$ and $mij(2)$ projections represent spectral variance (Burcaw et al., 2015; Chakwizira, Westin, et al., 2023; Lasič et al., 2006, 2024; Lundell & Lasič, 2020; Nilsson et al., 2017; Stepišnik, 1993; Stepišnik et al., 2006).

The eigenvectors of $mij(p)$ could be used to provide alternative SPAS versions. In contrast to the low-pass filtered encoding spectrum used in Eqs. (20) and (21), the spectral moments involve high-pass filtering. However, the SPAS axes derived from the low-pass and high-pass filtered encoding spectra are expected to be nearly equivalent, provided that the high- and low-frequency axes of the SPAS are swapped (Supplementary Fig. S1). The equivalence or complementarity of the SPAS from such different approaches is a consequence of low- and high-frequency contributions summing up to the full b-tensor.

The centroid encoding frequency has been used in oscillating gradient experiments with LTE (Arbabi et al., 2020) and adopted for tensor-valued encoding (Narvaez et al., 2022), defined as

Similarly, the “projected” centroid frequency could be defined for b-tensor projections along vectors $u$ as

Eq. (23) could be applied along the SPAS axes to yield the corresponding “projected” centroid frequencies. Note that $ωuc$ are not defined along the $mij(p)$ eigenvectors with zero eigenvalues (see Eq. (21)). Importantly, for the SPAS, $ωu=SPASc$ provides a measure of the spread of encoding frequencies, which could be used to gauge SA. Reduced variation in $ωuc$ reflects reduced SA, where in the limit of spectrally isotropic encoding (SA = 0), the variation in $ωuc$ vanishes and $ωuc=ωc$⁠. SA generally decreases with increasing degrees of motion nulling, for example, by repeating encoding over several periods, resulting in no encoding at zero frequencies (Fig. 1). This is similar to the narrowing of the frequency window for oscillating gradients (Callaghan & Stepišnik, 1996) or using double rotation of the q-vector (H. Jiang et al., 2023). The limit case of SA = 0 is practically unattainable since infinite encoding frequencies would be required, corresponding to an infinite number of q-trajectory excursions that would eventually even out any directional differences in encoding power spectra.

The STE was obtained using the Matlab (The MathWorks, Natick, MA) package NOW (https://github.com/jsjol/NOW) for numerical optimization of tensor-valued encoding waveforms under given hardware constraints (Sjölund et al., 2015). Matlab code for generating gradient waveforms used in this work is available at https://github.com/samo-lasic/Lasic_SPAS_ImagingNeuroscience2025. To avoid concomitant field effects and obtain velocity compensation, we repeated the STE waveforms after the inversion RF pulse (see Fig. 1A; Szczepankiewicz, Westin, et al., 2021). The SPAS was obtained as eigenvectors of the low-frequency filtered STE as described in Section 2.4 and illustrated in Figure 1B. The low-frequency band was determined from the spectral trace yielding 1/3 of the total encoding power (b-value). The LTE projections corresponding to SPAS are shown in panels B and C of Figure 1 as SPAS1, SPAS2, and SPAS3. Signals from these waveforms were geometrically averaged to obtain the geoSPAS dataset (Section 3.4).

The tuned LTEs were obtained as described in Section 2.3.2 by taking 1D projections (LTE) of the 3D waveform (STE) along 1000 uniformly distributed directions, and finding the 1D waveform that yields apparent diffusivity $Du$ value most similar to MD from STE for spheres of radius R = 2.5 µm and bulk diffusivity D₀ = 2 µm²/ms. The MD was calculated by Eqs. (11) and (15), and $Du$ by Eqs. (12) and (18). From 10% of best tuned LTE projections shown in Figure 2, the one minimizing the difference $|Du−MD|$ was chosen as the tuned LTE (tLTE), while the “optimized tuned” LTE (optLTE) (Lasič et al., 2023)) was chosen as the one with the lowest demand on the gradient system, defined as the one minimizing the product of maximum gradient magnitude and maximum gradient slew rate, $max{|gu|}⋅max{|g˙u|}$⁠.

The gradient waveforms were played in two equal 21 ms long encoding blocks separated by a 5 ms gap for the refocusing RF pulse. The waveforms were copied on both sides of the refocusing pulse, yielding velocity-compensated effective gradient waveforms of total duration 47 ms as shown in Figures 1C and 2B. Full velocity compensation ensures no encoding of incoherent ballistic flow, which can otherwise confound microstructural assessments due to pseudo-diffusion in vivo (Ahlgren et al., 2016). The centroid frequencies were $ωijc​/2π$ = 45 Hz for the STE and in the range 21–69 Hz for SPAS-LTEs. The exchange weighting times (Ning et al., 2018) were in the range Γ = 1.6–1.8 ms.

Simulations were implemented in Matlab (The MathWorks, Natick, MA). Code is available at https://github.com/samo-lasic/Lasic_SPAS_ImagingNeuroscience2025. Apparent diffusion coefficients (ADC) were calculated according to Eqs. (13) and (14) using expressions for diffusion restricted in spheres, cylinders, and planes (Lundell & Lasič, 2020; Stepišnik, 1993) with bulk diffusivity of D₀ = 2 µm²/ms.

Axisymmetric spheroids were modeled by the diffusion spectrum for spheres with two different restriction sizes. This approximation has been previously used and validated by simulations (Lasič et al., 2022; Lundell, Nilsson, Dyrby, et al., 2019; Nielsen et al., 2018). Diffusion in sticks of finite length and zero radial diffusivity was modeled using the diffusion spectrum for planar geometry. ADCs were obtained for 300 uniformly distributed rotations and average signals resulted from the arithmetic average of the corresponding 300 monoexponentially attenuated signals (Fig. 4A). Note that rotating the substrate is equivalent to rotating the encoding. This was performed for all 6 waveforms and for 100 restriction radii linearly spaced in the range R = 1–20 µm.

All experiments were approved by the Danish Animal Experiments Inspectorate (2018-15-0201-01551) following the European Communities Council Directive (2010/63/EU). Ex vivo and in vivo rat brains as well as phantom data were acquired on a Bruker BioSpec 7T (Bruker BioSpin, Ettlingen, Germany) MR Scanner with 0.66 T/m gradients. The individual protocols are described below.

Perfusion-fixed Sprague-Dawley adult rat brains stored at 4°C in paraformaldehyde (PFA) were transferred to a glass tube filled with a fluorocarbon fluid (Fluorinert, Sigma-Aldrich, Germany). After 2 to 3 h at room temperature (22 ± 1°C), the sample was secured in the Bruker mouse cryoprobe cradle and positioned at the center of the magnet. The sample temperature was stabilized at 27°C using both the surface temperature sensor control of the cryo-coil and a heated airflow with an additional temperature probe attached underneath the sample outside the field of view (SA Instruments, New York, USA). Temperature deviations recorded with the additional probe were within 0.1°C during the experiment. A mouse ¹H quadrature T/R cryoprobe surface array coil was used for measurements (Bruker, Germany).

One Sprague-Dawley female rat (n = 1; 180 g) was continuously anesthetized using 2–3% isoflurane in air and oxygen. The rat head was secured in a stereotactic frame using a bite bar and ear bars. Respiratory rate was monitored through a pillow placed underneath the rat’s body, and temperature was maintained at 37.5 ± 0.5°C using a rectal probe connected to warmed air blower (SA Instruments, New York, USA). An 86-mm diameter volume coil and 20 mm diameter single-loop surface coil (Bruker, Ettlingen, Germany) were used for transmission and reception, respectively.

A 1 cm plastic tube filled with demineralized water was used to test the calibration of the gradients for all the encoding waveforms and imaged with the same setup as for the ex vivo sample.

A multi-shot echo planar imaging sequence customized for diffusion encoding with general gradient waveforms (https://osf.io/t9vqn/) was used. Both in vivo and ex vivo scans used TR/TE = 3 s/53.136 ms, matrix size = 96 x 96 and 6 axial slices of 1 mm thickness. Adjustments included acquisition of a B0-field map and global shimming with MAPSHIM, along with EPI-ghost and trajectory corrections. Diffusion encodings were performed using 5 waveforms (STE, SPAS1, SPAS2, SPAS3, tLTE) with 12 rotations corresponding to the vertices of the icosahedron.

Ex vivo scans used 16 segments per image, 6 averages, FOV = 16 x 16 mm², in-plane resolution = 0.17 mm², scan time = 48.4 h, 10 logarithmically spaced b-values (123, 168, 257, 366, 550, 805, 1194, 1760, 2560, 3806 s/mm²).

In vivo scans used 4 segments per image, 3 averages, FOV = 25 x 25 mm², in-plane resolution = 0.26 mm², scan time = 3.65 h, 6 logarithmically spaced b-values (123, 238, 466, 952, 919, 3806 s/mm²).

Preprocessing steps included Marchenko–Pastur denoising implemented in Matlab (https://github.com/sunenj/MP-PCA-Denoising) (Does et al., 2019; Veraart et al., 2016) and extrapolation-based motion and eddy-current correction (Nilsson et al., 2015). Data analysis was performed with Matlab using a combination of custom scripts (https://github.com/samo-lasic/Lasic_SPAS_ImagingNeuroscience2025) and parts of the multidimensional diffusion MRI framework (https://github.com/markus-nilsson/md-dmri). Smoothing was applied per slice with 2D Gaussian smoothing kernel with standard deviation of 0.5. Signal attenuations from SPAS-LTEs were geometrically averaged, yielding the geoSPAS dataset. To estimate diffusion tensors, a subset of geoSPAS data acquired with 12 waveform directions and b-values in the range 123–805 s/mm² were used. Nonlinear least squares fitting was applied with symmetric positive-definite constraints to yield maps of signal without diffusion weighting (S₀), mean diffusivity (MD), and fractional anisotropy (FA). To generate maps of microscopic fractional anisotropy (µFA), data from STE and geoSPAS were averaged across 12 directions, and constrained nonlinear least squares fitting was applied using the gamma function as described previously (Lasič et al., 2014).

Model-free contrast maps associated with microscopic anisotropy (µA) and time-dependent diffusion (TDD) were obtained by subtraction of signals from different encodings (see Eq. (3)). Logarithms of direction averaged signals (⁠$log[〈S〉]$⁠) were taken at the maximum b-value of 3806 s/mm². The differences (⁠$Δlog[〈S〉]$⁠) between geoSPAS and STE were used to generate the µA contrast maps according to

and the differences between SPAS1 and SPAS3 were used for the TDD contrast maps according to

Subtracting logarithms of signals instead of raw signals yields contrast related to signal attenuation, effectively removing the b0 signal contribution. Regions of interests (ROIs) where drawn manually by reference to a rat brain atlas (Paxinos & Watson, 1996) in the cortical gray matter (CGM), choroid plexus (CP), dentate gyrus of the hippocampus (DGHC), and white matter (WM). Additional ROIs were generated by segmentation based on the ranges of µA and TDD contrasts (Supplementary Fig. S4). ROI-average signals versus b-value were normalized using S₀ extrapolation from the Gamma distribution fit of the multiexponential signal attenuations (Lasič et al., 2014). Full signal representations in the presence of restricted diffusion are derived in the Supplementary Materials.

STE waveforms (gradient, dephasing) and the spectral trace are shown in Figure 1A. The nonexistent zero-frequency component shows that the encoding is velocity compensated. The rotation relative to the laboratory frame (XYZ axes) is arbitrary and depends on the output from the optimization (Sjölund et al., 2015). Thus, using any of the XYZ projections for tuning purposes is generally not applicable, even though such approach can be justified as a proxy for tuning when specific q-trajectories are used (Lundell, Nilsson, Dyrby, et al., 2019). Spectral anisotropy of the STE is shown in Figure 1B, color coded with red, green, and blue, based on the relative amounts of encoding power from the low-, mid-, and high-frequency bands for each spectral projection. The crossovers between the frequency bands were determined from spectral trace yielding 1/3 and 2/3 of total encoding power (b-value), as depicted in panel B. The SPAS, shown in Figure 1B as red, green, and blue lines, corresponding to the eigenvectors of the low-pass filtered b-tensor, maximizes the spread of encoding power between the low-, mid-, and high-frequency bands. This is shown in Figure 1B and C by comparing the encoding power spectra for the SPAS1-SPAS3 LTE waveforms.

The tuning is shown in Figure 1B as contours outlining 10% of all LTE projections with smallest values of $|Du−MD|$ for spheres of R = 2.5 µm (solid lines) and R = 5 µm (dashed lines). For spheres of R = 2.5 µm, tuning is also shown in Figure 2A as the color-coded ratio $Du​/MD$⁠. The tLTE refers to the STE waveform projection with $Du=MD$⁠. As an alternative to the tLTE, the “optimized tuned” LTE (optLTE) is included in Figure 2. In this case, the tuning was slightly compromised to reduce the maximum gradient magnitude from 652 mT/m to 583 mT/m and slew rate from 2.6 T/ms/m to 1.3 T/ms/m, respectively. Due to the smooth nature of the encoding and diffusion spectra, the goal of matching MDs is equivalent to matching the envelopes of encoding spectra from LTE and the spectral trace from STE. The robustness of such tuning approach can be appreciated by comparing the tuning contours for two restriction sizes (solid and dashed lines), shown in Figure 1B, and noting that the tuned LTE projection shown in Figure 2A approximately corresponds to the intersection between the two contour lines. The tuning landscapes shown in Figure 3 illustrate how tuning varies with spherical angle across different sphere sizes, showing how closely the relative tuning approaches the ideal value of 1. For larger spheres, tuning remains close to 1 across the landscape, reflecting reduced sensitivity to TDD as diffusion approaches the Gaussian regime. Smaller spheres exhibit greater modulation, with sharper contours, while for very small structures, tuning becomes less relevant as both MD and $Du$ approach 0 (c.f. Fig. 4). Notably, the position of the tuning contours is largely independent of sphere size, indicating robust tuning across a wide range of restriction sizes.

The competing contrasts due to µA and TDD were simulated for diffusion within powders of cylinders, spheroids, and sticks of finite length (see Fig. 4). Signals versus b-value are shown for selected restriction sizes, for which both contrasts are visible. The encoding power gradually shifts from lower to higher frequencies for SPAS1 to SPAS3 encodings (Fig. 1), resulting in decreasing signals (Fig. 4A) and increasing ADC values (Fig. 4C) due to TDD. The tLTE and optLTE yield similar signals, which are between SPAS1 and SPAS3 and approximately match the geometrically averaged SPAS signals (geoSPAS). Note that the geoSPAS signals represent ideal matching of spectral traces and thus ideal tuning. While signal attenuations from LTEs are multi-exponential due to orientation dispersion, signals from STE are approximately mono-exponential due to mono-dispersed restriction size. The initial slope and the deviation from mono-exponential signal attenuation are related to the mean apparent diffusivity and the diffusion variance, respectively, shown in panels C and D. The signal difference between geoSPAS and STE at high b-values indicates the effect of µA unconfounded by restricted diffusion. The size dependence of TDD and µA is shown in Figure 4B. For cylinders, the µA contrast is maximized for small sizes, when diffusion appears Gaussian, similar as for example in hexagonal liquid crystals (Lundell, Nilsson, Dyrby, et al., 2019). With increasing cylinder radius, the radial ADC increases and the µA contrast decreases. The TDD contrast peaks at intermediate sizes, depending on encoding times, and fades away for large sizes when diffusion again appears Gaussian (compare panels B and C). For sticks, the size dependence of the TDD contrast is similar as in the case of cylinders, but the µA contrast exhibits the opposite trend. For shorter sticks, the axial ADC decreases toward the radial ADC of zero, while it is maximized for large sizes, leading to maximum µA contrast for longer sticks. A more complex behavior can be seen for the spheroids. As in the case of cylinders, both contrasts decrease for large sizes as radial and axial ADCs approach bulk diffusivity (⁠$D0$⁠). The TDD contrast is dominating due to the low compartment anisotropy and restriction along all axes (ratio 1:3). Two different restriction lengths result in a wider TDD peak compared with the cylinders. The µA contrast peaks at intermediate sizes and vanishes for small and large sizes, when diffusion appears Gaussian and ADCs are, respectively, zero and bulk (⁠$D0$⁠). SA in STE leads to small deviations from mono-exponential decays. This is reflected also in Figure 4D as non-zero values of the diffusion variance $VD$ (black solid lines) at intermediate restriction sizes for all three substrates. However, the maximum $VD$ from STE was only about 3% of the maximum $VD$ from LTE across all the substrates. The effects of SA for restricted diffusion are analyzed in more detail in the Supplementary Materials.

Reference experiments were conducted on water, confirming that the signals from all encodings overlap at equal b-values (Supplementary Figs. S2 and S3). Ex vivo results are shown in Figures 5 to 8 and in vivo results are shown in Figures 9 and 10. Note that the ex vivo results were obtained with a cryo-coil optimized for the mouse brain. This setup maximized overall SNR, but the coil’s sensitivity profile tapers off in deeper brain regions.

Signal attenuations from geometrically averaged SPAS LTEs (geoSPAS) yield S₀, MD, and FA parameter maps shown in Figure 5A–C. Rotationally averaged geoSPAS and STE data yield µFA maps unconfounded by restricted diffusion (Fig. 5D). The µFA² maps approximately correspond to the µA-based contrast maps shown in Figure 5F, obtained from high b-value geoSPAS and STE signals as a subtraction map (Eq. (24)). The protocol unfolds the µA contrast from the TDD contrast, which is shown in Figure 5E as a subtraction map of high b-value SPAS1 and SPAS3 signals (Eq. (25)). The relationship between TDD and µA contrasts is shown in the color-coded TDD-µA joint contrast maps in Figure 6. A closer look at the regions with cortical gray matter (CGM), choroid plexus (CP), dentate gyrus of the hippocampus (DGHC), and the white matter (WM) is shown in Figure 7. The relationship between µA and TDD contrasts can be visualized as a voxel-wise scatter plot in Figure 7A for the different regions outlined in Figure 7B. Overall, the µA contrast is pronounced in the white matter and TDD contrast is higher in the CP and in the DGHC. The normalized direction and ROI-average signals are shown in Figure 7C. The normalized ROI signals consistently decrease for SPAS1 to SPAS3 as expected from theoretical predictions in the presence of restricted diffusion (Fig. 4). Importantly, the signals from the tuned LTE (tLTE) coincide well with the signals from the geoSPAS. In CGM, both contrasts are relatively low. However, an interesting heterogeneity is visible within CGM with µFA exhibiting radial patterns (Fig. 5D) and TDD showing a bright tangential layer (Fig. 5F).

In WM, µA is large and TDD is relatively low but larger than in CGM. The CP and the DGHC exhibit intermediate µA and elevated TDD contrast, with pronounced TDD in the CP. The STE is close to mono-exponential in the CGM and WM, and more multi-exponential in the CP and DGHC, indicating larger tissue heterogeneity in those regions. The pronounced deviation from mono-exponential STE signals in the CP is reflected also by the larger spread of the TDD contrast in this region. The LTE signals in the CP resemble the simulated scenario with spheroidal restrictions (Fig. 4A), where TDD is pronounced due to restrictions along all three axes, which could be related to the high density of cuboidal epithelial cells, fibroblasts, or macrophages (Saunders et al., 2023).

The effects of TDD on the µFA parameter maps from the fixed tissue when using the geoSPAS signals (geometric averaging) or not tuned SPAS LTE signals are shown in Figure 8. This figure demonstrates the extent of µFA bias that arises when waveforms with different TDD sensitivities are used, highlighting the importance of tuning to achieve unbiased µFA estimation. The geoSPAS yields the “correct” µFA map not confounded by restricted diffusion. In comparison, using SPAS3 and SPAS1 data results in reduced and elevated µFA, respectively. Note that the fitting model used to estimate µFA assumes that STE and LTE yield equal MD values, an assumption that was not fulfilled for the maps in Figure 8A and C. The decreased/increased µFA for SPAS3/SPAS1 is due to increased LTE signal differences relative to the STE signal. The µFA difference maps reflect the effects of TDD (D-E). The fitting noise propagates from µFA with SPAS3 to the differences in D. The µFA difference from SPAS1 and geoSPAS (E) is more accurate, and this map resembles the TDD-based contrast maps in Figure 5F obtained by subtracting logarithms of direction-averaged signals. However, additional contrast is visible on the µFA difference (∆µFA²) in Figure 8E (slice 2) between the DGHC and the WM, and the layers in CGM are more pronounced than in Figure 5F. Additional µFA parameter maps can be found in the Supplementary Figure S5.

In vivo results are shown in Figures 9 and 10. Compared with the fixed tissue, MD was higher in vivo while FA and µFA remained in similar ranges (Fig. 9). Notably, the µA-based contrast was more pronounced in vivo (Fig. 9F), whereas the TDD-based contrast was generally lower (Fig. 9E). The relationship between both contrasts is shown in the TDD-µA joint contrast map in Supplementary Figure S6, corresponding to the fixed tissue case (Fig. 6). Similar µA-TDD contrasts were found in vivo for the CGM and WM, but relatively lower µA in the DGHC (Fig. 10). Due to the larger coverage of the coil in the in vivo setup, the cerebellar cortex (CBC) was visible in vivo, where relatively low µA and the highest TDD contrasts were observed. The generally lower TDD contrast in vivo is reflected also on the µFA parameter maps and their differences for SPAS LTEs (Supplementary Fig. S7), which were less pronounced compared with the fixed tissue (Fig. 8; Supplementary Fig. S5).

This work considers the effects of TDD in tensor-valued diffusion encoding, addressing the need for tuning gradient waveforms for equal sensitivity to restricted diffusion. The need for a systematic approach to tuning in tensor-valued encoding was first highlighted in Lundell, Nilsson, Dyrby, et al. (2019). The theoretical work presented in Lundell and Lasič (2020) identified the spectral trace as the “tuning lever” determining MD. The spectral trace allowed estimating MDs from direction-averaged signals for a wide range of encodings with STE and LTE waveforms with varying degrees of motion compensation and restriction encoding, enabling estimation of cell sizes in heart tissues (Lasič et al., 2022). However, previous work did not propose general strategies for tuning different b-tensors or identifying waveform projections with maximally different restriction encoding.

We have proposed two strategies for tuning b-tensors of different shapes and thereby obtain unbiased assessment of µA in the presentence of TDD. By tuning, we refer to the matching of apparent mean diffusivities (see Eq. (15)). One b-tensor is considered as a tuning reference and the other as the tensor being tuned to the reference. We remark that in the presence of restricted diffusion, µA may depend on the encoding timing parameters (Ianuş et al., 2017; Lasič et al., 2020; Lundell et al., 2021; Teh et al., 2021). However, the first order effects of restricted diffusion, as governed by tuning, remain constant between LTE and STE (Lasič et al., 2022). We have demonstrated the case of tuning LTE to STE (Figs. 1 and 2). Furthermore, we have introduced the concept of spectral anisotropy (SA) and defined the associated spectral principal axis system (SPAS) of encoding tensors. We have demonstrated how the SPAS can be used to generate independent contrast maps due to µA and TDD within a single multidimensional diffusion-encoding protocol.

The first tuning approach selects waveform projections such that the apparent mean diffusivity (MD) matches that of the reference encoding tensor (tLTE in Fig. 2). This approach provides a relative tuning tailored to different diffusion spectra. The smooth nature of the diffusion spectrum and smooth variation of encoding power distribution from tensor projections warrant robustness of this approach, reflected by the rather constant peak positions on the tuning landscapes for restricted diffusion in a wide range of sizes (Fig. 3). This robustness holds even for weaker diffusion dispersion, for example, in white matter, where the diffusion spectrum is flatter at low frequencies. In such cases, the tuning contours exhibit less directional variation, but the tuning remains determined by the stable position of the contours rather than their depth. While absolute tuning changes with the diffusion spectrum, the key factor is the relative tuning between the tensors. This is akin to tuning an orchestra, where the absolute pitch may shift, but the instruments remain tuned relative to each other. Importantly, tuning with spherical geometry can be generalized to other shapes, as long as the isotropic diffusion spectrum $λ¯iso(ω)$ in Eq. (15) approximates the diffusion spectrum in spheres of appropriate size. This tuning approach can also be applied to non-STE waveforms, as the tuning is determined by the encoding spectral trace $s(ω)$⁠. Furthermore, it is independent of the gradient moments and can, therefore, be applied to waveforms with arbitrary flow compensation properties. The flexibility of this approach allows for a degree of optimization, where some tuning could be sacrificed to better suit hardware constraints (optLTE vs. tLTE in Fig. 2). We speculate that this approach of deriving the tLTE could also be tailored for incoherent flow, which is characterized by an exponential velocity autocorrelation function (Kennan et al., 1994) and thus by a Lorentzian diffusion spectrum. Such adaptation could be useful for probing time-dependent anisotropic perfusion.

The second tuning approach aims to match the spectral traces of different encoding tensors precisely, making the tuning independent of the diffusion spectrum. This is achieved through geometric averaging of signals (geoSPAS) from any orthogonal set of encoding projections (Eqs. (15) and (19)). The SPAS provides a unique orthogonal set of LTE projections (SPAS1–SPAS3) with maximally different sensitivities to restricted diffusion, inherent in a spectrally anisotropic reference tensor (STE). In this context, a high SA of the encoding is valuable for distinguishing TDD contrast alongside µA contrast or the µFA map.

The competing effects of µA and TDD contrasts in different restricted geometries are illustrated in Figure 4. Signals decrease from SPAS1 to SPAS3 encodings, while the tuned encodings (tLTE and optLTE) yield intermediate signals approximately matching the geoSPAS signals. At higher b-values (2^nd cumulant), the difference between STE and geoSPAS reflects µA, while the difference between the SPAS-LTEs reflects TDD. The simulations presented in Figure 4 and the previous Monte-Carlo simulations (see supplementary information in Lundell, Nilsson, Dyrby, et al., 2019) suggest that SA in STE has a relatively small effect on the second cumulant, while tuning (first cumulant) is crucial for accurate µA assessment. As previously observed, the µA contrast depends on restriction scale relative to encoding time (Ianuş et al., 2017; Lasič et al., 2020; Lundell et al., 2021; Teh et al., 2021). However, it is crucial that LTE is tuned to STE for TDD to not confound the µA assessment (Lasič et al., 2022). Further analysis of SA and restricted diffusion effects in tensor-valued encoding is available in the Supplementary Materials section.

The multidimensional SPAS protocol yielded parameter maps (MD, FA, µFA) as well as µA- and TDD-based contrast maps via signal subtraction for the fixed tissue (Fig. 5) and in vivo (Fig. 9). While the parameter maps were generated based on higher-order cumulants from the joint fit of STE and LTE data, the µA- and TDD-based contrast maps were derived in a model-free manner by direct signal subtraction at high b-values. Although these results demonstrate the feasibility of generating µA- and TDD-based contrasts in both ex vivo and in vivo settings, it remains to be determined how reproducible these contrasts are across different samples and experimental conditions. It was crucial that the µFA parameter maps and the µA-based contrast maps use the geoSPAS signals along with the STE signals to ensure tuning in terms of equal MD values and thus reduce the confounding effects of TDD. While the geoSPAS signal only accounts for first order attenuation terms, it overlaps well with the tLTE signals (Figs. 7 and 10), suggesting that higher order terms are similarly represented by the different LTE waveforms. The µA- and TDD-based contrast maps display two independent contrasts, which can be combined to create color-coded TDD-µA joint contrast maps (Fig. 6; Supplementary Fig. S6), or used as a model-free approach for tissue segmentation (Supplementary Fig. S4). These joint contrast maps provide an integrative view of spatial relationships between cell size (TDD) and anisotropy (µA), facilitating identification of tissue regions with unique microstructural properties. Because the MD values from geoSPAS and STE are equal, the µA-based contrast maps, generated by subtracting $log[〈S〉]$⁠, are proportional to the differences in diffusion variance (or “mean kurtosis”). In contrast, the TDD-based contrast maps, created by subtracting the signals from SPAS1 and SPAS3 encodings, reflect both MD and diffusion variance differences due to TDD (see Eq. (3)). It should be noted that the µA and TDD maps are generated exclusively from images with the highest b-value. These maps could easily be improved by averaging more high-b data, rather than acquiring additional low-b data. In addition, using only high b-values helps reduce artifacts such as CSF partial volume effects and Gibbs ringing, contributing to cleaner contrast.

The importance of tuning for unbiased µFA maps is illustrated in Figure 8 and in Supplementary Figures S5 and S7. Note that just as MD, also µFA is “apparent” and depends on the TDD, but tuning of tensor-valued encoding is crucial for unbiased µFA estimation and separation of the isotropic and anisotropic sources of diffusion variance. Without this consideration, TDD could lead to severely biased results and interpretations particularly in preclinical diffusion-weighted imaging experiments with strong gradients that combine tensor-valued encodings such as STE and LTE with very different sensitivities to restricted diffusion (He et al., 2021). With the SPAS protocol, the “corrupted” µFA maps could be differentiated and thus potentially provide an additional valuable contrast (Fig. 8D, E; Supplementary Figs. S5D–F and S7D–F).

Elevated TDD has been observed with oscillating gradients in the ex vivo mouse brain, particularly in the granular layers of the cerebellum and hippocampus (Aggarwal et al., 2012, 2020). A similar trend has been observed in the cerebellum of the ex vivo monkey brain (Lundell et al., 2015; Lundell, Nilsson, Dyrby, et al., 2019). OGSE experiments in mice have shown similar results, with higher MD ex vivo and reduced TDD effects in vivo, consistent with our observations (Wu et al., 2014). The underlying mechanisms for the in vivo versus ex vivo differences may be related to neuronal beadings—a morphological feature suggested to drive diffusion changes during hypoxia and stroke (Budde & Frank, 2010). Interestingly, the decreased MD and increased TDD in fixed tissue resemble patterns observed in an in vivo stroke model, where TDD changes were recently demonstrated using tensor-valued encoding with differently tuned LTE waveforms (Zhou et al., 2024). This suggests that hypoxia and fixation may induce similar microstructural changes. Additionally, TDD effects have been shown to vary across cortical regions, such as the differences between the soma-rich occipital cortex and frontal regions (Lundell, Nilsson, Dyrby, et al., 2019). In this work, we also observed elevated TDD in specific cellular layers (Fig. 5). The comparison of the choroid plexus signals and simulation with spheroidal restrictions suggest a potentially high density of cuboidal epithelial cells, fibroblasts, or macrophages (Saunders et al., 2023). We speculate that the radial patterns visible within cortical gray matter on the µFA maps (Fig. 5D) could be related to cortical columns with high neurite densities or are associated with blood vessels, as suggested by the radial patterns observed in the cortex on the MD map (Fig. 5B). A bright tangential layer on the TDD map (Fig. 5F) possibly corresponds to the soma-rich layer 4 (Narayanan et al., 2017) with cell radii ~5 µm (Meitzen et al., 2011).

Alternative approaches to probe microscopic anisotropy rely on double diffusion encoding, which is a special case of employing planar tensor encoding (PTE) along with LTE (Jespersen et al., 2013). Tuning between LTE and PTE is in this case ensured by using long mixing times between the two encoding blocks, assuming equal timing is used in the two blocks. With double diffusion encoding, long mixing times could in principle also allow to independently vary tuning and SA for additional specificity to restricted anisotropic diffusion (Lundell et al., 2018). A related approach employing double oscillating gradients can be used to map time-dependent µA (Ianuş et al., 2017, 2018; Nielsen et al., 2018). However, such approaches rely on specific encoding protocols with limited degrees for optimizations, and importantly, cannot assess the relationship between µA and TDD within a single multidimensional protocol with minimal confounding factors.

Some degree of SA is generally unavoidable in tensor-valued encoding, which uses asynchronous gradient waveforms along orthogonal axes. This needs to be considered in tuning. Due to SA, STE may no longer be rotationally invariant (isotropic encoding), but depend on the relative angle between the encoding and the main diffusion axis (Jespersen et al., 2019; Lasič et al., 2021; Lundell et al., 2018; Lundell, Nilsson, Dyrby, et al., 2019; Lundell & Lasič, 2020). While this effect has not been observed in the human brain in vivo with encoding times of about 75 ms (Szczepankiewicz, Lasič, et al., 2019), it has been demonstrated in a fixed monkey brain on a preclinical scanner with encoding times of about 25 ms (Lasič et al., 2021) and in prostate cancer mice models with encoding times of about 15 ms (Szczepankiewicz et al., 2023). When TDD is significant, multiple rotations of the STE are required for the powder-averaged signal to be truly isotropic.

Identifying the SPAS is instrumental for taking SA into consideration. In addition to the low-pass filtering approach for defining the SPAS, we proposed alternative SPAS definitions based on encoding spectral moments (see Section 2.4 and Supplementary Materials). The spectral moments also relate to the generalization of centroid encoding frequencies from LTE (Arbabi et al., 2020) to tensor-valued encoding, which may be useful for gauging SA.

The proposed approach does not aim to maximize sensitivity to TDD per se, but rather exploits the available SA of the reference encoding. While we have here used the SPAS to obtain encodings with maximally different sensitivities to TDD, one could instead minimize such differences by intermediate rotations of the SPAS, thus providing yet another means of tuning, where a single projection could be used instead of employing geometric averaging of signals. Such approach resembles an approximate tuning used previously, where a waveform along a single axis of STE could be used for an approximately tuned LTE (Lundell, Nilsson, Dyrby, et al., 2019). We have since identified the spectral trace as the key property for tuning b-tensors (Lundell et al., 2021), which was instrumental for the herein proposed tuning approaches. Our analysis elucidates the role of SA in the geometric averaging approach to tuning. It also suggests that using reference tensors like STE with axially symmetric SA, resembling the q-MAS trajectory (Topgaard, 2013) used in Lundell, Nilsson, Dyrby, et al. (2019), could enable faster acquisitions of “tuned” signals via weighted geometric averaging of only two instead of three SPAS-LTEs.

An important aspect of our tuning approach is the use of b-tensors that are derived from the reference b-tensor, inheriting q-trajectory properties, which can be optimized for different hardware constraints (Sjölund et al., 2015) including gradient moment nulling (Lasič et al., 2020; Szczepankiewicz, Teh, et al., 2021) or minimization of concomitant field effects (Szczepankiewicz, Westin, et al., 2019). While it is tempting to consider additional LTEs designed to boost TDD sensitivity instead of employing derived LTEs, such an approach would require a separate protocol, potentially prolong scan time, and potentially compromise direct comparison with µA, if the average restriction encoding is not matched. We expect that the proposed SPAS-based approach could be adopted for any pair of b-tensors provided that the tensors being tuned are derived from the reference tensors via affine transformations, that is. by scaling the original q-trajectories.

There are some limitations of the proposed approach. The method relies on the spectral-domain analysis in the Gaussian phase approximation, neglecting effects of intra-compartmental kurtosis (Henriques et al., 2020; Jespersen et al., 2019) and exchange (Chakwizira, Westin, et al., 2023; Chakwizira, Zhu, et al., 2023; Lasič et al., 2024; Olesen et al., 2022; Pfeuffer et al., 1998). The accuracy of such approach would need to be tested with Monte-Carlo simulations for different substrates (Chakwizira, Westin, et al., 2023; Lundell, Nilsson, Dyrby, et al., 2019; Lundell, Nilsson, Szczepankiewicz, et al., 2019). Only simple restricted diffusion models for monodispersed compartments with a single size and shape were used in simulations and for deriving tuned projections. Additional simulations could be used to address precision and optimize the protocol in various experimental settings. Biophysical models comprising heterogeneous restricted compartments with exchange (Chakwizira, Westin, et al., 2023; Nilsson et al., 2018), hindered extracellular diffusion spaces (Novikov et al., 2011, 2014; Novikov, Fieremans, et al., 2018), or time-dependent perfusion (pseudo-diffusion) (Kennan et al., 1994) could be considered in future studies. Depending on the nature of TDD effects, modulating the first order effects, as in our approach, or the higher order terms such as kurtosis could have complementary value. Although the proposed method does not account for exchange, we note that in our experiments, the exchange weighting times (Ning et al., 2018) were relatively short and within a narrow range (Γ = 1.6-1.8 ms), approximately 20 times shorter than those used in studies specifically targeting exchange and TDD (Chakwizira, Westin, et al., 2023; Chakwizira, Zhu, et al., 2023). Since the tuned encoding inherits encoding features from the reference, potentially confounding effects of exchange from different tensor projections could be minimized by striving for reference b-tensors with minimal directional variation in exchange sensitivity.

This work is relevant to several studies of microscopic anisotropy employing tensor-valued encoding, including q-trajectory imaging and diffusional variance decomposition, applied, for example, on brain tumors (Brabec et al., 2022; Li et al., 2021; Nilsson et al., 2018; Szczepankiewicz et al., 2016), various neurological conditions (Afzali et al., 2022; Andersen et al., 2020; Lampinen, Zampeli, et al., 2020; Syed Nasser et al., 2022), and in various organs such as kidney (Nery et al., 2019), prostate (Langbein et al., 2021; Nilsson et al., 2021), breast (Cho et al., 2022), and heart (Teh et al., 2023). Considering b-tensor tuning and TDD is particularly necessary in preclinical settings employing shorter encoding times (He et al., 2021; Rios-Carrillo et al., 2023; Yon et al., 2020). Our findings support the potential of this method for studying microstructural tissue changes in various animal disease models, including conditions such as inflammation, stroke, and tumors. Recently, the ability of tensor-valued encoding to probe TDD in the human brain has been demonstrated (Johnson, Irfanoglu, et al., 2023; Johnson, Ross, et al., 2023), employing centroid encoding frequencies within the range of 6.5–21 Hz, reaching a maximum of only 11 Hz at b = 3000 s/mm². The encoding frequencies employed on the human scanner were notably lower than the frequencies in our study, which ranged from 21 to 69 Hz. This observation underscores the feasibility of applying SPAS-based protocols in human studies, not only for mapping microscopic anisotropy but also for detecting TDD. This is particularly relevant in conjunction with the emerging specialized gradient coils, such as MAGNUS (Foo et al., 2020), which offer exciting possibilities for microstructure imaging. The proposed method could help unveil microstructural features and play an important role in refining biophysical models of white matter (Assaf et al., 2004, 2008; Coelho et al., 2022; Jelescu et al., 2016; Nilsson et al., 2013; Novikov, Fieremans, et al., 2018), gray matter (Afzali et al., 2020; Jelescu et al., 2022; Novikov et al., 2011; Olesen et al., 2022; Palombo et al., 2020), or tumor tissues (X. Jiang et al., 2017; Nilsson et al., 2018; Reynaud, 2017; Reynaud et al., 2016b, 2016a).

In this work, we addressed the confounding effects of time-dependent diffusion (TDD) due to restricted diffusion in tensor-valued encoding and proposed strategies for tuning b-tensors to enable an unbiased assessment of microscopic anisotropy (µA). We introduced two approaches for tuning b-tensors to minimize the confounding effects of restricted diffusion. The first approach involves identifying encoding tensor projections that yield equal mean diffusivities, ensuring reliable tuning across various diffusion spectra for unbiased µA assessments. The second approach utilizes geometric averaging of signals from multiple LTEs, providing more robust tuning that is independent of the underlying diffusion spectra.

Additionally, the same encoding waveforms used in the geometric averaging approach for probing µA can also be used to probe an independent contrast due to TDD. This is facilitated by introducing the concepts of spectral anisotropy (SA) and the spectral principal axis system (SPAS) in tensor-valued encoding. Projections along the SPAS axes generate encodings with maximally different sensitivities to TDD, enabling the creation of two independent contrasts—µA and TDD—within a single multidimensional diffusion encoding protocol.

In the proposed framework, SA is not seen as a limitation that breaks the rotational invariance of STE, but rather as a valuable encoding feature that provides an additional independent contrast. The approach is flexible, accommodating free gradient waveforms that can be optimized to suit various hardware constraints. SPAS encodings provide complementary information to the b-tensor anisotropy, which is necessary to study orientationally dispersed, time-dependent anisotropic diffusion.

A minimal experimental protocol using SPAS encodings at a single b-value could serve as a model-free method to generate contrast maps specific to µA and TDD based on signal differences. The methodology shows potential for studying tissue microstructure changes in various disease models and may enable detection of both µA and TDD in the human brain. These methodological advances could refine biophysical models and enhance our understanding of tissue microstructure.

Code used for data analysis and simulations is available on GitHub at https://github.com/samo-lasic/Lasic_SPAS_ImagingNeuroscience2025.

Samo Lasič and Henrik Lundell: Conceptualization, Methodology, Simulation Software, Writing—Original Draft, Visualization. Samo Lasič, Markus Nilsson, and Nathalie Just: Data Analysis Software. Nathalie Just, Henrik Lundell, and Matthew Budde: Experiments. Samo Lasič, Nathalie Just, Markus Nilsson, Filip Szczepankiewicz, Matthew Budde, and Henrik Lundell: Writing—Review & Editing.

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 804746), VR (Swedish Research Council) (grant numbers 2020–04549, 2021-04844), Swedish Cancer Foundation (grant numbers 2022/2414, 22 0592 JIA, and 22 2011 Pj), eSSENCE (grant number 10:5), NIH (National Institutes of Health) (grant numbers R01NS109090, R01NS125781, R01MH074794, and P41EB015902).

S.L., M.N., F.S., and H.L. are inventors of related IP and may financially benefit from its use. All other authors declare that the research was conducted without any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary material for this article is available with the online version here: https://doi.org/10.1162/IMAG.a.35.