AxCaliber with time-dependent zeppelin¶

In this example, we use the AxCaliber method (Assaf et al. 2008) presented in a previous example while including the diffusion time-dependency of the extra-axonal space (Novikov et al. 2014, Burcaw et al. 2015).

The signal representation is

\begin{equation} E_{\textrm{AxCaliber}}^{\textrm{Temporal}}= \underbrace{(1-f_r)\overbrace{E_h(b,\textbf{n},\lambda_\parallel,D_{\infty},A)}^{\textrm{GPA-Zeppelin}}}_{\textrm{Extra-Axonal}} + \underbrace{f_r \overbrace{\Gamma(\alpha,\beta)}^{\textrm{Gamma Distribution}}*_{\mathbb{R}}\overbrace{E_r(\cdot|D_\perp)}^{\textrm{Cylinder}}}_{\textrm{Intra-Axonal}} \end{equation}

where the cylinder was initially represented using the Callaghan model, but was later replaced with the more general Van Gelderen model (Huang et al. 2015, De Santis et al. 2016). The initial Ball, which accounted for the extra-axonal signal, is replaced by the time-dependent Zeppelin (De Santis et al, 2016). The formulation of the time-dependent Zeppelin is as follows:

\begin{equation} E_h(b,\textbf{n},\lambda_\parallel,D_{\infty},A)=\exp(-b\textbf{n}^T(\textbf{R}\textbf{D}^h_{\textrm{diag}}\textbf{R}^T)\textbf{n})\quad\textrm{with}\quad \textbf{D}^h_{\textrm{diag}}= \begin{pmatrix} \lambda_\parallel & 0 & 0 \\ 0 & D_{\infty}+A\frac{\ln(\Delta/\delta)+3/2}{\Delta-\delta/3} & 0 \\ 0 & 0 & D_{\infty}+A\frac{\ln(\Delta/\delta)+3/2}{\Delta-\delta/3} \end{pmatrix} \end{equation}

here the perpendicular diffusivity is now dependent on $\Delta$ and $\delta$, with $D_{\infty}$ the bulk diffusion constant and $A$ a characteristic coefficient for extra-axonal hindrence.

Our goal here is to reproduce the Monte Carlo simulations presented in (De Santis et al, 2016), where they investigated the impact of including extra-axonal hindrence through the temporal zeppelin on axonal density and diameter estimation.

AxCaliber model with time-dependent zeppelin¶

Using Dmipy, we can construct the time-dependent version of AxCaliber by calling G3TemporalZeppelin instead of the previously used G1Ball.
As before, the Gamma-distributed Gaussian-Phase cylinders are instantiated by first calling the cylinder model, and then using it as input for the GammaDistributed model.

Then put the models together in a MultiCompartmentModel and display the parameters. Note that this model has a lot of parameters.

Visualize the model:

The data we will fit is measured only along 2 gradients perpendicular to the fiber axis (along x and y axis). This means there is no data that provides information to fit parallel diffusivity "lambda_par" or the axon orientation "mu".

So first, we make sure that the intra- and extra-axonal orientations are aligned and that their parallel diffusivities are equal by using the set_equal_parameter() :

Then to not fit these superfluous parameters, we fix the orientation mu along the z-axis and give lambda_par an arbitrary value (it will not affect the results).

Notice that we are now only fitting the volume fractions, the Gamma distribution parameters $\alpha$ and $\beta$, the bulk diffusivity $\lambda_\infty$ and the coefficient for extra-axonal hindrance $A$.

Load Monte Carlo simulated data¶

Using Camino, we previously generated a substrate of Gamma-distributed cylinders and simulated the perpendicular acquisition scheme that was used by (De Santis et al. 2016).

Unlike De Santis et al who simulated STEAM acquistion scheme, we simulate a PGSE scheme considering their parameters: $\delta = 17$ ms, $\Delta = 48,60,80,100,120,140,160,180,195$ ms, 4 b-values $(500, 1000, 2000$ and $4000$ $s/mm^2$) and 2 unweighted images for each $\Delta$.

To estimate the robustness of the model's parameter estimates, we simulate 20 Rician noise instances on the data at SNR=30.

Fit time-dependent AxCaliber model to simulated data¶

Visualize fitted parameters¶

We use the gaussian_kde function to show the distribution of the fitted parameters. For each figure, the red vertical line is the ground truth while the green vertical one represents the mean value of the estimated parameter.

Notice that the estimated mean axon diameter either defaults to "Stick-like" behavior or finds axon diameter much too large. The mean is therefore mostly outlier driven.

Except for the axon diameter estimation, our results are comparable to De Santis et al results (of the 3rd simulation presented in table 1). Indeed, the mean of the estmated axonal density $f_r$ is close to the ground truth, however, with this scheme, the axonal diameter is badly estimated.

As for the characteristic coefficient A, the estimated value is very close to De Santis et al, while its variance is a bit larger. And the mean of the estimated bulk diffision constant $D_{\infty}$ is 3 times bigger. Further study is needed to explain the differences between De Santis et al's findings and ours.

References¶

• Assaf, Yaniv, et al. "AxCaliber: a method for measuring axon diameter distribution from diffusion MRI." Magnetic resonance in medicine 59.6 (2008): 1347-1354.
• Novikov, Dmitry S., et al. "Revealing mesoscopic structural universality with diffusion." Proceedings of the National Academy of Sciences 111.14 (2014): 5088-5093.
• Huang, Susie Y., et al. "The impact of gradient strength on in vivo diffusion MRI estimates of axon diameter." NeuroImage 106 (2015): 464-472.
• Burcaw, Lauren M., Els Fieremans, and Dmitry S. Novikov. "Mesoscopic structure of neuronal tracts from time-dependent diffusion." NeuroImage 114 (2015): 18-37.
• De Santis, Silvia, Derek K. Jones, and Alard Roebroeck. "Including diffusion time dependence in the extra-axonal space improves in vivo estimates of axonal diameter and density in human white matter." NeuroImage 130 (2016): 91-103.