diff --git a/.ipynb_checkpoints/multi-area-model-checkpoint.ipynb b/.ipynb_checkpoints/multi-area-model-checkpoint.ipynb
index ab7ae2504ea3bb423670a84480ecde42ea785b26..e198e1b40beb8f45adf43cca81399c92105a1689 100644
--- a/.ipynb_checkpoints/multi-area-model-checkpoint.ipynb
+++ b/.ipynb_checkpoints/multi-area-model-checkpoint.ipynb
@@ -25,14 +25,16 @@
     "    * [2.2. Predict firing rates from theory](#section_2_2)\n",
     "    * [2.3. Extract interarea connectivity](#section_2_3)\n",
     "    * [2.4. Run the simulation](#section_2_4)\n",
-    "* [S3. Simulation results analysis and data processing](#section_3)\n",
-    "* [S4. Simulation results visualization](#section_4) \n",
-    "    * [4.1. Instantaneous firing rate and mean firing rate](#section_4_1)\n",
-    "    * [4.2. Raster plot of spiking activity for single area](#section_4_2)\n",
-    "    * [4.3. Population-averaged firing rate](#section_4_3)\n",
-    "    * [4.4. Average pairwise correlation coefficients of spiking activity](#section_4_4)\n",
-    "    * [4.5. Irregularity of spiking activity](#section_4_5)\n",
-    "    * [4.6. Time series of population- and area-averaged firing rates](#section_4_6)"
+    "* [S3. Simulation results validation and connection extraction](#section_3)\n",
+    "* [S4. Data loading and processing](#section_4)\n",
+    "* [S5. Simulation results visualization](#section_5) \n",
+    "    * [5.1. Instantaneous firing rate and mean firing rate across all populations](#section_5_1)\n",
+    "    * [5.2. Raster plot of spiking activity for single area](#section_5_2)\n",
+    "    * [5.3. Population-averaged firing rate](#section_5_3)\n",
+    "    * [5.4 Time-averaged population rates](#section_5_4)\n",
+    "    * [5.5. Average pairwise correlation coefficients of spiking activity](#section_5_5)\n",
+    "    * [5.6. Irregularity of spiking activity](#section_5_6)\n",
+    "    * [5.7. Time series of population- and area-averaged firing rates](#section_5_7)"
    ]
   },
   {
@@ -54,6 +56,24 @@
     "## S0. Configuration <a class=\"anchor\" id=\"section_0\"></a>"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "06e7f341-c226-4782-a913-4868027d7d06",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Create config file\n",
+    "with open('config.py', 'w') as fp:\n",
+    "    fp.write(\n",
+    "'''import os\n",
+    "base_path = os.path.abspath(\".\")\n",
+    "data_path = os.path.abspath(\"simulations\")\n",
+    "jobscript_template = \"python {base_path}/run_simulation.py {label}\"\n",
+    "submit_cmd = \"bash -c\"\n",
+    "''')"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -74,25 +94,7 @@
     "# Import the MultiAreaModel class\n",
     "from multiarea_model import MultiAreaModel\n",
     "from multiarea_model import Analysis\n",
-    "from config import base_path"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "2dd47c64",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# Create config file\n",
-    "with open('config.py', 'w') as fp:\n",
-    "    fp.write(\n",
-    "'''import os\n",
-    "base_path = os.path.abspath(\".\")\n",
-    "data_path = os.path.abspath(\"simulations\")\n",
-    "jobscript_template = \"python {base_path}/run_simulation.py {label}\"\n",
-    "submit_cmd = \"bash -c\"\n",
-    "''')"
+    "from config import base_path, data_path"
    ]
   },
   {
@@ -301,7 +303,10 @@
   {
    "cell_type": "markdown",
    "id": "1fd58841",
-   "metadata": {},
+   "metadata": {
+    "jp-MarkdownHeadingCollapsed": true,
+    "tags": []
+   },
    "source": [
     "### 2.1. Insantiate a multi-area model <a class=\"anchor\" id=\"section_2_1\"></a>"
    ]
@@ -455,7 +460,7 @@
     "tags": []
    },
    "source": [
-    "## S3. Simulation results analysis <a class=\"anchor\" id=\"section_3\"></a>"
+    "## S3. Simulation results validation and connection extraction <a class=\"anchor\" id=\"section_3\"></a>"
    ]
   },
   {
@@ -558,19 +563,21 @@
   {
    "cell_type": "markdown",
    "id": "57ff902c-d6ce-4f96-9e4f-8e3e7166ab66",
-   "metadata": {},
+   "metadata": {
+    "jp-MarkdownHeadingCollapsed": true,
+    "tags": []
+   },
    "source": [
-    "## S4. Data processing and simulation results visualization <a class=\"anchor\" id=\"section_4\"></a>"
+    "## S4. Data loading and processing <a class=\"anchor\" id=\"section_4\"></a>"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "6806564b-d5d4-47d4-afa9-5da0df83e025",
+   "id": "f5b58845-4d1a-430f-83f4-402fdf918aef",
    "metadata": {},
    "outputs": [],
    "source": [
-    "from config import data_path\n",
     "label_spikes = M.simulation.label\n",
     "label = M.simulation.label"
    ]
@@ -578,10 +585,8 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "7997d893-252d-4295-a22a-1e510f8424ae",
-   "metadata": {
-    "tags": []
-   },
+   "id": "6607a73d-1c74-4848-9603-081ad0e7cae8",
+   "metadata": {},
    "outputs": [],
    "source": [
     "\"\"\"\n",
@@ -604,6 +609,7 @@
     "    Default value is None and leads to loading of data for all\n",
     "    simulated areas.\n",
     "\"\"\"\n",
+    "# Instantiate an analysis class and load spike data\n",
     "A = Analysis(network=M, \n",
     "             simulation=M.simulation, \n",
     "             data_list=['spikes'],\n",
@@ -613,350 +619,300 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "da58921f-713b-424c-8fa1-80d9755558f3",
+   "id": "c471e9c8-b1e4-43e4-a6a1-443b8b8963be",
    "metadata": {},
    "outputs": [],
    "source": [
-    "\"\"\"\n",
-    "Loads simulation data of the requested type either from hdf5 files.\n",
-    "\n",
-    "Parameters\n",
-    "----------\n",
-    "\n",
-    "data_list : list\n",
-    "    list of observables to be loaded. Can contain 'spikes' and 'vm'\n",
-    "\"\"\"\n",
-    "A.load_data(data_list=['spikes'])"
+    "data = np.loadtxt(M.simulation.data_dir + '/recordings/' + M.simulation.label + \"-spikes-1-0.dat\", skiprows=3)\n",
+    "tsteps, spikecount = np.unique(data[:,1], return_counts=True)\n",
+    "rate = spikecount / M.simulation.params['dt'] * 1e3 / np.sum(M.N_vec)"
    ]
   },
   {
-   "cell_type": "markdown",
-   "id": "38ddd973",
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1870cf34-ee62-4614-bc25-c36bc9a7377c",
    "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
     "tags": []
    },
+   "outputs": [],
    "source": [
-    "### 4.1. Instantaneous firing rate and mean firing rate <a class=\"anchor\" id=\"section_4_1\"></a>"
+    "\"\"\"\n",
+    "Calculate time-averaged population rates and store them in member pop_rates.\n",
+    "If the rates had previously been stored with the same\n",
+    "parameters, they are loaded from file.\n",
+    "\n",
+    "Parameters\n",
+    "----------\n",
+    "t_min : float, optional\n",
+    "    Minimal time in ms of the simulation to take into account\n",
+    "    for the calculation. Defaults to 500 ms.\n",
+    "t_max : float, optional\n",
+    "    Maximal time in ms of the simulation to take into account\n",
+    "    for the calculation. Defaults to the simulation time.\n",
+    "compute_stat : bool, optional\n",
+    "    If set to true, the mean and variance of the population rate\n",
+    "    is calculated. Defaults to False.\n",
+    "    Caution: Setting to True slows down the computation.\n",
+    "areas : list, optional\n",
+    "    Which areas to include in the calculcation.\n",
+    "    Defaults to all loaded areas.\n",
+    "pops : list or {'complete'}, optional\n",
+    "    Which populations to include in the calculation.\n",
+    "    If set to 'complete', all populations the respective areas\n",
+    "    are included. Defaults to 'complete'.\n",
+    "\"\"\"\n",
+    "A.create_pop_rates()"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "56e368b6-72a2-43fb-b02e-f70ad6770e40",
+   "id": "50b7df89",
    "metadata": {},
    "outputs": [],
    "source": [
-    "data = np.loadtxt(M.simulation.data_dir + '/recordings/' + M.simulation.label + \"-spikes-1-0.dat\", skiprows=3)\n",
-    "tsteps, spikecount = np.unique(data[:,1], return_counts=True)\n",
-    "rate = spikecount / M.simulation.params['dt'] * 1e3 / np.sum(M.N_vec)"
+    "\"\"\"\n",
+    "Calculate synchrony as the coefficient of variation of the population rate\n",
+    "and store in member synchrony. Uses helper function synchrony.\n",
+    "If the synchrony has previously been stored with the\n",
+    "same parameters, they are loaded from file.\n",
+    "\n",
+    "\n",
+    "Parameters\n",
+    "----------\n",
+    "t_min : float, optional\n",
+    "    Minimal time in ms of the simulation to take into account\n",
+    "    for the calculation. Defaults to 500 ms.\n",
+    "t_max : float, optional\n",
+    "    Maximal time in ms of the simulation to take into account\n",
+    "    for the calculation. Defaults to the simulation time.\n",
+    "areas : list, optional\n",
+    "    Which areas to include in the calculcation.\n",
+    "    Defaults to all loaded areas.\n",
+    "pops : list or {'complete'}, optional\n",
+    "    Which populations to include in the calculation.\n",
+    "    If set to 'complete', all populations the respective areas\n",
+    "    are included. Defaults to 'complete'.\n",
+    "resolution : float, optional\n",
+    "    Resolution of the population rate. Defaults to 1 ms.\n",
+    "\"\"\"\n",
+    "A.create_synchrony()"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "bea30fc8",
+   "id": "d43b493c",
    "metadata": {},
    "outputs": [],
    "source": [
-    "fig, ax = plt.subplots()\n",
-    "ax.plot(tsteps, rate)\n",
-    "ax.plot(tsteps, np.average(rate)*np.ones(len(tsteps)), label='mean')\n",
-    "ax.set_title('Instantaneous firing rate across all populations')\n",
-    "ax.set_xlabel('time (ms)')\n",
-    "ax.set_ylabel('Firing rate (spikes / s)')\n",
-    "ax.set_xlim(0, sim_params['t_sim'])\n",
-    "ax.set_ylim(0, 50)\n",
-    "ax.legend()"
+    "\"\"\"\n",
+    "Calculate poulation-averaged LvR (see Shinomoto et al. 2009) and\n",
+    "store as member pop_LvR. Uses helper function LvR.\n",
+    "\n",
+    "Parameters\n",
+    "----------\n",
+    "t_min : float, optional\n",
+    "    Minimal time in ms of the simulation to take into account\n",
+    "    for the calculation. Defaults to 500 ms.\n",
+    "t_max : float, optional\n",
+    "    Maximal time in ms of the simulation to take into account\n",
+    "    for the calculation. Defaults to the simulation time.\n",
+    "areas : list, optional\n",
+    "    Which areas to include in the calculcation.\n",
+    "    Defaults to all loaded areas.\n",
+    "pops : list or {'complete'}, optional\n",
+    "    Which populations to include in the calculation.\n",
+    "    If set to 'complete', all populations the respective areas\n",
+    "    are included. Defaults to 'complete'.\n",
+    "\"\"\"\n",
+    "A.create_pop_LvR()"
    ]
   },
   {
-   "cell_type": "markdown",
-   "id": "ae19bcc3",
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "401ece2d-47c8-4775-80ae-92a8e432520c",
    "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
     "tags": []
    },
-   "source": [
-    "### 4.2 Raster plot of spiking activity for single area <a class=\"anchor\" id=\"section_4_2\"></a>\n",
-    "Raster plot of spiking activity of 3% of the neurons in area V1 (A), V2 (B), and FEF (C). Blue: excitatory neurons, red: inhibitory neurons. (D-F) Spiking statistics across all 32 areas for the respective populations shown as area-averaged box plots. Crosses: medians, boxes: interquartile range (IQR), whiskers extend to the most extremeobservat ions within 1.5×IQR beyond the IQR."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "c47e82ec-32f3-4de1-a32e-8f6b500787e0",
-   "metadata": {},
    "outputs": [],
    "source": [
-    "areas = ['V1', 'V2', 'FEF']\n",
-    "labels = ['A', 'B', 'C']"
+    "\"\"\"\n",
+    "Calculate time series of population- and area-averaged firing rates.\n",
+    "Uses ah.pop_rate_time_series.\n",
+    "If the rates have previously been stored with the\n",
+    "same parameters, they are loaded from file.\n",
+    "\n",
+    "\n",
+    "Parameters\n",
+    "----------\n",
+    "t_min : float, optional\n",
+    "    Minimal time in ms of the simulation to take into account\n",
+    "    for the calculation. Defaults to 500 ms.\n",
+    "t_max : float, optional\n",
+    "    Maximal time in ms of the simulation to take into account\n",
+    "    for the calculation. Defaults to the simulation time.\n",
+    "areas : list, optional\n",
+    "    Which areas to include in the calculcation.\n",
+    "    Defaults to all loaded areas.\n",
+    "pops : list or {'complete'}, optional\n",
+    "    Which populations to include in the calculation.\n",
+    "    If set to 'complete', all populations the respective areas\n",
+    "    are included. Defaults to 'complete'.\n",
+    "kernel : {'gauss_time_window', 'alpha_time_window', 'rect_time_window'}, optional\n",
+    "    Specifies the kernel to be convolved with the spike histogram.\n",
+    "    Defaults to 'binned', which corresponds to no convolution.\n",
+    "resolution: float, optional\n",
+    "    Width of the convolution kernel. Specifically it correponds to:\n",
+    "    - 'binned' : bin width of the histogram\n",
+    "    - 'gauss_time_window' : sigma\n",
+    "    - 'alpha_time_window' : time constant of the alpha function\n",
+    "    - 'rect_time_window' : width of the moving rectangular function\n",
+    "\"\"\"\n",
+    "A.create_rate_time_series()\n",
+    "A.save()"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "4066f042-995f-4987-bac2-7d7c61addbd0",
-   "metadata": {},
-   "outputs": [],
+   "cell_type": "markdown",
+   "id": "bb71c922",
+   "metadata": {
+    "jp-MarkdownHeadingCollapsed": true,
+    "tags": []
+   },
    "source": [
-    "# spike data  \n",
-    "spike_data = {}\n",
-    "for area in areas:\n",
-    "    spike_data[area] = {}\n",
-    "    for pop in M.structure[area]:\n",
-    "        spike_data[area][pop] = np.load(os.path.join(M.simulation.data_dir,\n",
-    "                                                     'recordings',\n",
-    "                                                     '{}-spikes-{}-{}.npy'.format(label_spikes, area, pop)))\n"
+    "## S5. Simulation results visualziation <a class=\"anchor\" id=\"section_5\"></a>"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "1da18fee",
-   "metadata": {},
-   "outputs": [],
+   "cell_type": "markdown",
+   "id": "38ddd973",
+   "metadata": {
+    "tags": []
+   },
    "source": [
-    "# \"\"\"\n",
-    "# Create raster display of a single area with populations stacked onto each other. Excitatory neurons in blue, inhibitory neurons in red.\n",
-    "\n",
-    "# Parameters\n",
-    "# ----------\n",
-    "# area : string {area}\n",
-    "#     Area to be plotted.\n",
-    "# frac_neurons : float, [0,1]\n",
-    "#     Fraction of cells to be considered.\n",
-    "# t_min : float, optional\n",
-    "#     Minimal time in ms of spikes to be shown. Defaults to 0 ms.\n",
-    "# t_max : float, optional\n",
-    "#     Minimal time in ms of spikes to be shown. Defaults to simulation time.\n",
-    "# output : {'pdf', 'png', 'eps'}, optional\n",
-    "#     If given, the function stores the plot to a file of the given format.\n",
-    "\n",
-    "# \"\"\"\n",
-    "# t_min = 0.\n",
-    "# t_max = 500.\n",
-    "\n",
-    "# # Draw V1\n",
-    "# area = 'V1'\n",
-    "# frac_neurons = 1.\n",
-    "# A.single_dot_display(area,  frac_neurons, t_min, t_max)\n",
-    "\n",
-    "# # Draw V2\n",
-    "# area = 'V2'\n",
-    "# frac_neurons = 1.\n",
-    "# A.single_dot_display(area,  frac_neurons, t_min, t_max)\n",
-    "\n",
-    "# # Draw FEF\n",
-    "# area = 'FEF'\n",
-    "# frac_neurons = 1.\n",
-    "# A.single_dot_display(area,  frac_neurons, t_min, t_max)"
+    "### 5.1. Instantaneous firing rate and mean firing rate across all populations <a class=\"anchor\" id=\"section_5_1\"></a>"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "73ca1021-0f0a-45bb-a3d5-51381d1357c3",
+   "id": "bea30fc8",
    "metadata": {},
    "outputs": [],
    "source": [
-    "for area, label in zip(areas, labels):\n",
-    "    label_pos = [-0.2, 1.01]\n",
-    "    pl.text(label_pos[0], label_pos[1], r'\\bfseries{}' + label + ': ' + area,\n",
-    "            fontdict={'fontsize': 10, 'weight': 'bold',\n",
-    "                      'horizontalalignment': 'left', 'verticalalignment':\n",
-    "                      'bottom'}, transform=axes[label].transAxes)\n",
-    "print(\"Raster plots\")\n",
-    "\n",
-    "t_min = 3000.\n",
-    "t_max = 3500.\n",
-    "\n",
-    "icolor = myred\n",
-    "ecolor = myblue\n",
-    "\n",
-    "frac_neurons = 0.03\n",
-    "\n",
-    "for i, area in enumerate(areas):\n",
-    "    ax = axes[labels[i]]\n",
-    "\n",
-    "    if area in spike_data:\n",
-    "        n_pops = len(spike_data[area])\n",
-    "        # Determine number of neurons that will be plotted for this area (for\n",
-    "        # vertical offset)\n",
-    "        offset = 0\n",
-    "        n_to_plot = {}\n",
-    "        for pop in M.structure[area]:\n",
-    "            n_to_plot[pop] = int(M.N[area][pop] * frac_neurons)\n",
-    "            offset = offset + n_to_plot[pop]\n",
-    "        y_max = offset + 1\n",
-    "        prev_pop = ''\n",
-    "        yticks = []\n",
-    "        yticklocs = []\n",
-    "        for jj, pop in enumerate(M.structure[area]):\n",
-    "            if pop[0:-1] != prev_pop:\n",
-    "                prev_pop = pop[0:-1]\n",
-    "                yticks.append('L' + population_labels[jj][0:-1])\n",
-    "                yticklocs.append(offset - 0.5 * n_to_plot[pop])\n",
-    "            ind = np.where(np.logical_and(\n",
-    "                spike_data[area][pop][:, 1] <= t_max, spike_data[area][pop][:, 1] >= t_min))\n",
-    "            pop_data = spike_data[area][pop][ind]\n",
-    "            pop_neurons = np.unique(pop_data[:, 0])\n",
-    "            neurons_to_ = np.arange(np.min(spike_data[area][pop][:, 0]), np.min(\n",
-    "                spike_data[area][pop][:, 0]) + n_to_plot[pop], 1)\n",
-    "\n",
-    "            if pop.find('E') > (-1):\n",
-    "                pcolor = ecolor\n",
-    "            else:\n",
-    "                pcolor = icolor\n",
-    "\n",
-    "            for kk in range(n_to_plot[pop]):\n",
-    "                spike_times = pop_data[pop_data[:, 0] == neurons_to_[kk], 1]\n",
-    "\n",
-    "                _ = ax.plot(spike_times, np.zeros(len(spike_times)) +\n",
-    "                            offset - kk, '.', color=pcolor, markersize=1)\n",
-    "            offset = offset - n_to_plot[pop]\n",
-    "        y_min = offset\n",
-    "        ax.set_xlim([t_min, t_max])\n",
-    "        ax.set_ylim([y_min, y_max])\n",
-    "        ax.set_yticklabels(yticks)\n",
-    "        ax.set_yticks(yticklocs)\n",
-    "        ax.set_xlabel('Time (s)', labelpad=-0.1)\n",
-    "        ax.set_xticks([t_min, t_min + 250., t_max])\n",
-    "        ax.set_xticklabels([r'$3.$', r'$3.25$', r'$3.5$'])"
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(tsteps, rate)\n",
+    "ax.plot(tsteps, np.average(rate)*np.ones(len(tsteps)), label='mean')\n",
+    "ax.set_title('Instantaneous firing rate across all populations')\n",
+    "ax.set_xlabel('time (ms)')\n",
+    "ax.set_ylabel('Firing rate (spikes / s)')\n",
+    "ax.set_xlim(0, sim_params['t_sim'])\n",
+    "ax.set_ylim(0, 50)\n",
+    "ax.legend()"
    ]
   },
   {
    "cell_type": "markdown",
-   "id": "019d805e",
+   "id": "ae19bcc3",
    "metadata": {
     "tags": []
    },
    "source": [
-    "### 4.3 Population-averaged firing rate <a class=\"anchor\" id=\"section_4_3\"></a>"
+    "### 5.2 Raster plot of spiking activity for single area <a class=\"anchor\" id=\"section_5_2\"></a>\n",
+    "Fig. 3 (A-C) Raster plot of spiking activity of 3% of the neurons in area V1 (A), V2 (B), and FEF (C). Blue: excitatory neurons, red: inhibitory neurons. (D-F) Spiking statistics across all 32 areas for the respective populations shown as area-averaged box plots. Crosses: medians, boxes: interquartile range (IQR), whiskers extend to the most extremeobservat ions within 1.5×IQR beyond the IQR."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "c05412f6-c842-415f-888a-b7604b795912",
+   "id": "1da18fee",
    "metadata": {},
    "outputs": [],
    "source": [
     "\"\"\"\n",
-    "Calculate time-averaged population rates and store them in member pop_rates.\n",
-    "If the rates had previously been stored with the same\n",
-    "parameters, they are loaded from file.\n",
+    "Create raster display of a single area with populations stacked onto each other. Excitatory neurons in blue, inhibitory neurons in red.\n",
     "\n",
     "Parameters\n",
     "----------\n",
+    "area : string {area}\n",
+    "    Area to be plotted.\n",
+    "frac_neurons : float, [0,1]\n",
+    "    Fraction of cells to be considered.\n",
     "t_min : float, optional\n",
-    "    Minimal time in ms of the simulation to take into account\n",
-    "    for the calculation. Defaults to 500 ms.\n",
+    "    Minimal time in ms of spikes to be shown. Defaults to 0 ms.\n",
     "t_max : float, optional\n",
-    "    Maximal time in ms of the simulation to take into account\n",
-    "    for the calculation. Defaults to the simulation time.\n",
-    "compute_stat : bool, optional\n",
-    "    If set to true, the mean and variance of the population rate\n",
-    "    is calculated. Defaults to False.\n",
-    "    Caution: Setting to True slows down the computation.\n",
-    "areas : list, optional\n",
-    "    Which areas to include in the calculcation.\n",
-    "    Defaults to all loaded areas.\n",
-    "pops : list or {'complete'}, optional\n",
-    "    Which populations to include in the calculation.\n",
-    "    If set to 'complete', all populations the respective areas\n",
-    "    are included. Defaults to 'complete'.\n",
+    "    Minimal time in ms of spikes to be shown. Defaults to simulation time.\n",
+    "output : {'pdf', 'png', 'eps'}, optional\n",
+    "    If given, the function stores the plot to a file of the given format.\n",
+    "\n",
     "\"\"\"\n",
-    "A.create_pop_rates(t_min=0)\n",
-    "print(\"Computing population rates done!\")"
+    "t_min = 0.\n",
+    "t_max = 500.\n",
+    "areas = ['V1', 'V2', 'FEF']\n",
+    "frac_neurons = 1.\n",
+    "for area in areas:\n",
+    "    A.single_dot_display(area,  frac_neurons, t_min, t_max)"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "721d1f03-df25-468d-8075-a807025a9c58",
+   "cell_type": "markdown",
+   "id": "019d805e",
+   "metadata": {
+    "tags": []
+   },
+   "source": [
+    "### 5.3 Population-averaged firing rate <a class=\"anchor\" id=\"section_5_3\"></a>\n",
+    "Fig 3. (D) Population-averaged firing rates"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "473d0882-8e45-4330-bfa2-2c7e1af0dac4",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "# stationary firing rates\n",
-    "fn = os.path.join(data_path, label, 'Analysis', 'pop_rates.json')\n",
-    "with open(fn, 'r') as f:\n",
-    "    pop_rates = json.load(f)"
+    "### 5.4 Time-averaged population rates <a class=\"anchor\" id=\"section_5_4\"></a>"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "9ba5ca35-7f90-47c8-a057-7cb02ee7be02",
+   "id": "721d1f03-df25-468d-8075-a807025a9c58",
    "metadata": {},
    "outputs": [],
    "source": [
-    "def set_boxplot_props(d):\n",
-    "    for i in range(len(d['boxes'])):\n",
-    "        if i % 2 == 0:\n",
-    "            d['boxes'][i].set_facecolor(icolor)\n",
-    "            d['boxes'][i].set_color(icolor)\n",
-    "        else:\n",
-    "            d['boxes'][i].set_facecolor(ecolor)\n",
-    "            d['boxes'][i].set_color(ecolor)\n",
-    "    pl.setp(d['whiskers'], color='k')\n",
-    "    pl.setp(d['fliers'], color='k', markerfacecolor='k', marker='+')\n",
-    "    pl.setp(d['medians'], color='none')\n",
-    "    pl.setp(d['caps'], color='k')\n",
-    "    pl.setp(d['means'], marker='x', color='k',\n",
-    "            markerfacecolor='k', markeredgecolor='k', markersize=3.)\n",
-    "    \n",
-    "print(\"plotting Population rates\")\n",
-    "\n",
-    "rates = np.zeros((len(M.area_list), 8))\n",
-    "for i, area in enumerate(M.area_list):\n",
-    "    for j, pop in enumerate(M.structure[area][::-1]):\n",
-    "        rate = pop_rates[area][pop][0]\n",
-    "        if rate == 0.0:\n",
-    "            rate = 1e-5\n",
-    "        if area == 'TH' and j > 3:  # To account for missing layer 4 in TH\n",
-    "            rates[i][j + 2] = rate\n",
-    "        else:\n",
-    "            rates[i][j] = rate\n",
-    "\n",
-    "\n",
-    "rates = np.transpose(rates)\n",
-    "masked_rates = np.ma.masked_where(rates < 1e-4, rates)\n",
-    "\n",
-    "ax = axes['D']\n",
-    "d = ax.boxplot(np.transpose(rates), vert=False,\n",
-    "               patch_artist=True, whis=1.5, showmeans=True)\n",
-    "set_boxplot_props(d)\n",
-    "\n",
-    "ax.plot(np.mean(rates, axis=1), np.arange(\n",
-    "    1., len(M.structure['V1']) + 1., 1.), 'x', color='k', markersize=3)\n",
-    "ax.set_yticklabels(population_labels[::-1], size=8)\n",
-    "ax.set_yticks(np.arange(1., len(M.structure['V1']) + 1., 1.))\n",
-    "ax.set_ylim((0., len(M.structure['V1']) + .5))\n",
+    "\"\"\"\n",
+    "Plot overview over time-averaged population rates encoded in colors\n",
+    "with areas along x-axis and populations along y-axis.\n",
     "\n",
-    "x_max = 220.\n",
-    "ax.set_xlim((-1., x_max))\n",
-    "ax.set_xlabel(r'Rate (spikes/s)', labelpad=-0.1)\n",
-    "ax.set_xticks([0., 50., 100.])"
+    "Parameters\n",
+    "----------\n",
+    "area_list : list, optional\n",
+    "   Specifies with areas are plotted in which order.\n",
+    "   Default to None, leading to plotting of  all areas ordered by architectural type.\n",
+    "output : {'pdf', 'png', 'eps'}, optional\n",
+    "    If given, the function stores the plot to a file of the given format.\n",
+    "\"\"\"\n",
+    "A.show_rates()"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "06a595de",
    "metadata": {
-    "jp-MarkdownHeadingCollapsed": true,
     "tags": []
    },
    "source": [
-    "### 4.4 Average pairwise correlation coefficients of spiking activity <a class=\"anchor\" id=\"section_4_4\"></a>"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "84d1689c",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "compute_corrcoeff.py"
+    "### 5.5 Average pairwise correlation coefficients of spiking activity <a class=\"anchor\" id=\"section_5_5\"></a>\n",
+    "Fig 5. (E) Average pairwise correlation coefficients of spiking activity"
    ]
   },
   {
@@ -1015,42 +971,11 @@
    "cell_type": "markdown",
    "id": "a3847e67",
    "metadata": {
-    "jp-MarkdownHeadingCollapsed": true,
     "tags": []
    },
    "source": [
-    "### 4.5 Irregularity of spiking activity <a class=\"anchor\" id=\"section_4_5\"></a>\n",
-    "Irregularity is measured by revised local variation LvR averaged across neurons"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "3c41c7d1-c39a-4c56-bda7-daa515cbaef7",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "\"\"\"\n",
-    "Calculate poulation-averaged LvR (see Shinomoto et al. 2009) and\n",
-    "store as member pop_LvR. Uses helper function LvR.\n",
-    "\n",
-    "Parameters\n",
-    "----------\n",
-    "t_min : float, optional\n",
-    "    Minimal time in ms of the simulation to take into account\n",
-    "    for the calculation. Defaults to 500 ms.\n",
-    "t_max : float, optional\n",
-    "    Maximal time in ms of the simulation to take into account\n",
-    "    for the calculation. Defaults to the simulation time.\n",
-    "areas : list, optional\n",
-    "    Which areas to include in the calculcation.\n",
-    "    Defaults to all loaded areas.\n",
-    "pops : list or {'complete'}, optional\n",
-    "    Which populations to include in the calculation.\n",
-    "    If set to 'complete', all populations the respective areas\n",
-    "    are included. Defaults to 'complete'.\n",
-    "\"\"\"\n",
-    "A.create_pop_cv_isi()"
+    "### 5.6 Irregularity of spiking activity <a class=\"anchor\" id=\"section_5_6\"></a>\n",
+    "Fig 5. (F) Irregularity measured by revised local variation LvR averaged across neurons"
    ]
   },
   {
@@ -1120,55 +1045,11 @@
    "cell_type": "markdown",
    "id": "90ae8f4c",
    "metadata": {
-    "jp-MarkdownHeadingCollapsed": true,
     "tags": []
    },
    "source": [
-    "### 4.6 Time series of population- and area-averaged firing rates <a class=\"anchor\" id=\"section_4_6\"></a>\n",
-    "Area-averaged firing rates, shown as raw binned spike histograms with 1ms bin width (gray) and convolved histograms, with aGaussian kernel (black) of optimal width"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "94b0b0c4-d70b-4c49-8b5d-e1ca75f0ccf4",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "\"\"\"\n",
-    "Calculate time series of population- and area-averaged firing rates.\n",
-    "Uses ah.pop_rate_time_series.\n",
-    "If the rates have previously been stored with the\n",
-    "same parameters, they are loaded from file.\n",
-    "\n",
-    "\n",
-    "Parameters\n",
-    "----------\n",
-    "t_min : float, optional\n",
-    "    Minimal time in ms of the simulation to take into account\n",
-    "    for the calculation. Defaults to 500 ms.\n",
-    "t_max : float, optional\n",
-    "    Maximal time in ms of the simulation to take into account\n",
-    "    for the calculation. Defaults to the simulation time.\n",
-    "areas : list, optional\n",
-    "    Which areas to include in the calculcation.\n",
-    "    Defaults to all loaded areas.\n",
-    "pops : list or {'complete'}, optional\n",
-    "    Which populations to include in the calculation.\n",
-    "    If set to 'complete', all populations the respective areas\n",
-    "    are included. Defaults to 'complete'.\n",
-    "kernel : {'gauss_time_window', 'alpha_time_window', 'rect_time_window'}, optional\n",
-    "    Specifies the kernel to be convolved with the spike histogram.\n",
-    "    Defaults to 'binned', which corresponds to no convolution.\n",
-    "resolution: float, optional\n",
-    "    Width of the convolution kernel. Specifically it correponds to:\n",
-    "    - 'binned' : bin width of the histogram\n",
-    "    - 'gauss_time_window' : sigma\n",
-    "    - 'alpha_time_window' : time constant of the alpha function\n",
-    "    - 'rect_time_window' : width of the moving rectangular function\n",
-    "\"\"\"\n",
-    "A.create_rate_time_series(t_max=1000.)\n",
-    "# M.analysis.save()"
+    "### 5.7 Time series of area-averaged firing rates <a class=\"anchor\" id=\"section_5_7\"></a>\n",
+    "Area-averagedfiring rates, shown as raw binned spike histograms with 1ms bin width (gray) and convolvedhistograms, with aGaussia nkernel (black) of optimal width"
    ]
   },
   {
diff --git a/multi-area-model.ipynb b/multi-area-model.ipynb
index 75bd71314fc0ca988c786a9b6b1b3b6b91ac9243..e198e1b40beb8f45adf43cca81399c92105a1689 100644
--- a/multi-area-model.ipynb
+++ b/multi-area-model.ipynb
@@ -564,6 +564,7 @@
    "cell_type": "markdown",
    "id": "57ff902c-d6ce-4f96-9e4f-8e3e7166ab66",
    "metadata": {
+    "jp-MarkdownHeadingCollapsed": true,
     "tags": []
    },
    "source": [
@@ -783,7 +784,10 @@
   {
    "cell_type": "markdown",
    "id": "bb71c922",
-   "metadata": {},
+   "metadata": {
+    "jp-MarkdownHeadingCollapsed": true,
+    "tags": []
+   },
    "source": [
     "## S5. Simulation results visualziation <a class=\"anchor\" id=\"section_5\"></a>"
    ]