forked from mavel101/AxonRadiusMapping
-
Notifications
You must be signed in to change notification settings - Fork 0
/
getAxonRadius.m
126 lines (94 loc) · 4.43 KB
/
getAxonRadius.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
function [r, beta] = getAxonRadius(delta, Delta, g, y, model)
% "getAxonRadius": Estimate the axon radius from multi-shell diffusion
% data, of which the lowest b-value is sufficiently high to suppress
% the extra-axonal compartment. For in vivo human MRI, we recommend
% b=6000s/mm2.
%
% [r, beta] = getAxonRadius(delta, Delta, g, y, model)
% output:
% - r: effective MR radius
% - beta: [1x1] The slope of the signal scaling as a function
% of 1/sqrt. ~ f/sqrt(Da)
% input:
% - delta: gradient duration [ms]
% - Delta: gradient seperation [ms]
% - g: gradient strenght [mT/m]
% - y: powder-averaged diffusion-weighted signal,
% normalized to y(b=0) = 1. If the dot comparment
% cannot be ignored, it need to be subtracted from y prior to fitting.
% - model: Van Gelderen or Neumann (the latter only applied in long pulse regimes, but typically applies)
%
% Author: Jelle Veraart (jelle.veraart@nyulangone.org)
% Copyright (c) 2019 New York University
if ~exist('model', 'var')
model = 'VanGelderen';
end
options = optimset('lsqnonlin');
options = optimset(options,'Jacobian','on','TolFun',1e-12,'TolX',1e-12,'MaxIter',10000,'Display','off');
gyroMagnRatio = 267.513*10^(-6);
q = g*gyroMagnRatio;
b = (q.*delta).^2.*(Delta - delta/3);
% Random initial guess seems to work well for Neuman, but not for van
% Gelderen model, so we use the Neumann results as initial guesses for
% the van Gelderen model
start = [sqrt(4*pi)*rand(), 1.5+2*rand()];
pars = lsqnonlin(@(x)residuals(x, [delta(:), Delta(:), g(:)], y, 'Neumann'),start,[0 0],[sqrt(4*pi) 5],options);
if strcmp(model,'VanGelderen')
pars = lsqnonlin(@(x)residuals(x, [delta(:), Delta(:), g(:)], y, 'VanGelderen'),pars,[0 0],[sqrt(4*pi) 5],options);
end
beta = pars(1);
r = pars(2);
end
function [E, J] = residuals(pars, X, y, model)
delta = X(:, 1);
Delta = X(:, 2);
g = X(:, 3);
[shat, dshat] = AxonDiameterFWD(delta, Delta, g, pars, model);
E = shat(:) - y(:);
J = dshat;
end
function [s, ds] = AxonDiameterFWD(delta, Delta, g, pars, model)
% read parameters with fixed D0
D0 = 2.5;
gyroMagnRatio = 267.513*10^(-6);
beta = pars(1);
r = pars(2);
% compute b-values using Tjeskall-Tanner sequence
q = g*gyroMagnRatio;
b = (q.*delta).^2.*(Delta - delta/3);
% Forward model
switch model
case 'Neumann'
s = beta*exp(-(7/48)*q.^2.*delta*r^4/D0)./sqrt(b);
ds_dbeta = exp(-(7/48)*q.^2.*delta*r^4/D0) ./ sqrt(b);
ds_dr = beta ./ sqrt(b) .* exp(-(7/48)*q.^2.*delta*r^4/D0) .*(-(7/12)*q.^2.*delta*r^3/D0);
case 'VanGelderen'
[Svg, dSvg] = vg(delta, Delta, q, r, D0);
s = beta*exp(-Svg) ./ sqrt(b);
ds_dbeta = exp(-Svg) ./ sqrt(b);
ds_dr = beta*exp(-Svg) ./ sqrt(b) .* -dSvg;
end
ds = [ds_dbeta, ds_dr];
end
function [s, ds] = vg(delta, Delta, q, r, D0)
td=r^2/D0;
bardelta=delta/td; dbardelta = -2*delta*D0 / r^3;
barDelta=Delta/td; dbarDelta = -2*Delta*D0 / r^3;
N=15;
b = [1.8412 5.3314 8.5363 11.7060 14.8636 18.0155 21.1644 24.3113 27.4571 30.6019 ...
33.7462 36.8900 40.0334 43.1766 46.3196 49.4624 52.6050 55.7476 58.8900 62.0323];
s = 0; ds = 0;
for k=1:N
s = s + (2/(b(k)^6*(b(k)^2-1)))*(-2 + 2*b(k)^2*bardelta + ...
2*(exp(-b(k)^2*bardelta)+exp(-b(k)^2*barDelta)) - ...
exp(-b(k)^2*(bardelta+barDelta)) - exp(-b(k)^2*(barDelta-bardelta)));
ds = ds + (2/(b(k)^6*(b(k)^2-1)))*( ...
2*b(k)^2*dbardelta + ...
2*exp(-b(k)^2*bardelta) - b(k)^2*dbardelta + ...
2*exp(-b(k)^2*barDelta) - b(k)^2*dbarDelta - ...
exp(-b(k)^2*(bardelta+barDelta)) - b(k)^2*(dbardelta + dbarDelta) - ...
exp(-b(k)^2*(barDelta-bardelta)) - b(k)^2*(dbarDelta - dbardelta));
end
s = s.*D0.*q.^2.*td^3;
ds = ds.*D0.*q.^2.*td^3 + 6*s.*q.^2.*r^5 / D0^2;
end