Modularity - michaelandric/steadystate GitHub Wiki

This is to assess network structures or graphs. Graphs are sets of nodes and edges. Here, voxels are "nodes" and the connection between them is an "edge." So, we can represent and analyze the brain network as a graph. The analysis measures the network's partitioning into modules - or communities. A community is a densely connected set of nodes.

Steps

I do this analysis in R with the bct package. Typically, I invoke the R code on the command line via a script that sets variables for the code. The exception is when I run blondel_community to get the actual modularity value. That is done with RInvoke_blondel.py. See below (Measure network modularity) for why.

1. Get the time series in text files.

Can use a simple AFNI maskdump

2. Correlate every voxel's time series against every other voxel's

Do this in R with fcorr. The fcorr script will first call AFNI 3dDeconvolve to 'clean' the time series - e.g., regress out motion effects. That 'clean' time series is then dumped from AFNI format to a text file. This text file is then used in the fcorr function: A full N x N correlation matrix is built.

3. Threshold the N x N correlation matrix

Using this threshold script to threshold (for now) the Pearson's correlations at 0.5. This parameter can be easily iterated. The limit right now is space; the thresholded matrix is the same size as the original correlation matrix. For a ~24,000^2 matrix, this is about 2 GB. Note: The threshold will change modularity. A higher threshold gives fewer connections between voxels. This will increase modularity - instead of global, connections become more "local" (in the sense of a community) as the connectivity threshold increases.

4. Convert matrix format to use with the "Blondel" (aka "Louvain") algorithm.

Done with this simple code that calls the blondel_convert function.

5. Measure network modularity

Done with the blondel_community function. Instead of RInvoke.sh, this is called on the command line with RInvoke_blondel.py. The R script writes out a "tree" file. This contains the index (voxel number) in the first column. The second column gives - with a numeric identifier - which community that voxel belongs. All hierarchical levels are in the same file. This file is parsed in the next step to give a tree for each hierarchical level.

Why this is done with a special RInvoke (RInvoke_blondel.py):
When run via the command line, an .Rout file is generated. Currently, 'blondel_community' prints out - in the R terminal - the modularity value and increase in modularity at each hierarchical level. The only value that can be written out by blondel_community is the final modularity value at the highest hierarchical level. So, the .Rout file contains information needed in the next step: the number of hierarchical levels in the tree. This is why RInvoke_blondel.py takes the blondel script and copies it into the data directory - with a name that's particular to the condition. Otherwise, the .Rout file would be overwritten each time it's run and lose the needed information.

6. Print out the community structure at each hierarchical level

Done with the blondel_hierarchy function. This will return a "tree" file for each hierarchical level. These "tree" files contain two columns. The first column is the index, i.e., the voxel number, 1 ... n where n is the total number of voxels. The second column is the numeric community label. The community numbers themselves are only identifiers. They do not index any feature other than their membership.

7. Analyze group-level modularity effect across conditions

This is a repeated measures ANOVA to determine whether modularity changes as a function of condition. I use this R code. Currently, this code prints out the stats - along with pairwise comparisons - to the R terminal. So, I run it via the command line and keep the .Rout file as record.