docs/tutorials/ball-and-stick-1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ball and stick 1: Basic cell"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the first part of a several page tutorial that will ultimately build a ring network of mutlicompartment neurons running on a parallel machine.\n",
    "\n",
    "Best software engineering practice is to carefully design both the model and code. A well-designed system is easier to debug, maintain, and extend.\n",
    "\n",
    "This tutorial will take a functional bottom-up approach, building key components and refactoring as needed for a better software design instead of describing a complex design and filling in the pieces.\n",
    "\n",
    "This part of the tutorial builds a two-compartment neuron consisting of a soma and dendrite. This representation is known as a <i>ball-and-stick</i> model. After building the cell, we will attach some instrumentation to it to simulate and monitor its dynamics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load NEURON"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We of course begin by loading NEURON's main submodule <tt>h</tt>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from neuron import h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As good practice, we'll load some unit definitions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from neuron.units import ms, mV, µm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll load the standard run library to give us high-level simulation control functions (e.g. running a simulation for a given period of time):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h.load_file(\"stdrun.hoc\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we're running in Jupyter, we should allow interactive graphics inline so that we can explore our <tt>PlotShape</tt>s interactively; skip this line if you're not using Jupyter (it'll cause a syntax error):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining the cell morphology"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create the sections"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A ball-and-stick cell by definition consists of two parts: the soma (ball) and a dendrite (stick). We could define two Sections at the top level as in the previous tutorial, but that wouldn't give us an easy way to create multiple cells. Instead, let's define a <tt>BallAndStick</tt> neuron class. The basic boilerplate for defining any class in Python looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class BallAndStick:\n",
    "    def __init__(self):\n",
    "        \"\"\"anything that should be done every time one of these is created goes here\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our case in particular, every time we say that there is a <tt>BallAndStick</tt> cell, we should create the soma and the dendrite sections:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class BallAndStick:\n",
    "    def __init__(self):\n",
    "        self.soma = h.Section(name=\"soma\", cell=self)\n",
    "        self.dend = h.Section(name=\"dend\", cell=self)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Any variables that describe properties of the cell must get stored as attributes of <tt>self</tt>. This is why we write <tt>self.soma</tt> instead of <tt>soma</tt>. Temporary variables, on the other hand, need not be prefixed with self and will simply stop existing when the initialization function ends.\n",
    "\n",
    "You will also note that we have introduced a new keyword argument for <tt>Section</tt>, namely the cell attribute. This will always be set to <tt>self</tt> (i.e. the current cell) to allow each Section to know what cell it belongs to."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that we can check the topology via:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h.topology()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No sections were displayed. Why?\n",
    "\n",
    "The explanation is that we haven't actually created any such cells yet; we've only defined a rule for them. Let's go ahead and create our first cell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_cell = BallAndStick()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We again check the topology and see that the sections have in fact been created:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h.topology()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly there is a soma and a dendrite, but the cell identifier before each is not very friendly.\n",
    "\n",
    "We can specify how the cell displays using the `__repr__` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class BallAndStick:\n",
    "    def __init__(self):\n",
    "        self.soma = h.Section(name=\"soma\", cell=self)\n",
    "        self.dend = h.Section(name=\"dend\", cell=self)\n",
    "\n",
    "    def __repr__(self):\n",
    "        return \"BallAndStick\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h.topology()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our sections display based on the rule they were created with, not the new rule, so we need to recreate the cell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_cell = BallAndStick()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we see something friendlier:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h.topology()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There's a problem though; do you see it?\n",
    "\n",
    "Every cell using this rule would print out the same cell identifier. So if we create a second cell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_other_cell = BallAndStick()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and look at the topology:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h.topology()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is no way to tell the two somas apart and see which goes with which cell. To fix this, we'll pass in an identifier <tt>gid</tt> when we create the cell and have `__repr__` incorporate that into the name:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class BallAndStick:\n",
    "    def __init__(self, gid):\n",
    "        self._gid = gid\n",
    "        self.soma = h.Section(name=\"soma\", cell=self)\n",
    "        self.dend = h.Section(name=\"dend\", cell=self)\n",
    "\n",
    "    def __repr__(self):\n",
    "        return \"BallAndStick[{}]\".format(self._gid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Important: the cell name (returned by `__repr__`) will be read when the Section is created, so any and all data that function needs -- here the `gid` -- must be stored before creating any Section objects."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's create our two cells:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_cell = BallAndStick(0)\n",
    "my_other_cell = BallAndStick(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now when we look at the topology, we can tell the Sections belonging to the two cells apart:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h.topology()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can tell which sections are the soma and which are dendrites. We can see which go with cell 0 and which go with cell 1, but there is nothing indicating that the dendrite is connected to the soma. This is because we have not told NEURON about any such relationships, so let's do so:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Connect the sections"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We attach `self.dend` to the `self.soma` using the `connect` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class BallAndStick:\n",
    "    def __init__(self, gid):\n",
    "        self._gid = gid\n",
    "        self.soma = h.Section(name=\"soma\", cell=self)\n",
    "        self.dend = h.Section(name=\"dend\", cell=self)\n",
    "        self.dend.connect(self.soma)\n",
    "\n",
    "    def __repr__(self):\n",
    "        return \"BallAndStick[{}]\".format(self._gid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As before, we must recreate the cells now that we've changed the rule:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_cell = BallAndStick(0)\n",
    "my_other_cell = BallAndStick(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this is not equivalent to attaching the `soma` to the `dend`; instead it means that the dendrite begins where the soma ends. We can see that in the topology:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h.topology()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For now, we can get rid of `my_other_cell`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "del my_other_cell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h.topology()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A `soma` can have many dendrites attached to it, but any dendrite only begins at one specific location."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What if we didn't want to attach the dendrite to the end of the soma (position 1)? We could explicitly specify the connection location via, e.g. `self.dend.connect(self.soma(0.5))` which would mean the dendrite was attached to the center of the soma. (Recall: segments are described using normalized positions.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The 1 end of the soma is said to be the parent of the `dend` section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define stylized geometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's set the length and the width for both sections. We'll make the soma have length and diameter of 12.6157 microns, the dendrite have length 200 microns and diameter 1 micron."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class BallAndStick:\n",
    "    def __init__(self, gid):\n",
    "        self._gid = gid\n",
    "        self.soma = h.Section(name=\"soma\", cell=self)\n",
    "        self.dend = h.Section(name=\"dend\", cell=self)\n",
    "        self.dend.connect(self.soma)\n",
    "        self.soma.L = self.soma.diam = 12.6157 * µm\n",
    "        self.dend.L = 200 * µm\n",
    "        self.dend.diam = 1 * µm\n",
    "\n",
    "    def __repr__(self):\n",
    "        return \"BallAndStick[{}]\".format(self._gid)\n",
    "\n",
    "\n",
    "my_cell = BallAndStick(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you're wondering why that number was chosen for the soma length and diameter: it is because it makes the surface area (which doesn't include end faces) approximately 500 μm<sup>2</sup>:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_cell.soma(0.5).area()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(NEURON only returns areas of segments which is why we asked for the center `soma` segment; since there is only one segment, the area here is the same as the whole Section area.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the surface area of a cylinder with equal length and diameter\n",
    "$$\n",
    "\\pi d \\ell = \\pi d^2 = 4 \\pi \\left (\\frac{d}{2} \\right) ^2 = 4 \\pi r^2\n",
    "$$\n",
    "is the same as the surface area of a sphere with the same diameter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is, if we're only interested in electrophysiology modeling, we can substitute a cylindrical soma with equal length and diameter for a spherical soma with the same diameter, as we've done here. (The volume, however, is of course different. So this substitution does not work if we're modeling diffusion or accumulation of ions.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the rest of this tutorial page, we'll focus solely on one cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h.topology()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now it is time to see what the cell looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "h.PlotShape(False).plot(plt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also make an interactive shape plot in a separate window using NEURON's built-in graphics, via:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# enable NEURON's graphics\n",
    "from neuron import gui\n",
    "\n",
    "# here: True means show using NEURON's GUI; False means do not do so, at least not at first\n",
    "ps = h.PlotShape(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Either way, you will notice that this looks like a line instead of a ball and stick. Why? It's because NEURON by default does not display diameters. This behavior is useful when we need to see the structure of small dendrites, and in NEURON 7.7, it's the only supported option for Jupyter notebooks with h.PlotShape (although gui2.PlotShape provides similar functionality)... but when using NEURON's built-in graphics, we can use the <tt>show</tt> method to show diamters via:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Why does this look like a line instead of a ball and stick? It's because NEURON by default does not display diameters. This behavior is useful when we need to see the structure of small dendrites, but for now, let's show the diameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ps = h.PlotShape(True)\n",
    "ps.show(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Jupyter, we can rotate images by clicking and dragging; we can zoom by right-clicking and dragging. When using NEURON's built-in graphics which appear in separate windows, right-click and select \"3D Rotate\", then drag to rotate. For this simple morphology, there is not anything else more to see."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Note: If you want to build your own visualization tool using matplotlib, mayavi, etc, ensure 3D points exist with `h.define_shape()`, then loop over all the sections with `h.allsec()` and read the morphology using `sec.x3d(i)` etc and `sec.diam3d(i)` for `i` in 0, .., `sec.n3d() - 1`. Less efficiently, the (x, y, z; diam) values for a whole section may be read by `sec.psection()['morphology']['pt3d']`.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Specify biophysics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our cell needs biophysical mechanisms in the membrane. We start by setting axial resistance and membrane capacitance. (Recall: NEURON's default axial resistance is appropriate for squid but low for mammalian models.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class BallAndStick:\n",
    "    def __init__(self, gid):\n",
    "        self._gid = gid\n",
    "        self.soma = h.Section(name=\"soma\", cell=self)\n",
    "        self.dend = h.Section(name=\"dend\", cell=self)\n",
    "        self.dend.connect(self.soma)\n",
    "        self.all = self.soma.wholetree()  # <-- NEW\n",
    "        self.soma.L = self.soma.diam = 12.6157 * µm\n",
    "        self.dend.L = 200 * µm\n",
    "        self.dend.diam = 1 * µm\n",
    "        for sec in self.all:  # <-- NEW\n",
    "            sec.Ra = 100  # Axial resistance in Ohm * cm                    # <-- NEW\n",
    "            sec.cm = 1  # Membrane capacitance in micro Farads / cm^2     # <-- NEW\n",
    "\n",
    "    def __repr__(self):\n",
    "        return \"BallAndStick[{}]\".format(self._gid)\n",
    "\n",
    "\n",
    "my_cell = BallAndStick(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've added a new variable `self.all` which is a list of all the sections in the cell. The <tt>wholetree</tt> method of a Section returns a list of all the sections it is attached to -- i.e. the corresponding neuron. This will help us iterate over them to -- in this case -- specify axial resistance and membrane capacitance, but can also be used for any other biophysics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is getting a little complicated. Let's split `__init__` into several functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class BallAndStick:\n",
    "    def __init__(self, gid):\n",
    "        self._gid = gid\n",
    "        self._setup_morphology()\n",
    "        self._setup_biophysics()\n",
    "\n",
    "    def _setup_morphology(self):\n",
    "        self.soma = h.Section(name=\"soma\", cell=self)\n",
    "        self.dend = h.Section(name=\"dend\", cell=self)\n",
    "        self.all = [self.soma, self.dend]\n",
    "        self.dend.connect(self.soma)\n",
    "        self.soma.L = self.soma.diam = 12.6157 * µm\n",
    "        self.dend.L = 200 * µm\n",
    "        self.dend.diam = 1\n",
    "\n",
    "    def _setup_biophysics(self):\n",
    "        for sec in self.all:\n",
    "            sec.Ra = 100  # Axial resistance in Ohm * cm\n",
    "            sec.cm = 1  # Membrane capacitance in micro Farads / cm^2\n",
    "\n",
    "    def __repr__(self):\n",
    "        return \"BallAndStick[{}]\".format(self._gid)\n",
    "\n",
    "\n",
    "my_cell = BallAndStick(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It may not seem like we gained a lot by doing this, but it makes the code a little more self-documenting (e.g. we can see at a glance which parts have to do with defining the morphology) and it divides it into pieces that could be reused in other cells."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll put Hodgkin-Huxley (`hh`) kinetics in the soma and specify some parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class BallAndStick:\n",
    "    def __init__(self, gid):\n",
    "        self._gid = gid\n",
    "        self._setup_morphology()\n",
    "        self._setup_biophysics()\n",
    "\n",
    "    def _setup_morphology(self):\n",
    "        self.soma = h.Section(name=\"soma\", cell=self)\n",
    "        self.dend = h.Section(name=\"dend\", cell=self)\n",
    "        self.dend.connect(self.soma)\n",
    "        self.all = self.soma.wholetree()\n",
    "        self.soma.L = self.soma.diam = 12.6157 * µm\n",
    "        self.dend.L = 200 * µm\n",
    "        self.dend.diam = 1 * µm\n",
    "\n",
    "    def _setup_biophysics(self):\n",
    "        for sec in self.all:\n",
    "            sec.Ra = 100  # Axial resistance in Ohm * cm\n",
    "            sec.cm = 1  # Membrane capacitance in micro Farads / cm^2\n",
    "        self.soma.insert(\"hh\")  # <-- NEW\n",
    "        for seg in self.soma:  # <-- NEW\n",
    "            seg.hh.gnabar = (\n",
    "                0.12  # Sodium conductance in S/cm2                  # <-- NEW\n",
    "            )\n",
    "            seg.hh.gkbar = (\n",
    "                0.036  # Potassium conductance in S/cm2               # <-- NEW\n",
    "            )\n",
    "            seg.hh.gl = 0.0003  # Leak conductance in S/cm2                    # <-- NEW\n",
    "            seg.hh.el = (\n",
    "                -54.3 * mV\n",
    "            )  # Reversal potential                           # <-- NEW\n",
    "\n",
    "    def __repr__(self):\n",
    "        return \"BallAndStick[{}]\".format(self._gid)\n",
    "\n",
    "\n",
    "my_cell = BallAndStick(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Note: here we loop over all segments in the soma, even though we only defined one segment. This gives us more general code, that will still work if we change the number of segments later. Always write your model implementations to be independent of the discretization.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let's insert a passive (leak) current in the dendrite:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class BallAndStick:\n",
    "    def __init__(self, gid):\n",
    "        self._gid = gid\n",
    "        self._setup_morphology()\n",
    "        self._setup_biophysics()\n",
    "\n",
    "    def _setup_morphology(self):\n",
    "        self.soma = h.Section(name=\"soma\", cell=self)\n",
    "        self.dend = h.Section(name=\"dend\", cell=self)\n",
    "        self.dend.connect(self.soma)\n",
    "        self.all = self.soma.wholetree()\n",
    "        self.soma.L = self.soma.diam = 12.6157 * µm\n",
    "        self.dend.L = 200 * µm\n",
    "        self.dend.diam = 1 * µm\n",
    "\n",
    "    def _setup_biophysics(self):\n",
    "        for sec in self.all:\n",
    "            sec.Ra = 100  # Axial resistance in Ohm * cm\n",
    "            sec.cm = 1  # Membrane capacitance in micro Farads / cm^2\n",
    "        self.soma.insert(\"hh\")\n",
    "        for seg in self.soma:\n",
    "            seg.hh.gnabar = 0.12  # Sodium conductance in S/cm2\n",
    "            seg.hh.gkbar = 0.036  # Potassium conductance in S/cm2\n",
    "            seg.hh.gl = 0.0003  # Leak conductance in S/cm2\n",
    "            seg.hh.el = -54.3 * mV  # Reversal potential\n",
    "        # Insert passive current in the dendrite                       # <-- NEW\n",
    "        self.dend.insert(\"pas\")  # <-- NEW\n",
    "        for seg in self.dend:  # <-- NEW\n",
    "            seg.pas.g = 0.001  # Passive conductance in S/cm2        # <-- NEW\n",
    "            seg.pas.e = -65 * mV  # Leak reversal potential             # <-- NEW\n",
    "\n",
    "    def __repr__(self):\n",
    "        return \"BallAndStick[{}]\".format(self._gid)\n",
    "\n",
    "\n",
    "my_cell = BallAndStick(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we weren't sure about the units for a given mechanism’s parameter, use <a class=\"reference external\" href=\"https://neuron.yale.edu/neuron/static/py_doc/programming/gui/widgets.html#units\" title=\"(in NEURON v7.5)\"><code class=\"xref py py-func docutils literal\"><span class=\"pre\">units()</span></code></a>. Pass in a string with the paramater name, an underscore, and then the mechanism name. e.g."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(h.units(\"gnabar_hh\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use `psection` to see what is present where:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for sec in h.allsec():\n",
    "    print(\"%s: %s\" % (sec, \", \".join(sec.psection()[\"density_mechs\"].keys())))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(A quick aside about `psection`: it's great for quickly getting information when interactively exploring a model because it provides a lot of data in one pass; for the same reason, however, other solutions that only extract specific desired information are more efficient for automatic explorations.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Instrumentation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have now created our neuron. Let's stimulate it and visualize its dynamics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stimulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll inject a current pulse into the distal (1) end of the dendrite starting 5 ms after the simulation starts, with a duration of 1 ms, and an amplitude of 0.1 nA. First, let's define and position the current clamp object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stim = h.IClamp(my_cell.dend(1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want, we can check the segment the current clamp is inserted into:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stim.get_segment()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that if we forget what the names of the attributes are, we can check the `dir`. Here we do that and ignore Python-created attributes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\", \".join(item for item in dir(stim) if not item.startswith(\"__\")))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The parameters we need to set for a current clamp are `delay` (measured in ms), `dur` (measured in ms), and `amp` (in nA):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stim.delay = 5\n",
    "stim.dur = 1\n",
    "stim.amp = 0.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Recording"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll start out recording the membrane potential at the center of the soma and the time in two NEURON Vectors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "soma_v = h.Vector().record(my_cell.soma(0.5)._ref_v)\n",
    "t = h.Vector().record(h._ref_t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run the simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We initialize membrane potential everywhere to -65 mV:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h.finitialize(-65 * mV)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we run until time 25 ms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h.continuerun(25 * ms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As in the scripting neuron basics part of the tutorial, we initialize `bokeh` graphics:"
   ]
  },
  {
   "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": [
    "(If you prefer to use `matplotlib` graphics, you can adapt the code below using the examples in the scripting neuron basics part of the tutorial.)"
   ]
  },
  {
   "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, soma_v, line_width=2)\n",
    "plt.show(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have a cell that responds to input, so that's good, let's explore the role of some parameters to see what they do and see if we can get an action potential."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation studies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Role of current amplitude"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's do a set of simulations, plotted in the same figure, where we vary the amplitude of the current in a `for` loop:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = plt.figure(x_axis_label=\"t (ms)\", y_axis_label=\"v (mV)\")\n",
    "amps = [0.075 * i for i in range(1, 5)]\n",
    "colors = [\"green\", \"blue\", \"red\", \"black\"]\n",
    "for amp, color in zip(amps, colors):\n",
    "    stim.amp = amp\n",
    "    h.finitialize(-65 * mV)\n",
    "    h.continuerun(25 * ms)\n",
    "    f.line(t, list(soma_v), line_width=2, legend_label=\"amp=%g\" % amp, color=color)\n",
    "plt.show(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `for` loop syntax as used here is looping over amplitude, color pairs. The only potential \"gotcha\" here is the need to make a copy of the values by passing in `list(soma_v)` instead of `soma_v`. If this copying was omitted, only the last set of values would be plotted."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting both the dendrite and the soma"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To plot the dendrite potential, we need to record it in a NEURON Vector:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dend_v = h.Vector().record(my_cell.dend(0.5)._ref_v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we do not need to reassign the existing time and soma membrane potential recording vectors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The running and plotting code is very similar to the above; the only difference is one additional call to the `line` method to plot the `dend_v` using a dashed line without a separate legend entry:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = plt.figure(x_axis_label=\"t (ms)\", y_axis_label=\"v (mV)\")\n",
    "amps = [0.075 * i for i in range(1, 5)]\n",
    "colors = [\"green\", \"blue\", \"red\", \"black\"]\n",
    "for amp, color in zip(amps, colors):\n",
    "    stim.amp = amp\n",
    "    h.finitialize(-65 * mV)\n",
    "    h.continuerun(25 * ms)\n",
    "    f.line(t, list(soma_v), line_width=2, legend_label=\"amp=%g\" % amp, color=color)\n",
    "    f.line(t, list(dend_v), line_width=2, line_dash=\"dashed\", color=color)\n",
    "\n",
    "plt.show(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that when the soma membrane potential is sufficiently low, it is possible for the dendrite to be more depolarized than the soma, but the peak membrane potential in the leaky dendrite is significantly reduced from the peak in the soma during an action potential. (This is because in our model all the voltage-gated channels that would depolarize the cell are only present in the soma.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The role of <tt>nseg</tt>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see the effects of nseg, the number of segments of the dendrite, on the signal through the dendrite."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll do this by modifying the above code to run all of our simulations for both `nseg = 1` (the default) and for `nseg = 101`. We'll use thin lines for the high resolution (`nseg = 101`) case and thick lines for the low resolution case:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = plt.figure(x_axis_label=\"t (ms)\", y_axis_label=\"v (mV)\")\n",
    "amps = [0.075 * i for i in range(1, 5)]\n",
    "colors = [\"green\", \"blue\", \"red\", \"black\"]\n",
    "for amp, color in zip(amps, colors):\n",
    "    stim.amp = amp\n",
    "    for my_cell.dend.nseg, width in [(1, 2), (101, 1)]:\n",
    "        h.finitialize(-65)\n",
    "        h.continuerun(25)\n",
    "        f.line(\n",
    "            t,\n",
    "            list(soma_v),\n",
    "            line_width=width,\n",
    "            legend_label=\"amp={:.3g}\".format(amp) if my_cell.dend.nseg == 1 else \"\",\n",
    "            color=color,\n",
    "        )\n",
    "        f.line(t, list(dend_v), line_width=width, line_dash=\"dashed\", color=color)\n",
    "plt.show(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This code has only two changes to the previous code: (1) addition of the `for` loop for `dend.nseg` values and width with accompanying indentation changes, (2) a modification to the legend formula to only generate a legend for the `nseg=1` case since the colors are unchanged. (The calls to `f.line` are split onto multiple lines but this is only for readability; there is no difference to Python if it is on one line or on multiple lines.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we see with the high-resolution simulation that the soma peaks should be reduced and delayed and the dendrite peaks increased relative to what was seen in the `nseg=1` case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "Modify the above example, increasing the low resolution case from `nseg=1` to find a low value of `nseg` that gives negligible error relative to the high-resolution `nseg=101` case. Hint: a non-quantitative way to proceed would be to find values such that all the thin and thick curves of a given color overlap. (Note that since we're plotting from the center of the dendrite, `nseg` should be odd, otherwise the center will fall on the segment boundaries and be ill-defined.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Advanced exercise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the approximate absolute error for your low `nseg` solution that graphically overlaps the high `nseg` solution. Compare this to the error for the `nseg=1` case. (Since the true solution is unknown, approximate it with the high `nseg` solution.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}