diff --git a/multi-area-model.ipynb b/multi-area-model.ipynb
index 8edb0c0b628a855ff3a4687c0a1572cd09495fd2..631cbcf2ccf93bc5e0e54bbd5627142d01e6455c 100644
--- a/multi-area-model.ipynb
+++ b/multi-area-model.ipynb
@@ -1,5 +1,26 @@
 {
  "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "b1331599",
+   "metadata": {},
+   "source": [
+    "# Down-scaled multi-area model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "96517739",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import os"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -10,16 +31,38 @@
     "!pip install nested_dict dicthash"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "0c6b6a3c",
+   "metadata": {},
+   "source": [
+    "Create config file."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "61e83d47",
+   "id": "72c170e4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "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,
+   "id": "2784f76b",
    "metadata": {},
    "outputs": [],
    "source": [
-    "import numpy as np\n",
-    "import os\n",
-    "\n",
     "from multiarea_model import MultiAreaModel"
    ]
   },
@@ -28,10 +71,11 @@
    "id": "2cedd26b",
    "metadata": {},
    "source": [
-    "# Down-scaled model\n",
-    "\n",
     "Neurons and indegrees are both scaled down to 1%.\n",
-    "Can usually be simulated on a local machine."
+    "Can usually be simulated on a local machine.\n",
+    "\n",
+    "**Warning: This will not yield reasonable dynamical results from the\n",
+    "network and is only meant to demonstrate the simulation workflow.**"
    ]
   },
   {
@@ -46,11 +90,10 @@
   },
   {
    "cell_type": "markdown",
-   "id": "241db524",
+   "id": "d53f1eab",
    "metadata": {},
    "source": [
-    "**Warning: This will not yield reasonable dynamical results from the\n",
-    "network and is only meant to demonstrate the simulation workflow.**"
+    "Specify model and simulation parameters."
    ]
   },
   {
@@ -60,7 +103,6 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "d = {}\n",
     "conn_params = {'replace_non_simulated_areas': 'het_poisson_stat',\n",
     "               'g': -11.,\n",
     "               \n",
@@ -88,110 +130,109 @@
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "d409be95",
+   "cell_type": "markdown",
+   "id": "de4a6703",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "M = MultiAreaModel(network_params, simulation=True,\n",
-    "                   sim_spec=sim_params,\n",
-    "                   theory=True,\n",
-    "                   theory_spec=theory_params)\n",
-    "p, r = M.theory.integrate_siegert()\n",
-    "print(\"Mean-field theory predicts an average \"\n",
-    "      \"rate of {0:.3f} spikes/s across all populations.\".format(np.mean(r[:, -1])))"
+    "Instantiate a multi-area model, predict firing rates from theroy, and run the simulation. "
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "15778e9c",
+   "id": "d409be95",
    "metadata": {},
    "outputs": [],
    "source": [
-    "M.simulation.simulate()"
+    "M = MultiAreaModel(network_params, simulation=True,\n",
+    "                   sim_spec=sim_params,\n",
+    "                   theory=True,\n",
+    "                   theory_spec=theory_params)"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "cb8e3edd",
+   "id": "918d907f",
    "metadata": {},
    "outputs": [],
    "source": [
-    "data = np.loadtxt(M.simulation.label + \"-spikes-1-0.dat\", skiprows=3)"
+    "p, r = M.theory.integrate_siegert()\n",
+    "print(\"Mean-field theory predicts an average \"\n",
+    "      \"rate of {0:.3f} spikes/s across all populations.\".format(np.mean(r[:, -1])))"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "84dc3bde",
+   "id": "15778e9c",
    "metadata": {},
    "outputs": [],
    "source": [
-    "data"
+    "M.simulation.simulate()"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "9590223b",
+   "cell_type": "markdown",
+   "id": "8726a93d",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "tsteps, spikecount = np.unique(data[:,1], return_counts=True)"
+    "Load spike data."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "0c3fc302",
+   "id": "cb8e3edd",
    "metadata": {},
    "outputs": [],
    "source": [
-    "import matplotlib.pyplot as plt"
+    "data = np.loadtxt(M.simulation.data_dir + '/recordings/' + M.simulation.label + \"-spikes-1-0.dat\", skiprows=3)"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "9f21dd40",
+   "cell_type": "markdown",
+   "id": "8793e033",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "fig, ax = plt.subplots()\n",
-    "ax.plot(tsteps, spikecount)"
+    "Compute instantaneous rate per neuron across all populations."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "41028d2c",
+   "id": "9590223b",
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.show()"
+    "tsteps, spikecount = np.unique(data[:,1], return_counts=True)\n",
+    "rate = spikecount / M.simulation.params['dt'] * 1e3 / np.sum(M.N_vec)"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "c6091457",
+   "cell_type": "markdown",
+   "id": "38ddd973",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "len(spikecount), np.average(spikecount)"
+    "Plot instantaneous and mean rate."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "ad2bd7b2",
+   "id": "bea30fc8",
    "metadata": {},
    "outputs": [],
    "source": [
-    "np.average(spikecount[2000:])"
+    "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 rate across all populations')\n",
+    "ax.set_xlabel('time (ms)')\n",
+    "ax.set_ylabel('rate (spikes / s)')\n",
+    "ax.set_xlim(0, sim_params['t_sim'])\n",
+    "ax.set_ylim(0, 50)\n",
+    "ax.legend()"
    ]
   }
  ],