Multi-compartmental_neuron_model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8a4cf380-de07-4e30-8fcb-f40663c3c98b",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "[![Open In Colab](./colab-badge.png)](https://colab.research.google.com/github/subhacom/moose-notebooks/blob/main/Multi-compartmental_neuron_model.ipynb) [![Binder](./binder_logo.png)](https://mybinder.org/v2/gh/subhacom/moose-notebooks/HEAD)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c500cd1b-7e78-4654-af6e-9e7f9bb55e4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Only required on colab!\n",
    "# !pip install pymoose"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bcd481a9-1241-4dc0-b01f-708c8927c067",
   "metadata": {},
   "source": [
    "# Neuronal cable with multiple compartments\n",
    "Electrically speaking, a single compartment represents a single point. A standalone compartment cannot accommodate voltage gradients along neurites. Real neurons have complex tree-like structure. To take the geometry of a neuron into account, we need to put together many compartments to imitate the structure (morphology) of the neuron."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d30722df-abcb-45b5-ba2f-f9ee73f888f2",
   "metadata": {},
   "source": [
    "## Two-compartment model\n",
    "Let us start with the simple case of a soma with an axon. We can model this with two compartments. Any current injected into the soma will raise its voltage from the resting membrane potential $Em$, but the axonal compartment will remain at $Em$ initially. Thus electric current will flow from the soma into the axon. This current will face resistance when flowing through the cytoplasm along the length of the neuron. This is the cytoplasmic resistance or axial resistance, and is represented by the field `Ra` of `Compartment` objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "483ac745-be6d-4add-977d-1f284b5229b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import moose"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b50d5d6e-7bca-4f4a-9cf2-5ba0d7d5d917",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = moose.Neutral('/data')\n",
    "model = moose.Neutral('/model')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0e88ad8-c914-419d-bb5a-2b1ee611e06d",
   "metadata": {},
   "outputs": [],
   "source": [
    "soma = moose.Compartment('/model/soma')\n",
    "axon = moose.Compartment('/model/axon')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef33f266-9958-41f9-8e3a-b55c4d411158",
   "metadata": {},
   "source": [
    "We have connect the soma to the axon somehow so the simulator can compute how much electric current is flowing from the soma into the axon (or vice versa). This requires that the voltage at the soma as well as the resistance between the soma and the axon be known when updating the axon. The computation at the soma too needs to know the voltage at the axon to figure out the current flowing out of the soma. This two way communication is achieved by a special kind of message called a shared message. \n",
    "If we want to use the soma's `Ra` for computing the current flowing from the soma into the axon, we must connect the `raxial` shared field of `soma` to `axial` dest field of `axon`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8ef257e-7497-49b3-8560-d1681d6f1484",
   "metadata": {},
   "outputs": [],
   "source": [
    "moose.connect(soma, 'raxial', axon, 'axial')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47e9783d-5a11-43af-bad2-7e6f10393bbd",
   "metadata": {},
   "source": [
    "We can connect a pulse generator to the soma in order to deliver a current injection from 10 ms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d55007f7-4a7a-44b1-8584-0224fefcc760",
   "metadata": {},
   "outputs": [],
   "source": [
    "pulse = moose.PulseGen('/model/pulse')\n",
    "moose.connect(pulse, 'output', soma, 'injectMsg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71585322-3761-4f1a-bc49-cab24258aee0",
   "metadata": {},
   "outputs": [],
   "source": [
    "pulse.firstDelay = 10e-3\n",
    "pulse.firstLevel = 10e-6\n",
    "pulse.firstWidth = 1e9  # hold a constant current injection"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1df1fe6-9abb-4106-9c26-75665c8b8140",
   "metadata": {},
   "source": [
    "We shall record the membrane potential at both the soma and the axon. Therefore two `Table`s must be created and connected to the two compartments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00d43823-95c1-4781-b490-cce571bc8bcc",
   "metadata": {},
   "outputs": [],
   "source": [
    "vm_soma = moose.Table('/data/vm_soma')\n",
    "vm_axon = moose.Table('/data/vm_axon')\n",
    "moose.connect(vm_soma, 'requestOut', soma, 'getVm')\n",
    "moose.connect(vm_axon, 'requestOut', axon, 'getVm')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6790ef63-566f-4249-9725-23a71925ab79",
   "metadata": {},
   "source": [
    "Now we can set the properties of the two compartments. For simplicity, let us keep the same values for both."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96cb0b76-ce68-4c00-9948-9fa59962e06f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set properties of the soma compartment\n",
    "soma.initVm = -70e-3\n",
    "soma.Em = -70e-3\n",
    "soma.Rm = 1e4\n",
    "soma.Cm = 1e-7\n",
    "\n",
    "# Set properties of the axonal compartment\n",
    "axon.initVm = -70e-3\n",
    "axon.Em = -70e-3\n",
    "axon.Rm = 1e4\n",
    "axon.Cm = 1e-7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b872a2a0-ea00-439e-9b93-b24b1989d1e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "runtime = 50e-3\n",
    "moose.reinit()\n",
    "moose.start(runtime)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6578b07-42e9-44a0-a0d9-c1d696987777",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(np.arange(len(vm_soma.vector)) * vm_soma.dt * 1e3, vm_soma.vector * 1e3, label='soma')\n",
    "plt.plot(np.arange(len(vm_axon.vector)) * vm_axon.dt * 1e3, vm_axon.vector * 1e3, '--', label='axon')\n",
    "plt.xlabel('Time (ms)')\n",
    "plt.ylabel('Voltage (mV)')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "777e85cb-beb9-43d1-8143-ce4f029e4906",
   "metadata": {},
   "source": [
    "The two membrane potentials seem to perfectly overlap! If we check the axial resistance of the soma, it becomes clear that the default `Ra` is too low compared to the membrane resistance. So almost all the current injected into the soma is flowing into the axon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2792c091-1582-4469-a255-73d013ece4df",
   "metadata": {},
   "outputs": [],
   "source": [
    "soma.Ra"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4199f48-fe42-4ded-b97e-f6dbdd5dc859",
   "metadata": {},
   "source": [
    "Let us bring the axial resistance to the same order as the membrane resistance, so the current gets split more evenly through the soma membrane (`Rm`) and the cytoplasm (`Ra`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7345d1d-75ea-47b9-9dbe-7cb102d2a9a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "soma.Ra = 1e4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a13f9a2d-f979-4f96-9bd9-cf6296ce4096",
   "metadata": {},
   "outputs": [],
   "source": [
    "runtime = 50e-3\n",
    "moose.reinit()\n",
    "moose.start(runtime)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ecaa21c-64b9-4f8e-9e5c-1d3a22e32d75",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(np.arange(len(vm_soma.vector)) * vm_soma.dt * 1e3, vm_soma.vector * 1e3, label='soma')\n",
    "plt.plot(np.arange(len(vm_axon.vector)) * vm_axon.dt * 1e3, vm_axon.vector * 1e3, '--', label='axon')\n",
    "plt.xlabel('Time (ms)')\n",
    "plt.ylabel('Voltage (mV)')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "669a85bb-991a-4ed9-ab36-002c10668eab",
   "metadata": {},
   "source": [
    "## Making a longer cable with many compartments\n",
    "We can investigate how a voltage deflection at the soma propagates through a neuronal cable by making a longer cable. We can put the repeating code inside a for loop. But before that, let us clean up the system. We can delete all the components we created so far and start on a clean slate:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "056c1204-f6d6-4b04-8709-d6a12ec2b5b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "moose.delete(data)\n",
    "moose.delete(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b023a95-6fef-4b9a-8c2c-67c97faddd46",
   "metadata": {},
   "source": [
    "The `delete` function deletes a moose object and all its children, along with any connections they have with other objects. This way we got rid of the earlier model and tables to record data. We re-create this structure below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80465efd-b959-4ecb-9613-de07996d043d",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = moose.Neutral('/data')\n",
    "model = moose.Neutral('/model')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12f62cfd-f414-4f30-b59a-73312c3b1dcd",
   "metadata": {},
   "outputs": [],
   "source": [
    "comps = []\n",
    "vm_tabs = []\n",
    "ncomp = 10\n",
    "for ii in range(ncomp):\n",
    "    comp = moose.Compartment(f'/model/comp_{ii}')\n",
    "    comp.initVm = -70e-3\n",
    "    comp.Em = -70e-3  \n",
    "    comp.Rm = 1e5   # 100 kOhm\n",
    "    comp.Cm = 1e-7  # 100 nF\n",
    "    comp.Ra = 1e4   # 10 kOhm\n",
    "    if len(comps) > 0:\n",
    "        moose.connect(comps[-1], 'raxial', comp, 'axial')\n",
    "    comps.append(comp)\n",
    "    # Recording\n",
    "    vm_tab = moose.Table(f'/data/vm_{ii}')\n",
    "    moose.connect(vm_tab, 'requestOut', comp, 'getVm')\n",
    "    vm_tabs.append(vm_tab)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4a6dbe0-fad4-48bf-9606-d77f0f8b5a4a",
   "metadata": {},
   "source": [
    "Create a pulse generator and connect it to the first compartment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b34fb390-7477-4afa-b87d-5bcbc4fdeeb6",
   "metadata": {},
   "outputs": [],
   "source": [
    "pulse = moose.PulseGen('/model/pulse')\n",
    "moose.connect(pulse, 'output', comps[0], 'injectMsg')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45707ae5-3595-489a-91d8-ccd70e231962",
   "metadata": {},
   "source": [
    "Set the properties of the pulse generator. We are using 10 ms delay and 1 uA current injection of 1 ms duration. To avoid repetition of the pulse during the simulation, we set the second delay (`delay[1]`) to a large value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1005eda2-836d-4a78-9f6e-af08a39f7dbb",
   "metadata": {},
   "outputs": [],
   "source": [
    "pulse.delay[0] = 10e-3\n",
    "pulse.level[0] = 1e-6\n",
    "pulse.width[0] = 1e-3\n",
    "pulse.delay[1] = 1e9   # very long time until repeating the pulse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6fe6b852-85ca-476e-aed0-e8f3db3f572c",
   "metadata": {},
   "outputs": [],
   "source": [
    "runtime = 50e-3\n",
    "moose.reinit()\n",
    "moose.start(runtime)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d2873c7-0420-4670-b88c-cd9ccb6495f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "for tab in vm_tabs:\n",
    "    plt.plot(np.arange(len(tab.vector)) * tab.dt * 1e3, tab.vector * 1e3, label=tab.name)\n",
    "plt.xlabel('Time (ms)')\n",
    "plt.ylabel('Voltage (mV)')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f441f75-06d4-4d90-b7c7-3798231e337d",
   "metadata": {},
   "source": [
    "We increased the values of `Rm` and `Ra` above so that the delay in signal propagation becomes clearly visible from the plots. Both the onset of depolarization and time to reach the peak voltage increase as we go further away from the current injection site (`vm_0` is the voltage of `comp_0`, where the pulse is injected).\n",
    "\n",
    "Also notice how the curves get smoother. This is lowpass filtering, where sharp changes in a signal get more attenuated and slower changes are passed with less attenuation. An RC circuit, as is our passive compartment model, acts as a lowpass filter."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c81a9b78-f46c-405e-81bd-fe180df19e5d",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "Modify the values of `Rm`, `Cm`, and `Ra` and try to explain the effects on signal propagation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd873f4d-9cd5-4e59-b74d-ef600d4237bf",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}