OSCR

Multiscale parcellation of dynamic causal models of the brain.

Code ↔ Paper

9 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 9 matches · 2 of them tie a paragraph to a whole file, not to given lines: weak matches, whose lines are not tinted
  1. [1] § Results › Empirical network › Conditional independence and Markov blanket ↔ Code/DEMO2_MB.m, lines 37–58 · score 0.76 · stretched beta prior, Bayes factor, posterior distribution, partial correlation, hyperparameter, analytic
  2. [2] § Materials and Methods › Empirical network ↔ Code/spm_dcm_J_mod.m, lines 1–72 · score 0.68 · Bayesian model reduction, local connectivity, redundant, linearized, effective connectivity, voxels
  3. [3] § Results › Empirical network › Conditional independence and Markov blanket ↔ Code/DEMO2_MB.m, lines 37–58 · score 0.64 · stretched beta prior, Bayes factors, partial correlation, hypothesis, log
  4. [4] § Results › Empirical network › Conditional independence and Markov blanket ↔ Code/spm_mb_ui_mod.m, lines 1–60 · score 0.59 · Markov blanket states, neuronal states, internal states, Sparsity, MB, partitioned
  5. [5] § Results › Empirical network › Conditional independence and Markov blanket ↔ Code/spm_mb_ui_mod.m, lines 1–60 · score 0.58 · Identifying Markov blanket, internal states, blanket states, blocks, couplings, causal
  6. [6] § Results › Empirical network › Conditional independence and Markov blanket ↔ Code/plot_HJ.m, the whole file · a weak match · score 0.57 · Cross PT, MB inducing, B1, B2, PT1, INT1
  7. [7] § Materials and Methods › Empirical network ↔ Code/spm_mb_ui_mod.m, lines 324–374 · score 0.57 · maximum intensity projection, variance explained, eigenmode, coupling, matrix, partitioned
  8. [8] § Materials and Methods › Empirical network ↔ Code/compute_J.m, the whole file · a weak match · score 0.57 · fMRI, visual motion, effective connectivity, variance, matrix
  9. [9] § Materials and Methods › Simulated networks ↔ Code/spm_dcm_J_mod.m, lines 1–72 · score 0.53 · Bayesian model reduction, log evidence, bmr, posterior, causal, spm

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 · 1,268 lines · 57 KB · AGPL-3.0 · 3 matches

  1. function [MB] = spm_mb_ui_mod(action,varargin)
  2. % This is a minimally MODIFIED version of spm_mb_ui.m in SPM12
  3. % Modified by Tahereh Zarghami, 2025
  5. % VOI extraction of adjusted data and Markov Blanket decomposition
  6. % FORMAT [MB] = spm_mb_ui('specify',SPM)
  7. % FORMAT [MB] = spm_mb_ui('blocking',MB)
  8. % FORMAT [MB] = spm_mb_ui('results' ,MB)
  9. %
  10. % SPM - structure containing generic analysis details
  11. %
  12. % MB.contrast - contrast name
  13. % MB.name - MB name
  14. % MB.c - contrast weights
  15. % MB.X - contrast subspace
  16. % MB.Y - whitened and adjusted data
  17. % MB.X0 - null space of contrast
  18. %
  19. % MB.XYZ - locations of voxels (mm)
  20. % MB.xyz - seed voxel location (mm)
  21. % MB.VOX - dimension of voxels (mm)
  22. %
  23. % MB.V - canonical vectors (data)
  24. % MB.v - canonical variates (data)
  25. % MB.W - canonical vectors (design)
  26. % MB.w - canonical variates (design)
  27. % MB.C - canonical contrast (design)
  28. %
  29. % MB.chi - Chi-squared statistics testing D >= i
  30. % MB.df - d.f.
  31. % MB.p - p-values
  32. %
  33. % also saved in MB_*.mat in the SPM working directory
  34. %
  35. % FORMAT [MB] = spm_cva_ui('results',MB)
  36. % Display the results of a MB analysis
  37. %__________________________________________________________________________
  38. %
  39. % This routine uses the notion of Markov blankets and the renormalisation
  40. % group to evaluate the coupling among neuronal systems at increasing
  41. % spatial scales. The underlying generative model is based upon the
  42. % renormalisation group: a working definition of renormalization involves
  43. % three elements: vectors of random variables, a course-graining operation
  44. % and a requirement that the operation does not change the functional form
  45. % of the Lagrangian. In our case, the random variables are neuronal states;
  46. % the course graining operation corresponds to the grouping (G) into a
  47. % particular partition and adiabatic reduction (R) - that leaves the
  48. % functional form of the dynamics unchanged.
  49. %
  50. % Here, the grouping operator (G) is based upon coupling among states as
  51. % measured by the Jacobian. In brief, the sparsity structure of the
  52. % Jacobian is used to recursively identify Markov blankets around internal
  53. % states to create a partition of states - at any level - into particles;
  54. % where each particle comprises external and blanket states. The ensuing
  55. % reduction operator (R) eliminates the internal states and retains the slow
  56. % eigenmodes of the blanket states. These then constitute the (vector)
  57. % states at the next level and the process begins again.
  58. %
  59. % This routine starts using a simple form of dynamic causal modelling
  60. % applied to the principal eigenvariate of local parcels (i.e., particles)
  61. % of voxels with compact support. The Jacobian is estimated using a
  62. % linearised dynamic causal (state space) model, where observations are
  63. % generated by applying a (e.g., haemodynamic) convolution operator to
  64. % hidden (e.g., neuronal) states. This estimation uses parametric empirical
  65. % Bayes (PEB: spm_PEB). The ensuing estimates of the Jacobian (i.e.,
  66. % effective connectivity) are reduced using Bayesian model reduction (BMR:
  67. % spm_dcm_BMR_all) within a bespoke routine (spm_dcm_J).
  68. %
  69. % The Jacobian is then partitioned using the course graining operator into
  70. % particles or parcels (using spm_markov blanket). The resulting partition
  71. % is then reduced by eliminating internal states and retaining slow
  72. % eigenmodes with the largest (real) eigenvalues (spm_A_reduce). The
  73. % Jacobian of the reduced states is then used to repeat the process -
  74. % recording the locations of recursively coarse-grained particles - until
  75. % there is a single particle.
  76. %
  77. % The result of this recursive decomposition (i.e., renormalisation)
  78. % affords a characterisation of directed coupling, as characterised by a
  79. % complex Jacobian; namely, a multivariate coupling matrix, describing the
  80. % coupling between eigenmodes of Markov blankets at successive scales. This
  81. % can be regarded as a recursive parcellation scheme based upon effective
  82. % connectivity and a generative (dynamic causal) model of multivariate
  83. % (neuronal) timeseries.
  84. %
  85. % The following lists the various results options. please see main body of
  86. % this script for a description of the (graphical) output
  87. %
  88. % display the results in terms of particular partitions and eigenmodes
  89. %--------------------------------------------------------------------------
  90. % spm_mb_ui('results',MB,'anatomy');
  91. %
  92. % characterise connectivity at the smallest scale
  93. %--------------------------------------------------------------------------
  94. % spm_mb_ui('results',MB,'distance');
  95. %
  96. % characterise scaling behaviour in terms of scaling exponent
  97. %--------------------------------------------------------------------------
  98. % spm_mb_ui('results',MB,'scaling');
  99. %
  100. % characterise intrinsic coupling in terms of transfer functions
  101. %--------------------------------------------------------------------------
  102. % spm_mb_ui('results',MB,'kernels');
  103. %
  104. % display the results in terms of particular partitions and eigenmodes
  105. %--------------------------------------------------------------------------
  106. % spm_mb_ui('results',MB,'dynamics');
  107. %
  108. % characterise extrinsic coupling with a connectogram
  109. %--------------------------------------------------------------------------
  110. % spm_mb_ui('results',MB,'connectogram');
  111. %
  112. % characterise extrinsic coupling in terms of cross covariance functions
  113. %--------------------------------------------------------------------------
  114. % spm_mb_ui('results',MB,'connectivity');
  115. %
  116. % characterise intrinsic coupling in terms of dissipative flow
  117. %--------------------------------------------------------------------------
  118. % spm_mb_ui('results',MB,'eigenmodes');
  119. %
  120. % characterise eigenmodes in terms of design or inputs
  121. %--------------------------------------------------------------------------
  122. % spm_mb_ui('results',MB,'responses');
  123. %
  124. % input effects as active states at base level
  125. %--------------------------------------------------------------------------
  126. % spm_mb_ui('results',MB,'inputs');
  127. %__________________________________________________________________________
  128. % Copyright (C) 2008-2014 Wellcome Trust Centre for Neuroimaging
  130. % Karl Friston
  131. % $Id: spm_mb_ui.m 7808 2020-03-31 11:18:26Z karl $
  133. OPT.d = 32; % maximum connection length (mm)
  134. OPT.np = 1024; % number of parcels (particles)
  135. OPT.xyz = [0;0;0]; % centre of (spherical) VOI (mm)
  136. OPT.spec = 128; % radius of (spherical) VOI (mm)
  137. OPT.Ic = 1; % contrast (for adjusted data)
  138. OPT.T = 1; % threshold for adiabatic reduction
  139. OPT.N = 8; % max modes for adiabatic reduction
  141. switch lower(action)
  143. case 'specify'
  144. %==================================================================
  145. % M B : S P E C I F Y
  146. %==================================================================
  147. if isempty(varargin)
  148. [SPM,sts] = spm_select(1,'mat','Select SPM',[],[],'^SPM.*\.mat$');
  149. if ~sts, return; end
  150. else
  151. SPM = varargin{1};
  152. end
  153. if ischar(SPM)
  154. SPM = load(SPM);
  155. SPM = SPM.SPM;
  156. end
  158. %-Contrast specification
  159. %------------------------------------------------------------------
  160. con = SPM.xCon(OPT.Ic).name;
  161. c = SPM.xCon(OPT.Ic).c;
  162. c = full(c);
  164. %-Extract required data from results files
  165. %==================================================================
  167. %-Get (normalised) explanatory variables (data)
  168. %------------------------------------------------------------------
  169. cd(SPM.swd)
  171. XYZ = SPM.xVol.XYZ;
  172. F = spm_get_data(SPM.xCon(1).Vspm,XYZ);
  173. XYZ = XYZ(:,logical(F > 1));
  174. R = spm_get_data(SPM.VResMS,XYZ); % computation order swapped with Y
  175. cd("..\attention\functional") % I added
  176. Y = spm_get_data(SPM.xY.VY, XYZ);
  177. Y = bsxfun(@times,Y,1./sqrt(R));
  178. XYZmm = SPM.xVol.M(1:3,1:3)*XYZ;
  179. XYZmm = bsxfun(@plus,XYZmm,SPM.xVol.M(1:3,4));
  181. %-Remove serial correlations and get design (note X := W*X)
  182. %------------------------------------------------------------------
  183. Y = SPM.xX.W*Y;
  184. X = SPM.xX.xKXs.X;
  185. M = SPM.xVol.M(1:3,1:3); %-voxels to mm
  186. VOX = diag(sqrt(diag(M'*M))'); %-voxel size
  188. %-Null-space
  189. %------------------------------------------------------------------
  190. X0 = [];
  191. try, X0 = [X0 blkdiag(SPM.xX.K.X0(:,1:min(16,end)))]; end %-drift terms
  192. try, X0 = [X0 spm_detrend(SPM.xGX.gSF)]; end %-global estimate
  193. X0 = full(spm_svd([X0, (X - X*c*pinv(c))])); %-null space of c
  194. Y = Y - X0*(pinv(X0)*Y);
  196. % exogenous inputs (decimated)
  197. %------------------------------------------------------------------
  198. for i = 1:numel(SPM.Sess.U)
  199. U(:,i) = SPM.Sess.U(i).u(33:end,:);
  200. name{i} = SPM.Sess.U(i).name{1};
  201. end
  202. Dy = spm_dctmtx(size(Y,1),size(Y,1));
  203. Du = spm_dctmtx(size(U,1),size(Y,1));
  204. Dy = Dy*sqrt(size(Y,1)/size(U,1));
  205. U = Dy*(Du'*U);
  207. %-First level partition: compact support
  208. %==================================================================
  209. spm_figure('GetWin','Markov blanket');clf
  210. subplot(2,2,1), spm_mip(std(Y),XYZmm,VOX);
  211. axis image, axis off, title('Variance of voxels','FontSize',16)
  213. % local SVD
  214. %------------------------------------------------------------------
  215. V = spm_mvb_U(Y,'compact',X0,XYZ,[],OPT.np);
  216. Y = Y*V;
  217. nv = size(Y,2);
  218. x = cell(3,nv);
  219. x(2,:) = num2cell(1:nv); % indices of partition
  220. z = num2cell(1:nv); % indices of states
  223. % distance matrix (mm)
  224. %------------------------------------------------------------------
  225. [R,S] = spm_mb_distance(V,XYZmm);
  228. subplot(2,2,2), spm_mip(abs(V)*std(Y)',XYZmm,VOX);
  229. axis image, axis off, title('Variance of modes','FontSize',16)
  230. subplot(2,2,4), spm_mip(logical(V)*std(Y)',XYZmm,VOX);
  231. axis image, axis off, title('Variance of particles','FontSize',16)
  232. subplot(2,2,3), imagesc(R), axis image, axis off
  233. title('Disttance between particles','FontSize',16), drawnow
  235. % self inhibition based on spatial scale
  236. %------------------------------------------------------------------
  237. I = -S^(-3/2)*16;
  239. % first order Jacobian
  240. %==================================================================
  241. [J,K,VAR] = spm_dcm_J_mod(Y,U,X0,SPM.xY.RT,R,I,OPT.d); % VAR added
  243. % Populate first level of market blanket structure
  244. %==================================================================
  245. clear MB
  247. %-Assemble results
  248. %------------------------------------------------------------------
  249. MB{1}.swd = SPM.swd; % results directory
  250. MB{1}.dt = SPM.xY.RT; % sampling interval
  251. MB{1}.con = con; % contrast name
  252. MB{1}.name = name; % input name
  253. MB{1}.XYZ = XYZmm; % locations of voxels (mm)
  254. MB{1}.VOX = VOX; % dimension of voxels (mm)
  256. MB{1}.J = J; % Jacobian (states)
  257. MB{1}.K = K; % Jacobian (inputs)
  258. MB{1}.VAR = VAR; % Variance %%% I added
  259. MB{1}.z = z; % indices of states
  260. MB{1}.x = x; % indices of states
  261. MB{1}.s = num2cell(-abs(diag(J))); % Lyapunov exponents
  262. MB{1}.U = U; % exogenous inputs
  263. MB{1}.V = V; % support in voxel space
  264. MB{1}.W = sparse(logical(abs(V))); % support in voxel space
  265. MB{1}.Y = Y; % time-series
  267. case 'blocking'
  268. %==================================================================
  269. % M B : B L O C K I N G
  270. %==================================================================
  271. MB = varargin{1};
  272. N = 4; % maximum number of scales
  273. m = ones(N,1); % # internal states
  275. % preclude long range coupling
  276. %------------------------------------------------------------------
  277. % R = cell(N,1);
  278. % R{1} = spm_mb_distance(MB{1}.V,MB{1}.XYZ);
  279. % R{1}(R{1} > 48) = 0;
  281. %% group renormalization (course graining or blocking)
  282. %------------------------------------------------------------------
  283. for i = 1:N
  285. % Markov blanket (particular) partition
  286. %--------------------------------------------------------------
  287. x = spm_Markov_blanket(MB{i}.J,MB{i}.z,m(i));
  289. % break if no states
  290. %--------------------------------------------------------------
  291. if size(x,2) < 1, break, end
  293. % dimension reduction (eliminating internal states)
  294. %--------------------------------------------------------------
  295. [J,z,v,s] = spm_A_reduce(MB{i}.J,x,OPT.T,OPT.N);
  296. MB{i + 1}.J = J;
  297. MB{i + 1}.z = z;
  298. MB{i + 1}.s = s;
  299. MB{i + 1}.x = x;
  300. MB{i + 1}.W = sparse(logical(abs(MB{i}.V)));
  302. % eigenmodes and variates (V and Y)
  303. %--------------------------------------------------------------
  304. for j = 1:numel(z)
  305. u = [MB{i + 1}.x{1:2,j}]; % blanket states
  306. w = MB{i + 1}.z{j}; % new states
  307. MB{i + 1}.V(:,w) = MB{i}.V(:,u)*v{j};
  308. MB{i + 1}.Y(:,w) = MB{i}.Y(:,u)*v{j};
  309. end
  311. % functional anatomy of states
  312. %--------------------------------------------------------------
  313. spm_figure('getwin',sprintf('Markov level %i',i)); clf;
  314. spm_mb_anatomy(MB{i},MB{1}.XYZ,MB{1}.VOX);
  316. % break if a single particle or parcel
  317. %--------------------------------------------------------------
  318. if size(MB{i}.x,2) < 2, break, end
  320. end
  324. %% Save results
  325. %==================================================================
  327. %-Save
  328. %------------------------------------------------------------------
  329. save(fullfile(MB{1}.swd,['MB_' MB{1}.con]),'MB');
  333. case 'results'
  334. %==================================================================
  335. % M B : R E S U L T S
  336. %==================================================================
  338. %-Get MB and ananlysis
  339. %------------------------------------------------------------------
  340. MB = varargin{1};
  341. analysis = varargin{2};
  342. if ischar(MB)
  343. MB = load(MB);
  344. MB = MB.MB;
  345. end
  347. switch analysis
  349. case('anatomy')
  351. % This figure illustrates the partition of states at the
  352. % first and subsequent scales. The upper right panel shows
  353. % all the constituent particles as a maximum intensity
  354. % projection, where the spatial support of each particle has
  355. % been colour-coded according to the variance explained by
  356. % its eigenmode. The upper middle panel shows the
  357. % corresponding adjacency matrix or coupling among particles.
  358. % The coloured circles encode the identity of each particle.
  359. % In this instance, the particles have been arranged in order
  360. % of descending dissipation (i.e., the principal eigenvalue
  361. % of each particles Jacobian). The upper right panel shows
  362. % these eigenvalues above the corresponding particle (encoded
  363. % by coloured dots) in terms of rate constants (i.e., the
  364. % negative inverse of the eigenvalue of each particle). The
  365. % remaining panels show the first 12 particles as maximum
  366. % intensity projections. The colour of the background
  367. % corresponds to the colour that designates each particle. In
  368. % this first level, each particle has a single eigenstate.
  369. % The numbers in brackets above each maximum intensity
  370. % projection correspond to the number of internal, active and
  371. % sensory states, respectively, where the active and sensory
  372. % states comprise blanket states. At this lowest level, every
  373. % eigenstate is a sensory state because it can influence, and
  374. % be influenced by,the eigenstates of other particles.
  376. % functional anatomy (space)
  377. %----------------------------------------------------------
  378. VOX = MB{1}.VOX;
  379. XYZ = MB{1}.XYZ;
  380. for i = 1:numel(MB)
  381. spm_figure('getwin',sprintf('Markov level %i',i)); clf;
  382. spm_mb_anatomy(MB{i},XYZ,VOX)
  383. end
  385. case('dynamics')
  387. % Intrinsic timescales in the brain: This figure reports
  388. % intrinsic timescales at intermediate scales. The left
  389. % column shows the eigenmodes in terms of their principal
  390. % frequency, i.e., the largest complex eigenvalue (divided by
  391. % 2*pi). The right column shows the corresponding eigenmodes
  392. % in terms of their principal time constants, i.e., the
  393. % reciprocal of the largest negative real part. These two
  394. % characterisations (principal frequency and time constant)
  395. % speak to different aspects of intrinsic timescales; both of
  396. % which contribute to the shape of an eigenstate’s
  397. % correlation function of time. The first quantifies the
  398. % frequency of solenoidal flow, while the second reflects the
  399. % rate of decay associated with the dissipative flow.
  401. % functional anatomy (time)
  402. %----------------------------------------------------------
  403. VOX = MB{1}.VOX;
  404. XYZ = MB{1}.XYZ;
  406. spm_figure('getwin','Timing'); clf;
  407. N = numel(MB) - 1;
  408. for i = 1:N
  410. % largest real and imaginary parts of each particle
  411. %------------------------------------------------------
  412. sr = zeros(size(MB{i}.V,2),1);
  413. si = zeros(size(MB{i}.V,2),1);
  414. for j = 1:numel(MB{i}.s)
  415. try
  416. z = MB{i}.z{j};
  417. sr(z) = max(real(MB{i}.s{j}));
  418. si(z) = max(imag(MB{i}.s{j}));
  419. end
  420. end
  422. % imaginary part (frequency)
  423. %------------------------------------------------------
  424. si = si/(2*pi);
  425. subplot(N,2,2*(i - 1) + 1),cla,hold off
  426. spm_mip(logical(abs(MB{i}.V))*si,XYZ,VOX);
  427. axis image, axis off
  428. str{1} = sprintf('Frequency - imaginary');
  429. str{2} = sprintf('%2.1f to %1.3f Hz',min(si),max(si));
  430. title(str,'FontSize',14,'FontWeight','bold')
  432. % real part (time constants)
  433. %------------------------------------------------------
  434. sr = (-1./sr);
  435. subplot(N,2,2*(i - 1) + 2),cla,hold off
  436. spm_mip(logical(abs(MB{i}.V))*sr,XYZ,VOX);
  437. axis image, axis off
  438. str{1} = sprintf('Dissipation - real (scale %i)',i);
  439. str{2} = sprintf('%2.1f to %2.1f sec',min(sr),max(sr));
  440. title(str,'FontSize',14,'FontWeight','bold')
  442. end
  444. case('distance')
  446. % Local connectivity: This figure reports some of the
  447. % statistical characteristics of dynamical coupling among
  448. % particles at the first level. The upper left panel plots
  449. % each connection in terms of the real part of the
  450. % corresponding Jacobian in Hz, against the distance spanned
  451. % by the connection (i.e., Euclidean distance between the
  452. % centres of the two particles). The upper left panel plots
  453. % the log-coupling (real part) against log-distance. The
  454. % lower panel plots the strength of reciprocal connections
  455. % against each other, to illustrate the relative proportions
  456. % of recurrent excitatory and inhibitory coupling. The
  457. % rarefied region in the centre of this scatterplot reflects
  458. % the fact that connections wIf youith small coupling strengths
  459. % have been eliminated during Bayesian model reduction.
  461. % spatial distance and (first-level) coupling
  462. %----------------------------------------------------------
  463. spm_figure('getwin','Distance rules'); clf;
  464. col = {'r','b'};
  466. % get distance (from averge location, xyz)
  467. %----------------------------------------------------------
  468. J = MB{1}.J;
  469. nv = length(J);
  470. D = abs(MB{1}.V);
  471. D = bsxfun(@rdivide,D,sum(D));
  472. xyz = MB{1}.XYZ*D;
  473. xyz(1,:) = abs(xyz(1,:));
  475. % characterise excitation and inhibitory coupling separately
  476. %----------------------------------------------------------
  477. u = 1/1024;
  478. for e = 1:2
  479. D = [];
  480. E = [];
  481. for i = 1:nv
  482. for j = 1:nv
  483. d = xyz(:,i) - xyz(:,j);
  484. d = sqrt(d'*d);
  485. if e > 1
  486. if J(i,j) < -u && d
  487. D(end + 1,1) = sqrt(d'*d);
  488. E(end + 1,1) = abs(J(i,j));
  489. end
  490. else
  491. if J(i,j) > +u && d
  492. D(end + 1,1) = sqrt(d'*d);
  493. E(end + 1,1) = abs(J(i,j));
  494. end
  495. end
  496. end
  497. end
  499. % coupling (cumulative)
  500. %----------------------------------------------------------
  501. n = 16;
  502. c = zeros(n,1);
  503. k = zeros(n,1);
  504. d = linspace(8,max(D),n + 1);
  505. for i = 1:n
  506. j = D > d(i) & D < d(i + 1);
  507. c(i,1) = mean(E(j)); % mean coupling
  508. k(i,1) = sum(j); % numer of edges
  509. end
  510. d = spm_vec(d(1:n));
  511. j = k > mean(k)/4;
  512. d = d(j);
  513. c = c(j);
  514. k = k(j);
  516. % regression
  517. %----------------------------------------------------------
  518. [F,df,B] = spm_ancova([ones(size(d)) log(d)],[],log(c));
  520. subplot(2,2,1)
  521. plot(D,E,['.' col{e}],'MarkerSize',1), hold on
  522. plot(d,c,['o' col{e}],'MarkerSize',8)
  523. plot(d,exp(B(1))*d.^B(2),col{e},'LineWidth',2)
  524. title('Distance scaling rule','FontSize',16)
  525. xlabel('Distance (mm)'), ylabel('coupling strength (Hz}')
  526. axis square, box off, set(gca,'YLim',[0, 1])
  528. subplot(2,2,2)
  529. plot(log(d),log(c),['o' col{e}],'MarkerSize',8), hold on
  530. plot(log(d),B(1) + B(2)*log(d),col{e},'LineWidth',1)
  531. if e == 1
  532. str = sprintf('Log scaling (%2.2f',B(2));
  533. else
  534. str = [str sprintf(', %2.2f)',B(2))];
  535. end
  536. title(str,'FontSize',16)
  537. xlabel('log distance'), ylabel('log coupling')
  538. axis square, box off
  539. legend('excitatory',' ','inhibitory',' '),legend(gca,'boxoff')
  541. % subplot(2,2,3)
  542. % p = k./(d.^2);
  543. % p = p/sum(p);
  544. % plot(d,p,col{e}), hold on
  545. % title('Density','FontSize',16)
  546. % xlabel('distance (mm)'), ylabel('Density')
  547. % axis square, box off
  549. end
  551. % assymetry
  552. %----------------------------------------------------------
  553. D = [];
  554. E = [];
  555. for i = 1:nv
  556. for j = 1:nv
  557. if abs(J(i,j)) > exp(-8) && j ~= i
  558. D(end + 1,1) = J(j,i);
  559. E(end + 1,1) = J(i,j);
  560. end
  561. end
  562. end
  564. subplot(2,1,2)
  565. hold off, plot(E,D,'.r','MarkerSize',1)
  566. title('Reciprocal coupling','FontSize',16)
  567. xlabel('coupling (Hz)'), ylabel('coupling (Hz)')
  568. axis square, box off, axis([-1 1 -1 1])
  570. % percentages of different recurrent connections
  571. %----------------------------------------------------------
  572. R(1) = sum(D > 0 & E > 0);
  573. R(2) = sum(D < 0 & E > 0);
  574. R(3) = sum(D > 0 & E < 0);
  575. R(4) = sum(D < 0 & E < 0);
  576. R = 100*R/sum(R);
  578. text( 1/4, 1/4,sprintf('%2.0f p.c.',R(1)),'FontWeight','bold');
  579. text( 1/4,-1/4,sprintf('%2.0f p.c.',R(2)),'FontWeight','bold');
  580. text(-1/4, 1/4,sprintf('%2.0f p.c.',R(3)),'FontWeight','bold');
  581. text(-1/4,-1/4,sprintf('%2.0f p.c.',R(4)),'FontWeight','bold');
  584. case('scaling')
  586. % Scale invariance: This figure illustrates scaling behaviour
  587. % across the scales of the particular decomposition. The
  588. % upper panel plots the real part of the eigenvalues of each
  589. % particle against its spatial scale; namely the Calliper
  590. % width of the particle’s eigenmode. This is replicated for
  591. % each of the four scales, denoted by the different colours
  592. % (green, pink, cyan and puce, respectively). The expected
  593. % values are shown as encircled large dots. The lower left
  594. % panel plots the logarithms of these temporal and spatial
  595. % expectations against each other. The resulting regression
  596. % slope corresponds to the scaling exponent. The light grey
  597. % circles correspond to what would have been seen at higher
  598. % and lower scales. The lower right plot shows the same
  599. % regression in terms of the implicit time constant, as a
  600. % function of spatial scale expressed in millimetres.
  602. % time constants over scales
  603. %----------------------------------------------------------
  604. spm_figure('getwin','Scaling'); clf;
  606. N = numel(MB);
  607. SS = cell(N,1);
  608. SR = cell(N,1);
  609. for i = 1:N
  610. for j = 1:numel(MB{i}.s)
  612. % get spatial scale (standard deviation)
  613. %--------------------------------------------------
  614. S = sum(logical(abs(MB{i}.V(:,MB{i}.z{j})))).^(1/3);
  615. SS{i} = [SS{i}; mean(S)];
  617. % get temporal scale (eigenvalue)
  618. %--------------------------------------------------
  619. S = real(MB{i}.s{j});
  620. SR{i} = [SR{i}; mean(S)];
  621. end
  622. end
  625. % Expectations and graphics
  626. %----------------------------------------------------------
  627. for i = 1:N
  628. j = isfinite(SS{i});
  629. SS{i} = SS{i}(j);
  630. ss(i,1) = mean(SS{i});
  631. j = isfinite(SR{i});
  632. SR{i} = SR{i}(j);
  633. sr(i,1) = mean(SR{i});
  634. end
  636. subplot(2,1,1), hold off
  637. col = spm_MB_col(N);
  638. for i = 1:N
  639. plot(SS{i},SR{i},'.','MarkerSize', 8,'Color',col{i}), hold on
  640. plot(ss(i),sr(i),'.','MarkerSize',32,'Color',col{i})
  641. plot(ss(i),sr(i),'o','MarkerSize',32,'Color',col{i})
  642. end
  643. title('Scale invariance','FontSize',16)
  644. xlabel('spatial scale (mm)'), ylabel('eigenvalue (Hz)')
  645. axis square, box off, hold off
  647. % regression (Time)
  648. %----------------------------------------------------------
  649. ii = (1:N)'; % scale
  650. [F,df,Br] = spm_ancova([ones(N,1) ii], [],-log(-sr));
  651. [F,df,Bs] = spm_ancova([ones(N,1) ii], [], log( ss));
  652. [F,df,A ] = spm_ancova([ones(N,1) log(ss)],[],-log(-sr));
  653. A = [A(1), Br(2)/Bs(2)];
  655. fprintf('Beta (exp) time: %2.2f (%2.2f)\n',Br(2),exp(Br(2)));
  656. fprintf('Beta (exp) space: %2.2f (%2.2f)\n',Bs(2),exp(Bs(2)));
  657. fprintf('Alpha: %2.2f\n',A(2));
  659. % extrapolate
  660. %----------------------------------------------------------
  661. bs = exp(Bs(2)*(-4:4:8) + Bs(1))';
  662. br = exp(A(1) + A(2)*log(bs));
  663. S = linspace(min(bs),max(bs),64);
  665. subplot(2,2,3), cla, hold on
  666. plot(log(bs),log(br),'.','MarkerSize',32,'Color',[1,1,1]*.8)
  667. plot(log(bs),log(br),'o','MarkerSize',32,'Color',[1,1,1]*.8)
  668. for i = 1:N
  669. plot(log(SS{i}),-log(-SR{i}),'.','MarkerSize',4,'Color',col{i})
  670. plot(log(ss(i)),-log(-sr(i)),'.','MarkerSize',32,'Color',col{i})
  671. plot(log(ss(i)),-log(-sr(i)),'o','MarkerSize',32,'Color',col{i})
  672. hold on
  673. end
  674. plot(log(S),A(1) + A(2)*log(S))
  676. title(sprintf('Log scaling - %2.2f',A(2)),'FontSize',16)
  677. xlabel('log distance (mm)'), ylabel('log time constant')
  678. axis square, box off, hold off
  680. subplot(2,2,4), hold off
  681. plot(S,exp(A(1) + A(2)*log(S))), hold on
  682. plot(bs,br,'.','MarkerSize',32,'Color',[1,1,1]*.8)
  683. plot(bs,br,'o','MarkerSize',32,'Color',[1,1,1]*.8)
  684. title('Scaling behaviour','FontSize',16)
  685. xlabel('distance (mm)'), ylabel('time constant (sec)')
  686. axis square, box off
  688. % some canonical scales
  689. %----------------------------------------------------------
  690. disp('mm')
  691. disp(num2str(bs,'%2.2e'))
  692. disp('Seconds')
  693. disp(num2str(br,'%2.2e'))
  694. for i = 1:numel(bs)
  695. hold on, plot([bs(i),bs(i)],[0 br(i)],'r-.')
  696. end
  698. case('kernels')
  700. % Transfer functions: This figure characterises the dynamics
  701. % at successive scales in terms of transfer functions, as
  702. % quantified by the complex eigenvalues (c.f., a pole-zero
  703. % map). The left column shows the transfer functions over
  704. % frequency for all particles with a complex eigenvalue (at
  705. % successive scales). These eigenvalues are shown in the
  706. % right column in the complex number plane. As we ascend from
  707. % one scale to the next, the real part of the eigenvalue
  708. % approaches zero from the left and the number of eigenvalues
  709. % falls with the coarse-graining, i.e., number of particles.
  710. % The complex part of the eigenvalues corresponds to the peak
  711. % frequency of the associated transfer functions, while the
  712. % dispersion around this peak decreases as the real part
  713. % approaches zero. The emergence of spectral peaks in the
  714. % transfer functions inherit from the complex part of the
  715. % eigenvalues that emerge under asymmetric coupling with
  716. % solenoidal flow.
  718. % frequencies (Hz)
  719. %----------------------------------------------------------
  720. w = linspace(0,1/8,128);
  722. % transfer functions over scales
  723. %----------------------------------------------------------
  724. spm_figure('getwin','Transfer functions'); clf;
  725. for i = 2:min(4,numel(MB))
  727. % loop over particles
  728. %------------------------------------------------------
  729. nz = numel(MB{i}.z);
  730. col = spm_MB_col(nz);
  731. for p = 1:nz
  733. % loop over eigenmodes
  734. %--------------------------------------------------
  735. for j = 1:numel(MB{i}.s{p})
  737. % condition unstable eigenmodes
  738. %----------------------------------------------
  739. s = MB{i}.s{p}(j);
  740. if real(s) > 0
  741. s = 1j*imag(s) - 1/64;
  742. end
  744. % Transfer functions (if imaginary)
  745. %----------------------------------------------
  746. if imag(s) ~= 0
  748. % Transfer function
  749. %------------------------------------------
  750. S = 1./(1j*2*pi*w - s);
  752. subplot(3,2,(i - 2)*2 + 1), hold on
  753. plot(w,abs(S),'color',col{p})
  755. % Pole-zero map
  756. %------------------------------------------
  757. subplot(3,2,(i - 2)*2 + 2), hold on
  758. plot(MB{i}.s{p}(j),'.','MarkerSize',32,'color',col{p})
  761. end
  762. end
  763. end
  765. subplot(3,2,(i - 2)*2 + 1)
  766. title(sprintf('Scale %i',i - 1),'fontsize',16)
  767. xlabel('frequency (Hz)'), ylabel('transfer tunctions')
  768. axis square, box off, spm_axis tight, hold off
  770. subplot(3,2,(i - 2)*2 + 2)
  771. plot([0,0],[-1,1],':')
  772. title('Eigenvalues','fontsize',16)
  773. xlabel('real (Hz)'), ylabel('imaginary (Hz)')
  774. axis square, box off, axis([-1 1 -1/8 1/8]), hold off
  776. end
  779. case('connectivity')
  781. % Dynamic coupling: This figure characterises the coupling
  782. % between the two eigenstates, which (unless specified)
  783. % correspond to the strongest coupling at this highest scale.
  784. % This coupling is mediated by the corresponding element of
  785. % the (complex) Jacobian, circled in red in the upper middle
  786. % panel. The flanking panels on the left and right show the
  787. % corresponding eigenmodes in voxel space. The middle row
  788. % shows the auto-covariance functions of the two eigenstates,
  789. % illustrating serial correlations that can last for many
  790. % seconds. The lower panels report the cross-covariance
  791. % function between the two eigenstates, over 256 seconds
  792. % (lower left panel) and focusing on cross-covariances over
  793. % 32 seconds (lower right panel). Asymmetrical
  794. % cross-covariance (and implicitly cross-correlation)
  795. % function reflects the solenoidal coupling and the implicit
  796. % breaking of detailed balance the particular decomposition
  797. % accommodates.
  800. % functional anatomy
  801. %----------------------------------------------------------
  802. VOX = MB{1}.VOX;
  803. XYZ = MB{1}.XYZ;
  804. spm_figure('getwin',' extrinsic connectivity'); clf;
  805. spm_mb_anatomy(MB{end},XYZ,VOX)
  808. % Jacobian and cross covariance function
  809. %----------------------------------------------------------
  810. dfdx = full(MB{end}.J);
  811. [ccf,pst] = spm_ssm2ccf(dfdx);
  813. % get pair of eigenvariates
  814. %----------------------------------------------------------
  815. try
  816. i = varargin{3};
  817. j = varargin{4};
  818. catch
  819. A = dfdx;
  820. A = A - diag(diag(A));
  821. [d,i] = max(real(A(:)));
  822. [i,j] = ind2sub(size(A),i);
  823. end
  824. cf = ccf(:,i,j);
  825. lim = [min(cf) max(cf)];
  827. % if no correlations than plot conjugate modes
  828. %----------------------------------------------------------
  829. if max(abs(lim)) < 1e-6
  830. [i,j] = sort(imag(spm_vec(MB{end}.s)),'descend');
  831. i = j(1);
  832. j = j(end);
  833. end
  835. % graphics
  836. %----------------------------------------------------------
  837. subplot(3,2,3), plot(pst,ccf(:,i,i))
  838. title(sprintf('covariance: state %i',i),'FontSize',16)
  839. xlabel('lag (sec)'),ylabel('covariance')
  840. axis square, box off, spm_axis tight
  842. subplot(3,2,4), plot(pst,ccf(:,j,j))
  843. title(sprintf('covariance: state %i',j),'FontSize',16)
  844. xlabel('lag (sec)'),ylabel('covariance')
  845. axis square, box off, spm_axis tight
  847. subplot(3,2,5), plot(pst,ccf(:,i,j),[0,0],lim,':b')
  848. title(sprintf('cross-covariance: states %i,%i',i,j),'FontSize',16)
  849. xlabel('lag (sec)'),ylabel('covariance')
  850. if max(abs(lim)) < 1e-6, set(gca,'YLim',[-1,1]), end
  851. axis square, box off, spm_axis tight
  853. % blow-up around zero lag
  854. %----------------------------------------------------------
  855. [d,m] = max(abs(ccf(:,i,j)));
  856. m = pst(m);
  858. subplot(3,2,6), plot(pst,ccf(:,i,j),[0,0],lim,':b',[m,m],lim,'-.r')
  859. title(sprintf('cross-covariance: states %i,%i',i,j),'FontSize',16)
  860. xlabel('lag (sec)'),ylabel('covariance')
  861. if max(abs(lim)) < 1e-6, set(gca,'YLim',[-1,1]), end
  862. axis square, box off, spm_axis tight, set(gca,'XLim',[-1,1]*32)
  865. % modes
  866. %----------------------------------------------------------
  867. subplot(3,3,1),cla,hold off
  868. spm_mip(abs(MB{end}.V(:,i)),XYZ,VOX);
  869. axis image, axis off
  870. title(sprintf('Eigenmode %i',i),'FontSize',16)
  872. subplot(3,3,3),cla,hold off
  873. spm_mip(abs(MB{end}.V(:,j)),XYZ,VOX);
  874. axis image, axis off
  875. title(sprintf('Eigenmode %i',j),'FontSize',16)
  877. subplot(3,3,2), hold on, plot(j,i,'or','MarkerSize',16)
  879. case('eigenmodes')
  881. % Dissipative and solenoidal dynamics: This figure unpacks
  882. % the intrinsic coupling at the final (usually, single
  883. % particle) level. At this level, there are can be no
  884. % coupling between particles and (by construction) the
  885. % dynamics is completely characterised in terms of the
  886. % eigenstates that comprise the particle. In turn, these are
  887. % completely characterised by their complex eigenvalues;
  888. % namely, the intrinsic complex coupling. The upper panels
  889. % show the dissipative and solenoidal (kinetic) energy of the
  890. % eigenstates that comprise the particle. The corresponding
  891. % eigenmodes are shown in the subsequent panels as maximum
  892. % intensity projections (of their absolute values).
  894. % functional anatomy
  895. %----------------------------------------------------------
  896. spm_figure('getwin',' extrinsic coupling'); clf;
  897. VOX = MB{1}.VOX;
  898. XYZ = MB{1}.XYZ;
  901. % Jacobian
  902. %----------------------------------------------------------
  903. J = full(MB{end}.J);
  904. Y = full(MB{end}.Y);
  906. % complex form
  907. %----------------------------------------------------------
  908. J = diag(J);
  909. D = -real(J); % dissipative
  910. Q = (imag(J).^2)./D; % solenoidal
  912. % show dynamics
  913. %==========================================================
  914. try, dt = MB{1}.dt; catch, dt = 1; end
  915. try, i = varargin{3}; catch, [d,i] = max(Q); end
  918. subplot(4,2,1), bar(D)
  919. xlabel('eigenmode'),ylabel('kinetic energy'), axis square, box off
  920. title('Dissipative energy','FontSize',16)
  922. subplot(4,2,2), bar(Q)
  923. xlabel('eigenmode'),ylabel('kinetic energy'), axis square, box off
  924. title('Solenoidal energy','FontSize',16)
  926. % pst = (1:size(Y,1))*dt;
  927. % subplot(3,2,4), plot(pst,abs(Y(:,i)),'-',pst,angle(Y(:,i)),':')
  928. % xlabel('time'),ylabel('state'), axis square, box off
  929. % legend({'amplitude','angle'}), legend(gca,'boxoff')
  930. % title('Dissipative dynamics','FontSize',16)
  932. % mode
  933. %----------------------------------------------------------
  934. for i = 1:min(9,numel(D))
  935. subplot(4,3,i + 3),cla,hold off
  936. spm_mip(abs(MB{end}.V(:,i)),XYZ,VOX);
  937. axis image, axis off
  938. title(sprintf('Eigenmode %i',i),'FontSize',12,'FontWeight','bold')
  939. end
  942. case('responses')
  944. % Induced responses: This figure illustrates the expression
  945. % of experimental or condition-specific effects at different
  946. % scales of the particular decomposition. The top panel is an
  947. % unusual form of statistical parametric mapping; namely, and
  948. % image of the F statistic, testing for the significance of
  949. % an effect of any of any exogenous inputs. Each row of the F
  950. % statistic map corresponds to a scale and comprises the F
  951. % statistic for each successive particle at that scale. This
  952. % map shows that the effect of (some linear mixture of)
  953. % exogenous inputs can be detected at several scales, as
  954. % evidenced by the dark bars. unless specified, the most
  955. % significant eigenmode is shown on the lower left in voxel
  956. % space. It's expression over time (in terms of its real
  957. % value) is depicted in the middle panel (blue line), with
  958. % the best fitting prediction based upon exogenous input
  959. % (green line). This prediction is a contrast (i.e., linear
  960. % mixture) of the input functions shown in the design matrix
  961. % on the lower right. The coefficients of this contrast are
  962. % shown below the design matrix. The last column of the
  963. % design matrix is simply a column of ones.
  965. % functional anatomy
  966. %----------------------------------------------------------
  967. VOX = MB{1}.VOX;
  968. XYZ = MB{1}.XYZ;
  969. spm_figure('getwin','responses'); clf;
  971. % ANOVA
  972. %----------------------------------------------------------
  973. U = MB{1}.U;
  974. n = size(U,2);
  975. U(:,n + 1) = 1;
  976. c = spm_speye(n + 1,n);
  978. % search over scles and eigenmodes
  979. %----------------------------------------------------------
  980. for j = 1:numel(MB)
  981. for i = 1:size(MB{j}.Y,2)
  982. Y = real(MB{j}.Y(:,i));
  983. F(j,i) = spm_ancova(U,[],Y,c);
  984. end
  985. end
  987. % most significant mode
  988. %----------------------------------------------------------
  989. [Fj,j] = max(F,[],2);
  990. [Fi,i] = max(Fj);
  991. try, i = varargin{3}; end
  992. j = j(i);
  993. Y = real(MB{i}.Y(:,j));
  994. [Fi,df,beta] = spm_ancova(U,[],Y,c);
  995. p = 1 - spm_Fcdf(Fi,df);
  996. df = round(df);
  998. % show dynamics
  999. %==========================================================
  1000. try, dt = MB{1}.dt; catch, dt = 1; end
  1001. pst = (1:size(Y,1))*dt;
  1002. str = 'Eigenvariate and prediction';
  1003. str = sprintf('%s: F(%d,%d) = %2.1f, p = %2.3f',str,df(1),df(2),Fi,p);
  1005. subplot(3,1,1), imagesc(1 - F)
  1006. xlabel('eigenmode'),ylabel('scale'), spm_axis tight, box off
  1007. title('F-statistic','FontSize',16), hold on
  1008. plot(j,i,'or','MarkerSize',8)
  1010. subplot(3,1,2), plot(pst,Y,':',pst,U*beta)
  1011. xlabel('time'),ylabel('real part'), spm_axis tight, box off
  1012. title(str,'FontSize',16)
  1014. % mode
  1015. %----------------------------------------------------------
  1016. subplot(3,2,5),cla,hold off
  1017. spm_mip(abs(MB{i}.V(:,j)),XYZ,VOX);
  1018. axis image, axis off
  1019. title(sprintf('Eigenmode %i',j),'FontSize',16)
  1021. % inputs and contribution
  1022. %----------------------------------------------------------
  1023. subplot(6,2,10), imagesc(U)
  1024. xlabel('input'),ylabel('time'), axis square, box off
  1025. title('Inputs','FontSize',16)
  1027. subplot(6,2,12), bar(beta)
  1028. xlabel('input'),ylabel('contribution'), axis square, box off
  1029. spm_axis tight
  1032. case('connectogram')
  1034. % Extrinsic connectivity: This figure illustrates asymmetric
  1035. % extrinsic (between particle) coupling at the final scale
  1036. % supplied in MB. The upper panels reproduce the results in
  1037. % 'anatomy', while the lower panel is a connectogram
  1038. % illustrating the coupling among eigenstates that constitute
  1039. % the particles at this scale. The width of each connector
  1040. % reflects the strength of the coupling; after dividing the
  1041. % strength into five bins and eliminating the lowest bin. The
  1042. % colour of the dots corresponds to the colour of the
  1043. % particle in the upper right panel. The colour of the
  1044. % connectors corresponds to the source of the strongest
  1045. % (reciprocal) connection. The coupling strength corresponds
  1046. % to the real part of the Jacobian, in Hz.
  1048. % functional anatomy of penultimate levels
  1049. %----------------------------------------------------------
  1050. VOX = MB{1}.VOX;
  1051. XYZ = MB{1}.XYZ;
  1052. i = numel(MB);
  1053. spm_figure('getwin',sprintf('connectogram level %i',i)); clf;
  1054. spm_mb_anatomy(MB{i},XYZ,VOX,0)
  1057. % assign colours to connections (afferent
  1058. %----------------------------------------------------------
  1059. MB = MB{i};
  1060. nz = numel(MB.z); % number of vector states
  1061. col = spm_MB_col(nz); % colours
  1062. COL = [];
  1063. for i = 1:nz
  1064. for j = 1:numel(MB.s{i})
  1065. COL(end + 1,:) = col{i}';
  1066. end
  1067. end
  1069. % connectivity of undirected connectivity
  1070. %----------------------------------------------------------
  1071. subplot(2,1,2)
  1072. spm_circularGraph(MB.J,'colormap',COL)
  1073. axis square off
  1074. title('Connectogram','FontSize',16)
  1076. case('inputs')
  1078. % Induced responses over space and time: this figure
  1079. % characterises induced responses in terms of first order
  1080. % Volterra kernels (i.e., impulse response functions) of
  1081. % particles at the first (finest) scale of course graining.
  1082. % Each row corresponds to the inputs specified. The left
  1083. % column shows the expression of these inputs over particles
  1084. % (weighted by the absolute value of their eigenmodes). This
  1085. % effect is the variance attributable to each input (i.e.,
  1086. % square of the corresponding kernel, summed over time), shown
  1087. % in the left row. These kernels are shown for the 32
  1088. % particles with the greatest (absolute) magnitude.
  1090. % functional anatomy (1st order kernels)
  1091. %----------------------------------------------------------
  1092. spm_figure('getwin','input effects'); clf;
  1093. VOX = MB{1}.VOX;
  1094. XYZ = MB{1}.XYZ;
  1095. nu = size(MB{1}.U,2);
  1097. dfdx = MB{1}.J;
  1098. dfdu = MB{1}.K;
  1099. pst = (1:64)*4;
  1100. ker = spm_ssm2ker(dfdx,dfdu,[],pst);
  1101. for i = 1:nu
  1103. % spatial expression
  1104. %-------------------------------------------------------
  1105. subplot(nu,2,2*(i - 1) + 1)
  1106. K = sum(ker(:,:,i).^2,1);
  1108. spm_mip(abs(MB{1}.V*K(:)),XYZ,VOX);
  1109. axis image, axis off
  1110. title(MB{1}.name{i},'FontSize',16)
  1112. % temporal expression
  1113. %-------------------------------------------------------
  1114. [d,j] = sort(abs(ker(1,:,i)),'descend');
  1116. subplot(nu,2,2*(i - 1) + 2)
  1117. plot(pst,ker(:,j(1:64),i))
  1118. xlabel('time (secs)'),ylabel('response'), axis square, box off
  1119. title('Neuronal responses','FontSize',16)
  1121. end
  1123. otherwise
  1124. error('Unknown action.');
  1126. end
  1129. otherwise
  1130. error('Unknown action.');
  1132. end
  1134. %% subroutines
  1135. %==========================================================================
  1137. % functional anatomy - distance (R)
  1138. %==========================================================================
  1139. function [R,S] = spm_mb_distance(V,XYZ)
  1140. % FORMAT [R,S] = spm_mb_distance(V,XYZ)
  1141. % V - pattern vectors
  1142. % XYZ - locations
  1143. %
  1144. % R - distance between patterns
  1145. % S - spatial dispersion of patterns (sd; mm)
  1146. %__________________________________________________________________________
  1148. % distance matrix (mm)
  1149. %--------------------------------------------------------------------------
  1150. nv = size(V,2);
  1151. R = zeros(nv,nv);
  1152. D = abs(V);
  1153. D = bsxfun(@rdivide,D,sum(D));
  1154. xyz = XYZ*D;
  1155. S = (XYZ.^2)*D - (XYZ*D).^2;
  1156. S = sqrt(mean(S(:)));
  1158. % interhemispheric coupling
  1159. %--------------------------------------------------------------------------
  1160. xyz(1,:) = abs(xyz(1,:));
  1162. for i = 1:nv
  1163. for j = i:nv
  1164. d = xyz(:,i) - xyz(:,j);
  1165. d = sqrt(d'*d);
  1166. R(i,j) = d;
  1167. R(j,i) = d;
  1168. end
  1169. end
  1172. % functional anatomy
  1173. %==========================================================================
  1174. function spm_mb_anatomy(MB,XYZ,VOX,N)
  1175. %__________________________________________________________________________
  1178. % set-up
  1179. %--------------------------------------------------------------------------
  1180. if nargin < 4, N = 12; end
  1182. nz = numel(MB.z); % number of particles
  1183. ns = size(MB.W,2); % number of states (n - 1)
  1184. nx = size(MB.J,2); % number of states (n)
  1185. col = spm_MB_col(nz); % colours
  1187. % functional anatomy of states
  1188. %--------------------------------------------------------------------------
  1189. subplot(3,3,1),cla,hold off
  1190. c = var(abs(MB.Y));
  1191. spm_mip(logical(abs(MB.V))*c',XYZ,VOX);
  1192. axis image, axis off
  1193. str{1} = sprintf('%d states, %d particles',ns,nz);
  1194. str{2} = sprintf('(%d eigenstates)',nx);
  1195. title(str,'FontSize',16)
  1197. subplot(3,3,2),cla,hold off
  1198. if nx > 128
  1199. imagesc(~abs(MB.J))
  1200. else
  1201. J = real(MB.J);
  1202. imagesc(max(J,-1/8))
  1203. end
  1204. axis square
  1205. title({'Jacobian';''},'FontSize',16)
  1206. xlabel('States'), hold on
  1208. subplot(3,3,3),cla,hold off
  1209. bar(-1./real(spm_vec(MB.s)))
  1210. axis square, set(gca,'XLim',[0,nx + 1])
  1211. title('Time constants','FontSize',16)
  1212. xlabel('states'),ylabel('seconds')
  1214. % label vector states (states of particles previous level)
  1215. %--------------------------------------------------------------------------
  1216. xi = 1;
  1217. for i = 1:nz
  1218. for j = 1:numel(MB.s{i})
  1219. hold on, subplot(3,3,2), plot(xi,0,'.','Color',col{i},'MarkerSize',24)
  1220. hold on, subplot(3,3,3), plot(xi,0,'.','Color',col{i},'MarkerSize',24)
  1221. xi = xi + 1;
  1222. end
  1223. end
  1226. % functional anatomy of particles
  1227. %--------------------------------------------------------------------------
  1228. nn = min(nz,N);
  1229. for j = 1:nn
  1231. Z = any(MB.W(:,MB.x{1,j}),2)*2; % active states
  1232. Z = any(MB.W(:,MB.x{2,j}),2)*1 + Z; % sensory states
  1233. Z = any(MB.W(:,MB.x{3,j}),2)*3 + Z; % internal states
  1236. % get subplot layout
  1237. %----------------------------------------------------------------------
  1238. if nn == 1
  1239. ni = 3; nj = 1;
  1240. elseif nn == 2
  1241. ni = 3; nj = 2;
  1242. elseif nn == 3
  1243. ni = 3; nj = 3;
  1244. elseif nn == 4
  1245. ni = 3; nj = 2;
  1246. elseif nn == 5
  1247. ni = 4; nj = 2;
  1248. elseif nn == 6
  1249. ni = 4; nj = 3;
  1250. else
  1251. ni = 4; nj = 4;
  1252. end
  1254. subplot(ni,nj,j + nj)
  1255. [z,i] = max(Z);
  1256. if z < 3, Z(i) = 3; end
  1257. mip = spm_mip(Z,XYZ,VOX);
  1258. for k = 1:3
  1259. c = col{j}(k)/2;
  1260. MIP(:,:,k) = (1 - mip/64)*(1 - c) + c;
  1261. end
  1262. image(MIP), axis image, axis off
  1263. str{1} = sprintf('Particle %d of %d (%d,%d,%d)',j,nz,...
  1264. numel(MB.x{3,j}),numel(MB.x{1,j}),numel(MB.x{2,j}));
  1265. str{2} = sprintf('%d eigenmodes',numel(MB.z{j}));
  1266. title(str,'FontWeight','Bold')
  1267. end

spm_mb_ui_mod.m at commit cdab728, under AGPL-3.0 · at the source

Overview

Authors: Tahereh S Zarghami1
ORCID iDs: Tahereh S Zarghami
  1. Bio-Electric Department, School of Electrical and Computer Engineering, University of Tehran, Tehran, Iran
Institutions: University of Tehran (Iran)
Journal: Imaging neuroscience (Cambridge, Mass.), volume 4, article IMAG.a.1234
Dates: received 26 June 2025; accepted 8 April 2026; published online 22 May 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1162/imag.a.1234 · PMID 42206220 · PMCID PMC13206500 · OpenAlex W4411325862
Open access: diamond, a free copy (OpenAlex)
Preprint: doi:10.1101/2025.06.14.659698
Status: code verified
Categories: computational modeling (no new data) (modality), computational (subfield)
Methods: Statistics, Machine learning
Keywords: dynamic causal modeling, multiscale analysis, Bayesian model reduction, structure learning, minimum cut, graph theory, Markov blanket, scale invariance
Topic: Functional Brain Connectivity Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Citations: not cited yet (Europe PMC); 196 references in the paper

Abstract

The hierarchical organization of the brain’s distributed network has received growing interest from the neuroscientific community, largely because of its potential to enhance our understanding of human cognition and behavior, in health and disease. This interest is motivated by the hypothesis that near-critical brain dynamics enable multiscale integration and segregation of neural dynamics. While most multiscale connectivity analyses focus on structural and functional networks, characterizing the effective connectome across multiple scales has been somewhat overlooked—primarily for computational reasons. The difficulty of estimating large cyclic causal models, together with the scarcity of theoretical frameworks for systematically moving between scales, has hindered progress in this direction. This technical note introduces a top–down multiscale parcellation scheme for dynamic causal models, with application to neuroimaging data. The method is based on Bayesian model comparison, as a generalization of the well-known ΔBIC method. To facilitate computation, recent developments in linear dynamic causal modeling (DCM) and Bayesian model reduction (BMR) are deployed. Specifically, a naïve version of BMR is introduced, enabling the parcellation scheme to scale to hundreds or thousands of regions. Notably, the derivations reveal an analytical relationship between reduced model evidence and minimum cut problem in graph theory. This duality puts the tools of graph theory at the service of model evidence optimization and significance testing. The proposed method was applied to simulated and empirical causal models to establish face and construct validity. Consequently, the large empirical causal network, inferred from a neuroimaging dataset, exhibited log–log scaling trends, suggestive of scale invariance in multiple dynamical measures. Future generalizations of this technique and its potential applications in systems and clinical neuroscience are discussed.

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 9 matches between paragraphs and lines of code.

tszarghami/Parcellation

License: AGPL-3.0
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: cdab728b3b1cf4541279f2d1165fdc0e171c780c, 14 May 2026
Languages: MATLAB (22)
Size: 46 files, 22 scripts
Software Heritage: not archived
Found in: “Data and Code Availability”
Holds: README, license file
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: SPM (8 files), Statistics and Machine Learning Toolbox (3 files)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
24 files
At the source: github.com/tszarghami/Parcellation

The paper's code and data availability statement is in the Data section.

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;
  • 22 scripts, each with its path and the digest of its content;
  • 9 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.

Data and Code Availability

MATLAB code to reproduce the results presented in this paper is available at: https://github.com/tszarghami/Parcellation. Linearized DCM, PEB, and BMR procedures are implemented using MATLAB routines (spm_dcm_J.m, spm_dcm_peb.m, spm_dcm_bmr.m) in SPM12 software package (https://www.fil.ion.ucl.ac.uk/spm/). Spectral clustering is implemented using the Compressive Spectral Clustering Toolbox (http://cscbox.gforge.inria.fr/), and graph randomization is performed using functions from the Brain Connectivity Toolbox (https://sites.google.com/site/bctnet/). The attention-to-visual-motion fMRI dataset analyzed in this study is publicly available at: https://www.fil.ion.ucl.ac.uk/spm/data/attention/. Appendices are provided as Supplementary Material (https://doi.org/10.1162/IMAG.a.1234#supplementary-data).

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

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 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 1 author, 8 keywords, 182 references.

Cite

This paper

Zarghami, T. S. (2026). Multiscale parcellation of dynamic causal models of the brain. Imaging neuroscience (Cambridge, Mass.), 4, IMAG.a.1234. https://doi.org/10.1162/imag.a.1234

BibTeX

@article{zarghami2026multiscale,
author = {Zarghami, Tahereh S},
title = {{Multiscale parcellation of dynamic causal models of the brain}},
journal = {Imaging neuroscience (Cambridge, Mass.)},
year = {2026},
month = may,
volume = {4},
pages = {IMAG.a.1234},
publisher = {MIT Press},
issn = {2837-6056},
doi = {10.1162/imag.a.1234},
url = {https://doi.org/10.1162/imag.a.1234},
pmid = {42206220},
pmcid = {PMC13206500}
}

RIS

TY - JOUR
AU - Zarghami, Tahereh S
TI - Multiscale parcellation of dynamic causal models of the brain
T2 - Imaging neuroscience (Cambridge, Mass.)
J2 - Imaging Neurosci (Camb)
PY - 2026
DA - 2026/05/22
VL - 4
SP - IMAG.a.1234
SN - 2837-6056
PB - MIT Press
DO - 10.1162/imag.a.1234
UR - https://doi.org/10.1162/imag.a.1234
LA - en
ER -

CSL-JSON

{
"id": "10.1162/imag.a.1234",
"type": "article-journal",
"title": "Multiscale parcellation of dynamic causal models of the brain",
"container-title": "Imaging neuroscience (Cambridge, Mass.)",
"author": [
{
"family": "Zarghami",
"given": "Tahereh S"
}
],
"container-title-short": "Imaging Neurosci (Camb)",
"volume": "4",
"page": "IMAG.a.1234",
"DOI": "10.1162/imag.a.1234",
"PMID": "42206220",
"PMCID": "PMC13206500",
"ISSN": "2837-6056",
"publisher": "MIT Press",
"URL": "https://doi.org/10.1162/imag.a.1234",
"language": "en",
"issued": {
"date-parts": [
[
2026,
5,
22
]
]
}
}

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.1038/s41467-026-74466-2 [code]
Neuromorphic hierarchical modular reservoirs.
Journal: Nature communications
In common: computational, 8 references
[2] doi:10.1007/s12021-025-09759-w [code]
Limitations of Variational Laplace-Based Dynamic Causal Modelling for Multistable Cortical Circuits.
Journal: Neuroinformatics
In common: SPM, Statistics and Machine Learning Toolbox, computational, 6 references
[3] doi:10.7554/elife.103097 [code]
Canonical neurodevelopmental trajectories of structural and functional manifolds.
Journal: eLife
In common: Statistics and Machine Learning Toolbox, 6 references
[4] doi:10.1038/s42003-026-10132-z [code]
Hippocampal criticality tracks cognitive demand and shifts with memory impairment.
Journal: Communications biology
In common: 6 references
[5] doi:10.1038/s41467-026-73347-y [code]
A week in the life of the human brain reveals stable states punctuated by chaotic-like transitions.
Journal: Nature communications
In common: Statistics and Machine Learning Toolbox, computational, computational modeling (no new data), 4 references
[6] doi:10.1038/s41593-026-02205-3 [code]
Competitive interactions shape mammalian brain network dynamics and computation.
Journal: Nature neuroscience
In common: SPM, Statistics and Machine Learning Toolbox, computational modeling (no new data), 3 references
[7] doi:10.1162/imag.a.1250 [code]
Dynamics-informed priors (DIP) for neural mass modelling.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: SPM, Statistics and Machine Learning Toolbox, 3 references
[8] doi:10.1002/advs.202523009 [code]
Personalized Network-Guided Neuromodulation Enhances Human Working Memory.
Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany)
In common: SPM, Statistics and Machine Learning Toolbox, 3 references
[9] doi:10.1371/journal.pbio.3003916 [code]
Arousal-driven critical roaming reproduces human functional connectivity dynamics.
Journal: PLoS biology
In common: Statistics and Machine Learning Toolbox, 4 references
[10] doi:10.1038/s43856-026-01707-2 [code]
Decreased amyloid-related structure-function coupling in preclinical Alzheimer's disease.
Journal: Communications medicine
In common: SPM, Statistics and Machine Learning Toolbox, 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.