Impact of the acquisition protocol on the sensitivity to demyelination and axonal loss of clinically feasible DWI techniques: a simulation study

Abstract

Objective

To evaluate: (a) the specific effect that the demyelination and axonal loss have on the DW signal, and (b) the impact of the sequence parameters on the sensitivity to damage of two clinically feasible DWI techniques, i.e. DKI and NODDI.

Methods

We performed a Monte Carlo simulation of water diffusion inside a novel synthetic model of white matter in the presence of axonal loss and demyelination, with three compartments with permeable boundaries between them. We compared DKI and NODDI in their ability to detect and assess the damage, using several acquisition protocols. We used the F test statistic as an index of the sensitivity for each DWI parameter to axonal loss and demyelination, respectively.

Results

DKI parameters significantly changed with increasing axonal loss, but, in most cases, not with demyelination; all the NODDI parameters showed sensitivity to both the damage processes (at p < 0.01). However, the acquisition protocol strongly affected the sensitivity to damage of both the DKI and NODDI parameters and, especially for NODDI, the parameter absolute values also.

Discussion

This work is expected to impact future choices for investigating white matter microstructure in focusing on specific stages of the disease, and for selecting the appropriate experimental framework to obtain optimal data quality given the purpose of the experiment.

Appendices

Appendix A: Intrinsic diffusivity in the intra-myelin space

In our synthetic model of WM, we set a priori the value of water intrinsic diffusivity, into the synthetic myelin, equal to the value set for the other compartments i.e. ${D}_{\mathrm{int}}(\mathrm{sim})=2.02\times {10}^{-9}{\mathrm{m}}^{2}/\mathrm{s}$, since there are no suitable measurements of the intrinsic diffusivity of water in the real myelin. Thus, in this Appendix, we will evaluate the plausibility of our hypothesis by comparing the apparent diffusivity calculated along the radial direction in the synthetic myelin, with the apparent diffusivity measured along the radial direction in the real myelin. In detail, Andrews et al. conducted an ex-vivo study on the excised frog sciatic nerve and measured the directional apparent diffusion coefficients in the intra-myelin space [57], reporting ${\mathrm{ADC}}_{\mathrm{m}\parallel }=3.7 \times {10}^{-10} {\mathrm{m}}^{2}/\mathrm{s}$ and ${\mathrm{ADC}}_{\mathrm{m}\perp }=1.3 \times {10}^{-10} {\mathrm{m}}^{2}/\mathrm{s}$ along, respectively, the axial and the radial direction of the axon. In our simulation, by setting the intrinsic diffusivity in the intra-myelin space equal to ${D}_{\mathrm{int}}(\mathrm{sim})=2.02\times {10}^{-9}{\mathrm{m}}^{2}/\mathrm{s}$, we obtained ${\mathrm{ADC}}_{\mathrm{m}\parallel }(\mathrm{sim})=0.7 \times {10}^{-10} {\mathrm{m}}^{2}/\mathrm{s}$ and ${\mathrm{ADC}}_{\mathrm{m}\perp }(\mathrm{sim})=0.06 \times {10}^{-10} {\mathrm{m}}^{2}/\mathrm{s}$; the synthetic myelinated fibers range in radius between 0.07–2.5 µm, with an average radius equal to $0.7 \, {\mu m}$. Given that the frog sciatic nerve contains myelinated fibers ranging in radius between 0.5–10 µm [91], we have to consider, at least, the differences in terms of fiber radius, to compare the directional apparent diffusivities in the myelin obtained in this simulation with those reported by Andrews et al. [57]. In the case of Brownian motion, i.e. diffusion motion in the absence of a concentration gradient, the average distance l a water molecule diffuses during the time Δ is $l=\sqrt{2{D\Delta }}$ [7]. In the intra-myelin space of a given axon A₁, along the radial direction, we can approximate the average distance l with half the difference between the outer R₁ and the inner radius r₁ of A_1, following the same procedure used by Stanisz et al. in [59]:

$$l = \frac{{R_{1} - r_{1} }}{2} = \frac{{R_{1} \left( {1 - g} \right)}}{2},$$

(5)

where g is the g-ratio = 0.7 in the healthy condition. Thus, the ratio of the apparent diffusion coefficients ${\mathrm{ADC}}_{\mathrm{m}\perp }$ in the intra-myelin space, along the radial direction, of two axons A₁ and A₂ with outer radii R₁ and R₂, is:

$$\frac{{{\text{ADC}}_{{{\text{m}} \bot 1}} }}{{{\text{ADC}}_{{{\text{m}} \bot 2}} }} = \frac{{R_{1}^{2} }}{{R_{2}^{2} }}.$$

(6)

There is no information about the average radius of the axons populating the specific sample analysed by Andrews et al. [57]. However, the typical distribution of the axonal radii of the frog sciatic nerve can be approximated with a Gamma function with a positive skewness [16], so we may presume that the average radius ranges from 0.5 to 5 µm, i.e., respectively, the inferior limit and midpoint of the radius interval [0.5–10 µm] reported in the literature for the sciatic frog nerve fibers [91]. By using Eq. 2 with ${\mathrm{ADC}}_{\mathrm{m}\perp 1}=0.06\times {10}^{-10} {\mathrm{m}}^{2}/\mathrm{s}$, i.e. our estimated value of the radial intra-myelin ADC, ${R}_{1}=0.7\,{\mu m}$, i.e. the average radius of the synthetic fibers, ${{\mathrm{ADC}}_{\mathrm{m}\perp }}_{2}=1.3\times {10}^{-10} {\mathrm{m}}^{2}/\mathrm{s}$, i.e. the intra-myelin ADC measured ex-vivo in a frog sciatic nerve by Andrews et al. in [57], we obtain an average radius of the frog sciatic nerve fibers equal to ${R}_{2}=3.3 \, {\mu m}$, that is physically plausible (0.5 µm < ${R}_{2}$<5 µm).

Appendix B: Boundary permeabilities

The synthetic tissue of WM presented in this simulation study consists of three compartments, i.e. the intra-axonal space, the intra-myelin space, and the extra-cellular space. The boundaries separating the intra-axonal and the intra-myelin spaces, and those separating the intra-myelin and the extra-cellular spaces are permeable: we had to characterize them considering that, to our knowledge, the only useful measurements of water permeability, reported in the literature, are affected by both the neuron membrane and myelin sheet.

In detail, when a walker encounters a membrane, it is either reflected or transmitted, with a membrane transmission probability P. P depends on the membrane biochemical characteristics translating to a permeability value $k$. To ensure no net flux across the membrane [32], the relation between k and P, from inside to outside and vice versa is:

$$P_{{{\text{in}}\left[ {{\text{out}}} \right] \to {\text{out}}\left[ {{\text{in}}} \right]}} = k/\left( {1/4\cdot\nu_{{{\text{in}}\left[ {{\text{out}}} \right]}} } \right),$$

(7)

where

$$\nu_{{{\text{in}}\left[ {{\text{out}}} \right]}} = \sqrt {6D_{{{\text{in}}\left[ {{\text{out}}} \right]}} /\delta t} ,$$

(8)

$\nu$ is the average speed of walkers, $\delta t$ is the unit time-step, and D is the intrinsic diffusivity. In this simulation, $\nu$ is obtained from Eq. 2, with $\delta t=15.6 {\mu s}$ and D=$2.02\times {10}^{-9}{\mathrm{m}}^{2}/\mathrm{s}$ as described in the section Materials and Methods; P is obtained from Eq. 1, given $\nu$ and the permeability k, that is the only free parameter of the simulation strictly linked to the boundary characterization.

As described above, in our synthetic tissue there are two kinds of membranes to characterize: k_a and P_a relate to the membrane separating the axon from the myelin; k_outer and P_outer refer to the boundary separating the myelin from the extra-cellular space. Notably, in this simulation, P_outer represents the probability, for a water molecule, to cross the whole myelin sheet, thus, we set k_outer to obtain physically plausible values for P_outer.

Stanisz et al. conducted a study on the bovine optic nerve [59] and reported an overall value of ${k}_{\mathrm{a}-\mathrm{outer}(1)}=0.9\pm 0.2\times {10}^{-5}\,\,\mathrm{m}/\mathrm{s}$ for an axonal membrane encompassing the myelin sheath. We set $0.7\times {10}^{-5}\,\mathrm{m}/\mathrm{s}$ as the lower limit of k_a. Similarly, an overall value of${k}_{\mathrm{a}-\mathrm{outer}(2)}=4.7\pm 3\times {10}^{-5}\,\mathrm{m}/\mathrm{s}$, was reported by Boss et al., in a cell culture study [58], as the water membrane permeability of the mice neurons: we set $5\times {10}^{-5}\,\mathrm{m}/\mathrm{s}$ as the upper limit of k_a. Initially, we set the values of k_a and k_outer equal for all the myelinated axons. We performed several simulation trials by setting k_a values in the range [$0.7 - 5.0]\times {10}^{-5}\,\mathrm{m}/\mathrm{s}$, and ${k}_{\mathrm{outer}}<{k}_{\mathrm{a}},$ since the myelin sheath is presumably less permeable than a simple cell membrane. We evaluated the impact of P_a and P_outer on DW signal changes induced by increasing the g-ratio, to simulate demyelination, using a typical sequence for DTI analysis. We expected that, with a decrease in myelin, the bulk apparent diffusion coefficient (derived from DTI) along the axial direction AD should not significantly change, whereas the one along the radial direction RD should significantly increase, as widely reported in DWI studies on shiverer mice CNS presenting demyelination but no axonal injury or inflammation [60,61,62,63]. On the contrary, in silico, when k_outer values are fixed to be equal for all the axons and lower than$0.7\times {10}^{-5}\,\mathrm{m}/\mathrm{s}$, an increase of g-ratio causes AD to decrease and RD to slightly decrease, albeit not significantly, regardless of the choice of k_a and the thickness of the myelin. This means that, if the myelin sheath is modeled as totally impermeable (${k}_{\mathrm{outer}}=0 \,\mathrm{m}/\mathrm{s}$) or slightly permeable (${k}_{\mathrm{outer}}<0.7\times {10}^{-5} \,\mathrm{m}/\mathrm{s}$), the trends of AD and RD obtained in simulation with increasing g-ratio do not match those reported in the literature [60,61,62,63].

We obtained the expected trends for AD and RD by modelling P_outer to exponentially decrease with myelin thickness and setting ${k}_{\mathrm{a}}=5\times {10}^{-5}\,\mathrm{m}/\mathrm{s}$ equal for all the axons. This kind of model for P_outer is conceived by considering that the true myelin sheath is a greatly extended plasma membrane of a glial cell (oligodendrocyte in WM) wrapped around the axon in a spiral fashion: therefore, a water molecule moves through it without being transported through aquaporins, crosses N cell membranes with a probability P_outer exponentially decreasing with myelin thickness:

$$P_{{\text{outer }}} = P_{{\text{g}}}^{N} ,{\text{ with }}P_{{\text{g}}} = P_{{\text{a}}} ,$$

(9)

where ${P}_{\mathrm{g}}$ is the membrane crossing probability of a glial cell: in this model we assume that ${P}_{g}$ is the same as for the axonal membrane crossing probability P_a. To realize such a dependence of P_outer on the myelin thickness, we set k*_outer for one particular axon A* and obtained the relative value of P*_outer from the Eq. 1. Then, we extracted the number N* of membrane turns wrapping around A* by using the Eq. 3. Finally, we calculated Nⁱ, and then Pⁱ_outer for all other axons Aⁱ, by considering the ratio in myelin thickness between A* and Aⁱ.

After implementing the ${P}_{\mathrm{outer}}$ model_, further simulation trials were performed by changing ${k}_{\mathrm{a}}$ in the range $[0.9 - 5]\times {10}^{-5}\,\mathrm{m}/\mathrm{s}$. Lower ${k}_{\mathrm{a}}$ values determine AD decrements, and non-significant RD increments, when the g-ratio increases, hence we definitely set ${k}_{\mathrm{a}}=5\times {10}^{-5}\,\mathrm{m}/\mathrm{s}$.

Finally, ${k*}_{\mathrm{outer}}$ was set to the value $1.5\times {10}^{-8}\,\mathrm{m}/\mathrm{s}$ to set a balance between its effects on the outer membrane permeability and the pre-exchange lifetime in myelin, which increase and decrease with increasing ${k*}_{\mathrm{outer}}$ values, respectively.

Technically, the pre-exchange lifetime ${\tau }_{\mathrm{m}}$ is the time required for 63% (e⁻¹) of the intra-cellular water to exchange [38]. In this regard, many studies conducted on animal models [92,93,94] have reported ${\tau }_{\mathrm{m}}$ values in the range $[43-150] \,\mathrm{ms}$, depending on the distribution of the axonal diameters: ${\tau }_{\mathrm{m}}=43 \,\mathrm{ms}$ was measured by Harkins et al. [92] in a rat spinal cord portion with a mean axonal radius of $(0.6\mp 0.2) \, {\mu m}$ and a myelin volume fraction of 0.51. In our simulation, the mean axonal radius is $0.7 \, {\mu m}$, and the myelin volume fraction is 0.3, therefore, to compare the value of the pre-exchange lifetime in real myelin (as reported in [92]) with our simulation, we have to take at least into account the difference in myelin volume fraction since the mean axonal radii are very similar.

Considering the linear dependence between the pre-exchange lifetime in a given space and volume fraction of the same space [38], we expect the following relation:

$$\frac{{{\text{MVF}}_{r} }}{{{\text{MVF}}_{s} }} = \frac{{\tau_{r} }}{{\tau_{s} }},$$

(10)

where MVF_r and ${\tau }_{r}$, and MVF_s and ${\tau }_{s}$, are, respectively, the myelin volume fraction and pre-exchange lifetime in real myelin and in simulation. From Eq. 4, given ${\mathrm{MVF}}_{r}=0.51$, ${\mathrm{MVF}}_{s}=0.3$_, and ${\tau }_{r}=43 \,\mathrm{ms}$, we obtain ${\tau }_{s}=25 \,\mathrm{ms}$. In this simulation, given all the free parameters (diffusivity and permeabilities) set as discussed above, we obtained a pre-exchange lifetime in the synthetic myelin equal to ${\tau }_{m}=10 \,\mathrm{ms}$: this value is in the correct order of magnitude and is physically plausible, considering the different experimental conditions in which the histological ex-vivo values [92] were obtained.

Appendix C: MR diffusion signal

The signal loss for each walker is computed considering the total time $\Delta {t}_{\mathrm{a}},\Delta {t}_{\mathrm{m}}, \Delta {t}_{\mathrm{ec}}$ spent during the TE, respectively, in the intra-axonal, intra-myelin, and extra-cellular space. Dephasing depends on the net displacement during the diffusion time Δ (assuming the short-pulse approximation). Ensemble effects are computed by integrating dephasing in time and summing over all the walkers:

$$\frac{S}{{S_{0} }} = \frac{{\mathop \sum \nolimits_{{{\text{walkers}}}} \cos \theta \left( \Delta \right) \cdot e^{{ - \left( {\frac{{\Delta t_{a} }}{{T2_{o} }} + \frac{{\Delta t_{m} }}{{T2_{m} }} + \frac{{\Delta t_{ec} }}{{T2_{o} }}} \right)}} }}{{\mathop \sum \nolimits_{{{\text{walkers}}}} e^{{ - \left( {\frac{{\Delta t_{a} }}{{T2_{o} }} + \frac{{\Delta t_{m} }}{{T2_{m} }} + \frac{{\Delta t_{ec} }}{{T2_{o} }}} \right)}} }}\;{\text{with}}\;\Delta t_{a} + \Delta t_{m} + \Delta t_{{{\text{ec}}}} = {\text{TE}},$$

(11)

$cos\theta \left(\Delta \right)$ is the attenuation of the signal for one molecule in the case of a dephasing $\theta$ due to the net displacement during the diffusion time Δ; T2_o is the T2 relaxation time in the intra-axonal and extra-cellular space, while T2_m relates to the intra-myelin compartment. In vivo [95] and in vitro [96] studies on the human brain and ex vivo studies on animal models [64, 96,97,99] support T2 in the range $[10-20] \,\mathrm{ms}$ for the myelin water. We set ${T2}_{\mathrm{m}}=15 \,\mathrm{ms}$ and ${T2}_{\mathrm{o}}=78 \,\mathrm{ms},$ as Whittall et al. [64] measured, respectively, the intra-myelin space of CC and in the whole WM.

Appendix D: Implementation of the axonal loss and axonal debris processes

In neurodegenerative lesions, and in particular, in the demyelinating diseases, the major contribution to the axonal loss is given by the Wallerian Degeneration (WD) occurring when a nerve fiber is cut or crushed and the part of the axon farther from the neuron cell body degenerates [100,101,102,103]. Irrespective of the specific mechanism underlying axonal loss, the result is a decrease in axonal density together with the formation of axonal debris. As reported later, both these microstructural changes of tissue substantially affect water diffusion.

A decrease in axonal density is modeled considering an experimental observation reported in post-mortem [66, 68] and in-vivo studies [65, 69,70,71] on MS lesions: in all the analysed CNS areas, a selective death of axons with a smaller radius, roughly lower than $1 {\mu \mathrm{m}}$, occurs in both acute and chronic lesions. Thus, a fraction p_loss of axonal loss translates to a random elimination of ${p}_{loss}\bullet {N}_{h}$ axons from their distribution in healthy conditions, with ${N}_{h}$ the initial number of axons. Selective death of smaller axons is implemented by extracting the outer radii of the axons to be deleted from a sigmoidal distribution function with the inflection point at 1 μm, and further selecting the axons with a radius smaller than 1 μm, if still present in the substrate.

Regarding the axonal debris, histological studies [104, 105] described what happens immediately after the axonal degeneration. Widespread granular disintegration of axoplasmic microtubules and neurofilaments occurs in the 24–72-h interval following axonal degeneration. Amorphous granular breakdown products from nerve fibers progressively leave the intra-cellular space in 48 h and amorphous granular deposits of increasing prominence appear within the space around the axonal membranes of the lost axon. The presence of such deposits in the lesioned area have been still observed 30–34-days after nerve degeneration. In a one-year follow-up diffusion study of stroke patients, Yu et al. [106] reported at the second week a sharply decreased DTI-derived FA and AD, and increased RD; then, from the second week to the third month, MD slightly increased accompanied by a decrease in FA, an increase in RD, and no change in AD; finally, all diffusion indices remained at a relatively stable level after three months. The authors hypothesized that such changes of DTI indices might be ascribed to WD followed by the formation of axonal debris. In 2012, Qin et al. [105] conducted an in-vivo diffusion study combined with a histological validation, on a pathological animal model (felis catus) presenting WD. Their results confirmed the trends for DTI indices and also the hypothesis of Yu et al. [106]. They found quick decreases in FA, and AD, and an increase in RD from the second day to the eighth day after axonal injury, which was accompanied by progressive axonal disintegration. Afterward, from the eighth day to the sixtieth day, a slow clearance of granular deposits of axonal debris within the space around the axonal membranes, caused FA, MD, and AD to slightly increase, while RD remained unchanged [105]. In detail, after 60 days, the AD value was found still 0.83 times lower than its value in a healthy condition. Another diffusion study on WD, by Thomalla et al. [107], reported an AD value, in the lesioned tissue, to be 0.86 times lower than its value in the healthy condition, 2 weeks after the damage.

Given all these findings we modeled axonal debris having in mind these considerations:

• 1. We hypothesized to scan the tissue when the amorphous granular deposits are in the extra-cellular space, around the axonal membrane of the lost axon, as mentioned before and observed in [104]. Thus, our model represents the lesioned tissue, at least 2 weeks after the axonal degeneration, when AD becomes relatively stable;

• 2. Since the granular deposits are similar to small organelles, the diffusion motion, occurring in the extra-cellular space in the presence of the amorphous granular deposits, may be considered as hindered (Gaussian-type);

• 3. Compared to a healthy condition, the presence of the granular deposits make AD and FA decrease, RD increase, and MD slightly increase;

• 4. Given the high and homogeneous density of the axons inside the substrate, in the healthy condition, we consider the extra-cellular space as approximately equally partitioned among the axons. Thus, for a given fraction p_loss of lost axons, we assume that a reduction in diffusivity, due to the presence of axonal debris, occurs in the same fraction p_loss of the extracellular space. Indeed, the axonal debris only affects the extracellular space closely surrounding each lost axon.

We modeled the presence of axonal debris in the extra-cellular space, by decreasing the diffusivity from $D=2.02\times {10}^{-3}{\mathrm{mm}}^{2}/\mathrm{s}$ (in the extra-cellular space in the healthy condition) to a lower value ${D}_{\mathrm{debris}}$ in a fractional volume of extra-cellular space equal to ${p}_{loss}$. Accordingly, the average diffusivity $\stackrel{-}{D}({p}_{\mathrm{loss}})$ in the extra-cellular space is given by:

$$\overline{D}\left( {p_{{{\text{loss}}}} } \right) = p_{{{\text{loss}}}} \times D_{{{\text{debris}}}} + \left( {1 - p_{{{\text{loss}}}} } \right) \times D .$$

(12)

We evaluated the effect of different ${D}_{\mathrm{debris}}$ on the AD and RD trend when ${p}_{\mathrm{loss}}$ increases. We expected a significant AD decrement and RD increment with axonal loss, as reported in DWI studies on axonal injury [100,101,102,103,104,105]. If ${D}_{\mathrm{debris}}=0,$ an increase in ${p}_{\mathrm{loss}}$ determines an increase in both AD and RD. Positive values of ${D}_{\mathrm{debris}}$ result in decreasing AD, and in a progressive reduction of the RD increment, with increasing ${p}_{\mathrm{loss}}$. Further increase in the ${D}_{\mathrm{debris}}$ value leads to a point in which RD does not significantly change, and then decreases, with increasing ${p}_{\mathrm{loss}}$. The expected AD and RD trends are obtained with approximately ${D}_{\mathrm{debris}}=1.8\times {10}^{-3}{\mathrm{mm}}^{2}/\mathrm{s}$. Setting such a value implies that when axonal loss affects the whole substrate (${p}_{\mathrm{loss}}=1)$, $\stackrel{-}{D}({p}_{\mathrm{loss}})=0.89\mathrm{ D}$, a physically plausible value, considering the results reported in the literature (in [107], AD value in the lesioned tissue was found to be 0.86 times lower than its value in the healthy condition). Figure 6 shows some examples of substrates with different degrees of demyelination and axonal loss.