OSCR

Axon Diameter Mapping From Myelin Water Diffusion MRI.

Code ↔ Paper

16 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 16 matches
  1. [1] § Methods › Diffusion Simulations in Cylindrical Shells With Caliber Variations and Undulations ↔ demo_4_bead_undulation_RD_AD.m, lines 18–156 · score 0.91 · undulation wavelength, bead width, undulation amplitude, axon skeleton, caliber variations, axon length
  2. [2] § Methods › Diffusion Simulations in Cylindrical Shells With Caliber Variations and Undulations ↔ demo_5_bead_undulation_ADM.m, lines 65–107 · score 0.91 · undulation wavelength, bead width, undulation amplitude, axon skeleton, caliber variations, axon length
  3. [3] § Results › Diffusion Simulations of ADM ↔ demo_3_cylinder_shell_ADM.m, lines 493–588 · score 0.74 · log linear, Gamma distribution, axon diameter, inner radius, histology, radii
  4. [4] § Methods › Diffusion Simulations of RD: Two-Dimensional Ring ↔ demo_6_GPA.m, lines 12–102 · score 0.70 · zero thickness, layer thickness, finite thickness, cylindrical surface, cylindrical shell, myelin layer
  5. [5] § Methods › Diffusion Simulations of ADM: Three-Dimensional Shells ↔ demo_3_cylinder_shell_ADM.m, lines 175–246 · score 0.70 · volume weighted, axial diffusivity, narrow pulse solution, gradient directions, wide pulse, spherical
  6. [6] § Methods › Diffusion Simulations of ADM: Three-Dimensional Shells ↔ demo_5_bead_undulation_ADM.m, lines 289–360 · score 0.70 · volume weighted, axial diffusivity, narrow pulse solution, gradient directions, wide pulse, spherical
  7. [7] § Methods › Diffusion Simulations of ADM: Three-Dimensional Shells ↔ demo_3_cylinder_shell_ADM.m, lines 493–588 · score 0.67 · log linear, inner diameter, myelin layers, histological, ratio, outer
  8. [8] § Methods › Diffusion Simulations of ADM: Three-Dimensional Shells ↔ demo_3_cylinder_shell_ADM.m, lines 175–246 · score 0.66 · Rician noise floor, noise ratio, SNR, signals, ADM, Shells
  9. [9] § Methods › Diffusion Simulations of ADM: Three-Dimensional Shells ↔ demo_5_bead_undulation_ADM.m, lines 289–360 · score 0.66 · Rician noise floor, noise ratio, SNR, signals, ADM, Shells
  10. [10] § Methods › Diffusion Simulations in Cylindrical Shells With Caliber Variations and Undulations ↔ lib/rms/main_PGSE_stringbead.cu, lines 484–524 · score 0.65 · periodic boundary, rejection sampling, walkers, membranes, layer
  11. [11] § Methods › Diffusion Simulations in Cylindrical Shells With Caliber Variations and Undulations ↔ demo_4_bead_undulation_RD_AD.m, lines 18–156 · score 0.64 · caliber variations, random walker, intrinsic diffusivity, frustums, undulations, CUDA
  12. [12] § Methods › Diffusion Simulations of RD: One-Dimensional Circle ↔ demo_1_cylindrical_surface_RD.m, lines 9–45 · score 0.60 · Monte Carlo simulations, cylindrical surface, intrinsic diffusivity, walkers, radius
  13. [13] § Results › Diffusion Simulations of RD ↔ demo_6_GPA.m, lines 12–102 · score 0.59 · radius ratio, finite thickness, cylindrical surface, cylindrical shell, narrow pulse, wide pulse
  14. [14] § Methods › Diffusion Simulations of RD: Two-Dimensional Ring ↔ demo_1_cylindrical_surface_RD.m, lines 9–45 · score 0.58 · Monte Carlo simulations, random walker, intrinsic diffusivity, cylindrical, radius
  15. [15] § Methods › Diffusion Simulations of ADM: Three-Dimensional Shells ↔ demo_3_cylinder_shell_ADM.m, lines 18–108 · score 0.56 · permeable membranes, random walkers, intrinsic diffusivity, myelin layer, ADM, Shells
  16. [16] § Methods › Diffusion Simulations of RD: Two-Dimensional Ring ↔ demo_2_cylinder_shell_RD.m, lines 18–105 · score 0.56 · random walker, inner radius, intrinsic diffusivity, myelin layer, CUDA, thickness

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

MATLAB · 590 lines · 18 KB · MIT · 5 matches

  1. % ********** Setup the directory on your computer **********
  2. % demo 3: MC simulations of SMT in coaxial cylindrical shells of finite
  3. % thickness, permeable membrane
  4. clear
  5. restoredefaultpath
  6. filePath = matlab.desktop.editor.getActiveFilename;
  7. root0 = fileparts(filePath);
  8. addpath(genpath(fullfile(root0,'lib')));
  10. root = fullfile(root0,'data');
  11. root_cuda = fullfile(root0,'lib','rms');
  13. % project name
  14. projname = 'coaxial_cyliner_shell_ADM';
  16. mkdir(fullfile(root,projname));
  18. %% Generate packing, permeable membrane
  20. % inner radius, um
  21. r = 0.1:0.1:5;
  23. % thickness of each myelin layer, um
  24. % myelin layer thickness lm should be larger than the step size sqrt(4*D0*dt)
  25. lm = 12/1e3;
  27. % # myelin layers
  28. C0 = 0.35/lm;
  29. C1 = 0.006/lm;
  30. C2 = 0.024/lm;
  31. Nm = round(C0 + C1*2*r + C2*log(2*r));
  33. % permeability, um/ms
  34. kappa = 6.7e-3;
  36. % pulse width, ms
  37. Td = 6;
  39. % diffusion time, ms
  40. TD = 13;
  42. % b-value, ms/um2
  43. bval = [50 350 800 1500 2400 3450 4750 6000]/1e3;
  45. % gradient direction, 64 directions per b-shell
  46. % requirement: MRtrix3
  47. bvec = dirgen(64);
  49. % simulation parameters
  50. dt = 1e-6; % time step, ms
  51. TN = ceil(max(TD+Td)/dt)+100; % # steps
  52. NPar = 1e4; % # random walkers
  53. D0 = 0.8; % intrinsic diffusivity, um2/ms
  54. threadpb = 256; % thread per block for CUDA
  56. seed = 0;
  57. for i = 1:numel(r)
  58. seed = seed + 1;
  59. target = fullfile(root,projname,sprintf('CoCyl_%04u',seed));
  60. mkdir(target);
  62. % field of view, um
  63. res = 2*(r(i)+Nm(i)*lm)*1.05;
  65. % save simulation parameters
  66. fileID = fopen(fullfile(target,'simParamInput.txt'),'w');
  67. fprintf(fileID,sprintf('%g\n%u\n%u\n%g\n%g\n%u\n%g\n%g\n%u\n%g\n',...
  68. dt, TN, NPar, D0, kappa, threadpb, ...
  69. r(i)/res, lm/res, Nm(i), res));
  70. fclose(fileID);
  72. % save diffusion time and pulse width
  73. NDelta = numel(TD);
  74. DELdel = zeros(NDelta,2);
  75. for iii = 1:numel(TD)
  76. DELdel(iii,:) = [TD(iii), Td(iii)];
  77. end
  78. DELdel = DELdel.';
  79. DELdel = DELdel(:);
  80. fid = fopen(fullfile(target,'gradient_NDelta.txt'),'w');
  81. fprintf(fid,sprintf('%u\n',NDelta));
  82. fclose(fid);
  84. fid = fopen(fullfile(target,'gradient_DELdel.txt'),'w');
  85. fprintf(fid,sprintf('%.8f\n',DELdel));
  86. fclose(fid);
  88. % save b-table
  89. ig = 0;
  90. btab = zeros(numel(bval)*size(bvec,1),4);
  91. for jjj = 1:numel(bval)
  92. bvalj = bval(jjj);
  93. for kkk = 1:size(bvec,1)
  94. ig = ig+1;
  95. bveck = bvec(kkk,:);
  96. btab(ig,:) = [bvalj bveck];
  97. end
  98. end
  99. btab = btab.';
  100. btab = btab(:);
  101. fid = fopen(fullfile(target,'gradient_Nbtab.txt'),'w');
  102. fprintf(fid,sprintf('%u\n',numel(bval)*size(bvec,1)));
  103. fclose(fid);
  105. fid = fopen(fullfile(target,'gradient_btab.txt'),'w');
  106. fprintf(fid,sprintf('%.8f\n',btab));
  107. fclose(fid);
  108. end
  110. %% Create a shell script to run the codes
  112. fileID = fopen(fullfile(root,projname,'job.sh'),'w');
  113. fprintf(fileID,'#!/bin/bash\n');
  114. for j = 1:50
  115. target = fullfile(root,projname,sprintf('CoCyl_%04u',j));
  116. fprintf(fileID,sprintf('cd %s\n',target));
  117. fprintf(fileID,sprintf('cp -a %s .\n',fullfile(root_cuda,'main_PGSE_perm_cuda')));
  118. fprintf(fileID,'./main_PGSE_perm_cuda\n');
  119. end
  120. fclose(fileID);
  122. % Open the terminal window in the project folder and run "sh job.sh"
  124. % You may need to open the terminal in the root_cuda folder and compile the
  125. % CUDA code using "nvcc main_PGSE_perm.cu -o main_PGSE_perm_cuda"
  127. %% Plot simulation results
  128. Nr = 50; % # cylinder radius
  129. Nbval = 8; % # b-value
  130. Nbvec = 64; % # gradient direction
  132. % direction dependent diffusion signals
  133. S_dir = zeros(Nr, Nbval, Nbvec);
  134. for i = 1:Nr
  135. rms = simul3DcoCyl_cuda_pgse_bvec(fullfile(root,projname,sprintf('CoCyl_%04u',i)));
  136. sigi = rms.sig;
  137. bval = unique(rms.bval);
  138. sigi = reshape(sigi,Nbvec,[]);
  139. S_dir(i,:,:) = sigi.';
  140. end
  142. % effective radius of the n-th order, um
  143. % n = 0, <r> mean radius
  144. % n = 1, <r^2>/<r> volume weighted averaged radius
  145. % n = 2, (<r^3>/<r>)^(1/2) effective radius for narrow-pulse
  146. % n = 4, (<r^5>/<r>)^(1/4) effective radius for wide-pulse
  147. n = [0 1 2 4];
  148. r_mom = zeros(Nr, numel(n));
  149. vol = zeros(Nr, 1); % ~volume, um2
  150. Nm = zeros(Nr, 1); % # myelin layer
  151. for i = 1:Nr
  152. rms = simul3DcoCyl_cuda_pgse_bvec(fullfile(root,projname,sprintf('CoCyl_%04u',i)));
  154. % cylindrical shell radius at the middle thickness, um
  155. r = rms.rCir + (1:rms.Nm - 1/2)*rms.lm;
  157. % inner radius, um
  158. ri = rms.rCir;
  160. % outer radius, um
  161. ro = rms.rCir + rms.Nm*rms.lm;
  163. vol(i) = pi*(ro^2 - ri^2);
  164. Nm(i) = rms.Nm;
  165. for k = 1:numel(n)
  166. ni = n(k);
  167. if ni==0
  168. r_mom(i,k) = sum(r);
  169. else
  170. r_mom(i,k) = sum(r.^(ni+1));
  171. end
  172. end
  173. end
  175. %% MyeCaliber
  176. % Gamma distribution for axon radius
  177. ri_mean = (0.25:0.05:1).'; % mean, um
  178. ri_var = (ri_mean/2).^2; % variance, um2
  179. b = ri_var./ri_mean; % scale parameter, um
  180. a = ri_mean./b; % shape parameter
  181. x = (0.1:0.1:5).'; % sampled axon radius, um
  183. Nbval = 8; % # b-value
  184. Nbvec = 64; % # gradient direction per b-shell
  185. SNR = [Inf 50 20 10]; % signal-to-noise ratio for Rician noise
  187. % fitted parameters: radius, axial diffusivity, their standard deviations
  188. r_fit = zeros(numel(ri_mean), numel(SNR), 2);
  189. D_fit = zeros(numel(ri_mean), numel(SNR), 2);
  190. r_std = zeros(numel(ri_mean), numel(SNR), 2);
  191. D_std = zeros(numel(ri_mean), numel(SNR), 2);
  193. % spherical mean signal with noise for the fitting
  194. S_fit = zeros(numel(ri_mean), numel(SNR), Nbval);
  196. % model selection
  197. % narrow: narrow-pulse solution (Canales-Rodriguez et al. 2005)
  198. % wide: wide-pulse solution
  199. pulsetype = {'narrow', 'wide'};
  201. tic;
  202. for j = 1:numel(ri_mean)
  203. % direction-dependent diffusion signals from axons with radii in Gamma
  204. % distribution
  205. Ni = gampdf(x, a(j), b(j));
  206. Si = sum(Ni.*vol.*S_dir)/sum(Ni.*vol);
  207. Si = squeeze(Si);
  208. for i = 1:numel(SNR)
  209. % apply Rician noise to direction-dependent diffusion signals
  210. sigma = 1/SNR(i);
  211. Sj = abs( Si + sigma*randn(Nbval, Nbvec) + 1j*sigma*randn(Nbval, Nbvec) );
  213. % calculate spherical mean signal
  214. Sj = mean(Sj, 2);
  216. % Rician noise floor correction
  217. Sj = sqrt(max(Sj.^2-sigma^2, 0));
  219. % model fitting
  220. S_fit(j,i,:) = Sj;
  221. for k = 1:numel(pulsetype)
  222. [r_fit(j,i,k), D_fit(j,i,k), r_std(j,i,k), D_std(j,i,k)] = myelinADMfit(Sj, bval, rms.del, rms.Del, pulsetype{k});
  223. end
  224. end
  225. end
  226. toc;
  228. % effective radius of the n-th order with axon radius distribution, um
  229. % n = 0, <r> mean radius
  230. % n = 1, <r^2>/<r> volume weighted averaged radius
  231. % n = 2, (<r^3>/<r>)^(1/2) effective radius for narrow-pulse
  232. % n = 4, (<r^5>/<r>)^(1/4) effective radius for wide-pulse
  233. n = [0 1 2 4];
  234. r_eff = zeros(numel(ri_mean), numel(n));
  235. for j = 1:numel(ri_mean)
  236. Ni = gampdf(x, a(j), b(j));
  237. for i = 1:numel(n)
  238. ni = n(i);
  239. if ni==0
  240. r_eff(j,i) = sum(Ni.*r_mom(:,i))/sum(Ni.*Nm);
  241. else
  242. I = find(n==0);
  243. r_eff(j,i) = ( sum(Ni.*r_mom(:,i))/sum(Ni.*r_mom(:,I)) ).^(1/ni);
  244. end
  245. end
  246. end
  248. %% Canales-Rodriguez et al. 2025, full solution fitting
  250. % build training data for random forest regression
  251. rCR = 0.01:0.01:3; % axon radius, um
  252. DCR = 0.01:0.01:1.6; % intrinsic diffusivity, um2/ms
  253. % normalized spherical mean signal
  254. SCR = zeros(numel(bval),numel(rCR),numel(DCR));
  255. option = 'fast';
  256. tic;
  257. if strcmpi(option, 'slow')
  258. for i = 1:numel(bval)
  259. SCR(i,:,:) = myelinADMmodelCR(rCR(:), bval(i), rms.del, rms.Del, DCR(:).','slow');
  260. end
  261. else
  262. for i = 1:numel(bval)
  263. for j = 1:numel(DCR)
  264. SCR(i,:,j) = myelinADMmodelCR(rCR(:), bval(i), rms.del, rms.Del, DCR(j),'fast');
  265. end
  266. end
  267. end
  268. toc;
  270. % random forest regression for each SNR
  271. X = permute(SCR, [2, 3, 1]);
  272. X = reshape(X, [], 8);
  273. Y = repmat(rCR(:),numel(DCR),1);
  274. r_fit_CR = zeros(numel(ri_mean), numel(SNR));
  275. tic;
  276. for i = 1:numel(SNR)
  277. % apply Rician noise
  278. sigma = 1/SNR(i);
  279. Xi = abs( X + sigma*randn([size(X), Nbvec]) + 1j*sigma*randn([size(X), Nbvec]) );
  280. Xi = mean(Xi, 3);
  282. % Rician noise floor correction
  283. Xi = sqrt(max(Xi.^2-sigma^2, 0));
  285. % model fitting
  286. Mdl = TreeBagger(10,Xi,Y,'Method','regression','OOBPrediction','On');
  287. for j = 1:numel(ri_mean)
  288. r_fit_CR(j,i) = predict(Mdl,squeeze(S_fit(j,i,:)).');
  289. end
  290. end
  291. toc;
  293. %% Plot normalized spherical mean signal
  294. figure('unit','inch','position',[0 0 18 4]);
  295. cmap = colormap('lines');
  296. mk = {'o','v','s'};
  297. for i = 1:numel(SNR)
  298. subplot(1,numel(SNR),i)
  299. clear h lgtxt
  300. hold on;
  301. list = [1 numel(ri_mean)];
  302. for j = 1:numel(list)
  303. Si = S_fit(list(j),i,:);
  304. Si = squeeze(Si);
  305. h(j) = plot(1./sqrt(bval), Si, mk{j}, 'markersize', 6, 'linewidth', 1, 'color', cmap(j,:));
  306. bval_plot = 1./linspace(0.01, 5, 1000).^2;
  307. S_WP = myelinADMmodel(r_fit(list(j),i,2), bval_plot, rms.del, rms.Del, D_fit(list(j),i,2), 'wide');
  308. plot(1./sqrt(bval_plot), S_WP, '-' , 'linewidth', 1, 'color', cmap(j,:));
  309. lgtxt{j} = sprintf('$2\\bar{r}_i=$%.1f $\\mu$m', ri_mean(list(j))*2);
  310. end
  311. h(3) = plot(-1, -1, 'k-', 'linewidth', 1);
  312. lgtxt{3} = 'fitting';
  313. xlim([0 2]);
  314. ylim([0 1]);
  315. pbaspect([1 1 1]);
  316. xlabel('$1/\sqrt{b}$, $\mu$m$\cdot$ms$^{-1/2}$', 'interpreter', 'latex', 'fontsize', 20)
  317. if i==1, ylabel('$\bar{S}(b)$', 'interpreter', 'latex', 'fontsize', 20); end
  318. title(sprintf('SNR=%u',SNR(i)),'interpreter','latex','fontsize',20)
  319. box on; grid on;
  320. legend(h, lgtxt,'interpreter','latex','fontsize',16,'box','off','location','southeast')
  321. end
  323. %% Resolution limit
  324. D0 = rms.Din; % intrinsic diffusivity, um2/ms
  325. Da = rms.Din; % axial diffusivity, um2/ms
  326. za = 1.64; % z-score at alpha = 0.05
  327. Nav = 64; % # gradient direction per b-shell
  328. r_min_NP = zeros(numel(SNR), 1)-1; % resolution limit of narrow-pulse
  329. r_min_WP = zeros(numel(SNR), 1)-1; % resolution limit of wide-pulse
  330. for j = 1:numel(SNR)
  331. r_min_NP(j) = myelinADMrmin(max(bval), rms.del, rms.Del, D0, Da, za, SNR(j), Nav, 'narrow');
  332. r_min_WP(j) = myelinADMrmin(max(bval), rms.del, rms.Del, D0, Da, za, SNR(j), Nav, 'wide');
  333. end
  335. %% Plot figure, Canales-Rodriguez et al. 2005, narrow pulse solution
  336. figure('unit','inch','position',[0 0 18 4]);
  337. cmap = colormap('lines');
  338. mk = {'o','v','s','^','d'};
  339. clear h hx lgtxt
  340. for i = 1:numel(n)
  341. subplot(1,4,i)
  342. hold on;
  343. for j = numel(SNR):-1:1
  344. plotstd(2*r_eff(:,i), 2*r_fit(:,j,1), 2*r_std(:,j,1), cmap(j,:), 0.3, 'area');
  345. end
  346. end
  348. for i = 1:numel(n)
  349. subplot(1,4,i)
  350. hold on;
  351. for j = 1:numel(SNR)
  352. if j==3
  353. h(j) = plot(r_eff(:,i)*2, r_fit(:,j,1)*2, mk{j}, 'color', cmap(j,:)*0.85,'linewidth',1);
  354. else
  355. h(j) = plot(r_eff(:,i)*2, r_fit(:,j,1)*2, mk{j}, 'color', cmap(j,:),'linewidth',1);
  356. end
  357. if ~isinf(SNR(j))
  358. if j==3
  359. hx = xline(r_min_NP(j)*2); set(hx, 'color', cmap(j,:)*0.85, 'linestyle', '--','linewidth',1);
  360. else
  361. hx = xline(r_min_NP(j)*2); set(hx, 'color', cmap(j,:), 'linestyle', '--','linewidth',1);
  362. end
  363. end
  364. lgtxt{j} = sprintf('SNR=%u', SNR(j));
  365. end
  366. xlim([0 4]);
  367. ylim([0 4]);
  368. yticks(0:4);
  369. hr = refline(1); set(hr, 'color', 'k');
  370. hl = plot(-1, -1, 'k--', 'linewidth', 1); lgtxt{numel(SNR)+1} = '2$r_{\rm min,NP}$';
  371. hg = legend([h hl], lgtxt, 'interpreter', 'latex', 'location', 'northwest', 'box', 'on', 'fontsize', 12);
  372. box on;
  373. pbaspect([1 1 1])
  374. switch n(i)
  375. case 0
  376. xlabel('2$\langle r\rangle$, $\mu$m', 'interpreter', 'latex', 'fontsize',20);
  377. case 1
  378. xlabel('2$r_{\rm vwa}$, $\mu$m', 'interpreter', 'latex', 'fontsize',20);
  379. case 2
  380. xlabel('2$r_{\rm eff, NP}$, $\mu$m', 'interpreter', 'latex', 'fontsize',20);
  381. case 4
  382. xlabel('2$r_{\rm eff, WP}$, $\mu$m', 'interpreter', 'latex', 'fontsize',20);
  383. end
  385. if i==1
  386. ylabel('fitted 2$r$, $\mu$m', 'interpreter', 'latex', 'fontsize',20);
  387. end
  388. end
  390. %% Plot figure, wide-pulse solution
  391. figure('unit','inch','position',[0 0 18 4]);
  392. cmap = colormap('lines');
  393. mk = {'o','v','s','^','d'};
  394. clear h hx lgtxt
  395. for i = 1:numel(n)
  396. subplot(1,4,i)
  397. hold on;
  398. for j = numel(SNR):-1:1
  399. plotstd(2*r_eff(:,i), 2*r_fit(:,j,2), 2*r_std(:,j,2), cmap(j,:), 0.3, 'area');
  400. end
  401. end
  403. for i = 1:numel(n)
  404. subplot(1,4,i)
  405. hold on;
  406. for j = 1:numel(SNR)
  407. if j==3
  408. h(j) = plot(2*r_eff(:,i), 2*r_fit(:,j,2), mk{j}, 'color', cmap(j,:)*0.85,'linewidth',1);
  409. else
  410. h(j) = plot(2*r_eff(:,i), 2*r_fit(:,j,2), mk{j}, 'color', cmap(j,:),'linewidth',1);
  411. end
  412. if ~isinf(SNR(j))
  413. if j==3
  414. hx = xline(2*r_min_WP(j)); set(hx, 'color', cmap(j,:)*0.85, 'linestyle', '--','linewidth',1);
  415. else
  416. hx = xline(2*r_min_WP(j)); set(hx, 'color', cmap(j,:), 'linestyle', '--','linewidth',1);
  417. end
  418. end
  419. lgtxt{j} = sprintf('SNR=%u', SNR(j));
  420. end
  421. xlim([0 4]);
  422. ylim([0 4]);
  423. yticks(0:4);
  424. hr = refline(1); set(hr, 'color', 'k');
  425. hl = plot(-1, -1, 'k--', 'linewidth', 1); lgtxt{numel(SNR)+1} = '2$r_{\rm min,WP}$';
  426. if i==1
  427. hg = legend([h hl], lgtxt, 'interpreter', 'latex', 'location', 'southeast', 'box', 'on', 'fontsize', 12);
  428. else
  429. hg = legend([h hl], lgtxt, 'interpreter', 'latex', 'location', 'northwest', 'box', 'on', 'fontsize', 12);
  430. end
  431. box on;
  432. pbaspect([1 1 1])
  433. switch n(i)
  434. case 0
  435. xlabel('2$\langle r\rangle$, $\mu$m', 'interpreter', 'latex', 'fontsize',20);
  436. case 1
  437. xlabel('2$r_{\rm vwa}$, $\mu$m', 'interpreter', 'latex', 'fontsize',20);
  438. case 2
  439. xlabel('2$r_{\rm eff, NP}$, $\mu$m', 'interpreter', 'latex', 'fontsize',20);
  440. case 4
  441. xlabel('2$r_{\rm eff, WP}$, $\mu$m', 'interpreter', 'latex', 'fontsize',20);
  442. end
  443. if i==1
  444. ylabel('fitted 2$r$, $\mu$m', 'interpreter', 'latex', 'fontsize',20);
  445. end
  446. end
  448. %% Plot figure, Canales Rodriguez et al. 2005, narrow-pulse full solution
  449. figure('unit','inch','position',[0 0 18 4]);
  450. cmap = colormap('lines');
  451. mk = {'o','v','s','^','d'};
  452. clear h hx lgtxt
  453. for i = 1:numel(n)
  454. subplot(1,4,i)
  455. hold on;
  456. for j = 1:numel(SNR)
  457. h(j) = plot(2*r_eff(:,i), 2*r_fit_CR(:,j), mk{j}, 'color', cmap(j,:),'linewidth',1);
  458. if ~isinf(SNR(j))
  459. if j==3
  460. hx = xline(2*r_min_NP(j)); set(hx, 'color', cmap(j,:)*0.85, 'linestyle', '--','linewidth',1);
  461. else
  462. hx = xline(2*r_min_NP(j)); set(hx, 'color', cmap(j,:), 'linestyle', '--','linewidth',1);
  463. end
  464. end
  465. lgtxt{j} = sprintf('SNR=%u', SNR(j));
  466. end
  467. xlim([0 4]);
  468. ylim([0 4]);
  469. hr = refline(1); set(hr, 'color', 'k');
  470. hl = plot(-1, -1, 'k--', 'linewidth', 1); lgtxt{numel(SNR)+1} = '2$r_{\rm min,NP}$';
  471. if i==1
  472. hg = legend([h hl], lgtxt, 'interpreter', 'latex', 'location', 'northwest', 'box', 'on', 'fontsize', 12);
  473. else
  474. hg = legend([h hl], lgtxt, 'interpreter', 'latex', 'location', 'northwest', 'box', 'on', 'fontsize', 12);
  475. end
  476. box on;
  477. pbaspect([1 1 1])
  478. switch n(i)
  479. case 0
  480. xlabel('2$\langle r\rangle$, $\mu$m', 'interpreter', 'latex', 'fontsize',20);
  481. case 1
  482. xlabel('2$r_{\rm vwa}$, $\mu$m', 'interpreter', 'latex', 'fontsize',20);
  483. case 2
  484. xlabel('2$r_{\rm eff, NP}$, $\mu$m', 'interpreter', 'latex', 'fontsize',20);
  485. case 4
  486. xlabel('2$r_{\rm eff, WP}$, $\mu$m', 'interpreter', 'latex', 'fontsize',20);
  487. end
  488. if i==1
  489. ylabel('fitted 2$r$, $\mu$m', 'interpreter', 'latex', 'fontsize',20);
  490. end
  491. end
  493. %% Plot histology information
  495. % inner radius, um
  496. r = 0.1:0.1:5;
  497. r = r(:);
  499. % thickness of myelin layer, um
  500. lm = 12/1e3;
  502. % # myelin layers
  503. C0 = 0.35/lm;
  504. C1 = 0.006/lm;
  505. C2 = 0.024/lm;
  506. Nm = round(C0 + C1*2*r + C2*log(2*r));
  508. % Total myelin thickness, um
  509. Lm = Nm*lm;
  511. % g-ratio: ratio of inner to outer radii
  512. gratio = r./(r+Lm);
  514. figure('unit','inch','position',[0 0 15 5])
  516. % plot inner diameter vs g-ratio
  517. subplot(131)
  518. hold on;
  519. hp = plot(2*r, gratio, 'o', 'linewidth', 1);
  521. xx = linspace(0, 10, 100);
  522. yy = C0 + C1*xx + C2*log(xx);
  523. yy = xx./(xx + 2*yy*lm);
  524. hh = plot(xx, yy, '-', 'linewidth', 1);
  526. legend([hh, hp], {'log-linear', 'numerical phantom'}, 'interpreter', 'latex', ...
  527. 'box', 'off', 'location', 'southeast', 'fontsize', 14)
  529. pbaspect([1 1 1]);
  530. xlabel('inner diameter 2$r_i$, $\mu$m', 'interpreter', 'latex', 'fontsize', 20);
  531. ylabel('$g$-ratio', 'interpreter', 'latex', 'fontsize', 20)
  532. xlim([0 10]);
  533. ylim([0 1]);
  534. yticks(0:0.2:1);
  535. box on;
  537. % plot inner diameter vs # myelin layers
  538. subplot(132)
  539. hold on;
  540. hp = plot(2*r, Nm, 'o', 'linewidth', 1);
  542. xx = linspace(0, 10, 100);
  543. yy = C0 + C1*xx + C2*log(xx);
  544. hh =plot(xx, yy, '-', 'linewidth', 1);
  546. legend([hh, hp], {'log-linear', 'numerical phantom'}, 'interpreter', 'latex', ...
  547. 'box', 'off', 'location', 'southeast', 'fontsize', 14)
  549. pbaspect([1 1 1]);
  550. xlabel('inner diameter 2$r_i$, $\mu$m', 'interpreter', 'latex', 'fontsize', 20);
  551. ylabel('\# myelin layers $N_m$', 'interpreter', 'latex', 'fontsize', 20)
  552. xlim([0 10]); ylim([25 40])
  553. box on;
  555. % plot axon diameter distribution in Gamma distribution
  556. di_mean = 2*(0.25:0.05:1).';
  557. di_var = (di_mean/2).^2;
  558. b = di_var./di_mean;
  559. a = di_mean./b;
  560. di = 2*(0.1:0.1:5); di = di(:);
  561. Ni = zeros(numel(di), numel(di_mean));
  562. for j = 1:numel(di_mean)
  563. Ni(:,j) = gampdf(di, a(j), b(j));
  564. end
  566. ax = subplot(133);
  567. hold on;
  568. plot([2.5 2.5], [0.9 1.9], 'k-', 'linewidth', 1.5)
  569. text(1.5, 1.4, '2$\bar{r}_i$', 'interpreter', 'latex', 'fontsize', 14);
  570. cmap = colormap('parula');
  571. list = 1:numel(di_mean);
  572. Pi = Ni./sum(Ni,1)/mean(diff(di));
  573. clear h lgtxt
  574. for i = 1:numel(list)
  575. h(i) = plot(di, Pi(:,list(i)), '.-', 'color', cmap(i*floor(size(cmap,1)/numel(list)),:));
  576. end
  577. pbaspect([1 1 1]);
  578. box on;
  579. xlabel('inner diameter 2$r_i$, $\mu$m', 'interpreter', 'latex', 'fontsize', 20);
  580. ylabel('PDF, $\mu$m$^{-1}$', 'interpreter', 'latex', 'fontsize', 20)
  581. for i = 1:numel(list)
  582. lgtxt{i} = sprintf('%.2f $\\mu$m', di_mean(list(i)));
  583. end
  584. legend(h, lgtxt, 'interpreter','latex','fontsize',12, 'box', 'off', ...
  585. 'NumColumns',2);
  586. xlim([0 10]);
  587. ylim([0 2]);
  588. yticks(0:5)

demo_3_cylinder_shell_ADM.m at commit 9de680c, under MIT · at the source

Overview

Authors: Hong-Hsi Lee1, Kwok-Shing Chan1, Dmitry S. Novikov2, Els Fieremans2, Susie Y. Huang1
ORCID iDs: Hong-Hsi Lee, Kwok-Shing Chan, Susie Y. Huang
  1. Department of Radiology, Athinoula A. Martinos Center for Biomedical Imaging, Massachusetts General Hospital, Charlestown, MA 02129 USA, and also with the Harvard Medical School, Boston, MA 02115 USA
  2. Center for Advanced Imaging Innovation and Research, New York, NY 10016 USA, and also with the Department of Radiology, New York University School of Medicine, New York, NY 10016 USA
Institutions: Massachusetts General Hospital (United States); Athinoula A. Martinos Center for Biomedical Imaging (United States); Center for Advanced Imaging Innovation and Research (United States)
Journal: IEEE transactions on medical imaging, volume 45, issue 6, pages 2883-2896
Dates: published online 1 June 2026; in print June 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1109/tmi.2026.3664328 · PMID 41686667 · PMCID PMC13387500 · OpenAlex W7128810266
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: structural MRI / diffusion (modality), human (organism), cellular / molecular (subfield)
Keywords: Diffusion MRI, axon diameter, myelinated axon, Monte Carlo simulations, microstructure imaging
MeSH: Axons*, Diffusion Magnetic Resonance Imaging*, Image Processing, Computer-Assisted*, Myelin Sheath*, Animals, Computer Simulation, Humans, Monte Carlo Method, Phantoms, Imaging, Water (* major topic)
Topic: Advanced Neuroimaging Techniques and Applications (Radiology, Nuclear Medicine and Imaging, Medicine), according to OpenAlex
Funding: NIBIB NIH HHS (P41 EB030006, U01 EB026996, P41 EB015896, P41 EB017183); National Institute of Health (P41EB015896, R01NS088040, R01NS118187, R21NS081230, P41EB030006, P41EB017183, DP5OD031854, U01EB026996, U24NS137077); NINDS NIH HHS (R01 NS088040, R01 NS118187, U24 NS137077, R21 NS081230); ZonMw (04520232330012); NIH HHS (DP5 OD031854); NWO/ZonMw through the Rubicon (04520232330012)
Citations: not cited yet (Europe PMC); 101 references in the paper

Abstract

Probing diffusion in myelin water using diffusion-weighted T1-/T2-selective MRI acquisitions enables noninvasive measurement of myelinated axon diameter. Its application for in vivo measurements requires numerical verification through diffusion simulations. Here, we propose the theory of myelin water diffusion as measured with a diffusion MRI pulse sequence with wide gradient pulses using the Gaussian phase approximation. We establish its applicability to axonal diameter mapping via Monte Carlo simulations in either infinitely thin cylindrical surfaces or concentric cylindrical shells of finite thickness, mimicking the micro-geometry of myelin sheaths. The estimated diameters are shown to be weighted more toward outer than inner calibers. Simulation results evaluate the theory of myelin water diffusion and axon diameter estimation using spherical mean diffusion signals, demonstrating its applicability at signal-to-noise ratio above 20 on the Connectome 2.0 MRI scanner equipped with maximum gradient strength of 500 mT/m and slew rate of 600 T/m/s. Measuring restricted diffusion of myelin water in-between myelin sheaths using diffusion MRI allows one to measure myelinated axon diameters in vivo. The protocol can potentially be adapted for clinically available high-gradient performance scanners.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 16 matches between paragraphs and lines of code.

Connectome20/MyCaliber

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 9de680c81bfa1180ab3af722c203496535e405c6, 25 February 2026
Languages: MATLAB (20), CUDA (3)
Size: 30 files, 23 scripts
Software Heritage: not archived
Found in: the text, “Diffusion Simulations in Cylindrical Shells With”
Holds: README, license file
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Statistics and Machine Learning Toolbox (5 files), Optimization Toolbox (3 files), Image Processing Toolbox (1 file), MRtrix3 (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
25 files
At the source: github.com/Connectome20/MyCaliber

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 23 scripts, each with its path and the digest of its content;
  • 16 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 2, 28 September 2026

  • Publisher: n/a → Institute of Electrical and Electronics Engineers

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 5 authors, 5 keywords, 10 MeSH terms, 6 funders, 91 references.

Cite

This paper

Lee, H.-H., Chan, K.-S., Novikov, D. S., Fieremans, E., & Huang, S. Y. (2026). Axon Diameter Mapping From Myelin Water Diffusion MRI. IEEE transactions on medical imaging, 45(6), 2883-2896. https://doi.org/10.1109/tmi.2026.3664328

BibTeX

@article{lee2026axon,
author = {Lee, Hong-Hsi and Chan, Kwok-Shing and Novikov, Dmitry S. and Fieremans, Els and Huang, Susie Y.},
title = {{Axon Diameter Mapping From Myelin Water Diffusion MRI}},
journal = {IEEE transactions on medical imaging},
year = {2026},
month = jun,
volume = {45},
number = {6},
pages = {2883--2896},
publisher = {Institute of Electrical and Electronics Engineers},
issn = {0278-0062},
doi = {10.1109/tmi.2026.3664328},
url = {https://doi.org/10.1109/tmi.2026.3664328},
pmid = {41686667},
pmcid = {PMC13387500}
}

RIS

TY - JOUR
AU - Lee, Hong-Hsi
AU - Chan, Kwok-Shing
AU - Novikov, Dmitry S.
AU - Fieremans, Els
AU - Huang, Susie Y.
TI - Axon Diameter Mapping From Myelin Water Diffusion MRI
T2 - IEEE transactions on medical imaging
J2 - IEEE Trans Med Imaging
PY - 2026
DA - 2026/06/01
VL - 45
IS - 6
SP - 2883
EP - 2896
SN - 0278-0062
PB - Institute of Electrical and Electronics Engineers
DO - 10.1109/tmi.2026.3664328
UR - https://doi.org/10.1109/tmi.2026.3664328
LA - en
ER -

CSL-JSON

{
"id": "10.1109/tmi.2026.3664328",
"type": "article-journal",
"title": "Axon Diameter Mapping From Myelin Water Diffusion MRI",
"container-title": "IEEE transactions on medical imaging",
"author": [
{
"family": "Lee",
"given": "Hong-Hsi"
},
{
"family": "Chan",
"given": "Kwok-Shing"
},
{
"family": "Novikov",
"given": "Dmitry S."
},
{
"family": "Fieremans",
"given": "Els"
},
{
"family": "Huang",
"given": "Susie Y."
}
],
"container-title-short": "IEEE Trans Med Imaging",
"volume": "45",
"issue": "6",
"page": "2883-2896",
"DOI": "10.1109/tmi.2026.3664328",
"PMID": "41686667",
"PMCID": "PMC13387500",
"ISSN": "0278-0062",
"publisher": "Institute of Electrical and Electronics Engineers",
"URL": "https://doi.org/10.1109/tmi.2026.3664328",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
1
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1002/hbm.70553 [code]
Axon Diameter Mapping in the Living Human Brain with Ultra-High-Gradient Diffusion MRI at 500 mT/m Gradient Strength.
Journal: Human brain mapping
In common: MRtrix3, Optimization Toolbox, Image Processing Toolbox, 1 other tool, structural MRI / diffusion, 23 references, 2 authors
[2] doi:10.1002/mrm.70490 [code]
Dependence of the Extra-Cellular Diffusion Coefficient on the Fractions of Neurites and Cell Bodies in Gray Matter.
Journal: Magnetic resonance in medicine
In common: structural MRI / diffusion, cellular / molecular, 22 references, author Hong‐Hsi Lee
[3] doi:10.1371/journal.pbio.3003861
Learning engages transient and sustained cellular mechanisms in the human brain.
Journal: PLoS biology
In common: structural MRI / diffusion, cellular / molecular, 10 references
[4] doi:10.1002/mrm.70336 [code]
Offline Reconstruction of Diffusion MRI Acquisitions for Comparison Between Complex PCA-Based and AI-Based Denoising.
Journal: Magnetic resonance in medicine
In common: MRtrix3, Optimization Toolbox, Image Processing Toolbox, 1 other tool, structural MRI / diffusion, 5 references
[5] doi:10.1038/s41598-026-39162-7 [code]
White matter microstructure differences in obstructive sleep apnea severity groups assessed by diffusion tensor metrics and biophysical modeling.
Journal: Scientific reports
In common: MRtrix3, structural MRI / diffusion, 7 references
[6] doi:10.1002/mrm.70378 [code]
Investigating the Sensitivity of the Diffusion MRI Signal to Magnetization Transfer and Permeability via Monte-Carlo Simulations.
Journal: Magnetic resonance in medicine
In common: structural MRI / diffusion, 8 references
[7] doi:10.1038/s41467-026-70018-w [code]
Geometry of the cumulant series in diffusion MRI.
Journal: Nature communications
In common: Statistics and Machine Learning Toolbox, structural MRI / diffusion, 7 references
[8] doi:10.1038/s41598-026-51531-w [code]
Multimodal age-dependent diffusion-MRI analysis of the neocortex in a rat model of cortical dysplasia.
Journal: Scientific reports
In common: MRtrix3, Image Processing Toolbox, Statistics and Machine Learning Toolbox, structural MRI / diffusion, 4 references
[9] doi: [code]
Diffusion-relaxation MRI as virtual histology: separable microstructural signatures of AD pathology in ex vivo human brain
Journal: Research square
In common: Image Processing Toolbox, Statistics and Machine Learning Toolbox, structural MRI / diffusion, cellular / molecular, 3 references
[10] doi:10.7554/elife.107661 [code]
In vivo mapping of striatal neurodegeneration in Huntington's disease with Soma and Neurite Density Imaging.
Journal: eLife
In common: Optimization Toolbox, Image Processing Toolbox, Statistics and Machine Learning Toolbox, structural MRI / diffusion, 3 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.