This repository contains all the code necessary to reproduce the analyses in this paper. In this paper, we explore the mesoscale organization of functional networks supporting brain computer interface learning.

• Required matlab toolboxes: eeglab, fieldtrip
• Required python toolboxes: os, numpy, scipy, sys, (NMF from https://github.com/akhambhati/Echobase/tree/master/Echobase/Network/Partitioning/Subgraph)
• Required R toolboxes: R.matlab, ggplot2, dplyr, coin, lmPer, car, aplpack, RColorBrewer, wesanderson, ez, plyr

Data is very clean, so minimal preprocessing is needed

• Run preproc_behav.m to get summary statistics for each subjects performance and learning rate (slope)
• Run preproc_wpli and preproc_wpli_pr.m to make FC data and phase randomized FC data.
  • This will call the function wrapper_wpli.m and wrapper_pr_wpli.m. This script will combine planar gradiometers using a helper function from fieldtrip, calculate the wPLI using custom helper function get_window_wpli.m for each trial, and save a matrix for each session. The function outputs any errors in the pipeline

• Run format_NMF_wpli.m. This saves a single connection matrix of FC and performance for each subject/band.
• For every subject/band combination, run NMF_pipe_wpli.py and NMF_pipe_wpli_pr.py. This was done on the cluster, where each dataset was a job. The performs the NMF decomposition, and saves the selected parameters (alpha, beta, and m), the error, as well as the loadings and temporal coefficients

• Run get_noise_sg.m. This will create an index of subgraph that are not regularized for empirical and UPR data. This index will be called in future scripts

• Run behavior.R. This script essentially just plots data loaded from a mat file. It also does a repeated measures anova on the subjects average performance

• Run exp_beh_wpli.m. All this script does is reformat some of the data so that it can be analyzed in R. It combines performance data with ranked loading data
• Run exp_beh_corr_wpli.R. This will generate figures and statistics for fig 2

• Run consensus_subgraphs.m. This calls thresh.mat. A helper function for thresholding matrices. This will create the consensus networks, where each edge is the number of subjects that have a given edge in the X% strongest connections of their subgraph. This also makes some files for plotting. It also saves the average consistency for each graph
  • Plots are made in gephi. You will need to run make_gml.py to make files for gephi
• Run consistent_subgraphs_pr.m. to repeat this process for null data.
• Run compare_consensistency.m. to compare consistency between null and empirical data
• Run consistent_edge_location.m. to get lobes of the most consistent edges
• Run consistent_check_norm.R. to get quantile-quantile plots of consistency

• Run temp_exp.m. This reformats data about the energy, and peak of temporal expression across subgraphs for analysis in R
• Run temp_exp.R. This makes the plot and statistics for figure 4.

• Run get_oc_parameters.m This will find parameters that will give low error
• Run optimal_control_MI.m This generates data for most selectivity analyses, as well as positive and negative controls
• Run density_slope.m to get density data. This is used to see if effects are above and beyond what is expected from density
• Run optimal_control_state_variations.m and optimal_control_magnitude.m This generates data for all the state specific controls
• Run optimal_control_alt_params.m This generates data for different parameter sets
• Run opt_energy.R to generate all plots and statistics in the paper

• NMF parameter distributions and comparisons: run NMF_parameters.ipynb
• Comparison of the Number of Triangles:
  • Run consistency_thresh.m to explore consistency at different thresholds
• Model Validation:
  • Run count_triangles.m
  • Run compare_triangles.R
  • Run make_states.m, which calls get_mtmfft_power.m to get states (z-scored power across channels)
  • Run model_vaidation.m
  • Run model_validation_all_sg.m
  • Run model_validation_stats.R

• Singular matrix: there are a couple sources for this. (1) your original data is not full rank (2) you data is transposed, make sure what you submit is nWin x nCon, the code does not check (3) Subgraph partitioning creates rank deficient matrices, lowering the max rank will help with this

• Jennifer Stiso