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

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ball and stick 3: Extensible network of cells"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the third part of a tutorial series where we build a multicompartment cell and evolve it into a network of cells running on a parallel machine. In this part, we take the functionality of the ring network we constructed in the previous page and encapsulate it into various classes so that the network is more extensible. We also begin parameterizing the model so that particular values are not hard-coded, but remain variable so that the model is flexible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading libraries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As before, we will begin by loading relevant NEURON libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from neuron import h, gui\n",
    "from neuron.units import ms, mV\n",
    "\n",
    "h.load_file(\"stdrun.hoc\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generic Cell class"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the last tutorial, we created a generic `Cell` class (actually, two versions) but we can expand this to make it more powerful. For example, let's make each Cell record its spike times, some membrane potential timeseries, and keep track of NetCons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Cell:\n",
    "    def __init__(self, gid, x, y, z, theta):\n",
    "        self._gid = gid\n",
    "        self._setup_morphology()\n",
    "        self.all = self.soma.wholetree()\n",
    "        self._setup_biophysics()\n",
    "        self.x = self.y = self.z = 0\n",
    "        h.define_shape()\n",
    "        self._rotate_z(theta)\n",
    "        self._set_position(x, y, z)\n",
    "\n",
    "        # everything below here in this method is NEW\n",
    "        self._spike_detector = h.NetCon(self.soma(0.5)._ref_v, None, sec=self.soma)\n",
    "        self.spike_times = h.Vector()\n",
    "        self._spike_detector.record(self.spike_times)\n",
    "\n",
    "        self._ncs = []\n",
    "\n",
    "        self.soma_v = h.Vector().record(self.soma(0.5)._ref_v)\n",
    "\n",
    "    def __repr__(self):\n",
    "        return \"{}[{}]\".format(self.name, self._gid)\n",
    "\n",
    "    def _set_position(self, x, y, z):\n",
    "        for sec in self.all:\n",
    "            for i in range(sec.n3d()):\n",
    "                sec.pt3dchange(\n",
    "                    i,\n",
    "                    x - self.x + sec.x3d(i),\n",
    "                    y - self.y + sec.y3d(i),\n",
    "                    z - self.z + sec.z3d(i),\n",
    "                    sec.diam3d(i),\n",
    "                )\n",
    "        self.x, self.y, self.z = x, y, z\n",
    "\n",
    "    def _rotate_z(self, theta):\n",
    "        \"\"\"Rotate the cell about the Z axis.\"\"\"\n",
    "        for sec in self.all:\n",
    "            for i in range(sec.n3d()):\n",
    "                x = sec.x3d(i)\n",
    "                y = sec.y3d(i)\n",
    "                c = h.cos(theta)\n",
    "                s = h.sin(theta)\n",
    "                xprime = x * c - y * s\n",
    "                yprime = x * s + y * c\n",
    "                sec.pt3dchange(i, xprime, yprime, sec.z3d(i), sec.diam3d(i))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the `BallAndStick` has a simple geometry, we could modify it to assume that all inputs go into a single location; we will call it the `.syn`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class BallAndStick(Cell):\n",
    "    name = \"BallAndStick\"\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.soma.L = self.soma.diam = 12.6157\n",
    "        self.dend.L = 200\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",
    "        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  # Reversal potential in mV\n",
    "        # Insert passive current in the dendrite\n",
    "        self.dend.insert(\"pas\")\n",
    "        for seg in self.dend:\n",
    "            seg.pas.g = 0.001  # Passive conductance in S/cm2\n",
    "            seg.pas.e = -65  # Leak reversal potential mV\n",
    "\n",
    "        # NEW: the synapse\n",
    "        self.syn = h.ExpSyn(self.dend(0.5))\n",
    "        self.syn.tau = 2 * ms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Due to the nature of `h.ExpSyn` decay, there is mathematically no difference between having two ExpSyn objects at the same point or one synapse where multiple inputs add linearly, so it suffices to have just the one as long as we're happy with all inputs going into `dend(0.5)`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make a Ring class"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Encapsulating code into discrete objects is not only conceptually useful for code management, but as we know with cell objects, it lets us make several instances of the object for use in a network. Thinking ahead, we may very well need several networks – each network configured differently. This allows scripting of several simulations <i>en masse</i>, either in a for loop that sequentially processes the networks, or it can be used with NEURON's <a href=\"https://neuron.yale.edu/neuron/static/py_doc/modelspec/programmatic/network/parcon.html#ParallelContext.subworlds\">subworlds</a> architecture in a parallel context."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Ring:\n",
    "    \"\"\"A network of *N* ball-and-stick cells where cell n makes an\n",
    "    excitatory synapse onto cell n + 1 and the last, Nth cell in the\n",
    "    network projects to the first cell.\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(\n",
    "        self, N=5, stim_w=0.04, stim_t=9, stim_delay=1, syn_w=0.01, syn_delay=5, r=50\n",
    "    ):\n",
    "        \"\"\"\n",
    "        :param N: Number of cells.\n",
    "        :param stim_w: Weight of the stimulus\n",
    "        :param stim_t: time of the stimulus (in ms)\n",
    "        :param stim_delay: delay of the stimulus (in ms)\n",
    "        :param syn_w: Synaptic weight\n",
    "        :param syn_delay: Delay of the synapse\n",
    "        :param r: radius of the network\n",
    "        \"\"\"\n",
    "        self._syn_w = syn_w\n",
    "        self._syn_delay = syn_delay\n",
    "        self._create_cells(N, r)\n",
    "        self._connect_cells()\n",
    "        # add stimulus\n",
    "        self._netstim = h.NetStim()\n",
    "        self._netstim.number = 1\n",
    "        self._netstim.start = stim_t\n",
    "        self._nc = h.NetCon(self._netstim, self.cells[0].syn)\n",
    "        self._nc.delay = stim_delay\n",
    "        self._nc.weight[0] = stim_w\n",
    "\n",
    "    def _create_cells(self, N, r):\n",
    "        self.cells = []\n",
    "        for i in range(N):\n",
    "            theta = i * 2 * h.PI / N\n",
    "            self.cells.append(\n",
    "                BallAndStick(i, h.cos(theta) * r, h.sin(theta) * r, 0, theta)\n",
    "            )\n",
    "\n",
    "    def _connect_cells(self):\n",
    "        for source, target in zip(self.cells, self.cells[1:] + [self.cells[0]]):\n",
    "            nc = h.NetCon(source.soma(0.5)._ref_v, target.syn, sec=source.soma)\n",
    "            nc.weight[0] = self._syn_w\n",
    "            nc.delay = self._syn_delay\n",
    "            source._ncs.append(nc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `_create_cells` method is basically the same as the `create_n_BallAndStick` function in the previous part of the tutorial; the only difference is that the cells are stored in `self._cells` instead of being returned. `_connect_cells` is shorter than the previous version because it can take advantage of the existing synapses and lists."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test the network"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's make a `Ring` object with 5 cells, render it using NEURON's built-in graphics, and run a simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ring = Ring(N=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now to check that it is constructed correctly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shape_window = h.PlotShape(True)\n",
    "shape_window.show(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks good so far; let's run the simulation and record time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = h.Vector().record(h._ref_t)\n",
    "h.finitialize(-65 * mV)\n",
    "h.continuerun(100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember that if we are running in Jupyter to make a plot appear inline we must:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now plot the trace of cell 0's soma:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.plot(t, ring.cells[0].soma_v)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cell 0 looks good. Let's look at the raster diagram:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "for i, cell in enumerate(ring.cells):\n",
    "    plt.vlines(list(cell.spike_times), i + 0.5, i + 1.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explore effects of parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare two simulations: one with the same parameters as above, which we'll plot in black, and one with half the synaptic weight, which we'll plot in red:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "for syn_w, color in [(0.01, \"black\"), (0.005, \"red\")]:\n",
    "    ring = Ring(N=5, syn_w=syn_w)\n",
    "    h.finitialize(-65 * mV)\n",
    "    h.continuerun(100 * ms)\n",
    "    for i, cell in enumerate(ring.cells):\n",
    "        plt.vlines(list(cell.spike_times), i + 0.5, i + 1.5, color=color)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In both simulations, the first spike occurs at 10.925 ms. After that, the red spikes lag the black ones by steadily increasing amounts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next part of the tutorial will translate this serial model into a parallel model. That part will not work in Jupyter and must be run from a terminal."
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}