diff --git a/code_essential/Build_Coop.m b/code_essential/Build_Coop.m
new file mode 100644
index 0000000..134ce93
--- /dev/null
+++ b/code_essential/Build_Coop.m
@@ -0,0 +1,45 @@
+function [M] = Build_Coop(CONST_NUM_NODES,NODE_TYPE,Neuron_type,mult_fact)
+
+% Building the cooperation matrix (first with only 1s on the diagonal)
+M = eye(CONST_NUM_NODES,CONST_NUM_NODES);
+
+% Setting the required elements of the diagonal to 0
+switch Neuron_type
+    case 'Sensory'
+        for n = 1:279
+            if NODE_TYPE(n) ~= 1 && NODE_TYPE(n) ~= 4 && NODE_TYPE(n) ~= 5
+                M(n,n) = mult_fact;
+            end
+        end
+
+    case 'Motor'
+        for n = 1:279
+            if NODE_TYPE(n) ~= 3 && NODE_TYPE(n) ~= 6
+                M(n,n) = mult_fact;
+            end
+        end
+        
+    case 'Inter'
+        for n = 1:279
+            if NODE_TYPE(n) ~= 2
+                M(n,n) = mult_fact;
+            end
+        end
+        
+    case 'Sensory-Inter'
+        for n = 1:279
+            if NODE_TYPE(n) == 3 || NODE_TYPE(n) == 6
+                M(n,n) = mult_fact;
+            end
+        end
+        
+    case 'All'
+        
+    otherwise
+        errordlg('!!!!!');
+end
+
+
+
+
+end
\ No newline at end of file
diff --git a/code_essential/Compute_Matrix.m b/code_essential/Compute_Matrix.m
new file mode 100644
index 0000000..7d278a4
--- /dev/null
+++ b/code_essential/Compute_Matrix.m
@@ -0,0 +1,16 @@
+function [MOI_deg] = Compute_Matrix(A,M,deg)
+
+    MOI_deg = 0;
+
+    % First term
+    for k = 1:deg
+        MOI_deg = MOI_deg + Compute_ck(k) * (M*A^k+A^k*M);
+    end
+    
+    % Second term
+    for k1 = 1:deg-1
+        for k2 = 1:deg-k1    
+            MOI_deg = MOI_deg - (Compute_ck(k1)*Compute_ck(k2)) * A^k1 * M * A^k2;
+        end
+    end
+end
\ No newline at end of file
diff --git a/code_essential/Compute_Slepian_Concentration.m b/code_essential/Compute_Slepian_Concentration.m
new file mode 100644
index 0000000..c66eed9
--- /dev/null
+++ b/code_essential/Compute_Slepian_Concentration.m
@@ -0,0 +1,16 @@
+function [Slep,Mu] = Compute_Slepian_Concentration(U,W,S)
+
+    % Trims the basis
+    U = U(:,1:W);
+    
+    C = U' * S * U;
+    
+    C = (C + C')/2;
+    
+    [V,Mu] = eig(C);
+    
+    [Mu,idx] = sort(diag(Mu),'descend');
+    V = V(:,idx);
+
+    Slep = U*V;
+end
\ No newline at end of file
diff --git a/code_essential/Compute_Slepian_ModifiedEmbedded.m b/code_essential/Compute_Slepian_ModifiedEmbedded.m
new file mode 100644
index 0000000..1a4f6d5
--- /dev/null
+++ b/code_essential/Compute_Slepian_ModifiedEmbedded.m
@@ -0,0 +1,17 @@
+function [Slep,Xi] = Compute_Slepian_ModifiedEmbedded(U,Lambda,W,S)
+
+    % Trims the basis
+    U = U(:,1:W);
+    Lambda = Lambda(1:W,1:W);
+    
+    MC = Lambda^(0.5) * U' * S * U * Lambda^(0.5);
+    
+    MC = (MC+MC')/2;
+    
+    [V,Xi] = eig(MC);
+    
+    [Xi,idx] = sort(diag(Xi),'ascend');
+    V = V(:,idx);
+
+    Slep = U*V;
+end
\ No newline at end of file
diff --git a/code_essential/Compute_Slepianstuff.m b/code_essential/Compute_Slepianstuff.m
new file mode 100644
index 0000000..ce8cc97
--- /dev/null
+++ b/code_essential/Compute_Slepianstuff.m
@@ -0,0 +1,19 @@
+function [V,Zeta,Mu,Xi] = Compute_Slepianstuff(MOI,S,L)
+
+    [V,Zeta] = eig((MOI+MOI')/2);
+    Zeta = diag(Zeta);
+    [Zeta,idx] = sort(Zeta,'descend');
+    V = V(:,idx);
+    Mu = (V'*S*V)./(V'*V);
+    Mu = diag(Mu);
+    Xi = (V'*L^(0.5)*S*L^(0.5)*V)./(V'*V);
+    Xi = diag(Xi);
+
+    for v = 1:size(MOI,1)
+        
+        [~,idx_norm]=max(abs(V(:,v)));
+        if sign(V(idx_norm(1),v))<0
+            V(:,v)=-V(:,v);
+        end
+    end
+end
\ No newline at end of file
diff --git a/code_essential/Compute_ck.m b/code_essential/Compute_ck.m
new file mode 100644
index 0000000..80742a3
--- /dev/null
+++ b/code_essential/Compute_ck.m
@@ -0,0 +1,5 @@
+function [c_k] = Compute_ck(k)
+
+    c_k = (factorial(2*k))/((2^(2*k))*(factorial(k)^2)*(2*k-1));
+
+end
\ No newline at end of file
diff --git a/code_essential/Script_NetworkNeuroscience_Main.m b/code_essential/Script_NetworkNeuroscience_Main.m
new file mode 100644
index 0000000..afdef45
--- /dev/null
+++ b/code_essential/Script_NetworkNeuroscience_Main.m
@@ -0,0 +1,136 @@
+%% 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
\ No newline at end of file
diff --git a/code_essential/Worm279dir.mat b/code_essential/Worm279dir.mat
new file mode 100644
index 0000000..c85498c
Binary files /dev/null and b/code_essential/Worm279dir.mat differ
diff --git a/code_essential/slepNormalize.m b/code_essential/slepNormalize.m
new file mode 100644
index 0000000..d3abec8
--- /dev/null
+++ b/code_essential/slepNormalize.m
@@ -0,0 +1,35 @@
+% Construct normalized (sparse) adjacency matrix
+% 
+% returns
+% A - normalized adjacency matrix
+% D - degree vector
+%
+% Dimitri Van De Ville, Sep 2014
+
+function [An,Dn]=slepNormalize(A,CONST_NORMALIZE)
+
+if ~issparse(A)
+    warning('Adjacency matrix A is not sparse.');
+end;
+
+msize=size(A,1);
+
+% normalize adjacency matrix
+if CONST_NORMALIZE
+    D=A*ones(msize,1);
+    n=D.^(-.5);
+    n=spdiags(n,0,spalloc(msize,msize,msize));
+    An=n*A*n;
+    Dn=spdiags(An*ones(msize,1),0,spalloc(msize,msize,msize));
+else
+    An=A;
+    Dn=spdiags(A*ones(msize,1),0,spalloc(msize,msize,msize));
+end;
+
+% construct D
+%Dn=sum(An);
+%Dn=spdiags(Dn.',0,spalloc(msize,msize,msize));
+
+%L = D - A;
+% determine the normalized Laplacian nL
+%nL = (full(D)^(-.5))*L*(full(D)^(-.5));
\ No newline at end of file