diff --git a/figures/Schmidt2018_dyn/compute_bold_signal.py b/figures/Schmidt2018_dyn/compute_bold_signal.py
index a10f5a59501fe13557fa38348b0ccf82385a91cd..ac72a790db177e98733d05980fb3da277e17d855 100644
--- a/figures/Schmidt2018_dyn/compute_bold_signal.py
+++ b/figures/Schmidt2018_dyn/compute_bold_signal.py
@@ -3,6 +3,12 @@ import os
 import sys
 
 
+"""
+Compute BOLD signal for a given area from the time series of
+population-averaged spike rates of a given simulation using the
+neuRosim package of R (see Schmidt et al. 2018 for more details).
+"""
+
 data_path = sys.argv[1]
 label = sys.argv[2]
 area = sys.argv[3]
diff --git a/figures/Schmidt2018_dyn/compute_corrcoeff.py b/figures/Schmidt2018_dyn/compute_corrcoeff.py
index 34ced53c77f54bd56c553768c3b6d9e78778c6c9..4bd7f889d3be4fed38bf7ac0258a0f83b4046858 100644
--- a/figures/Schmidt2018_dyn/compute_corrcoeff.py
+++ b/figures/Schmidt2018_dyn/compute_corrcoeff.py
@@ -6,6 +6,11 @@ import correlation_toolbox.helper as ch
 from multiarea_model import MultiAreaModel
 import sys
 
+"""
+Compute correlation coefficients for a subsample
+of neurons for the entire network from raw spike files of a given simulation.
+"""
+
 data_path = sys.argv[1]
 label = sys.argv[2]
 
@@ -31,7 +36,7 @@ M = MultiAreaModel({})
 
 spike_data = {}
 cc_dict = {}
-for area in ['V1', 'V2', 'FEF']:
+for area in M.area_list:
     cc_dict[area] = {}
     LvR_list = []
     N = []
diff --git a/figures/Schmidt2018_dyn/compute_cross_correlation.py b/figures/Schmidt2018_dyn/compute_cross_correlation.py
index fe607f781acb6cbb1b3d5d7dc79db703576f0807..905ebb98c08f08fb7769ed211b47abaeff6f8951 100644
--- a/figures/Schmidt2018_dyn/compute_cross_correlation.py
+++ b/figures/Schmidt2018_dyn/compute_cross_correlation.py
@@ -5,6 +5,11 @@ import numpy as np
 import os
 import sys
 
+""" 
+Compute the cross-correlation betwen two given areas from their
+time series of population-averaged spike rates of a given simulation.
+"""
+
 data_path = sys.argv[1]
 label = sys.argv[2]
 area1 = sys.argv[3]
diff --git a/figures/Schmidt2018_dyn/compute_functional_connectivity.py b/figures/Schmidt2018_dyn/compute_functional_connectivity.py
index 87eb4872b56285f193f01179fc5c688683291a83..14afe4cc9c8387e3fff99d1ffb7a42f2f33fd621 100644
--- a/figures/Schmidt2018_dyn/compute_functional_connectivity.py
+++ b/figures/Schmidt2018_dyn/compute_functional_connectivity.py
@@ -7,6 +7,12 @@ from multiarea_model import MultiAreaModel
 from scipy.spatial.distance import pdist
 from scipy.spatial.distance import squareform
 
+"""
+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]
diff --git a/figures/Schmidt2018_dyn/compute_granger_causality.py b/figures/Schmidt2018_dyn/compute_granger_causality.py
index b93d17d3f7b8c6b2b7b29fa9c94f3b85acb5bd78..e4575826e9bf083a7fc50a8da1fbd1b923e062cd 100644
--- a/figures/Schmidt2018_dyn/compute_granger_causality.py
+++ b/figures/Schmidt2018_dyn/compute_granger_causality.py
@@ -9,6 +9,13 @@ from multiarea_model.multiarea_helpers import create_mask
 from scipy.stats import levene
 from statsmodels.tsa.vector_ar.var_model import VAR
 
+
+"""
+Compute the conditional Granger causality to a given population of an
+area based on the population-averaged spike rates from a given
+simulation.
+"""
+
 data_path = sys.argv[1]
 label = sys.argv[2]
 area = sys.argv[3]
diff --git a/figures/Schmidt2018_dyn/compute_louvain_communities.py b/figures/Schmidt2018_dyn/compute_louvain_communities.py
index f3913fd04c00612c7b2591d2141d6bc14b01f592..5d558723fe6c27480351ee1d24fb00e916ceed9e 100644
--- a/figures/Schmidt2018_dyn/compute_louvain_communities.py
+++ b/figures/Schmidt2018_dyn/compute_louvain_communities.py
@@ -8,6 +8,13 @@ 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]
diff --git a/figures/Schmidt2018_dyn/compute_pop_LvR.py b/figures/Schmidt2018_dyn/compute_pop_LvR.py
index 9128ffd93572b0f19b6f2ce93a0ac98d6e53ac3c..cd0d0704838ec0003f2fa9c67b7681a66dec866b 100644
--- a/figures/Schmidt2018_dyn/compute_pop_LvR.py
+++ b/figures/Schmidt2018_dyn/compute_pop_LvR.py
@@ -6,6 +6,11 @@ from multiarea_model.analysis_helpers import pop_LvR
 from multiarea_model import MultiAreaModel
 import sys
 
+"""
+Compute LvR for the entire network from raw spike
+files of a given simulation.
+"""
+
 data_path = sys.argv[1]
 label = sys.argv[2]
 
diff --git a/figures/Schmidt2018_dyn/compute_pop_rates.py b/figures/Schmidt2018_dyn/compute_pop_rates.py
index 784f9a00e4345eab71b63ca2c998f0d4ede7457b..88a9945bb783e75dab5e77db4b010b04aa3d0f16 100644
--- a/figures/Schmidt2018_dyn/compute_pop_rates.py
+++ b/figures/Schmidt2018_dyn/compute_pop_rates.py
@@ -6,6 +6,11 @@ from multiarea_model.analysis_helpers import pop_rate
 from multiarea_model import MultiAreaModel
 import sys
 
+"""
+Compute stationary spike rates for the entire network from raw spike
+files of a given simulation.
+"""
+
 data_path = sys.argv[1]
 label = sys.argv[2]
 
@@ -27,7 +32,7 @@ M = MultiAreaModel({})
 
 spike_data = {}
 pop_rates = {}
-for area in ['V1', 'V2', 'FEF']:
+for area in M.area_list:
     pop_rates[area] = {}
     rate_list = []
     N = []
diff --git a/figures/Schmidt2018_dyn/compute_power_spectrum.py b/figures/Schmidt2018_dyn/compute_power_spectrum.py
index 3cba554db5ed1dff93dd42ba95203a95f4c1e07c..3555077596fa8b4eb44268dc425f68bbf0e4e568 100644
--- a/figures/Schmidt2018_dyn/compute_power_spectrum.py
+++ b/figures/Schmidt2018_dyn/compute_power_spectrum.py
@@ -6,6 +6,18 @@ from multiarea_model import MultiAreaModel
 from multiarea_model.analysis_helpers import centralize
 from scipy.signal import welch
 
+
+"""
+Compute the power spectrum time series for a given area of
+population-averaged spike rates of a given simulation.
+
+The spike rates can be based on three different methods:
+- binned spike histograms on all neurons ('full')
+- binned spike histograms on a subsample of 140 neurons ('subsample')
+- spike histograms convolved with a Gaussian kernel of optimal width
+  after Shimazaki et al. (2010)
+"""
+
 # Parameters for Welch Power Spectral density and spectrogram
 noverlap = 1000
 nperseg = 1024
diff --git a/figures/Schmidt2018_dyn/compute_rate_histogram.py b/figures/Schmidt2018_dyn/compute_rate_histogram.py
index c351b3360a9dfd494f3f43806ac68597a6a0542e..6e5dec02f045e299b728a1b1f104d95434050263 100644
--- a/figures/Schmidt2018_dyn/compute_rate_histogram.py
+++ b/figures/Schmidt2018_dyn/compute_rate_histogram.py
@@ -6,6 +6,12 @@ import sys
 from multiarea_model import MultiAreaModel
 from multiarea_model.analysis_helpers import pop_rate_distribution
 
+
+"""
+Compute histogram of spike rates over single neurons for a given area
+from raw spike files of a given simulation.
+"""
+
 assert(len(sys.argv) == 4)
 data_path = sys.argv[1]
 label = sys.argv[2]
diff --git a/figures/Schmidt2018_dyn/compute_rate_time_series.py b/figures/Schmidt2018_dyn/compute_rate_time_series.py
index 4a2ebfcc1d82a335955d4b72aa5a175ade5ab54e..c4956af202d6f155601762c99fbf623ac4ad373a 100644
--- a/figures/Schmidt2018_dyn/compute_rate_time_series.py
+++ b/figures/Schmidt2018_dyn/compute_rate_time_series.py
@@ -2,7 +2,6 @@ import json
 import neo
 import numpy as np
 import os
-import pandas as pd
 import quantities as pq
 
 from multiarea_model.analysis_helpers import pop_rate_time_series
@@ -10,6 +9,17 @@ from elephant.statistics import instantaneous_rate
 from multiarea_model import MultiAreaModel
 import sys
 
+"""
+Compute time series of population-averaged spike rates for a given
+area from raw spike files of a given simulation.
+
+Implements three different methods:
+- binned spike histograms on all neurons ('full')
+- binned spike histograms on a subsample of 140 neurons ('subsample')
+- spike histograms convolved with a Gaussian kernel of optimal width
+  after Shimazaki et al. (2010)
+"""
+
 assert(len(sys.argv) == 5)
 
 data_path = sys.argv[1]