Merge pull request #61 from didi-hou/functional_connectivity

Functional connectivity added to jupyter-notebook tutorial

figures/MAM2EBRAINS/M2E_compute_fc.py

@@ -0,0 +1,56 @@
+import correlation_toolbox.helper as ch
+import numpy as np
+import os
+import sys
+
+from multiarea_model import MultiAreaModel
+from scipy.spatial.distance import pdist
+from scipy.spatial.distance import squareform
+
+from M2E_compute_synaptic_input import compute_synaptic_input
+
+"""
+Compute the functional connectivity between all areas of a given
+simulation based on their time series of spiking rates or their
+estimated BOLD signal.
+"""
+
+# data_path = sys.argv[1]
+# label = sys.argv[2]
+# method = sys.argv[3]
+
+def compute_fc(M, data_path, label):
+    # compute synaptic input
+    for area in M.area_list:
+        compute_synaptic_input(M, data_path, label, area)
+
+    method = "synaptic_input"
+
+    load_path = os.path.join(data_path,
+                             label,
+                             'Analysis',
+                             method)
+    save_path = os.path.join(data_path,
+                             label,
+                             'Analysis')
+
+    # """
+    # Create MultiAreaModel instance to have access to data structures
+    # """
+    # M = MultiAreaModel({})
+
+    time_series = []
+    for area in M.area_list:
+        fn = os.path.join(load_path,
+                          '{}_{}.npy'.format(method, area))
+        si = np.load(fn)
+        if method == 'bold_signal':  # Cut off the long initial transient of the BOLD signal
+            si = si[5000:]
+        time_series.append(ch.centralize(si, units=True))
+
+    D = pdist(time_series, metric='correlation')
+    correlation_matrix = 1. - squareform(D)
+
+    np.save(os.path.join(save_path,
+                         'functional_connectivity_{}.npy'.format(method)),
+            correlation_matrix)

figures/MAM2EBRAINS/M2E_compute_louvain_communities.py

@@ -0,0 +1,88 @@
+# import community
+from community import community_louvain
+import csv
+import json
+import networkx as nx
+import numpy as np
+import os
+import sys
+
+from multiarea_model.multiarea_model import MultiAreaModel
+
+"""
+Determines communities in the functional connectivity of either the
+experimental fMRI data used in Schmidt et al. 2018 or of a given
+simulation (the functional connectivity being based either on spike
+rates or an estimated BOLD signal).
+"""
+
+# data_path = sys.argv[1]
+# label = sys.argv[2]
+# method = sys.argv[3]
+
+# """
+# Create MultiAreaModel instance to have access to data structures
+# """
+# M = MultiAreaModel({})
+
+def compute_communities(M, data_path, label):
+    method = "synaptic_input"
+
+    if label == 'exp':
+        load_path = ''
+
+        func_conn_data = {}
+
+        with open('./figures/Schmidt2018_dyn/Fig8_exp_func_conn.csv', 'r') as f:
+            myreader = csv.reader(f, delimiter='\t')
+            # Skip first 3 lines
+            next(myreader)
+            next(myreader)
+            next(myreader)
+            areas = next(myreader)
+            for line in myreader:
+                dict_ = {}
+                for i in range(len(line)):
+                    dict_[areas[i]] = float(line[i])
+                func_conn_data[areas[myreader.line_num - 5]] = dict_
+
+        FC = np.zeros((len(M.area_list),
+                       len(M.area_list)))
+        for i, area1 in enumerate(M.area_list):
+            for j, area2 in enumerate(M.area_list):
+                FC[i][j] = func_conn_data[area1][area2]
+
+    else:
+        load_path = os.path.join(data_path,
+                                 label,
+                                 'Analysis',
+                                 'functional_connectivity_{}.npy'.format(method))
+        FC = np.load(load_path)
+
+    # Set diagonal to 0
+    for i in range(FC.shape[0]):
+        FC[i][i] = 0.
+
+    G = nx.Graph()
+    for area in M.area_list:
+        G.add_node(area)
+
+    edges = []
+    for i, area in enumerate(M.area_list):
+        for j, area2 in enumerate(M.area_list):
+            edges.append((area, area2, FC[i][j]))
+    G.add_weighted_edges_from(edges)
+
+    # part = community.best_partition(G)
+    part = community_louvain.best_partition(G)
+
+    if label == 'exp':
+        fn = os.path.join('FC_exp_communities.json')
+    else:
+        fn = os.path.join(data_path,
+                          label,
+                          'Analysis',
+                          'FC_{}_communities.json'.format(method))
+
+    with open(fn, 'w') as f:
+        json.dump(part, f)

figures/MAM2EBRAINS/M2E_compute_synaptic_input.py

@@ -0,0 +1,100 @@
+import correlation_toolbox.helper as ch
+import correlation_toolbox.correlation_analysis as corr
+import json
+import numpy as np
+import os
+import sys
+
+from multiarea_model.multiarea_model import MultiAreaModel
+
+# data_path = sys.argv[1]
+# label = sys.argv[2]
+# area = sys.argv[3]
+
+def compute_synaptic_input(M, data_path, label, area):
+    load_path = os.path.join(data_path,
+                             label,
+                             'Analysis',
+                             'rate_time_series_full')
+    save_path = os.path.join(data_path,
+                             label,
+                             'Analysis',
+                             'synaptic_input')
+
+    with open(os.path.join(data_path, label, 'custom_params_{}'.format(label)), 'r') as f:
+        sim_params = json.load(f)
+    # T = sim_params['T']
+    T = M.simulation.params['t_sim']
+
+
+    # """
+    # Create MultiAreaModel instance to have access to data structures
+    # """
+    # connection_params = {'g': -11.,
+    #                      'cc_weights_factor': sim_params['cc_weights_factor'],
+    #                      'cc_weights_I_factor': sim_params['cc_weights_I_factor'],
+    #                      'K_stable': '../SchueckerSchmidt2017/K_prime_original.npy'}
+    # network_params = {'connection_params': connection_params}
+    # M = MultiAreaModel(network_params)
+
+
+    """
+    Synaptic filtering kernel
+    """
+    t = np.arange(0., 20., 1.)
+    tau_syn = M.params['neuron_params']['single_neuron_dict']['tau_syn_ex']
+    kernel = np.exp(-t / tau_syn)
+
+
+    """
+    Load rate time series
+    """
+    rate_time_series = {}
+    for source_area in M.area_list:
+        rate_time_series[source_area] = {}
+        for source_pop in M.structure[source_area]:
+            fn = os.path.join(load_path,
+                              'rate_time_series_full_{}_{}.npy'.format(source_area, source_pop))
+            dat = np.load(fn)
+            rate_time_series[source_area][source_pop] = dat
+
+
+    synaptic_input_list = []
+    N_list = []
+    for pop in M.structure[area]:
+        time_series = np.zeros(int((T - 500.)))
+        for source_area in M.area_list:
+            for source_pop in M.structure[source_area]:
+                weight = M.W[area][pop][source_area][source_pop]
+                time_series += (rate_time_series[source_area][source_pop] *
+                                abs(weight) *
+                                M.K[area][pop][source_area][source_pop])
+        syn_current = np.convolve(kernel, time_series, mode='same')
+        synaptic_input_list.append(syn_current)
+        N_list.append(M.N[area][pop])
+
+        fp = '_'.join(('synaptic_input',
+                       area,
+                       pop))
+        try:
+            os.mkdir(save_path)
+        except FileExistsError:
+            pass
+        np.save('{}/{}.npy'.format(save_path, fp), syn_current)
+
+    synaptic_input_list = np.array(synaptic_input_list)
+    area_time_series = np.average(synaptic_input_list, axis=0, weights=N_list)
+
+    fp = '_'.join(('synaptic_input',
+                   area))
+    np.save('{}/{}.npy'.format(save_path, fp), area_time_series)
+
+    par = {'areas': M.area_list,
+           'pops': 'complete',
+           'resolution': 1.,
+           't_min': 500.,
+           't_max': T}
+    fp = '_'.join(('synaptic_input',
+                   'Parameters.json'))
+    with open('{}/{}'.format(save_path, fp), 'w') as f:
+        json.dump(par, f)