diff --git a/figures/Schmidt2018_dyn/Fig2_bistability.py b/figures/Schmidt2018_dyn/Fig2_bistability.py
index bd16f49c5bc4206d94596abac957822d2e411bee..d0f1fe3dab3c4f65bdc223335366228e3be658e8 100644
--- a/figures/Schmidt2018_dyn/Fig2_bistability.py
+++ b/figures/Schmidt2018_dyn/Fig2_bistability.py
@@ -31,13 +31,21 @@ axes['C2'] = panel_factory.new_empty_panel(2, 2, r'', label_position=-0.25)
 # Simulation
 if LOAD_ORIGINAL_DATA:
     data = {}
-    data_labels = [('LA', '533d73357fbe99f6178029e6054b571b485f40f6'),
-                   ('HA', '0adda4a542c3d5d43aebf7c30d876b6c5fd1d63e'),
-                   ('LA_post', '33fb5955558ba8bb15a3fdce49dfd914682ef3ea')]
-    for key, label in data_labels:
-        fn = os.path.join(original_data_path, label, 'Analysis/pop_rates.json')
-        with open(fn, 'r') as f:
-            data[key] = json.load(f)
+    data_labels = [('533d73357fbe99f6178029e6054b571b485f40f6'),
+                   ('0adda4a542c3d5d43aebf7c30d876b6c5fd1d63e'),
+                   ('33fb5955558ba8bb15a3fdce49dfd914682ef3ea')]
+    data_path = original_data_path
+else:
+    from network_simulations import init_models
+    from config import data_path
+    models = init_models('Fig2')
+    data_labels = [M.simulation.label for M in models]
+
+keys = ['LA', 'HA', 'LA_post']
+for key, label in zip(keys, data_labels):
+    fn = os.path.join(data_path, label, 'Analysis/pop_rates.json')
+    with open(fn, 'r') as f:
+        data[key] = json.load(f)
 
     """
     Create MultiAreaModel instance to have access to data structures
diff --git a/figures/Schmidt2018_dyn/Fig3_ground_state_chi1.py b/figures/Schmidt2018_dyn/Fig3_ground_state_chi1.py
index 0e0781b002945daf9b297e120fdcfce38b1f1dee..9fc559ac510da94b154b8fea2e53133bcd25a74f 100644
--- a/figures/Schmidt2018_dyn/Fig3_ground_state_chi1.py
+++ b/figures/Schmidt2018_dyn/Fig3_ground_state_chi1.py
@@ -89,53 +89,61 @@ LOAD_ORIGINAL_DATA = True
 
 if LOAD_ORIGINAL_DATA:
     label = '33fb5955558ba8bb15a3fdce49dfd914682ef3ea'
+    data_path = original_data_path
+else:
+    from network_simulations import init_models
+    from config import data_path
+    models = init_models('Fig3')
+    label = models[0].simulation.label
 
-    """
-    Create MultiAreaModel instance to have access to data structures
-    """
-    M = MultiAreaModel({})
-
-    # spike data
-    spike_data = {}
-    for area in areas:
-        spike_data[area] = {}
-        for pop in M.structure[area]:
-            spike_data[area][pop] = np.load(os.path.join(original_data_path,
-                                                         label,
-                                                         'recordings',
-                                                         'spike_data_{}_{}.npy'.format(area, pop)))
-    # stationary firing rates
-    fn = os.path.join(original_data_path, label, 'Analysis', 'pop_rates.json')
-    with open(fn, 'r') as f:
-        pop_rates = json.load(f)
-
-    # time series of firing rates
-    rate_time_series = {}
-    for area in M.area_list:
-        fn = os.path.join(original_data_path, label,
-                          'Analysis',
-                          'rate_time_series',
-                          'rate_time_series_{}.npy'.format(area))
-        rate_time_series[area] = np.load(fn)
-
-    # time series of firing rates convolved with a kernel
-    rate_time_series_auto_kernel = {}
-    for area in M.area_list:
-        fn = os.path.join(original_data_path, label,
-                          'Analysis',
-                          'rate_time_series_auto_kernel',
-                          'rate_time_series_auto_kernel_{}.npy'.format(area))
-        rate_time_series_auto_kernel[area] = np.load(fn)
-
-    # local variance revised (LvR)
-    fn = os.path.join(original_data_path, label, 'Analysis', 'pop_LvR.json')
-    with open(fn, 'r') as f:
-        pop_LvR = json.load(f)
     
-    # correlation coefficients
-    fn = os.path.join(original_data_path, label, 'Analysis', 'corrcoeff.json')
-    with open(fn, 'r') as f:
-        corrcoeff = json.load(f)
+"""
+Create MultiAreaModel instance to have access to data structures
+"""
+M = MultiAreaModel({})
+
+# spike data
+spike_data = {}
+for area in areas:
+    spike_data[area] = {}
+    for pop in M.structure[area]:
+        spike_data[area][pop] = np.load(os.path.join(data_path,
+                                                     label,
+                                                     'recordings',
+                                                     '{}-spikes-{}-{}.npy'.format(label,
+                                                                                  area, pop)))
+# stationary firing rates
+fn = os.path.join(data_path, label, 'Analysis', 'pop_rates.json')
+with open(fn, 'r') as f:
+    pop_rates = json.load(f)
+
+# time series of firing rates
+rate_time_series = {}
+for area in M.area_list:
+    fn = os.path.join(data_path, label,
+                      'Analysis',
+                      'rate_time_series_full',
+                      'rate_time_series_full_{}.npy'.format(area))
+    rate_time_series[area] = np.load(fn)
+
+# time series of firing rates convolved with a kernel
+rate_time_series_auto_kernel = {}
+for area in M.area_list:
+    fn = os.path.join(data_path, label,
+                      'Analysis',
+                      'rate_time_series_auto_kernel',
+                      'rate_time_series_auto_kernel_{}.npy'.format(area))
+    rate_time_series_auto_kernel[area] = np.load(fn)
+
+# local variance revised (LvR)
+fn = os.path.join(data_path, label, 'Analysis', 'pop_LvR.json')
+with open(fn, 'r') as f:
+    pop_LvR = json.load(f)
+
+# correlation coefficients
+fn = os.path.join(data_path, label, 'Analysis', 'corrcoeff.json')
+with open(fn, 'r') as f:
+    corrcoeff = json.load(f)
     
 """
 Plotting
diff --git a/figures/Schmidt2018_dyn/Fig4_metastability.py b/figures/Schmidt2018_dyn/Fig4_metastability.py
index b40c113e9d8a6771c3a88529795e64d7d911fee6..755089af8eb2499f1d6137f2148a058fe711af5f 100644
--- a/figures/Schmidt2018_dyn/Fig4_metastability.py
+++ b/figures/Schmidt2018_dyn/Fig4_metastability.py
@@ -97,14 +97,14 @@ Load data
 """
 chi_list = [1.0, 1.8, 1.9, 2., 2.1, 2.5]
 
+"""
+Create MultiAreaModel instance to have access to data structures
+"""
+M = MultiAreaModel({})
+
 LOAD_ORIGINAL_DATA = True
 
 if LOAD_ORIGINAL_DATA:
-    """
-    Create MultiAreaModel instance to have access to data structures
-    """
-    M = MultiAreaModel({})
-
     labels = ['33fb5955558ba8bb15a3fdce49dfd914682ef3ea',
               '1474e1884422b5b2096d3b7a20fd4bdf388af7e0',
               '99c0024eacc275d13f719afd59357f7d12f02b77',
@@ -112,34 +112,42 @@ if LOAD_ORIGINAL_DATA:
               '08a3a1a88c19193b0af9d9d8f7a52344d1b17498',
               '5bdd72887b191ec22a5abcc04ca4a488ea216e32']
 
-    rate_time_series = {label: {} for label in labels}
-    rate_time_series_pops = {label: {} for label in labels}
-    for label in labels:
-        for area in M.area_list:
-            fn = os.path.join(original_data_path, label,
+    label_stat_rate = '99c0024eacc275d13f719afd59357f7d12f02b77'
+    data_path = original_data_path
+else:
+    from network_simulations import init_models
+    from config import data_path
+    models = init_models('Fig4')
+    labels = [M.simulation.label for M in models]
+    label_stat_rate = labels[2]  # chi=1.9
+
+rate_time_series = {label: {} for label in labels}
+rate_time_series_pops = {label: {} for label in labels}
+for label in labels:
+    for area in M.area_list:
+        fn = os.path.join(data_path, label,
+                          'Analysis',
+                          'rate_time_series_full',
+                          'rate_time_series_full_{}.npy'.format(area))
+        rate_time_series[label][area] = np.load(fn)
+        rate_time_series_pops[label][area] = {}
+        for pop in M.structure[area]:
+            fn = os.path.join(data_path, label,
                               'Analysis',
                               'rate_time_series_full',
-                              'rate_time_series_full_{}.npy'.format(area))
-            rate_time_series[label][area] = np.load(fn)
-            rate_time_series_pops[label][area] = {}
-            for pop in M.structure[area]:
-                fn = os.path.join(original_data_path, label,
-                                  'Analysis',
-                                  'rate_time_series_full',
-                                  'rate_time_series_full_{}_{}.npy'.format(area, pop))
-                rate_time_series_pops[label][area][pop] = np.load(fn)
-        
-        with open(os.path.join(original_data_path, label,
-                               'Analysis',
-                               'rate_time_series_full',
-                               'rate_time_series_full_Parameters.json')) as f:
-            rate_time_series[label]['Parameters'] = json.load(f)
-
-    # stationary firing rates
-    label = '99c0024eacc275d13f719afd59357f7d12f02b77'
-    fn = os.path.join(original_data_path, label, 'Analysis', 'pop_rates.json')
-    with open(fn, 'r') as f:
-        pop_rates = {label: json.load(f)}
+                              'rate_time_series_full_{}_{}.npy'.format(area, pop))
+            rate_time_series_pops[label][area][pop] = np.load(fn)
+
+    with open(os.path.join(data_path, label,
+                           'Analysis',
+                           'rate_time_series_full',
+                           'rate_time_series_full_Parameters.json')) as f:
+        rate_time_series[label]['Parameters'] = json.load(f)
+
+# stationary firing rates
+fn = os.path.join(data_path, label_stat_rate, 'Analysis', 'pop_rates.json')
+with open(fn, 'r') as f:
+    pop_rates = {label: json.load(f)}
 
 
 # Meanfield part: first initialize base class to compute initial rates
diff --git a/figures/Schmidt2018_dyn/Fig5_ground_state.py b/figures/Schmidt2018_dyn/Fig5_ground_state.py
index e5384ebc3c47d9733403ad8f5d04d280050cdfc4..66c0425c2be1a1cdf9a4069c919dbb59d1532aed 100644
--- a/figures/Schmidt2018_dyn/Fig5_ground_state.py
+++ b/figures/Schmidt2018_dyn/Fig5_ground_state.py
@@ -87,58 +87,67 @@ Load data
 """
 LOAD_ORIGINAL_DATA = True
 
+
 if LOAD_ORIGINAL_DATA:
     # use T=10500 simulation for spike raster plots
     label_spikes = '3afaec94d650c637ef8419611c3f80b3cb3ff539'
     # and T=100500 simulation for all other panels
     label = '99c0024eacc275d13f719afd59357f7d12f02b77'
-    
-    """
-    Create MultiAreaModel instance to have access to data structures
-    """
-    M = MultiAreaModel({})
-
-    # spike data
-    spike_data = {}
-    for area in areas:
-        spike_data[area] = {}
-        for pop in M.structure[area]:
-            spike_data[area][pop] = np.load(os.path.join(original_data_path,
-                                                         label_spikes,
-                                                         'recordings',
-                                                         'spike_data_{}_{}.npy'.format(area, pop)))
-    # stationary firing rates
-    fn = os.path.join(original_data_path, label, 'Analysis', 'pop_rates.json')
-    with open(fn, 'r') as f:
-        pop_rates = json.load(f)
-
-    # time series of firing rates
-    rate_time_series = {}
-    for area in M.area_list:
-        fn = os.path.join(original_data_path, label,
-                          'Analysis',
-                          'rate_time_series',
-                          'rate_time_series_{}.npy'.format(area))
-        rate_time_series[area] = np.load(fn)
-
-    # time series of firing rates convolved with a kernel
-    rate_time_series_auto_kernel = {}
-    for area in M.area_list:
-        fn = os.path.join(original_data_path, label,
-                          'Analysis',
-                          'rate_time_series_auto_kernel',
-                          'rate_time_series_auto_kernel_{}.npy'.format(area))
-        rate_time_series_auto_kernel[area] = np.load(fn)
-
-    # local variance revised (LvR)
-    fn = os.path.join(original_data_path, label, 'Analysis', 'pop_LvR.json')
-    with open(fn, 'r') as f:
-        pop_LvR = json.load(f)
-    
-    # correlation coefficients
-    fn = os.path.join(original_data_path, label, 'Analysis', 'corrcoeff.json')
-    with open(fn, 'r') as f:
-        corrcoeff = json.load(f)
+    data_path = original_data_path
+else:
+    from network_simulations import init_models
+    from config import data_path
+    models = init_models('Fig5')
+    label_spikes = models[0].simulation.label
+    label = models[1].simulation.label
+
+"""
+Create MultiAreaModel instance to have access to data structures
+"""
+M = MultiAreaModel({})
+
+# spike data
+spike_data = {}
+for area in areas:
+    spike_data[area] = {}
+    for pop in M.structure[area]:
+        spike_data[area][pop] = np.load(os.path.join(data_path,
+                                                     label_spikes,
+                                                     'recordings',
+                                                     '{}-spikes-{}-{}.npy'.format(label_spikes,
+                                                                                  area, pop)))
+# stationary firing rates
+fn = os.path.join(data_path, label, 'Analysis', 'pop_rates.json')
+with open(fn, 'r') as f:
+    pop_rates = json.load(f)
+
+# time series of firing rates
+rate_time_series = {}
+for area in M.area_list:
+    fn = os.path.join(data_path, label,
+                      'Analysis',
+                      'rate_time_series_full',
+                      'rate_time_series_full_{}.npy'.format(area))
+    rate_time_series[area] = np.load(fn)
+
+# time series of firing rates convolved with a kernel
+rate_time_series_auto_kernel = {}
+for area in M.area_list:
+    fn = os.path.join(data_path, label,
+                      'Analysis',
+                      'rate_time_series_auto_kernel',
+                      'rate_time_series_auto_kernel_{}.npy'.format(area))
+    rate_time_series_auto_kernel[area] = np.load(fn)
+
+# local variance revised (LvR)
+fn = os.path.join(data_path, label, 'Analysis', 'pop_LvR.json')
+with open(fn, 'r') as f:
+    pop_LvR = json.load(f)
+
+# correlation coefficients
+fn = os.path.join(data_path, label, 'Analysis', 'corrcoeff.json')
+with open(fn, 'r') as f:
+    corrcoeff = json.load(f)
     
 """
 Plotting
diff --git a/figures/Schmidt2018_dyn/Fig6_comparison_exp_spiking_data.py b/figures/Schmidt2018_dyn/Fig6_comparison_exp_spiking_data.py
index b4c5fa32d59f822c60b32968ad205867cd33efdb..2a0e77e04f8dcfeff56b36914a3169c0de42071d 100644
--- a/figures/Schmidt2018_dyn/Fig6_comparison_exp_spiking_data.py
+++ b/figures/Schmidt2018_dyn/Fig6_comparison_exp_spiking_data.py
@@ -65,30 +65,36 @@ if LOAD_ORIGINAL_DATA:
               '5bdd72887b191ec22a5abcc04ca4a488ea216e32',
               '3afaec94d650c637ef8419611c3f80b3cb3ff539',
               '99c0024eacc275d13f719afd59357f7d12f02b77']
+    data_path = original_data_path
+else:
+    from network_simulations import init_models
+    from config import data_path
+    models = init_models('Fig6')
+    labels = [M.simulation.label for M in models]
 
-    area = 'V1'
-
-    power_spectra = {chi: {} for chi in chi_list}
-    for chi, label in zip(chi_list, labels):
-        fp = os.path.join(original_data_path,
-                          label,
-                          'Analysis',
-                          'power_spectrum_subsample')
-        power_spectra[chi] = {'f': np.load(os.path.join(fp,
-                                                        'power_spectrum_subsample_freq.npy')),
-                              'power': np.load(os.path.join(fp,
-                                                            'power_spectrum_subsample_V1.npy'))}
-    rate_histograms = {chi: {} for chi in chi_list}
-    for chi, label in zip(chi_list, labels):
-        fp = os.path.join(original_data_path,
-                          label,
-                          'Analysis',
-                          'rate_histogram')
-
-        rate_histograms[chi] = {'bins': np.load(os.path.join(fp,
-                                                             'rate_histogram_bins.npy')),
-                                'vals': np.load(os.path.join(fp,
-                                                             'rate_histogram_V1.npy'))}
+area = 'V1'
+
+power_spectra = {chi: {} for chi in chi_list}
+for chi, label in zip(chi_list, labels):
+    fp = os.path.join(data_path,
+                      label,
+                      'Analysis',
+                      'power_spectrum_subsample')
+    power_spectra[chi] = {'f': np.load(os.path.join(fp,
+                                                    'power_spectrum_subsample_freq.npy')),
+                          'power': np.load(os.path.join(fp,
+                                                        'power_spectrum_subsample_V1.npy'))}
+rate_histograms = {chi: {} for chi in chi_list}
+for chi, label in zip(chi_list, labels):
+    fp = os.path.join(data_path,
+                      label,
+                      'Analysis',
+                      'rate_histogram')
+
+    rate_histograms[chi] = {'bins': np.load(os.path.join(fp,
+                                                         'rate_histogram_bins.npy')),
+                            'vals': np.load(os.path.join(fp,
+                                                         'rate_histogram_V1.npy'))}
 
 """
 Load experimental data
diff --git a/figures/Schmidt2018_dyn/Fig7_temporal_hierarchy.py b/figures/Schmidt2018_dyn/Fig7_temporal_hierarchy.py
index 980a66d6bf142ebce437241433c2eede7e39cc51..2b9bac7efbb1b9635b4ff2f3c5cafd8cb295fa9f 100644
--- a/figures/Schmidt2018_dyn/Fig7_temporal_hierarchy.py
+++ b/figures/Schmidt2018_dyn/Fig7_temporal_hierarchy.py
@@ -72,44 +72,53 @@ pl.text(-0.05, 1.0, r'\bfseries{}' + 'G',
 """
 Load data
 """
-LOAD_ORIGINAL_DATA = True
 
+"""
+Create MultiAreaModel instance to have access to data structures
+"""
+M = MultiAreaModel({})
+
+LOAD_ORIGINAL_DATA = True
 if LOAD_ORIGINAL_DATA:
-    """
-    Create MultiAreaModel instance to have access to data structures
-    """
-    M = MultiAreaModel({})
     label = '99c0024eacc275d13f719afd59357f7d12f02b77'
-    rate_time_series = {}
-    for area in M.area_list:
-        fn = os.path.join(original_data_path, label,
-                          'Analysis',
-                          'rate_time_series_full',
-                          'rate_time_series_full_{}.npy'.format(area))
-        rate_time_series[area] = np.load(fn)
+    data_path = original_data_path
+else:
+    from network_simulations import init_models
+    from config import data_path
+    models = init_models('Fig7')
+    label = models[0].simulation.label
 
-    fn = os.path.join(original_data_path, label,
+
+rate_time_series = {}
+for area in M.area_list:
+    fn = os.path.join(data_path, label,
                       'Analysis',
                       'rate_time_series_full',
-                      'rate_time_series_full_Parameters.json')
-    with open(fn, 'r') as f:
-        rate_time_series['Parameters'] = json.load(f)
-
-    cross_correlation = {}
-    for area in M.area_list:
-        cross_correlation[area] = {}
-        for area2 in M.area_list:
-            fn = os.path.join(original_data_path, label,
-                              'Analysis',
-                              'cross_correlation',
-                              'cross_correlation_{}_{}.npy'.format(area, area2))
-            cross_correlation[area][area2] = np.load(fn)
-
-    fn = os.path.join(original_data_path, label,
-                      'Analysis',
-                      'cross_correlation',
-                      'cross_correlation_time.npy')
-    cross_correlation['time'] = np.load(fn)
+                      'rate_time_series_full_{}.npy'.format(area))
+    rate_time_series[area] = np.load(fn)
+
+fn = os.path.join(data_path, label,
+                  'Analysis',
+                  'rate_time_series_full',
+                  'rate_time_series_full_Parameters.json')
+with open(fn, 'r') as f:
+    rate_time_series['Parameters'] = json.load(f)
+
+cross_correlation = {}
+for area in M.area_list:
+    cross_correlation[area] = {}
+    for area2 in M.area_list:
+        fn = os.path.join(data_path, label,
+                          'Analysis',
+                          'cross_correlation',
+                          'cross_correlation_{}_{}.npy'.format(area, area2))
+        cross_correlation[area][area2] = np.load(fn)
+
+fn = os.path.join(data_path, label,
+                  'Analysis',
+                  'cross_correlation',
+                  'cross_correlation_time.npy')
+cross_correlation['time'] = np.load(fn)
 
 
 # Correlation peak
diff --git a/figures/Schmidt2018_dyn/Fig8_interactions.py b/figures/Schmidt2018_dyn/Fig8_interactions.py
index 528dadf6291826326356a0f10d7e587d1b044485..841af43f6b773f1ac7e082da1587e68dfc728563 100644
--- a/figures/Schmidt2018_dyn/Fig8_interactions.py
+++ b/figures/Schmidt2018_dyn/Fig8_interactions.py
@@ -85,7 +85,7 @@ Load data
 """
 Create MultiAreaModel instance to have access to data structures
 """
-conn_params = {'g': -16.,
+conn_params = {'g': -11.,
                'fac_nu_ext_TH': 1.2,
                'fac_nu_ext_5E': 1.125,
                'fac_nu_ext_6E': 1.41666667,
@@ -128,11 +128,9 @@ Simulation data
 """
 LOAD_ORIGINAL_DATA = True
 
-if LOAD_ORIGINAL_DATA:
-    tmin = 500.
-    tmax = 10000.
+cc_weights_factor = [1.0, 1.4, 1.5, 1.6, 1.7, 1.75, 1.8, 2., 2.1, 2.5, 1.9]
 
-    cc_weights_factor = [1.0, 1.4, 1.5, 1.6, 1.7, 1.75, 1.8, 1.9, 2., 2.1, 2.5]
+if LOAD_ORIGINAL_DATA:
     labels = ['33fb5955558ba8bb15a3fdce49dfd914682ef3ea',
               '783cedb0ff27240133e3daa63f5d0b8d3c2e6b79',
               '380856f3b32f49c124345c08f5991090860bf9a3',
@@ -140,42 +138,50 @@ if LOAD_ORIGINAL_DATA:
               'c1876856b1b2cf1346430cf14e8d6b0509914ca1',
               'a30f6fba65bad6d9062e8cc51f5483baf84a46b7',
               '1474e1884422b5b2096d3b7a20fd4bdf388af7e0',
-              '99c0024eacc275d13f719afd59357f7d12f02b77',
               'f18158895a5d682db5002489d12d27d7a974146f',
               '08a3a1a88c19193b0af9d9d8f7a52344d1b17498',
-              '5bdd72887b191ec22a5abcc04ca4a488ea216e32']
-
-    sim_FC = {}
-    for label in labels:
-        fn = os.path.join(original_data_path,
-                          label,
-                          'Analysis',
-                          'functional_connectivity_synaptic_input.npy')
-        sim_FC[label] = np.load(fn)
-
-    sim_FC_bold = {}
-    for label in ['99c0024eacc275d13f719afd59357f7d12f02b77']:
-        fn = os.path.join(original_data_path,
-                          label,
-                          'Analysis',
-                          'functional_connectivity_bold_signal.npy')
-        sim_FC_bold[label] = np.load(fn)
-
-    label = '99c0024eacc275d13f719afd59357f7d12f02b77'
-    fn = os.path.join(original_data_path,
+              '5bdd72887b191ec22a5abcc04ca4a488ea216e32',
+              '99c0024eacc275d13f719afd59357f7d12f02b77']
+    data_path = original_data_path
+    label_plot = labels[-1]  # chi=1.9
+else:
+    from network_simulations import init_models
+    from config import data_path
+    models = init_models('Fig8')
+    labels = [M.simulation.label for M in models]
+
+
+sim_FC = {}
+for label in labels:
+    fn = os.path.join(data_path,
+                      label,
+                      'Analysis',
+                      'functional_connectivity_synaptic_input.npy')
+    sim_FC[label] = np.load(fn)
+
+sim_FC_bold = {}
+for label in [label_plot]:
+    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, ind_MDP)
+                      'functional_connectivity_bold_signal.npy')
+    sim_FC_bold[label] = np.load(fn)
+
+label = label_plot
+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, ind_MDP)
 
     
 def zero_diagonal(matrix):
@@ -231,14 +237,14 @@ def matrix_plot(ax, matrix, index, vlim, pos=None):
 Plotting
 """
 ax = axes['A']
-label = '99c0024eacc275d13f719afd59357f7d12f02b77'
+label = label_plot
 
 
 matrix_plot(ax, zero_diagonal(sim_FC[label]),
             part_sim_index, 1., pos=(0, 0))
 
 ax = axes['B']
-label = '99c0024eacc275d13f719afd59357f7d12f02b77'
+label = label_plot
 
 matrix_plot(ax, zero_diagonal(sim_FC_bold[label]),
             part_sim_index, 1., pos=(0, 0))
@@ -266,7 +272,7 @@ ax.plot(cc_weights_factor[1:], cc_list[1:], '.', ms=10,
 ax.plot(cc_weights_factor[0], cc_list[0], '^', ms=5,
         markeredgecolor='none', label='Sim. vs. Exp.', color='k')
 
-label = '99c0024eacc275d13f719afd59357f7d12f02b77'
+label = label_plot
 cc_bold = np.corrcoef(zero_diagonal(sim_FC_bold[label]).flatten(),
                       zero_diagonal(exp_FC).flatten())[0][1]
 ax.plot([1.9], cc_bold, '.',
diff --git a/figures/Schmidt2018_dyn/Fig9_laminar_interactions.py b/figures/Schmidt2018_dyn/Fig9_laminar_interactions.py
index 0bab2ddd11b38e78c1b9caf81ed2563fbfbc734e..d5477caebca10c7bbf7788ad78c3c237df8c7f89 100644
--- a/figures/Schmidt2018_dyn/Fig9_laminar_interactions.py
+++ b/figures/Schmidt2018_dyn/Fig9_laminar_interactions.py
@@ -39,39 +39,45 @@ with open(os.path.join(datapath, 'viscortex_processed_data.json'), 'r') as f:
     proc = json.load(f)
 arch_types = proc['architecture_completed']
 
+conn_params = {'g': -16.,
+               'fac_nu_ext_TH': 1.2,
+               'fac_nu_ext_5E': 1.125,
+               'fac_nu_ext_6E': 1.41666667,
+               'av_indegree_V1': 3950.,
+               'K_stable': '../SchueckerSchmidt2017/K_prime_original.npy'}
+network_params = {'N_scaling': 1.,
+                  'K_scaling': 1.,
+                  'connection_params': conn_params}
+M = MultiAreaModel(network_params)
+
 
 LOAD_ORIGINAL_DATA = True
 
 if LOAD_ORIGINAL_DATA:
-    conn_params = {'g': -16.,
-                   'fac_nu_ext_TH': 1.2,
-                   'fac_nu_ext_5E': 1.125,
-                   'fac_nu_ext_6E': 1.41666667,
-                   'av_indegree_V1': 3950.,
-                   'K_stable': '../SchueckerSchmidt2017/K_prime_original.npy'}
-    network_params = {'N_scaling': 1.,
-                      'K_scaling': 1.,
-                      'connection_params': conn_params}
-    M = MultiAreaModel(network_params)
-
     label = '99c0024eacc275d13f719afd59357f7d12f02b77'
-
-    gc = {}
-    for area in M.area_list:
-        gc[area] = {}
-        for pop in M.structure[area]:
-            fn = os.path.join(original_data_path,
-                              label,
-                              'Analysis',
-                              'granger_causality',
-                              'granger_causality_{}_{}.json'.format(area, pop))
-            with open(fn, 'r') as f:
-                gc[area][pop] = json.load(f)
-
-    with open('Fig9_{}_significant_channels.json'.format(label), 'r') as f:
-        significant_channels = json.load(f)
-    for typ in significant_channels:
-        significant_channels[typ] = [tuple(pair) for pair in significant_channels[typ]]
+    data_path = original_data_path
+else:
+    from network_simulations import init_models
+    from config import data_path
+    models = init_models('Fig7')
+    label = models[0].simulation.label
+
+gc = {}
+for area in M.area_list:
+    gc[area] = {}
+    for pop in M.structure[area]:
+        fn = os.path.join(data_path,
+                          label,
+                          'Analysis',
+                          'granger_causality',
+                          'granger_causality_{}_{}.json'.format(area, pop))
+        with open(fn, 'r') as f:
+            gc[area][pop] = json.load(f)
+
+with open('Fig9_{}_significant_channels.json'.format(label), 'r') as f:
+    significant_channels = json.load(f)
+for typ in significant_channels:
+    significant_channels[typ] = [tuple(pair) for pair in significant_channels[typ]]
 
 """
 Bottom row
diff --git a/figures/Schmidt2018_dyn/Snakefile b/figures/Schmidt2018_dyn/Snakefile
index b79a3fe4375571159c3807c1290d2fcc82d3f506..531cc6ec8a4c2158082afe62174f39a002073fd6 100644
--- a/figures/Schmidt2018_dyn/Snakefile
+++ b/figures/Schmidt2018_dyn/Snakefile
@@ -6,12 +6,40 @@ area_list = ['V1', 'V2', 'VP', 'V3', 'V3A', 'MT', 'V4t', 'V4', 'VOT', 'MSTd',
              'STPa', '46', 'AITd', 'TH']
 population_list = ['23E', '23I', '4E',  '4I', '5E', '5I', '6E', '6I']
 
-LOAD_ORIGINAL_DATA = True
-
-ORIGINAL_SIMULATIONS = {'all': ['533d73357fbe99f6178029e6054b571b485f40f6',
+LOAD_ORIGINAL_DATA = False
+
+ORIGINAL_SIM_LABELS = {'all': ['533d73357fbe99f6178029e6054b571b485f40f6',
+                               '0adda4a542c3d5d43aebf7c30d876b6c5fd1d63e',
+                               '33fb5955558ba8bb15a3fdce49dfd914682ef3ea',
+                               '3afaec94d650c637ef8419611c3f80b3cb3ff539',
+                               '783cedb0ff27240133e3daa63f5d0b8d3c2e6b79',
+                               '380856f3b32f49c124345c08f5991090860bf9a3',
+                               '5a7c6c2d6d48a8b687b8c6853fb4d98048681045',
+                               'c1876856b1b2cf1346430cf14e8d6b0509914ca1',
+                               'a30f6fba65bad6d9062e8cc51f5483baf84a46b7',
+                               '1474e1884422b5b2096d3b7a20fd4bdf388af7e0',
+                               'f18158895a5d682db5002489d12d27d7a974146f',
+                               '08a3a1a88c19193b0af9d9d8f7a52344d1b17498',
+                               '5bdd72887b191ec22a5abcc04ca4a488ea216e32'],
+                       'Fig1': None,
+                       'Fig2': ['533d73357fbe99f6178029e6054b571b485f40f6',
                                 '0adda4a542c3d5d43aebf7c30d876b6c5fd1d63e',
-                                '33fb5955558ba8bb15a3fdce49dfd914682ef3ea',
-                                '3afaec94d650c637ef8419611c3f80b3cb3ff539',
+                                '33fb5955558ba8bb15a3fdce49dfd914682ef3ea'],
+                       'Fig3': '33fb5955558ba8bb15a3fdce49dfd914682ef3ea',
+                       'Fig4': ['33fb5955558ba8bb15a3fdce49dfd914682ef3ea',
+                                '1474e1884422b5b2096d3b7a20fd4bdf388af7e0',
+                                '99c0024eacc275d13f719afd59357f7d12f02b77',
+                                'f18158895a5d682db5002489d12d27d7a974146f',
+                                '08a3a1a88c19193b0af9d9d8f7a52344d1b17498',
+                                '5bdd72887b191ec22a5abcc04ca4a488ea216e32'],
+                       'Fig5': ['3afaec94d650c637ef8419611c3f80b3cb3ff539',
+                                '99c0024eacc275d13f719afd59357f7d12f02b77'],
+                       'Fig6': ['33fb5955558ba8bb15a3fdce49dfd914682ef3ea',
+                                '5bdd72887b191ec22a5abcc04ca4a488ea216e32',
+                                '99c0024eacc275d13f719afd59357f7d12f02b77',
+                                '3afaec94d650c637ef8419611c3f80b3cb3ff539'],
+                       'Fig7': ['99c0024eacc275d13f719afd59357f7d12f02b77'],
+                       'Fig8': ['33fb5955558ba8bb15a3fdce49dfd914682ef3ea',
                                 '783cedb0ff27240133e3daa63f5d0b8d3c2e6b79',
                                 '380856f3b32f49c124345c08f5991090860bf9a3',
                                 '5a7c6c2d6d48a8b687b8c6853fb4d98048681045',
@@ -20,43 +48,19 @@ ORIGINAL_SIMULATIONS = {'all': ['533d73357fbe99f6178029e6054b571b485f40f6',
                                 '1474e1884422b5b2096d3b7a20fd4bdf388af7e0',
                                 'f18158895a5d682db5002489d12d27d7a974146f',
                                 '08a3a1a88c19193b0af9d9d8f7a52344d1b17498',
-                                '5bdd72887b191ec22a5abcc04ca4a488ea216e32'],
-                        'Fig1': None,
-                        'Fig2': ['533d73357fbe99f6178029e6054b571b485f40f6',
-                                 '0adda4a542c3d5d43aebf7c30d876b6c5fd1d63e',
-                                 '33fb5955558ba8bb15a3fdce49dfd914682ef3ea'],
-                        'Fig3': '33fb5955558ba8bb15a3fdce49dfd914682ef3ea',
-                        'Fig4': ['33fb5955558ba8bb15a3fdce49dfd914682ef3ea',
-                                 '1474e1884422b5b2096d3b7a20fd4bdf388af7e0',
-                                 '99c0024eacc275d13f719afd59357f7d12f02b77',
-                                 'f18158895a5d682db5002489d12d27d7a974146f',
-                                 '08a3a1a88c19193b0af9d9d8f7a52344d1b17498',
-                                 '5bdd72887b191ec22a5abcc04ca4a488ea216e32'],
-                        'Fig5': ['3afaec94d650c637ef8419611c3f80b3cb3ff539',
-                                 '99c0024eacc275d13f719afd59357f7d12f02b77'],
-                        'Fig6': ['33fb5955558ba8bb15a3fdce49dfd914682ef3ea',
-                                 '5bdd72887b191ec22a5abcc04ca4a488ea216e32',
-                                 '99c0024eacc275d13f719afd59357f7d12f02b77',
-                                 '3afaec94d650c637ef8419611c3f80b3cb3ff539'],
-                        'Fig7': ['99c0024eacc275d13f719afd59357f7d12f02b77'],
-                        'Fig8': ['33fb5955558ba8bb15a3fdce49dfd914682ef3ea',
-                                 '783cedb0ff27240133e3daa63f5d0b8d3c2e6b79',
-                                 '380856f3b32f49c124345c08f5991090860bf9a3',
-                                 '5a7c6c2d6d48a8b687b8c6853fb4d98048681045',
-                                 'c1876856b1b2cf1346430cf14e8d6b0509914ca1',
-                                 'a30f6fba65bad6d9062e8cc51f5483baf84a46b7',
-                                 '1474e1884422b5b2096d3b7a20fd4bdf388af7e0',
-                                 'f18158895a5d682db5002489d12d27d7a974146f',
-                                 '08a3a1a88c19193b0af9d9d8f7a52344d1b17498',
-                                 '5bdd72887b191ec22a5abcc04ca4a488ea216e32',
-                                 '99c0024eacc275d13f719afd59357f7d12f02b77'],
-                        'Fig8': '99c0024eacc275d13f719afd59357f7d12f02b77',
-                        'Fig9': '99c0024eacc275d13f719afd59357f7d12f02b77'}
+                                '5bdd72887b191ec22a5abcc04ca4a488ea216e32',
+                                '99c0024eacc275d13f719afd59357f7d12f02b77'],
+                       'Fig9': ['99c0024eacc275d13f719afd59357f7d12f02b77']}
 
 
 if LOAD_ORIGINAL_DATA:
     DATA_DIR = original_data_path
-    SIMULATIONS = ORIGINAL_SIMULATIONS
+    SIM_LABELS = ORIGINAL_SIM_LABELS_LABELS
+else:
+    from network_simulations import create_label_dict
+    from config import data_path
+    DATA_DIR = data_path
+    SIM_LABELS = create_label_dict()
 
 rule all:
     input:
@@ -74,7 +78,7 @@ include: './Snakefile_preprocessing'
 
 rule Fig2_bistability:
     input:
-        expand(os.path.join(DATA_DIR, '{simulation}', 'Analysis/pop_rates.json'), simulation=SIMULATIONS['Fig2'])
+        expand(os.path.join(DATA_DIR, '{simulation}', 'Analysis/pop_rates.json'), simulation=SIM_LABELS['Fig2'])
     output:
         'Fig2_bistability.eps'
     shell:
@@ -82,14 +86,14 @@ rule Fig2_bistability:
     
 rule Fig3_ground_state_chi1:
     input:
-        os.path.join(DATA_DIR, SIMULATIONS['Fig3'], 'Analysis/pop_rates.json'),
-        os.path.join(DATA_DIR, SIMULATIONS['Fig3'], 'Analysis/pop_LvR.json'),
-        os.path.join(DATA_DIR, SIMULATIONS['Fig3'], 'Analysis/corrcoeff.json'),
-        expand(os.path.join(DATA_DIR, SIMULATIONS['Fig3'], 'recordings', '-'.join((SIMULATIONS['Fig3'], 'spikes-{area}-{pop}.npy'))),
+        os.path.join(DATA_DIR, SIM_LABELS['Fig3'], 'Analysis/pop_rates.json'),
+        os.path.join(DATA_DIR, SIM_LABELS['Fig3'], 'Analysis/pop_LvR.json'),
+        os.path.join(DATA_DIR, SIM_LABELS['Fig3'], 'Analysis/corrcoeff.json'),
+        expand(os.path.join(DATA_DIR, SIM_LABELS['Fig3'], 'recordings', '-'.join((SIM_LABELS['Fig3'], 'spikes-{area}-{pop}.npy'))),
                area=['V1', 'V2', 'FEF'], pop=population_list),
-        expand(os.path.join(DATA_DIR, SIMULATIONS['Fig3'], 'Analysis', 'rate_time_series_full', 'rate_time_series_full_{area}.npy'),
+        expand(os.path.join(DATA_DIR, SIM_LABELS['Fig3'], 'Analysis', 'rate_time_series_full', 'rate_time_series_full_{area}.npy'),
                area=area_list),
-        expand(os.path.join(DATA_DIR, SIMULATIONS['Fig3'], 'Analysis', 'rate_time_series_auto_kernel', 'rate_time_series_auto_kernel_{area}.npy'), area=area_list),
+        expand(os.path.join(DATA_DIR, SIM_LABELS['Fig3'], 'Analysis', 'rate_time_series_auto_kernel', 'rate_time_series_auto_kernel_{area}.npy'), area=area_list),
     output:
         'Fig3_ground_state_chi1.eps'
     shell:
@@ -104,7 +108,7 @@ rule Fig4_theory_calculations:
 rule Fig4_metastability:
     input:
         expand(os.path.join(DATA_DIR, '{simulation}', 'Analysis', 'rate_time_series_full', 'rate_time_series_full_V1.npy'),
-               simulation=SIMULATIONS['Fig4']),
+               simulation=SIM_LABELS['Fig4']),
         expand('Fig4_theory_data/results_{cc_weights_factor}.npy',
                cc_weights_factor=[1.0, 1.8, 1.9, 2., 2.1, 2.5])
     output:
@@ -114,14 +118,14 @@ rule Fig4_metastability:
     
 rule Fig5_ground_state:
     input:
-        os.path.join(DATA_DIR, SIMULATIONS['Fig5'][1], 'Analysis', 'pop_rates.json'),
-        os.path.join(DATA_DIR, SIMULATIONS['Fig5'][1], 'Analysis', 'pop_LvR.json'),
-        os.path.join(DATA_DIR, SIMULATIONS['Fig5'][1], 'Analysis', 'corrcoeff.json'),
-        expand(os.path.join(DATA_DIR, SIMULATIONS['Fig5'][0], 'recordings', '-'.join((SIMULATIONS['Fig5'][0], 'spikes-{area}-{pop}.npy'))),
+        os.path.join(DATA_DIR, SIM_LABELS['Fig5'][1], 'Analysis', 'pop_rates.json'),
+        os.path.join(DATA_DIR, SIM_LABELS['Fig5'][1], 'Analysis', 'pop_LvR.json'),
+        os.path.join(DATA_DIR, SIM_LABELS['Fig5'][1], 'Analysis', 'corrcoeff.json'),
+        expand(os.path.join(DATA_DIR, SIM_LABELS['Fig5'][0], 'recordings', '-'.join((SIM_LABELS['Fig5'][0], 'spikes-{area}-{pop}.npy'))),
                area=['V1', 'V2', 'FEF'], pop=population_list),
-        expand(os.path.join(DATA_DIR, SIMULATIONS['Fig5'][1], 'Analysis', 'rate_time_series_full', 'rate_time_series_full_{area}.npy'),
+        expand(os.path.join(DATA_DIR, SIM_LABELS['Fig5'][1], 'Analysis', 'rate_time_series_full', 'rate_time_series_full_{area}.npy'),
                area=area_list),
-        expand(os.path.join(DATA_DIR, SIMULATIONS['Fig5'][1], 'Analysis', 'rate_time_series_auto_kernel', 'rate_time_series_auto_kernel_{area}.npy'), area=area_list),
+        expand(os.path.join(DATA_DIR, SIM_LABELS['Fig5'][1], 'Analysis', 'rate_time_series_auto_kernel', 'rate_time_series_auto_kernel_{area}.npy'), area=area_list),
     output:
         'Fig5_ground_state.eps'
     shell:
@@ -130,11 +134,11 @@ rule Fig5_ground_state:
 rule Fig6_comparison_exp_spiking_data:
     input:
         expand(os.path.join(DATA_DIR, '{simulation}', 'Analysis', 'power_spectrum_subsample', 'power_spectrum_subsample_V1.npy'),
-               simulation=SIMULATIONS['Fig6']),
+               simulation=SIM_LABELS['Fig6']),
         expand(os.path.join(DATA_DIR, '{simulation}', 'Analysis', 'rate_histogram', 'rate_histogram_V1.npy'),
-               simulation=SIMULATIONS['Fig6']),
+               simulation=SIM_LABELS['Fig6']),
         expand(os.path.join(DATA_DIR,'{simulation}', 'Analysis', 'rate_histogram', 'rate_histogram_V1.npy'),
-               simulation=SIMULATIONS['Fig6'][:2] + SIMULATIONS['Fig6'][-1:]),
+               simulation=SIM_LABELS['Fig6'][:2] + SIM_LABELS['Fig6'][-1:]),
         os.path.join(chu2014_path, 'Analysis', 'spike_data.npy'),
         os.path.join(chu2014_path, 'Analysis', 'neuron_depths.npy'),
         os.path.join(chu2014_path, 'Analysis', 'power_spectrum_freq.npy'),
@@ -158,9 +162,9 @@ rule Fig6_comparison_exp_spiking_data:
 rule Fig7_temporal_hierarchy:
     input:
         expand(os.path.join(DATA_DIR, '{simulation}', 'Analysis', 'rate_time_series_full', 'rate_time_series_full_{area}.npy'),
-               simulation=SIMULATIONS['Fig7'], area=area_list),
+               simulation=SIM_LABELS['Fig7'], area=area_list),
         expand(os.path.join(DATA_DIR, '{simulation}', 'Analysis', 'cross_correlation', 'cross_correlation_{area1}_{area2}.npy'),
-               simulation=SIMULATIONS['Fig7'], area1=area_list, area2=area_list),
+               simulation=SIM_LABELS['Fig7'], area1=area_list, area2=area_list),
     output:
         'Fig7_temporal_hierarchy.eps'
     shell:
@@ -169,13 +173,13 @@ rule Fig7_temporal_hierarchy:
 rule Fig8_interactions:
     input:
         expand(os.path.join(DATA_DIR, '{simulation}', 'Analysis', 'functional_connectivity_synaptic_input.npy'),
-               simulation=SIMULATIONS['Fig8']),
+               simulation=SIM_LABELS['Fig8']),
         expand(os.path.join(DATA_DIR, '{simulation}', 'Analysis', 'FC_synaptic_input_communities.json'),
-               simulation=SIMULATIONS['Fig8'][-1]),
+               simulation=SIM_LABELS['Fig8'][-1]),
         'Fig8_exp_func_conn.csv',
         'FC_exp_communities.json',
         expand(os.path.join(DATA_DIR, '{simulation}', 'Analysis', 'functional_connectivity_bold_signal.npy'),
-               simulation=SIMULATIONS['Fig8'][-1])
+               simulation=SIM_LABELS['Fig8'][-1])
     output:
         'Fig8_interactions.eps'
     shell:
@@ -184,14 +188,14 @@ rule Fig8_interactions:
 rule Fig9_laminar_interactions:
     input:
         expand(os.path.join(DATA_DIR, '{simulation}', 'Analysis', 'granger_causality', 'granger_causality_{area}_{pop}.json'),
-               simulation=SIMULATIONS['Fig9'], area=area_list, pop=population_list),
-        'Fig9_{}_significant_channels.json'.format(SIMULATIONS['Fig9']),
-        'Fig9_{}_HL_interactions.eps'.format(SIMULATIONS['Fig9']),
-        'Fig9_{}_HZ_interactions.eps'.format(SIMULATIONS['Fig9']),
-        'Fig9_{}_LH_interactions.eps'.format(SIMULATIONS['Fig9']),
-        'Fig9_{}_HL_paths.eps'.format(SIMULATIONS['Fig9']),
-        'Fig9_{}_HZ_paths.eps'.format(SIMULATIONS['Fig9']),
-        'Fig9_{}_LH_paths.eps'.format(SIMULATIONS['Fig9']),
+               simulation=SIM_LABELS['Fig9'], area=area_list, pop=population_list),
+        'Fig9_{}_significant_channels.json'.format(SIM_LABELS['Fig9'][0]),
+        'Fig9_{}_HL_interactions.eps'.format(SIM_LABELS['Fig9'][0]),
+        'Fig9_{}_HZ_interactions.eps'.format(SIM_LABELS['Fig9'][0]),
+        'Fig9_{}_LH_interactions.eps'.format(SIM_LABELS['Fig9'][0]),
+        'Fig9_{}_HL_paths.eps'.format(SIM_LABELS['Fig9'][0]),
+        'Fig9_{}_HZ_paths.eps'.format(SIM_LABELS['Fig9'][0]),
+        'Fig9_{}_LH_paths.eps'.format(SIM_LABELS['Fig9'][0]),
     output:
         'Fig9_laminar_interactions.eps'
     shell:
diff --git a/figures/Schmidt2018_dyn/Snakefile_preprocessing b/figures/Schmidt2018_dyn/Snakefile_preprocessing
index ceca19d08e75c2eb23af9265c562e788a503bf3e..ff3b593706866c82a758e52f4daa558d4ade8280 100644
--- a/figures/Schmidt2018_dyn/Snakefile_preprocessing
+++ b/figures/Schmidt2018_dyn/Snakefile_preprocessing
@@ -85,7 +85,7 @@ rule process_chu2014_data:
         os.path.join(chu2014_path, 'Analysis', 'rate_histogram_high.npy'),
         os.path.join(chu2014_path, 'Analysis', 'rate_histogram_full.npy')
     shell:
-        'python3 process_chu2014_data.py {} {}'.format(DATA_DIR, SIMULATIONS['Fig6'][-1])
+        'python3 process_chu2014_data.py {} {}'.format(DATA_DIR, SIM_LABELS['Fig6'][-1])
 
 rule cross_correlation:
     input:
diff --git a/figures/Schmidt2018_dyn/network_simulations.py b/figures/Schmidt2018_dyn/network_simulations.py
new file mode 100644
index 0000000000000000000000000000000000000000..3ab10ce1ba1e058f02704ec427e855f7ed4d2101
--- /dev/null
+++ b/figures/Schmidt2018_dyn/network_simulations.py
@@ -0,0 +1,216 @@
+import copy
+
+from multiarea_model import MultiAreaModel
+from start_jobs import start_job
+from config import submit_cmd, jobscript_template
+
+"""
+This script provides the code to execute all simulations presented in
+Schmidt M, Bakker R, Shen K, Bezgin B, Diesmann M & van Albada SJ
+(2018) A multi-scale layer-resolved spiking network model of
+resting-state dynamics in macaque cortex.
+
+Needs to be simulated with sufficient
+resources, for instance on a compute cluster.
+"""
+
+# Common parameter settings
+input_params = {'rate_ext': 10.}
+neuron_params = {'V0_mean': -150.,
+                 'V0_sd': 50.}
+sim_params = {'num_processes': 20,  # Needs to be adapted to the HPC system used
+              'local_num_threads': 24,  # Needs to be adapted to the HPC system used
+              'recording_dict': {'record_vm': False}}
+
+"""
+Simulation with  kappa = 1. leading to the low-activity fixed point
+shown in Fig. 2.
+"""
+d = {}
+
+sim_params.update({'t_sim': 10500.})
+conn_params = {'g': -16.,
+               'fac_nu_ext_TH': 1.2,
+               'fac_nu_ext_5E': 1.,
+               'fac_nu_ext_6E': 1.,
+               'av_indegree_V1': 3950.}
+network_params = {'N_scaling': 1.,
+                  'K_scaling': 1.,
+                  'connection_params': conn_params,
+                  'neuron_params': neuron_params,
+                  'input_params': input_params}
+
+params_LA = (network_params, copy.deepcopy(sim_params))
+
+"""
+Simulation with kappa = 1.125 leading to the high-activity fixed point
+shown in Fig. 2.
+"""
+sim_params.update({'t_sim': 10500.})
+conn_params = {'g': -16.,
+               'fac_nu_ext_TH': 1.2,
+               'fac_nu_ext_5E': 1.125,
+               'fac_nu_ext_6E': 1.41666667,
+               'av_indegree_V1': 3950.}
+network_params = {'N_scaling': 1.,
+                  'K_scaling': 1.,
+                  'connection_params': conn_params,
+                  'neuron_params': neuron_params}
+
+params_HA = (network_params, copy.deepcopy(sim_params))
+
+
+params_stab = {}
+"""
+Simulation with kappa = 1.125, chi=1
+Presented in Fig. 2, 4, and 8.
+"""
+sim_params.update({'t_sim': 10500.})
+conn_params = {'g': -11.,
+               'fac_nu_ext_TH': 1.2,
+               'fac_nu_ext_5E': 1.125,
+               'fac_nu_ext_6E': 1.41666667,
+               'av_indegree_V1': 3950.,
+               'K_stable': '../SchueckerSchmidt2017/K_prime_original.npy'}
+network_params = {'N_scaling': 1.,
+                  'K_scaling': 1.,
+                  'connection_params': conn_params,
+                  'neuron_params': neuron_params}
+params_stab[1.0] = (network_params, copy.deepcopy(sim_params))
+
+"""
+Simulation with kappa = 1.125, chi=1.9
+Presented in Fig. 4 and all following.
+
+One simulation with t_sim=10500. ms to plot spike raster plots (to
+reduce the computational load when loading spikes) and one with
+t_sim=100500. ms for all quantitative measures.
+"""
+sim_params.update({'t_sim': 10500.})
+conn_params = {'g': -11.,
+               'fac_nu_ext_TH': 1.2,
+               'fac_nu_ext_5E': 1.125,
+               'fac_nu_ext_6E': 1.41666667,
+               'av_indegree_V1': 3950.,
+               'K_stable': '../SchueckerSchmidt2017/K_prime_original.npy',
+               'cc_weights_factor': 1.9,
+               'cc_weights_I_factor': 2.}
+network_params = {'N_scaling': 1.,
+                  'K_scaling': 1.,
+                  'connection_params': conn_params,
+                  'neuron_params': neuron_params}
+params_stab['1.9_spikes'] = (network_params, copy.deepcopy(sim_params))
+
+sim_params.update({'t_sim': 100500.})
+conn_params = {'g': -11.,
+               'fac_nu_ext_TH': 1.2,
+               'fac_nu_ext_5E': 1.125,
+               'fac_nu_ext_6E': 1.41666667,
+               'av_indegree_V1': 3950.,
+               'K_stable': '../SchueckerSchmidt2017/K_prime_original.npy',
+               'cc_weights_factor': 1.9,
+               'cc_weights_I_factor': 2.}
+network_params = {'N_scaling': 1.,
+                  'K_scaling': 1.,
+                  'connection_params': conn_params,
+                  'neuron_params': neuron_params}
+params_stab[1.9] = (network_params, copy.deepcopy(sim_params))
+
+
+"""
+Simulations with kappa = 1.125
+and varying chi, presented in Fig. 4 and 8.
+"""
+cc_weights_factor_list = [1.4, 1.5, 1.6, 1.7, 1.75, 1.8, 2., 2.1, 2.5]
+for cc_weights_factor in cc_weights_factor_list:
+    sim_params.update({'t_sim': 10500.})
+    conn_params = {'g': -11.,
+                   'fac_nu_ext_TH': 1.2,
+                   'fac_nu_ext_5E': 1.125,
+                   'fac_nu_ext_6E': 1.41666667,
+                   'av_indegree_V1': 3950.,
+                   'K_stable': '../SchueckerSchmidt2017/K_prime_original.npy',
+                   'cc_weights_factor': cc_weights_factor,
+                   'cc_weights_I_factor': 2.}
+
+    network_params = {'N_scaling': 1.,
+                      'K_scaling': 1.,
+                      'connection_params': conn_params,
+                      'neuron_params': neuron_params}
+
+    params_stab[cc_weights_factor] = (network_params, copy.deepcopy(sim_params))
+
+"""
+Collect all parameter dictionaries in one place
+"""
+NEW_SIM_PARAMS = {'all': [params_LA,
+                          params_HA]
+                  + [params_stab[chi] for chi in params_stab],
+                  'Fig1': None,
+                  'Fig2': [params_LA,
+                           params_HA,
+                           params_stab[1.0]],
+                  'Fig3': [params_stab[1.0]],
+                  'Fig4': [params_stab[chi] for chi in [1., 1.8, 1.9, 2., 2.1, 2.5]],
+                  'Fig5': [params_stab['1.9_spikes'],
+                           params_stab[1.9]],
+                  'Fig6': [params_stab[chi] for chi in [1., 2.5, '1.9_spikes', 1.9]],
+                  'Fig7': [params_stab[1.9]],
+                  'Fig8': [params_stab[chi] for chi in [1., 1.4, 1.5, 1.6, 1.7, 1.75,
+                                                        1.8, 2., 2.1, 2.5, 1.9]],
+                  'Fig9': [params_stab[1.9]]}
+
+"""
+Collect all labels in one dictionary
+"""
+
+
+def init_model(par):
+    return MultiAreaModel(par[0],
+                          simulation=True,
+                          sim_spec=par[1])
+
+
+def init_models(figure):
+    """
+    Create instances of the models required for the given figure.
+    """
+    models = []
+    for par in NEW_SIM_PARAMS[figure]:
+        M = init_model(par)
+        models.append(M)
+    return models
+
+
+def create_label_dict():
+    M_LA = init_model(NEW_SIM_PARAMS['all'][0])
+    M_HA = init_model(NEW_SIM_PARAMS['all'][1])
+    M_stab = {chi: init_model(params_stab[chi]) for chi in params_stab}
+    NEW_SIM_LABELS = {'all': [M_LA.simulation.label,
+                              M_HA.simulation.label]
+                      + [M_stab[chi].simulation.label for chi in M_stab.keys()],
+                      'Fig1': None,
+                      'Fig2': [M_LA.simulation.label,
+                               M_HA.simulation.label,
+                               M_stab[1.0].simulation.label],
+                      'Fig3': [M_stab[1.0].simulation.label],
+                      'Fig4': [M_stab[chi].simulation.label for chi in [1., 1.8, 1.9,
+                                                                        2., 2.1, 2.5]],
+                      'Fig5': [M_stab['1.9_spikes'].simulation.label,
+                               M_stab[1.9].simulation.label],
+                      'Fig6': [M_stab[chi].simulation.label for chi in [1., 2.5, 1.9,
+                                                                        '1.9_spikes']],
+                      'Fig7': [M_stab[1.9].simulation.label],
+                      'Fig8': [M_stab[chi].simulation.label for chi in [1., 1.4, 1.5, 1.6,
+                                                                        1.7, 1.75, 1.8, 2.,
+                                                                        2.1, 2.5, 1.9]],
+                      'Fig9': [M_stab[1.9].simulation.label]}
+    return NEW_SIM_LABELS
+
+
+def run_simulation(figure):
+    models = init_models(figure)
+    for M in models:
+        start_job(M.simulation.label, submit_cmd, jobscript_template)
+    
+    
diff --git a/figures/Schmidt2018_dyn/run_network_simulations.py b/figures/Schmidt2018_dyn/run_network_simulations.py
deleted file mode 100644
index a04938c2e3d2bc6bf4b095aaa252d37a5d2553f8..0000000000000000000000000000000000000000
--- a/figures/Schmidt2018_dyn/run_network_simulations.py
+++ /dev/null
@@ -1,141 +0,0 @@
-import os
-
-from multiarea_model import MultiAreaModel
-from start_jobs import start_job
-from config import submit_cmd, jobscript_template
-from config import base_path
-
-"""
-This script provides the code to execute all simulations presented in
-Schmidt M, Bakker R, Shen K, Bezgin B, Diesmann M & van Albada SJ
-(2018) A multi-scale layer-resolved spiking network model of
-resting-state dynamics in macaque cortex.
-
-Needs to be simulated with sufficient
-resources, for instance on a compute cluster.
-"""
-
-# Common parameter settings
-input_params = {'rate_ext': 10.}
-neuron_params = {'V0_mean': -150.,
-                 'V0_sd': 50.}
-sim_params = {'num_processes': 20,  # Needs to be adapted to the HPC system used
-              'local_num_threads': 24,  # Needs to be adapted to the HPC system used
-              'recording_dict': {'record_vm': False}}
-
-
-"""
-Simulation with  kappa = 1. leading to the low-activity fixed point
-shown in Fig. 2.
-"""
-d = {}
-
-sim_params.update({'t_sim': 10500.})
-conn_params = {'g': -16.,
-               'fac_nu_ext_TH': 1.2,
-               'fac_nu_ext_5E': 1.,
-               'fac_nu_ext_6E': 1.,
-               'av_indegree_V1': 3950.}
-network_params = {'N_scaling': 1.,
-                  'K_scaling': 1.,
-                  'connection_params': conn_params,
-                  'neuron_params': neuron_params,
-                  'input_params': input_params}
-
-M_LA = MultiAreaModel(network_params, simulation=True,
-                      sim_spec=sim_params)
-
-# start_job(M_LA.simulation.label, submit_cmd, jobscript_template)
-
-
-"""
-Simulation with kappa = 1.125 leading to the high-activity fixed point
-shown in Fig. 2.
-"""
-sim_params.update({'t_sim': 10500.})
-conn_params = {'g': -16.,
-               'fac_nu_ext_TH': 1.2,
-               'fac_nu_ext_5E': 1.125,
-               'fac_nu_ext_6E': 1.41666667,
-               'av_indegree_V1': 3950.}
-network_params = {'N_scaling': 1.,
-                  'K_scaling': 1.,
-                  'connection_params': conn_params,
-                  'neuron_params': neuron_params}
-
-M_HA = MultiAreaModel(network_params, simulation=True,
-                      sim_spec=sim_params)
-p, r_HA = M_HA.theory.integrate_siegert()
-# start_job(M_HA.simulation.label, submit_cmd, jobscript_template)
-
-
-"""
-Simulation with kappa = 1.125, chi=1
-Presented in Fig. 2, 4, and 8.
-"""
-sim_params.update({'t_sim': 10500.})
-conn_params = {'g': -11.,
-               'fac_nu_ext_TH': 1.2,
-               'fac_nu_ext_5E': 1.125,
-               'fac_nu_ext_6E': 1.41666667,
-               'av_indegree_V1': 3950.,
-               'K_stable': '../SchueckerSchmidt2018/K_prime_original.npy'}
-network_params = {'N_scaling': 1.,
-                  'K_scaling': 1.,
-                  'connection_params': conn_params,
-                  'neuron_params': neuron_params}
-
-M_stab = MultiAreaModel(network_params, simulation=True,
-                        sim_spec=sim_params)
-p, r_stab = M_stab.theory.integrate_siegert()
-# start_job(M_stab.simulation.label, submit_cmd, jobscript_template)
-
-
-"""
-Simulation with kappa = 1.125, chi=1.9
-Presented in Fig. 4 and all following.
-"""
-sim_params.update({'t_sim': 100500.})
-conn_params = {'g': -11.,
-               'fac_nu_ext_TH': 1.2,
-               'fac_nu_ext_5E': 1.125,
-               'fac_nu_ext_6E': 1.41666667,
-               'av_indegree_V1': 3950.,
-               'K_stable': '../SchueckerSchmidt2018/K_prime_original.npy',
-               'cc_weights_factor': 1.9,
-               'cc_weights_I_factor': 2.}
-network_params = {'N_scaling': 1.,
-                  'K_scaling': 1.,
-                  'connection_params': conn_params,
-                  'neuron_params': neuron_params}
-
-M_stab = MultiAreaModel(network_params, simulation=True,
-                        sim_spec=sim_params)
-p, r_stab = M_stab.theory.integrate_siegert()
-# start_job(M_stab.simulation.label, submit_cmd, jobscript_template)
-
-
-"""
-Simulations with kappa = 1.125
-and varying chi, presented in Fig. 4 and 8.
-"""
-for cc_weights_factor in [1.4, 1.5, 1.6, 1.7, 1.75, 1.8, 2., 2.1, 2.5]:
-    sim_params.update({'t_sim': 10500.})
-    conn_params = {'g': -11.,
-                   'fac_nu_ext_TH': 1.2,
-                   'fac_nu_ext_5E': 1.125,
-                   'fac_nu_ext_6E': 1.41666667,
-                   'av_indegree_V1': 3950.,
-                   'K_stable': '../SchueckerSchmidt2018/K_prime_original.npy',
-                   'cc_weights_factor': cc_weights_factor,
-                   'cc_weights_I_factor': 2.}
-
-    network_params = {'N_scaling': 1.,
-                      'K_scaling': 1.,
-                      'connection_params': conn_params,
-                      'neuron_params': neuron_params}
-
-    M_stab = MultiAreaModel(network_params, simulation=True,
-                            sim_spec=sim_params)
-    p, r_stab = M_stab.theory.integrate_siegert()
-    # start_job(M_stab.simulation.label, submit_cmd, jobscript_template)