network_ring.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Qe5XtIdg7VcU"
   },
   "source": [
    "## Are you ready to install Arbor?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Uncomment for Google Collab\n",
    "# !pip install pandas seaborn arbor\n",
    "# !mv /usr/local/lib/python3.7/site-packages/arbor /usr/local/lib/python3.7/dist-packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "iL1fXlbz7Whp",
    "outputId": "685a9459-6779-4c17-d15b-fbd859b398d7"
   },
   "outputs": [],
   "source": [
    "import arbor, pandas, seaborn\n",
    "from math import sqrt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WEdy_Twg9ZIs"
   },
   "source": [
    "# A ring network"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "uv9gW2to7txz"
   },
   "source": [
    "In this example, a small network of cells, arranged in a ring, will be created and the simulation distributed over multiple threads or GPUs if available.\n",
    "\n",
    "Concepts covered in this example:\n",
    "* Building a basic `arbor.cell` with a synapse site and spike generator.\n",
    "* Building a `arbor.recipe` with a network of interconnected cells.\n",
    "* Running the simulation and extract the results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "JaCLO1my8Z9K"
   },
   "source": [
    "## The cell\n",
    "\n",
    "![tutorial_network_ring_morph.svg](https://docs.arbor-sim.org/en/latest/_images/tutorial_network_ring_morph.svg)\n",
    "_A 4-segment cell with a soma (pink) and a branched dendrite (light blue)._\n",
    "\n",
    "**Step (1)** shows a function that creates a simple cell with a dendrite. We construct the following morphology and label the soma and dendrite:\n",
    "\n",
    "```python\n",
    "def make_cable_cell(gid):\n",
    "    # (1) Build a segment tree\n",
    "    tree = arbor.segment_tree()\n",
    "\n",
    "    # Soma (tag=1) with radius 6 μm, modelled as cylinder of length 2*radius\n",
    "    s = tree.append(\n",
    "        arbor.mnpos, arbor.mpoint(-12, 0, 0, 6), arbor.mpoint(0, 0, 0, 6), tag=1\n",
    "    )\n",
    "\n",
    "    # (b0) Single dendrite (tag=3) of length 50 μm and radius 2 μm attached to soma.\n",
    "    b0 = tree.append(s, arbor.mpoint(0, 0, 0, 2), arbor.mpoint(50, 0, 0, 2), tag=3)\n",
    "\n",
    "    # Attach two dendrites (tag=3) of length 50 μm to the end of the first dendrite.\n",
    "    # (b1) Radius tapers from 2 to 0.5 μm over the length of the dendrite.\n",
    "    tree.append(\n",
    "        b0,\n",
    "        arbor.mpoint(50, 0, 0, 2),\n",
    "        arbor.mpoint(50 + 50 / sqrt(2), 50 / sqrt(2), 0, 0.5),\n",
    "        tag=3,\n",
    "    )\n",
    "    # (b2) Constant radius of 1 μm over the length of the dendrite.\n",
    "    tree.append(\n",
    "        b0,\n",
    "        arbor.mpoint(50, 0, 0, 1),\n",
    "        arbor.mpoint(50 + 50 / sqrt(2), -50 / sqrt(2), 0, 1),\n",
    "        tag=3,\n",
    "    )\n",
    "\n",
    "    # Associate labels to tags\n",
    "    labels = arbor.label_dict()\n",
    "    labels[\"soma\"] = \"(tag 1)\"\n",
    "    labels[\"dend\"] = \"(tag 3)\"\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "KQmM8K4q9KQ4"
   },
   "source": [
    "In **step (2)** we create a label for both the root and the site of the synapse. These locations will form the endpoints of the connections between the cells.\n",
    "\n",
    "![](https://docs.arbor-sim.org/en/latest/_images/tutorial_network_ring_synapse_site.svg)\n",
    "_We’ll create labels for the root (red) and a synapse_site (black)._\n",
    "\n",
    "```python\n",
    "    # (2) Mark location for synapse at the midpoint of branch 1 (the first dendrite).\n",
    "    labels['synapse_site'] = '(location 1 0.5)'\n",
    "    # Mark the root of the tree.\n",
    "    labels['root'] = '(root)'\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "lzunVL9X-tJk"
   },
   "source": [
    "After we’ve created a basic `arbor.decor`, **step (3)** places a synapse with an exponential decay (`'expsyn'`) on the `'synapse_site'`. The synapse is given the label `'syn'`, which is later used to form `arbor.connection` objects terminating at the cell.\n",
    "\n",
    "**Step (4)** places a spike detector at the 'root'. The detector is given the label 'detector', which is later used to form `arbor.connection` objects originating from the cell.\n",
    "\n",
    "```python\n",
    "    # (3) Create a decor and a cable_cell\n",
    "    decor = arbor.decor()\n",
    "\n",
    "    # Put hh dynamics on soma, and passive properties on the dendrites.\n",
    "    decor.paint('\"soma\"', arbor.density(\"hh\"))\n",
    "    decor.paint('\"dend\"', arbor.density(\"pas\"))\n",
    "\n",
    "    # (4) Attach a single synapse.\n",
    "    decor.place('\"synapse_site\"', arbor.synapse(\"expsyn\"), \"syn\")\n",
    "\n",
    "    # Attach a spike detector with threshold of -10 mV.\n",
    "    decor.place('\"root\"', arbor.spike_detector(-10), \"detector\")\n",
    "\n",
    "    cell = arbor.cable_cell(tree, labels, decor)\n",
    "\n",
    "    return cell\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "27_aD44W_NbM"
   },
   "source": [
    "Tying it all together:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "id": "vldZLHrX9FNX"
   },
   "outputs": [],
   "source": [
    "def make_cable_cell(gid):\n",
    "    # (1) Build a segment tree\n",
    "    tree = arbor.segment_tree()\n",
    "\n",
    "    # Soma (tag=1) with radius 6 μm, modelled as cylinder of length 2*radius\n",
    "    s = tree.append(\n",
    "        arbor.mnpos, arbor.mpoint(-12, 0, 0, 6), arbor.mpoint(0, 0, 0, 6), tag=1\n",
    "    )\n",
    "\n",
    "    # (b0) Single dendrite (tag=3) of length 50 μm and radius 2 μm attached to soma.\n",
    "    b0 = tree.append(s, arbor.mpoint(0, 0, 0, 2), arbor.mpoint(50, 0, 0, 2), tag=3)\n",
    "\n",
    "    # Attach two dendrites (tag=3) of length 50 μm to the end of the first dendrite.\n",
    "    # (b1) Radius tapers from 2 to 0.5 μm over the length of the dendrite.\n",
    "    tree.append(\n",
    "        b0,\n",
    "        arbor.mpoint(50, 0, 0, 2),\n",
    "        arbor.mpoint(50 + 50 / sqrt(2), 50 / sqrt(2), 0, 0.5),\n",
    "        tag=3,\n",
    "    )\n",
    "    # (b2) Constant radius of 1 μm over the length of the dendrite.\n",
    "    tree.append(\n",
    "        b0,\n",
    "        arbor.mpoint(50, 0, 0, 1),\n",
    "        arbor.mpoint(50 + 50 / sqrt(2), -50 / sqrt(2), 0, 1),\n",
    "        tag=3,\n",
    "    )\n",
    "\n",
    "    # Associate labels to tags\n",
    "    labels = arbor.label_dict()\n",
    "    labels[\"soma\"] = \"(tag 1)\"\n",
    "    labels[\"dend\"] = \"(tag 3)\"\n",
    "\n",
    "    # (2) Mark location for synapse at the midpoint of branch 1 (the first dendrite).\n",
    "    labels[\"synapse_site\"] = \"(location 1 0.5)\"\n",
    "    # Mark the root of the tree.\n",
    "    labels[\"root\"] = \"(root)\"\n",
    "\n",
    "    # (3) Create a decor and a cable_cell\n",
    "    decor = arbor.decor()\n",
    "\n",
    "    # Put hh dynamics on soma, and passive properties on the dendrites.\n",
    "    decor.paint('\"soma\"', arbor.density(\"hh\"))\n",
    "    decor.paint('\"dend\"', arbor.density(\"pas\"))\n",
    "\n",
    "    # (4) Attach a single synapse.\n",
    "    decor.place('\"synapse_site\"', arbor.synapse(\"expsyn\"), \"syn\")\n",
    "\n",
    "    # Attach a spike detector with threshold of -10 mV.\n",
    "    decor.place('\"root\"', arbor.spike_detector(-10), \"detector\")\n",
    "\n",
    "    cell = arbor.cable_cell(tree, labels, decor)\n",
    "\n",
    "    return cell"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "3OyFuxfxBaAn"
   },
   "source": [
    "## The recipe\n",
    "\n",
    "To create a model with multiple connected cells, we need to use a `arbor.recipe`. The recipe is where the different cells and the connections between them are defined.\n",
    "\n",
    "Step **(5)** shows a class definition for a recipe with multiple cells. Instantiating the class requires the desired number of cells as input. We are connecting the cells **(8)**, returning a configurable number of cells **(6)** and returning a new cell per `gid` **(7)**.\n",
    "\n",
    "Step **(8)** creates an `arbor.connection` between consecutive cells. If a cell has gid `gid`, the previous cell has a gid `(gid-1)%self.ncells`. The connection has a weight of 0.01 (inducing a conductance of 0.01 μS in the target mechanism `expsyn`) and a delay of 5 ms. The first two arguments to `arbor.connection` are the **source** and **target** of the connection.\n",
    "\n",
    "The **source** is a `arbor.cell_global_label` object containing a cell index `gid`, the source label corresponding to a valid detector label on the cell and an optional selection policy (for choosing a single detector out of potentially many detectors grouped under the same label - remember, in this case the number of detectors labeled 'detector' is 1). The `arbor.cell_global_label` can be initialized with a `(gid, label)` tuple, in which case the selection policy is the default `arbor.selection_policy.univalent`; or a `(gid, (label, policy))` tuple.\n",
    "\n",
    "The **target** is a `arbor.cell_local_label` object containing a cell index `gid`, the target label corresponding to a valid synapse label on the cell and an optional selection policy (for choosing a single synapse out of potentially many synapses grouped under the same label - remember, in this case the number of synapses labeled 'syn' is 1). The `arbor.cell_local_label` can be initialized with a `label` string, in which case the selection policy is the default `arbor.selection_policy.univalent`; or a `(label, policy)` tuple. The `gid` of the target cell doesn't need to be explicitly added to the connection, it is the argument to the `arbor.recipe.connections_on` method.\n",
    "\n",
    "Step **(9)** attaches an `arbor.event_generator` on the 0th target (synapse) on the 0th cell; this means it is connected to the `\"synapse_site\"` on cell 0. This initiates the signal cascade through the network. The `arbor.explicit_schedule` in instantiated with a list of times in milliseconds, so here a single event at the 1 ms mark is emitted. Note that this synapse is connected twice, once to the event generator, and once to another cell.\n",
    "\n",
    "Step **(10)** places a `probe` at the `\"root\"` of each cell.\n",
    "\n",
    "Step **(11)** instantiates the recipe with 4 cells.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "id": "odQw1y-kEQ-W"
   },
   "outputs": [],
   "source": [
    "# (5) Create a recipe that generates a network of connected cells.\n",
    "class ring_recipe(arbor.recipe):\n",
    "    def __init__(self, ncells):\n",
    "        # The base C++ class constructor must be called first, to ensure that\n",
    "        # all memory in the C++ class is initialized correctly.\n",
    "        arbor.recipe.__init__(self)\n",
    "        self.ncells = ncells\n",
    "        self.props = arbor.neuron_cable_properties()\n",
    "\n",
    "    # (6) The num_cells method that returns the total number of cells in the model\n",
    "    # must be implemented.\n",
    "    def num_cells(self):\n",
    "        return self.ncells\n",
    "\n",
    "    # (7) The cell_description method returns a cell\n",
    "    def cell_description(self, gid):\n",
    "        return make_cable_cell(gid)\n",
    "\n",
    "    # The kind method returns the type of cell with gid.\n",
    "    # Note: this must agree with the type returned by cell_description.\n",
    "    def cell_kind(self, gid):\n",
    "        return arbor.cell_kind.cable\n",
    "\n",
    "    # (8) Make a ring network. For each gid, provide a list of incoming connections.\n",
    "    def connections_on(self, gid):\n",
    "        src = (gid - 1) % self.ncells\n",
    "        w = 0.01  # 0.01 μS on expsyn\n",
    "        d = 5  # ms delay\n",
    "        return [arbor.connection((src, \"detector\"), \"syn\", w, d)]\n",
    "\n",
    "    # (9) Attach a generator to the first cell in the ring.\n",
    "    def event_generators(self, gid):\n",
    "        if gid == 0:\n",
    "            sched = arbor.explicit_schedule([1])  # one event at 1 ms\n",
    "            weight = 0.1  # 0.1 μS on expsyn\n",
    "            return [arbor.event_generator(\"syn\", weight, sched)]\n",
    "        return []\n",
    "\n",
    "    # (10) Place a probe at the root of each cell.\n",
    "    def probes(self, gid):\n",
    "        return [arbor.cable_probe_membrane_voltage('\"root\"')]\n",
    "\n",
    "    def global_properties(self, kind):\n",
    "        return self.props\n",
    "\n",
    "\n",
    "# (11) Instantiate recipe\n",
    "ncells = 4\n",
    "recipe = ring_recipe(ncells)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "kl6w1y6xETtk"
   },
   "source": [
    "## The execution\n",
    "\n",
    "To create a simulation, we must create an `arbor.context` and\n",
    "`arbor.domain_decomposition`.\n",
    "\n",
    "Step **(12)** creates a default execution context, and uses the\n",
    "`arbor.partition_load_balance` to create a default domain decomposition.\n",
    "You can print the objects to see what defaults they produce on your\n",
    "system.\n",
    "\n",
    "Step **(13)** sets all spike generators to record using the\n",
    "`arbor.spike_recording.all` policy. This means the timestamps of the\n",
    "generated events will be kept in memory. Be default, these are\n",
    "discarded.\n",
    "\n",
    "In addition to having the timestamps of spikes, we want to extract the\n",
    "voltage as a function of time.\n",
    "\n",
    "Step **(14)** sets the probes (step **10**) to measure at a certain\n",
    "schedule. This is sometimes described as attaching a `sampler` to a\n",
    "`probe`. `arbor.simulation.sample` expects a `probe id` and the\n",
    "desired schedule (here: a recording frequency of 10 kHz, or a `dt` of\n",
    "0.1 ms). Note that the probe id is a separate index from those of\n",
    "`connection` endpoints; probe ids correspond to the index of the list\n",
    "produced by `arbor.recipe.probes` on cell `gid`.\n",
    "\n",
    "`arbor.simulation.sample` returns a handle to the samples\n",
    "that will be recorded. We store these handles for later use.\n",
    "\n",
    "Step **(15)** executes the simulation for a duration of 100 ms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "nLNsaf96EgNL",
    "outputId": "671d4930-d7d7-442d-b9d9-d342c117efd1"
   },
   "outputs": [],
   "source": [
    "# (12) Create an execution context using all locally available threads and simulation\n",
    "sim = arbor.simulation(recipe)\n",
    "\n",
    "# (13) Set spike generators to record\n",
    "sim.record(arbor.spike_recording.all)\n",
    "\n",
    "# (14) Attach a sampler to the voltage probe on cell 0. Sample rate of 10 sample every ms.\n",
    "handles = [sim.sample((gid, 0), arbor.regular_schedule(0.1)) for gid in range(ncells)]\n",
    "\n",
    "# (15) Run simulation for 100 ms\n",
    "sim.run(100)\n",
    "print(\"Simulation finished\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "JdMAKtkQEu-J"
   },
   "source": [
    "## The results\n",
    "\n",
    "Step **(16)** prints the timestamps of the spikes.\n",
    "\n",
    "Step **(17)** generates a plot of the sampling data.\n",
    "`arbor.simulation.samples` takes a ``handle`` of the probe we wish to examine. It returns a list\n",
    "of ``(data, meta)`` terms: ``data`` being the time and value series of the probed quantity; and\n",
    "``meta`` being the location of the probe. The size of the returned list depends on the number of\n",
    "discrete locations pointed to by the handle, which in this case is 1, so we can take the first element.\n",
    "(Recall that in step **(10)** we attached a probe to the ``\"root\"``, which describes one location.\n",
    "It could have described a `locset`.)\n",
    "\n",
    "Since we have created ``ncells`` cells, we have ``ncells`` traces. We should be seeing phase shifted traces, as the action potential propagated through the network.\n",
    "\n",
    "We plot the results using pandas and seaborn:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 646
    },
    "id": "HqpJ01fsEuf6",
    "outputId": "8f3c03af-8807-4c9a-8210-33cde75e30aa"
   },
   "outputs": [],
   "source": [
    "# (16) Print spike times\n",
    "print(\"spikes:\")\n",
    "for sp in sim.spikes():\n",
    "    print(\" \", sp)\n",
    "\n",
    "# (17) Plot the recorded voltages over time.\n",
    "print(\"Plotting results ...\")\n",
    "df_list = []\n",
    "for gid in range(ncells):\n",
    "    samples, meta = sim.samples(handles[gid])[0]\n",
    "    df_list.append(\n",
    "        pandas.DataFrame(\n",
    "            {\"t/ms\": samples[:, 0], \"U/mV\": samples[:, 1], \"Cell\": f\"cell {gid}\"}\n",
    "        )\n",
    "    )\n",
    "\n",
    "df = pandas.concat(df_list, ignore_index=True)\n",
    "seaborn.relplot(data=df, kind=\"line\", x=\"t/ms\", y=\"U/mV\", hue=\"Cell\", ci=None)"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "network_ring.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}