diff --git a/figures/MAM2EBRAINS/M2E_compute_louvain_communities.py b/figures/MAM2EBRAINS/M2E_compute_louvain_communities.py
deleted file mode 100644
index 629e7bd..0000000
--- a/figures/MAM2EBRAINS/M2E_compute_louvain_communities.py
+++ /dev/null
@@ -1,82 +0,0 @@
-# import community
-from community import community_louvain
-import csv
-import json
-import networkx as nx
-import numpy as np
-import os
-
-def compute_communities(M, data_path, label):
-    """
-    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).
-
-    Parameters:
-        - M (MultiAreaModel): The M object containing the area list.
-        - data_path (str): The path to the data directory.
-        - label (str): The label for the data.
-
-    Returns:
-        None
-    """
-    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_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)
diff --git a/figures/MAM2EBRAINS/M2E_visualize_fc.py b/figures/MAM2EBRAINS/M2E_visualize_fc.py
index 089f24f..cb38740 100644
--- a/figures/MAM2EBRAINS/M2E_visualize_fc.py
+++ b/figures/MAM2EBRAINS/M2E_visualize_fc.py
@@ -15,7 +15,6 @@
 sys.path.append('./figures/Schmidt2018')
 
 from M2E_compute_fc import compute_fc
-from M2E_compute_louvain_communities import compute_communities
 
 cmap = pl.cm.coolwarm
 cmap = cmap.from_list('mycmap', [myblue, 'white', myred], N=256)
@@ -94,7 +93,6 @@ def visualize_fc(M, data_path):
     
     # Compute functional connectivity
     compute_fc(M, data_path, label)
-    compute_communities(M, data_path, label)
     
     """
     Figure layout
@@ -161,21 +159,13 @@ def visualize_fc(M, data_path):
         sim_FC[label] = np.load(fn)
 
     label = M.simulation.label
-    fn = os.path.join(data_path,
-                      label,
-                      'Analysis',
-                      'FC_synaptic_input_communities.json')
-    with open(fn, 'r') as f:
-        part_sim = json.load(f)
-    part_sim_list = [part_sim[area] for area in M.area_list]
-    part_sim_index = np.argsort(part_sim_list, kind='mergesort')
-    # Manually position MDP in between the two clusters for visual purposes
-    ind_MDP = M.area_list.index('MDP')
-    ind_MDP_index = np.where(part_sim_index == ind_MDP)[0][0]
-    part_sim_index = np.append(part_sim_index[:ind_MDP_index], part_sim_index[ind_MDP_index+1:])
-    new_ind_MDP_index = np.where(np.array(part_sim_list)[part_sim_index] == 0.)[0][-1]
-    part_sim_index = np.insert(part_sim_index, new_ind_MDP_index+1, ind_MDP)
+    areas_reordered = ['V1', 'V2', 'VP', 'V4t', 'V4', 'VOT', 'MSTd', 'PITv', 
+                   'PITd', 'CITv', 'CITd', 'AITv', 'AITd', 'MDP', 'V3', 'V3A',
+                   'MT', 'PIP', 'PO', 'DP', 'MIP', 'VIP', 'LIP', 'MSTI', 
+                   'FEF', 'TF', 'FST', '7a', 'STPp', 'STPa', '46', 'TH']
 
+    part_sim = {area: M.area_list.index(area) for area in areas_reordered if area in M.area_list}
+    part_sim_index = list(part_sim.values())
 
     """
     Plotting