7-Reduced-Dynamical-Models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "# Week 10: Reduced Dynamical Neuron Models\n",
    "\n",
    "Last week we discussed dynamical systems theory for 1-D and 2-D systems. We learned about the following key concepts:\n",
    "\n",
    "- state variables, parameters, differential equations\n",
    "- stable and unstable equilibria\n",
    "- topological equivalence and local linearization\n",
    "- limit cycles and saddle equilibria\n",
    "- bifurcations\n",
    "\n",
    "We'll conclude our investigation of 2-D systems with a discussion of a bifurcation that's not found in 1-D systems, and then look at some simple 2-D neuron models that can recapitulate many of the behaviors of more complex, realistic models.\n",
    "\n",
    "This notebook has a review of the material covered in lecture, followed by three exercises at the end.\n",
    "\n",
    "### Additional resources\n",
    "\n",
    "Many of the figures in this notebook are from Chapters 4 and 8 in *Dynamical systems in neuroscience*, by Eugene Izhikevich (2006).\n",
    "\n",
    "To see some really cool visualizations of dynamical systems, check out this [video textbook](https://www.youtube.com/watch?v=fufF25O4vD8&list=PL8erL0pXF3JZqdlYIfTTyibOqSqwzRdVV&index=1) by Robert Ghrist."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "# Lecture Review\n",
    "\n",
    "## Bifurcations in 2-D systems\n",
    "\n",
    "Just as we saw in 1-D dynamical systems, when an equilibrium has one or more zero-valued eigenvalues, it indicates that the system may be undergoing a **bifurcation**: a qualitative shift in dynamics caused by the loss or gain of equilibria.\n",
    "\n",
    "2-D systems can bifurcate in all the ways that 1-D systems can - saddle nodes, pitchforks, etc - but also in new and more exciting ways."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Heteroclinic orbits\n",
    "\n",
    "As a refresher, here's our simple persistent sodium + potassium model:\n",
    "\n",
    "\\begin{align}\n",
    "C_\\mathrm{m} \\dot{V} & = g_L(E_L - V) + \\overline{g}_\\mathrm{Na}m_\\infty(V)(E_\\mathrm{Na} - V) + \\overline{g}_\\mathrm{K}n(E_\\mathrm{K} - V) + I_e \\\\\n",
    "\\dot{n} & = \\frac{n_\\infty(V) - n}{\\tau_n(V)} \\\\\n",
    "\\end{align}\n",
    "\n",
    "When the injected current is low, there are three equilibria: a node, a saddle, and an unstable focus. The saddle and node are connected by **heteroclinic orbits**:\n",
    "\n",
    "<img src=\"https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l12_i4.26.png\" alt=\"heteroclinic orbits\" style=\"width: 400px;\"/>\n",
    "\n",
    "This system will always wind up back at the stable node. Why?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Homoclinic orbits\n",
    "\n",
    "It's also possible, though rare, for an orbit to loop back on the unstable equilibrium where it started. For example, with slightly higher $I_\\mathrm{e}$:\n",
    "\n",
    "<img src=\"https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l12_i4.27.png\" alt=\"homoclinic orbit\" style=\"width: 400px;\"/>\n",
    "\n",
    "If we increase $I_\\mathrm{e}$ further, we get a *saddle homoclinic orbit* bifurcation. We won't go into detail on this kind of bifurcation, but more can be found in Chapter 6 of Izhikevich (2006)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Saddle-node bifurcation\n",
    "\n",
    "What happens as the two equilibria get closer together? The stable node and saddle will eventually meet, forming a *saddle-node equilibrium*:\n",
    "\n",
    "<img src=\"https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l12_i4.28.png\" alt=\"saddle-node equilibrium\" style=\"width: 400px;\"/>\n",
    "\n",
    "Further increase of $I_\\mathrm{e}$ causes this equilibrium to disappear, forming a limit cycle around the unstable focus:\n",
    "\n",
    "<img src=\"https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l12_i4.30.png\" alt=\"limit cycle\" style=\"width: 400px;\"/>\n",
    "\n",
    "Hopefully this reminds you of the 1-D saddle-node bifurcation!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Saddle-node spiking\n",
    "\n",
    "If we apply a ramp current to this model, we will pass through the bifurcation, and the system will abruptly transition from rest to spiking:\n",
    "\n",
    "<img src=\"https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l12_i4.29.png\" alt=\"saddle-node spiking\" style=\"width: 500px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Attractor ruins\n",
    "\n",
    "In the phase plane, this transition looks like this:\n",
    "\n",
    "<img src=\"https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l12_i4.31.png\" alt=\"saddle-node spiking phase plane\" style=\"width: 500px;\"/>\n",
    "\n",
    "Just like the 1-D case, the 2-D saddle-node bifurcation leaves *attractor ruins* where the system moves slowly. As a consequence, a system undergoing saddle-node bifurcation can spike very slowly around the rheobase. This is called **Type I excitability**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Saddle-node bifurcation diagram\n",
    "\n",
    "The saddle-node bifurcation diagram for this system looks a lot like the 1-D system:\n",
    "\n",
    "<img src=\"https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l12_i4.32.png\" alt=\"saddle-node bifurcation diagram\" style=\"width: 500px;\"/>\n",
    "\n",
    "To generate this plot, we can simply replace the state-variable $n$ with its steady-state value and solve for the points where $\\dot{V} = 0$:\n",
    "\n",
    "$$\n",
    "0 = g_L(E_L - V) + \\overline{g}_\\mathrm{Na}m_\\infty(V)(E_\\mathrm{Na} - V) + \\overline{g}_\\mathrm{K}n_\\infty(V)(E_\\mathrm{K} - V) + I_e\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Hopf Bifurcation\n",
    "\n",
    "Our model can undergo a different kind of bifurcation if we change the threshold of the potassium current, such that the $n$-nullcline only intersects with the $V$-nullcline at a single point. When $I_\\mathrm{e}$ is low, this equilibrium is a *stable focus*:\n",
    "\n",
    "<img src=\"https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l12_i4.34a.png\" alt=\"low-threshold stable focus\" style=\"width: 300px;\"/>\n",
    "\n",
    "As the injected current increases, however, the system starts to oscillate around the equilibrium:\n",
    "\n",
    "<img src=\"https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l12_i4.34b.png\" alt=\"low-threshold bifurcation\" style=\"width: 300px;\"/>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Hopf-bifurcation spiking\n",
    "\n",
    "If we apply a ramp current to the low-threshold K model, we see a very different transition to spiking:\n",
    "\n",
    "<img src=\"https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l12_i4.33.png\" alt=\"Hopf spiking\" style=\"width: 500px;\"/>\n",
    "\n",
    "Note that the period of the spiking remains relatively constant, even though the amplitude increases."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Hopf bifurcation in the phase plane\n",
    "\n",
    "Oscillation is evidence of a limit cycle, and indeed this is what we see in the phase plane as the current increases. \n",
    "\n",
    "The stable equilibrium becomes unstable, and becomes surrounded by a limit cycle:\n",
    "\n",
    "<img src=\"https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l12_i4.34.png\" alt=\"low-threshold bifurcation\" style=\"width: 400px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Bifurcation plot\n",
    "\n",
    "Hopf bifurcations are not as easy to find, because the equilibrium does not appear or disappear. We have to simulate the model or determine when the real component of the Jacobian eigenvalues changes sign.\n",
    "\n",
    "<img src=\"https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l12_i4.35.png\" alt=\"Hopf spiking\" style=\"width: 500px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bifurcation and excitability types\n",
    "\n",
    "Because the Hopf bifurcation is characterized by a relative constant firing frequency, with no attractor ruins, a neuron that exhibits this kind of dynamics exhibits a discontinuous f-I plot, or **Type II excitability** (right):\n",
    "\n",
    "<img src=\"https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l12_excitability.png\" alt=\"type I and type II firing\" style=\"width: 500px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Simple Neuron Models\n",
    "\n",
    "One take-away from our dive into dynamical systems theory is that we don't necessarily need to model every detail of neuron biophysics.\n",
    "\n",
    "Instead, we might be able to approximate many aspects of neural behavior using models with only a few state variables.\n",
    "\n",
    "This would have significant advantages for the explanatory power of the models and for parameter inference."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Integrate-and-Fire\n",
    "\n",
    "Let's start with one of the simplest models, which represents neurons as a ohmic leak current and a spike threshold:\n",
    "\n",
    "<img src=\"https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l12_i8.1.png\" alt=\"LIF model\" style=\"width: 400px;\"/>\n",
    "\n",
    "The membrane dynamics are identical to our passive model:\n",
    "\n",
    "$$\n",
    "C_\\mathrm{m}\\dot{V} = I_\\mathrm{e} + g_L(E_L - V)\n",
    "$$\n",
    "\n",
    "However, when $V \\geq E_\\mathrm{thresh}$, the model is said to fire an action potential and the voltage is **reset** to $E_K$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The integrate-and-fire (IF) model can be linearly rescaled to this canonical form:\n",
    "\n",
    "$$\n",
    "\\dot{v} = b - v \\quad \\mathrm{if} \\; v = 1, v \\rightarrow 0\n",
    "$$\n",
    "\n",
    "The IF model has a number of neuron-like properties:\n",
    "\n",
    "- all-or-nothing spiking\n",
    "- well-defined threshold\n",
    "- relative refractory period\n",
    "- type I excitability\n",
    "\n",
    "From a practical standpoint, the IF model is very fast and easy to simulate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Quadratic Integrate-and-Fire\n",
    "\n",
    "The LIF model doesn't really have dynamics *per se* (neither a saddle-node nor a Hopf bifurcation), and its behavior around the spike threshold is decidedly not neuronal.\n",
    "\n",
    "The main reason is that it lacks the nonlinear \"n\" shape of the $V$-nullcline. We can achieve a similar topology while maintaining computational efficiency by replacing the $-v$ with $v^2$:\n",
    "\n",
    "$$\n",
    "\\dot{v} = b + v^2 \\quad \\mathrm{if} \\; v = v_\\mathrm{peak}, v \\rightarrow v_\\mathrm{reset}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### QIF bifurcation\n",
    "\n",
    "Unlike the LIF model, the QIF model has bifurcations. Indeed, $\\dot{v} = b + v^2$ is the *topological normal form* for saddle node bifurcations along the $b$ parameter. Depending on what the value of $v_\\mathrm{reset}$ is, there are several other kinds of possible bifurcations in this model:\n",
    "\n",
    "<img src=\"https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l12_i8.3.png\" alt=\"QIF bifurcations\" style=\"width: 400px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The QIF model exhibits some key properties that are seen in real neurons but not the LIF model:\n",
    "\n",
    "- acceleration toward spike initiation (soft threshold)\n",
    "- saddle-node bifurcations\n",
    "\n",
    "<img src=\"images/l12_i8.7.png\" alt=\"QIF model spiking\" style=\"width: 500px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Izhikevich Model\n",
    "\n",
    "The QIF model can be further improved (i.e., more able to generate real neuronal firing behaviors) by adding a recovery variable $u$:\n",
    "\n",
    "\\begin{align}\n",
    "\\dot{v} & = I + v^2 - u \\quad &\\mathrm{if} \\; v \\geq v_\\mathrm{peak}, \\\\\n",
    "\\dot{u} & = a(bv - u) \\quad & v \\rightarrow c, u \\rightarrow u + d\n",
    "\\end{align}\n",
    "\n",
    "Depending on the values of $a$ and $b$, the $u$ and $v$ nullclines can be made to cross at a single point (resonator, Hopf bifurcation) or at two points (integrator, saddle-node bifurcation):\n",
    "\n",
    "<img src=\"https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l12_i8.5.png\" alt=\"Izhikevich model phase portrait\" style=\"width: 500px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "With suitably tuned parameters, the Izhikevich model can generate a wide range of spiking patterns:\n",
    "\n",
    "<img src=\"https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l12_i8.8.png\" alt=\"Izhikevich model phase portrait\" style=\"width: 500px;\"/>\n",
    "\n",
    "See Figure 8.8 in Izhikevich (2006) for more patterns and a link to the code that generates them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Exercises\n",
    "\n",
    "## Simulating reduced models\n",
    "\n",
    "The big advantages of the LIF, Izhekevich, and other reduced models are:\n",
    "\n",
    "- they have few parameters $\\rightarrow$ high explanatory power\n",
    "- they have few state variables $\\rightarrow$ low computational demands\n",
    "\n",
    "This makes it feasible to simulate large networks of neurons fast enough to gain some insight into how various parameters, including connectivity within the network, influence larger-scale dynamics and neurocomputational properties."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dependencies\n",
    "\n",
    "To run the simulations in this notebook, you'll need to have the [brian2](https://briansimulator.org/) simulator installed. If you're running from binder or a docker image, this is already taken care of, but if you're using a local anaconda installation, you need to run the following cell to set up before running any of the following cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install brian2==2.5.0.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# load matplotlib inline mode\n",
    "%matplotlib inline\n",
    "\n",
    "# \"star imports\" are not generally considered good practice, but this does simply our demo somewhat\n",
    "from brian2 import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Leaky integrate and fire\n",
    "\n",
    "To create a model neuron in brian, we specify the differential equations directly. Variable names are inferred from the equations, but parameters like `tau` have to be defined as python variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_scope() # this clears all the models\n",
    "\n",
    "tau = 10 * ms\n",
    "eqs = '''\n",
    "dv/dt = (1-v)/tau : 1\n",
    "'''\n",
    "\n",
    "# create a neuron (or rather, a group of neurons with 1 member) and specify how to integrate numerically\n",
    "G = NeuronGroup(1, eqs, method='exact')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "To run the model, we call `run()` with the duration of the integration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"initial state:\", G.v[0])\n",
    "run(100 * ms)\n",
    "print(\"final state:\", G.v[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "If we want to see the time course of the voltage, we need to create a `StateMonitor`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_scope()\n",
    "G = NeuronGroup(1, eqs, method='exact')\n",
    "M = StateMonitor(G, 'v', record=True)\n",
    "run(30 * ms)\n",
    "plot(M.t/ms, M.v[0])\n",
    "xlabel('Time (ms)')\n",
    "ylabel('v');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Generating spikes\n",
    "\n",
    "To make our model fire spikes, we need to specify a threshold and a reset rule. These are supplied to the `NeuronGroup` function as additional arguments. We also need to specify an additional monitor to record spike times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_scope()\n",
    "\n",
    "tau = 10*ms\n",
    "eqs = '''\n",
    "dv/dt = (1-v)/tau : 1\n",
    "'''\n",
    "\n",
    "G = NeuronGroup(1, eqs, threshold='v>0.8', reset='v = 0', method='exact')\n",
    "statemon = StateMonitor(G, 'v', record=True)\n",
    "spikemon = SpikeMonitor(G)\n",
    "\n",
    "run(50*ms)\n",
    "\n",
    "plot(statemon.t/ms, statemon.v[0])\n",
    "for t in spikemon.t:\n",
    "    axvline(t/ms, ls='--', c='C1', lw=3)\n",
    "xlabel('Time (ms)')\n",
    "ylabel('v');\n",
    "\n",
    "print(\"spike times:\", spikemon.t[:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Time-varying driving current\n",
    "\n",
    "Let's see if we can simulate the response to an increasing ramp current. We can also specify parameters that change as a function of time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_scope()\n",
    "tau = 10*ms\n",
    "slope = 0.01 / ms\n",
    "\n",
    "eqs = '''\n",
    "dv/dt = (I-v)/tau : 1\n",
    "I = slope * t : 1\n",
    "'''\n",
    "\n",
    "G = NeuronGroup(1, eqs, threshold='v>0.8', reset='v = 0', method='euler')\n",
    "M = StateMonitor(G, variables=True, record=True)\n",
    "\n",
    "run(200 * ms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = subplots(nrows=2, ncols=1, sharex=True, figsize=(9, 6))\n",
    "axes[0].plot(M.t/ms, M.v[0])\n",
    "axes[1].plot(M.t/ms, M.I[0])\n",
    "axes[0].set_ylabel('V')\n",
    "axes[1].set_ylabel('I')\n",
    "axes[1].set_xlabel(\"Time (ms)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Quadratic Integrate and Fire\n",
    "\n",
    "Now let's simulate a QIF model. This is still a 1-D system, but the differential equation has changed. Also, because this model is capable of spiking (rather than us imputing a spike when the voltage crosses a threshold), we set the threshold much higher."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_scope()\n",
    "tau = 10*ms\n",
    "# driving current\n",
    "I = 1\n",
    "\n",
    "eqs = '''\n",
    "dv/dt = (I + v**2)/tau : 1\n",
    "'''\n",
    "\n",
    "G = NeuronGroup(1, eqs, threshold='v>40', reset='v = 0', method='euler')\n",
    "M = StateMonitor(G, variables=True, record=True)\n",
    "\n",
    "# initial state\n",
    "G.v = 0\n",
    "\n",
    "run(200 * ms)\n",
    "plot(M.t/ms, M.v[0])\n",
    "ylabel('V (mV)')\n",
    "xlabel(\"Time (ms)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Exercise 1\n",
    "\n",
    "This is the parameter space for the QIF model (Note that the `I` parameter in our model is $b$ in the figure).\n",
    "\n",
    "<img src=\"https://gracula.psyc.virginia.edu/public/courseware/comp-neurosci/images/l12_i8.3.png\" alt=\"QIF bifurcations\" style=\"width: 400px;\"/>\n",
    "\n",
    "**A** Which dynamical regime (i.e., which region in the parameter space) does the model in the code cell above correspond to? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**B** Copy the previous code cell, edit it to set the initial voltage to a different value, then run the code. What dynamical behavior do you observe? Does the behavior depend on the initial\n",
    "condition in this part of the parameter space?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**C**  Copy and paste the code above and edit the driving current `I` (and if necessary) the initial voltage to place the model in the *excitable* regime. To demonstrate excitability, show that the model spikes once if its initial condition is above the threshold."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**D** Adjust the parameter values (`I` and `Vreset`) to place the system in the the *bistable* regime. To demonstrate bistability, you need to show that the model spikes repetitively if starts above threshold, but remains at rest if it's below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Izhikevich model\n",
    "\n",
    "To implement the Izhikevich model, we need to add an additional state variable, $u$.\n",
    "\n",
    "\\begin{align}\n",
    "\\dot{v} & = I + v^2 - u \\quad &\\mathrm{if} \\; v \\geq v_\\mathrm{reset}, \\\\\n",
    "\\dot{u} & = a(bv - u) \\quad & v \\rightarrow c, u \\rightarrow u + d\n",
    "\\end{align}\n",
    "\n",
    "We're going to use an alternate form that allows us to separately specify the resting potential and threshold:\n",
    "\n",
    "\\begin{align}\n",
    "C\\dot{v} & = k(v - v_\\mathrm{rest})(v - v_\\mathrm{thresh}) - u + I\\quad &\\mathrm{if} \\; v \\geq v_\\mathrm{reset}, \\\\\n",
    "\\dot{u} & = a(b(v - v_\\mathrm{rest}) - u) \\quad & v \\rightarrow c, u \\rightarrow u + d\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_scope()\n",
    "# from Izhikevich (2006) Fig 8.6a\n",
    "C  = 100 * ms\n",
    "vr = -60\n",
    "vt = -40\n",
    "k  = 0.7\n",
    "a  = 0.03 / ms\n",
    "b  = -2\n",
    "c  = -50\n",
    "d  = 100\n",
    "vpeak = 35\n",
    "I = TimedArray((np.arange(100) > 7) * 70, dt=10 * ms)\n",
    "eqs = '''\n",
    "dv/dt = (k * (v - vr) * (v - vt) - u + I(t))/C : 1\n",
    "du/dt = a * (b * (v - vr) - u) : 1\n",
    "'''\n",
    "\n",
    "resetters = '''\n",
    "v = c\n",
    "u += d\n",
    "'''\n",
    "\n",
    "G = NeuronGroup(1, eqs, threshold='v>vpeak', reset=resetters, method='euler')\n",
    "M = StateMonitor(G, variables=True, record=True)\n",
    "G.v = vr\n",
    "\n",
    "run(1000 * ms)\n",
    "plot(M.t/ms, M.v[0])\n",
    "ylabel('V (mV)')\n",
    "xlabel(\"Time (ms)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Exercise 2\n",
    "\n",
    "**A**. Draw or plot the phase portrait for the Izhikevich model with $I = 0$ using the parameter values above. Recall that the phase portrait\n",
    "shows the vector field and the nullclines. I recommend finding the nullclines first. To find the $v$-nullcline, set $dv/dt = 0$ and solve for $v$ as a function of $u$. To find the $u$-nullcline, set $du/dt = 0$ and solve for $v$ as a function of $u$. You do not need to indicate the fixed points or their stability.\n",
    "\n",
    "Set the x-axis limits of your plot from -80 to 0 for $v$ and plot the nullclines. The min and max of the $u$-nullcline will give you the limits to use for the y-axis. When creating your meshgrid, use `np.linspace` to pick ~20 values between the upper and lower limits for $u$ and $v$.\n",
    "\n",
    "If you are not sure how to use these functions (`np.linspace`, `np.meshgrid`, `plt.quiver`), take a look at the examples in notebook 6."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**B**. What kind of bifurcation does this particular set of parameter values produce as $I$ is increased? *Note*: you do not need to\n",
    "simulate anything to answer this question. You should be able to visualize how the fixed points will change as $I$ is increased."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**C**. Find a new set of $b$ and $v_{rest}$ values that will alter the intersection of the $u$ and $v$ nullclines such that they only intersect once in the range of $v \\in (-80, 0)$. Plot the phase portrait. What kind of bifurcation should the system exhibit now?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**D**. Copy the simulation code from above and paste it below this cell, setting the values of $b$ and $v_{rest}$ accordingly. Test\n",
    "if your prediction is correct (you will need show how the system behaves on both sides of the bifurcation)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3\n",
    "\n",
    "Copy the simulation code from above and adjust the parameters to make the model burst (i.e. fire a series of closely spaced spikes and then stop for some period of time). (Hint: consult chapter 8 in Izhikevich)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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.4"
  },
  "rise": {
   "scroll": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}