docs/tutorials/scripting-neuron-basics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scripting NEURON basics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The objectives of this part of the tutorial are to get familiar with basic operations of NEURON using Python. In this worksheet we will:\n",
    "\n",
    "* Create a passive cell membrane in NEURON.\n",
    "* Create a synaptic stimulus onto the neuron.\n",
    "* Modify parameters of the membrane and stimulus.\n",
    "* Visualize results with bokeh."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is NEURON?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The NEURON simulation environment is a powerful engine for performing simulations of neurons and biophysical neural networks. It permits the construction of biologically realistic membranes with active and passive ion channels, combined with virtual connectivity and electrophysiology tools to drive and measure neuron and network behaviors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1: Import the <tt>neuron</tt> module into Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Any code that is not part of Python's <a class=\"reference external\" href=\"https://docs.python.org/3/library/functions.html#built-in-funcs\" title=\"(in Python v3.6)\"><span>Built-in Functions</span></a> must be imported. The Python interface to NEURON goes through the <tt>neuron</tt> module, especially the <tt>neuron.h</tt> submodule. The <tt>neuron</tt> module has additional submodules, including <tt>neuron.rxd</tt> for reaction-diffusion dynamics, <tt>neuron.gui2</tt> for Jupyter-compatible PlotShape graphs, and <tt>neuron.gui</tt> for Interviews-based GUI tools. The <tt>neuron.h</tt> submodule also allows loading files, executing code, and calling functions written in HOC, an older scripting language supported by NEURON. This allows the continued use of HOC libraries in Python code."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To import <tt>neuron</tt>, we <i>could</i> use:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import neuron"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the above succeeded, it produces no output (in recent versions of NEURON), so how can we know what version of NEURON we have? Simple: ask for the `__version__`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(neuron.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are only a limited number of functions avaiable directly from the <tt>neuron</tt> module. In practice, we usually want to directly import the submodules we need; i.e. <b>do something like the below instead</b>:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from neuron import h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When using NEURON, you will always want the <tt>h</tt> submodule. You may or may not need to import the additional submodules mentioned above. If you do, they can be imported separately or loaded in one line with a comma separated list, as in:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from neuron import h, rxd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NB: When importing <tt>h</tt>, etc like this, there is usually no need for importing <tt>neuron</tt> separately."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>NEURON assumes certain default units (concentration in mM, time in ms, voltage in mV), but units can be specified explicitly by importing unit definitions from <tt>neuron.units</tt>. Even if you are using the default units, being explicit makes your code more readable by others. For example:</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from neuron.units import ms, mV"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the above gives you an error, then you are using a version of NEURON older than 7.7. Update before proceeding."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: Create a cell"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A <a class=\"reference external\" href=\"https://neuron.yale.edu/neuron/static/py_doc/modelspec/programmatic/topology/geometry.html#geometry\" title=\"(in NEURON v7.5)\"><span class=\"xref std std-ref\">Section</span></a> is the basic morphological building-block in NEURON. We typically think of a Section as an unbranched cable, but it can also be used to represent a soma. Thus a simple model neuron with only a soma can be created as in:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "soma = h.Section(name=\"soma\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is no output, so how can we tell that we successfully created a Section?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aside 1: NEURON's <tt>h.topology</tt> function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>NEURON's <tt>h.topology()</tt> function displays the topological structure of the entire model, indicating which sections are connected to which sections, where they are connected, and how many <i>segments</i> each section is divided into.</p>\n",
    "<p>If you're following along with our example, there's not much to see yet since there is only one section, but it does demonstrate that the soma has been created and has one segment (one dash is shown):</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h.topology()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The <tt>h.topology()</tt> function displays its data to screen and returns <tt>1.0</tt> indicating success (this function always succeeds). Note: This function is only for displaying data; other methods must be used to store the data in a variable for programmatic analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aside 2: The <tt>psection</tt> method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>Every NEURON section has a <tt>psection</tt> method (think: <b>p</b>roperties of the <b>section</b>) that returns a Python dictionary providing a structured data representation of the properties of the section.</p><p>For example, we can query the soma via:</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "soma.psection()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results tell us the soma is a cylinder with length 100 microns, diameter 500 microns, axial resistivity 35.4 ohm*cm, and specific membrance capacitance 1 μF/cm<sup>2</sup>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Note:</b> calling this method does not itself print anything to the screen. Instead it returns a dictionary. We see the contents of the dictionary only because we are running interactively; from a script, nothing would be printed unless we explicitly printed it with <tt>print</tt>, or, better pretty-printed it with <tt>pprint.pprint</tt>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since this is a dictionary, we can extract any properties we want using square brackets. For example, the length of the section is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "soma.psection()[\"morphology\"][\"L\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All of these values can be individually accessed in more efficient ways, but <tt>psection</tt> provides an overview of the full properties of the section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, the length of the <tt>soma</tt> is more efficiently available (and settable) via:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "soma.L"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<tt>soma.psection()['morphology']['diam']</tt> is a list (of length 1 here), with each entry corresponding to the value for each segment. Similarly for <tt>soma.psection()['cm']</tt>, etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Important:</b> You may have noticed that the default diameter is 500 μm, which is excessively large for mammalian neurons. It's the default because it's appropriate for the squid giant axons studied by Hodgkin and Huxley. NEURON also uses squid-relevant values for axial resistivity (<tt>soma.Ra</tt>) and temperature (<tt>h.celsius</tt>). These should all be adjusted for mammalian models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3: Set the cell's morphological properties"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we're simulating a soma, the default length of 100 μm and diameter of 500 μm are inappropriate. Let's set the length (<tt>L</tt>) and diameter (<tt>diam</tt>) to 20 μm instead:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "soma.L = 20\n",
    "soma.diam = 20"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In many models, you will have cells consisting of many connected sections. In brief, this can be done using the Section's <tt>connect</tt> method. That will be described in a subsequent part of the tutorial. For now though, we consider only the soma."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aside 3: Python's <tt>dir</tt> function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also probe objects with Python’s built-in <a class=\"reference external\" href=\"https://docs.python.org/3/library/functions.html#dir\" title=\"(in Python v3.6)\"><code class=\"xref py py-func docutils literal\"><span class=\"pre\">dir()</span></code></a> function. Let’s see what it says about soma."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dir(soma)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tells us all of the Python methods and variables associated with the object. Any methods with two leading and trailing underscores are reserved by Python. The other items in the list are additional members of soma that we can call. To see all of the functions, variables, etc available through NEURON's <tt>h</tt> submodule, try:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import textwrap\n",
    "\n",
    "print(textwrap.fill(\", \".join(dir(h))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(The <tt>', '.join(...)</tt> tells Python to build a string out of the list returned by <tt>dir</tt> where the items are separated from each other with a comma and a space. The <tt>textwrap.fill(...)</tt> tells Python to split long lines into multiple lines, by default a maximum of 70 characters long.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aside 4: Getting more help"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to probing objects with <a class=\"reference external\" href=\"https://docs.python.org/3/library/functions.html#dir\" title=\"(in Python v3.6)\"><code class=\"xref py py-func docutils literal\"><span class=\"pre\">dir()</span></code></a>, help from docstrings is available using <a class=\"reference external\" href=\"https://docs.python.org/3/library/functions.html#help\" title=\"(in Python v3.6)\"><code class=\"xref py py-func docutils literal\"><span class=\"pre\">help()</span></code></a>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, from <tt>dir(soma)</tt> above, we know that there is a <tt>connect</tt> method available. Let's inquire about that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(soma.connect)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When running interactively in Jupyter, the same information is available in a window that can be popped out by prefacing the method/function/etc with a question mark; e.g."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "?soma.connect"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Biophysical mechanisms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NEURON comes with a few built in biophysical mechanisms that can be added to a model:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<table border=\"1\" class=\"docutils\">\n",
    "<colgroup>\n",
    "<col width=\"26%\" />\n",
    "<col width=\"74%\"  />\n",
    "</colgroup>\n",
    "<tbody valign=\"top\">\n",
    "<tr class=\"row-odd\"><td><tt>pas</tt></td>\n",
    "<td>Passive (“leak”) channel.</td>\n",
    "</tr>\n",
    "<tr class=\"row-even\"><td><tt>extracellular</tt></td>\n",
    "<td>For simulating effects of nonzero extracellular\n",
    "potential, as may happen with leaky patch clamps,\n",
    "or detailed propertes of the myelin sheath.</td>\n",
    "</tr>\n",
    "<tr class=\"row-odd\"><td><tt>hh</tt></td>\n",
    "<td>Hodgkin-Huxley sodium, potassium, and leakage channels.</td>\n",
    "</tr>\n",
    "</tbody>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thousands of additional mechanisms (for A-currents, etc) are available as MOD files as part of published model codes on <a href=\"http://modeldb.yale.edu\">ModelDB</a>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4: Insert ion channels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A section's <tt>insert</tt> method is used to insert density mechanisms (i.e. anything where we don't want to specify every single instance separately). Let's insert Hodgkin-Huxley channels into the <tt>soma</tt>'s membrane. We do this by passing 'hh' as the mechanism type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "soma.insert(\"hh\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The section is returned so that multiple insertions can be chained together if desired."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We note that Hodgkin-Huxley channel kinetics are based on the squid giant axon. If that's not your model organism, then for your actual modeling projects, you'll want to use other kinetics, either by downloading them from online resources like <a href=\"http://modeldb.yale.edu\">ModelDB</a> or by writing them yourself in NMODL or NeuroML."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aside 5: Sections and segments"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A NEURON Section is considered a piece of cable. Depending on the resolution desired, it may be necessary to divide the cable into a number of segments where voltage varies linearly between centers of adjacent segments. The number of segments within a section is given by the variable, `nseg`. The total ionic current across the segment membrane is approximately the area of the segment multiplied by the ionic current density at the center of the segment. To access a part of the section, specify a value between 0 and 1, where 0 is typically the end closest to the soma and 1 is the distal end. Because nseg divides the cable into equal-length parts, it should be an odd number so that to address the middle of the cable, (0.5), gives the middle segment."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To summarize, we access sections by their name and segments by some location on the section.\n",
    "\n",
    "* Section: ``section``\n",
    "* Segment: ``section(loc)``"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the Python `type` function can tell us what a variable is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"type(soma) = {}\".format(type(soma)))\n",
    "print(\"type(soma(0.5)) = {}\".format(type(soma(0.5))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aside 6: Accessing segment variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Segment variables follow the idiom:\n",
    "\n",
    "    section(loc).var\n",
    "\n",
    "And for mechanisms on the segment:\n",
    "\n",
    "    section(loc).mech.var\n",
    "\n",
    "or\n",
    "\n",
    "    section(loc).var_mech\n",
    "\n",
    "The first form is preferred."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mech = soma(0.5).hh\n",
    "print(dir(mech))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(mech.gkbar)\n",
    "print(soma(0.5).hh.gkbar)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 5: Insert a stimulus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's insert a current clamp (an <a class=\"reference external\" href=\"https://neuron.yale.edu/neuron/static/py_doc/modelspec/programmatic/mechanisms/mech.html#IClamp\" title=\"(in NEURON v7.5)\"><code class=\"xref py py-class docutils literal\"><span class=\"pre\">IClamp</span></code></a> object) into the center of the soma to induce some membrane dynamics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iclamp = h.IClamp(soma(0.5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An <tt>IClamp</tt> is a Point Process. Point processes are point sources of current. When making a new PointProcess, you pass the segment to which it will bind.\n",
    "\n",
    "Again, with the <tt>dir</tt> function, we can validate that <tt>iclamp</tt> is an object and contains some useful parameters. Let's look at some of those parameters. We use a <a href=\"https://www.pythonforbeginners.com/basics/list-comprehensions-in-python\">list comprehension</a> to ignore those elements of the <tt>dir</tt> that start with double underscores (and are thus Python magic methods and <i>not</i> functions/variables intended to be used directly). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print([item for item in dir(iclamp) if not item.startswith(\"__\")])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In particular, we notice three key properties of a current clamp: <tt>amp</tt> -- the amplitude (in nA), <tt>delay</tt> -- the time the current clamp switches on (in ms), and <tt>dur</tt> -- how long (in ms) the current clamp stays on. Let's set these values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iclamp.delay = 2\n",
    "iclamp.dur = 0.1\n",
    "iclamp.amp = 0.9"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use <tt>psection</tt> to get a representation of the soma model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "soma.psection()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 6: Set up recording variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cell should be configured to run a simulation. However, we need to indicate which variables we wish to record; these will be stored in a NEURON Vector (<tt>h.Vector</tt> object). For now, we will record the membrane potential, which is <tt>soma(0.5).v</tt> and the corresponding time points (<tt>h.t</tt>). References to variables are available by preceding the last part of the variable name with a `_ref_`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v = h.Vector().record(soma(0.5)._ref_v)  # Membrane potential vector\n",
    "t = h.Vector().record(h._ref_t)  # Time stamp vector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 7: Run the simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default, the NEURON <tt>h</tt> module provides the low level <a href=\"https://www.neuron.yale.edu/neuron/static/py_doc/simctrl/programmatic.html#fadvance\"><tt>fadvance</tt></a> function for advancing one time step. For higher-level simulation control specification, we load NEURON's <tt>stdrun</tt> library:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h.load_file(\"stdrun.hoc\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then initialize our simulation such that our cell has a resting membrane potential of -65 mV:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h.finitialize(-65 * mV)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now continue the simulation from the current time (0) until 40 ms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h.continuerun(40 * ms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(For those who are interested: we initialized to a resting membrane potential of -65 mV because that's the default reversal potential for the <tt>hh</tt> channels, the only channel (set) inserted in this model.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Strictly speaking, we didn't need to specify the units here -- recall they were defined above in the <tt>from neuron.units import ms, mV</tt> -- as they are the defaults assumed by NEURON, but it is good practice to be explicitly clear.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 8: Plot the results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using bokeh\n",
    "<p>When working in Jupyter with an active internet connection, it is often convenient to use the <tt>bokeh</tt> module for plotting, as it provides interactive graphs that can be panned, zoomed, and saved from the Jupyter notebook.</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's load <tt>bokeh</tt> and tell it to output to the Jupyter notebook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bokeh.io import output_notebook\n",
    "import bokeh.plotting as plt\n",
    "\n",
    "output_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we plot membrane potential vs time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = plt.figure(x_axis_label=\"t (ms)\", y_axis_label=\"v (mV)\")\n",
    "f.line(t, v, line_width=2)\n",
    "plt.show(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using matplotlib"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"https://matplotlib.org/\"><tt>matplotlib</tt></a> is a mature non-Javascript based graphics library. While it does offer an interactive Jupyter mode, this interactivity blocks subsequent Python execution until interactive-mode is cancelled."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Jupyter, to ensure <tt>matplotlib</tt> graphs appear inline, use:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(When not using Jupyter, skip the above step -- it would be a syntax error in pure Python. Your graphs will appear in separate windows instead.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As with <tt>bokeh</tt>, we must import the <tt>matplotlib</tt> module before using it. In particular, we will load the <tt>pyplot</tt> submodule and give it the shorter name <tt>plt</tt>:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The <tt>matplotlib</tt> equivalent to the above <tt>bokeh</tt> example is then:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(t, v)\n",
    "plt.xlabel(\"t (ms)\")\n",
    "plt.ylabel(\"v (mV)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 9: Saving and loading results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### CSV\n",
    "<p>The <tt>csv</tt> (comma separated variables) file format is widely used for data interchange, and can be used to transfer data to MATLAB, Excel, etc without writing any special conversion code.</p><p>(Many Python distributions provide the <tt>pandas</tt> module which can do the same using a slightly simpler interface, but the code below works in all versions of Python... an example of reading CSV using pandas is shown below.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Python provides the <a href=\"https://docs.python.org/3/library/csv.html\"><tt>csv</tt></a> module to simplify reading and writing <tt>csv</tt> files. We load it via:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4>Writing</h4>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"data.csv\", \"w\") as f:\n",
    "    csv.writer(f).writerows(zip(t, v))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should now have a <tt>data.csv</tt> file. Try loading it in your favorite spreadsheet.\n",
    "\n",
    "To write additional variables (as columns) to the CSV file, just add them inside the <tt>zip</tt>; e.g. <tt>zip(t, v1, v2, ca1)</tt>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4>Reading</h4>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"data.csv\") as f:\n",
    "    reader = csv.reader(f)\n",
    "    tnew, vnew = zip(*[[float(val) for val in row] for row in reader if row])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The argument to the <tt>zip</tt> is a nested list comprehension; the zip and the asterisk together effectively transpose the data turning it from a list of (t, v) pairs into a list of t values and a list of v values. For loading more variables, the right hand side of the last line is unchanged; all that changes is that the variables need to be listed on the left; e.g. <tt>tnew, vnew, canew = zip(...)</tt>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot our newly loaded data (here with matplotlib) to see that it is the same as before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(tnew, vnew)\n",
    "plt.xlabel(\"t (ms)\")\n",
    "plt.ylabel(\"v (mV)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Reading and plotting using <tt>plotnine</tt> and <tt>pandas</tt>:</b>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The modules <tt>plotnine</tt> and <tt>pandas</tt> are key components of the Python data science ecosystem. We can use them to work with our NEURON data, but they may have to be installed separately."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we'll load the modules. It is common to use shortened names for the modules, but this is not necessary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import plotnine as p9\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can load the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv(\"data.csv\", header=None, names=[\"t\", \"v\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(If the CSV file had a header row identifying the columns, then <tt>pd.read_csv</tt> would have handled the names automatically and we would not have had to specify the last two arguments above.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now plot with <tt>plotnine</tt>'s version of <tt>ggplot</tt>. This function provides an implementation of Wilkinson's <a href=\"http://doi.org/10.1007/0-387-28695-0\">Grammar of Graphics</a>; as such, the interface is essentially identical to the R function of the same name. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = (p9.ggplot(data, p9.aes(x=\"t\", y=\"v\")) + p9.geom_path()).draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### JSON\n",
    "<p>JSON is used for structured data interchange. It is a newer but widely used format and libraries for reading and writing JSON exist for most programming languages.</p><p>Python provides the <a href=\"https://docs.python.org/3/library/json.html\">json</a> module to simplify reading and writing JSON files. We load it via:</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import json"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4>Writing</h4>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"data.json\", \"w\") as f:\n",
    "    json.dump({\"t\": list(t), \"v\": list(v)}, f, indent=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we built a dictionary with keys <tt>t</tt> and <tt>v</tt>, and stored their values as a list. Since JSON is a language-independent format, it does not have a concept of NEURON Vectors, which is why we had to create a list copy of them before saving. The <tt>indent=4</tt> argument is optional, but indents the output to make it more human-readable (at the cost of a larger file size)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4>Reading</h4>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"data.json\") as f:\n",
    "    data = json.load(f)\n",
    "tnew = data[\"t\"]\n",
    "vnew = data[\"v\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As for <tt>csv</tt>, we plot the newly loaded data to show that it is the same as the original:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(tnew, vnew)\n",
    "plt.xlabel(\"t (ms)\")\n",
    "plt.ylabel(\"v (mV)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pickles\n",
    "<p>Pickles are a Python-specific data exchange format.</p><p>Python provides the <a href=\"https://docs.python.org/3/library/pickle.html\">pickle</a> module to read and write pickled files. We load it via:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4>Writing</h4>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"data.p\", \"wb\") as f:\n",
    "    pickle.dump({\"t\": t, \"v\": v}, f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is slightly cleaner than the JSON solution above because it is Python specific and therefore able to explicitly encode NEURON Vector objects."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4>Reading</h4>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"data.p\", \"rb\") as f:\n",
    "    data = pickle.load(f)\n",
    "tnewp = data[\"t\"]\n",
    "vnewp = data[\"v\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pickles in Python 3 are by default binary files, so we have to specify write and read flags of <tt>wb</tt> and <tt>rb</tt> respectively. We use a different variable name here than before simply to indicate that what is loaded in has type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(tnewp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and in particular is a NEURON <tt>Vector</tt>:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tnewp.hname()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unlike the other solutions provided, which construct regular Python lists:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(tnew)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a minor distinction though, as we've already seen <tt>list(vec)</tt> copies a Vector <tt>vec</tt> into a new list. Using `h.Vector(old_list)` makes a NEURON Vector that is a copy of `old_list`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As before, plotting suggests that we have successfully loaded in the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(tnewp, vnewp)\n",
    "plt.xlabel(\"t (ms)\")\n",
    "plt.ylabel(\"v (mV)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}