5 - simulator.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bff6608b-d5a8-4b7b-9b93-4587910698c9",
   "metadata": {
    "id": "4qlKKIdjfjPU",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# My First Neuron Simulator 𖡬\n",
    "\n",
    "![Johannes Blenke](https://images.unsplash.com/photo-1662242742653-13f0f76de431?ixlib=rb-4.0.3&ixid=MnwxMjA3fDB8MHxwaG90by1wYWdlfHx8fGVufDB8fHx8&auto=format&fit=crop&w=1740&q=80)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b45d1add",
   "metadata": {
    "id": "Wjt1ql5YVgts",
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Like how libraries provide packaged functionality, we can create a **Class** that holds _attributes_ (or _properties_) and _methods_ (and _functions_).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c0e92f1",
   "metadata": {
    "id": "Wjt1ql5YVgts",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "This helps **abstract** a lot of detail by **ENCAPSULATING** _implementation_ details and allowing us to re-use blocks of related code easily (i.e. collated pieces of logic forming a concept / object).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6146e61c",
   "metadata": {
    "id": "Wjt1ql5YVgts",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "You should aim to abstract away as much detail as necessary, but not more.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e38e925",
   "metadata": {
    "id": "Wjt1ql5YVgts",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "> The process of re-visiting old code and changing it (_hopefully an improvement!_) is known as **refactorisation** (noun) or **refactoring** (verb).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7358d7dc",
   "metadata": {
    "id": "Wjt1ql5YVgts",
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "You can define an object by creating a **CLASS** definition using the syntax\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb3fc2be",
   "metadata": {
    "id": "Wjt1ql5YVgts",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "```\n",
    "class <MyClassName>:\n",
    "  \"\"\"some comment about the class\"\"\"\n",
    "\n",
    "  def __init__(self):\n",
    "    <things to do when object is created>\n",
    "\n",
    "  def <other class method>(self, a):\n",
    "    \"\"\"a short comment about the method\"\"\"\n",
    "    self.a  = a # this class now has property a!\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e311f57d-52e3-411f-a40d-8dd3e0237a5a",
   "metadata": {
    "id": "Wjt1ql5YVgts",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "> Note how the first parameter for a class method must be `self`. This provides the method with access to _other_ class methods and properties. It's also a way to store arguments provided to methods as object properties.\n",
    "\n",
    "> `__init__` is a special method name (hence it starts _and_ ends with `__`). It is called when a class is _initialised_ or created (also known as **instantiation**) using `obj = MyClassName()`). `obj` is the **instantiated OBJECT** of the **class**. Hence, \"Object Orientated Programming\", or OOP.\n",
    "\n",
    "> Classes define the attributes and methods, Objects are what you play around with when you run your program.\n",
    "\n",
    "> There are a number of other special methods that classes have that can be overwritten!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f7ce666-617b-4a28-bf0d-bbe87d6a67da",
   "metadata": {
    "id": "UwheocVbF5UP",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Encapsulate the neuron logic so far in a single entity 🎁\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c160b39d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# download the synapses file so we can import it in the next cell\n",
    "!wget https://github.com/ChrisCurrin/imbizo_2024_problem_solving/raw/main/synapses.py synapses.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a57a9ec-e0e7-40a7-b894-064f452026c5",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from synapses import I_AMPA, I_GABA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44297228-f9f0-4c04-a20d-3257c74e13ca",
   "metadata": {
    "id": "tHZ-tL6Gfh35",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# create 'constants' with uppercase variable names\n",
    "DEFAULT_V_REST = -65  # resting membrane potential (mV)\n",
    "DEFAULT_C_M = 1  # membrane capacitance (nF)\n",
    "DEFAULT_R_M = 100  # membrane resistance (MOhm)\n",
    "\n",
    "\n",
    "# define a 'point' neuron's core properties\n",
    "class MyFirstNeuron:\n",
    "    \"\"\"A point neuron that can receive spike times for GABAA and AMPA synapses,\n",
    "    as well as constant external current. Note this only models subthreshold\n",
    "    membrane potentials (no action potentials are generated).\n",
    "\n",
    "    Note: although a point neuron, the length and radius of the neuron must be\n",
    "    provided (in um).\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(\n",
    "        self,\n",
    "        length,\n",
    "        radius,\n",
    "        name=\"my first neuron\",\n",
    "        R_m=DEFAULT_R_M,\n",
    "        C_m=DEFAULT_C_M,\n",
    "        V_0=DEFAULT_V_REST,\n",
    "    ):\n",
    "        print(f\"creating neuron {name}\")\n",
    "        # initialise variables\n",
    "        self.name = name  # name now belongs, and is a property of \"MyFirstNeuron\"\n",
    "        self.length = length * 1e-5  # um to dm\n",
    "        self.r = radius * 1e-5  # um to dm\n",
    "        self.R_m = R_m\n",
    "        self.C_m = C_m\n",
    "        self.V_rest = V_0\n",
    "        self.I_e = 0\n",
    "\n",
    "        # surface area of electrode placement (cm2)\n",
    "        self.A = (np.pi * self.length * self.r) * 100\n",
    "        # volume of sphere (dm3)\n",
    "        self.volume = 4 / 3 * np.pi * (self.length / 2) * self.r**2\n",
    "        # surface area of open cylinder (cm2)\n",
    "        self.SA = (self.length * np.pi * self.r + 2 * np.pi * self.r**2) * 100\n",
    "\n",
    "        # create lists, even if they start empty\n",
    "        self.V_t = []  # we are going to record the voltage\n",
    "        self.inh_spks = []  # list for spikes that arrive at GABA synapses\n",
    "        self.exc_spks = []  # list for spikes that arrive at AMPA synapses\n",
    "\n",
    "        # initialise variables for a simulation\n",
    "        self.reinit()\n",
    "\n",
    "    def set_inh_input(self, inh_spks):\n",
    "        self.inh_spks = inh_spks\n",
    "\n",
    "    def set_exc_input(self, exc_spks):\n",
    "        self.exc_spks = exc_spks\n",
    "\n",
    "    def set_external_input(self, I_e):\n",
    "        # current is in μA\n",
    "        self.I_e = I_e\n",
    "\n",
    "    def _I_syn(self, t, V_t):\n",
    "        \"\"\"Returns the sum of all synaptic currents.\n",
    "        Note that the method starts with an underscore \"_\" to indicate that\n",
    "        it should be for internal (to the class) use only!\n",
    "        \"\"\"\n",
    "        return I_GABA(t, V_t, self.inh_spks) + I_AMPA(t, V_t, self.exc_spks)\n",
    "\n",
    "    def step(self, t, dt):\n",
    "        V_t = self.V_t[-1]  # get last value\n",
    "        I_leak = (V_t - self.V_rest) / self.R_m\n",
    "        dV = (-I_leak - self._I_syn(t, V_t) + self.I_e * 1e-5 / self.A) * dt / self.C_m\n",
    "        self.V_t.append(V_t + dV)\n",
    "\n",
    "    def reinit(self, V_0=None):\n",
    "        \"\"\"Re-initialise the membrane potential list, starting at V_0 (or V_rest\n",
    "        if no argument is provided).\n",
    "        \"\"\"\n",
    "        if V_0 is None:  # None is a special value that provides a good placeholder\n",
    "            V_0 = self.V_rest  # use original value\n",
    "        self.V_t = [V_0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80ba1efb-9a5d-4ffe-bd03-0a1b062eda9d",
   "metadata": {
    "id": "rnwdjN-GEVKu",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "> Note that methods that start with a single underscore _discourage_ use outside the class. E.g. `_I_syn`\n",
    "\n",
    "> Double underscore makes the method private (difficult to 'see' outside the class and child classes).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35be2b61-5107-405a-be8d-09b4224f16a7",
   "metadata": {
    "id": "CFL3IYqYF-ff",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Add some logic that allows us to 'step through' time. 🚶\n",
    "\n",
    "Objects in the simulator move through time through their `step` method.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e95898e-ad3b-4ae7-9091-d4e6a017a279",
   "metadata": {
    "code_folding": [],
    "id": "oh5cDnIr_s6t",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "class Simulator:\n",
    "    \"\"\"A Simulator can be used to step through time by calling each of its\n",
    "    encapsulated neuron's ``MyFirstNeuron.step`` method.\n",
    "    The duration and the time step can be chosen when call ``run`` or\n",
    "    ``continuefor`` (which does does initialise objects)\n",
    "    The time is recorded in the ``t_list`` variable for retrieval later.\n",
    "\n",
    "    Note that only one Simulator ever exists.\n",
    "    We call this a Singleton class because efforts to construct a new object\n",
    "    instead retrieve the existing object and so there only exists a single object.\n",
    "    Here, this is implemented here by\n",
    "    - overriding the __new__ method (takes in `cls` - shorthand for 'class' -\n",
    "    instead of `self`)\n",
    "    - having a _sim class variable (it is at the top of the class definition)\n",
    "    other approaches can be found\n",
    "    https://python-3-patterns-idioms-test.readthedocs.io/en/latest/Singleton.html\n",
    "    \"\"\"\n",
    "\n",
    "    _sim = None\n",
    "\n",
    "    def __new__(cls):\n",
    "        if not cls._sim:\n",
    "            print(\"new Simulator\")\n",
    "            cls._sim = super(Simulator, cls).__new__(cls)\n",
    "            cls._sim.neurons = []\n",
    "        else:\n",
    "            print(\"existing Simulator\")\n",
    "        return cls._sim\n",
    "\n",
    "    def __init__(self):\n",
    "        self.init()\n",
    "\n",
    "    def init(self):\n",
    "        self.t_list = [0]\n",
    "        # reinitialise all neurons (for repeated simulations)\n",
    "        for neuron in self.neurons:\n",
    "            neuron.reinit()\n",
    "\n",
    "    def add_neuron(self, neuron):\n",
    "        # ensure only neurons with unique names are in the list\n",
    "        for nrn in self.neurons:\n",
    "            if neuron.name == nrn.name:  # found a previous version\n",
    "                sim.neurons.remove(nrn)  # remove from simulator\n",
    "                del nrn  # delete object\n",
    "        self.neurons.append(neuron)\n",
    "\n",
    "    def _run(self, duration, time_step):\n",
    "        \"\"\"private method to run the simulation from t_list[-1] until duration.\"\"\"\n",
    "        print()\n",
    "        start = self.t_list[-1]\n",
    "        stop = self.t_list[-1] + duration + time_step\n",
    "        for t in np.arange(start, stop, time_step):\n",
    "            self.t_list.append(t)\n",
    "            for neuron in self.neurons:\n",
    "                neuron.step(t, time_step)\n",
    "            # if i%log_every==0:\n",
    "            #   print(\"#\"*int(100*i/num_iter) + f\"> {t} ms\", end='\\r')\n",
    "        print()\n",
    "\n",
    "    def run(self, duration=100, time_step=0.1):\n",
    "        \"\"\"Run a simulation over time, calling the step method of each neuron in the\n",
    "        neurons list.\n",
    "\n",
    "        Values at t_0 are defined in all contained objects' reinit() methods\n",
    "        The first time step is therefore t_1 = dt.\n",
    "        \"\"\"\n",
    "        num_iter = int(duration / time_step)\n",
    "        print(\n",
    "            f\"running a simulation for {duration} ms \"\n",
    "            f\"using a dt of {time_step} ms.\"\n",
    "            f\"\\nThere will be {num_iter} iterations for each of \"\n",
    "            f\"{len(self.neurons)} neurons\"\n",
    "        )\n",
    "        self.init()\n",
    "        self._run(duration, time_step)\n",
    "        print(\"done!\")\n",
    "\n",
    "    def continuefor(self, duration=100, time_step=0.1):\n",
    "        num_iter = int(duration / time_step)\n",
    "        print(\n",
    "            f\"continuing a simulation from {self.t_list[-1]} ms \"\n",
    "            f\"for {duration} ms \"\n",
    "            f\"using a dt of {time_step} ms.\"\n",
    "            f\"\\nThere will be {num_iter} iterations for each of \"\n",
    "            f\"{len(self.neurons)} neurons\"\n",
    "        )\n",
    "        self._run(duration, time_step)\n",
    "        print(\"done!\")\n",
    "\n",
    "    # define a class method to reset the simulator\n",
    "    def clean(self):\n",
    "        while len(self.neurons) > 0:\n",
    "            neuron = self.neurons.pop()\n",
    "            del neuron\n",
    "        self.init()\n",
    "        return self\n",
    "\n",
    "    # create an alias for the clean method\n",
    "    reset = clean"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6df3b435",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Simulate a Neuron! 🎉\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "efc635c7-2ae9-43ac-8d8b-608a1b822043",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 386
    },
    "id": "s3f96XNg_3PF",
    "outputId": "ae6f8d53-c8a4-411c-f343-8c0b1ae0ac39",
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "T = 1000  # duration (ms)\n",
    "dt = 0.1  # time step (ms)\n",
    "\n",
    "# create a simulator\n",
    "sim = Simulator()\n",
    "\n",
    "# create a neuron with input\n",
    "nrn = MyFirstNeuron(23, 11, name=\"my 1st neuron!\")\n",
    "nrn.set_inh_input(np.random.randint(0, T, size=10))\n",
    "nrn.set_exc_input(np.random.randint(0, T, size=10))\n",
    "\n",
    "# TODO: create ANOTHER neuron with external input\n",
    "# ...\n",
    "\n",
    "# adding all neurons we want to update\n",
    "sim.add_neuron(nrn)\n",
    "\n",
    "# run a simulation of a neuron\n",
    "sim.run(T, dt)\n",
    "\n",
    "# visualise\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# create a figure with an axis for plotting\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "# plot voltage over iterations\n",
    "ax.plot(sim.t_list, nrn.V_t)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "319ea83b-1462-4ebd-9b09-5afb5e6ad578",
   "metadata": {
    "id": "FW13rGMfr0PQ",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## **Tasks** ✏️\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "194e2a98-3b3c-4f29-97a2-57fa8346e280",
   "metadata": {
    "id": "UHawHpxYxJKL",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Task 1:** Run the cell below with different time steps.\n",
    "\n",
    "- How does this affect the voltage trace?\n",
    "\n",
    "**Task 2:** Create another neuron with additional external input and compare the traces\n",
    "\n",
    "- So the same excitatory and inhibitory input as \"my 1st neuron!\" but **also** with external input.\n",
    "- We will use this neuron in future traces\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6a695cf-679f-4722-8725-ec03287fef9d",
   "metadata": {
    "colab": {
     "background_save": true,
     "base_uri": "https://localhost:8080/",
     "height": 434
    },
    "id": "pIqdvVssCnIG",
    "outputId": "27232229-c451-422f-a526-a590143eeab1",
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# @title Task Solution Space { display-mode: \"form\" }\n",
    "T = 1000  # duration (ms)\n",
    "\n",
    "### --------------------------------------------\n",
    "### TASK 1: change dt\n",
    "# TODO: edit dt below\n",
    "dt = 0.1  # time step (ms)\n",
    "### --------------------------------------------\n",
    "\n",
    "# create a simulator\n",
    "sim = Simulator()\n",
    "\n",
    "# create a neuron with input\n",
    "nrn = MyFirstNeuron(23, 11, \"my 1st neuron!\")\n",
    "exc_input = np.random.randint(0, T, size=10)\n",
    "inh_input = np.random.randint(0, T, size=10)\n",
    "\n",
    "nrn.set_inh_input(inh_input)\n",
    "nrn.set_exc_input(exc_input)\n",
    "\n",
    "### --------------------------------------------\n",
    "### TASK 2: create ANOTHER neuron with external input\n",
    "nrn_with_ext_input = MyFirstNeuron(23, 11, \"my 1st neuron with external input\")\n",
    "# TODO: EDIT BELOW >>>\n",
    "\n",
    "### <<< EDIT ABOVE\n",
    "### --------------------------------------------\n",
    "\n",
    "# adding all neurons we want to update\n",
    "sim.add_neuron(nrn)\n",
    "sim.add_neuron(nrn_with_ext_input)\n",
    "\n",
    "# run a simulation of a neuron\n",
    "sim.run(T, dt)\n",
    "\n",
    "# visualise\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# create a figure with an axis for plotting\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "# plot voltage over iterations\n",
    "names = []\n",
    "for n in sim.neurons:\n",
    "    ax.plot(sim.t_list, n.V_t)\n",
    "    names.append(n.name)\n",
    "\n",
    "# this is a *list comprehension*\n",
    "names = [n.name for n in sim.neurons]\n",
    "ax.legend(names, frameon=False)\n",
    "\n",
    "# add some axis labels\n",
    "ax.set_ylabel(\"Voltage [$V_m$] (mV)\")\n",
    "ax.set_xlabel(\"Time (ms)\")\n",
    "\n",
    "# don't forget a title\n",
    "ax.set_title(\"How does external input affect voltage?\")\n",
    "\n",
    "# make the borders prettier\n",
    "ax.spines[\"right\"].set_visible(False)\n",
    "ax.spines[\"top\"].set_visible(False)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e03f9bd-48f5-45fb-9d4e-77334e80574e",
   "metadata": {
    "id": "8zUXuf54lDWY",
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Task 3**\n",
    "\n",
    "Instead of the **analytical** solution (can compute any time `t`) for `g_syn`, implement the **numerical** solution (need to know `t-1` to calculate `t`).\n",
    "\n",
    "$\\begin{aligned}\n",
    "\\frac{dg_{syn}}{dt} &= -\\frac{g_{syn}}{\\tau_{syn}} \\\\\n",
    "g_{syn} &= g_{syn} + g_{syn_{\\rm{max}}}   (\\textrm{when } t = t_{s})\n",
    "\\end{aligned}$\n",
    "\n",
    "Further questions\n",
    "\n",
    "1. What are the advantages + disadvantages between the 2 solutions?\n",
    "1. Derive the numerical solution from the analytical solution (ignore the $\\sum$)\n",
    "1. Integrate the numerical solution to get the analytical solution\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2e63aa5",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# @title Task 3 Scaffold\n",
    "\n",
    "\n",
    "class Synapse:\n",
    "    \"\"\"Implemented synaptic conductance using Object Orientated Programming and\n",
    "    numerical solution to g_syn\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, spks, g_max=50, tau=5, E=0):\n",
    "        pass\n",
    "\n",
    "    def step(self, t, dt, V_t):\n",
    "        \"\"\"\n",
    "        Parameters:\n",
    "            t (float): The time point at which to calculate the current (ms).\n",
    "            dt (float): The time step (ms).\n",
    "            V_t (float): The membrane potential at time t (mV).\n",
    "\n",
    "        Returns:\n",
    "            float: The current generated by the synapse at time t (mA).\n",
    "\n",
    "        Algorithm:\n",
    "            1. decay any conductance\n",
    "            2. add conductance for incoming spikes\n",
    "            3. set synaptic current\n",
    "        \"\"\"\n",
    "        # 1. decay any conductance\n",
    "\n",
    "        # 2. add conductance for incoming spikes\n",
    "\n",
    "        # 3. set synaptic current\n",
    "\n",
    "        self.I_syn = 0\n",
    "\n",
    "\n",
    "# Usage\n",
    "\n",
    "\n",
    "class NeuronWithNumericalSynapse(MyFirstNeuron):\n",
    "    def __init__(self, *args, name=\"numerical synapses\", **kwargs):\n",
    "        super().__init__(*args, name=name, **kwargs)\n",
    "\n",
    "    def set_inh_input(self, inh_input, g_max=50, tau=8, E=-70):\n",
    "        self.inh_synapse = Synapse(inh_input, g_max=g_max, tau=tau, E=E)  # <-- changed\n",
    "\n",
    "    def set_exc_input(self, exc_spks, g_max=50, tau=4, E=0):\n",
    "        self.exc_synapse = Synapse(exc_spks, g_max=g_max, tau=tau, E=E)  # <-- changed\n",
    "\n",
    "    def _I_syn(self, t, dt, V_t):  # <-- method changed (call in `step` changed)\n",
    "        self.inh_synapse.step(t, dt, V_t)\n",
    "        self.exc_synapse.step(t, dt, V_t)\n",
    "        return self.inh_synapse.I_syn + self.exc_synapse.I_syn"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0efb33ab-792d-459b-a2c7-5b8672fab3fe",
   "metadata": {
    "id": "WhwrRDnw-PBy",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The Leaky Integrate-and-Fire Neuron 🔥\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a1bc750-2ad1-4ebd-a9ef-a4b55d286e01",
   "metadata": {
    "id": "IXpTUVkT-Ors",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Thus far, we have implemented the change in voltage given inputs, but there's no concept of an action potential going on!\n",
    "\n",
    "The simplest approach is to define a _threshold_ $V_{thresh}$ for the neuron that once reached at time $t_i$\n",
    "\n",
    "1. $V_m$ goes to `0` (or `20`) at $t_{i+1}$\n",
    "1. After hyperpolarises to $V_{reset}$ at $t_{i+2}$\n",
    "\n",
    "Because we want to add extra functionality, we'll use this opportunity to introduce another 2 (of the 3) tenets of programming. **Inheritance** and **Polymorphism** (the other tenet was **encapsulation**).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc514d5a-9121-4dec-9248-208a582936d6",
   "metadata": {
    "id": "TMTZyPVwKhhS",
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Inheritance 🧬\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb242c19-c28d-4d0c-85bf-a518f1c8972f",
   "metadata": {
    "id": "Lch7V1p0MdCE",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We already have a bunch of logic in `MyFirstNeuron` that has the same applicability to a leaky integrate-and-fire (LIF) neuron (setting initial values in `__init__`, adding inputs, summing synaptic inputs, and re-initialising the neuron's voltage).\n",
    "\n",
    "We can _inherit_ this logic using\n",
    "\n",
    "```\n",
    "class <NewClass>(<ClassToInheritFrom>):\n",
    "```\n",
    "\n",
    "> note that previously we _inherited from object_.\n",
    "\n",
    "All methods that are in `<ClassToInheritFrom>` (e.g. `MyFirstNeuron`) are accessible in this new class!\n",
    "\n",
    "```\n",
    "class LIFNeuron(MyFirstNeuron):\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a09ff5be-1595-4b3f-9af3-aaf9409cfaca",
   "metadata": {
    "id": "RaniDxe-Qnif",
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "But what if we want to _change_ a method.\n",
    "\n",
    "For example, the `step` method in `MyFirstNeuron` doesn't have the logic to handle our (artificial) action potentials -the \"fire\" 🔥 in \"integrate-and-fire\". It only has the sub-threshold \"integrate\" part.\n",
    "\n",
    "> Q: Where does the \"leaky\" come from?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f096800-566d-49f6-83cb-0fd4be962843",
   "metadata": {
    "id": "i5gzV5i0RdKr",
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Polymorphism 𐩯\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11e695b0",
   "metadata": {
    "id": "QvDPEhL3QsII",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "To change this, we create a `step` method in this new class with the new logic we want. Using the same name to implement new functionality is known as **polymorphism**! Python determines which `step` method you want by the object that calls it (i.e. is it a `MyFirstNeuron` or `LIFNeuron`).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c48bb42e-19f9-4a62-82ba-c02a12313a19",
   "metadata": {
    "id": "QvDPEhL3QsII",
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In fact, the way we have implemented things, we do need to explicitly _repeat_ some logic. We could go back to `MyFirstNeuron` and **refactor** `step` to make things smoother by re-using the same logic. But we're only repeating logic once so we'll accept this for now. If this happens again, go back and optimise to reduce sources of error and the number of places you would need to alter something.\n",
    "\n",
    "```\n",
    "def step(self, t, dt):\n",
    "  V_t = self.V_t[-1] # get last value\n",
    "  if V_t == 0:\n",
    "    # neuron has spiked at previous time step\n",
    "    V_t = self.V_rest\n",
    "  elif V_t >= V_thresh: # what if you change the order of the if statements?\n",
    "    # voltage has reached threshold, spike!\n",
    "    V_t = 0\n",
    "  else:\n",
    "    # this is repeated from MyFirstNeuron\n",
    "    I_leak = (V_t - self.V_rest)/self.R_m\n",
    "    dV = (-I_leak - self.__I_syn(t, V_t) + self.I_e/self.A) * dt/self.C_m\n",
    "    V_t = V_t + dV\n",
    "  self.V_t.append(V_t)\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7749aa0c-7685-4b9c-a359-e608d4187b48",
   "metadata": {
    "id": "JicSGsIJRjyy",
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Polymorphism with Inheritance 𐩯xx\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a95220db",
   "metadata": {
    "id": "LFWDxqUNRVcT",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We also have new parameters, `V_thresh` and `V_reset`, we want to introduce that only needs to be defined once when we create (**instantiate**) our neuron.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c2d15f9",
   "metadata": {
    "id": "LFWDxqUNRVcT",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We overwrite `__init__` by creating a method with the same name, but then we **call** the **parent** (`MyFirstNeuron`) using the super special `super()` function, along with the method we want: `super().__init__()`.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d58676a6",
   "metadata": {
    "id": "LFWDxqUNRVcT",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Because our parent method has parameters it accepts, we need to pass these along the hierarchy. We take in arbitrary _arguments_ (e.g. `radius`) and _keyword arguments_ (e.g. `R_m=1`) that other methods want with `*args` and `**kwargs`, represpectively.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3d2bb11-7ddf-41fc-9301-2b6804506353",
   "metadata": {
    "id": "LFWDxqUNRVcT",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "> We can explicitly accept arguments and pass them on too, e.g. maybe we want a different default value!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36c2355c-037b-4535-bd0a-aab9d6bfd4de",
   "metadata": {
    "id": "EIx1XJlv5OId",
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Another thing: [Logging 🪵](https://docs.python.org/3/library/logging.html)\n",
    "\n",
    "We have been using `print` to display text, but a better (more informative and robust) way is using the `logging` module.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e9f9571-42b5-4f6e-a2d0-6180707d7ed4",
   "metadata": {
    "colab": {
     "background_save": true,
     "base_uri": "https://localhost:8080/"
    },
    "id": "XLLDq2xh5ydZ",
    "outputId": "45481d68-44cf-489f-ef07-ccd5676ea8e0",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import logging\n",
    "\n",
    "logging.basicConfig(\n",
    "    format=\"%(asctime)s %(name)-15s %(levelname)-8s %(message)s\",\n",
    "    datefmt=\"%H:%M:%S\",\n",
    "    level=logging.INFO,\n",
    ")\n",
    "\n",
    "logger = logging.getLogger(\"MyLogger\")\n",
    "\n",
    "logger.setLevel(logging.WARNING)\n",
    "\n",
    "logger.debug(\"Only show when debugging\")\n",
    "logger.info(\"replace print statements with this instead\")\n",
    "logger.warning(\"this displays a warning\")\n",
    "logger.error(\"there was an error\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34b7c6d2-426b-43be-acb8-7e021147af3e",
   "metadata": {
    "id": "y827oGSuRgCw",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Implementation of Leaky Integrate-and-Fire Neuron\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a50bf43d-5524-4843-ad32-565fa50e21e9",
   "metadata": {
    "colab": {
     "background_save": true
    },
    "id": "lDh9lxxn98ph",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "class LIFNeuron(MyFirstNeuron):\n",
    "    \"\"\"\n",
    "    Leaky integrate-and-fire neuron that inherits from MyFirstNeuron.\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(\n",
    "        self,\n",
    "        *args,\n",
    "        name=\"LIF\",  # replace default\n",
    "        V_thresh=-55,\n",
    "        V_reset=-68,  # new arguments\n",
    "        **kwargs,\n",
    "    ):  # other keyword arguments\n",
    "        logging.debug(f\"other arguments passed = {args}\")\n",
    "        logging.debug(f\"other keyword arguments passed = {kwargs}\")\n",
    "        super().__init__(*args, name=name, **kwargs)\n",
    "        self.V_thresh = V_thresh\n",
    "        self.V_reset = V_reset\n",
    "\n",
    "    def step(self, t, dt):\n",
    "        V_t = self.V_t[-1]  # get last value\n",
    "        if V_t == 0:\n",
    "            # neuron has spiked at previous time step\n",
    "            V_t = self.V_reset\n",
    "        elif V_t >= self.V_thresh:  # what if you change the order of the if statements?\n",
    "            # voltage has reached threshold, spike!\n",
    "            V_t = 0\n",
    "        else:\n",
    "            # this is repeated from MyFirstNeuron,\n",
    "            # but we could also call super().step(t, dt) with some consideration\n",
    "            # (the parent step method also appends to V_t list)\n",
    "            I_leak = (V_t - self.V_rest) / self.R_m\n",
    "            dV = (\n",
    "                (-I_leak - self._I_syn(t, V_t) + self.I_e * 1e-5 / self.A)\n",
    "                * dt\n",
    "                / self.C_m\n",
    "            )\n",
    "            V_t = V_t + dV\n",
    "        self.V_t.append(V_t)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80517a4a-f679-4338-8816-a17e02eb78a0",
   "metadata": {
    "id": "Fg10qLgVFVLY",
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Create a shortcut for visualisations\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84b46671-bb8c-4ead-aeda-59b2e4abb499",
   "metadata": {
    "colab": {
     "background_save": true
    },
    "id": "WlSbtab1FXdU",
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# visualise\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.cbook import flatten\n",
    "\n",
    "\n",
    "def plot_v(sim, attrs=None, title=\"\", **kwargs):\n",
    "    # note that list parameters should have defaults as None in the method and\n",
    "    #   defined in the body, even if it's empty\n",
    "    if attrs is None:\n",
    "        attrs = []\n",
    "\n",
    "    # create a figure with an axis for plotting\n",
    "    fig, axs = plt.subplots(\n",
    "        nrows=1 + len(attrs),\n",
    "        squeeze=False,\n",
    "        sharex=True,\n",
    "        gridspec_kw={\"height_ratios\": [5] + [1] * len(attrs)},\n",
    "    )\n",
    "    ax = axs[0, 0]\n",
    "\n",
    "    # plot voltage over iterations\n",
    "    for neuron in sim.neurons:\n",
    "        ax.plot(sim.t_list, neuron.V_t, **kwargs)\n",
    "        for i, attr in enumerate(attrs):\n",
    "            try:\n",
    "                axs[i + 1, 0].plot(sim.t_list, getattr(neuron, attr))\n",
    "            except AttributeError:\n",
    "                logger.debug(\n",
    "                    f\"attribute {attr} did not exist in {neuron.name}. Continuing...\"\n",
    "                )\n",
    "                logger.debug(\n",
    "                    f\"\"\"do this next time `val = getattr(neuron, attr, None)`\n",
    "                              [see `help(getattr)`]\"\"\"\n",
    "                )\n",
    "\n",
    "    names = [n.name for n in sim.neurons]\n",
    "    ax.legend(names, loc=\"upper left\", bbox_to_anchor=(1, 1), frameon=False)\n",
    "\n",
    "    # add some axis labels\n",
    "    ax.set_ylabel(\"$V_m$ (mV)\")  # <-- math symbols in-between $ $\n",
    "    axs[-1, -1].set_xlabel(\"Time (ms)\")\n",
    "\n",
    "    # don't forget a title\n",
    "    ax.set_title(title)\n",
    "\n",
    "    # make the borders prettier\n",
    "    for ax in flatten(axs):\n",
    "        ax.spines[\"right\"].set_visible(False)\n",
    "        ax.spines[\"top\"].set_visible(False)\n",
    "\n",
    "    return fig, axs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b72346ba-fdef-4bb2-a2c9-4f7e3e5d6940",
   "metadata": {
    "colab": {
     "background_save": true,
     "base_uri": "https://localhost:8080/",
     "height": 434
    },
    "id": "BGuB-c_JJIoK",
    "outputId": "1899840f-863d-48cd-a717-939263b91749",
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# create a simulator\n",
    "sim = Simulator()\n",
    "sim.clean()\n",
    "\n",
    "T = 1000  # duration (ms)\n",
    "dt = 0.1  # time step (ms)\n",
    "\n",
    "# create a neuron with input\n",
    "lif_nrn = LIFNeuron(23, 11)\n",
    "lif_nrn.set_inh_input(np.random.randint(0, T, size=10))\n",
    "lif_nrn.set_exc_input(np.random.randint(0, T, size=10))\n",
    "\n",
    "sim.add_neuron(lif_nrn)\n",
    "\n",
    "# run a simulation of a neuron\n",
    "sim.run(T, dt)\n",
    "\n",
    "plot_v(sim, title=\"Let there be spikes!\")\n",
    "print(f\"there were {np.sum(np.array(lif_nrn.V_t)>=0)} spikes\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f462d2f-f69b-4913-adfa-8d89d1e39fdb",
   "metadata": {
    "id": "9eGzo8DP5VOz",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### How does _external input current_ affect the _firing rate_ of a neuron?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e3cb425-a8f2-492f-b28e-b9c64ba09d1d",
   "metadata": {
    "colab": {
     "background_save": true,
     "base_uri": "https://localhost:8080/"
    },
    "id": "iwIaHGNa7LNp",
    "outputId": "2c7dce96-4835-405e-e023-4f42b2466300",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "sim = Simulator()\n",
    "# get rid of previous neurons\n",
    "sim.clean()\n",
    "\n",
    "# create neurons each with their own external current value\n",
    "for current in np.arange(0, 0.5, 0.05):\n",
    "    nrn = LIFNeuron(23, 11, name=f\"{current:.2f}\")\n",
    "    nrn.set_external_input(current)\n",
    "    sim.add_neuron(nrn)\n",
    "sim.run(T, dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e950f29-4bf5-4756-88b7-293dd44cba1c",
   "metadata": {
    "cellView": "form",
    "colab": {
     "background_save": true,
     "base_uri": "https://localhost:8080/",
     "height": 295
    },
    "id": "YcvzeRzcEsOU",
    "outputId": "3988f5de-5620-4ada-edf8-e45707b9e294",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# @title Plot\n",
    "fig, ax = plot_v(sim, title=\"Let there be spikes!\", alpha=0.5, lw=1)\n",
    "ax[0, 0].get_legend().set_title(\"External\\nCurrent (μA)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "075eac49-68ab-4044-a8a9-ba9051acb094",
   "metadata": {
    "colab": {
     "background_save": true
    },
    "id": "F_W775u4CXQp",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# convert Voltage for each neuron to an array\n",
    "v = np.empty(shape=(len(sim.t_list), len(sim.neurons)))\n",
    "for i, nrn in enumerate(sim.neurons):\n",
    "    v[:, i] = nrn.V_t\n",
    "# calculate the instantaneous firing rate\n",
    "ifr = np.sum(v >= 0, axis=0) / (sim.t_list[-1] / 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a139fa23-57e9-4bf6-ac7f-6b33690c335b",
   "metadata": {
    "colab": {
     "background_save": true,
     "base_uri": "https://localhost:8080/",
     "height": 312
    },
    "id": "Smtlt98vHl9S",
    "outputId": "e0e06126-6a39-4c24-8c98-30aec6e91008",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot([float(nrn.name) for nrn in sim.neurons], ifr)\n",
    "ax.set_xlabel(\"Current (μA)\")\n",
    "ax.set_ylabel(\"Firing rate (Hz)\")\n",
    "ax.set_title(\"I-F 'curve'\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "662b09f7-537d-4536-8794-61a28462d12b",
   "metadata": {
    "id": "06LikNcZLeNj",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### **Homework/Project 2**: How does _**excitatory** and **inhibitory** input_ affect the _firing rate_ of a neuron?\n",
    "\n",
    "Similar to above for external current, determine\n",
    "\n",
    "- how does excitatory drive affect firing rate?\n",
    "\n",
    "  - you should produce an input-output (I-O) curve where excitation is on the x-axis and firing rate is on the y-axis\n",
    "\n",
    "- how does inhibition affect this I-O relationship?\n",
    "\n",
    "- how does timing of the excitatory + inhibitory inputs affect the relationship?\n",
    "\n",
    "  - what if they both have the same input frequency (\"Balanced input\")\n",
    "  - what if they have the same frequency, but are offset?\n",
    "  - what does noise add or take away from this I-O function?\n",
    "\n",
    "- create a circuit with 2 excitatory neurons and 1 inhibitory neuron. Connect them to each other. \\\n",
    "  _hint when one fires, add a spike time to its connected pair_\n",
    "  - Increase the number of neurons (try keep a 80:20 percentage of excitatory: inhibitory neurons).\n",
    "  - Implement Hebbian learning\n"
   ]
  },
  {
   "cell_type": "raw",
   "id": "adcbf2fc",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}