notebook/models/08-spiking_neural_systems.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "This is an introduction to spiking neural systems with Julia's DifferentialEquations package.\nWe will cover four different models: leaky integrate-and-fire, Izhikevich, adaptive exponential\n integrate-and-fire and Hodgkin-Huxley model. Let's get started with the leaky integrate-and-fire (LIF) model.\n\n## The Leaky-Integrate-and-Fire Model\nThe LIF model is an extension of the integrate-and-fire (IF) model. While the IF\nmodel simply integrates input until it fires, the LIF model integrates input but\nalso decays towards an equilibrium potential. This means that inputs that arrive\nin quick succession have a much higher chance to make the cell spike as opposed\nto inputs that are further apart in time. The LIF is a more realistic neuron\nmodel than the IF because it is known from real neurons that the timing of\ninputs is extremely relevant for their spiking.\n\nThe LIF model has four parameters, `gL, EL, C, Vth, I` and we define it in the `lif(u, p, t)` function."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "using DifferentialEquations\nusing Plots\ngr()\n\nfunction lif(u,p,t);\n    gL, EL, C, Vth, I = p\n    (-gL*(u-EL)+I)/C\nend"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "Our system is described by one differential equation: `(-gL*(u-EL)+I)/C`, where\n`u` is the voltage, `I` is the input, `gL` is the leak conductance, `EL` is the\nequilibrium potential of the leak conductance and `C` is the membrane capacitance.\nGenerally, any change of the voltage is slowed down (filtered) by the membrane\ncapacitance. That's why we divide the whole equation by `C`. Without any external\ninput, the voltage always converges towards `EL`. If `u` is larger than `EL`,\n`u` decreases until it is at `EL`. If `u` is smaller than `EL`, `u` increases\nuntil it is at `EL`. The only other thing that can change the voltage is the\nexternal input `I`.\n\nOur `lif` function requires a certain parameter structure because it will need\nto be compatible with the `DifferentialEquations` interface. The input signature\nis `lif(u, p, t)` where `u` is the voltage, `p` is the collection of the parameters\nthat describe the equation and `t` is time. You might wonder why time does not\nshow up in our equation, although we need to calculate the change in voltage\nwith respect to time. The ODE solver will take care of time for us. One of\nthe advantages of the ODE solver as opposed to calculating the change of\nu in a for loop is that many ODE solver algorithm can dynamically adjust the\ntime step in a way that is efficient and accurate.\n\nOne crucial thing is still missing however. This is supposed to be a model of\nneural spiking, right? So we need a mechanism that recognizes the spike and\nhyperpolarizes `u` in response. For this purpose we will use callbacks.\nThey can make discontinuous changes to the model when certain conditions are met."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "function thr(u,t,integrator)\n    integrator.u > integrator.p[4]\nend\n\nfunction reset!(integrator)\n    integrator.u = integrator.p[2]\nend\n\nthreshold = DiscreteCallback(thr,reset!)\ncurrent_step= PresetTimeCallback([2,15],integrator -> integrator.p[5] += 210.0)\ncb = CallbackSet(current_step,threshold)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "Our condition is `thr(u,t,integrator)` and the condition kicks in when `integrator.u > integrator.p[4]` where `p[4]` is our threshold parameter `Vth`. Our effect of the condition is `reset!(integrator)`. It sets `u` back to the equilibrium potential `p[2]`. We then wrap both the condition and the effect into a `DiscreteCallback` called threshold. There is one more particularly handy callback called `PresetTimeCallback`. This one increases the input `p[5]` at `t=2` and `t=15` by `210.0`. Both callbacks are then combined into a `CallbackSet`. We are almost done to simulate our system we just need to put numbers on our initial voltage and parameters."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "u0 = -75\ntspan = (0.0, 40.0)\n# p = (gL, EL, C, Vth, I)\np = [10.0, -75.0, 5.0, -55.0, 0]\n\nprob = ODEProblem(lif, u0, tspan, p, callback=cb)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "Our initial voltage is `u0 = - 75`, which will be the same as our equilibrium potential, so we start at a stable point. Then we define the timespan we want to simulate. The time scale of the LIF as it is defined conforms roughly to milliseconds. Then we define our parameters as `p = [10.0, -75.0, 5.0, -55.0, 0]`. Remember that `gL, EL, C, Vth, I = p`. Finally we wrap everything into a call to `ODEProblem`. Can't forget the `CallbackSet`. With that our model is defined. Now we just need to solve it with a quick call to `solve`."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "sol = solve(prob)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "First of all the `solve` output tells us if solving the system generally worked. In this case we know it worked because the return code (`retcode`) say `Success`. Then we get the numbers for the timestep and the solution to `u`. The raw numbers are not super interesting to let's plot our solution."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "plot(sol)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "We see that the model is resting at `-75` while there is no input. At `t=2` the input increases by `210` and the model starts to spike. Spiking does not start immediately because the input first has to charge the membrane capacitance. Notice how once spiking starts it very quickly becomes extremely regular. Increasing the input again at `t=15` increases firing as we would expect but it is still extremely regular. This is one of the features of the LIF. The firing frequency is regular for constant input and a linear function of the input strength. There are ways to make LIF models less regular. For example we could use certain noise types at the input. We could also simulate a large number of LIF models and connect them synaptically. Instead of going into those topics, we will move on to the Izhikevich model, which is known for its ability to generate a large variety of spiking dynamics during constant inputs.\n\n## The Izhikevich Model\n[The Izhikevich model](https://www.izhikevich.org/publications/spikes.htm) is a two-dimensional model of neuronal spiking. It was derived from a bifurcation analysis of a cortical neuron. Because it is two-dimensional it can generate much more complex spike dynamics than the LIF model. The kind of dynamics depend on the four parameters and the input `a, b, c, d, I = p`. All the concepts are the same as above, expect for some minor changes to our function definitions to accomodate for the second dimension."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "#Izhikevichch Model\nusing DifferentialEquations\nusing Plots\n\nfunction izh!(du,u,p,t);\n    a, b, c, d, I = p\n\n    du[1] = 0.04*u[1]^2+5*u[1]+140-u[2]+I\n    du[2] = a*(b*u[1]-u[2])\nend"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "This is our Izhikevich model. There are two important changes here. First of all, note the additional input parameter `du`. This is a sequence of differences. `du[1]` corresponds to the voltage (the first dimension of the system) and `du[2]`corresponds to the second dimension. This second dimension is called `u` in the original Izhikevich work amnd it makes the notation a little annoying. In this tutorial I will generally stick to Julia and `DifferentialEquations` conventions as opposed to conventions of the specific models and `du` is commonly used. We will never define `du` ourselves outside of the function but the ODE solver will use it internally. The other change here is the `!` after our function name. This signifies that du will be preallocated before integration, which saves a lot of allocation time. In other words, du will be changed in place. Now we just need our callbacks to take care of spikes and increase the input."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "function thr(u,t,integrator)\n    integrator.u[1] >= 30\nend\n\nfunction reset!(integrator)\n    integrator.u[1] = integrator.p[3]\n    integrator.u[2] += integrator.p[4]\nend\n\nthreshold = DiscreteCallback(thr,reset!)\ncurrent_step= PresetTimeCallback(50,integrator -> integrator.p[5] += 10)\ncb = CallbackSet(current_step,threshold)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "One key feature of the Izhikevich model is that each spike increases our second dimension `u[2]` by a preset amount `p[4]`. Between spikes `u[2]` decays to a value that depends on `p[1]` and `p[2]` and the equilibrium potential `p[3]`. Otherwise the code is not too different from the LIF model. We will again need to define our parameters and we are ready to simulate."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "p = [0.02, 0.2, -50, 2, 0]\nu0 = [-65, p[2]*-65]\ntspan = (0.0, 300)\n\nprob = ODEProblem(izh!, u0, tspan, p, callback=cb)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "sol = solve(prob);\nplot(sol, vars=1)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "This spiking type is called chattering. It fires with intermittent periods of silence. Note that the input starts at `t=50` and remain constant for the duration of the simulation. One of mechanisms that sustains this type of firing is the spike induced hyperpolarization coming from our second dimension, so let's look at this variable."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "plot(sol, vars=2)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "Our second dimension `u[2]` increases with every spike. When it becomes too large and the system cannot generate another spike until `u[2]`has decayed to a value small enough that spiking can resume. This process repeats. In this model spiking is no longer regular like it was in the LIF. Here we have two frequencies, the frequency during the spiking state and the frequency between spiking states. The LIF model was dominated by one single frequency that was a function of the input strength. Let's see if we can generate another spiking type by changing the parameters."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "p = [0.02, 0.2, -65, 8, 0]\nu0 = [-65, p[2]*-65]\ntspan = (0.0, 300)\n\nprob = ODEProblem(izh!, u0, tspan, p, callback=cb)\nsol = solve(prob);\nplot(sol, vars=1)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "This type is called regularly spiking and we created it just by lowering `p[3]` and increasing `p[4]`. Note that the type is called regularly spiking but it is not instantaneously regular. The instantenous frequency is higher in the beginning. This is called spike frequency adaptation and is a common property of real neurons. There are many more spike types that can be generated. Check out the [original Izhikevich work](https://www.izhikevich.org/publications/spikes.htm) and create your own favorite neuron!\n\n## Hodgkin-Huxley Model\nThe Hodgkin-Huxley (HH) model is our first so called biophysically realistic model. This means that all parameters and mechanisms of the model represent biological mechanisms. Specifically, the HH model simulates the ionic currents that depolarize and hyperpolarize a neuron during an action potential. This makes the HH model four-dimensional. Let's see how it looks."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "using DifferentialEquations\nusing Plots\n\n# Potassium ion-channel rate functions\nalpha_n(v) = (0.02 * (v - 25.0)) / (1.0 - exp((-1.0 * (v - 25.0)) / 9.0))\nbeta_n(v) = (-0.002 * (v - 25.0)) / (1.0 - exp((v - 25.0) / 9.0))\n\n# Sodium ion-channel rate functions\nalpha_m(v) = (0.182*(v + 35.0)) / (1.0 - exp((-1.0 * (v + 35.0)) / 9.0))\nbeta_m(v) = (-0.124 * (v + 35.0)) / (1.0 - exp((v + 35.0) / 9.0))\n\nalpha_h(v) = 0.25 * exp((-1.0 * (v + 90.0)) / 12.0)\nbeta_h(v) = (0.25 * exp((v + 62.0) / 6.0)) / exp((v + 90.0) / 12.0)\n\nfunction HH!(du,u,p,t);\n    gK, gNa, gL, EK, ENa, EL, C, I = p\n    v, n, m, h = u\n\n    du[1] = ((gK * (n^4.0) * (EK - v)) + (gNa * (m ^ 3.0) * h * (ENa - v)) + (gL * (EL - v)) + I) / C\n    du[2] = (alpha_n(v) * (1.0 - n)) - (beta_n(v) * n)\n    du[3] = (alpha_m(v) * (1.0 - m)) - (beta_m(v) * m)\n    du[4] = (alpha_h(v) * (1.0 - h)) - (beta_h(v) * h)\nend"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "We have three different types of ionic conductances. Potassium, sodium and the leak. The potassium and sodium conducance are voltage gated. They increase or decrease depending on the voltage. In ion channel terms, open channels can transition to the closed state and closed channels can transition to the open state. It's probably easiest to start with the potassium current described by `gK * (n^4.0) * (EK - v)`. Here gK is the total possible conductance that we could reach if all potassium channels were open. If all channels were open, `n` would equal 1 which is usually not the case. The transition from open state to closed state is modeled in `alpha_n(v)` while the transition from closed to open is in `beta_n(v)`. Because potassium conductance is voltage gated, these transitions depend on `v`. The numbers in `alpha_n; beta_n` were calculated by Hodgkin and Huxley based on their extensive experiments on the squid giant axon. They also determined, that `n` needs to be taken to the power of 4 to correctly model the amount of open channels.\n\nThe sodium current is not very different but it two gating variables, `m; h` instead of one. The leak conductance gL has no gating variables because it is not voltage gated. Let's move on to the parameter. If you want all the details on the HH model you can find a great description [here](https://neuronaldynamics.epfl.ch/online/Ch2.S2.html)."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "current_step= PresetTimeCallback(100,integrator -> integrator.p[8] += 1)\n\n# n, m & h steady-states\nn_inf(v) = alpha_n(v) / (alpha_n(v) + beta_n(v))\nm_inf(v) = alpha_m(v) / (alpha_m(v) + beta_m(v))\nh_inf(v) = alpha_h(v) / (alpha_h(v) + beta_h(v))\n\np = [35.0, 40.0, 0.3, -77.0, 55.0, -65.0, 1, 0]\nu0 = [-60, n_inf(-60), m_inf(-60), h_inf(-60)]\ntspan = (0.0, 1000)\n\nprob = ODEProblem(HH!, u0, tspan, p, callback=current_step)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "For the HH model we need only one callback. The PresetTimeCallback that starts our input current. We don't need to reset the voltage when it reaches threshold because the HH model has its own repolarization mechanism. That is the potassium current, which activates at large voltages and makes the voltage more negative. The three functions `n_inf; m_inf; h_inf` help us to find good initial values for the gating variables. Those functions tell us that the steady-state gating values should be for the initial voltage. The parameters were chosen in a way that the properties of the model roughly resemble that of a cortical pyramidal cell instead of the giant axon Hodgkin and Huxley were originally working on."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "sol = solve(prob);\nplot(sol, vars=1)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "That's some good regular voltage spiking. One of the cool things about a biophysically realistic model is that the gating variables tell us something about the mechanisms behind the action potential. You might have seen something like the following plot in a biology textbook."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "plot(sol, vars=[2,3,4], tspan=(105.0,130.0))"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "That's it for the tutorial on spiking neural systems."
      ],
      "metadata": {}
    }
  ],
  "nbformat_minor": 2,
  "metadata": {
    "language_info": {
      "file_extension": ".jl",
      "mimetype": "application/julia",
      "name": "julia",
      "version": "1.4.2"
    },
    "kernelspec": {
      "name": "julia-1.4",
      "display_name": "Julia 1.4.2",
      "language": "julia"
    }
  },
  "nbformat": 4
}