Resolution limit of cylinder diameter estimation by diffusion MRI: The impact of gradient waveform and orientation dispersion.

Diffusion MRI has been proposed as a non-invasive technique for axonal diameter mapping. However, accurate estimation of small diameters requires strong gradients, which is a challenge for the transition of the technique from preclinical to clinical MRI scanners, since these have weaker gradients. In this work, we develop a framework to estimate the lower bound for accurate diameter estimation, which we refer to as the resolution limit. We analyse only the contribution from the intra-axonal space and assume that axons can be represented by impermeable cylinders. To address the growing interest in using techniques for diffusion encoding that go beyond the conventional single diffusion encoding (SDE) sequence, we present a generalised analysis capable of predicting the resolution limit regardless of the gradient waveform. Using this framework, waveforms were optimised to minimise the resolution limit. The results show that, for parallel cylinders, the SDE experiment is optimal in terms of yielding the lowest possible resolution limit. In the presence of orientation dispersion, diffusion encoding sequences with square-wave oscillating gradients were optimal. The resolution limit for standard clinical MRI scanners (maximum gradient strength 60-80 mT/m) was found to be between 4 and 8 μm, depending on the noise levels and the level of orientation dispersion. For scanners with a maximum gradient strength of 300 mT/m, the limit was reduced to between 2 and 5 μm.

double diffusion encoding

diffusion MRI

full width at half‐maximum

Gaussian phase distribution

intra‐voxel incoherent motion

oscillating gradient spin echo

single diffusion encoding

signal‐to‐noise ratio

q‐trajectory encoding.

INTRODUCTION

Axons in the white matter serve as the backbone of the brain network. The information transmission through this network is determined by the conduction velocity along the axons and the axon density, which both depend on the axon diameter.1, 2, 3 Non‐invasive methods to determine the axon diameter and the axon density are thus important for mapping the network of the brain.

Most axons have diameters between 0.2 and 20 μm.4 Large axons facilitate rapid communication, e.g. for processing of sensorimotor stimuli, and are found in structures such as the corticospinal tract and the midbody of the corpus callosum. However, large axons occupy much space, yielding low axon density, and demand much energy per bit of information transmitted.5 Smaller axons, with a diameter of 0.7 μm, minimise the energy cost per bit.5 Not surprisingly, smaller axons are the most prevalent in the brain and fewer than 1% of its axons are larger than 3 μm.6 In the optic nerve, most axons have diameters below 2 μm, with a peak of 0.7 μm.5

Diffusion magnetic resonance imaging (dMRI) may enable model‐based estimation of compartment sizes and densities.7, 8 For example, Assaf et al. investigated fixed nerves and demonstrated that the average displacement of water molecules due to diffusion is limited to approximately 2 μm in the direction perpendicular to a coherent nerve fiber structure.9 By modelling the axons as cylinders and extra‐axonal water as undergoing Gaussian diffusion, the full axon diameter distribution has been recovered from dMRI data.10, 11

Techniques for axon diameter mapping have been developed in systems capable of producing magnetic field gradients of up to 1500 mT/m,10 but human MRI scanners feature gradients with much lower amplitudes. Conventional clinical systems can deliver gradients of up to 80 mT/m, and custom systems such as the Connectom system can reach as high as 300 mT/m.12 The gradient amplitude is important, because it defines the so‐called resolution limit in diffusion MRI.13, 14 The emergence of this limit is obvious in q‐space diffusion MRI, where sizes are obtained from the width of the so‐called ensemble average propagator.9, 15, 16 This propagator is obtained by means of an inverse Fourier transform of the signal‐versus‐q curve.15 However, limited gradient performance leads to limited support in terms of high q‐values, and thus the resulting propagator is convolved with a low‐pass kernel with a width defined by the inverse of the gradient amplitude.14 As the true size goes towards zero, the size estimated from the width of the estimated propagator remains at the width of the kernel. Weaker gradients result in a wider kernel, and thereby a poorer resolution in terms of a higher value of the resolution limit. The gradient amplitude can thus be compared with the wavelength of light in optical microscopy, which defines the resolution in terms of the Abbe diffraction limit.17, 18 Coincidentally, this limit prevents accurate quantification of axon diameters below approximately 0.4 μm with conventional light microscopy.19

The resolution limit is important, not only in q‐space dMRI but also for model‐based recovery of the axon diameter.8 Approaches such as CHARMED, AxCaliber, and ActiveAx estimate the axon diameter by solving an inverse problem in which axons are modelled as straight cylinders.10, 20, 21 Specificity to the axon diameter is assumed to be obtained from the signal attenuation of intra‐axonal water; however, the sensitivity of the MR signal to small cylinder diameters is low, because a small change in the diameter produces a negligible change in the measured signal. Hence, small cylinder diameters are challenging to estimate accurately (for example, see Figure 1a–d of Dyrby et al.22). In other words, cylinders with a diameter below the resolution limit are indistinguishable from virtual cylinders with a diameter of zero. Preliminary results assuming parallel cylinders indicated that the resolution limit is approximately 6 μm for gradient amplitudes of 60 mT/m and 3 μm for amplitudes of 300 mT/m.23 More realistic cases including orientation dispersion may result in even lower sensitivity to the diameter24 and cause further complications for solving the inverse problem and interpreting its solution.25, 26, 27

Most diameter mapping studies have employed the Stejskal–Tanner experiment,28 here referred to as the single diffusion encoding (SDE) experiment following the nomenclature in Shemesh et al.29 Diffusion‐encoding techniques that go beyond SDE have recently been proposed as potential solutions to reduce the resolution limit. Such techniques have generally been adapted from the fields of porous materials research, and include the double diffusion encoding (DDE) sequence,30, 31, 32, 33, 34 oscillating diffusion encoding (ODE), also known as oscillating gradient spin echo (OGSE) or modulated‐gradient NMR,24, 35, 36, 37, 38, 39 and non‐pulsed and non‐parametric gradient waveforms,40, 41 which we refer to as q‐trajectory encoding (QTE).42 In combination with improved gradient hardware for clinical MRI,43 gradient waveforms beyond SDE may enable non‐invasive recovery of the axon diameter.

In this work, we introduce an analytical framework to predict the resolution limit for any gradient waveform. Prior approaches in this direction were fully numerical, and confined to the SDE and ODE sequences.24 We analyse three cases: the first where axons are parallel, the second where there is full axonal orientation dispersion, and the third where there is partial alignment. We limit our investigation to estimation of the diameter from intra‐axonal water diffusion, assuming axons can be modelled by impermeable and straight cylinders. Provided this assumption holds, which can be debated,27 our results can be used directly in the analysis of measurements on intra‐axonal metabolites.44 For water measurements, contributions from extra‐axonal components must be incorporated in the analysis. Such components are often assumed to exhibit Gaussian diffusion,20, 21 which may not be congruent with the physics of extracellular diffusion.45 Accounting for time‐dependent diffusion outside axons is likely needed for accurate mapping of axonal characteristics,46, 47 but investigating its impact on the resolution limit was beyond the scope of the present study. Our analysis nevertheless contributes with a lower bound on the resolution limit of the intra‐axonal component.

THEORY

We express the attenuation of the signal due to diffusion perpendicular to the main axis of a cylinder as S(b|d), where b refers to the diffusion encoding strength (b‐value) and d to the cylinder diameter. We will use this notation to derive the resolution limit ( ), which we define as the diameter where the signal attenuation is indistinguishable from the case where the diameter goes towards zero, We will derive the value of for three cases: parallel cylinders, randomly ordered cylinders, and finally for any level of orientation dispersion. We assume all cylinders to be equal in size. We limit the analysis to one‐dimensional (1D) waveforms, like those applied in SDE and ODE, but the analysis is applicable also to 1D aspects of multi‐dimensional diffusion encoding, such as DDE or QTE. For all waveforms, note that g(t) must fulfil the condition in order to form an echo at the centre of k‐space. Here and throughout the analysis, we will assume g(t) to describe the effective gradient waveform after effects of RF pulses have been accounted for.

Defining the resolution limit

To define the resolution limit, we begin by defining the difference in signal (ΔS) between cylinders with a diameter approaching zero and those with a diameter corresponding to the resolution limit, We formally define the resolution limit as a hypothesis test of whether the observed ΔS is statistically higher than zero. Assuming that, due to noise, ΔS is normally distributed,48 we phrase the test in terms of a requirement on the z‐score z(ΔS)⩾z , where z is the z‐threshold for significance level α and and σ is the standard deviation of the signal due to noise, defined from the signal‐to‐noise ratio (SNR) and the signal amplitude at b=0(S 0) according to SNR=S 0/σ. Moreover, n is the number of signal averages. Deriving d from a model fitted to data acquired with multiple values of b would improve precision and be similar to increasing n, although it would be less efficient than sampling with just b = 0 and .49

Altogether, this yields a requirement for the normalised signal at the resolution limit: and If the attenuation ΔS/S 0 is less than , it will be indistinguishable from zero, i.e. the attenuation would be identical to that from a cylinder with a diameter of zero. In other words, the diameter would be below the resolution limit. When fitting models to measurements in systems with structures having sizes below the resolution limit, the diameter should not be a free model parameter. Exploiting the resolution limit thus allows for model simplifications, such as using ‘stick’ diffusion tensors with zero radial diffusivity and non‐zero axial diffusivity to model diffusion inside thin axons.26, 50

In this study, we set z to 1.64 for a one‐sided test at the 5% significance level. For a SNR of 50 (σ=S 0/SNR) and n=10, we obtain . For completeness, note that for SNR = 30 and n=1. Throughout this work, we will use the level = 1% as a reference. This is a reasonable lower limit for in vivo measurements on a clinical MRI system. Although smaller effects are detectable in principle, in practice they may be difficult to separate from effects not accounted for in the model, for example residual eddy currents,51 or effects that may be difficult to model accurately for non‐SDE waveforms, such as intra‐voxel incoherent motion (IVIM).52 In other words, effects smaller than 1% may be statistically significant but practically irrelevant.

In order to derive the resolution limit, we consider the attenuation for diffusion encoded in a direction perpendicular to a cylinder, given by where D ⊥ is the apparent radial diffusivity, which depends on d as well as on the timing of the gradient waveform used for the diffusion encoding. The approximation of the exponential is valid where the attenuation factor b D ⊥ is small, which is true per definition at the resolution limit. Since , the expression for ΔS in Equation 3 is reduced to

Parallel cylinders

We begin by deriving the resolution limit for parallel cylinders, first for the SDE sequence and then for the case of arbitrary gradient waveforms.

Single diffusion encoding sequence

Three parameters define an SDE experiment: the duration and leading‐edge separation of the gradient pulses (δ and Δ, respectively), and the gradient amplitude (g). How should these three parameters be selected in order to minimise ? The question is equivalent to maximising ΔS≈b D(δ,Δ|d), where D(·) depends on the timing variables δ and Δ, and γ is the gyromagnetic ratio. Since b, and thus ΔS, increases monotonically with g, the gradient should assume its maximal value in order to minimise . In order to select the timing parameters that minimise , we express D ⊥ using the Gaussian phase distribution (GPD) approximation,53, 54 according to where k(α,β) is defined in the Appendix, D 0 is the free diffusivity of the intra‐cylinder water, and From Equations 8, 9, and 10, we obtain According to numerical computations, this expression is maximised when δ=Δ, or expressed in unitless variables when α=β. Moreover, close to the resolution limit, where d is small, α≫1. Under these conditions, we can approximate D ⊥ as55, 56 Hence This expression can be used to estimate cylinder diameters, assuming intra‐axonal‐specific data are acquired, and with prior knowledge of D 0. Potential errors in the assumed value of D 0 are not critical, since errors of up to 50 % in D 0 give at most 10–15% errors in d. The expression in Equation 14 also gives the resolution limit for parallel cylinders and SDE, according to For , g=80 mT/m, and δ=40 ms, we obtain for the high‐SNR case where and when . Since D 0 decreases with temperature, investigations of cold fixed tissue may be beneficial to reduce the resolution limit.

Spectral domain analysis of restricted diffusion

In the previous section, we derived the resolution limit for the SDE sequence. However, gradient waveforms other than SDE may offer improved sensitivity to small cylinders. In order to investigate arbitrary gradient waveforms, we analyse the diffusion process in the spectral domain,57 where the effect of diffusion encoding on the normalised signal (S/S 0) is given by where D(ω) is the diffusion spectrum and q(ω) is the Fourier transform of q(t), defined by For completeness, we note that where f is the frequency measured in Hertz (2πf=ω). For convenience, we will use both ω and f. The relation in Equation 18 is also known as Parseval's identity. For free diffusion, D(ω)=D 0, and thus Equation 16 evaluates to .

For restricted diffusion in a cylinder, as well as for diffusion between parallel plates or in spherical geometries, the diffusion encoding spectrum D(ω) can be described by a sum of Lorenzian functions57: where C and b are coefficients that depend on the geometry and are defined for a cylinder geometry in the Appendix. The Appendix also shows the derivation of this specific form of Equation  19 from the expression in Equation (36) of Stepisnik.57

For low frequencies, D(ω) can be approximated by a second‐order polynomial47, 57, 58 where k=7/1536. This approximation yields reasonably accurate signal predictions as long as the diffusion encoding spectrum has negligible power for frequencies above a cut‐off frequency f 0. The difference between the true spectrum and the approximation is less than 20% as long as , which can also be seen in Figure 1. This inequality can be used to define f 0 according to For and d=2μm, we obtain f 0≈500 Hz. For reference, note that the maximal frequency of a sine or cosine wave that utilises the maximal gradient amplitude is given by , which equates to approximately 400 Hz for a high‐performance clinical scanner with  mT/m/ms and  mT/m. As a consequence, the approximation in Equation 20 can be used in Equation 16 to predict the signal with reasonable accuracy for many of the waveforms that are useful in practice.

Figure 1

Diffusion and encoding spectra, shown by black lines and blue areas. The dashed line shows the low‐frequency approximation. The spectra were generated for d=3μm and , and were normalised by the bulk diffusivity (D 0). Insets show gradient waveforms (g) generated for an SDE experiment with δ=10 ms and Δ=15 ms in panel A, with a sine wave with f=150 Hz in panel B, and a cosine wave (f=150 Hz) in panel C. The b‐values for these waveforms were 0.5, 0.05, and 0.017 ms/μm2

Resolution limit for parallel axons and general waveforms

A key result of the present analysis concerns the implications of the low‐frequency approximation for the resolution limit. Following Equations 7, 16, and 20, we obtain An important feature of this relation is that it separates the effects of the geometry (d 4) from the diffusion encoding (the integral part). By utilising Parseval's identity and Equation 22, we can simplify the integral further: We use this result to define an important entity: where since We refer to the entity V as the ‘spectral encoding variance’, since it is defined from the second moment of q(ω). The unit of V is 1/s2.

The product b V captures all features of the gradient waveforms required to predict the signal attenuation perpendicular to the cylinder. The b‐value captures the attenuation and V the time‐dependence, since we can approximate D ⊥ for an arbitrary gradient waveform by The signal attenuation is thus approximated by The expression above can be shown to be equivalent to Equation (119) of Grebenkov,59 which was derived for diffusion in the motional narrowing regime. Previous work on the motional narrowing regime in the context of spin relaxation has found similar equations relating the square of the gradient field, the compartment size to fourth order, and the inverse of the bulk diffusion coefficient.60, 61 In that context, motional narrowing occurs when spins diffuse rapidly through an inhomogeneous field. In our context, the motional narrowing regime occurs when the gradient varies slowly compared with the time it takes to traverse the confinement. Formally, this can be expressed as the case where the gradient waveform contains no energy above f 0. At this cut‐off frequency, the root‐mean‐square displacement due to diffusion is of the order of the compartment size ( ). This approach to the motional narrowing regime is an extension to that described by Hurlimann et al.,62 who analysed three different temporal regimes of the constant gradient experiment.

To obtain another key result, we extend the expression in Equation 24 according to where T is the duration of the gradient waveform and η is a time‐invariant factor that depends on the waveform as For SDE, T=δ+Δ, and by assuming negligible ramp times we obtain and Note the similarity of the expression for V and the denominator in Equation 13. Combining Equations 13 and 32 gives Equation 27, which shows that the same solution is obtained for both time‐domain and frequency‐domain approaches in the SDE case.

To evaluate the resolution limit in the case of parallel cylinders and arbitrary gradient waveforms, we use Equations 8 and 27 to define Together with the definition of the resolution limit in Equation 5, we obtain a general expression for the resolution limit in the case of parallell cylinders (par): In order to examine how to optimise the gradient waveform to minimise , we first note that k, D 0, and γ are independent of the gradient waveform. Moreover, note that , since the SNR is reduced as more T 2 relaxation takes place for long encoding gradients (the echo time, TE, is given by T+T 0, where T 0 is the time required for imaging gradients and RF pulses). Utilising Equation 29, we now obtain a simplified expression, This equation yields two important results. First, the resolution limit is minimised if T=T 2, since this value of T minimises . The optimal duration of the diffusion encoding is thus equal to the value of T 2 (approximately 80 ms for white matter at B 0=3T). Second, η should be maximised for optimal resolution. Note that η obtains its maximal value of unity if and only if the gradients are at full amplitude during the whole period of T(Equation 30), while still fulfilling Equation 2. This can be obtained with SDE if δ=Δ (see Equation 31). We have thus reproduced the result behind Equation 15. A square gradient waveform would be equally as good as SDE, but constraints such as limited slew rates reduce the value of η proportional to the time required for slewing, which is greater for gradients with multiple pulses. Hence, this result shows that SDE yields a value of lower than what is possible with DDE and ODE, due to the effects of limited slew rates. DDE with short mixing times increases η, and it is thus not surprising that short rather than long mixing time DDE experiments are preferred for size estimations.63 Waveforms with multiple oscillations may result in encoding spectra with power at frequencies above f 0(Equation 21). In that case, the low‐frequency assumption would not hold, but the resulting signal attenuation would be lower than predicted, yielding a poorer resolution (higher ).

Orientation dispersion

Most, if not all, regions of the brain feature axonal bundles with different orientations.64, 65, 66 We therefore proceed to analyse the resolution limit in the case of full orientation dispersion (cylinders oriented along all directions with equal probability). In this scenario, the signal is reduced due to increased attenuation from diffusion weighting along the fibers, which gives67, 68, 69 where and Assuming D ⊥≈0 close to the resolution limit, h(A) becomes independent of d.

The resolution limit for dispersed cylinders is now given by rearranging Equation 5 so that and thus Since h(A)⩽1, this shows that the resolution is worse in the dispersed case ( ).

Optimisation for unlimited slew rate

To gain some intuition on how to choose an optimal waveform for obtaining the best resolution in the case with complete orientation dispersion, consider a waveform composed of m identical pulsed gradient pairs, i.e. a square wave. For now, we assume that the slew rate is infinite. The total b‐value is then given by b=m b 0, where and δ=T/2m, assuming that Δ=δ for each pulsed gradient pair. For a train of diffusion encoding pulses, we obtain70 For completeness, we note that, for an SDE experiment with maximal duration of the gradients, we have m=1, δ=Δ=T/2, and thus as expected.

In order to maximise ΔS, and thus minimize , it is sufficient to consider the terms affected by the gradient waveform, Note that, for a perfect square wave with η=1, we get Hence, the term b V is here independent of the number of oscillations. Maximisation of ΔS is thus done by maximising h(A) for a fixed value of T. Curiously, h(A) is at its maximum when b=0, which we obtain when . To understand this result, we first note that the resolution limit is in part determined by the SNR, which in turn is reduced by high b due to partial alignment of the dispersed cylinders and the encoding direction. Increasing m has no effect on the radial attenuation component, but decreases b, which is beneficial for the SNR and thus reduces the resolution limit. In practice, however, when , due to limited slew rates. An optimal waveform will thus need to balance the two objectives of maximising η while minimising b.

Intermediate orientation dispersion

In the intermediate case, where the level of orientation dispersion is somewhere between full coherence and full dispersion, we found it challenging to derive an analytical expression for ΔS that is also illuminating. The most straightforward way we found was to extend h(A), assuming D ⊥≈0, so that where A is defined in Equation 38 and C d is the orientation dispersion factor, defined in the interval from zero to unity, which represents full coherence and full dispersion, respectively. This factor is defined by where κ is the orientation dispersion factor in the Watson distribution.26 In analogy with the derivation for the complete orientation dispersion case in Equation 40, we thus describe the resolution limit in the intermediate (int) case by Since h(A,0)=1 for orientation coherence, the expression agrees with Equation 34. With full dispersion and high b‐values, h(A,1)≈h(A), which agrees with Equation 40. For low b‐values, for example due to oscillating gradients, and thus , so that . This result is in agreement with the previous notion that, with dispersion, oscillating gradients are preferred, since the attenuation resulting from axial diffusion is thus minimised.

METHODS

We first verified the calculations of D ⊥(d) for SDE, DDE, and ODE by using analytical expressions and Monte Carlo random‐walk simulations. We also verified the calculations of using numerical calculations for both parallel and fully dispersed cylinders. Second, we optimised gradient waveforms in order to minimise . Finally, we investigated the capability to recover cylinder diameters using numerical simulations, and tested whether the theoretically predicted values of agreed with the simulation results.

Model verification

We suggest that the frequency‐domain approximation in Equation 20 can be used to calculate D ⊥ from d for any gradient waveform, as long as the condition in Equation 21 is fulfilled. We tested this assertion by comparing D ⊥ calculated from Equation 20 with values predicted for SDE by Equation 10 using the GPD approximation. These verifications were performed for , and d in the range 0–10μm. Note that we have previously compared the GPD approximation with results from Monte Carlo simulations in cylinders.71

For a wide range of waveforms including both SDE and DDE, we compared the value of D ⊥ predicted from the frequency‐domain approximation in Equation 20 with the value obtained from Monte Carlo simulations, using an approach described previously.71 In short, random walkers (n=5×104) were placed on a grid within a circle with reflective boundaries and a diameter of d. The step size was set to Δx=0.08μm, and accordingly Δt=Δx 2/2D 0=1.6μs. The phase ϕ accumulated by particle i from a gradient waveform, assumed to have been applied along x, is theoretically given by where g (t) is the gradient and r (t) is the position of the particle. In the simulations, ϕ was obtained according to where m=T/Δt is the number of time steps. The MR signal was given by , where ⟨·⟩ refers to ensemble averaging and |·| yields the magnitude of a complex signal. For moderate attenuations, i.e. for low values of b, the normalised signal can be approximated according to54 and thus The values of D ⊥ obtained by simulations were compared with the theoretical predictions for different gradient waveforms.

In addition, we computed values of using numerical calculations described by Drobnjak et al.,24 and compared the results with our analytical results. Numerical calculations employed the matrix formalism implemented in the MISST software package,72 and were used to predict the diffusion‐weighted signal for a range of gradients {60, 80, 150, and 300} mT/m, SNR ∈{10, 20, and 50}, and diameters sampled finely in the range d∈[0,10]μm. We then used the same approach as in Drobnjak et al.24 to calculate and numerically find the smallest for which Equation 5 is satisfied. We performed calculations for each of the combinations of and SNR, and for both the SDE case (m=1) and a square wave with m=2, 3, and 4 identical pulsed gradient pairs. Infinite slew rates were assumed in all cases.

Optimising gradient waveforms

Gradient waveforms were optimised in order to maximise ΔS in the case of parallel cylinders and in the presence of full orientation dispersion. For this purpose, we expressed the gradient waveform in terms of a cosine series, We chose the cosine basis, since this yielded lower b‐values than the corresponding sine basis waveforms, and thus a lower value of the resolution limit in the dispersed case. The coefficients c were optimised in order to maximise ΔS, after which the resolution limit was calculated. However, before ΔS was calculated, we convolved the gradient waveform with a Gaussian kernel with a standard deviation computed to ensure slew rates below 200 mT/m/ms (see Appendix). This procedure ensured that the resulting waveform can be played out on a clinical scanner. In these optimisations, we assumed T=80 ms and hardware parameters corresponding to high‐end clinical MRI scanners and the Connectome scanner, i.e. and , respectively. We noted that the optimisations tended to yield square‐wave functions. We therefore also explicitly generated square waveforms with different frequencies.

To determine the resolution limit for an optimised waveform, the optimisation procedure was repeated for different values of d until the value of ΔS/S 0 reached 1%.

Recovery of the diameter

In order to assess the ability to recover the diameter from protocols optimised to minimise the resolution limit, we simulated the MR signal using the Monte Carlo approach described above for an SDE sequence with δ = 40 ms and Δ = 40 ms. After that, we added Gaussian noise to the MR signal (corresponding to SNR = 200, i.e. ) and estimated the diameter from the relationship in Equation 7, with D ⊥(d) given by the frequency‐domain approximation in Equation 20. Note that, for the SDE timing parameters used, the approximation is valid for d<10μm, since the spectrum of the waveform contains negligible energy above 20 Hz. In the estimation, we assumed prior knowledge of the correct value of D 0of 2 μm2/ms. The procedure of generating noise and estimating the diameter was repeated 3000 times for diameters between 0 and 10μm.

RESULTS

Figure 2 shows values of D ⊥(d) for diameters between 0 and 10 μm. For the SDE case (panel A), the low‐frequency approximation in Equation 20 agreed well with both the analytical GPD‐based approximation in Equation 10 and the Monte Carlo simulations. For the gradient waveform composed of a f=80 Hz sine wave (panel B), however, the approximation and the Monte Carlo simulations agreed only for d⩽5μm. The discrepancy for higher diameters was expected, according to the limit specified in Equation 21, which shows that the approximation should be valid for d⩽5 μm, when f=80 Hz and .

Figure 2

Low‐frequency approximation versus Monte Carlo results. The plots show the apparent radial diffusivity (D) predicted from the low‐frequency approximation (blue solid line), results from the Monte Carlo simulations (circles), and corresponding gradient waveforms (g(t)) and b‐values (inset). For the SDE case, the diffusivity predicted by the GPD approximation is also shown. The dashed line and the grey area show the mean and the 95% interval of the estimated diffusivities, assuming SNR = 1000. For small diameters (d), the low‐frequency approximation is in good agreement with the Monte Carlo simulations, whereas they start to diverge for larger diameters. In panel C, where the gradient waveform was a sine wave with f=80 Hz, the approximation agreed with the Monte Carlo simulations for d<5μm, as expected from Equation  21

Figure 3 shows the effect on the signal difference (ΔS) resulting from varying the frequency of a square‐wave oscillating gradient waveform, for cylinders where d=3.6μm. As the frequency increases, the b‐value decreases whereas D ⊥ increases. The net effect, however, is that ΔS is reduced for higher frequencies, illustrating that waveforms other than SDE lead to worse performance in terms of obtaining a minimal resolution limit for parallel cylinders. In the presence of orientation dispersion, however, the signal difference was maximised at a frequency of approximately 100 Hz. Note that here we assumed a slew rate of 200 mT/m/ms.

Figure 3

Impact of frequency on attenuation. The left panel shows ΔS, which is the signal difference between systems where d=3.6μm and d=0, for square‐wave oscillating gradients with varying frequency. Results are shown for both parallel and dispersed fibers. The two panels to the right shows the gradient waveforms that maximised ΔS for the parallel case (top) and the dispersed case (bottom)

Figure 4 shows a comparison of the resolution limit computed by our theoretical approach (Equations 34 and 40), and from the numerical approach described by Drobnjak et al.24 Investigations were performed under varying and SNR, assuming infinite slew rates. Panel A shows that numerical and analytical results were in excellent agreement for all gradient strengths and all SNR levels, for parallel cylinders. Panel B shows the results for fully dispersed cylinders. Here, the resolution limits depend on the number of oscillations m. Numerical and analytical results were well aligned except at low SNR levels, where our analytical approach underestimated the resolution limit. The discrepancies decreased at higher SNR, and for SNR=50 the match was excellent.

Figure 4

Resolution limit ( ) compared between analytical and numerical approaches. Panel A shows results for parallel cylinders under varying SNR levels, maximal gradient amplitudes ( ), and for a varying number of square‐wave oscillations (m) assuming infinite slew rates. The two approaches show excellent agreement. Panel B shows the corresponding result for the case with complete orientation dispersion. Resolution limits were higher than for the parallel case, as expected. The analytical approach underestimates the resolution limit for low SNR values, but the two approaches show good agreement for higher SNR levels. In all cases, the x‐axis was scaled according to

Optimisation

Figure 5 shows gradient waveforms optimised to minimise . In the case of parallel cylinders, the optimisation resulted in SDE waveforms, regardless of . This is not surprising, since SDE maximises η(Equation 30). In the case of fully dispersed cylinders, the optimal waveforms comprised a train of square pulses. For the case with stronger gradients, the waveforms became triangular, due to limitations in the available slew rate.

Figure 5

Waveforms optimised for best possible resolution. In the case of parallel cylinders, the optimisation resulted in SDE waveforms. For the dispersed case, waveforms included multiple gradient pulses. With a 80 mT/m system, a square‐wave oscillating waveform emerged. For the 300 mT/m system, the limited slew rate prevented the emergence of a square wave, and thus a triangular wave appeared instead

Table 1 shows the resolution limit in the cases of parallel and dispersed cylinders. The resolution limit was 20% higher in the presence of orientation dispersion for the 80 mT/m gradient system. For the 300 mT/m system, the corresponding number was 40%. For completeness, we also investigated the resolution limit for a preclinical system with 1000 mT/m gradients (slew rate 5000 T/m/s and T=30ms), and found it to be of the order of 1 μm.

Table 1

Resolution limits ( ) for waveforms optimised so that ΔS=1%, with a waveform amplitude of and a duration of T=80 ms at field strength B 0=3T and T=50 ms at B 0=7T. The unit of the encoding strength b is ms/μm2

g_max	B 0	dmin(par)	b	dmin(disp)	b
80 mT/m	3T	3.3 μm	20	3.4 μm	0.1
300 mT/m	3T	1.7 μm	260	2.6 μm	1.4
60 mT/m	7T	3.5 μm	2.4	3.7 μm	0.1
80 mT/m	7T	3.0 μm	4.5	3.2 μm	0.1

In the intermediate case between coherence and full dispersion, the resolution limit depends on the specific level of orientation dispersion. For the level of dispersion recorded for the corpus callosum (FWHM = 34° , see Leergaard et al.73), we find for the 300 mT/m system, which is approximately 15–20% higher than the case for parallel fibres. In other words, even for an orientation dispersion as small as the one in the corpus callosum, otherwise known for its high orientation coherence, the resolution is degraded.

At 7T, the SNR is higher while the T 2 relaxation time is shorter (approximately 50 ms), thus requiring the use of shorter gradient pulses. Whether there is a benefit at 7T compared with 3T depends on the gradient system. If 80 mT/m is available at both 3T and 7T, the higher field strength has an advantage (Table 1).

Recovery

The ability to recover cylinder diameters from a simple experiment with two b‐values (0 and ) was assessed by numerical simulations. As shown in Figure 6 for a SDE‐like experiment in a system with parallell cylinders, diameters were correctly recovered above the resolution limit at . However, below this limit, cylinders of different sizes were indistinguishable.

Figure 6

Recovery of cylinder diameters (d). Red lines show the 80% confidence interval. The resolution limit is shown as the blue vertical line, at m. The data represent a case with parallel cylinders and SDE encoding with δ=Δ=40ms

DISCUSSION

In this work, we assessed the minimal cylinder diameter that can be recovered from diffusion MRI of water within the cylinders. We label this diameter ‘the resolution limit’, in analogy with optical microscopy, where the diffraction limit defines the best optical resolution of the microscope. The theory presented herein allows the resolution limit to be assessed for arbitrary gradient waveforms, and also shows how to optimise the gradient waveform in order to minimise the resolution limit. In the very simple case of parallel cylinders, we demonstrated that the SDE experiment is preferred over any other waveform, which is in agreement with the findings of Drobnjak et al.,24 who used numerical simulations to compare SDE with ODE. Other studies have reported DDE to be beneficial for size estimation,34, 63 which stands in apparent contrast to our findings. However, those studies compare an optimised DDE with a suboptimal SDE, in the sense that the values of η(Equation 30) were not maximised for SDE and DDE separately (with T as the longest time common to the two sequences). Although DDE sequences are not intrinsically beneficial for size estimation, DDE has other unique applications, for example to study diffusion diffraction,74 water exchange,75, 76, 77, 78, 79 and intra‐voxel incoherent flow of blood.52 DDE is also useful for estimation of microscopic anisotropy,80, 81 although this is also possible with QTE, which may be experimentally advantageous compared with DDE.42, 82, 83, 84

The theory developed herein can be used to define three spatiotemporal regimes of the diffusion encoding. In the first regime, the diffusion is completely restricted and . In the second regime, the diffusion is partly restricted and , where f 0 is the maximal relevant frequency component of the encoding spectrum. In the third regime, the diffusion is weakly restricted and . The motional narrowing regime encompasses both the first and second regimes, whereas the third regime is a mixture of the motional narrowing and the free diffusion regimes described by Hurlimann et al.62 These regimes can be applied to simplify modelling of restricted diffusion. In the first regime, we can assume D ⊥=0, and in the second, D ⊥ can be reliably calculated using only V based on the low‐frequency approximation. In the third regime, knowledge of the full encoding spectrum is required to predict D ⊥, and thereby specific details such as the number of oscillations of the waveform become relevant.

According to our analysis, the analytically computed resolution limit is in perfect agreement with numerical results in the completely and partly restricted regimes (Figure 4). However, in the weakly restricted regime, the low‐frequency assumption does not hold. As a result, the analytically calculated resolution limit is underestimated, especially at a large number of oscillations. Nevertheless, even for these scenarios, the curves from the analytical and numerical approaches synchronously follow the same trend, determined by b V (Figure 4). Note that, in clinical scanners, gradient waveforms have little encoding power above 100 Hz, which yields d=4–5μm as the boundary between the partially and weakly restricted regimes.

The present analysis concerns diffusion in cylinders, and can be used to analyse axon diameter estimation by methods that model axons as straight cylinders.8, 20, 21 Our analysis suggests that gradient amplitudes in currently available clinical systems are insufficient to quantify the axon diameter accurately. This result is in agreement with Sepehrband et al.,85 who showed that the axon diameter estimated by such models depends on up to 1350 mT/m, which indicates that diameters probed by weaker gradients reflect the available amplitude rather than the underlying diameters. Moreover, in most of the white matter, axons do not run in straight and parallel courses, but rather in tortuous and non‐parallel configurations.27, 64 Even in a structure such as the corpus callosum, where axons are unusually coherent, there is substantial orientation dispersion.65, 66 In the presence of orientation dispersion, higher b‐values reduce the effective SNR, since the intra‐axonal water signal is attenuated proportional to the relative alignment of the axon to the direction of the diffusion encoding. Reducing the b‐value thus improves SNR, but if this is done by just reducing g, the resolution gets worse. A key result of the present study offers an alternative. Close to the resolution limit, the most important determinant of the signal attenuation is the integral of the squared gradient (b V ). This entity can be kept constant, while b and V can be varied independently by the use of square‐wave oscillating gradients. In order to improve the SNR and thus improve the resolution, waveforms should, in the presence of dispersion, use more oscillations and hence lower b‐values than for the parallel case. This finding is in agreement with the results of Drobnjak et al.24 We also see that, with orientation dispersion, the relative benefit of using strong gradients to achieve high b‐values diminishes (Table 1).

We wish to highlight three limitations concerning the extrapolation of the present theoretical results to practical estimation of axon diameters in white matter. First, we investigate a simplified case where there is only intra‐axonal water. However, we believe the resolution limits reported herein to be applicable also to multi‐component systems, at least as lower limits, since adding complexity such as partial volume effects to the model will only make it more difficult to estimate the diameter accurately, not easier. The second limitation concerns the relationship between the axon diameter and the structure of the extra‐cellular space,47 which we neglected in the present analysis by considering the intra‐axonal component only. Including the extra‐cellular space in the model may, however, contribute with information on structural dimensions on its own.47, 86 Due to its weaker frequency dependence (|ω| versus ω 2), the resolution limit for estimating structural dimensions of the extracellular versus intra‐axonal space may differ, and, for low frequencies, effects of time‐dependent diffusion in the extra‐cellular space may dominate over those in the intra‐axonal space. Third and finally, the present analysis neglects exchange between water environments. However, exchange can probably be safely neglected, since our previous results have shown exchange times in white matter that are of the order of seconds or longer,78 which is much longer than the time‐scales during which effects of restricted diffusion can be observed.

Finally, note that the diffraction limit in optical microscopy, introduced in 1873, was recently broken.87 It took 135 years. Perhaps breaking the resolution limit in diffusion MRI can be done a little faster?