% Main code to generate functional hierarchies
%
% To obtain a functional hierarchy as in the SciAdv paper, run this code 5
% times with different values of 'indicedata' and average the results using
% the code 'rankfast.m'
%
% Jorge Mejias, 2016
function []=main(indicedata)
%add path for the GC toolbox:
%addpath('/scratch/jmp20/mvgc_v1.0');
%startup;%matlabpool('open', 'local', 10);
%initialize random number generator and load data:
seed1=round(938190+12*indicedata);rng(seed1);
%-----------Data to obtain from core-nets.org:
load('subgraphData30.mat'); %this must include:
%FLN (as flnMat), rank-ordered,
%SLN (as SLNmat), rank-ordered,
%a list of area names rank-ordered (as areaList)
load('subgraphWiring30.mat');
%area-to area distances, rank-ordered and given in mm.
%------------------Now the simulations:
Areas=1:30;
Nareas=length(Areas);
ROI=[1 2 3 4 6 8 9 13]; %Kennedy-Fries Neuron 2015 selection
Nareas2=length(ROI);
%prepare the parameters and compute the main data:
par=parameters(Areas,ROI,flnMat,slnMat,wiring);
[mDAI,X2,X5,DAI,f,F,fshuf,Fshuf]=hierarchy(indicedata,Areas,par,ROI,areaList,flnMat,slnMat);
A=mDAI;
B=F;
C=f;
D=Fshuf;
E=fshuf;
%let's get the gamma and alpha power along the cortex:
gammac=zeros(Nareas,30);alphac=gammac;
for j=1:30
for i=1:Nareas
%gamma power:
[~,~,fg,zg]=analysis(par,X2(i,:,j),30);
%alpha power:
[~,~,fa,za]=analysis(par,X5(i,:,j),3);
gammac(i,j)=zg;alphac(i,j)=za;
end
end
%save the data
output=sprintf('output%d.mat',indicedata);
save(output,'A','B','C','D','E','alphac','gammac','DAI','par');