AxCaliber¶
To estimate the axon diameter distribution itself, (Assaf et al. 2008) proposed a method called AxCaliber. In AxCaliber, a distribution of cylinders with different diameters are fitted to the diffusion signal. In particular, a Gamma distribution is used. A requirement of this technique is that the signal is measured perpendicular to the axon bundle direction for different gradient strengths and diffusion times. The signal representation in AxCaliber is given as
\begin{equation} E_{\textrm{AxCaliber}}= \underbrace{(1-f_r)\overbrace{E_h(\lambda_{\textrm{h}})}^{\textrm{Ball}}}_{\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 (Assaf et al. 2008, Huang et al. 2015, De Santis et al. 2016).
Advantages:
- Estimates the axon diameter distribution.
Limitations:
- Tissue orientation must be known for the signal to be measured perpendicular,
- Does not take axon dispersion into account,
- Does not include the effects of extra axonal restriction (Burcaw et al. 2015).
AxCaliber model with Gaussian Phase Approximation¶
The process of setting up the AxCaliber model with Gamma-distributed cylinder diameters is the same as setting up an axon dispersed model. We import the needed Ball and Cylinder models, and then distribute the cylinder using the DD1GammaDistributed function from distribute_models.
In [1]:
from dmipy.distributions import distribute_models
from dmipy.signal_models import gaussian_models, cylinder_models
ball = gaussian_models.G1Ball()
cylinder = cylinder_models.C4CylinderGaussianPhaseApproximation()
gamma_cylinder = distribute_models.DD1GammaDistributed(models=[cylinder])
Then put the models together in a MultiCompartmentModel and display the parameters
In [2]:
from dmipy.core import modeling_framework
axcaliber_gamma = modeling_framework.MultiCompartmentModel(models=[ball, gamma_cylinder])
axcaliber_gamma.parameter_cardinality
Out[2]:
OrderedDict([('DD1GammaDistributed_1_DD1Gamma_1_alpha', 1),
('G1Ball_1_lambda_iso', 1),
('DD1GammaDistributed_1_C4CylinderGaussianPhaseApproximation_1_mu',
2),
('DD1GammaDistributed_1_DD1Gamma_1_beta', 1),
('DD1GammaDistributed_1_C4CylinderGaussianPhaseApproximation_1_lambda_par',
1),
('partial_volume_0', 1),
('partial_volume_1', 1)])
Visualize the model:
In [3]:
from IPython.display import Image
axcaliber_gamma.visualize_model_setup(view=False, cleanup=False, with_parameters=True)
Image('Model Setup.png')
Out[3]:
The data we will fit is measured only along 2 gradient orientations that are perpendicular to the axon 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". To not fit these superfluous parameters we fix mu along the z-axis and give lambda_par an arbitrary value (it will not affect the results).
In [3]:
axcaliber_gamma.set_fixed_parameter(
'DD1GammaDistributed_1_C4CylinderGaussianPhaseApproximation_1_lambda_par', 1.7e-9)
axcaliber_gamma.set_fixed_parameter(
'DD1GammaDistributed_1_C4CylinderGaussianPhaseApproximation_1_mu', [0, 0])
Load AxCaliber-compatible data¶
For the example data, we kindly make use of a cat spinal cord dataset hosted at the White Matter Microscopy Database (https://osf.io/yp4qg/). The data includes 2D Axcaliber data, 3D multi-shell data and the corresponding histology values for axon diameter, intra-axonal volume fraction (and others).
In [4]:
from dmipy.data import saved_data
scheme_spinal_cord, data_spinal_cord = saved_data.duval_cat_spinal_cord_2d()
This data was used by Duval et al. 'Validation of quantitative MRI metrics using full slice histology with automatic axon segmentation', ISMRM 2016. Reference at Cohen-Adad et al. White Matter Microscopy Database. http://doi.org/10.17605/OSF.IO/YP4QG
In [5]:
import matplotlib.pyplot as plt
%matplotlib inline
plt.imshow(data_spinal_cord.signal[:,:,0,0])
plt.title('b0 of spinal cord data')
plt.axis('off');
Printing the acquisition scheme information we can see the Delta-shells of the acquisition.
In [6]:
scheme_spinal_cord.print_acquisition_info
Acquisition scheme summary total number of measurements: 1791 number of b0 measurements: 72 number of DWI shells: 418 shell_index |# of DWIs |bvalue [s/mm^2] |gradient strength [mT/m] |delta [ms] |Delta[ms] |TE[ms] 0 |96 |191 |0 |3.0 |7.0 |36.152 1 |4 |627 |402 |3.0 |7.0 |36.152 2 |4 |681 |419 |3.0 |7.0 |36.152 3 |4 |738 |437 |3.0 |7.0 |36.152 4 |4 |797 |454 |3.0 |7.0 |36.152 5 |4 |858 |471 |3.0 |7.0 |36.152 6 |4 |922 |488 |3.0 |7.0 |36.152 7 |4 |988 |505 |3.0 |7.0 |36.152 8 |4 |1056 |522 |3.0 |7.0 |36.152 9 |4 |1126 |539 |3.0 |7.0 |36.152 10 |4 |1199 |557 |3.0 |7.0 |36.152 11 |4 |1274 |574 |3.0 |7.0 |36.152 12 |4 |1351 |591 |3.0 |7.0 |36.152 13 |4 |1431 |608 |3.0 |7.0 |36.152 14 |4 |1512 |625 |3.0 |7.0 |36.152 15 |4 |1596 |642 |3.0 |7.0 |36.152 16 |4 |1683 |659 |3.0 |7.0 |36.152 17 |4 |1771 |677 |3.0 |7.0 |36.152 18 |4 |1862 |694 |3.0 |7.0 |36.152 19 |4 |1955 |711 |3.0 |7.0 |36.152 20 |4 |2051 |728 |3.0 |7.0 |36.152 21 |4 |2148 |745 |3.0 |7.0 |36.152 22 |4 |2248 |762 |3.0 |7.0 |36.152 23 |3 |2350 |779 |3.0 |7.0 |36.152 24 |3 |2455 |797 |3.0 |7.0 |36.152 25 |3 |2562 |814 |3.0 |7.0 |36.152 26 |3 |2671 |831 |3.0 |7.0 |36.152 27 |3 |2782 |848 |3.0 |7.0 |36.152 28 |20 |31 |0 |3.0 |40.0 |57.288 29 |4 |149 |77 |3.0 |40.0 |57.288 30 |4 |223 |94 |3.0 |40.0 |57.288 31 |4 |311 |111 |3.0 |40.0 |57.288 32 |4 |415 |128 |3.0 |40.0 |57.288 33 |4 |533 |145 |3.0 |40.0 |57.288 34 |4 |666 |162 |3.0 |40.0 |57.288 35 |4 |813 |179 |3.0 |40.0 |57.288 36 |4 |976 |197 |3.0 |40.0 |57.288 37 |4 |1153 |214 |3.0 |40.0 |57.288 38 |4 |1345 |231 |3.0 |40.0 |57.288 39 |4 |1551 |248 |3.0 |40.0 |57.288 40 |4 |1773 |265 |3.0 |40.0 |57.288 41 |4 |2009 |282 |3.0 |40.0 |57.288 42 |4 |2260 |299 |3.0 |40.0 |57.288 43 |4 |2526 |317 |3.0 |40.0 |57.288 44 |4 |2806 |334 |3.0 |40.0 |57.288 45 |4 |3101 |351 |3.0 |40.0 |57.288 46 |4 |3411 |368 |3.0 |40.0 |57.288 47 |4 |3736 |385 |3.0 |40.0 |57.288 48 |4 |4076 |402 |3.0 |40.0 |57.288 49 |4 |4430 |419 |3.0 |40.0 |57.288 50 |4 |4799 |437 |3.0 |40.0 |57.288 51 |4 |5183 |454 |3.0 |40.0 |57.288 52 |4 |5582 |471 |3.0 |40.0 |57.288 53 |4 |5995 |488 |3.0 |40.0 |57.288 54 |4 |6423 |505 |3.0 |40.0 |57.288 55 |4 |6866 |522 |3.0 |40.0 |57.288 56 |4 |7323 |539 |3.0 |40.0 |57.288 57 |4 |7796 |557 |3.0 |40.0 |57.288 58 |4 |8283 |574 |3.0 |40.0 |57.288 59 |4 |8785 |591 |3.0 |40.0 |57.288 60 |4 |9301 |608 |3.0 |40.0 |57.288 61 |4 |9833 |625 |3.0 |40.0 |57.288 62 |4 |10379 |642 |3.0 |40.0 |57.288 63 |4 |10940 |659 |3.0 |40.0 |57.288 64 |4 |11516 |677 |3.0 |40.0 |57.288 65 |4 |12106 |694 |3.0 |40.0 |57.288 66 |4 |12711 |711 |3.0 |40.0 |57.288 67 |4 |13332 |728 |3.0 |40.0 |57.288 68 |4 |13966 |745 |3.0 |40.0 |57.288 69 |4 |14616 |762 |3.0 |40.0 |57.288 70 |3 |15280 |779 |3.0 |40.0 |57.288 71 |3 |15959 |797 |3.0 |40.0 |57.288 72 |3 |16653 |814 |3.0 |40.0 |57.288 73 |3 |17362 |831 |3.0 |40.0 |57.288 74 |3 |18085 |848 |3.0 |40.0 |57.288 75 |12 |10 |0 |8.0 |12.0 |46.152 76 |4 |78 |42 |8.0 |12.0 |46.152 77 |4 |153 |59 |8.0 |12.0 |46.152 78 |4 |254 |77 |8.0 |12.0 |46.152 79 |4 |379 |94 |8.0 |12.0 |46.152 80 |4 |530 |111 |8.0 |12.0 |46.152 81 |4 |706 |128 |8.0 |12.0 |46.152 82 |4 |907 |145 |8.0 |12.0 |46.152 83 |4 |1133 |162 |8.0 |12.0 |46.152 84 |4 |1384 |179 |8.0 |12.0 |46.152 85 |4 |1661 |197 |8.0 |12.0 |46.152 86 |4 |1962 |214 |8.0 |12.0 |46.152 87 |4 |2289 |231 |8.0 |12.0 |46.152 88 |4 |2641 |248 |8.0 |12.0 |46.152 89 |4 |3017 |265 |8.0 |12.0 |46.152 90 |4 |3419 |282 |8.0 |12.0 |46.152 91 |4 |3846 |299 |8.0 |12.0 |46.152 92 |4 |4299 |317 |8.0 |12.0 |46.152 93 |4 |4776 |334 |8.0 |12.0 |46.152 94 |4 |5278 |351 |8.0 |12.0 |46.152 95 |4 |5806 |368 |8.0 |12.0 |46.152 96 |4 |6359 |385 |8.0 |12.0 |46.152 97 |4 |6937 |402 |8.0 |12.0 |46.152 98 |4 |7539 |419 |8.0 |12.0 |46.152 99 |4 |8167 |437 |8.0 |12.0 |46.152 100 |4 |8820 |454 |8.0 |12.0 |46.152 101 |4 |9499 |471 |8.0 |12.0 |46.152 102 |4 |10202 |488 |8.0 |12.0 |46.152 103 |4 |10931 |505 |8.0 |12.0 |46.152 104 |4 |11684 |522 |8.0 |12.0 |46.152 105 |4 |12463 |539 |8.0 |12.0 |46.152 106 |4 |13267 |557 |8.0 |12.0 |46.152 107 |4 |14096 |574 |8.0 |12.0 |46.152 108 |4 |14950 |591 |8.0 |12.0 |46.152 109 |4 |15830 |608 |8.0 |12.0 |46.152 110 |4 |16734 |625 |8.0 |12.0 |46.152 111 |4 |17664 |642 |8.0 |12.0 |46.152 112 |4 |18618 |659 |8.0 |12.0 |46.152 113 |4 |19598 |677 |8.0 |12.0 |46.152 114 |4 |20603 |694 |8.0 |12.0 |46.152 115 |4 |21633 |711 |8.0 |12.0 |46.152 116 |4 |22688 |728 |8.0 |12.0 |46.152 117 |4 |23769 |745 |8.0 |12.0 |46.152 118 |4 |24874 |762 |8.0 |12.0 |46.152 119 |3 |26004 |779 |8.0 |12.0 |46.152 120 |3 |27160 |797 |8.0 |12.0 |46.152 121 |3 |28341 |814 |8.0 |12.0 |46.152 122 |3 |29547 |831 |8.0 |12.0 |46.152 123 |3 |30778 |848 |8.0 |12.0 |46.152 124 |12 |13 |0 |8.0 |15.0 |46.152 125 |4 |103 |42 |8.0 |15.0 |46.152 126 |4 |203 |59 |8.0 |15.0 |46.152 127 |4 |336 |77 |8.0 |15.0 |46.152 128 |4 |502 |94 |8.0 |15.0 |46.152 129 |4 |701 |111 |8.0 |15.0 |46.152 130 |4 |933 |128 |8.0 |15.0 |46.152 131 |4 |1199 |145 |8.0 |15.0 |46.152 132 |4 |1498 |162 |8.0 |15.0 |46.152 133 |4 |1829 |179 |8.0 |15.0 |46.152 134 |4 |2195 |197 |8.0 |15.0 |46.152 135 |4 |2593 |214 |8.0 |15.0 |46.152 136 |4 |3025 |231 |8.0 |15.0 |46.152 137 |4 |3489 |248 |8.0 |15.0 |46.152 138 |4 |3987 |265 |8.0 |15.0 |46.152 139 |4 |4518 |282 |8.0 |15.0 |46.152 140 |4 |5083 |299 |8.0 |15.0 |46.152 141 |4 |5681 |317 |8.0 |15.0 |46.152 142 |4 |6311 |334 |8.0 |15.0 |46.152 143 |4 |6975 |351 |8.0 |15.0 |46.152 144 |4 |7672 |368 |8.0 |15.0 |46.152 145 |4 |8403 |385 |8.0 |15.0 |46.152 146 |4 |9166 |402 |8.0 |15.0 |46.152 147 |4 |9963 |419 |8.0 |15.0 |46.152 148 |4 |10793 |437 |8.0 |15.0 |46.152 149 |4 |11656 |454 |8.0 |15.0 |46.152 150 |4 |12553 |471 |8.0 |15.0 |46.152 151 |4 |13482 |488 |8.0 |15.0 |46.152 152 |4 |14445 |505 |8.0 |15.0 |46.152 153 |4 |15440 |522 |8.0 |15.0 |46.152 154 |4 |16469 |539 |8.0 |15.0 |46.152 155 |4 |17532 |557 |8.0 |15.0 |46.152 156 |4 |18628 |574 |8.0 |15.0 |46.152 157 |4 |19756 |591 |8.0 |15.0 |46.152 158 |4 |20918 |608 |8.0 |15.0 |46.152 159 |4 |22113 |625 |8.0 |15.0 |46.152 160 |4 |23342 |642 |8.0 |15.0 |46.152 161 |4 |24603 |659 |8.0 |15.0 |46.152 162 |4 |25898 |677 |8.0 |15.0 |46.152 163 |4 |27225 |694 |8.0 |15.0 |46.152 164 |4 |28586 |711 |8.0 |15.0 |46.152 165 |4 |29981 |728 |8.0 |15.0 |46.152 166 |4 |31409 |745 |8.0 |15.0 |46.152 167 |4 |32869 |762 |8.0 |15.0 |46.152 168 |3 |34363 |779 |8.0 |15.0 |46.152 169 |3 |35890 |797 |8.0 |15.0 |46.152 170 |3 |37451 |814 |8.0 |15.0 |46.152 171 |3 |39044 |831 |8.0 |15.0 |46.152 172 |3 |40670 |848 |8.0 |15.0 |46.152 173 |12 |19 |0 |8.0 |20.0 |46.152 174 |4 |145 |42 |8.0 |20.0 |46.152 175 |4 |285 |59 |8.0 |20.0 |46.152 176 |4 |472 |77 |8.0 |20.0 |46.152 177 |4 |705 |94 |8.0 |20.0 |46.152 178 |4 |985 |111 |8.0 |20.0 |46.152 179 |4 |1312 |128 |8.0 |20.0 |46.152 180 |4 |1685 |145 |8.0 |20.0 |46.152 181 |4 |2105 |162 |8.0 |20.0 |46.152 182 |4 |2571 |179 |8.0 |20.0 |46.152 183 |4 |3085 |197 |8.0 |20.0 |46.152 184 |4 |3645 |214 |8.0 |20.0 |46.152 185 |4 |4251 |231 |8.0 |20.0 |46.152 186 |4 |4904 |248 |8.0 |20.0 |46.152 187 |4 |5604 |265 |8.0 |20.0 |46.152 188 |4 |6350 |282 |8.0 |20.0 |46.152 189 |4 |7144 |299 |8.0 |20.0 |46.152 190 |4 |7984 |317 |8.0 |20.0 |46.152 191 |4 |8870 |334 |8.0 |20.0 |46.152 192 |4 |9803 |351 |8.0 |20.0 |46.152 193 |4 |10783 |368 |8.0 |20.0 |46.152 194 |4 |11810 |385 |8.0 |20.0 |46.152 195 |4 |12883 |402 |8.0 |20.0 |46.152 196 |4 |14002 |419 |8.0 |20.0 |46.152 197 |4 |15168 |437 |8.0 |20.0 |46.152 198 |4 |16381 |454 |8.0 |20.0 |46.152 199 |4 |17642 |471 |8.0 |20.0 |46.152 200 |4 |18948 |488 |8.0 |20.0 |46.152 201 |4 |20301 |505 |8.0 |20.0 |46.152 202 |4 |21700 |522 |8.0 |20.0 |46.152 203 |4 |23146 |539 |8.0 |20.0 |46.152 204 |4 |24640 |557 |8.0 |20.0 |46.152 205 |4 |26179 |574 |8.0 |20.0 |46.152 206 |4 |27766 |591 |8.0 |20.0 |46.152 207 |4 |29398 |608 |8.0 |20.0 |46.152 208 |4 |31078 |625 |8.0 |20.0 |46.152 209 |4 |32805 |642 |8.0 |20.0 |46.152 210 |4 |34578 |659 |8.0 |20.0 |46.152 211 |4 |36397 |677 |8.0 |20.0 |46.152 212 |4 |38263 |694 |8.0 |20.0 |46.152 213 |4 |40176 |711 |8.0 |20.0 |46.152 214 |4 |42136 |728 |8.0 |20.0 |46.152 215 |4 |44142 |745 |8.0 |20.0 |46.152 216 |4 |46195 |762 |8.0 |20.0 |46.152 217 |3 |48294 |779 |8.0 |20.0 |46.152 218 |3 |50440 |797 |8.0 |20.0 |46.152 219 |3 |52634 |814 |8.0 |20.0 |46.152 220 |3 |54873 |831 |8.0 |20.0 |46.152 221 |3 |57159 |848 |8.0 |20.0 |46.152 222 |8 |3 |0 |8.0 |25.0 |47.288 223 |4 |67 |25 |8.0 |25.0 |47.288 224 |4 |187 |42 |8.0 |25.0 |47.288 225 |4 |368 |59 |8.0 |25.0 |47.288 226 |4 |608 |77 |8.0 |25.0 |47.288 227 |4 |909 |94 |8.0 |25.0 |47.288 228 |4 |1269 |111 |8.0 |25.0 |47.288 229 |4 |1690 |128 |8.0 |25.0 |47.288 230 |4 |2171 |145 |8.0 |25.0 |47.288 231 |4 |2712 |162 |8.0 |25.0 |47.288 232 |4 |3313 |179 |8.0 |25.0 |47.288 233 |4 |3974 |197 |8.0 |25.0 |47.288 234 |4 |4696 |214 |8.0 |25.0 |47.288 235 |4 |5478 |231 |8.0 |25.0 |47.288 236 |4 |6319 |248 |8.0 |25.0 |47.288 237 |4 |7221 |265 |8.0 |25.0 |47.288 238 |4 |8182 |282 |8.0 |25.0 |47.288 239 |4 |9205 |299 |8.0 |25.0 |47.288 240 |4 |10287 |317 |8.0 |25.0 |47.288 241 |4 |11429 |334 |8.0 |25.0 |47.288 242 |4 |12631 |351 |8.0 |25.0 |47.288 243 |4 |13893 |368 |8.0 |25.0 |47.288 244 |4 |15216 |385 |8.0 |25.0 |47.288 245 |4 |16599 |402 |8.0 |25.0 |47.288 246 |4 |18041 |419 |8.0 |25.0 |47.288 247 |4 |19544 |437 |8.0 |25.0 |47.288 248 |4 |21107 |454 |8.0 |25.0 |47.288 249 |4 |22731 |471 |8.0 |25.0 |47.288 250 |4 |24414 |488 |8.0 |25.0 |47.288 251 |4 |26157 |505 |8.0 |25.0 |47.288 252 |4 |27960 |522 |8.0 |25.0 |47.288 253 |4 |29823 |539 |8.0 |25.0 |47.288 254 |4 |31748 |557 |8.0 |25.0 |47.288 255 |4 |33731 |574 |8.0 |25.0 |47.288 256 |4 |35775 |591 |8.0 |25.0 |47.288 257 |4 |37879 |608 |8.0 |25.0 |47.288 258 |4 |40043 |625 |8.0 |25.0 |47.288 259 |4 |42268 |642 |8.0 |25.0 |47.288 260 |4 |44552 |659 |8.0 |25.0 |47.288 261 |4 |46896 |677 |8.0 |25.0 |47.288 262 |4 |49300 |694 |8.0 |25.0 |47.288 263 |4 |51765 |711 |8.0 |25.0 |47.288 264 |4 |54291 |728 |8.0 |25.0 |47.288 265 |4 |56875 |745 |8.0 |25.0 |47.288 266 |4 |59520 |762 |8.0 |25.0 |47.288 267 |3 |62225 |779 |8.0 |25.0 |47.288 268 |3 |64990 |797 |8.0 |25.0 |47.288 269 |3 |67817 |814 |8.0 |25.0 |47.288 270 |3 |70702 |831 |8.0 |25.0 |47.288 271 |3 |73647 |848 |8.0 |25.0 |47.288 272 |8 |4 |0 |8.0 |30.0 |52.288 273 |4 |82 |25 |8.0 |30.0 |52.288 274 |4 |229 |42 |8.0 |30.0 |52.288 275 |4 |450 |59 |8.0 |30.0 |52.288 276 |4 |744 |77 |8.0 |30.0 |52.288 277 |4 |1112 |94 |8.0 |30.0 |52.288 278 |4 |1554 |111 |8.0 |30.0 |52.288 279 |4 |2069 |128 |8.0 |30.0 |52.288 280 |4 |2657 |145 |8.0 |30.0 |52.288 281 |4 |3319 |162 |8.0 |30.0 |52.288 282 |4 |4055 |179 |8.0 |30.0 |52.288 283 |4 |4864 |197 |8.0 |30.0 |52.288 284 |4 |5748 |214 |8.0 |30.0 |52.288 285 |4 |6704 |231 |8.0 |30.0 |52.288 286 |4 |7734 |248 |8.0 |30.0 |52.288 287 |4 |8837 |265 |8.0 |30.0 |52.288 288 |4 |10014 |282 |8.0 |30.0 |52.288 289 |4 |11266 |299 |8.0 |30.0 |52.288 290 |4 |12590 |317 |8.0 |30.0 |52.288 291 |4 |13988 |334 |8.0 |30.0 |52.288 292 |4 |15459 |351 |8.0 |30.0 |52.288 293 |4 |17004 |368 |8.0 |30.0 |52.288 294 |4 |18623 |385 |8.0 |30.0 |52.288 295 |4 |20315 |402 |8.0 |30.0 |52.288 296 |4 |22081 |419 |8.0 |30.0 |52.288 297 |4 |23920 |437 |8.0 |30.0 |52.288 298 |4 |25832 |454 |8.0 |30.0 |52.288 299 |4 |27820 |471 |8.0 |30.0 |52.288 300 |4 |29879 |488 |8.0 |30.0 |52.288 301 |4 |32013 |505 |8.0 |30.0 |52.288 302 |4 |34220 |522 |8.0 |30.0 |52.288 303 |4 |36500 |539 |8.0 |30.0 |52.288 304 |4 |38856 |557 |8.0 |30.0 |52.288 305 |4 |41283 |574 |8.0 |30.0 |52.288 306 |4 |43784 |591 |8.0 |30.0 |52.288 307 |4 |46359 |608 |8.0 |30.0 |52.288 308 |4 |49007 |625 |8.0 |30.0 |52.288 309 |4 |51731 |642 |8.0 |30.0 |52.288 310 |4 |54526 |659 |8.0 |30.0 |52.288 311 |4 |57395 |677 |8.0 |30.0 |52.288 312 |4 |60338 |694 |8.0 |30.0 |52.288 313 |4 |63354 |711 |8.0 |30.0 |52.288 314 |4 |66446 |728 |8.0 |30.0 |52.288 315 |4 |69609 |745 |8.0 |30.0 |52.288 316 |4 |72846 |762 |8.0 |30.0 |52.288 317 |3 |76156 |779 |8.0 |30.0 |52.288 318 |3 |79540 |797 |8.0 |30.0 |52.288 319 |3 |82999 |814 |8.0 |30.0 |52.288 320 |3 |86531 |831 |8.0 |30.0 |52.288 321 |3 |90135 |848 |8.0 |30.0 |52.288 322 |8 |5 |0 |8.0 |35.0 |57.288 323 |4 |97 |25 |8.0 |35.0 |57.288 324 |4 |271 |42 |8.0 |35.0 |57.288 325 |4 |533 |59 |8.0 |35.0 |57.288 326 |4 |881 |77 |8.0 |35.0 |57.288 327 |4 |1316 |94 |8.0 |35.0 |57.288 328 |4 |1838 |111 |8.0 |35.0 |57.288 329 |4 |2447 |128 |8.0 |35.0 |57.288 330 |4 |3144 |145 |8.0 |35.0 |57.288 331 |4 |3927 |162 |8.0 |35.0 |57.288 332 |4 |4797 |179 |8.0 |35.0 |57.288 333 |4 |5754 |197 |8.0 |35.0 |57.288 334 |4 |6799 |214 |8.0 |35.0 |57.288 335 |4 |7930 |231 |8.0 |35.0 |57.288 336 |4 |9149 |248 |8.0 |35.0 |57.288 337 |4 |10454 |265 |8.0 |35.0 |57.288 338 |4 |11846 |282 |8.0 |35.0 |57.288 339 |4 |13327 |299 |8.0 |35.0 |57.288 340 |4 |14893 |317 |8.0 |35.0 |57.288 341 |4 |16546 |334 |8.0 |35.0 |57.288 342 |4 |18287 |351 |8.0 |35.0 |57.288 343 |4 |20114 |368 |8.0 |35.0 |57.288 344 |4 |22030 |385 |8.0 |35.0 |57.288 345 |4 |24031 |402 |8.0 |35.0 |57.288 346 |4 |26120 |419 |8.0 |35.0 |57.288 347 |4 |28295 |437 |8.0 |35.0 |57.288 348 |4 |30558 |454 |8.0 |35.0 |57.288 349 |4 |32909 |471 |8.0 |35.0 |57.288 350 |4 |35345 |488 |8.0 |35.0 |57.288 351 |4 |37869 |505 |8.0 |35.0 |57.288 352 |4 |40480 |522 |8.0 |35.0 |57.288 353 |4 |43177 |539 |8.0 |35.0 |57.288 354 |4 |45963 |557 |8.0 |35.0 |57.288 355 |4 |48835 |574 |8.0 |35.0 |57.288 356 |4 |51794 |591 |8.0 |35.0 |57.288 357 |4 |54840 |608 |8.0 |35.0 |57.288 358 |4 |57972 |625 |8.0 |35.0 |57.288 359 |4 |61194 |642 |8.0 |35.0 |57.288 360 |4 |64501 |659 |8.0 |35.0 |57.288 361 |4 |67895 |677 |8.0 |35.0 |57.288 362 |4 |71375 |694 |8.0 |35.0 |57.288 363 |4 |74943 |711 |8.0 |35.0 |57.288 364 |4 |78600 |728 |8.0 |35.0 |57.288 365 |4 |82342 |745 |8.0 |35.0 |57.288 366 |4 |86171 |762 |8.0 |35.0 |57.288 367 |3 |90087 |779 |8.0 |35.0 |57.288 368 |3 |94090 |797 |8.0 |35.0 |57.288 369 |3 |98182 |814 |8.0 |35.0 |57.288 370 |3 |102359 |831 |8.0 |35.0 |57.288 371 |3 |106623 |848 |8.0 |35.0 |57.288 372 |8 |6 |0 |8.0 |40.0 |62.288 373 |4 |113 |25 |8.0 |40.0 |62.288 374 |4 |314 |42 |8.0 |40.0 |62.288 375 |4 |615 |59 |8.0 |40.0 |62.288 376 |4 |1017 |77 |8.0 |40.0 |62.288 377 |4 |1519 |94 |8.0 |40.0 |62.288 378 |4 |2122 |111 |8.0 |40.0 |62.288 379 |4 |2826 |128 |8.0 |40.0 |62.288 380 |4 |3630 |145 |8.0 |40.0 |62.288 381 |4 |4534 |162 |8.0 |40.0 |62.288 382 |4 |5539 |179 |8.0 |40.0 |62.288 383 |4 |6644 |197 |8.0 |40.0 |62.288 384 |4 |7851 |214 |8.0 |40.0 |62.288 385 |4 |9157 |231 |8.0 |40.0 |62.288 386 |4 |10564 |248 |8.0 |40.0 |62.288 387 |4 |12071 |265 |8.0 |40.0 |62.288 388 |4 |13678 |282 |8.0 |40.0 |62.288 389 |4 |15387 |299 |8.0 |40.0 |62.288 390 |4 |17196 |317 |8.0 |40.0 |62.288 391 |4 |19105 |334 |8.0 |40.0 |62.288 392 |4 |21115 |351 |8.0 |40.0 |62.288 393 |4 |23225 |368 |8.0 |40.0 |62.288 394 |4 |25436 |385 |8.0 |40.0 |62.288 395 |4 |27748 |402 |8.0 |40.0 |62.288 396 |4 |30159 |419 |8.0 |40.0 |62.288 397 |4 |32671 |437 |8.0 |40.0 |62.288 398 |4 |35283 |454 |8.0 |40.0 |62.288 399 |4 |37998 |471 |8.0 |40.0 |62.288 400 |4 |40811 |488 |8.0 |40.0 |62.288 401 |4 |43725 |505 |8.0 |40.0 |62.288 402 |4 |46739 |522 |8.0 |40.0 |62.288 403 |4 |49854 |539 |8.0 |40.0 |62.288 404 |4 |53071 |557 |8.0 |40.0 |62.288 405 |4 |56387 |574 |8.0 |40.0 |62.288 406 |4 |59803 |591 |8.0 |40.0 |62.288 407 |4 |63320 |608 |8.0 |40.0 |62.288 408 |4 |66937 |625 |8.0 |40.0 |62.288 409 |4 |70657 |642 |8.0 |40.0 |62.288 410 |4 |74475 |659 |8.0 |40.0 |62.288 411 |4 |78394 |677 |8.0 |40.0 |62.288 412 |4 |82413 |694 |8.0 |40.0 |62.288 413 |4 |86533 |711 |8.0 |40.0 |62.288 414 |4 |90755 |728 |8.0 |40.0 |62.288 415 |4 |95076 |745 |8.0 |40.0 |62.288 416 |4 |99497 |762 |8.0 |40.0 |62.288 417 |3 |104018 |779 |8.0 |40.0 |62.288 418 |3 |108640 |797 |8.0 |40.0 |62.288 419 |3 |113365 |814 |8.0 |40.0 |62.288 420 |3 |118188 |831 |8.0 |40.0 |62.288 421 |3 |123112 |848 |8.0 |40.0 |62.288
The 2D data only has data measured in the xy-plane, where we only use the [1, 0, 0] dummy direction for b0 measurements. We can see the unique gradient directions as follows:
In [7]:
import numpy as np
unique_directions = np.unique(scheme_spinal_cord.gradient_directions, axis=0)
unique_directions
Out[7]:
array([[-0.70711 , -0.70711 , 0. ],
[-0.70711 , 0.70711 , 0. ],
[-0.707107, -0.707107, 0. ],
[-0.707107, 0.707107, 0. ],
[ 0.707107, -0.707107, 0. ],
[ 0.707107, 0.707107, 0. ],
[ 0.70711 , -0.70711 , 0. ],
[ 0.70711 , 0.70711 , 0. ],
[ 1. , 0. , 0. ]])
It is valuable to check if actually the signal between the two orthogonal directions is the same:
In [8]:
dir1mask = np.all(scheme_spinal_cord.gradient_directions == unique_directions[0], axis=1)
dir2mask = np.all(scheme_spinal_cord.gradient_directions == unique_directions[1], axis=1)
test_voxel = data_spinal_cord.signal[20, 20, 0]
data_dir1 = test_voxel[dir1mask]
data_dir2 = test_voxel[dir2mask]
Delta_dir1 = scheme_spinal_cord.Delta[dir1mask]
Delta_dir2 = scheme_spinal_cord.Delta[dir2mask]
G = scheme_spinal_cord.gradient_strengths
colors = ['r', 'g', 'b', 'orange']
unique_Deltas = np.unique(scheme_spinal_cord.shell_Delta)[2:-1:2]
for i, Delta_ in enumerate(unique_Deltas):
mask_Delta1 = Delta_dir1 == Delta_
mask_Delta2 = Delta_dir2 == Delta_
plt.plot(G[dir1mask][mask_Delta1], data_dir1[mask_Delta1], c=colors[i], label='dir1 Delta='+str(Delta_))
plt.plot(G[dir2mask][mask_Delta2], data_dir2[mask_Delta2], c=colors[i], ls='--', label='dir2 Delta='+str(Delta_))
plt.legend()
plt.xlabel('Gradient Strength [T/m]')
plt.ylabel('Signal [-]');
It can be seen that the two measured perpendicular directions experience slightly different levels of diffusion restriction; the dir2 signal profiles attenuate faster than the dir1 ones. This could implies that that some kind of anisotropic axon dispersion could be happening, or that the axon bundle is slightly tilted towards dir2.
Dmipy allows us to add axon dispersion into the AxCaliber model to compensate for this, but we'll ignore that for this example and fit the signal as-is.
Fit AxCaliber to spinal cord data¶
It is known that a Gamma distribution is tricky to fit since there are several parameter combinations that produce similar distributions. We will use the stochastic MIX optimization (Farooq et al. 2016), which is slow but is likely to robustly give the global minimum.
NOTE: Fitting AxCaliber to this dataset with MIX takes over 7 hours!¶
In [9]:
axcaliber_gamma_fit = axcaliber_gamma.fit(
scheme_spinal_cord, data_spinal_cord.signal,
mask=data_spinal_cord.mask, solver='mix', maxiter=100)
Using parallel processing with 8 workers. Setup MIX optimizer in 4.57763671875e-05 seconds Fitting of 968 voxels complete in 30356.3031921 seconds. Average of 31.3598173473 seconds per voxel.
Visualize fitted parameter maps¶
We can visualize the fitted parameters as before:
In [10]:
fitted_parameters = axcaliber_gamma_fit.fitted_parameters
fig, axs = plt.subplots(2, 3, figsize=[15, 10])
axs = axs.ravel()
counter = 0
for name, values in fitted_parameters.items():
cf = axs[counter].imshow(values.squeeze().T, origin=True, interpolation='nearest')
axs[counter].set_title(name)
axs[counter].set_axis_off()
fig.colorbar(cf, ax=axs[counter], shrink=0.5)
counter += 1
We can see that the maps for alpha and beta are not super smooth. This could be because we're actually fitting two signals experiencing different restriction levels at the same time, but hard to say without fitting each signal separately. We can visualize the distributions themselves as follows:
In [15]:
axcaliber_alpha = axcaliber_gamma_fit.fitted_parameters['DD1GammaDistributed_1_DD1Gamma_1_alpha']
axcaliber_beta = axcaliber_gamma_fit.fitted_parameters['DD1GammaDistributed_1_DD1Gamma_1_beta']
from dmipy.distributions.distributions import DD1Gamma
gamma = DD1Gamma()
x = 15 # some x-pos in the data
y = 25 # some y-pos in the data
alpha = axcaliber_alpha[x, y, 0]
beta = axcaliber_beta[x, y, 0]
radii, Pgamma = gamma(alpha=alpha, beta=beta)
plt.plot(2 * radii * 1e6, Pgamma)
plt.xlabel('cylinder diameter [$\mu$m]')
plt.ylabel('Gamma probability density')
plt.title('Estimated Gamma distribution');
Comparison between estimated and histology axon diameter¶
The mean axon diameter can be estimated from the axon diameter distribution as diameter = 2$\alpha\beta$. We compare the histological axon diameter with the estimated ones:
In [16]:
mean_diameter = 2 * axcaliber_alpha * axcaliber_beta * 1e6 # in microns
mean_diameter_histology = data_spinal_cord.histology.h1_axonEquivDiameter # was already in microns
In [17]:
fig, [ax1, ax2] = plt.subplots(1, 2, figsize=[10, 5])
ax1plot = ax1.imshow(mean_diameter.squeeze().T,
origin=True, vmax=6, interpolation='nearest')
ax1.set_title('Estimated Mean Axon Diameter')
ax1.set_axis_off()
plt.colorbar(ax1plot, ax=ax1, shrink=0.75)
ax2plot = ax2.imshow(mean_diameter_histology.T,
origin=True, vmax=6, interpolation='nearest')
ax2.set_title('Histology Mean Axon Diameter')
ax2.set_axis_off()
plt.colorbar(ax2plot, ax=ax2, shrink=0.75)
Out[17]:
<matplotlib.colorbar.Colorbar at 0x7f56346bea10>
We can see that our estimated mean axon diameter map is completely smooth, aside from a few outliers in areas which actually lie outside of the white matter area (see right histology image). The estimated diameters are overestimated in general. We find between 4-6 $\mu$m, while the histological values are between 0-5 $\mu$m.
Using the Seaborn package, we can nicely visualize the correlation between our estimated and the histology diameters. We separate visualize this correlation of axon diameters between 1-2$\mu$m, 2-3$\mu$m and 3+$\mu$m, and remove the exceptional outliers.
In [18]:
mean_diameter_mask12 = np.all([mean_diameter_histology > 1,
mean_diameter_histology < 2,
mean_diameter.squeeze() < 10], axis=0)
mean_diameter_mask23 = np.all([mean_diameter_histology > 2,
mean_diameter_histology < 3,
mean_diameter.squeeze() < 10], axis=0)
mean_diameter_mask3p = np.all([mean_diameter_histology > 3,
mean_diameter.squeeze() < 10], axis=0)
Axon diameter correlation 1-2 $\mu$m - not significant¶
In [19]:
import seaborn as sns
mask = mean_diameter_mask12
im = sns.jointplot(mean_diameter[mask].squeeze(),
mean_diameter_histology[mask].squeeze(), kind="reg")
im.set_axis_labels("Estimated Axon Diameter", "Histology Axon Diameter", fontsize=15)
Out[19]:
<seaborn.axisgrid.JointGrid at 0x7f56139c7550>
Axon diameter correlation 2-3 $\mu$m - not significant¶
In [20]:
mask = mean_diameter_mask23
im = sns.jointplot(mean_diameter[mask].squeeze(),
mean_diameter_histology[mask].squeeze(), kind="reg", color='r')
im.set_axis_labels("Estimated Axon Diameter", "Histology Axon Diameter", fontsize=15)
Out[20]:
<seaborn.axisgrid.JointGrid at 0x7f5634a06a10>
Axon diameter correlation 3+ $\mu$m - Significant¶
In [21]:
mask = mean_diameter_mask3p
im = sns.jointplot(mean_diameter[mask].squeeze(),
mean_diameter_histology[mask].squeeze(), kind="reg", color='g')
im.set_axis_labels("Estimated Axon Diameter", "Histology Axon Diameter", fontsize=15)
Out[21]:
<seaborn.axisgrid.JointGrid at 0x7f56138db490>
The results show that for our current approach to the data set we can only estimate axon diameters that correlate with the histology values for axon diameters of 3+$\mu$m.
Comparison between estimated and histology volume fraction¶
Similarly, we can also visualize the estimated intra-axonal volume fraction between our model and the histology values.
In [22]:
axcaliber_fraction = axcaliber_gamma_fit.fitted_parameters['partial_volume_1'].squeeze()
histology_fraction = data_spinal_cord.histology.h4_fr
In [23]:
fig, [ax1, ax2] = plt.subplots(1, 2, figsize=[10, 5])
ax1plot = ax1.imshow(axcaliber_fraction.T,
origin=True,interpolation='nearest', vmax=0.55)
ax1.set_title('Estimated Volume Fraction')
ax1.set_axis_off()
plt.colorbar(ax1plot, ax=ax1, shrink=0.75)
ax2plot = ax2.imshow(histology_fraction.T,
origin=True, interpolation='nearest', vmax=0.55)
ax2.set_title('Histology Volume Fraction')
ax2.set_axis_off()
plt.colorbar(ax2plot, ax=ax2, shrink=0.75)
Out[23]:
<matplotlib.colorbar.Colorbar at 0x7f561330e450>
In [24]:
im = sns.jointplot(axcaliber_fraction[mean_diameter_histology > 1],
histology_fraction[mean_diameter_histology > 1], kind="reg", color='y')
im.set_axis_labels("Estimated Volume Fraction", "Histology Volume Fraction", fontsize=15);
The results show that estimating intra-axonal volume fractions that correlate with histology is easier than doing the same for axon diameter. However, the model does systematically underestimate the value compared to histology.
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.
- 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.
- Farooq, Hamza, et al. "Microstructure imaging of crossing (MIX) white matter fibers from diffusion MRI." Scientific reports 6 (2016): 38927.
- Huang, Susie Y., et al. "The impact of gradient strength on in vivo diffusion MRI estimates of axon diameter." NeuroImage 106 (2015): 464-472.