a function for calculating the small world propensity of a given network - assumes that matrix is undirected (symmeteric) and if not, creates a symmetric matrix which is used for the calculations
0001 function [SWP,delta_C,delta_L] = smallWorldPropensity(A, varargin) 0002 0003 % a function for calculating the small world propensity of 0004 % a given network - assumes that matrix is undirected (symmeteric) and if 0005 % not, creates a symmetric matrix which is used for the calculations 0006 0007 %NOTE: This code requires the Bioinformatics Toolbox to be installed 0008 % (uses graphallshortestpaths.m) 0009 0010 %Inputs: 0011 % A the connectivity matrix, weighted or binary 0012 % varargin a string corresponding to the method of clustering 0013 % to be used, where 'O' is Onnela, 'Z' is Zhang, 0014 % 'B' is Barrat, 'bin' is binary (default is Onnela). 0015 % If the user specifies binary analysis, a 0016 % weighted matrix will be converted to a binary matrix 0017 % before proceeding. 0018 % 0019 0020 %Outputs: 0021 % SWP the small world propensity of the matrix 0022 % delta_C the fractional deviation from the expected culstering coefficient of a 0023 % random network 0024 % delta_L the fractional deviation from the expected path length of a 0025 % random network 0026 0027 %written by Eric Bridgeford and modified by Sarah F. Muldoon 0028 0029 % Reference: Muldoon, Bridgeford, and Bassett (2015) "Small-World Propensity in Weighted, 0030 % Real-World Networks" http://arxiv.org/abs/1505.02194 0031 % Code from http://commdetect.weebly.com/ 0032 % Danielle S. Bassett - UPenn 0033 % Some edits by Byungchan Kim - Washington University 0034 0035 if isempty(varargin) 0036 varargin{1} = 'O'; 0037 end 0038 0039 if sum(sum(A)) > 0 0040 0041 bin_matrix = 0; 0042 if strcmp(varargin{1},'bin') == 1 0043 bin_matrix = 1; 0044 A = A > 0; 0045 end 0046 0047 %check to see if matrix is symmeteric 0048 symcheck=abs(A-A'); 0049 if sum(sum(symcheck)) > 0 0050 % adjust the input matrix to symmeterize 0051 disp('Input matrix is not symmetric. Symmetrizing.') 0052 W = symm_matrix(A, bin_matrix); 0053 else 0054 W=A; 0055 end 0056 0057 %calculate the number of nodes 0058 n = length(W); 0059 %compute the weighted density of the network 0060 dens_net = sum(sum(W))/(max(max(W))*n*(n-1)); 0061 0062 %compute the average degree of the unweighted network, to give 0063 %the approximate radius 0064 numb_connections = length(find(W>0)); 0065 avg_deg_unw = numb_connections/n; 0066 avg_rad_unw = avg_deg_unw/2; 0067 avg_rad_eff = ceil(avg_rad_unw); 0068 0069 0070 %compute the regular and random matrix for the network W 0071 W_reg = regular_matrix_generator(W, avg_rad_eff); 0072 W_rand = randomize_matrix(W); 0073 0074 %compute all path length calculations for the network 0075 reg_path = avg_path_matrix(1./W_reg); %path of the regular network 0076 rand_path = avg_path_matrix(1./W_rand); %path of the random netowork 0077 net_path = avg_path_matrix(1./W); %path of the network 0078 0079 A = (net_path - rand_path); 0080 if A < 0 0081 A = 0; 0082 end 0083 diff_path = A/ (reg_path - rand_path); 0084 if net_path == Inf || rand_path == Inf || reg_path == Inf 0085 diff_path = 1; 0086 end 0087 if diff_path > 1 0088 diff_path = 1; 0089 end 0090 0091 0092 %compute all clustering calculations for the network 0093 reg_clus = avg_clus_matrix(W_reg,varargin{1}); 0094 rand_clus = avg_clus_matrix(W_rand,varargin{1}); 0095 net_clus = avg_clus_matrix(W,varargin{1}); 0096 0097 B = (reg_clus - net_clus); 0098 if B < 0 0099 B = 0; 0100 end 0101 0102 diff_clus = B / (reg_clus - rand_clus); 0103 if isnan(reg_clus) || isnan(rand_clus) || isnan(net_clus) 0104 diff_clus = 1; 0105 end 0106 if diff_clus > 1 0107 diff_clus = 1; 0108 end 0109 0110 %calculate small world value, the root sum of the squares of 0111 %diff path and diff clus 0112 SWP = 1 - (sqrt(diff_clus^2 + diff_path^2)/sqrt(2)); 0113 delta_C=diff_clus; 0114 delta_L=diff_path; 0115 end 0116 end 0117 0118 0119 %XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 0120 %XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 0121 %XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 0122 %XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 0123 %the above code requires the following funcitons 0124 0125 function [Clus] = avg_clus_matrix(W, met) 0126 %a function to compute the average clusteirng coefficient for a 0127 %input matrix M 0128 0129 %Inputs: 0130 % W a matrix, weighted or unweighted 0131 % met a string, to represent the method to be used for computing 0132 % the clustering coefficient 0133 % possible strings: 'O' (Onnela), 'Z' (Zhang), 'B' (Barrat), 0134 % 'bin' (binary) 0135 % default if none is chosen is Onnela 0136 0137 %Outputs: 0138 % Clus the average clustering coefficient 0139 0140 %written by Eric Bridgeford 0141 0142 n = length(W); 0143 [C] = clustering_coef_matrix(W, met); 0144 Clus = nanmean(C); 0145 end 0146 0147 0148 function [Len] = avg_path_matrix(M) 0149 0150 %a function to compute the average path length of a given matrix 0151 %using the graphallshortestpaths built-in matlab function 0152 0153 %written by Eric Bridgeford 0154 0155 n = length(M); 0156 M = sparse(M); 0157 D = graphallshortestpaths(M); 0158 0159 %checks if a node is disconnected from the system, and replaces 0160 %its value with 0 0161 for i = 1:n 0162 for j = 1:n 0163 if isinf(D(i,j)) == 1 0164 D(i,j) = 0; 0165 end 0166 end 0167 end 0168 0169 Len = mean(mean(D)); 0170 end 0171 0172 0173 function [C] = clustering_coef_matrix(W, met) 0174 0175 %a modification of the clustering coefficient function provided 0176 %in the brain connectivity toolbox 0177 0178 %improved definition of Onnela Clustering Coefficient, as well as 0179 %implementation of function for Zhang and Barrat clustering values 0180 0181 %Reference: 0182 % Onnela et al., Phys. Rev. E71, 065103(R)(2005) 0183 % B.Zhang and S. Horvath, Stat. App. Genet. Mol. Biol.4, 17(2005) 0184 % Barrat et al., Proc. Natl. Acad. Sci. U.S.A.101, 3747(2004) 0185 % Watts and Strogatz (1998) Nature 393:440-442 0186 0187 %Inputs: 0188 % W the weighted or unweighted connectivity matrix 0189 % met a string, to represent the method to be used for computing 0190 % the clustering coefficient 0191 % possible strings: 'O' (Onnela), 'Z' (Zhang), 'B' (Barrat), 'bin' 0192 % (binary) 0193 % default if none is chosen is Onnela 0194 0195 0196 %code originally written by Mika Rubinov, UNSW, 2007-2010 0197 %modified/written by Eric Bridgeford 0198 0199 if met == 'O' 0200 K=sum(W~=0,2); 0201 W = double(W); 0202 W2 = W/max(max(W)); 0203 cyc3=diag(W2.^(1/3)^3); 0204 K(cyc3==0)=inf; %if no 3-cycles exist, make C=0 (via K=inf) 0205 C=cyc3./(K.*(K-1)); 0206 end 0207 if met == 'bin' 0208 G = double(W>0); 0209 n=length(G); 0210 C=zeros(n,1); 0211 for u=1:n 0212 V=find(G(u,:)); 0213 k=length(V); 0214 if k>=2; %degree must be at least 2 0215 S=G(V,V); 0216 C(u)=sum(S(:))/(k^2-k); 0217 end 0218 end 0219 end 0220 0221 if met == 'Z' 0222 K=sum(W~=0,2); 0223 W = double(W); 0224 W2 = W/max(max((W))); 0225 cyc3=diag((W2)^3); 0226 denom = zeros(length(W),1); 0227 for i = 1:length(W) 0228 denom(i) = (sum(W2(i,:))^2-sum(W2(i,:).^2)); 0229 end 0230 C = cyc3./denom; 0231 end 0232 0233 if met == 'B' 0234 A = double(W>0); 0235 C = zeros(length(W),1); 0236 for i = 1:length(W) 0237 sum1 = 0; 0238 for j = 1:length(W) 0239 for k = 1:length(W) 0240 sum1 = ((W(i,j)+W(i,k))/2)*A(i,j)*A(j,k)*A(i,k)+sum1; 0241 end 0242 end 0243 C(i) = 1/(sum(W(i,:))*(sum(A(i,:))-1))*sum1; 0244 end 0245 end 0246 0247 0248 end 0249 0250 0251 function A_rand=randomize_matrix(A); 0252 0253 %This code creates a random undirected network from the connectivity 0254 %distribution of an undirected adjacency matrix, ie, the intital matrix 0255 %must be symmetric. 0256 0257 % INPUTS: 0258 % A: an undirected adjacency matrix (symmetric) with no self connections 0259 0260 % OUTPUTS: 0261 % A_rand: a comparable random network with same number of nodes and 0262 % connectivity distribution 0263 0264 % written by Sarah F. Muldoon 0265 0266 num_nodes=length(A); 0267 A_rand=zeros(num_nodes); 0268 mask=triu(ones(num_nodes),1); 0269 grab_indices=find(mask > 0); 0270 0271 orig_edges=A(grab_indices); 0272 num_edges=length(orig_edges); 0273 0274 rand_index=randperm(num_edges); 0275 randomized_edges=orig_edges(rand_index); 0276 0277 edge=1; 0278 for i=1:num_nodes-1 0279 for j=i+1:num_nodes 0280 A_rand(i,j)=randomized_edges(edge); 0281 A_rand(j,i)=randomized_edges(edge); 0282 edge=edge+1; 0283 end 0284 end 0285 end 0286 0287 0288 function M = regular_matrix_generator(G,r) 0289 %generates a regular matrix, with weights obtained form the 0290 %original adjacency matrix representation of the network 0291 0292 % note that all inputs should be symmeterized prior to forming a regular 0293 % matrix, since otherwise half of the connnections will be trashed. This 0294 % can be accomplished with the built in symm_matrix function, however, 0295 % the function is not performed here so that users can use their own 0296 % symmeterization procedure. 0297 0298 %Inputs: 0299 % G the adjacency matrix for the given network; must be symmmeterized 0300 % r the approximate radius of the regular network 0301 0302 %Outputs: 0303 % M the regular matrix for the given network, where all 0304 % weights are sorted such that the inner radius has the 0305 % highest weights randomly distributed across the nodes, 0306 % and so on 0307 0308 %written by Eric W. Bridgeford 0309 0310 n = length(G); 0311 G = triu(G); 0312 %reshape the matrix G into an array, B 0313 B = reshape(G,[length(G)^2,1]); 0314 %sorts the array in descending order 0315 B = sort(B,'descend'); 0316 %computes the number of connections and adds zeros if 0317 %numel(G) < 2*n*r 0318 num_els =ceil(numel(G)/(2*n)); 0319 num_zeros = 2*n*num_els - numel(G); 0320 %adds zeros to the remaineder of the list, so length(B) = 2*n*r 0321 B = cat(1,B,zeros(num_zeros,1)); 0322 %reshapes B into a matrix, where the values descend top to 0323 %bottom, as well as left to right. The greatest value in each 0324 %column is less than the smallest value of the column to its left. 0325 B = reshape(B,[n],[]); 0326 0327 M = zeros(length(G)); 0328 0329 %distributes the connections into a regular network, M, where 0330 %the innermost radius represents the highest values 0331 for i = 1:length(G) 0332 for z = 1:r 0333 a = randi([1,n]); 0334 0335 %random integer chosen to take a value from B 0336 while (B(a,z) == 0 && z ~= r) || (B(a,z) == 0 && z == r && ~isempty(find(B(:,r),1))) 0337 a = randi([1,n]); 0338 end 0339 %finds the two nodes a distance of z from the origin node 0340 %and places the entries from the matrix B 0341 y_coor_1 = mod(i+z-1,length(G))+1; 0342 %y_coor_2 = mod(i-z-1,length(G))+1; 0343 M(i,y_coor_1) = B(a,z); 0344 M(y_coor_1,i) = B(a,z); 0345 %removes the weights from the matrix B so they cannot be 0346 %reused 0347 B(a,z) = 0; 0348 0349 end 0350 end 0351 0352 end 0353 0354 0355 function [W] = symm_matrix(A, bin_key) 0356 0357 % a function to symmetrize an input matrix. The procedure 0358 % by which items are symmetrized such that: 0359 % in the binary case: 0360 % if a(i,j) || a(j,i) == 1 for i,j in A 0361 % w(i,j) && w(j,i) == 1 0362 % in the weighted case: 0363 % if (a(i,j) || a(j,i) > 0 for i,j in A 0364 % w(i,j) && w(j,i) == (a(i,j) + a(j,i) )/ 2 0365 0366 % Inputs: 0367 % A: The binary or weighted input matrix 0368 % bin_key: the key to indicate whether weighted or binary analysis 0369 % will take place 0370 % 1 indicates binarized, 0 indicates weighted 0371 0372 % Outputs 0373 % W: The symmeterized matrix 0374 0375 % if binary analysis is specified, let binary symmeterization take place 0376 0377 % written by Eric W. Bridgeford 0378 0379 W = zeros(length(A)); 0380 0381 if bin_key == 1 0382 A = A > 0; % verify that the input matrix is binary 0383 for i = 1:length(A) 0384 for j = i:length(A) 0385 if A(i,j) || A(j,i) 0386 W(i,j) = 1; 0387 W(j,i) = 1; 0388 end 0389 end 0390 end 0391 else 0392 for i = 1:length(A) 0393 for j = i:length(A) 0394 if A(i,j) || A(j,i) 0395 val = (A(i,j) + A(j,i)) / 2; 0396 W(i,j) = val; 0397 W(j,i) = val; 0398 end 0399 end 0400 end 0401 end 0402 end 0403 0404 0405 0406 0407 0408