diff --git a/NEST-Elephant-Tutorial.ipynb b/NEST-Elephant-Tutorial.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..34d51f050eac8bd31bd90bf3b0aa4e91991b9d4a
--- /dev/null
+++ b/NEST-Elephant-Tutorial.ipynb
@@ -0,0 +1,1072 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "7820e086",
+   "metadata": {},
+   "source": [
+    "# Showcasing the interface between NEST and Elephant"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "88a9985a-a535-4762-aff3-d6b3efcb1637",
+   "metadata": {
+    "jp-MarkdownHeadingCollapsed": true,
+    "tags": []
+   },
+   "source": [
+    "See <path/to/ebrains-lab-setup/> to set up a notebook in lab.ebrains.eu\n",
+    "\n",
+    "## NEST and Elephant - overview\n",
+    "\n",
+    "NEST is a simulation software that can output a variety of data. One aim of NEST is to create biologically realistic electrophysiological output in a computer simulation, rather than a real world experiment.  Elephant is an analysis toolkit for electrophysiological data. So it can take real world data or data generated by simulators like NEST.\n",
+    "\n",
+    "Why use NEST? Why use Elephant? // should we elaborate on the benefits of the two? is there something from a neuroscientists perspective that is really attractive in using these tools?\n",
+    "\n",
+    "\n",
+    "\n",
+    "The output generated by NEST can be imported into Elephant with a few lines of code.\n",
+    "\n",
+    "\n",
+    "Here we will show an overview of what lines you need to include in both NEST and Elephant so they work together, and additionally show an entire NEST and Elephant script to show a realistic analysis using Elephant of output from NEST.\n",
+    "\n",
+    "* Rather than starting from scratch, both NEST and Elephant have examples you can work with. Explore the documentation for NEST and Elephant to learn about all the possibilities.\n",
+    "    \n",
+    "    * Available notebooks for PyNEST examples can be downloaded from online documentation (https://nest-simulator.readthedocs.io/en/v3.2/examples/index.html) \n",
+    "\n",
+    "    * Elephant provides notebooks to their tutorial series in the online documentation: https://elephant.readthedocs.io/en/latest/tutorials.html\n",
+    "    \n",
+    "    * For EBRAINS users: For ease of access, you can clone those repositories and have them available in the root of your directory on lab.ebrains.eu. \n",
+    "   \n",
+    "\n",
+    "## Quick overview of your NEST script and Elephant script \n",
+    "\n",
+    "Here is a brief summary of the lines of code required to generate the output in NEST and import into Elephant.\n",
+    "\n",
+    "### Setup in NEST (create output that can be read by Elephant)\n",
+    "\n",
+    "In your NEST script, you need to\n",
+    "\n",
+    "* create the type of recorder you want output from.\n",
+    "    \n",
+    "        espikes = nest.Create(\"spike_recorder\")        \n",
+    "        ispikes = nest.Create(\"spike_recorder\")\n",
+    "        \n",
+    "* Set the recorder to output to text file (ascii) in .dat or .gdf format. See documentation for more info: https://nest-simulator.readthedocs.io/en/v3.2/models/recording_backend_ascii.html\n",
+    "    \n",
+    "        espikes.set(label=\"brunel-py-ex\", record_to=\"ascii\", file_extension=\"gdf\")        \n",
+    "        ispikes.set(label=\"brunel-py-in\", record_to=\"ascii\", file_extension=\"gdf\")\n",
+    "        \n",
+    "* Connect your recorder in your network and run the simulation. Your output will appear after simulation.\n",
+    "     \n",
+    "        nest.Connect(nodes_ex, espikes, syn_spec=\"excitatory\")        \n",
+    "        nest.Connect(nodes_in, ispikes, syn_spec=\"excitatory\")        \n",
+    "        nest.Simulate(1000.0)\n",
+    "\n",
+    "### Setup in Elephant (import data from NEST)\n",
+    "\n",
+    "In your Elephant script, you need to\n",
+    "\n",
+    "* Remove the header lines from the output files. // this will change when neo is updated (since NEST 3, this has to be done)\n",
+    "     \n",
+    "        for file in gdf_files:\n",
+    "           !tail -n +4 {file}>{'no_header_'+file} # start at line 4 and save\n",
+    "           \n",
+    "* Import from neo, NESTIO\n",
+    "\n",
+    "        from neo.io.nestio import NestIO\n",
+    "        \n",
+    "* Prepare the file(s) so it's usable by Elephant\n",
+    "\n",
+    "        file = ('no_header_'+gdf_files[0])\n",
+    "        IO = NestIO(filenames=file)\n",
+    "        seg = IO.read_segment(\n",
+    "                     gid_list=[], # empty list to retrieve spiketrains from all neurons\n",
+    "                     t_start=0 * pq.ms,\n",
+    "                     t_stop=1000 * pq.ms,\n",
+    "                     id_column_gdf=0,  # column 0 contains the neuron ID\n",
+    "                     time_column_gdf=1, # column 1 contains the times in ms\n",
+    "                     )\n",
+    "\n",
+    "\n",
+    "Now let's take a look at a complete NEST and Elephant script to fully see what this looks like."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "cb2aef7e-90d0-4933-a488-6759e30a2f2d",
+   "metadata": {
+    "tags": []
+   },
+   "source": [
+    "# Full example - NEST script that outputs spike data and is analyzed using Elephant  // delete output files brunel-py* before running!)\n",
+    "\n",
+    "\n",
+    "* Simulating a Brunel model (random balanced network) with delta synapses, and outputting spike data to file.\n",
+    "\n",
+    "\n",
+    "* Calculate instantaneous firing rates and investigate spike train correlations"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c3989f0e-40d0-4e37-b184-f12240c8fc5b",
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "%matplotlib inline"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c291a98d-3adb-4b47-9157-4418219e53bf",
+   "metadata": {},
+   "source": [
+    "\n",
+    "## Random balanced network (delta synapses)\n",
+    "\n",
+    "This script simulates an excitatory and an inhibitory population on\n",
+    "the basis of the network used in [1]_\n",
+    "\n",
+    "This example of a random balanced network has been adapted from https://nest-simulator.readthedocs.io/en/v3.2/auto_examples/brunel_delta_nest.html.\n",
+    "\n",
+    "When connecting the network, customary synapse models are used, which\n",
+    "allow for querying the number of created synapses. Using spike\n",
+    "recorders, the average firing rates of the neurons in the populations\n",
+    "are established. The building as well as the simulation time of the\n",
+    "network are recorded.\n",
+    "\n",
+    "### References\n",
+    "\n",
+    ".. [1] Brunel N (2000). Dynamics of sparsely connected networks of excitatory and\n",
+    "       inhibitory spiking neurons. Journal of Computational Neuroscience 8,\n",
+    "       183-208.\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e1ff4fee-33df-42c0-a524-ad4c616dbb46",
+   "metadata": {},
+   "source": [
+    "Import all necessary modules "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fcafe876-b966-45c6-9244-27bb06362e6d",
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "import nest\n",
+    "\n",
+    "nest.ResetKernel()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e880eaba-e4fc-4c64-9c1f-5c4d16ade155",
+   "metadata": {},
+   "source": [
+    "Assigning the simulation parameters to variables.\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ad72b701-310f-49b9-8f1b-67b247722e9d",
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "dt = 0.1  # the resolution in ms\n",
+    "simtime = 1000.0  # Simulation time in ms\n",
+    "delay = 1.5  # synaptic delay in ms"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f2b8afb8-7c40-434b-86ec-5513fdb43fd1",
+   "metadata": {},
+   "source": [
+    "Definition of the parameters crucial for asynchronous irregular firing of\n",
+    "the neurons.\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0725fd8d-dc2a-4eb0-80d1-d18e22309cbd",
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "g = 5.0  # ratio inhibitory weight/excitatory weight\n",
+    "eta = 2.0  # external rate relative to threshold rate\n",
+    "epsilon = 0.1  # connection probability"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "28c5bbf5-ef74-4036-bcf4-1e673f39a39c",
+   "metadata": {},
+   "source": [
+    "Definition of the number of neurons in the network and the number of neurons\n",
+    "recorded from\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8c428035-04b4-4f1d-b298-bd35c62c0635",
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "order = 500 # originally 2500\n",
+    "NE = 4 * order  # number of excitatory neurons\n",
+    "NI = 1 * order  # number of inhibitory neurons\n",
+    "N_neurons = NE + NI  # number of neurons in total\n",
+    "N_rec = 50  # record from 50 neurons"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0e05bfe1-1b80-4472-9a0a-25df0658a1a1",
+   "metadata": {},
+   "source": [
+    "Definition of connectivity parameters\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "96bc58df-8e91-4726-aa33-a62d28b8aed3",
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "CE = int(epsilon * NE)  # number of excitatory synapses per neuron\n",
+    "CI = int(epsilon * NI)  # number of inhibitory synapses per neuron\n",
+    "C_tot = int(CI + CE)  # total number of synapses per neuron"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d284755e-9b79-49c2-9514-86c639ea79ed",
+   "metadata": {},
+   "source": [
+    "Initialization of the parameters of the integrate and fire neuron and the\n",
+    "synapses. The parameters of the neuron are stored in a dictionary.\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ee8fe368-9724-40ec-a878-37c9e1e102d3",
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "tauMem = 20.0  # time constant of membrane potential in ms\n",
+    "theta = 20.0  # membrane threshold potential in mV\n",
+    "neuron_params = {\"C_m\": 1.0,\n",
+    "                 \"tau_m\": tauMem,\n",
+    "                 \"t_ref\": 2.0,\n",
+    "                 \"E_L\": 0.0,\n",
+    "                 \"V_reset\": 0.0,\n",
+    "                 \"V_m\": 0.0,\n",
+    "                 \"V_th\": theta}\n",
+    "J = 0.1  # postsynaptic amplitude in mV\n",
+    "J_ex = J  # amplitude of excitatory postsynaptic potential\n",
+    "J_in = -g * J_ex  # amplitude of inhibitory postsynaptic potential"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5d6ec304-6ed4-4e13-8015-80273ba9354a",
+   "metadata": {},
+   "source": [
+    "Definition of threshold rate, which is the external rate needed to fix the\n",
+    "membrane potential around its threshold, the external firing rate and the\n",
+    "rate of the poisson generator which is multiplied by the in-degree CE and\n",
+    "converted to Hz by multiplication by 1000.\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "711121df-9762-498f-9c1e-d746b9da6ec7",
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "nu_th = theta / (J * CE * tauMem)\n",
+    "nu_ex = eta * nu_th\n",
+    "p_rate = 1000.0 * nu_ex * CE"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "85c6a71a-a393-4b35-adba-81eb4937adb4",
+   "metadata": {},
+   "source": [
+    "Configuration of the simulation kernel by the previously defined time\n",
+    "resolution used in the simulation. \n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "371f239f-cf84-4926-bfca-a5b65d5feb22",
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "nest.resolution = dt\n",
+    "nest.overwrite_files = True\n",
+    "\n",
+    "print(\"Building network\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "fbf20f51-8073-48fc-bd7f-19178fd263b9",
+   "metadata": {},
+   "source": [
+    "Creation of the nodes using ``Create``. We store the returned handles in\n",
+    "variables for later reference. Here the excitatory and inhibitory, as well\n",
+    "as the poisson generator and two spike recorders. The spike recorders will\n",
+    "later be used to record excitatory and inhibitory spikes. Properties of the\n",
+    "nodes are specified via ``params``, which expects a dictionary.\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f0320f7b-0695-4198-b139-0ede14326d4b",
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "nodes_ex = nest.Create(\"iaf_psc_delta\", NE, params=neuron_params)\n",
+    "nodes_in = nest.Create(\"iaf_psc_delta\", NI, params=neuron_params)\n",
+    "noise = nest.Create(\"poisson_generator\", params={\"rate\": p_rate})\n",
+    "espikes = nest.Create(\"spike_recorder\")\n",
+    "ispikes = nest.Create(\"spike_recorder\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1b759958-32f0-4e80-a33a-6a54117746a1",
+   "metadata": {},
+   "source": [
+    "Configuration of the spike recorders recording excitatory and inhibitory\n",
+    "spikes by sending parameter dictionaries to ``set``. Setting the property\n",
+    "`record_to` to *\"ascii\"* ensures that the spikes will be recorded to a file,\n",
+    "whose name starts with the string assigned to the property `label`.\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6fe48bbf-1301-44bc-9df5-86b4ef98cf39",
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "espikes.set(label=\"brunel-py-ex\", record_to=\"ascii\", file_extension=\"gdf\")\n",
+    "ispikes.set(label=\"brunel-py-in\", record_to=\"ascii\", file_extension=\"gdf\")\n",
+    "\n",
+    "print(\"Connecting devices\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f42c16da-3236-4b19-bf41-d7f0b1bc6347",
+   "metadata": {},
+   "source": [
+    "Definition of a synapse using ``CopyModel``, which expects the model name of\n",
+    "a pre-defined synapse, the name of the customary synapse and an optional\n",
+    "parameter dictionary. The parameters defined in the dictionary will be the\n",
+    "default parameter for the customary synapse. Here we define one synapse for\n",
+    "the excitatory and one for the inhibitory connections giving the\n",
+    "previously defined weights and equal delays.\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d3bc0a07-e8d9-4bed-95ba-a4e27e8da154",
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "nest.CopyModel(\"static_synapse\", \"excitatory\",\n",
+    "               {\"weight\": J_ex, \"delay\": delay})\n",
+    "nest.CopyModel(\"static_synapse\", \"inhibitory\",\n",
+    "               {\"weight\": J_in, \"delay\": delay})"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0d7284b8-eea5-4a52-b171-cd91b2fd25d5",
+   "metadata": {},
+   "source": [
+    "Connecting the previously defined poisson generator to the excitatory and\n",
+    "inhibitory neurons using the excitatory synapse. Since the poisson\n",
+    "generator is connected to all neurons in the population the default rule\n",
+    "(# ``all_to_all``) of ``Connect`` is used. The synaptic properties are inserted\n",
+    "via ``syn_spec`` which expects a dictionary when defining multiple variables\n",
+    "or a string when simply using a pre-defined synapse.\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "efbbdd70-944a-4ba0-87a4-2c8b189f3604",
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "nest.Connect(noise, nodes_ex, syn_spec=\"excitatory\")\n",
+    "nest.Connect(noise, nodes_in, syn_spec=\"excitatory\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e05da954-362b-496e-8d0e-fdb5f1fdb6e7",
+   "metadata": {},
+   "source": [
+    "Connecting the first ``N_rec`` nodes of the excitatory and inhibitory\n",
+    "population to the associated spike recorders using excitatory synapses.\n",
+    "Here the same shortcut for the specification of the synapse as defined\n",
+    "above is used.\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c31c394b-2fc0-406b-a9a2-60be5707608e",
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "nest.Connect(nodes_ex[:N_rec], espikes, syn_spec=\"excitatory\")\n",
+    "nest.Connect(nodes_in[:N_rec], ispikes, syn_spec=\"excitatory\")\n",
+    "\n",
+    "print(\"Connecting network\")\n",
+    "\n",
+    "print(\"Excitatory connections\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "291bca8b-9145-43ae-8aa9-fde64d02a29b",
+   "metadata": {},
+   "source": [
+    "Connecting the excitatory population to all neurons using the pre-defined\n",
+    "excitatory synapse. Beforehand, the connection parameter are defined in a\n",
+    "dictionary. Here we use the connection rule ``fixed_indegree``,\n",
+    "which requires the definition of the indegree. Since the synapse\n",
+    "specification is reduced to assigning the pre-defined excitatory synapse it\n",
+    "suffices to insert a string.\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b74d5425-1723-4b16-97ce-ce3fb6287888",
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "conn_params_ex = {'rule': 'fixed_indegree', 'indegree': CE}\n",
+    "nest.Connect(nodes_ex, nodes_ex + nodes_in, conn_params_ex, \"excitatory\")\n",
+    "\n",
+    "print(\"Inhibitory connections\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ff1deac7-31fb-4541-88c4-0e7139e963b4",
+   "metadata": {},
+   "source": [
+    "Connecting the inhibitory population to all neurons using the pre-defined\n",
+    "inhibitory synapse. The connection parameters as well as the synapse\n",
+    "parameters are defined analogously to the connection from the excitatory\n",
+    "population defined above.\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bedcb3c7-13fb-4524-a572-014dfff95f56",
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "conn_params_in = {'rule': 'fixed_indegree', 'indegree': CI}\n",
+    "nest.Connect(nodes_in, nodes_ex + nodes_in, conn_params_in, \"inhibitory\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "2a0308e3-bb90-48e7-8857-860f9938e1fe",
+   "metadata": {},
+   "source": [
+    "Simulation of the network.\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e666a00a-8dcd-434b-a05d-9d651a040f00",
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "print(\"Simulating\")\n",
+    "\n",
+    "nest.Simulate(simtime)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "20424db7-4ecd-48a3-8e44-9811b0dfd156",
+   "metadata": {},
+   "source": [
+    "Storage of the time point after the simulation of the network in a variable.\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e97e2135-e995-4d16-a777-89c31f2107dc",
+   "metadata": {
+    "tags": []
+   },
+   "source": [
+    "## Processing NEST data with Elephant\n",
+    "This script loads the simulation data generated with NEST in `brunel_alpha_nest.ipynb` and performs anaysis with Elephant (Electrophysiology Analysis Toolkit)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "181f1d0b-dbc1-450f-aea7-1beb555e9bf9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "!pip install viziphant"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e9c93f69-4c03-4235-8b93-657799488646",
+   "metadata": {},
+   "source": [
+    "---\n",
+    "### Hotfix: delete descriptive meta-data header prior to data\n",
+    "\n",
+    "**with NEST 3.0 the following changes occured:** \\\n",
+    "*In addition, most backends now write the name of the recorded variable for each column as a descriptive meta-data header prior to writing any data.\\\n",
+    "see: (https://nest-simulator.readthedocs.io/en/v3.1/release_notes/v3.0/features/recording_simulations.html?highlight=meta-data)*\n",
+    "\n",
+    "In order to use the `neo.io.nestio` the meta-data header for each column needs to be removed.\n",
+    "\n",
+    "NestIO was developed with NEST 2.10.0, see `nestio.py`:\\\n",
+    "*Class for reading output files from NEST simulations\n",
+    "( http://www.nest-simulator.org/ ).* \\\n",
+    "**Tested with NEST 2.10.0** \\\n",
+    "*Depends on: numpy, quantities* \\\n",
+    "*Supported: Read* \\\n",
+    "*Authors: Julia Sprenger, Maximilian Schmidt, Johanna Senk*\n",
+    "\n",
+    "Until the IO for NEST in Neo is updated, the following hotfix is applied to the `.gdf`-files:\\\n",
+    "(*once Neo is updated, simply remove* `'no_header_'`)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ed65043b-f734-443f-b9b5-610bdf3e8ec0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gdf_files=['brunel-py-ex-2502-0.gdf','brunel-py-in-5203-0.gdf']\n",
+    "#gdf_folder=('brunel_alpha_nest_output'+'/')\n",
+    "\n",
+    "for file in gdf_files:\n",
+    "    !tail -n +4 {file}>{'no_header_'+file} # start at line 4 and save"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c40a36f2-bcba-4397-ae4d-af586417eee1",
+   "metadata": {},
+   "source": [
+    "---"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ffe00b5a-2486-4438-9c69-9702e55857ed",
+   "metadata": {},
+   "source": [
+    "## Load simulation data \n",
+    "\n",
+    "In this example the data from the `.gdf` files is loaded using the nestIO from Neo package.\n",
+    "\n",
+    "see also: https://neo.readthedocs.io/"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6928e585-a812-40e6-ad7e-1873d9ddf683",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import quantities as pq\n",
+    "from neo.io.nestio import NestIO\n",
+    "\n",
+    "file = ('no_header_'+gdf_files[0])\n",
+    "IO = NestIO(filenames=file)\n",
+    "seg = IO.read_segment(\n",
+    "                     gid_list=[], # empty list to retrieve spiketrains from all neurons\n",
+    "                     t_start=0 * pq.ms,\n",
+    "                     t_stop=1000 * pq.ms,\n",
+    "                     id_column_gdf=0,  # column 0 contains the neuron ID\n",
+    "                     time_column_gdf=1, # column 1 contains the times in ms\n",
+    "                     )"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "118dde25-b640-4f0e-948e-81c75f27b14e",
+   "metadata": {},
+   "source": [
+    "## 1. Elephant: Calculate instantaneous firing rates\n",
+    "\n",
+    "see also https://elephant.readthedocs.io/\n",
+    "\n",
+    "The mean firing rate is just the temporal average of the number of spikes per time in the spiking activity of the neuron for the\n",
+    "whole trial. We can have a time-varying estimate of the firing rate along the duration of the trial by computing the instantaneous rate.\n",
+    "\n",
+    "For this estimation, Elephant's `statistics` module provides the `instantaneous_rate` function. The function performs kernel\n",
+    "convolution: each spike that ocurred during the trial is blurred with a kernel function. The bandwidth of the kernel function is\n",
+    "controlled by the parameter `sigma`. Different types of kernel and bandwidths can be selected for more or less smooth estimates.\n",
+    "\n",
+    "Let's plot the instantaneous rate for one neuron that we are analyzing. Let's use a Gaussian kernel with two different `sigma`\n",
+    "parameters: `20 ms` or `100 ms`. Use `.1 ms` as `sampling_period`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "11792f80-bbff-4c0f-b213-f8b439552dbc",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from elephant import statistics\n",
+    "from elephant import kernels\n",
+    "\n",
+    "spiketrain_id=0\n",
+    "# Define kernel 1\n",
+    "kernel_1 = kernels.GaussianKernel(sigma=20 * pq.ms)\n",
+    "# calculate instantaneous firing rate 1\n",
+    "inst_rate_1 = statistics.instantaneous_rate(\n",
+    "                                    seg.spiketrains[spiketrain_id],\n",
+    "                                    sampling_period=.1 * pq.ms,\n",
+    "                                    kernel=kernel_1\n",
+    "                                    )\n",
+    "\n",
+    "# Define kernel 2\n",
+    "kernel_2 = kernels.GaussianKernel(sigma=100 * pq.ms)\n",
+    "# calculate instantaneous firing rate 2\n",
+    "inst_rate_2 = statistics.instantaneous_rate(\n",
+    "                                    seg.spiketrains[spiketrain_id],\n",
+    "                                    sampling_period=.1 * pq.ms,\n",
+    "                                    kernel=kernel_2\n",
+    "                                    )"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b2f0fd19-f830-439c-9c73-890bb0395ac4",
+   "metadata": {},
+   "source": [
+    "### Plot analysis results\n",
+    "\n",
+    "Execute the code below to plot the spike train together with the rate profile that we calculated.\n",
+    "The spike times will be added to the bottom of each plot. Note that the parameters of the rate computations\n",
+    "are retained in the resulting rate variables, which in turn are Neo `AnalogSignal` objects."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f819f76f-ce68-45b9-9bd0-0ef1ad44861f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "fig, axes = plt.subplots(2, 1, sharey=True, figsize=(12,6))\n",
+    "for plot, rates in enumerate([inst_rate_1, inst_rate_2]):\n",
+    "    axes[plot].plot(rates.times.rescale(pq.ms), rates, linewidth=2)\n",
+    "    axes[plot].eventplot(seg.spiketrains[spiketrain_id].magnitude, color='black', linelengths=5)\n",
+    "    axes[plot].set_ylabel(\"Firing rate (Hz)\")\n",
+    "    axes[plot].set_title(f\"$\\sigma$ = {rates.annotations['kernel']['sigma']}\")\n",
+    "axes[-1].set_xlabel(\"Time (ms)\");"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7b0b28fd-62b9-4731-906d-f9fae841a09d",
+   "metadata": {},
+   "source": [
+    "We can see that the smoothness of the rate curve is affected by the sigma parameter. Elephant can try to estimate the best kernel bandwidth for a Gaussian kernel by using 'auto' as the value of the kernel parameter. Let's try this option."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "49ea8b4f-3d88-424a-8edd-acb2e92c536e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "rates_auto = statistics.instantaneous_rate(\n",
+    "                                            seg.spiketrains[spiketrain_id],\n",
+    "                                            kernel='auto', \n",
+    "                                            sampling_period=.1 * pq.ms,\n",
+    "                                            )"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "469f5bd4-e9d7-48fd-815c-f1b37c477738",
+   "metadata": {},
+   "source": [
+    "Execute the code below to plot the instantaneous firing rate obtained from the automatic kernel bandwidth."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d5c13889-2c4e-4987-85ef-a68af212db58",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, axes = plt.subplots(figsize=(12,3))\n",
+    "axes.plot(rates_auto.times.rescale(pq.ms), rates_auto)\n",
+    "axes.eventplot(seg.spiketrains[spiketrain_id].magnitude, color='black', linelengths=5)\n",
+    "axes.set_ylabel(\"Firing rate [Hz]\")\n",
+    "axes.set_xlabel(\"Time (ms)\")\n",
+    "axes.set_title(f\"$\\sigma$ = {rates_auto.annotations['kernel']['sigma']}\");"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "604b8e36-7da0-4166-bb4e-996e8d597a38",
+   "metadata": {},
+   "source": [
+    "Finally, let's compute and plot the instantaneous firing rates for all spiketrains.\n",
+    "\n",
+    "You can pass the list with all spike trains for the instantaneous_rate function, to generate a 2D-matrix. Use `kernel_1` as the kernel parameter and .1 ms as sampling_period."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8a86e14d-d799-418e-bc47-566e0d7fb25e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "all_rates = statistics.instantaneous_rate(\n",
+    "                                            seg.spiketrains, \n",
+    "                                            kernel=kernel_1, \n",
+    "                                            sampling_period=.1 * pq.ms,\n",
+    "                                            )"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "931e7e5b-be74-45b6-9205-7d0980f84b84",
+   "metadata": {},
+   "source": [
+    "### Viziphant\n",
+    "\n",
+    "Now use the Viziphant function `plot_instantaneous_rates_colormesh` to plot the rates of all spiketrains.\n",
+    "\n",
+    "see also https://viziphant.readthedocs.io/"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e6773f36-7510-479b-a71f-8f77a5bbe3ce",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import viziphant\n",
+    "# Create the axes, setting the figure size\n",
+    "fig, axes = plt.subplots(figsize=(15,6))\n",
+    "# Plot the rates using the axes\n",
+    "viziphant.statistics.plot_instantaneous_rates_colormesh(all_rates, axes=axes);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "21a71b93-d0d5-4205-af6a-21fd08c51b29",
+   "metadata": {},
+   "source": [
+    "## 2. Elephant: Investigate the spike train correlations\n",
+    "\n",
+    "From the rasterplot that we generated above, it is difficult to assess if there are correlations between the spike trains.\n",
+    "\n",
+    "To investigate that, Elephant provides functions to calculate and plot the cross-correlation matrix. This matrix quantifies the similarity for each pair of spike trains in the trial. We start by binning the spike trains. With this, we are obtaining the number of spikes that occurred during small intervals. Let's use a bin size of 3 ms."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4f438dd0-70d5-44d7-b3d0-8e9aed471841",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from elephant import conversion\n",
+    "binned_spiketrains = conversion.BinnedSpikeTrain(\n",
+    "                                                seg.spiketrains, \n",
+    "                                                binsize=3*pq.ms\n",
+    "                                                )"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "caf909fb-6d19-4497-9cb5-25af252f1d42",
+   "metadata": {},
+   "source": [
+    "The binning above showed several warnings that happened due to the machine error precision of floating point computation. Since now we're aware of it, let's filter them out throughout the rest of the notebook."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e977f0e7-187d-4669-bab0-f2428c0d362d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import warnings\n",
+    "warnings.filterwarnings(\"ignore\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c8c32ce7-a246-454d-93f2-de8aa1cfcaf5",
+   "metadata": {},
+   "source": [
+    "If we inspect the spike times of one neuron, we can see that if the neuron fired in a particular bin,\n",
+    "the `BinnedSpikeTrain` object will store the value 1."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6c033582-d2d3-48c0-a777-64d12e9e690e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(\"Timepoints in ms where spike occured\")\n",
+    "print(seg.spiketrains[spiketrain_id])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "25095160-62fb-45a0-9dcd-78dca20f563d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print (\"Index of the 3 ms bin that the spike occurred\")\n",
+    "print(seg.spiketrains[spiketrain_id].times / (3 * pq.ms))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bf026c31-43a4-499f-a0b1-46a0a4c150f7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(\"Total number of bins: \"+str(binned_spiketrains.n_bins))\n",
+    "print(binned_spiketrains.to_array()[0,:])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b6121418-a6ae-4d7f-81de-0e836c8715b3",
+   "metadata": {},
+   "source": [
+    "Then we calculate the cross-correlation matrix of the binned spikes using the `correlation_coefficient`\n",
+    "function from the `spike_train_correlation` module."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "65f496c7-f535-4587-b78b-e36346f02af1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from elephant import spike_train_correlation\n",
+    "cross_corr_matrix = spike_train_correlation.correlation_coefficient(\n",
+    "    binned_spiketrains)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "da769643-a05b-4ac9-8f5f-f47d253d4e22",
+   "metadata": {},
+   "source": [
+    "### Viziphant\n",
+    "We visualize the matrix using the plot_corrcoef function of Viziphant. To better visualize the coefficients, we will set correlation_range to 'auto', to use the color bar only in the range of the obtained coefficients. We will also not plot the values along the main diagonal, as those are equal to 1, by setting remove_diagonal to True."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5d18702e-5504-48a9-a90f-0e5422307065",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Create the axes, setting the figure size\n",
+    "fig, axes = plt.subplots(figsize=(10,5))\n",
+    "# Plot the correlation matrix\n",
+    "viziphant.spike_train_correlation.plot_corrcoef(\n",
+    "    cross_corr_matrix, axes=axes,correlation_range='auto',remove_diagonal=True)\n",
+    "# Set labels and title\n",
+    "axes.set_xlabel('Neuron')\n",
+    "axes.set_ylabel('Neuron')\n",
+    "axes.set_title(\"Correlation coefficient matrix\");"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a8c53593-b49b-4fdc-a4ac-edec08a3ef69",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "EBRAINS_release_v0.1_202109",
+   "language": "python",
+   "name": "spack_python_kernel_release_20210930"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.11"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}