Home‎ > ‎

 adjacency_plot_und ``` ADJACENCY_PLOT_UND Quick visualization tool [X,Y,Z] = ADJACENCY_PLOT(AIJ,COOR) takes adjacency matrix AIJ and node spatial coordinates COOR and generates three vectors that can be used for quickly plotting the edges in AIJ. If no coordinates are specified, then each node is assigned a position on a circle. COOR can, in general, be 2D or 3D. Example: >> load AIJ; % load your adjacency matrix >> load COOR; % load 3D coordinates for each node >> [x,y,z] = adjacency_plot_und(AIJ,COOR); % call function >> plot3(x,y,z); % plots network as a single line object If COOR were 2D, the PLOT3 command changes to a PLOT command. NOTE: This function is similar to MATLAB's GPLOT command. Richard Betzel, Indiana University, 2013 ``` agreement ``` AGREEMENT Agreement matrix from clusters D = AGREEMENT(CI) takes as input a set of vertex partitions CI of dimensions [vertex x partition]. Each column in CI contains the assignments of each vertex to a class/community/module. This function aggregates the partitions in CI into a square [vertex x vertex] agreement matrix D, whose elements indicate the number of times any two vertices were assigned to the same class. In the case that the number of nodes and partitions in CI is large (greater than ~1000 nodes or greater than ~1000 partitions), the script can be made faster by computing D in pieces. The optional input BUFFSZ determines the size of each piece. Trial and error has found that BUFFSZ ~ 150 works well. Inputs, CI, set of (possibly) degenerate partitions BUFFSZ, optional second argument to set buffer size Outputs: D, agreement matrix Richard Betzel, Indiana University, 2012 ``` agreement_weighted ``` AGREEMENT_WEIGHTED Weights agreement matrix D = AGREEMENT_WEIGHTED(CI,WTS) is identical to AGREEMENT, with the exception that each partitions contribution is weighted according to the corresponding scalar value stored in the vector WTS. As an example, suppose CI contained partitions obtained using some heuristic for maximizing modularity. A possible choice for WTS might be the Q metric (Newman's modularity score). Such a choice would add more weight to higher modularity partitions. NOTE: Unlike AGREEMENT, this script does not have the input argument BUFFSZ. Inputs: CI, set of partitions WTS, relative weight of importance of each paritition Outputs: D, weighted agreement matrix Richard Betzel, Indiana University, 2013 ``` align_matrices ``` ALIGN_MATRICES Alignment of two matrices [Mreordered,Mindices,cost] = align_matrices(M1,M2,dfun,flag) This function aligns two matrices relative to one another by reordering the nodes in M2. The function uses a version of simulated annealing. Inputs: M1 = first connection matrix (square) M2 = second connection matrix (square) dfun = distance metric to use for matching: 'absdff' = absolute difference 'sqrdff' = squared difference 'cosang' = cosine of vector angle Mreordered = reordered connection matrix M2 Mindices = reordered indices cost = distance between M1 and Mreordered Connection matrices can be weighted or binary, directed or undirected. They must have the same number of nodes. M1 can be entered in any node ordering. Note that in general, the outcome will depend on the initial condition (the setting of the random number seed). Also, there is no good way to determine optimal annealing parameters in advance - these parameters will need to be adjusted "by hand" (particularly H, Texp, T0, and Hbrk). For large and/or dense matrices, it is highly recommended to perform exploratory runs varying the settings of 'H' and 'Texp' and then select the best values. Based on extensive testing, it appears that T0 and Hbrk can remain unchanged in most cases. Texp may be varied from 1-1/H to 1-10/H, for example. H is the most important parameter - set to larger values as the problem size increases. Good solutions can be obtained for matrices up to about 100 nodes. It is advisable to run this function multiple times and select the solution(s) with the lowest 'cost'. If the two matrices are related it may be very helpful to pre-align them by reordering along their largest eigenvectors: [v,~] = eig(M1); v1 = abs(v(:,end)); [a1,b1] = sort(v1); [v,~] = eig(M2); v2 = abs(v(:,end)); [a2,b2] = sort(v2); [a,b,c] = overlapMAT2(M1(b1,b1),M2(b2,b2),'dfun',1); Setting 'Texp' to zero cancels annealing and uses a greedy algorithm instead. Yusuke Adachi, University of Tokyo, 2010 Olaf Sporns, Indiana University, 2010 ``` assortativity_bin ``` ASSORTATIVITY_BIN Assortativity coefficient r = assortativity(CIJ,flag); The assortativity coefficient is a correlation coefficient between the degrees of all nodes on two opposite ends of a link. A positive assortativity coefficient indicates that nodes tend to link to other nodes with the same or similar degree. Inputs: CIJ, binary directed/undirected connection matrix flag, 0, undirected graph: degree/degree correlation 1, directed graph: out-degree/in-degree correlation 2, directed graph: in-degree/out-degree correlation 3, directed graph: out-degree/out-degree correlation 4, directed graph: in-degree/in-degree correlation Outputs: r, assortativity coefficient Notes: The function accepts weighted networks, but all connection weights are ignored. The main diagonal should be empty. For flag 1 the function computes the directed assortativity described in Rubinov and Sporns (2010) NeuroImage. Reference: Newman (2002) Phys Rev Lett 89:208701 Foster et al. (2010) PNAS 107:10815�10820 Olaf Sporns, Indiana University, 2007/2008 Vassilis Tsiaras, University of Crete, 2009 Murray Shanahan, Imperial College London, 2012 Mika Rubinov, University of Cambridge, 2012 ``` assortativity_wei ``` ASSORTATIVITY_WEI Assortativity coefficient r = assortativity_wei(CIJ,flag); The assortativity coefficient is a correlation coefficient between the strengths (weighted degrees) of all nodes on two opposite ends of a link. A positive assortativity coefficient indicates that nodes tend to link to other nodes with the same or similar strength. Inputs: CIJ, weighted directed/undirected connection matrix flag, 0, undirected graph: strength/strength correlation 1, directed graph: out-strength/in-strength correlation 2, directed graph: in-strength/out-strength correlation 3, directed graph: out-strength/out-strength correlation 4, directed graph: in-strength/in-strength correlation Outputs: r, assortativity coefficient Notes: The main diagonal should be empty. For flag 1 the function computes the directed assortativity described in Rubinov and Sporns (2010) NeuroImage. Reference: Newman (2002) Phys Rev Lett 89:208701 Foster et al. (2010) PNAS 107:10815-10820 Olaf Sporns, Indiana University, 2007/2008 Vassilis Tsiaras, University of Crete, 2009 Murray Shanahan, Imperial College London, 2012 Mika Rubinov, University of Cambridge, 2012 ``` backbone_wu ``` BACKBONE_WU Backbone [CIJtree,CIJclus] = backbone_wu(CIJ,avgdeg) The network backbone contains the dominant connections in the network and may be used to aid network visualization. This function computes the backbone of a given weighted and undirected connection matrix CIJ, using a minimum-spanning-tree based algorithm. input: CIJ, connection/adjacency matrix (weighted, undirected) avgdeg, desired average degree of backbone output: CIJtree, connection matrix of the minimum spanning tree of CIJ CIJclus, connection matrix of the minimum spanning tree plus strongest connections up to an average degree 'avgdeg' NOTE: nodes with zero strength are discarded. NOTE: CIJclus will have a total average degree exactly equal to (or very close to) 'avgdeg'. NOTE: 'avgdeg' backfill is handled slightly differently than in Hagmann et al 2008. Reference: Hidalgo et al. (2007) Science 317, 482. Hagmann et al. (2008) PLoS Biol Olaf Sporns, Indiana University, 2007/2008/2010/2012 ``` betweenness_bin ``` BETWEENNESS_BIN Node betweenness centrality BC = betweenness_bin(A); Node betweenness centrality is the fraction of all shortest paths in the network that contain a given node. Nodes with high values of betweenness centrality participate in a large number of shortest paths. Input: A, binary (directed/undirected) connection matrix. Output: BC, node betweenness centrality vector. Note: Betweenness centrality may be normalised to the range [0,1] as BC/[(N-1)(N-2)], where N is the number of nodes in the network. Reference: Kintali (2008) arXiv:0809.1906v2 [cs.DS] (generalization to directed and disconnected graphs) Mika Rubinov, UNSW/U Cambridge, 2007-2012 ``` betweenness_wei ``` BETWEENNESS_WEI Node betweenness centrality BC = betweenness_wei(L); Node betweenness centrality is the fraction of all shortest paths in the network that contain a given node. Nodes with high values of betweenness centrality participate in a large number of shortest paths. Input: L, Directed/undirected connection-length matrix. Output: BC, node betweenness centrality vector. Notes: The input matrix must be a connection-length matrix, typically obtained via a mapping from weight to length. For instance, in a weighted correlation network higher correlations are more naturally interpreted as shorter distances and the input matrix should consequently be some inverse of the connectivity matrix. Betweenness centrality may be normalised to the range [0,1] as BC/[(N-1)(N-2)], where N is the number of nodes in the network. Reference: Brandes (2001) J Math Sociol 25:163-177. Mika Rubinov, UNSW/U Cambridge, 2007-2012 ``` breadth ``` BREADTH Auxiliary function for breadthdist.m [distance,branch] = breadth(CIJ,source); Implementation of breadth-first search. Input: CIJ, binary (directed/undirected) connection matrix source, source vertex Outputs: distance, distance between 'source' and i'th vertex (0 for source vertex) branch, vertex that precedes i in the breadth-first search tree (-1 for source vertex) Notes: Breadth-first search tree does not contain all paths (or all shortest paths), but allows the determination of at least one path with minimum distance. The entire graph is explored, starting from source vertex 'source'. Olaf Sporns, Indiana University, 2002/2007/2008 ``` breadthdist ``` BREADTHDIST Reachability and distance matrices [R,D] = breadthdist(CIJ); The binary reachability matrix describes reachability between all pairs of nodes. An entry (u,v)=1 means that there exists a path from node u to node v; alternatively (u,v)=0. The distance matrix contains lengths of shortest paths between all pairs of nodes. An entry (u,v) represents the length of shortest path from node u to node v. The average shortest path length is the characteristic path length of the network. Input: CIJ, binary (directed/undirected) connection matrix Outputs: R, reachability matrix D, distance matrix Note: slower but less memory intensive than "reachdist.m". Algorithm: Breadth-first search. Olaf Sporns, Indiana University, 2002/2007/2008 ``` charpath ``` CHARPATH Characteristic path length, global efficiency and related statistics lambda = charpath(D); lambda = charpath(D); [lambda,efficiency] = charpath(D); [lambda,efficiency,ecc,radius,diameter] = charpath(D,diagonal_dist,infinite_dist); The network characteristic path length is the average shortest path length between all pairs of nodes in the network. The global efficiency is the average inverse shortest path length in the network. The nodal eccentricity is the maximal path length between a node and any other node in the network. The radius is the minimal eccentricity, and the diameter is the maximal eccentricity. Input: D, distance matrix diagonal_dist optional argument include distances on the main diagonal (default: diagonal_dist=0) infinite_dist optional argument include infinite distances in calculation (default: infinite_dist=1) Outputs: lambda, network characteristic path length efficiency, network global efficiency ecc, nodal eccentricity radius, network radius diameter, network diameter Notes: The input distance matrix may be obtained with any of the distance functions, e.g. distance_bin, distance_wei. Characteristic path length is defined here as the mean shortest path length between all pairs of nodes, for consistency with common usage. Note that characteristic path length is also defined as the median of the mean shortest path length from each node to all other nodes. Infinitely long paths (i.e. paths between disconnected nodes) are included in computations by default. This behavior may be modified with via the infinite_dist argument. Olaf Sporns, Indiana University, 2002/2007/2008 Mika Rubinov, U Cambridge, 2010/2015 ``` clique_communities ``` CLIQUE_COMMUNITIES Overlapping community structure via clique percolation M = clique_communities(A, cq_thr) The optimal community structure is a subdivision of the network into groups of nodes which have a high number of within-group connections and a low number of between group connections. This algorithm uncovers overlapping community structure in binary undirected networks via the clique percolation method. Inputs: A, Binary undirected connection matrix. cq_thr, Clique size threshold (integer). Larger clique size thresholds potentially result in larger communities. Output: M, Overlapping community-affiliation matrix Binary matrix of size CxN [communities x nodes] Algorithms: Bron–Kerbosch algorithm for detection of maximal cliques. Dulmage-Mendelsohn decomposition for detection of components (implemented in get_components.m) Note: This algorithm can be slow and memory intensive in large matrices. The algorithm requires the function get_components.m Reference: Palla et al. (2005) Nature 435, 814-818. Mika Rubinov, Janelia HHMI, 2017 ``` clustering_coef_bd ``` CLUSTERING_COEF_BD Clustering coefficient C = clustering_coef_bd(A); The clustering coefficient is the fraction of triangles around a node (equiv. the fraction of node's neighbors that are neighbors of each other). Input: A, binary directed connection matrix Output: C, clustering coefficient vector Reference: Fagiolo (2007) Phys Rev E 76:026107. Mika Rubinov, UNSW, 2007-2010 ``` clustering_coef_bu ``` CLUSTERING_COEF_BU Clustering coefficient C = clustering_coef_bu(A); The clustering coefficient is the fraction of triangles around a node (equiv. the fraction of node's neighbors that are neighbors of each other). Input: A, binary undirected connection matrix Output: C, clustering coefficient vector Reference: Watts and Strogatz (1998) Nature 393:440-442. Mika Rubinov, UNSW, 2007-2010 ``` clustering_coef_wd ``` CLUSTERING_COEF_WD Clustering coefficient C = clustering_coef_wd(W); The weighted clustering coefficient is the average "intensity" (geometric mean) of all triangles associated with each node. Input: W, weighted directed connection matrix (all weights must be between 0 and 1) Output: C, clustering coefficient vector Reference: Fagiolo (2007) Phys Rev E 76:026107. Note: All weights must be between 0 and 1. This may be achieved using the weight_conversion.m function, W_nrm = weight_conversion(W, 'normalize'); Mika Rubinov, UNSW/U Cambridge, 2007-2015 ``` clustering_coef_wu ``` CLUSTERING_COEF_WU Clustering coefficient C = clustering_coef_wu(W); The weighted clustering coefficient is the average "intensity" (geometric mean) of all triangles associated with each node. Input: W, weighted undirected connection matrix (all weights must be between 0 and 1) Output: C, clustering coefficient vector Note: All weights must be between 0 and 1. This may be achieved using the weight_conversion.m function, W_nrm = weight_conversion(W, 'normalize'); Reference: Onnela et al. (2005) Phys Rev E 71:065103 Mika Rubinov, UNSW/U Cambridge, 2007-2015 ``` clustering_coef_wu_sign ``` CLUSTERING_COEF_WU_SIGN Multiple generalizations of the clustering coefficient [C_pos,C_neg,Ctot_pos,Ctot_neg] = clustering_coef_wu_sign(W,coef_type); The weighted clustering coefficient is the average weight or intensity of all triangles associated with each node. Inputs: W, Weighted undirected connection matrix corr_type, Desired type of clustering coefficient. Options: 1, (default) Onnela et al. formula, used in original clustering_coef_wu.m. Computed separately for positive & negative weights. 2, Zhang & Horvath formula, similar to Onnela formula except denominator of Onnela formula relies on binarizing the network whereas this denominator is based on weight value, which reduces the sensitivity of this measure to the weights directly connected to the node of interest. Computed separately for positive & negative weights. 3, Constantini & Perugini's generalization of the Zhang & Horvath formula. This formula takes both positive & negative weights into account simultaneously, & is particularly sensitive to non-redundancy in path information based on sign (i.e., when two weights are positive & one negative, or all three are negative, both of which indicate that the weight of the third path is not redundant information). Produces only one value. Outputs: C_pos/C_neg, Clustering coefficient vector for positive/negative weights. For the third option, only one vector is outputted (as C_pos). Ctot_pos/Ctot_neg, Mean clustering coefficient for positive and negative weights. References: Onnela et al. (2005) Phys Rev E 71:065103 Zhang & Horvath (2005) Stat Appl Genet Mol Biol 41:1544-6115 Costantini & Perugini (2014) PLOS ONE 9:e88669 Contributor: Jeff Spielberg, Boston University, 2014-2015 (script based on clustering_coef_wu.m) ``` community_louvain ``` COMMUNITY_LOUVAIN Optimal community structure M = community_louvain(W); [M,Q] = community_louvain(W,gamma); [M,Q] = community_louvain(W,gamma,M0); [M,Q] = community_louvain(W,gamma,M0,'potts'); [M,Q] = community_louvain(W,gamma,M0,'negative_asym'); [M,Q] = community_louvain(W,[],[],B); The optimal community structure is a subdivision of the network into nonoverlapping groups of nodes which maximizes the number of within- group edges, and minimizes the number of between-group edges. This function is a fast and accurate multi-iterative generalization of the Louvain community detection algorithm. This function subsumes and improves upon, modularity_louvain_und.m, modularity_finetune_und.m, modularity_louvain_dir.m, modularity_finetune_dir.m, modularity_louvain_und_sign.m and additionally allows to optimize other objective functions (includes built-in Potts-model Hamiltonian, allows for custom objective-function matrices). Inputs: W, directed/undirected weighted/binary connection matrix with positive and possibly negative weights. gamma, resolution parameter (optional) gamma>1, detects smaller modules 0<=gamma<1, detects larger modules gamma=1, classic modularity (default) M0, initial community affiliation vector (optional) B, objective-function type or custom objective matrix (optional) 'modularity', modularity (default) 'potts', Potts-model Hamiltonian (for binary networks) 'negative_sym', symmetric treatment of negative weights 'negative_asym', asymmetric treatment of negative weights B, custom objective-function matrix Note: see Rubinov and Sporns (2011) for a discussion of symmetric vs. asymmetric treatment of negative weights. Outputs: M, community affiliation vector Q, optimized community-structure statistic (modularity by default) Example: % Iterative community finetuning. % W is the input connection matrix. n = size(W,1); % number of nodes M = 1:n; % initial community affiliations Q0 = -1; Q1 = 0; % initialize modularity values while Q1-Q0>1e-5; % while modularity increases Q0 = Q1; % perform community detection [M, Q1] = community_louvain(W, [], M); end References: Blondel et al. (2008) J. Stat. Mech. P10008. Reichardt and Bornholdt (2006) Phys. Rev. E 74, 016110. Ronhovde and Nussinov (2008) Phys. Rev. E 80, 016109 Sun et al. (2008) Europhysics Lett 86, 28004. Rubinov and Sporns (2011) Neuroimage 56:2068-79. Mika Rubinov, U Cambridge 2015-2016 ``` consensus_und ``` CONSENSUS_UND Consensus clustering CIU = CONSENSUS(D,TAU,REPS) seeks a consensus partition of the agreement matrix D. The algorithm used here is almost identical to the one introduced in Lancichinetti & Fortunato (2012): The agreement matrix D is thresholded at a level TAU to remove an weak elements. The resulting matrix is then partitions REPS number of times using the Louvain algorithm (in principle, any clustering algorithm that can handle weighted matrixes is a suitable alternative to the Louvain algorithm and can be substituted in its place). This clustering produces a set of partitions from which a new agreement is built. If the partitions have not converged to a single representative partition, the above process repeats itself, starting with the newly built agreement matrix. NOTE: In this implementation, the elements of the agreement matrix must be converted into probabilities. NOTE: This implementation is slightly different from the original algorithm proposed by Lanchichinetti & Fortunato. In its original version, if the thresholding produces singleton communities, those nodes are reconnected to the network. Here, we leave any singleton communities disconnected. Inputs: D, agreement matrix with entries between 0 and 1 denoting the probability of finding node i in the same cluster as node j TAU, threshold which controls the resolution of the reclustering REPS, number of times that the clustering algorithm is reapplied Outputs: CIU, consensus partition References: Lancichinetti & Fortunato (2012). Consensus clustering in complex networks. Scientific Reports 2, Article number: 336. Richard Betzel, Indiana University, 2012 modified on 3/2014 to include "unique_partitions" ``` core_periphery_dir ``` CORE_PERIPHERY_DIR Core/periphery structure and core-ness statistic C = core_periphery_dir(W) [C,q] = core_periphery_dir(W,gamm,C0) The optimal core/periphery subdivision is a partition of the network into two non-overlapping groups of nodes, a core group and a periphery group, in a way that maximizes the number/weight of within core-group edges, and minimizes the number/weight of within periphery-group edges. The core-ness is a statistic which quantifies the goodness of the optimal core/periphery subdivision. Input: W directed (weighted or binary) connection matrix. gamma, core-ness resolution parameter (optional) gamma>1 detects small core/large periphery 0<=gamma<1 detects large core/small periphery default is gamma=1 Outputs: C, binary vector of optimal core structure C = 1 represents nodes in the core C = 0 represents nodes in the periphery q, maximized core-ness statistic Algorithm: A version of Kernighan-Lin algorithm for graph partitioning used in community detection (Newman, 2006) applied to optimize a core-structure objective described in Borgatti and Everett (2000). Reference: Borgatti and Everett (2000) Soc Networks 21:375–395. Newman (2006) Phys Rev E 74:036104, PNAS 23:8577-8582. Rubinov, Ypma et al. (2015) PNAS 112:10032-7 2015, Mika Rubinov, U Cambridge ``` cycprob ``` CYCPROB Cycle probability [fcyc,pcyc] = cycprob(Pq); Cycles are paths which begin and end at the same node. Cycle probability for path length d, is the fraction of all paths of length d-1 that may be extended to form cycles of length d. Input: Pq, 3D matrix, with Pq(i,j,q) = number of paths from 'i' to 'j' of length 'q' (produced by 'findpaths') Outputs: fcyc, fraction of all paths that are cycles for each path length 'q'. pcyc, probability that a non-cyclic path of length 'q-1' can be extended to form a cycle of length 'q', for each path length 'q', Olaf Sporns, Indiana University, 2002/2007/2008 ``` degrees_dir ``` DEGREES_DIR Indegree and outdegree [id,od,deg] = degrees_dir(CIJ); Node degree is the number of links connected to the node. The indegree is the number of inward links and the outdegree is the number of outward links. Input: CIJ, directed (binary/weighted) connection matrix Output: id, node indegree od, node outdegree deg, node degree (indegree + outdegree) Notes: Inputs are assumed to be on the columns of the CIJ matrix. Weight information is discarded. Olaf Sporns, Indiana University, 2002/2006/2008 ``` degrees_und ``` DEGREES_UND Degree deg = degrees_und(CIJ); Node degree is the number of links connected to the node. Input: CIJ, undirected (binary/weighted) connection matrix Output: deg, node degree Note: Weight information is discarded. Olaf Sporns, Indiana University, 2002/2006/2008 ``` density_dir ``` DENSITY_DIR Density kden = density_dir(CIJ); [kden,N,K] = density_dir(CIJ); Density is the fraction of present connections to possible connections. Input: CIJ, directed (weighted/binary) connection matrix Output: kden, density N, number of vertices K, number of edges Notes: Assumes CIJ is directed and has no self-connections. Weight information is discarded. Olaf Sporns, Indiana University, 2002/2007/2008 ``` density_und ``` DENSITY_UND Density kden = density_und(CIJ); [kden,N,K] = density_und(CIJ); Density is the fraction of present connections to possible connections. Input: CIJ, undirected (weighted/binary) connection matrix Output: kden, density N, number of vertices K, number of edges Notes: Assumes CIJ is undirected and has no self-connections. Weight information is discarded. Olaf Sporns, Indiana University, 2002/2007/2008 ``` diffusion_efficiency ``` DIFFUSION_EFFICIENCY Global mean and pair-wise diffusion efficiency [GEdiff,Ediff] = diffusion_efficiency(adj); The diffusion efficiency between nodes i and j is the inverse of the mean first passage time from i to j, that is the expected number of steps it takes a random walker starting at node i to arrive for the first time at node j. Note that the mean first passage time is not a symmetric measure -- mfpt(i,j) may be different from mfpt(j,i) -- and the pair-wise diffusion efficiency matrix is hence also not symmetric. Input: adj, Weighted/Unweighted, directed/undirected adjacency matrix Outputs: GEdiff, Mean Global diffusion efficiency (scalar) Ediff, Pair-wise diffusion efficiency (matrix) References: Goñi J, et al (2013) PLoS ONE Joaquin Goñi and Andrea Avena-Koenigsberger, IU Bloomington, 2012 ``` distance_bin ``` DISTANCE_BIN Distance matrix D = distance_bin(A); The distance matrix contains lengths of shortest paths between all pairs of nodes. An entry (u,v) represents the length of shortest path from node u to node v. The average shortest path length is the characteristic path length of the network. Input: A, binary directed/undirected connection matrix Output: D, distance matrix Notes: Lengths between disconnected nodes are set to Inf. Lengths on the main diagonal are set to 0. Algorithm: Algebraic shortest paths. Mika Rubinov, U Cambridge Jonathan Clayden, UCL 2007-2013 ``` distance_wei ``` DISTANCE_WEI Distance matrix (Dijkstra's algorithm) D = distance_wei(L); [D,B] = distance_wei(L); The distance matrix contains lengths of shortest paths between all pairs of nodes. An entry (u,v) represents the length of shortest path from node u to node v. The average shortest path length is the characteristic path length of the network. Input: L, Directed/undirected connection-length matrix. *** NB: The length matrix L isn't the weights matrix W (see below) *** Output: D, distance (shortest weighted path) matrix B, number of edges in shortest weighted path matrix Notes: The input matrix must be a connection-length matrix, typically obtained via a mapping from weight to length. For instance, in a weighted correlation network higher correlations are more naturally interpreted as shorter distances and the input matrix should consequently be some inverse of the connectivity matrix. The number of edges in shortest weighted paths may in general exceed the number of edges in shortest binary paths (i.e. shortest paths computed on the binarized connectivity matrix), because shortest weighted paths have the minimal weighted distance, but not necessarily the minimal number of edges. Lengths between disconnected nodes are set to Inf. Lengths on the main diagonal are set to 0. Algorithm: Dijkstra's algorithm. Mika Rubinov, UNSW/U Cambridge, 2007-2012. Rick Betzel and Andrea Avena, IU, 2012 ``` distance_wei_floyd ``` DISTANCE_WEI_FLOYD Distance matrix (Floyd-Warshall algorithm) [SPL,hops,Pmat] = distance_wei_floyd(D,transform) Computes the topological length of the shortest possible path connecting every pair of nodes in the network. Inputs: D, Weighted/unweighted directed/undirected connection *weight* OR *length* matrix. transform, If the input matrix is a connection *weight* matrix, specify a transform that map input connection weights to connection lengths. Two transforms are available. 'log' -> l_ij = -log(w_ij) 'inv' -> l_ij = 1/w_ij If the input matrix is a connection *length* matrix, do not specify a transform (or specify an empty transform argument). Outputs: SPL, Unweighted/Weighted shortest path-length matrix. If W is directed matrix, then SPL is not symmetric. hops, Number of edges in the shortest path matrix. If W is unweighted, SPL and hops are identical. Pmat, Elements {i,j} of this matrix indicate the next node in the shortest path between i and j. This matrix is used as an input argument for function 'retrieve_shortest_path.m', which returns as output the sequence of nodes comprising the shortest path between a given pair of nodes. Notes: There may be more than one shortest path between any pair of nodes in the network. Non-unique shortest paths are termed shortest path degeneracies, and are most likely to occur in unweighted networks. When the shortest-path is degenerate, The elements of matrix Pmat correspond to the first shortest path discovered by the algorithm. The input matrix may be either a connection weight matrix, or a connection length matrix. The connection length matrix is typically obtained with a mapping from weight to length, such that higher weights are mapped to shorter lengths (see above). Algorithm: Floyd–Warshall Algorithm Andrea Avena-Koenigsberger, IU, 2012 ``` diversity_coef_sign ``` DIVERSITY_COEF_SIGN Shannon-entropy based diversity coefficient [Hpos Hneg] = diversity_coef_sign(W,Ci); The Shannon-entropy based diversity coefficient measures the diversity of intermodular connections of individual nodes and ranges from 0 to 1. Inputs: W, undirected connection matrix with positive and negative weights Ci, community affiliation vector Output: Hpos, diversity coefficient based on positive connections Hneg, diversity coefficient based on negative connections References: Shannon CE (1948) Bell Syst Tech J 27, 379-423. Rubinov and Sporns (2011) NeuroImage. 2011-2012, Mika Rubinov, U Cambridge ``` edge_betweenness_bin ``` EDGE_BETWEENNESS_BIN Edge betweenness centrality EBC = edge_betweenness_bin(A); [EBC BC] = edge_betweenness_bin(A); Edge betweenness centrality is the fraction of all shortest paths in the network that contain a given edge. Edges with high values of betweenness centrality participate in a large number of shortest paths. Input: A, binary (directed/undirected) connection matrix. Output: EBC, edge betweenness centrality matrix. BC, node betweenness centrality vector. Note: Betweenness centrality may be normalised to the range [0,1] as BC/[(N-1)(N-2)], where N is the number of nodes in the network. Reference: Brandes (2001) J Math Sociol 25:163-177. Mika Rubinov, UNSW/U Cambridge, 2007-2012 ``` edge_betweenness_wei ``` EDGE_BETWEENNESS_WEI Edge betweenness centrality EBC = edge_betweenness_wei(L); [EBC BC] = edge_betweenness_wei(L); Edge betweenness centrality is the fraction of all shortest paths in the network that contain a given edge. Edges with high values of betweenness centrality participate in a large number of shortest paths. Input: L, Directed/undirected connection-length matrix. Output: EBC, edge betweenness centrality matrix. BC, nodal betweenness centrality vector. Notes: The input matrix must be a connection-length matrix, typically obtained via a mapping from weight to length. For instance, in a weighted correlation network higher correlations are more naturally interpreted as shorter distances and the input matrix should consequently be some inverse of the connectivity matrix. Betweenness centrality may be normalised to the range [0,1] as BC/[(N-1)(N-2)], where N is the number of nodes in the network. Reference: Brandes (2001) J Math Sociol 25:163-177. Mika Rubinov, UNSW/U Cambridge, 2007-2012 ``` edge_nei_overlap_bd ``` EDGE_NEI_OVERLAP_BD Overlap amongst neighbors of two adjacent nodes [EC,ec,degij] = edge_nei_bd(CIJ); This function determines the neighbors of two nodes that are linked by an edge, and then computes their overlap. Connection matrix must be binary and directed. Entries of 'EC' that are 'inf' indicate that no edge is present. Entries of 'EC' that are 0 denote "local bridges", i.e. edges that link completely non-overlapping neighborhoods. Low values of EC indicate edges that are "weak ties". If CIJ is weighted, the weights are ignored. Neighbors of a node can be linked by incoming, outgoing, or reciprocal connections. Inputs: CIJ, directed (binary/weighted) connection matrix Outputs: EC, edge neighborhood overlap matrix ec, edge neighborhood overlap per edge, in vector format degij, degrees of node pairs connected by each edge Reference: Easley and Kleinberg (2010) Networks, Crowds, and Markets. Cambridge University Press, Chapter 3 Olaf Sporns, Indiana University, 2012 ``` edge_nei_overlap_bu ``` EDGE_NEI_OVERLAP_BU Overlap amongst neighbors of two adjacent nodes [EC,ec,degij] = edge_nei_bu(CIJ); This function determines the neighbors of two nodes that are linked by an edge, and then computes their overlap. Connection matrix must be binary and directed. Entries of 'EC' that are 'inf' indicate that no edge is present. Entries of 'EC' that are 0 denote "local bridges", i.e. edges that link completely non-overlapping neighborhoods. Low values of EC indicate edges that are "weak ties". If CIJ is weighted, the weights are ignored. Inputs: CIJ, undirected (binary/weighted) connection matrix Outputs: EC, edge neighborhood overlap matrix ec, edge neighborhood overlap per edge, in vector format degij, degrees of node pairs connected by each edge Reference: Easley and Kleinberg (2010) Networks, Crowds, and Markets. Cambridge University Press, Chapter 3. Olaf Sporns, Indiana University, 2012 ``` efficiency_bin ``` EFFICIENCY_BIN Global efficiency, local efficiency. Eglob = efficiency_bin(A); Eloc = efficiency_bin(A,1); The global efficiency is the average of inverse shortest path length, and is inversely related to the characteristic path length. The local efficiency is the global efficiency computed on the neighborhood of the node, and is related to the clustering coefficient. Inputs: A, binary undirected or directed connection matrix local, optional argument local=0 computes global efficiency (default) local=1 computes local efficiency Output: Eglob, global efficiency (scalar) Eloc, local efficiency (vector) Algorithm: algebraic path count Reference: Latora and Marchiori (2001) Phys Rev Lett 87:198701. Fagiolo (2007) Phys Rev E 76:026107. Rubinov M, Sporns O (2010) NeuroImage 52:1059-69 Mika Rubinov, U Cambridge Jonathan Clayden, UCL 2008-2013 ``` efficiency_wei ``` EFFICIENCY_WEI Global efficiency, local efficiency. Eglob = efficiency_wei(W); Eloc = efficiency_wei(W,2); The global efficiency is the average of inverse shortest path length, and is inversely related to the characteristic path length. The local efficiency is the global efficiency computed on the neighborhood of the node, and is related to the clustering coefficient. Inputs: W, weighted undirected or directed connection matrix local, optional argument local=0 computes the global efficiency (default). local=1 computes the original version of the local efficiency. local=2 computes the modified version of the local efficiency (recommended, see below). Output: Eglob, global efficiency (scalar) Eloc, local efficiency (vector) Notes: The efficiency is computed using an auxiliary connection-length matrix L, defined as L_ij = 1/W_ij for all nonzero L_ij; This has an intuitive interpretation, as higher connection weights intuitively correspond to shorter lengths. The weighted local efficiency broadly parallels the weighted clustering coefficient of Onnela et al. (2005) and distinguishes the influence of different paths based on connection weights of the corresponding neighbors to the node in question. In other words, a path between two neighbors with strong connections to the node in question contributes more to the local efficiency than a path between two weakly connected neighbors. Note that the original weighted variant of the local efficiency (described in Rubinov and Sporns, 2010) is not a true generalization of the binary variant, while the modified variant (described in Wang et al., 2016) is a true generalization. For ease of interpretation of the local efficiency it may be advantageous to rescale all weights to lie between 0 and 1. Algorithm: Dijkstra's algorithm References: Latora and Marchiori (2001) Phys Rev Lett 87:198701. Onnela et al. (2005) Phys Rev E 71:065103 Fagiolo (2007) Phys Rev E 76:026107. Rubinov M, Sporns O (2010) NeuroImage 52:1059-69 Wang Y et al. (2016) Neural Comput 21:1-19. Mika Rubinov, U Cambridge/Janelia HHMI, 2011-2017 ``` eigenvector_centrality_und ``` EIGENVECTOR_CENTRALITY_UND Spectral measure of centrality v = eigenvector_centrality_und(CIJ) Eigenector centrality is a self-referential measure of centrality: nodes have high eigenvector centrality if they connect to other nodes that have high eigenvector centrality. The eigenvector centrality of node i is equivalent to the ith element in the eigenvector corresponding to the largest eigenvalue of the adjacency matrix. Inputs: CIJ, binary/weighted undirected adjacency matrix. Outputs: v, eigenvector associated with the largest eigenvalue of the adjacency matrix CIJ. Reference: Newman, MEJ (2002). The mathematics of networks. Contributors: Xi-Nian Zuo, Chinese Academy of Sciences, 2010 Rick Betzel, Indiana University, 2012 Mika Rubinov, University of Cambridge, 2015 ``` erange ``` ERANGE Shortcuts [Erange,eta,Eshort,fs] = erange(CIJ); Shorcuts are central edges which significantly reduce the characteristic path length in the network. Input: CIJ, binary directed connection matrix Outputs: Erange, range for each edge, i.e. the length of the shortest path from i to j for edge c(i,j) AFTER the edge has been removed from the graph. eta average range for entire graph. Eshort entries are ones for shortcut edges. fs fraction of shortcuts in the graph. Follows the treatment of 'shortcuts' by Duncan Watts Olaf Sporns, Indiana University, 2002/2007/2008 ``` evaluate_generative_model ``` EVALUATE_GENERATIVE_MODEL Generation and evaluation of synthetic networks [B,E,K] = EVALUATE_GENERATIVE_MODEL(A,Atgt,D,m,modeltype,modelvar,params) Generates synthetic networks and evaluates their energy function (see below) using the models described in the study by Betzel et al (2016) in Neuroimage. Inputs: A, binary network of seed connections Atgt, binary network against which synthetic networks are compared D, Euclidean distance/fiber length matrix m, number of connections that should be present in final synthetic network modeltype, specifies the generative rule (see below) modelvar, specifies whether the generative rules are based on power-law or exponential relationship ({'powerlaw'}|{'exponential}) params, either a vector (in the case of the geometric model) or a matrix (for all other models) of parameters at which the model should be evaluated. Outputs: B, m x number of networks matrix of connections E, energy for each synthetic network K, Kolmogorov-Smirnov statistics for each synthetic network. Full list of model types: (each model type realizes a different generative rule) 1. 'sptl' spatial model 2. 'neighbors' number of common neighbors 3. 'matching' matching index 4. 'clu-avg' average clustering coeff. 5. 'clu-min' minimum clustering coeff. 6. 'clu-max' maximum clustering coeff. 7. 'clu-diff' difference in clustering coeff. 8. 'clu-prod' product of clustering coeff. 9. 'deg-avg' average degree 10. 'deg-min' minimum degree 11. 'deg-max' maximum degree 12. 'deg-diff' difference in degree 13. 'deg-prod' product of degree Note: Energy is calculated in exactly the same way as in Betzel et al (2016). There are four components to the energy are KS statistics comparing degree, clustering coefficient, betweenness centrality, and edge length distributions. Energy is calculated as the maximum across all four statistics. Reference: Betzel et al (2016) Neuroimage 124:1054-64. Richard Betzel, Indiana University/University of Pennsylvania, 2015 ``` find_motif34 ``` FIND_MOTIF34 Motif legend Motif_matrices = find_motif34(Motif_id,Motif_class); Motif_id = find_motif34(Motif_matrix); This function returns all motif isomorphs for a given motif id and class (3 or 4). The function also returns the motif id for a given motif matrix 1. Input: Motif_id, e.g. 1 to 13, if class is 3 Motif_class, number of nodes, 3 or 4. Output: Motif_matrices, all isomorphs for the given motif 2. Input: Motif_matrix e.g. [0 1 0; 0 0 1; 1 0 0] Output Motif_id e.g. 1 to 13, if class is 3 Mika Rubinov, UNSW, 2007-2008 ``` findpaths ``` FINDPATHS Network paths [Pq,tpath,plq,qstop,allpths,util] = findpaths(CIJ,sources,qmax,savepths); Paths are sequences of linked nodes, that never visit a single node more than once. This function finds all paths that start at a set of source nodes, up to a specified length. Warning: very memory-intensive. Inputs: CIJ, binary (directed/undirected) connection matrix qmax, maximal path length sources, source units from which paths are grown savepths, set to 1 if all paths are to be collected in 'allpths' Outputs: Pq, 3D matrix, with P(i,j,q) = number of paths from 'i' to 'j' of length 'q'. tpath, total number of paths found (lengths 1 to 'qmax') plq, path length distribution as a function of 'q' qstop, path length at which 'findpaths' is stopped allpths, a matrix containing all paths up to 'qmax' util, node use index Note that Pq(:,:,N) can only carry entries on the diagonal, as all "legal" paths of length N-1 must terminate. Cycles of length N are possible, with all vertices visited exactly once (except for source and target). 'qmax = N' can wreak havoc (due to memory problems). Note: Weights are discarded. Note: I am certain that this algorithm is rather inefficient - suggestions for improvements are welcome. Olaf Sporns, Indiana University, 2002/2007/2008/2010 ``` findwalks ``` FINDWALKS Network walks [Wq,twalk,wlq] = findwalks(CIJ); Walks are sequences of linked nodes, that may visit a single node more than once. This function finds the number of walks of a given length, between any two nodes. Input: CIJ binary (directed/undirected) connection matrix Outputs: Wq 3D matrix, Wq(i,j,q) is the number of walks from 'i' to 'j' of length 'q'. twalk total number of walks found wlq walk length distribution as function of 'q' Notes: Wq grows very quickly for larger N,K,q. Weights are discarded. Algorithm: algebraic path count Olaf Sporns, Indiana University, 2002/2007/2008 ``` flow_coef_bd ``` FLOW_COEF_BD Node-wise flow coefficients [hc,HC,total_flo] = flow_coef_bd(CIJ) Computes the flow coefficient for each node and averaged over the network, as described in Honey et al. (2007) PNAS. The flow coefficient is similar to betweenness centrality, but works on a local neighborhood. It is mathematically related to the clustering coefficient (cc) at each node as, fc+cc <= 1. input: CIJ, connection/adjacency matrix (binary, directed) output: fc, flow coefficient for each node FC, average flow coefficient over the network total_flo, number of paths that "flow" across the central node Reference: Honey et al. (2007) Proc Natl Acad Sci U S A Olaf Sporns, Indiana University, 2007/2010/2012 ``` gateway_coef_sign ``` GATEWAY_COEF_SIGN Gateway coefficient [Gpos,Gneg] = gateway_coef_sign(W,Ci,centtype); Gateway coefficient is a variant of participation coefficient. Similar to participation coefficient, gateway coefficient measures the diversity of intermodular connections of individual nodes, but this is weighted by how critical these connections are to intermodular connectivity (e.g., if a node is the only connection between it's module and another module, it will have a higher gateway coefficient). Inputs: W, undirected connection matrix with positive and negative weights Ci, community affiliation vector centtype, centrality measure to use 1 = Node Strength 2 = Betweenness Centrality Output: Gpos, gateway coefficient for positive weights Gneg, gateway coefficient for negative weights Reference: Vargas ER, Wahl LM. Eur Phys J B (2014) 87:1-10. Jeff Spielberg, Boston University ``` generate_fc ``` GENERATE_FC Generation of synthetic functional connectivity matrices [FCpre,pred_data,Fcorr] = generate_fc(SC,beta,ED,{'SPLwei_log','SIwei_log'},FC) [FCpre,pred_data] = generate_fc(SC,beta,[],{'SPLwei_log','SIwei_log'}) Uses a vector beta of regression coefficients from the model FC = pred_var*beta to predict FC. pred_var are structural-based network measures derived from the structural connectivity network. Inputs: SC, Weighted/unweighted undirected NxN Structural Connectivity matrix. beta, Regression coefficients (vector). These may be obtained as an output parameter from function predict_fc.m ED, Euclidean distance matrix or upper triangular vector of the matrix (optional) pred_var, Set of M predictors. These can be given as an KxM array where K = ((N*(N-1))/2) and M is the number of predictors. Alternatively, pred_var can be a cell with the names of network measures to be used as predictors. Accepted network measure names are: SPLbin - Shortest-path length (binary) SPLwei_inv - Shortest-path length computed with an inv transform SPLwei_log - Shortest-path length computed with a log transform SPLdist - Shortest-path length computed with no transform SIbin - Search Information of binary shortest-paths SIwei_inv - Search Information of shortest-paths computed with an inv transform SIwei_log - Search Information of shortest-paths computed with a log transform SIdist - Search Information of shortest-paths computed with no transform T - Path Transitivity deltaMFPT - Column-wise z-scored mean first passage time neighOverlap - Neighborhood Overlap MI - Matching Index Predictors must be specified in the order that matches the given beta values. model, Specifies the order of the regression model used within matlab's function regstats.m. 'model' can be any option accepted by matlab's regstats.m function (e.g.'linear', 'interaction' 'quadratic', etc). If no model is specified, 'linear' is the default. FC, Functional connections. FC can be a NxN symmetric matrix or a ((N*(N-1))/2) x 1 vector containing the upper triangular elements of the square FC matrix (excluding diagonal elements). This argument is optional and only used to compute the correlation between the predicted FC and empirical FC. Outputs: FCpre, Predicted NxN Functional Connectivity matrix pred_data, KxM array of predictors. FCcorr, Pearson Correlation between FCpred and FC Reference: Goñi et al. (2014) PNAS 111: 833–838 Andrea Avena-Koenigsberger, Joaquin Goñi and Olaf Sporns; IU Bloomington, 2016 ``` generative_model ``` GENERATIVE_MODEL Run generative model code B = GENERATIVE_MODEL(A,D,m,modeltype,modelvar,params) Generates synthetic networks using the models described in the study by Betzel et al (2016) in Neuroimage. Inputs: A, binary network of seed connections D, Euclidean distance/fiber length matrix m, number of connections that should be present in final synthetic network modeltype, specifies the generative rule (see below) modelvar, specifies whether the generative rules are based on power-law or exponential relationship ({'powerlaw'}|{'exponential}) params, either a vector (in the case of the geometric model) or a matrix (for all other models) of parameters at which the model should be evaluated. epsilon, the baseline probability of forming a particular connection (should be a very small number {default = 1e-5}). Output: B, m x number of networks matrix of connections Full list of model types: (each model type realizes a different generative rule) 1. 'sptl' spatial model 2. 'neighbors' number of common neighbors 3. 'matching' matching index 4. 'clu-avg' average clustering coeff. 5. 'clu-min' minimum clustering coeff. 6. 'clu-max' maximum clustering coeff. 7. 'clu-diff' difference in clustering coeff. 8. 'clu-prod' product of clustering coeff. 9. 'deg-avg' average degree 10. 'deg-min' minimum degree 11. 'deg-max' maximum degree 12. 'deg-diff' difference in degree 13. 'deg-prod' product of degree Example usage: load demo_generative_models_data % get number of bi-directional connections m = nnz(A)/2; % get cardinality of network n = length(A); % set model type modeltype = 'neighbors'; % set whether the model is based on powerlaw or exponentials modelvar = [{'powerlaw'},{'powerlaw'}]; % choose some model parameters params = [-2,0.2; -5,1.2; -1,1.5]; nparams = size(params,1); % generate synthetic networks B = generative_model(Aseed,D,m,modeltype,modelvar,params); % store them in adjacency matrix format Asynth = zeros(n,n,nparams); for i = 1:nparams; a = zeros(n); a(B(:,i)) = 1; a = a + a'; Asynth(:,:,i) = a; end Reference: Betzel et al (2016) Neuroimage 124:1054-64. Richard Betzel, Indiana University/University of Pennsylvania, 2015 ``` get_components ``` GET_COMPONENTS Connected components [comps,comp_sizes] = get_components(adj); Returns the components of an undirected graph specified by the binary and undirected adjacency matrix adj. Components and their constitutent nodes are assigned the same index and stored in the vector, comps. The vector, comp_sizes, contains the number of nodes beloning to each component. Inputs: adj, binary and undirected adjacency matrix Outputs: comps, vector of component assignments for each node comp_sizes, vector of component sizes Note: disconnected nodes will appear as components of size 1 J Goni, University of Navarra and Indiana University, 2009/2011 ``` grid_communities ``` GRID_COMMUNITIES Outline communities along diagonal [X Y INDSORT] = GRID_COMMUNITIES(C) takes a vector of community assignments C and returns three output arguments for visualizing the communities. The third is INDSORT, which is an ordering of the vertices so that nodes with the same community assignment are next to one another. The first two arguments are vectors that, when overlaid on the adjacency matrix using the PLOT function, highlight the communities. Example: >> load AIJ; % load adjacency matrix >> [C,Q] = modularity_louvain_und(AIJ); % get community assignments >> [X,Y,INDSORT] = fcn_grid_communities(C); % call function >> imagesc(AIJ(INDSORT,INDSORT)); % plot ordered adjacency matrix >> hold on; % hold on to overlay community visualization >> plot(X,Y,'r','linewidth',2); % plot community boundaries Inputs: C, community assignments Outputs: X, x coor Y, y coor INDSORT, indices Richard Betzel, Indiana University, 2012 ``` gtom ``` GTOM Generalized topological overlap measure gt = gtom(adj,numSteps); The m-th step generalized topological overlap measure (GTOM) quantifies the extent to which a pair of nodes have similar m-th step neighbors. Mth-step neighbors are nodes that are reachable by a path of at most length m. This function computes the the M x M generalized topological overlap measure (GTOM) matrix for number of steps, numSteps. Inputs: adj, adjacency matrix (binary,undirected) numSteps, number of steps Outputs: gt, GTOM matrix NOTE: When numSteps is equal to 1, GTOM is identical to the topological overlap measure (TOM) from reference [2]. In that case the 'gt' matrix records, for each pair of nodes, the fraction of neighbors the two nodes share in common, where "neighbors" are one step removed. As 'numSteps' is increased, neighbors that are furter out are considered. Elements of 'gt' are bounded between 0 and 1. The 'gt' matrix can be converted from a similarity to a distance matrix by taking 1-gt. References: [1] Yip & Horvath (2007) BMC Bioinformatics 2007, 8:22 [2] Ravasz et al (2002) Science 297 (5586), 1551. J Goni, University of Navarra and Indiana University, 2009/2011 ``` jdegree ``` JDEGREE Joint degree distribution [J,J_od,J_id,J_bl] = jdegree(CIJ); This function returns a matrix in which the value of each element (u,v) corresponds to the number of nodes that have u outgoing connections and v incoming connections. Input: CIJ, directed (weighted/binary) connection matrix Outputs: J, joint degree distribution matrix (shifted by one) J_od, number of vertices with od>id. J_id, number of vertices with id>od. J_bl, number of vertices with id=od. Note: Weights are discarded. Olaf Sporns, Indiana University, 2002/2006/2008 ``` kcore_bd ``` KCORE_BD K-core [CIJkcore,kn,peelorder,peellevel] = kcore_bd(CIJ,k); The k-core is the largest subnetwork comprising nodes of degree at least k. This function computes the k-core for a given binary directed connection matrix by recursively peeling off nodes with degree lower than k, until no such nodes remain. input: CIJ, connection/adjacency matrix (binary, directed) k, level of k-core output: CIJkcore, connection matrix of the k-core. This matrix only contains nodes of degree at least k. kn, size of k-core peelorder, indices in the order in which they were peeled away during k-core decomposition peellevel, corresponding level - nodes at the same level have been peeled away at the same time 'peelorder' and 'peellevel' are similar the the k-core sub-shells described in Modha and Singh (2010). References: e.g. Hagmann et al. (2008) PLoS Biology Olaf Sporns, Indiana University, 2007/2008/2010/2012 ``` kcore_bu ``` KCORE_BU K-core [CIJkcore,kn,peelorder,peellevel] = kcore_bu(CIJ,k); The k-core is the largest subnetwork comprising nodes of degree at least k. This function computes the k-core for a given binary undirected connection matrix by recursively peeling off nodes with degree lower than k, until no such nodes remain. input: CIJ, connection/adjacency matrix (binary, undirected) k, level of k-core output: CIJkcore, connection matrix of the k-core. This matrix only contains nodes of degree at least k. kn, size of k-core peelorder, indices in the order in which they were peeled away during k-core decomposition peellevel, corresponding level - nodes at the same level were peeled away at the same time 'peelorder' and 'peellevel' are similar the the k-core sub-shells described in Modha and Singh (2010). References: e.g. Hagmann et al. (2008) PLoS Biology Olaf Sporns, Indiana University, 2007/2008/2010/2012 ``` kcoreness_centrality_bd ``` KCORENESS_CENTRALITY_BD K-coreness centrality [coreness,kn] = kcoreness_centrality_bd(CIJ) The k-core is the largest subgraph comprising nodes of degree at least k. The coreness of a node is k if the node belongs to the k-core but not to the (k+1)-core. This function computes k-coreness of all nodes for a given binary directed connection matrix. input: CIJ, connection/adjacency matrix (binary, directed) output: coreness, node coreness. kn, size of k-core References: e.g. Hagmann et al. (2008) PLoS Biology Olaf Sporns, Indiana University, 2007/2008/2010/2012 ``` kcoreness_centrality_bu ``` KCORENESS_CENTRALITY_BU K-coreness centrality [coreness,kn] = kcoreness_centrality_bu(CIJ) The k-core is the largest subgraph comprising nodes of degree at least k. The coreness of a node is k if the node belongs to the k-core but not to the (k+1)-core. This function computes the coreness of all nodes for a given binary undirected connection matrix. input: CIJ, connection/adjacency matrix (binary, undirected) output: coreness, node coreness. kn, size of k-core References: e.g. Hagmann et al. (2008) PLoS Biology Olaf Sporns, Indiana University, 2007/2008/2010/2012 ``` latmio_dir ``` LATMIO_DIR Lattice with preserved in/out degree distribution [Rlatt,Rrp,ind_rp,eff] = latmio_dir(R,ITER,D); This function "latticizes" a directed network, while preserving the in- and out-degree distributions. In weighted networks, the function preserves the out-strength but not the in-strength distributions. Input: R, directed (binary/weighted) connection matrix ITER, rewiring parameter (each edge is rewired approximately ITER times) D, distance-to-diagonal matrix Output: Rlatt, latticized network in original node ordering Rrp, latticized network in node ordering used for latticization ind_rp, node ordering used for latticization eff, number of actual rewirings carried out References: Maslov and Sneppen (2002) Science 296:910 Sporns and Zwi (2004) Neuroinformatics 2:145 Mika Rubinov, UNSW, 2007-2010 Olaf Sporns, IU, 2012 ``` latmio_dir_connected ``` LATMIO_DIR_CONNECTED Lattice with preserved in/out degree distribution [Rlatt,Rrp,ind_rp,eff] = latmio_dir_connected(R,ITER,D); This function "latticizes" a directed network, while preserving the in- and out-degree distributions. In weighted networks, the function preserves the out-strength but not the in-strength distributions. The function also ensures that the randomized network maintains connectedness, the ability for every node to reach every other node in the network. The input network for this function must be connected. Input: R, directed (binary/weighted) connection matrix ITER, rewiring parameter (each edge is rewired approximately ITER times) D, distance-to-diagonal matrix Output: Rlatt, latticized network in original node ordering Rrp, latticized network in node ordering used for latticization ind_rp, node ordering used for latticization eff, number of actual rewirings carried out References: Maslov and Sneppen (2002) Science 296:910 Sporns and Zwi (2004) Neuroinformatics 2:145 Mika Rubinov, UNSW, 2007-2010 Olaf Sporns, Indiana University, 2012 ``` latmio_und ``` LATMIO_UND Lattice with preserved degree distribution [Rlatt,Rrp,ind_rp,eff] = latmio_und(R,ITER,D); This function "latticizes" an undirected network, while preserving the degree distribution. The function does not preserve the strength distribution in weighted networks. Input: R, undirected (binary/weighted) connection matrix ITER, rewiring parameter (each edge is rewired approximately ITER times) D, distance-to-diagonal matrix Output: Rlatt, latticized network in original node ordering Rrp, latticized network in node ordering used for latticization ind_rp, node ordering used for latticization eff, number of actual rewirings carried out References: Maslov and Sneppen (2002) Science 296:910 Sporns and Zwi (2004) Neuroinformatics 2:145 2007-2012 Mika Rubinov, UNSW Jonathan Power, WUSTL Olaf Sporns, IU ``` latmio_und_connected ``` LATMIO_UND_CONNECTED Lattice with preserved degree distribution [Rlatt,Rrp,ind_rp,eff] = latmio_und_connected(R,ITER,D); This function "latticizes" an undirected network, while preserving the degree distribution. The function does not preserve the strength distribution in weighted networks. The function also ensures that the randomized network maintains connectedness, the ability for every node to reach every other node in the network. The input network for this function must be connected. Input: R, undirected (binary/weighted) connection matrix ITER, rewiring parameter (each edge is rewired approximately ITER times) D, distance-to-diagonal matrix Output: Rlatt, latticized network in original node ordering Rrp, latticized network in node ordering used for latticization ind_rp, node ordering used for latticization eff, number of actual rewirings carried out References: Maslov and Sneppen (2002) Science 296:910 Sporns and Zwi (2004) Neuroinformatics 2:145 2007-2012 Mika Rubinov, UNSW Jonathan Power, WUSTL Olaf Sporns, IU ``` link_communities ``` LINK_COMMUNITIES Optimal overlapping community structure M = link_communities(W) M = link_communities(W,'complete'); The optimal community structure is a subdivision of the network into groups of nodes which have a high number of within-group connections and a low number of between group connections. This algorithm uncovers overlapping community structure via hierarchical clustering of network links. This algorith is generalized for weighted/directed/fully-connected networks. Input: W, directed (weighted or binary) connection matrix. type_clustering, type of hierarchical clustering (optional) 'single' single-linkage (default) 'complete' complete-linkage Output: M, nodal community-affiliation matrix binary matrix of size CxN [communities x nodes] NB: The algorithm can be slow and memory intensive. Reference: Ahn, Bagrow and Lehmann (2010) Nature 466, 761–764. Mika Rubinov, U Cambridge, 2014-2015 ``` local_assortativity_wu_sign ``` LOCAL_ASSORTATIVITY_WU_SIGN Local Assortativity [loc_assort_pos,loc_assort_neg] = local_assortativity_wu_sign(W); Local Assortativity measures the extent to which nodes are connected to nodes of similar strength (vs. higher or lower strength). Adapted from Thedchanamoorthy et al. (2014)'s formula to allow weighted/signed networks (node degree replaced with node strength). Note, output values sum to total assortativity. Inputs: W, undirected connection matrix with positive and negative weights Output: loc_assort_pos, local assortativity from positive weights loc_assort_neg, local assortativity from negative weights Reference: Thedchanamoorthy G, Piraveenan M, Kasthuriratna D, Senanayake U. Proc Comp Sci (2014) 29:2449-2461. Jeff Spielberg, Boston University ``` make_motif34lib ``` MAKE_MOTIF34LIB Auxiliary motif library function make_motif34lib; This function generates the motif34lib.mat library required for all other motif computations. Mika Rubinov, UNSW, 2007-2010 ``` makeevenCIJ ``` MAKEEVENCIJ Synthetic modular small-world network CIJ = makeevenCIJ(N,K,sz_cl); This function generates a random, directed network with a specified number of fully connected modules linked together by evenly distributed remaining random connections. Inputs: N, number of vertices (must be power of 2) K, number of edges sz_cl, size of clusters (power of 2) Outputs: CIJ, connection matrix Notes: N must be a power of 2. A warning is generated if all modules contain more edges than K. Cluster size is 2^sz_cl; Olaf Sporns, Indiana University, 2005/2007 ``` makefractalCIJ ``` MAKEFRACTALCIJ Synthetic hierarchical modular network [CIJ,K] = makefractalCIJ(mx_lvl,E,sz_cl); This function generates a directed network with a hierarchical modular organization. All modules are fully connected and connection density decays as 1/(E^n), with n = index of hierarchical level. Inputs: mx_lvl, number of hierarchical levels, N = 2^mx_lvl E, connection density fall-off per level sz_cl, size of clusters (power of 2) Outputs: CIJ, connection matrix K, number of connections present in the output CIJ Olaf Sporns, Indiana University, 2005/2007 ``` makelatticeCIJ ``` MAKELATTICECIJ Synthetic lattice network CIJ = makelatticeCIJ(N,K); This function generates a directed lattice network without toroidal boundary counditions (i.e. no ring-like "wrapping around"). Inputs: N, number of vertices K, number of edges Outputs: CIJ, connection matrix Note: The lattice is made by placing connections as close as possible to the main diagonal, without wrapping around. No connections are made on the main diagonal. In/Outdegree is kept approx. constant at K/N. Olaf Sporns, Indiana University, 2005/2007 ``` makerandCIJ_dir ``` MAKERANDCIJ_DIR Synthetic directed random network CIJ = makerandCIJ_dir(N,K); This function generates a directed random network Inputs: N, number of vertices K, number of edges Output: CIJ, directed random connection matrix Note: no connections are placed on the main diagonal. Olaf Sporns, Indiana University, 2007/2008 ``` makerandCIJ_und ``` MAKERANDCIJ_UND Synthetic directed random network CIJ = makerandCIJ_und(N,K); This function generates an undirected random network Inputs: N, number of vertices K, number of edges Output: CIJ, undirected random connection matrix Note: no connections are placed on the main diagonal. Olaf Sporns, Indiana University, 2007/2008 ``` makerandCIJdegreesfixed ``` MAKERANDCIJDEGREESFIXED Synthetic directed random network CIJ = makerandCIJdegreesfixed(N,K); This function generates a directed random network with a specified in-degree and out-degree sequence. The function returns a flag, denoting whether the algorithm succeeded or failed. Inputs: in, indegree vector out, outdegree vector Output: CIJ, binary directed connectivity matrix flag, flag=1 if the algorithm succeeded; flag=0 otherwise Notes: Necessary conditions include: length(in) = length(out) = n sum(in) = sum(out) = k in(i), out(i) < n-1 in(i) + out(j) < n+2 in(i) + out(i) < n No connections are placed on the main diagonal Aviad Rubinstein, Indiana University 2005/2007 ``` makeringlatticeCIJ ``` MAKERINGLATTICECIJ Synthetic lattice network CIJ = makeringlatticeCIJ(N,K); This function generates a directed lattice network with toroidal boundary counditions (i.e. with ring-like "wrapping around"). Inputs: N, number of vertices K, number of edges Outputs: CIJ, connection matrix Note: The lattice is made by placing connections as close as possible to the main diagonal, with wrapping around. No connections are made on the main diagonal. In/Outdegree is kept approx. constant at K/N. Olaf Sporns, Indiana University, 2005/2007 ``` maketoeplitzCIJ ``` MAKETOEPLITZCIJ A synthetic directed network with Gaussian drop-off of connectivity with distance CIJ = maketoeprandCIJ(N,K,s) This function generates a directed network with a Gaussian drop-off in edge density with increasing distance from the main diagonal. There are toroidal boundary counditions (i.e. no ring-like "wrapping around"). Inputs: N, number of vertices K, number of edges s, standard deviation of toeplitz Output: CIJ, connection matrix Note: no connections are placed on the main diagonal. Olaf Sporns, Indiana University, 2005/2007 ``` matching_ind ``` MATCHING_IND Matching index [Min,Mout,Mall] = matching_ind(CIJ); For any two nodes u and v, the matching index computes the amount of overlap in the connection patterns of u and v. Self-connections and u-v connections are ignored. The matching index is a symmetric quantity, similar to a correlation or a dot product. Input: CIJ, connection/adjacency matrix Output: Min, matching index for incoming connections Mout, matching index for outgoing connections Mall, matching index for all connections Notes: Does not use self- or cross connections for comparison. Does not use connections that are not present in BOTH u and v. All output matrices are calculated for upper triangular only. Olaf Sporns, Indiana University, 2002/2007/2008 ``` matching_ind_und ``` MATCHING_IND_UND matching index M0 = MATCHING_IND_UND(CIJ) computes matching index for undirected graph specified by adjacency matrix CIJ. Matching index is a measure of similarity between two nodes' connectivity profiles (excluding their mutual connection, should it exist). Inputs: CIJ, undirected adjacency matrix Outputs: M0, matching index matrix. Richard Betzel, Indiana University, 2013 ``` mean_first_passage_time ``` MEAN_FIRST_PASSAGE_TIME Mean first passage time MFPT = mean_first_passage_time(adj) The first passage time (MFPT) from i to j is the expected number of steps it takes a random walker starting at node i to arrive for the first time at node j. The mean first passage time is not a symmetric measure: mfpt(i,j) may be different from mfpt(j,i). Input: adj, Weighted/Unweighted, directed/undirected adjacency matrix Output: MFPT, Pairwise mean first passage time matrix. References: Goñi J, et al (2013) PLoS ONE Joaquin Goñi, IU Bloomington, 2012 ``` mleme_constraint_model ``` MLEME_CONSTRAINT_MODEL Unbiased sampling of networks with soft constraints W0 = mleme_constraint_model(samp, W); W0 = mleme_constraint_model(samp, W, M); W0 = mleme_constraint_model(samp, W, M, Lo, Li, Lm); [W0, E0, P0, Delt0] = mleme_constraint_model(samp, W, M, Lo, Li, Lm, opts); This function returns an ensemble of unbiasedly sampled networks with weighted node-strength and module-weight constraints. These constraints are soft in that they are satisfied on average for the full network ensemble but not, in general, for each individual network. Inputs (for a network with n nodes, m modules and c constraints): samp, Number of networks to sample. W, (length n) square directed and weighted connectivity matrix. All weights must be nonnegative integers. Note that real-valued weights may be converted to integers with arbitrary precision through rescaling and rounding, e.g. W_int = round(10^precision * W_real). M, (length n) module affiliation vector. This vector is often obtained as the output of a community detection algorithm. The vector must contain nonnegative integers, with zeros specifying nodes which are not part of any community. This input may be left empty if there are no module constraints. Lo, (length n) out-strength constraint logical vector. This vector specifies out-strength constraints for each node. Alternatively, it is possible to specify 1 to constrain all out-strengths or 0 for no constraints. Empty or no input results in default behavour (no constraints). Lo, (length n) in-strength constraint logical vector. This vector specifies in-strength constraints for each node. Alternatively, it is possible to specify 1 to constrain all in-strengths or 0 for no constraints. Empty or no input results in default behavour (no constraints). Lm, (length m) module-weight constraint logical matrix. This matrix specifies module-weight constraints for all pairs of modules. Alternatively, it is possible to specify 2 to constrain all inter-module and intra-module weights, 1 to constrain all intra-module weights, or 0 for no constraints. Empty or no input results in default behavour (no constraints). opts, optional argument: pass optimization and display options with optimset. Default: optimset('MaxFunEvals', 1e6*c, 'MaxIter', 1e6, 'Display', 'iter'); Outputs: W0, an ensemble of sampled networks with constraints. E0, expected weights matrix. P0, probability matrix. Delt0, algorithm convergence error. Algorithm: Maximum-likelihood estimation of network probability distribution by numerical solution of systems of nonlinear equations, and sampling of individual networks directly from this distribution. Notes: Empirical connection weights are not preserved. Constraint errors are guaranteed to vanish in the limit of the full network ensemble. Examples: % get community structure of a weighted network W M = community_louvain(W, 2); % specify node and module constraints n = length(W); % number of nodes m = max(M); % number of modules Lo = true(n, 1); % out-strength constraints Li = true(n, 1); % in-strength constraints Lm = eye(m); % module-weight constraints % sample networks with the above constraints [W0, E0, P0, Delt0] = mleme_constraint_model(samp, W, M, Lo, Li, Lm); % equivalent formulation [W0, E0, P0, Delt0] = mleme_constraint_model(samp, W, M, 1, 1, 1); % alternative: sample networks with average weight constraints only [W0, E0, P0, Delt0] = mleme_constraint_model(samp, W); References: Squartini and Garlaschelli (2011) New J Phys 13:083001 Rubinov (2016) Nat Commun 7:13812 2016, Mika Rubinov, Janelia HHMI ``` modularity_dir ``` MODULARITY_DIR Optimal community structure and modularity Ci = modularity_dir(W); [Ci Q] = modularity_dir(W); The optimal community structure is a subdivision of the network into nonoverlapping groups of nodes in a way that maximizes the number of within-group edges, and minimizes the number of between-group edges. The modularity is a statistic that quantifies the degree to which the network may be subdivided into such clearly delineated groups. Inputs: W, directed weighted/binary connection matrix gamma, resolution parameter (optional) gamma>1, detects smaller modules 0<=gamma<1, detects larger modules gamma=1, classic modularity (default) Outputs: Ci, optimal community structure Q, maximized modularity Note: This algorithm is essentially deterministic. The only potential source of stochasticity occurs at the iterative finetuning step, in the presence of non-unique optimal swaps. However, the present implementation always makes the first available optimal swap and is therefore deterministic. References: Leicht and Newman (2008) Phys Rev Lett 100:118703. Reichardt and Bornholdt (2006) Phys Rev E 74:016110. 2008-2016 Mika Rubinov, UNSW Jonathan Power, WUSTL Dani Bassett, UCSB Xindi Wang, Beijing Normal University Roan LaPlante, Martinos Center for Biomedical Imaging ``` modularity_und ``` MODULARITY_UND Optimal community structure and modularity Ci = modularity_und(W); [Ci Q] = modularity_und(W,gamma); The optimal community structure is a subdivision of the network into nonoverlapping groups of nodes in a way that maximizes the number of within-group edges, and minimizes the number of between-group edges. The modularity is a statistic that quantifies the degree to which the network may be subdivided into such clearly delineated groups. Inputs: W, undirected weighted/binary connection matrix gamma, resolution parameter (optional) gamma>1, detects smaller modules 0<=gamma<1, detects larger modules gamma=1, classic modularity (default) Outputs: Ci, optimal community structure Q, maximized modularity Note: This algorithm is essentially deterministic. The only potential source of stochasticity occurs at the iterative finetuning step, in the presence of non-unique optimal swaps. However, the present implementation always makes the first available optimal swap and is therefore deterministic. References: Newman (2006) -- Phys Rev E 74:036104, PNAS 23:8577-8582. Reichardt and Bornholdt (2006) Phys Rev E 74:016110. 2008-2016 Mika Rubinov, UNSW Jonathan Power, WUSTL Dani Bassett, UCSB Xindi Wang, Beijing Normal University Roan LaPlante, Martinos Center for Biomedical Imaging ``` module_degree_zscore ``` MODULE_DEGREE_ZSCORE Within-module degree z-score Z=module_degree_zscore(W,Ci,flag); The within-module degree z-score is a within-module version of degree centrality. Inputs: W, binary/weighted, directed/undirected connection matrix Ci, community affiliation vector flag, 0, undirected graph (default) 1, directed graph: out-degree 2, directed graph: in-degree 3, directed graph: out-degree and in-degree Output: Z, within-module degree z-score. Reference: Guimera R, Amaral L. Nature (2005) 433:895-900. Mika Rubinov, UNSW, 2008-2010 ``` motif3funct_bin ``` MOTIF3FUNCT_BIN Frequency of functional class-3 motifs [f,F] = motif3funct_bin(A); *Structural motifs* are patterns of local connectivity in complex networks. In contrast, *functional motifs* are all possible subsets of patterns of local connectivity embedded within structural motifs. Such patterns are particularly diverse in directed networks. The motif frequency of occurrence around an individual node is known as the motif fingerprint of that node. The total motif frequency of occurrence in the whole network is correspondingly known as the motif fingerprint of the network. Input: A, binary directed connection matrix Output: F, node motif frequency fingerprint f, network motif frequency fingerprint Notes: 1. The function find_motif34.m outputs the motif legend. 2. There is a source of possible confusion in motif terminology. Motifs ("structural" and "functional") are most frequently considered only in the context of anatomical brain networks (Sporns and Kötter, 2004). On the other hand, motifs are not commonly studied in undirected networks, due to the paucity of local undirected connectivity patterns. References: Milo et al. (2002) Science 298:824-827 Sporns O, Kötter R (2004) PLoS Biol 2: e369 Mika Rubinov, UNSW/U Cambridge, 2007-2015 ``` motif3funct_wei ``` MOTIF3FUNCT_WEI Intensity and coherence of functional class-3 motifs [I,Q,F] = motif3funct_wei(W); *Structural motifs* are patterns of local connectivity in complex networks. In contrast, *functional motifs* are all possible subsets of patterns of local connectivity embedded within structural motifs. Such patterns are particularly diverse in directed networks. The motif frequency of occurrence around an individual node is known as the motif fingerprint of that node. The motif intensity and coherence are weighted generalizations of the motif frequency. The motif intensity is equivalent to the geometric mean of weights of links comprising each motif. The motif coherence is equivalent to the ratio of geometric and arithmetic means of weights of links comprising each motif. Input: W, weighted directed connection matrix (all weights must be between 0 and 1) Output: I, node motif intensity fingerprint Q, node motif coherence fingerprint F, node motif frequency fingerprint Notes: 1. The function find_motif34.m outputs the motif legend. 2. Average intensity and coherence are given by I./F and Q./F 3. All weights must be between 0 and 1. This may be achieved using the weight_conversion.m function, as follows: W_nrm = weight_conversion(W, 'normalize'); 4. There is a source of possible confusion in motif terminology. Motifs ("structural" and "functional") are most frequently considered only in the context of anatomical brain networks (Sporns and Kötter, 2004). On the other hand, motifs are not commonly studied in undirected networks, due to the paucity of local undirected connectivity patterns. References: Onnela et al. (2005), Phys Rev E 71:065103 Milo et al. (2002) Science 298:824-827 Sporns O, Kötter R (2004) PLoS Biol 2: e369 Mika Rubinov, UNSW/U Cambridge, 2007-2015 ``` motif3struct_bin ``` MOTIF3STRUCT_BIN Frequency of structural class-3 motifs [f,F] = motif3struct_bin(A); Structural motifs are patterns of local connectivity in complex networks. Such patterns are particularly diverse in directed networks. The motif frequency of occurrence around an individual node is known as the motif fingerprint of that node. The total motif frequency of occurrence in the whole network is correspondingly known as the motif fingerprint of the network. Input: A, binary directed connection matrix Output: F, node motif frequency fingerprint f, network motif frequency fingerprint Note: The function find_motif34.m outputs the motif legend. References: Milo et al. (2002) Science 298:824-827 Sporns O, Kötter R (2004) PLoS Biol 2: e369 Mika Rubinov, UNSW/U Cambridge, 2007-2015 ``` motif3struct_wei ``` MOTIF3STRUCT_WEI Intensity and coherence of structural class-3 motifs [I,Q,F] = motif3struct_wei(W); Structural motifs are patterns of local connectivity in complex networks. Such patterns are particularly diverse in directed networks. The motif frequency of occurrence around an individual node is known as the motif fingerprint of that node. The motif intensity and coherence are weighted generalizations of the motif frequency. The motif intensity is equivalent to the geometric mean of weights of links comprising each motif. The motif coherence is equivalent to the ratio of geometric and arithmetic means of weights of links comprising each motif. Input: W, weighted directed connection matrix (all weights must be between 0 and 1) Output: I, node motif intensity fingerprint Q, node motif coherence fingerprint F, node motif frequency fingerprint Notes: 1. The function find_motif34.m outputs the motif legend. 2. Average intensity and coherence are given by I./F and Q./F 3. All weights must be between 0 and 1. This may be achieved using the weight_conversion.m function, as follows: W_nrm = weight_conversion(W, 'normalize'); References: Onnela et al. (2005), Phys Rev E 71:065103 Milo et al. (2002) Science 298:824-827 Sporns O, Kötter R (2004) PLoS Biol 2: e369% Mika Rubinov, UNSW/U Cambridge, 2007-2015 ``` motif4funct_bin ``` MOTIF4FUNCT_BIN Frequency of functional class-4 motifs [f,F] = motif4funct_bin(A); *Structural motifs* are patterns of local connectivity in complex networks. In contrast, *functional motifs* are all possible subsets of patterns of local connectivity embedded within structural motifs. Such patterns are particularly diverse in directed networks. The motif frequency of occurrence around an individual node is known as the motif fingerprint of that node. The total motif frequency of occurrence in the whole network is correspondingly known as the motif fingerprint of the network. Input: A, binary directed connection matrix Output: F, node motif frequency fingerprint f, network motif frequency fingerprint Notes: 1. The function find_motif34.m outputs the motif legend. 2. There is a source of possible confusion in motif terminology. Motifs ("structural" and "functional") are most frequently considered only in the context of anatomical brain networks (Sporns and Kötter, 2004). On the other hand, motifs are not commonly studied in undirected networks, due to the paucity of local undirected connectivity patterns. References: Milo et al. (2002) Science 298:824-827 Sporns O, Kötter R (2004) PLoS Biol 2: e369 Mika Rubinov, UNSW/U Cambridge, 2007-2015 ``` motif4funct_wei ``` MOTIF4FUNCT_WEI Intensity and coherence of functional class-4 motifs [I,Q,F] = motif4funct_wei(W); *Structural motifs* are patterns of local connectivity in complex networks. In contrast, *functional motifs* are all possible subsets of patterns of local connectivity embedded within structural motifs. Such patterns are particularly diverse in directed networks. The motif frequency of occurrence around an individual node is known as the motif fingerprint of that node. The motif intensity and coherence are weighted generalizations of the motif frequency. The motif intensity is equivalent to the geometric mean of weights of links comprising each motif. The motif coherence is equivalent to the ratio of geometric and arithmetic means of weights of links comprising each motif. Input: W, weighted directed connection matrix (all weights must be between 0 and 1) Output: I, node motif intensity fingerprint Q, node motif coherence fingerprint F, node motif frequency fingerprint Notes: 1. The function find_motif34.m outputs the motif legend. 2. Average intensity and coherence are given by I./F and Q./F 3. All weights must be between 0 and 1. This may be achieved using the weight_conversion.m function, as follows: W_nrm = weight_conversion(W, 'normalize'); 4. There is a source of possible confusion in motif terminology. Motifs ("structural" and "functional") are most frequently considered only in the context of anatomical brain networks (Sporns and Kötter, 2004). On the other hand, motifs are not commonly studied in undirected networks, due to the paucity of local undirected connectivity patterns. References: Onnela et al. (2005), Phys Rev E 71:065103 Milo et al. (2002) Science 298:824-827 Sporns O, Kötter R (2004) PLoS Biol 2: e369% Mika Rubinov, UNSW/U Cambridge, 2007-2015 ``` motif4struct_bin ``` MOTIF4STRUCT_BIN Frequency of structural class-4 motifs [f,F] = motif4struct_bin(A); Structural motifs are patterns of local connectivity in complex networks. Such patterns are particularly diverse in directed networks. The motif frequency of occurrence around an individual node is known as the motif fingerprint of that node. The total motif frequency of occurrence in the whole network is correspondingly known as the motif fingerprint of the network. Input: A, binary directed connection matrix Output: F, node motif frequency fingerprint f, network motif frequency fingerprint Note: The function find_motif34.m outputs the motif legend. References: Milo et al. (2002) Science 298:824-827 Sporns O, Kötter R (2004) PLoS Biol 2: e369 Mika Rubinov, UNSW/U Cambridge, 2007-2015 ``` motif4struct_wei ``` MOTIF4STRUCT_WEI Intensity and coherence of structural class-4 motifs [I,Q,F] = motif4struct_wei(W); Structural motifs are patterns of local connectivity in complex networks. Such patterns are particularly diverse in directed networks. The motif frequency of occurrence around an individual node is known as the motif fingerprint of that node. The motif intensity and coherence are weighted generalizations of the motif frequency. The motif intensity is equivalent to the geometric mean of weights of links comprising each motif. The motif coherence is equivalent to the ratio of geometric and arithmetic means of weights of links comprising each motif. Input: W, weighted directed connection matrix (all weights must be between 0 and 1) Output: I, node motif intensity fingerprint Q, node motif coherence fingerprint F, node motif frequency fingerprint Notes: 1. The function find_motif34.m outputs the motif legend. 2. Average intensity and coherence are given by I./F and Q./F 3. All weights must be between 0 and 1. This may be achieved using the weight_conversion.m function, as follows: W_nrm = weight_conversion(W, 'normalize'); References: Onnela et al. (2005), Phys Rev E 71:065103 Milo et al. (2002) Science 298:824-827 Sporns O, Kötter R (2004) PLoS Biol 2: e369% Mika Rubinov, UNSW/U Cambridge, 2007-2015 ``` null_model_dir_sign ``` NULL_MODEL_DIR_SIGN Directed random graphs with preserved weight, degree and strength distributions W0 = null_model_dir_sign(W); W0 = null_model_dir_sign(W,bin_swaps); W0 = null_model_dir_sign(W,bin_swaps,wei_freq); [W0 R] = null_model_dir_sign(W,bin_swaps,wei_freq); This function randomizes an directed network with positive and negative weights, while preserving the degree and strength distributions. This function calls randmio_dir_signed.m Inputs: W, Directed weighted connection matrix bin_swaps, Average number of swaps of each edge in binary randomization. bin_swap=5 is the default (each edge rewired 5 times) bin_swap=0 implies no binary randomization wei_freq, Frequency of weight sorting in weighted randomization wei_freq must be in the range of: 0 < wei_freq <= 1 wei_freq=1 implies that weights are sorted at each step (default in older [<2011] versions of MATLAB) wei_freq=0.1 implies that weights are sorted at each 10th step (faster, default in newer versions of MATLAB) Output: W0, Randomized weighted connection matrix R, Correlation coefficients between strength sequences of input and output connection matrices Notes: The value of bin_swaps is ignored when binary topology is fully connected (e.g. when the network has no negative weights). Randomization may be better (and execution time will be slower) for higher values of bin_swaps and wei_freq. Higher values of bin_swaps may enable a more random binary organization, and higher values of wei_freq may enable a more accurate conservation of strength sequences. R are the correlation coefficients between positive and negative in-strength and out-strength sequences of input and output connection matrices and are used to evaluate the accuracy with which strengths were preserved. Note that correlation coefficients may be a rough measure of strength-sequence accuracy and one could implement more formal tests (such as the Kolmogorov-Smirnov test) if desired. Example usage: %Create random directed weights matrix W=randn(100); %Compute one instance of null model (slow execution time): %bin_swaps=5, rewire each binary edge 5 times on average %wei_freq=1, sort all edges at every step tic; [W0_slow R_slow]=null_model_dir_sign(W,5,1); R_slow, toc R_slow = 0.9795 0.9724 0.9772 0.9773 Elapsed time is 3.485388 seconds. %Compute another instance of of null model (fast execution time): %bin_swaps=5, rewire each binary edge 5 times on average %wei_freq=0.1, sort all edges at every 10th step (10=1/0.1) tic; [W0_fast R_fast]=null_model_dir_sign(W,5,0.1); R_fast, toc R_fast = 0.9655 0.9652 0.9717 0.9804 Elapsed time is 0.763831 seconds. Reference: Rubinov and Sporns (2011) Neuroimage 56:2068-79 2011-2015, Mika Rubinov, U Cambridge ``` null_model_und_sign ``` NULL_MODEL_UND_SIGN Random graphs with preserved weight, degree and strength distributions W0 = null_model_und_sign(W); W0 = null_model_und_sign(W,bin_swaps); W0 = null_model_und_sign(W,bin_swaps,wei_freq); [W0 R] = null_model_und_sign(W,bin_swaps,wei_freq); This function randomizes an undirected network with positive and negative weights, while preserving the degree and strength distributions. This function calls randmio_und_signed.m Inputs: W, Undirected weighted connection matrix bin_swaps, Average number of swaps of each edge in binary randomization. bin_swap=5 is the default (each edge rewired 5 times) bin_swap=0 implies no binary randomization wei_freq, Frequency of weight sorting in weighted randomization wei_freq must be in the range of: 0 < wei_freq <= 1 wei_freq=1 implies that weights are resorted at each step (default in older [<2011] versions of MATLAB) wei_freq=0.1 implies that weights are sorted at each 10th step (faster, default in newer versions of Matlab) Output: W0, Randomized weighted connection matrix R, Correlation coefficient between strength sequences of input and output connection matrices Notes: The value of bin_swaps is ignored when binary topology is fully connected (e.g. when the network has no negative weights). Randomization may be better (and execution time will be slower) for higher values of bin_swaps and wei_freq. Higher values of bin_swaps may enable a more random binary organization, and higher values of wei_freq may enable a more accurate conservation of strength sequences. R are the correlation coefficients between positive and negative strength sequences of input and output connection matrices and are used to evaluate the accuracy with which strengths were preserved. Note that correlation coefficients may be a rough measure of strength-sequence accuracy and one could implement more formal tests (such as the Kolmogorov-Smirnov test) if desired. Example usage: %Create random weights matrix W=tril(randn(100),-1); W=W+W.'; %Compute one instance of null model (slow execution time): %bin_swaps=5, rewire each binary edge 5 times on average %wei_freq=1, sort all edges at every step tic; [W0_slow R_slow]=null_model_und_sign(W,5,1); R_slow, toc R_slow = 0.9720 0.9746 Elapsed time is 2.112715 seconds. %Compute another instance of of null model (fast execution time): %bin_swaps=5, rewire each binary edge 5 times on average %wei_freq=0.1, sort all edges at every 10th step (10=1/0.1) tic; [W0_fast R_fast]=null_model_und_sign(W,5,0.1); R_fast, toc R_fast = 0.9623 0.9789 Elapsed time is 0.584797 seconds. Reference: Rubinov and Sporns (2011) Neuroimage 56:2068-79 2011-2015, Mika Rubinov, U Cambridge ``` pagerank_centrality ``` PAGERANK_CENTRALITY PageRank centrality r = pagerank_centrality(A, d, falff) The PageRank centrality is a variant of eigenvector centrality. This function computes the PageRank centrality of each vertex in a graph. Formally, PageRank is defined as the stationary distribution achieved by instantiating a Markov chain on a graph. The PageRank centrality of a given vertex, then, is proportional to the number of steps (or amount of time) spent at that vertex as a result of such a process. The PageRank index gets modified by the addition of a damping factor, d. In terms of a Markov chain, the damping factor specifies the fraction of the time that a random walker will transition to one of its current state's neighbors. The remaining fraction of the time the walker is restarted at a random vertex. A common value for the damping factor is d = 0.85. Inputs: A, adjacency matrix d, damping factor falff, initial page rank probability (non-negative) Outputs: r, vectors of page rankings Note: The algorithm will work well for smaller matrices (number of nodes around 1000 or less) References: [1]. GeneRank: Using search engine technology for the analysis of microarray experiments, by Julie L. Morrison, Rainer Breitling, Desmond J. Higham and David R. Gilbert, BMC Bioinformatics, 6:233, 2005. [2]. Boldi P, Santini M, Vigna S (2009) PageRank: Functional dependencies. ACM Trans Inf Syst 27, 1-23. Xi-Nian Zuo, Institute of Psychology, Chinese Academy of Sciences, 2011 Rick Betzel, Indiana University, 2012 ``` participation_coef ``` PARTICIPATION_COEF Participation coefficient P = participation_coef(W,Ci); Participation coefficient is a measure of diversity of intermodular connections of individual nodes. Inputs: W, binary/weighted, directed/undirected connection matrix Ci, community affiliation vector flag, 0, undirected graph (default) 1, directed graph: out-degree 2, directed graph: in-degree Output: P, participation coefficient Reference: Guimera R, Amaral L. Nature (2005) 433:895-900. 2008-2015 Mika Rubinov, UNSW/U Cambridge Alex Fornito, University of Melbourne ``` participation_coef_sign ``` PARTICIPATION_COEF_SIGN Participation coefficient [Ppos Pneg] = participation_coef_sign(W,Ci); Participation coefficient is a measure of diversity of intermodular connections of individual nodes. Inputs: W, undirected connection matrix with positive and negative weights Ci, community affiliation vector Output: Ppos, participation coefficient from positive weights Pneg, participation coefficient from negative weights Reference: Guimera R, Amaral L. Nature (2005) 433:895-900. 2011, Mika Rubinov, UNSW ``` partition_distance ``` PARTITION_DISTANCE Distance or similarity between community partitions This function quantifies information-theoretic distance (normalized variation of information) or similarity (normalized mutual information) between community partitions. VIn = partition_distance(Cx); VIn = partition_distance(Cx, Cy); [VIn, MIn] = partition_distance(Cx, Cy); Inputs: Cx, Community partition vector or matrix of n rows and p columns, n is the number of network nodes, and p is the number of input community partitions (in the case of vector input p=1). Cy (optional argument), Community partition vector or matrix of n rows and q columns. n is the number of nodes (must be equal to the number of nodes in Cq) and q is the number of input community partitions (may be different to the number of nodes in Cq). This argument may be omitted, in which case, the partition distance is computed between all pairwise partitions of Cx. Outputs: VIn, Normalized variation of information ([p, q] matrix) MIn, Normalized mutual information ([p, q] matrix) Notes: Mathematical definitions. VIn = [H(X) + H(Y) - 2MI(X, Y)]/log(n) MIn = 2MI(X, Y) / [H(X) + H(Y)] where H is the entropy and MI is the mutual information Reference: Meila M (2007) J Multivar Anal 98, 873-895. 2011-2017, Mika Rubinov, UNSW, Janelia HHMI ``` path_transitivity ``` PATH_TRANSITIVITY Transitivity based on shortest paths T = path_transitivity(W,transform) This function computes the density of local detours (triangles) that are available along the shortest-paths between all pairs of nodes. Inputs: W, unweighted/weighted undirected connection *weight* OR *length* matrix. transform, If the input matrix is a connection *weight* matrix, specify a transform that map input connection weights to connection lengths. Two transforms are available. 'log' -> l_ij = -log(w_ij) 'inv' -> l_ij = 1/w_ij If the input matrix is a connection *length* matrix, do not specify a transform (or specify an empty transform argument). Output: T, matrix of pairwise path transitivity. Olaf Sporns, Andrea Avena-Koenigsberger and Joaquin Goñi, IU Bloomington, 2014 References: Goñi et al (2014) PNAS doi: 10.1073/pnas.131552911 ``` predict_fc ``` PREDICT_FC Prediction of functional connectivity from structural connectivity [FCpre,FCcorr,beta,pred_data,R] = predict_fc(SC,FC,ED,pred_var,model) [FCpre,FCcorr,beta] = predict_fc(SC,FC,[],{'SPLwei_log','SIwei_log'},'quadratic') This function computes regression coefficients to predict FC from structural-based measures that are used as predictor variables. Inputs: SC, Weighted/unweighted undirected NxN Structural Connectivity matrix. FC, Functional connections. FC can be a NxN symmetric matrix or an ((N*(N-1))/2) x 1 vector containing the upper triangular elements of the square FC matrix (excluding diagonal elements). ED, Euclidean distance matrix or upper triangular vector of the matrix (optional) pred_var, Set of M predictors. These can be given as an KxM array where K = ((N*(N-1))/2) and M is the number of predictors. Alternatively, pred_var can be a cell with the names of network measures to be used as predictors. Accepted network measure names are: SPLbin - Shortest-path length (binary) SPLwei_inv - Shortest-path length computed with an inv transform SPLwei_log - Shortest-path length computed with a log transform SPLdist - Shortest-path length computed with no transform SIbin - Search Information of binary shortest-paths SIwei_inv - Search Information of shortest-paths computed with an inv transform SIwei_log - Search Information of shortest-paths computed with a log transform SIdist - Search Information of shortest-paths computed with no transform T - Path Transitivity deltaMFPT - Column-wise z-scored mean first passage time neighOverlap - Neighborhood Overlap MI - Matching Index If no predictors are specified, the defaults are {'SPLwei_log', 'SIwei_log'}. model, Specifies the order of the regression model used within matlab's function regstats.m. 'model' can be any option accepted by matlab's regstats.m function (e.g. 'linear', 'interaction', 'quadratic', etc.) If no model is specified, 'linear' is the default. Output: FCpre, Predicted NxN Functional Connectivity matrix FCcorr, Pearson Correlation between PCpred and FC beta, Regression Coefficients pred_data, KxM array of predictors. R, Output from regstats.m (e.g. 'beta', 'yhat', 'rsquare', 'adjrsquare', 'tstat', 'r', 'mse', 'standres'). References: Goñi et al (2014) PNAS, 833–838, doi: 10.1073/pnas.1315529111 Andrea Avena-Koenigsberger and Joaquin Goñi, IU Bloomington, 2014 ``` randmio_dir ``` RANDMIO_DIR Random graph with preserved in/out degree distribution R = randmio_dir(W, ITER); [R eff] = randmio_dir(W, ITER); This function randomizes a directed network, while preserving the in- and out-degree distributions. In weighted networks, the function preserves the out-strength but not the in-strength distributions. Input: W, directed (binary/weighted) connection matrix ITER, rewiring parameter (each edge is rewired approximately ITER times) Output: R, randomized network eff, number of actual rewirings carried out References: Maslov and Sneppen (2002) Science 296:910 2007-2012 Mika Rubinov, UNSW Olaf Sporns, IU ``` randmio_dir_connected ``` RANDMIO_DIR_CONNECTED Random graph with preserved in/out degree distribution R = randmio_dir_connected(W, ITER); [R eff] = randmio_dir_connected(W, ITER); This function randomizes a directed network, while preserving the in- and out-degree distributions. In weighted networks, the function preserves the out-strength but not the in-strength distributions. The function also ensures that the randomized network maintains connectedness, the ability for every node to reach every other node in the network. The input network for this function must be connected. Input: W, directed (binary/weighted) connection matrix ITER, rewiring parameter (each edge is rewired approximately ITER times) Output: R, randomized network eff, number of actual rewirings carried out References: Maslov and Sneppen (2002) Science 296:910 2007-2012 Mika Rubinov, UNSW Olaf Sporns, IU ``` randmio_dir_signed ``` RANDMIO_DIR_SIGNED Random directed graph with preserved signed in/out degree distribution R = randmio_dir(W, ITER); [R,eff] = randmio_dir(W, ITER); This function randomizes a directed weighted network with positively and negatively signed connections, while preserving the positively and negatively signed in- and out-degree distributions. In weighted networks, the function preserves the out-strength but not the in-strength distributions. Input: W, directed (binary/weighted) connection matrix ITER, rewiring parameter (each edge is rewired approximately ITER times) Output: R, randomized network eff, number of actual rewirings carried out Reference: Maslov and Sneppen (2002) Science 296:910 2011-2015 Dani Bassett, UCSB Olaf Sporns, Indiana U Mika Rubinov, U Cambridge ``` randmio_und ``` RANDMIO_UND Random graph with preserved degree distribution R = randmio_und(W,ITER); [R eff]=randmio_und(W, ITER); This function randomizes an undirected network, while preserving the degree distribution. The function does not preserve the strength distribution in weighted networks. Input: W, undirected (binary/weighted) connection matrix ITER, rewiring parameter (each edge is rewired approximately ITER times) Output: R, randomized network eff, number of actual rewirings carried out References: Maslov and Sneppen (2002) Science 296:910 2007-2012 Mika Rubinov, UNSW Jonathan Power, WUSTL Olaf Sporns, IU ``` randmio_und_connected ``` RANDMIO_UND_CONNECTED Random graph with preserved degree distribution R = randmio_und_connected(W,ITER); [R eff] = randmio_und_connected(W, ITER); This function randomizes an undirected network, while preserving the degree distribution. The function does not preserve the strength distribution in weighted networks. The function also ensures that the randomized network maintains connectedness, the ability for every node to reach every other node in the network. The input network for this function must be connected. Input: W, undirected (binary/weighted) connection matrix ITER, rewiring parameter (each edge is rewired approximately ITER times) Output: R, randomized network eff, number of actual rewirings carried out References: Maslov and Sneppen (2002) Science 296:910 2007-2012 Mika Rubinov, UNSW Jonathan Power, WUSTL Olaf Sporns, IU ``` randmio_und_signed ``` RANDMIO_UND_SIGNED Random graph with preserved signed degree distribution R = randmio_und_signed(W,ITER); [R,eff] = randmio_und_signed(W,ITER); This function randomizes an undirected network with positively and negatively signed connections, while preserving the positively and negatively signed degree distribution. The function does not preserve the strength distribution in weighted networks. Input: W, undirected (binary/weighted) connection matrix ITER, rewiring parameter (each edge is rewired approximately ITER times) Output: R, randomized network eff, number of actual rewirings carried out Reference: Maslov and Sneppen (2002) Science 296:910 2011-2015 Dani Bassett, UCSB Olaf Sporns, Indiana U Mika Rubinov, U Cambridge ``` randomize_graph_partial_und ``` RANDOMIZE_GRAPH_PARTIAL_UND Swap edges with preserved degree sequence A = RANDOMIZE_GRAPH_PARTIAL_UND(A,B,MAXSWAP) takes adjacency matrices A and B and attempts to randomize matrix A by performing MAXSWAP rewirings. The rewirings will avoid any spots where matrix B is nonzero. Inputs: A, undirected adjacency matrix B, edges to avoid MAXSWAP, number of rewirings Outputs: A, randomized matrix Richard Betzel, Indiana University, 2013 Notes: 1. Based on the script randmio_und.m. 2. Graph may become disconnected as a result of rewiring. Always important to check. 3. A can be weighted, though the weighted degree sequence will not be preserved. ``` randomizer_bin_und ``` RANDOMIZER_BIN_UND Random graph with preserved in/out degree distribution R = randomizer_bin_und(A,alpha); This function randomizes a binary undirected network, while preserving the degree distribution. The function directly searches for rewirable edge pairs (rather than trying to rewire edge pairs at random), and hence avoids long loops and works especially well in dense matrices. Inputs: A, binary undirected connection matrix alpha, fraction of edges to rewire Outputs: R, randomized network References: Maslov and Sneppen (2002) Science 296:910 Jonathan Power, WUSTL. 3/1/10. ``` reachdist ``` REACHDIST Reachability and distance matrices [R,D] = reachdist(CIJ); The binary reachability matrix describes reachability between all pairs of nodes. An entry (u,v)=1 means that there exists a path from node u to node v; alternatively (u,v)=0. The distance matrix contains lengths of shortest paths between all pairs of nodes. An entry (u,v) represents the length of shortest path from node u to node v. The average shortest path length is the characteristic path length of the network. Input: CIJ, binary (directed/undirected) connection matrix Outputs: R, reachability matrix D, distance matrix Note: faster but more memory intensive than "breadthdist.m". Algorithm: algebraic path count. Olaf Sporns, Indiana University, 2002/2007/2008 ``` rentian_scaling_2d ``` RENTIAN_SCALING_2D Rentian scaling for networks embedded in two dimensions. [N,E] = rentian_scaling_2d(A,XY,n,tol) Physical Rentian scaling (or more simply Rentian scaling) is a property of systems that are cost-efficiently embedded into physical space. It is what is called a "topo-physical" property because it combines information regarding the topological organization of the graph with information about the physical placement of connections. Rentian scaling is present in very large scale integrated circuits, the C. elegans neuronal network, and morphometric and diffusion-based graphs of human anatomical networks. Rentian scaling is determined by partitioning the system into cubes, counting the number of nodes inside of each cube (N), and the number of edges traversing the boundary of each cube (E). If the system displays Rentian scaling, these two variables N and E will scale with one another in loglog space. The Rent's exponent is given by the slope of log10(E) vs. log10(N), and can be reported alone or can be compared to the theoretical minimum Rent's exponent to determine how cost efficiently the network has been embedded into physical space. Note: if a system displays Rentian scaling, it does not automatically mean that the system is cost-efficiently embedded (although it does suggest that). Validation occurs when comparing to the theoretical minimum Rent's exponent for that system. INPUTS: A: MxM adjacency matrix. Must be unweighted, binary, and symmetric. XY: Matrix of node placement coordinates. Must be in the form of an Mx2 matrix [x y], where M is the number of nodes and x and y are column vectors of node coordinates. n: Number of partitions to compute. Each partition is a data point. You want a large enough number to adequately estimate the Rent's exponent. tol: This should be a small value (for example 1e-6). In order to mitigate the effects of boundary conditions due to the finite size of the network, we only allow partitions that are contained within the boundary of the network. This is achieved by first computing the volume of the convex hull of the node coordinates (V). We then ensure that the volume of the convex hull computed on the original node coordinates plus the coordinates of the randomly generated partition (Vnew) is within a given tolerance of the original (i.e. check abs(V - Vnew) < tol). Thus tol, should be a small value in order to make sure the partitions are contained largely within the boundary of the network, and thus the number of nodes and edges within the box are not skewed by finite size effects. OUTPUTS: N: nx1 vector of the number of nodes in each of the n partitions. E: nx1 vector of the number of edges crossing the boundary of each partition. Subsequent Analysis: Rentian scaling plots are created by: figure; loglog(E,N,'*'); To determine the Rent's exponent, p, we need to determine the slope of E vs. N in loglog space, which is the Rent's exponent. There are many ways of doing this with more or less statistical rigor. Robustfit in MATLAB is one such option: [b,stats] = robustfit(log10(N),log10(E)) Then the Rent's exponent is b(1,2) and the standard error of the estimation is given by stats.se(1,2). Note: n=5000 was used in Bassett et al. 2010 in PLoS CB. Reference: Danielle S. Bassett, Daniel L. Greenfield, Andreas Meyer-Lindenberg, Daniel R. Weinberger, Simon W. Moore, Edward T. Bullmore. Efficient physical embedding of topologically complex information processing networks in brains and computer circuits. PLoS Comput Biol, 2010, 6(4):e1000748. Modification History: 2010: Original (Dani Bassett) Dec 2016: Updated code so that both partition centers and partition sizes are chosen at random. Also added in a constraint on partition placement that prevents boxes from being located outside the edges of the network. This helps prevent skewed results due to boundary effects arising from the finite size of the network. (Lia Papadopoulos) ``` rentian_scaling_3d ``` RENTIAN_SCALING_3D Rentian scaling for networks embedded in three dimensions. [N,E] = rentian_scaling_3d(A,XYZ,n,tol) Physical Rentian scaling (or more simply Rentian scaling) is a property of systems that are cost-efficiently embedded into physical space. It is what is called a "topo-physical" property because it combines information regarding the topological organization of the graph with information about the physical placement of connections. Rentian scaling is present in very large scale integrated circuits, the C. elegans neuronal network, and morphometric and diffusion-based graphs of human anatomical networks. Rentian scaling is determined by partitioning the system into cubes, counting the number of nodes inside of each cube (N), and the number of edges traversing the boundary of each cube (E). If the system displays Rentian scaling, these two variables N and E will scale with one another in loglog space. The Rent's exponent is given by the slope of log10(E) vs. log10(N), and can be reported alone or can be compared to the theoretical minimum Rent's exponent to determine how cost efficiently the network has been embedded into physical space. Note: if a system displays Rentian scaling, it does not automatically mean that the system is cost-efficiently embedded (although it does suggest that). Validation occurs when comparing to the theoretical minimum Rent's exponent for that system. INPUTS: A: MxM adjacency matrix. Must be unweighted, binary, and symmetric. XYZ: Matrix of node placement coordinates. Must be in the form of an Mx3 matrix [x y z], where M is the number of nodes and x, y, z are column vectors of node coordinates. n: Number of partitions to compute. Each partition is a data point. You want a large enough number to adequately estimate the Rent's exponent. tol: This should be a small value (for example 1e-6). In order to mitigate the effects of boundary conditions due to the finite size of the network, we only allow partitions that are contained within the boundary of the network. This is achieved by first computing the volume of the convex hull of the node coordinates (V). We then ensure that the volume of the convex hull computed on the original node coordinates plus the coordinates of the randomly generated partition (Vnew) is within a given tolerance of the original (i.e. check abs(V - Vnew) < tol). Thus tol, should be a small value in order to make sure the partitions are contained largely within the boundary of the network, and thus the number of nodes and edges within the box are not skewed by finite size effects. OUTPUTS: N: nx1 vector of the number of nodes in each of the n partitions. E: nx1 vector of the number of edges crossing the boundary of each partition. Subsequent Analysis: Rentian scaling plots are created by: figure; loglog(E,N,'*'); To determine the Rent's exponent, p, we need to determine the slope of E vs. N in loglog space, which is the Rent's exponent. There are many ways of doing this with more or less statistical rigor. Robustfit in MATLAB is one such option: [b,stats] = robustfit(log10(N),log10(E)) Then the Rent's exponent is b(1,2) and the standard error of the estimation is given by stats.se(1,2). Note: n=5000 was used in Bassett et al. 2010 in PLoS CB. Reference: Danielle S. Bassett, Daniel L. Greenfield, Andreas Meyer-Lindenberg, Daniel R. Weinberger, Simon W. Moore, Edward T. Bullmore. Efficient physical embedding of topologically complex information processing networks in brains and computer circuits. PLoS Comput Biol, 2010, 6(4):e1000748. Modification History: 2010: Original (Dani Bassett) Dec 2016: Updated code so that both partition centers and partition sizes are chosen at random. Also added in a constraint on partition placement that prevents boxes from being located outside the edges of the network. This helps prevent skewed results due to boundary effects arising from the finite size of the network. (Lia Papadopoulos) ``` reorderMAT ``` REORDERMAT Reorder matrix for visualization [MATreordered,MATindices,MATcost] = reorderMAT(MAT,H,cost); This function reorders the connectivity matrix in order to place more edges closer to the diagonal. This often helps in displaying community structure, clusters, etc. Inputs: MAT, connection matrix H, number of reordering attempts cost, 'line' or 'circ', for shape of lattice (linear or ring lattice) MATreordered reordered connection matrix MATindices reordered indices MATcost cost of reordered matrix Olaf Sporns, Indiana University ``` reorder_matrix ``` REORDER_MATRIX Matrix reordering for visualization [Mreordered,Mindices,cost] = reorder_matrix(M1,cost,flag) This function rearranges the nodes in matrix M1 such that the matrix elements are squeezed along the main diagonal. The function uses a version of simulated annealing. Inputs: M1 = connection matrix (weighted or binary, directed or undirected) cost = 'line' or 'circ', for shape of lattice cost (linear or ring lattice) Mreordered = reordered connection matrix Mindices = reordered indices cost = distance between M1 and Mreordered Note that in general, the outcome will depend on the initial condition (the setting of the random number seed). Also, there is no good way to determine optimal annealing parameters in advance - these paramters will need to be adjusted "by hand" (particularly H, Texp, and T0). For large and/or dense matrices, it is highly recommended to perform exploratory runs varying the settings of 'H' and 'Texp' and then select the best values. Based on extensive testing, it appears that T0 and Hbrk can remain unchanged in most cases. Texp may be varied from 1-1/H to 1-10/H, for example. H is the most important parameter - set to larger values as the problem size increases. It is advisable to run this function multiple times and select the solution(s) with the lowest 'cost'. Setting 'Texp' to zero cancels annealing and uses a greedy algorithm instead. Yusuke Adachi, University of Tokyo 2010 Olaf Sporns, Indiana University 2010 ``` reorder_mod ``` REORDER_MOD Reorder connectivity matrix by modular structure On = reorder_mod(W,M); [On Wr] = reorder_mod(W,M); This function reorders the connectivity matrix by modular structure and may consequently be useful in visualization of modular structure. Inputs: W, connectivity matrix (binary/weighted undirected/directed) M, module affiliation vector Outputs: On, new node order Wr, reordered connectivity matrix Used in: Rubinov and Sporns (2011) NeuroImage; Zingg et al. (2014) Cell. 2011, Mika Rubinov, UNSW/U Cambridge ``` resource_efficiency_bin ``` RESOURCE_EFFICIENCY_BIN Resource efficiency and shortest-path probability [Eres,prob_SPL] = resource_efficiency_bin(adj,lambda,SPL,M) The resource efficiency between nodes i and j is inversly proportional to the amount of resources (i.e. number of particles or messages) required to ensure with probability 0 < lambda < 1 that at least one of them will arrive at node j in exactly SPL steps, where SPL is the length of the shortest-path between i and j. The shortest-path probability between nodes i and j is the probability that a single random walker starting at node i will arrive at node j by following (one of) the shortest path(s). Inputs: adj, Unweighted, undirected adjacency matrix lambda, Probability (0 < lambda < 1) set to NAN if computation of Eres is not desired SPL, Shortest-path length matrix (optional) M, Transition probability matrix (optional) Outputs: Eres, Resource efficiency matrix. prob_SPL, Shortest-path probability matrix Notes: Global measures for both Eres and prob_SPL are well defined and can be computed as the average across the off-diagonal elements: GEres = mean(Eres(~eye(N)>0)); Gprob_SPL = mean(rob_SPL(~eye(N)>0)); Reference: Goñi J, et al (2013) PLoS ONE Joaquin Goñi, IU Bloomington, 2012 ``` retrieve_shortest_path ``` RETRIEVE_SHORTEST_PATH Retrieval of shortest path This function finds the sequence of nodes that comprise the shortest path between a given source and target node. Inputs: s, Source node: i.e. node where the shortest path begins. t, Target node: i.e. node where the shortest path ends. hops, Number of edges in the path. This matrix may be obtained as the second output argument of the function "distance_wei_floyd.m". Pmat, Pmat is a matrix whose elements {k,t} indicate the next node in the shortest path between k and t. This matrix may be obtained as the third output of the function "distance_wei_floyd.m" Output: path, Nodes comprising the shortest path between nodes s and t. Andrea Avena-Koenigsberger and Joaquin Goñi, IU, 2012 ``` rich_club_bd ``` RICH_CLUB_BD Rich club coefficients (binary directed graph) R = rich_club_bd(CIJ) [R,Nk,Ek] = rich_club_bd(CIJ,klevel) The rich club coefficient, R, at level k is the fraction of edges that connect nodes of degree k or higher out of the maximum number of edges that such nodes might share. Input: CIJ, connection matrix, binary and directed klevel, optional input argument. klevel sets the maximum level at which the rich club coefficient will be calculated. If klevel is not included the the maximum level will be set to the maximum degree of CIJ. Output: R, vector of rich-club coefficients for levels 1 to klevel. Nk, number of nodes with degree>k Ek, number of edges remaining in subgraph with degree>k Reference: Colizza et al. (2006) Nat. Phys. 2:110. Martijn van den Heuvel, University Medical Center Utrecht, 2011 ``` rich_club_bu ``` RICH_CLUB_BU Rich club coefficients (binary undirected graph) R = rich_club_bu(CIJ) [R,Nk,Ek] = rich_club_bu(CIJ,klevel) The rich club coefficient, R, at level k is the fraction of edges that connect nodes of degree k or higher out of the maximum number of edges that such nodes might share. Input: CIJ, connection matrix, binary and undirected klevel, optional input argument. klevel sets the maximum level at which the rich club coefficient will be calculated. If klevel is not included the the maximum level will be set to the maximum degree of CIJ. Output: R, vector of rich-club coefficients for levels 1 to klevel. Nk, number of nodes with degree>k Ek, number of edges remaining in subgraph with degree>k Reference: Colizza et al. (2006) Nat. Phys. 2:110. Martijn van den Heuvel, University Medical Center Utrecht, 2011 ``` rich_club_wd ``` RICH_CLUB_WD Rich club coefficients curve (weighted directed graph) Rw = rich_club_wd(CIJ,varargin) The weighted rich club coefficient, Rw, at level k is the fraction of edge weights that connect nodes of degree k or higher out of the maximum edge weights that such nodes might share. Inputs: CIJ: weighted directed connection matrix k-level: (optional) max level of RC(k). (by default k-level quals the maximal degree of CIJ) Output: Rw: rich-club curve References: T Opsahl et al. Phys Rev Lett, 2008, 101(16) M van den Heuvel, O Sporns, J Neurosci 2011 31(44) Martijn van den Heuvel, University Medical Center Utrecht, 2011 ``` rich_club_wu ``` RICH_CLUB_WU Rich club coefficients curve (weighted undirected graph) Rw = rich_club_wu(CIJ,varargin) % rich club curve for weighted graph The weighted rich club coefficient, Rw, at level k is the fraction of edge weights that connect nodes of degree k or higher out of the maximum edge weights that such nodes might share. Inputs: CIJ: weighted directed connection matrix k-level: (optional) max level of RC(k). (by default k-level quals the maximal degree of CIJ) Output: Rw: rich-club curve References: T Opsahl et al. Phys Rev Lett, 2008, 101(16) M van den Heuvel, O Sporns, J Neurosci 2011 31(44) Martijn van den Heuvel, University Medical Center Utrecht, 2011 ``` rout_efficiency ``` ROUT_EFFICIENCY Mean, pair-wise and local routing efficiency [GErout,Erout,Eloc] = rout_efficiency(D,transform); The routing efficiency is the average of inverse shortest path length. The local routing efficiency of a node u is the routing efficiency computed on the subgraph formed by the neighborhood of node u (excluding node u). Inputs: D, Weighted/unweighted directed/undirected connection *weight* OR *length* matrix. transform, If the input matrix is a connection *weight* matrix, specify a transform that map input connection weights to connection lengths. Two transforms are available. 'log' -> l_ij = -log(w_ij) 'inv' -> l_ij = 1/w_ij If the input matrix is a connection *length* matrix, do not specify a transform (or specify an empty transform argument). Outputs: GErout, Mean global routing efficiency (scalar). Erout, Pair-wise routing efficiency (matrix). Eloc, Local efficiency (vector) Note: The input matrix may be either a connection weight matrix, or a connection length matrix. The connection length matrix is typically obtained with a mapping from weight to length, such that higher weights are mapped to shorter lengths (see above). Algorithm: Floyd–Warshall Algorithm References: Latora and Marchiori (2001) Phys Rev Lett Goñi et al (2013) PLoS ONE Avena-Koenigsberger et al (2016) Brain Structure and Function Andrea Avena-Koenigsberger and Joaquin Goñi, IU Bloomington, 2012 ``` score_wu ``` SCORE_WU S-score [CIJscore,sn] = score_wu(CIJ,s); The s-core is the largest subnetwork comprising nodes of strength at least s. This function computes the s-core for a given weighted undirected connection matrix. Computation is analogous to the more widely used k-core, but is based on node strengths instead of node degrees. input: CIJ, connection/adjacency matrix (weighted, undirected) s, level of s-core. Note: s can take on any fractional value output: CIJscore, connection matrix of the s-core. This matrix contains only nodes with a strength of at least s. sn, size of s-score Olaf Sporns, Indiana University, 2007/2008/2010/2012 ``` search_information ``` SEARCH_INFORMATION Search information SI = search_information(adj,transform,has_memory) Computes the amount of information (measured in bits) that a random walker needs to follow the shortest path between a given pair of nodes. Inputs: adj, Weighted/unweighted directed/undirected connection *weight* OR *length* matrix. transform, If the input matrix is a connection *weight* matrix, specify a transform that map input connection weights to connection lengths. Two transforms are available. 'log' -> l_ij = -log(w_ij) 'inv' -> l_ij = 1/w_ij If the input matrix is a connection *length* matrix, do not specify a transform (or specify an empty transform argument). has_memory, This flag defines whether or not the random walker "remembers" its previous step, which has the effect of reducing the amount of information needed to find the next state. If this flag is not set, the walker has no memory by default. Outputs: SI, pair-wise search information (matrix). Note that SI(i,j) may be different from SI(j,i), hense, SI is not a symmetric matrix even when adj is symmetric. References: Rosvall et al. (2005) Phys Rev Lett 94, 028701 Goñi et al (2014) PNAS doi: 10.1073/pnas.131552911 Andrea Avena-Koenigsberger and Joaquin Goñi, IU Bloomington, 2014 ``` strengths_dir ``` STRENGTHS_DIR In-strength and out-strength [is,os,str] = strengths_dir(CIJ); Node strength is the sum of weights of links connected to the node. The instrength is the sum of inward link weights and the outstrength is the sum of outward link weights. Input: CIJ, directed weighted connection matrix Output: is, node instrength os, node outstrength str, node strength (instrength + outstrength) Notes: Inputs are assumed to be on the columns of the CIJ matrix. Olaf Sporns, Indiana University, 2002/2006/2008 ``` strengths_und ``` STRENGTHS_UND Strength str = strengths_und(CIJ); Node strength is the sum of weights of links connected to the node. Input: CIJ, undirected weighted connection matrix Output: str, node strength Olaf Sporns, Indiana University, 2002/2006/2008 ``` strengths_und_sign ``` STRENGTHS_UND_SIGN Strength and weight [Spos Sneg] = strengths_und_sign(W); [Spos Sneg vpos vneg] = strengths_und_sign(W); Node strength is the sum of weights of links connected to the node. Inputs: W, undirected connection matrix with positive and negative weights Output: Spos/Sneg, nodal strength of positive/negative weights vpos/vneg, total positive/negative weight 2011, Mika Rubinov, UNSW ``` subgraph_centrality ``` SUBGRAPH_CENTRALITY Subgraph centrality of a network Cs = subgraph_centrality(CIJ) The subgraph centrality of a node is a weighted sum of closed walks of different lengths in the network starting and ending at the node. This function returns a vector of subgraph centralities for each node of the network. Inputs: CIJ, adjacency matrix (binary) Outputs: Cs, subgraph centrality Reference: Estrada and Rodriguez-Velasquez (2005) Phys Rev E 71, 056103 Estrada and Higham (2010) SIAM Rev 52, 696. Xi-Nian Zuo, Chinese Academy of Sciences, 2010 Rick Betzel, Indiana University, 2012 ``` threshold_absolute ``` THRESHOLD_ABSOLUTE Absolute thresholding W_thr = threshold_absolute(W, thr); This function thresholds the connectivity matrix by absolute weight magnitude. All weights below the given threshold, and all weights on the main diagonal (self-self connections) are set to 0. Inputs: W weighted or binary connectivity matrix thr weight treshold Output: W_thr thresholded connectivity matrix Mika Rubinov, UNSW, 2009-2010 ``` threshold_proportional ``` THRESHOLD_PROPORTIONAL Proportional thresholding W_thr = threshold_proportional(W, p); This function "thresholds" the connectivity matrix by preserving a proportion p (0