Code that generated the results published in (please, cite)

 * Girardi-Schappo, M., Fadaie, F., Lee, H.M., Caldairou, B., Sziklas, V., Crane, J.,
   Bernhardt, B.C., Bernasconi, A. and Bernasconi, N. (2021),
   Altered communication dynamics reflect cognitive deficits in temporal lobe epilepsy.
   Epilepsia, 62: 1022-1033. https://doi.org/10.1111/epi.16864



Spreading times source code
---------------------------

by Mauricio Girardi Schappo, PhD -- 2016-08-01 - 2017-08-24


This file contains: the structure of the program, instructions on how to run the simulation, and an abstract example of how to do so...

I) Structure of the program:

The program has 4 main procedures:
    0) Flipping RTLE adjacency matrices
    1) Simulation
    2) Output data processing
    3) Plotting

External dependencies:
    0) Flipping RTLE adjacency matrices
        a) aal data structure: /data/noel/noel7/mauricio/data/AAL_data/aal_cortex_map_olf294_fix.mat
        b) patient info table: /data/noel/noel7/mauricio/data/TLE_NB_BB.mat (provided by Neda)
    1) Simulation
        No external dependencies
    2) Output data processing
        a) aal data structure: /data/noel/noel7/mauricio/data/AAL_data/aal_cortex_map_olf294_fix.mat
        b) patient info table: /data/noel/noel7/mauricio/data/TLE_NB_BB.mat (provided by Neda)
        c) community files generated by ./func/calcCommunities.m for each case
        d) shortest path files generated by ./func/calcShortestPath.m for each case
    3) Plotting
        a) aal data structure: /data/noel/noel7/mauricio/data/AAL_data/aal_cortex_map_olf294_fix.mat

The source code is organized in three directories:
    0) . (this) -> this directory is reserved for the "main" functions of each of the procedures above
    1) ./src    -> simulation source code, compiling scripts and mex files
    2) ./func   -> miscellaneous functions that are used by the many steps of the program procedures listed above

Important files for each procedure:
    0) Flipping RTLE adjacency matrices
        ./flip_AAL_matrices.m                 -> main flipping file
        ./func/flipAALMatrix_rowsPerCols.m    -> main subject matrix flipping algorithm
        ./func/getContraLateralLabel.m        -> gets the numeric indices of the contra-lateral parcels that correspond to the provided parcel indices

    1) Simulation
        ./RunPhaseTransitionSingleMatrix.m  -> Runs the simulation for many activation thresholds (theta) for each input case;
                                                Output data: quantity of activated nodes versus theta for each input case
        ./RunSpreadingTimeSingleMatrix.m    -> Runs the simulation for 1-seed signal spreading, starting from each node, to every other nodes for each input case
                                                 Output data: spreading time matrix for each case
        ./RunTwoSeedsSpreading.m            -> Runs 2-seed spreading simulation for every seed combination
                                                 Output data: entropy, conformity, competitiveness, maximum activation time; all for each input case
        ./RunNetworkActivationTime.m        -> Runs signal spreading simulation for the given input configuration of seeds
                                                 Output data: network activation time for each case and the given seed configuration

    2) Output data processing
        ./filter_parse_comm_shortPath_process.m     -> main file for output data processing; however, I generally run sectors of this script separatelly on demand; it serves more as a reminder
                                                         for the command lines I use to call the two scripts below (parse and process SpreadingData).
        ./func/parseSpreadingData.m                 -> Parses all the raw output data from all the cases and generates a coherent data structure to be processed
        ./func/processSpreadingData.m               -> Calculates all quantities from the data structure generated above (t-tests included)

    3) Plotting
        ./plot_sim_data -> main plotting function; plots all the figures that are used in the paper


II) Running the simulation

The input parameters of the main simulation scripts described above (Run*.m) are self-evident.
    a) RunNetworkActivationTime(inputFile, useBinaryMatrix, Theta, tTotal, nSeeds, isCompetitive, seedType, seedInd, outDirName, outFileType)
        inputFile        -> .mat file containing "matrix" variable with the FA DTI adjacency of a given subject
        useBinaryMatrix  -> replaces all weights by 1 if true
        Theta            -> value of the threshold of the nodes
        tTotal           -> max amount of time steps to run the dynamics
        nSeeds           -> number of seeds to be used
        isCompetitive    -> if true, each seed generates a competing cascade
        seedType         -> 'random': randomly select seeds
                            'fixed': select fixed seedInd for every trial + a quantity of (nSeeds-numel(seedInd)) seeds in 1:size(A,1);
                            'mostconn': selects the nSeeds most connected nodes as seeds
        seedInd          -> double(1:nSim,1:nSeeds) index of 'fixed' seeds
                            if there is less rows than nSim, the remaining rows are a periodic repetition of seedInd(:,:)
                            if there is less columns than nSeeds, the remaining cols are going to be filled up with random seeds
        outDirName       -> name of the output directory
        outFileType      -> leave it as 'mat'

    b) RunPhaseTransitionSingleMatrix(inputFile, useBinaryMatrix, Theta0, Theta1, nTheta, nSim, tTotal, x0, icType, useWeightSampling, outDirName, outFileType)
        inputFile        -> .mat file containing "matrix" variable with the FA DTI adjacency of a given subject
        useBinaryMatrix  -> replaces all weights by 1 if true
        Theta0           -> initial value of the threshold
        Theta1           -> final value of the threshold
        nTheta           -> amount of theta to simulate between theta0 and theta1
        nSim             -> number of trials for each theta
        tTotal           -> max amount of time steps to run the dynamics
        x0               -> initial state of the nodes of the network
        icType           -> 'fixed': the same node is stimulated for every theta and sim
                            'mostconn': the same node (most connected) is stimulated for every theta and sim
                            'random': randomly choose one seed node
        useWeightSampling -> if true, a new weight sampling is generated for each simulation of each theta (leave it as FALSE to use original subject FA weights)
        outDirName       -> name of the output directory
        outFileType      -> leave it as 'mat'

    c) RunSpreadingTimeSingleMatrix(inputFile, useBinaryMatrix, Theta, nSim, tTotal, outDirName, outFileType)
        inputFile        -> .mat file containing "matrix" variable with the FA DTI adjacency of a given subject
        useBinaryMatrix  -> replaces all weights by 1 if true
        Theta            -> value of the threshold of the nodes
        nSim             -> number of trials for each theta
        tTotal           -> max amount of time steps to run the dynamics
        outDirName       -> name of the output directory
        outFileType      -> leave it as 'mat'

    d) RunTwoSeedsSpreading(inputFile, useBinaryMatrix, Theta, tTotal, isCompetitive, outDirName, outFileType)
        inputFile        -> .mat file containing "matrix" variable with the FA DTI adjacency of a given subject
        useBinaryMatrix  -> replaces all weights by 1 if true
        Theta            -> value of the threshold of the nodes
        tTotal           -> max amount of time steps to run the dynamics
        isCompetitive    -> if true, each seed generates a competing cascade
        outDirName       -> name of the output directory
        outFileType      -> leave it as 'mat'

In order to run in the queue for every case, there are shell scripts in the directory: /data/noel/noel7/mauricio/programs/matlab to do that.
The scripts are:
    a) qsub_netact.sh
    b) qsub_phase_trans_sim.sh
    c) qsub_spreading_time.sh
    d) qsub_twoseeds.sh


PLEASE, DO NOT MODIFY THESE SCRIPTS. If you want to run test simulations, do so by copying these scripts if you need to change any input parameter out of curiosity.


III) Example

    0) Organize the desired input files in an input directory (I will use 4 patients + 4 controls):
    mkdir ~/test_input
    cp /data/noel/noel7/mauricio/data/AALmatrix_N90/flipped_cols/TLE_0361_1_FA_AAL_lowres.mat ~/test_input
    cp /data/noel/noel7/mauricio/data/AALmatrix_N90/flipped_cols/TLE_0362_1_FA_AAL_lowres.mat ~/test_input
    cp /data/noel/noel7/mauricio/data/AALmatrix_N90/flipped_cols/TLE_0369_1_FA_AAL_lowres_flip.mat ~/test_input
    cp /data/noel/noel7/mauricio/data/AALmatrix_N90/flipped_cols/TLE_0380_1_FA_AAL_lowres_flip.mat ~/test_input
    cp /data/noel/noel7/mauricio/data/AALmatrix_N90/flipped_cols/TLE_301_1_FA_AAL_lowres.mat ~/test_input
    cp /data/noel/noel7/mauricio/data/AALmatrix_N90/flipped_cols/TLE_304_1_FA_AAL_lowres.mat ~/test_input
    cp /data/noel/noel7/mauricio/data/AALmatrix_N90/flipped_cols/TLE_305_2_FA_AAL_lowres.mat ~/test_input
    cp /data/noel/noel7/mauricio/data/AALmatrix_N90/flipped_cols/TLE_306_1_FA_AAL_lowres.mat ~/test_input

    1) Create an output directory:
    mkdir ~/test_output

    2) change directory and copy scripts qsub_*.sh
    cd /data/noel/noel7/mauricio/programs/matlab
    cp qsub_netact.sh qsub_netact_BC.sh
    cp qsub_phase_trans_sim.sh  qsub_phase_trans_sim_BC.sh
    cp qsub_spreading_time.sh qsub_spreading_time_BC.sh
    cp qsub_twoseeds.sh qsub_twoseeds_BC.sh

    3) add the following lines somewhere in the middle of the file
        a) qsub_netact_BC.sh (COMMENT OUT ALL THE OTHER LINES that start with ./run_net_act.sh ):
        ./run_net_act.sh ~/test_input 0 0.00 100 1 0 fixed "\"(1:90)'\"" ~/test_output/netact_N90_theta0.00 noel64.q
        ./run_net_act.sh ~/test_input 0 0.01 100 1 0 fixed "\"(1:90)'\"" ~/test_output/netact_N90_theta0.01 noel64.q
        ./run_net_act.sh ~/test_input 0 0.03 100 1 0 fixed "\"(1:90)'\"" ~/test_output/netact_N90_theta0.03 noel64.q
        ./run_net_act.sh ~/test_input 0 0.05 100 1 0 fixed "\"(1:90)'\"" ~/test_output/netact_N90_theta0.05 noel64.q
        
        b) qsub_phase_trans_sim_BC.sh (COMMENT OUT ALL THE OTHER LINES that start with ./run_phase_transition.sh ):
        ./run_phase_transition.sh ~/test_input 0 0 1 200 200 1000 random 0 ~/test_output/pt_nT200_nSim200_r_w noel64.q

        c) qsub_spreading_time_BC.sh (COMMENT OUT ALL THE OTHER LINES that start with ./run_spreading_time.sh ):
        ./run_spreading_time.sh ~/test_input 0 0.00 1 100 ~/test_output/sp_N90_theta0.00_w noel64.q
        ./run_spreading_time.sh ~/test_input 0 0.01 1 100 ~/test_output/sp_N90_theta0.01_w noel64.q
        ./run_spreading_time.sh ~/test_input 0 0.03 1 100 ~/test_output/sp_N90_theta0.03_w noel64.q
        ./run_spreading_time.sh ~/test_input 0 0.05 1 100 ~/test_output/sp_N90_theta0.05_w noel64.q

        d) qsub_twoseeds_BC.sh (COMMENT OUT ALL THE OTHER LINES that start with ./run_twoseeds_spread.sh ):
        ./run_twoseeds_spread.sh ~/test_input 0 0.00 100 0 ~/test_output/2s_N90_theta0.00_w_coop noel64.q
        ./run_twoseeds_spread.sh ~/test_input 0 0.00 100 1 ~/test_output/2s_N90_theta0.00_w_comp noel64.q
        ./run_twoseeds_spread.sh ~/test_input 0 0.01 100 0 ~/test_output/2s_N90_theta0.01_w_coop noel64.q
        ./run_twoseeds_spread.sh ~/test_input 0 0.01 100 1 ~/test_output/2s_N90_theta0.01_w_comp noel64.q

    4) submit the scripts to the queue:
    ./qsub_netact_BC.sh && ./qsub_phase_trans_sim_BC.sh && ./qsub_spreading_time_BC.sh && ./qsub_twoseeds_BC.sh

    5) open MATLAB and go to the /data/noel/noel7/mauricio/programs/matlab/spreading directory
    cd /data/noel/noel7/mauricio/programs/matlab/spreading
    matlab &

Now, inside MATLAB:

    6) add the directory 'func' to path:
    addpath('func')

    6) generate communities and shortest path data:
    load('/data/noel/noel7/mauricio/data/AAL_data/aal_labels_N90.mat') % inputs aal_labels variable
    calcCommunities('~/test_input/TLE*_FA_*.mat', 1:0.5:3.5, 10, aal_labels, true, '~/test_output/comm.mat');
    calcShortestPath('~/test_input/TLE*_FA_*.mat', true, '~/test_output/shortPath.mat');

    7) parse the raw spreading data:
    parseSpreadingData('~/test_output','~/test_output/si_data.mat');

    8) process the si_data.mat structure
    processSpreadingData('~/test_output/si_data.mat', '', '', { {'RTLE','LTLE'} }, { {'HS'} }, 1, 0, [], { 'AlphaTTest',0.01,'AlphaFDR',0.05,'SplitNegPos','on' }, 'sim_data_alpha0.01_spTT.mat', 0, 0);

Now there's a file ~/test_output/sim_data_alpha0.01_spTT.mat that can be used to plot

    9) outside of matlab, create output directory for the figures:
    mkdir ~/test_output/figs

    10) open script plot_sim_data.m
    open plot_sim_data.m

    11) uncomment line 13 and comment out line 14, such that:
    file: plot_sim_data.m
    13 | N = 90; % size of the parcellation of the main results
    14 | % N = 306; % size of the parcellation of the main results

    12) comment out lines 47--49 and add just below these lines the following:
    file: plot_sim_data.m
    47 | % outputDir = fullfile(outputMainDirN,['full_alpha0.01_spTT',odirSuffix,filesep]);
    48 | % si_n90 = load(fullfile(dataDir,'sirs','flipped_cols','N90','sim_data_alpha0.01_spTT.mat'));
    49 | % si_n306 = load(fullfile(dataDir,'sirs','flipped_cols','N306','sim_data_alpha0.01_spTT.mat'));
    50 | outputDir = '~/test_output/figs';
    51 | si_n90 = load('~/test_output/sim_data_alpha0.01_spTT.mat');
    52 | si_n306 = [];

    13) execute plot_sim_data.m
    plot_sim_data