Homec4science
Diffusion Guided Graph Spectral Embedding (master) Browse

guidedGSE/code_essential/Script_NetworkNeuroscience_Main.m
6e883b19735amaster

Last Change
Raw File

Script_NetworkNeuroscience_Main.m
View Options

%% This is a cleaned version of the NetSci 2018 code snippets with the
% purpose to write a journal version of the abstract...
%
% Initially written by Dimitri Van De Ville, and adjusted by Thomas Bolton
% Last modified on August 7th 2018
%% 1. Stuff to specify for specificity of the analysis
% A. What input matrix do we use ?
% Data kindly provided by Dr. Petra Vertes/Bullmore Lab/Cambridge
load('Worm279dir');
A=sparse(Worm279_matrix);
A = full(A);
A = (A+A')/2;
% Number of edges present in the loaded adjacency matrix (set to 2287 if
% full, 1961 is only chemical synapses, 514 if only gap junctions)
n_edges=2287;
% If only using gap junctions
% A=sparse(Worm279_ejunct_matrix_dir);
% n_edges=514;
%
% If only using chemical synapses
% A=sparse((Worm279_synapse_matrix_dir+Worm279_synapse_matrix_dir.')/2);
% n_edges=1961;
% B. What type of neurons do we want to focus on (select between 'Inter',
% 'Motor', 'Sensory', or 'Sensory-Inter'
Neuron_type = 'Sensory';
addpath(genpath(fullfile(pwd,'cbrewer')));
%% 2. Declaration of parameters, loading of variables
% Number of nodes in the C. elegans graph
CONST_NUM_NODES=279;
% Do we want to normalize the adjacency and degree matrices so that D = I ?
% This is done in the current paper version, so set to 1
CONST_NORMALIZE=1;
% Type of operator to use in the computations: either Laplacian matrix (1),
% or modularity matrix (2); currently using the Laplacian matrix, so set to
% 1
CONST_OPERATOR=1;
% Names of the different neuronal subtypes
STR_NEURON_TYPES=...
{ 'ADE|ADF|ADL|AFD|ALM|ALN|AQR|ASE|ASG|ASH|ASI|ASJ|ASK|PVM|PH|AW|AF|OLL|BAG|PLM|CEP|FL|PQR|AVM|URY|HOA|HOB|PDE|PVD|SPV', % sensory
'ADA|AI|AVG|AUA|DVA|DBV|DVC|DVE|DVF|LUA|AVF|PGA|PVT|PVU|PVQ|PVW|PVP|BD|RIC|RIS|RIH|RIP|RIR|RIA|RIB|RIF|PVR|PVS|PVX|RIG|PVC|AVA|AVB|AVD|AVE|AVH|AVJ|AVK|RMG|URB|MCML|MCMR', % inter
'AS01|AS02|AS03|AS04|AS05|AS06|AS07|AS08|AS09|AS10|AS11|MCL|MCR|DA|VB|VA|HS|DD|VC|DB|VD|SM|RMD|RME|RMF|RMH|URA|PDA|PDB', % motor
'IL1|IL2', % polymodal (S, I, M)
'OLQ|PLN|SDQ|URX|ALA', % polymodal (S, I)
'AVL|DVB|MI|PDC|PVN|PVV|PVY|PVZ|RID|RIM|RIV|SAA|SAB|SIA|SIB|CA',% polymodal (I, M)
'PCA|PCB|PCC|SPC|SPD|NSM' }; % polymodal (S, M)
% Types of neurons that we can have
NODE_TYPE_LABEL = { 'sensory','inter','motor', 'sensory-inter-motor', 'sensory-inter', 'inter-motor', 'sensory-motor' };
% Number of such types
NUM_NODE_TYPES=7;
% We scan each entry from Worm279_labels and if there is a match with the
% above table (between neuron name), then we add the neuron type label in
% NODE_TYPE
NODE_TYPE=zeros(CONST_NUM_NODES,1);
for iterN=1:CONST_NUM_NODES
for iterT=1:length(STR_NEURON_TYPES)
if regexp(Worm279_labels{iterN},sprintf('^(%s)\\w*',STR_NEURON_TYPES{iterT}))
if NODE_TYPE(iterN)~=0
fprintf('Double type assignment: node %d (%d-%d)\n',iterN,NODE_TYPE(iterN),iterT);
end
NODE_TYPE(iterN)=iterT;
end
end
end
% Labels, XY coordinates and birth time of the nodes
NODE_LABEL=Worm279_labels;
NODES_XY=Worm279_positions;
NODES_BIRTH = Worm279_birthtime;
% Creates normalized degree and adjacency matrices, so that I - A = L
[A,D]=slepNormalize(sparse(A),CONST_NORMALIZE);
L = eye(CONST_NUM_NODES,CONST_NUM_NODES)-(A+A')/2;
[U_Lap,Lambda_Lap] = eig(L);
[Lambda_Lap,idx_Lap] = sort(diag(Lambda_Lap),'ascend');
U_Lap = U_Lap(:,idx_Lap);
for u = 1:size(U_Lap,1)
[~,idx_norm]=max(abs(U_Lap(:,u)));
if sign(U_Lap(idx_norm(1),u))<0
U_Lap(:,u)=-U_Lap(:,u);
end
end
%% 3. Analysis using guided spectral analysis
% Building of the cooperation matrix
M = Build_Coop(CONST_NUM_NODES,NODE_TYPE,Neuron_type,0);
% Creating selectivity matrix (basically a logical version of M; used only
% in Mu/Xi computations)
S = M;
S(S > 0) = 1;
% Computation of the matrix to perform eigendecomposition on (full,
% linear/qiadratic approximations)
MOI = M - (L^(0.5))*M*(L^(0.5));
MOI_lin = (M*A+A*M)/2;
MOI_quad = (M*A+A*M)/2 + (M*A*A+A*A*M)/8 - (A*M*A)/4;
% Computation of the same matrix using a higher-order approximation
max_deg = 5;
for d = 2:max_deg
disp(d);
MOI_deg(:,:,d-1) = Compute_Matrix(A,M,d);
[V_deg(:,:,d-1),Zeta_deg(:,d-1),Mu_deg(:,d-1),Xi_deg(:,d-1)] = Compute_Slepianstuff(MOI_deg(:,:,d-1),S,L);
end
% Extraction of eigenvalues and eigenvectors; the called function also
% reorders them, and flips sign to have the largest magnitude element
% positive
[V,Zeta,Mu,Xi] = Compute_Slepianstuff(MOI,S,L);
[V_lin,Zeta_lin,Mu_lin,Xi_lin] = Compute_Slepianstuff(MOI_lin,S,L);
[V_quad,Zeta_quad,Mu_quad,Xi_quad] = Compute_Slepianstuff(MOI_quad,S,L);
% Computation of old Slepian criteria at different bandwidths
W = [100 150 200 279];
for w = 1:length(W)
[Slep_con{w}(:,:),Slep_Mu{w}] = Compute_Slepian_Concentration(U_Lap,W(w),S);
[Slep_mod{w}(:,:),Slep_Xi{w}] = Compute_Slepian_ModifiedEmbedded(U_Lap,diag(Lambda_Lap),W(w),S);
end
c4science · Help