6-Optional-Dynamical-Systems.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Week 9: Dynamical Systems (optional notebook)\n",
    "\n",
    "As first shown by Hodgkin and Huxley, the biophysics of neurons can be described as a system of ordinary differential equations. \n",
    "\n",
    "One route we can take with such models is to increase their realism by adding currents and compartments that correspond to the channels expressed by real neurons, and to their spatially elaborated morphology. This realism comes at the cost of more parameters and variables, which take more computational power to simulate and which can be difficult to constrain with experimental data.\n",
    "\n",
    "Another route, which we will take today, is to look for general principles of neural dynamics that are independent of the specific channels and their spatial distribution. In other words, what can we learn about neurons by studying them through the mathematical formalism of **dynamical systems**?\n",
    "\n",
    "### Before class\n",
    "\n",
    "Read chapters 3-4 of *Dynamical Systems in Neuroscience : The Geometry of Excitability and Bursting* (available as an e-book through the library or on Collab). Many of the figures in this notebook come from this text.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# load matplotlib inline mode\n",
    "%matplotlib inline\n",
    "\n",
    "# import some useful libraries\n",
    "import numpy as np                # numerical analysis linear algebra\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt   # plotting\n",
    "\n",
    "# set some style options\n",
    "mpl.rcParams['image.origin'] = 'lower'\n",
    "mpl.rcParams['image.aspect'] = 'auto'\n",
    "mpl.rcParams['image.cmap'] = 'jet'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Key terms\n",
    "\n",
    "DST is (or attempts to be) a general theory of how different kinds of dynamical systems behave.\n",
    "\n",
    "- **dynamical system**: a system of equations that describes how a system changes in time (e.g. ODEs)\n",
    "- **state variables**: time-varying quantities that describe the state of a system\n",
    "- **parameters**: quantities in the dynamical system that don't change\n",
    "- **state space**: the space defined by the state variables of the system. \n",
    "- **parameter space**: the space defined by the parameters of the system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Hodgkin-Huxley Model\n",
    "\n",
    "For example, the HH model is a dynamical system defined by these four differential equations:\n",
    "\n",
    "\\begin{align}\n",
    "C_\\mathrm{m} \\dot{V} & = g_L(E_L - V) + \\overline{g}_\\mathrm{Na}m^3h(E_\\mathrm{Na} - V) + \\overline{g}_\\mathrm{K}n^4(E_\\mathrm{K} - V) + I_e\\\\\n",
    "\\dot{m} & = \\frac{m_\\infty(V) - m}{\\tau_m(V)} \\\\\n",
    "\\dot{h} & = \\frac{h_\\infty(V) - h}{\\tau_h(V)} \\\\\n",
    "\\dot{n} & = \\frac{n_\\infty(V) - n}{\\tau_n(V)} \\\\\n",
    "\\end{align}\n",
    "\n",
    "What are the **state variables**, and what are the **parameters** in this model?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## One-dimensional neuron\n",
    "\n",
    "Let's start with the simplest possible dynamical system, with one state variable (and therefore a one-dimensional state space).\n",
    "\n",
    "We will base this loosely on real neuron dynamics by considering a **persistent sodium current** that does not inactivate.\n",
    "\n",
    "\\begin{align}\n",
    "C_\\mathrm{m} \\dot{V} & = g_L(E_L - V) + \\overline{g}_\\mathrm{Na}m(E_\\mathrm{Na} - V) + I_e\\\\\n",
    "\\dot{m} & = \\frac{m_\\infty(V) - m}{\\tau_m(V)} \\\\\n",
    "\\end{align}\n",
    "\n",
    "To further simplify things, let's assume that the kinetics of $m$ are very fast, with $\\tau_m(V) \\ll C_\\mathrm{m}R_\\mathrm{m}$. That means that $m \\approx m_\\infty(V)$, so we can drop the second equation to get\n",
    "\n",
    "$$\n",
    "C_\\mathrm{m} \\dot{V} = g_L(E_L - V) + \\overline{g}_\\mathrm{Na}m_\\infty(V)(E_\\mathrm{Na} - V) + I_e\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## One-dimensional dynamical system\n",
    "\n",
    "Now $\\dot{V}$ is just a function of $V$. Let's see what this function looks like:\n",
    "\n",
    "<img src=\"images/l11_i3.5.png\" alt=\"1-D Na_p model derivative\" style=\"width: 400px;\"/>\n",
    "\n",
    "As we discussed last week, this ODE has to be **integrated** to determine how $V(t)$ evolves through time. It's complicated enough that we can't solve it analytically, but instead have to use **numerical integration**. Here are some integrated trajectories:\n",
    "\n",
    "<img src=\"images/l11_i3.6.png\" alt=\"1-D Na_p model solutions\" style=\"width: 400px;\"/>\n",
    "\n",
    "Note that we get different solutions depending on the values of the initial condition $V(0)$ and the parameter $I_\\mathrm{e}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Geometric analysis\n",
    "\n",
    "We can also understand a lot about how a dynamical system will behave simply by looking at the graph of $\\dot{V}$.\n",
    "\n",
    "The key features we want to identify are the **equilibrium points** and whether they are **stable** or **unstable**.\n",
    "\n",
    "An equilibrium point is any point in the state space where there is no movement. In this case, where $\\dot{V} = 0$.\n",
    "\n",
    "Equilibrium points are also called **attractors** (but not every attractor is an equilibrium, as we'll see later).\n",
    "\n",
    "<img src=\"images/l11_i3.8a.png\" alt=\"equilibrium points\" style=\"width: 500px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Stability\n",
    "\n",
    "A **stable equilibrium** is one where where small perturbations of the state result in the system returning to the equilibrium. This is called **convergence**. \n",
    "\n",
    "An **unstable equilibrium** is one where small perturbations result in the system leaving the equilibrium, or **diverging**.\n",
    "\n",
    "<img src=\"images/l11_i3.8.png\" alt=\"equilibrium points\" style=\"width: 400px;\"/>\n",
    "\n",
    "What differentiates the stable from the unstable equilibrium points in this system?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Eigenvalues\n",
    "\n",
    "If we denote the ODE as $\\dot{V} = F(V)$, then the stability of the system around some value of $V$ is determined by the derivative of $F(V)$:\n",
    "\n",
    "$$\n",
    "\\lambda = F'(V)\n",
    "$$\n",
    "\n",
    "Negative slopes around an equilibrium ($\\lambda < 0$) imply stability; positive slopes ($\\lambda > 0$) imply instability.\n",
    "\n",
    "The value of $\\lambda$ at some equilibrium point $V$ is also called an **eigenvalue**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Attraction domains\n",
    "\n",
    "Unstable equilibria define the boundaries of **attraction domains** or **basins**.\n",
    "\n",
    "<img src=\"images/l11_i3.12.png\" alt=\"attraction domains\" style=\"width: 350px;\"/>\n",
    "\n",
    "There is always an unstable equilibrium between any two stable equilibria. Why?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Thresholds\n",
    "\n",
    "Unstable equilibria also define **thresholds**. \n",
    "\n",
    "If the system is perturbed from a stable equilibrium above the threshold, then it will converge on a different stable equilibrium (or diverge to infinity).\n",
    "\n",
    "<img src=\"images/l11_i3.13a.png\" alt=\"threshold derivative\" style=\"width: 350px;\"/>\n",
    "\n",
    "Spike initiation in biological neurons can be seen as an example of threshold crossing:\n",
    "\n",
    "<img src=\"images/l11_i3.14.png\" alt=\"threshold spike\" style=\"width: 350px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Phase plots\n",
    "\n",
    "We can understand a lot about the behavior of a dynamical system simply by looking at the slope of $F(V)$. The equilibrium points are where $F'(V) = 0$, and everywhere else the movement of the system can be represented by arrows. This is called a **phase plot** or **phase portrait**:\n",
    "\n",
    "<img src=\"images/l11_i3.18.png\" alt=\"phase plot\" style=\"width: 350px;\"/>\n",
    "\n",
    "Indeed, the dynamics of the system as a whole is much more a function of the **topology** of the phase plot than it is by the specific locations of the equilibria. Big changes in $F(V)$ may have minimal effects on the dynamics if the equilibria simply shift loctions; in contrast, small changes in $F(V)$ can completely transform the dynamics if they introduce or eliminate equilibria. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Local equivalence (Hartman-Grobman Theorem)\n",
    "\n",
    "We can take this idea a step further. If what really matters is the topology of the system, not the exact form of $F(V)$, then we can understand how a system will behave around an equilibrium by approximating it with a straight line!\n",
    "\n",
    "<img src=\"images/l11_i3.21.png\" alt=\"local equivalence\" style=\"width: 350px;\"/>\n",
    "\n",
    "That is, a nonlinear dynamical system $\\dot{V} = F(V)$ near an equilibrium point $V_\\mathrm{eq}$ is topologically equivalent to the linear dynamical system\n",
    "\n",
    "$$\n",
    "\\dot{V} = \\lambda(V - V_\\mathrm{eq})\n",
    "$$\n",
    "\n",
    "provided that $\\lambda$ is nonzero.\n",
    "\n",
    "In contrast to the full dynamical system, which may not have an analytical solution, the solution to the linear approximation is just\n",
    "\n",
    "$$\n",
    "V(t) = V_\\mathrm{eq} + e^{\\lambda t}(V_0 - V_\\mathrm{eq})\n",
    "$$\n",
    "\n",
    "Thus, if $\\lambda$ is negative, the system will exponentially converge to the equilibrium value. If $\\lambda$ is positive, it will exponentially diverge."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Bifurcations\n",
    "\n",
    "What happens when $\\lambda$ is zero? The simple answer is that the system cannot be approximated by $\\dot{V} = \\lambda(V - V_\\mathrm{eq})$, but depends instead on higher-order derivatives of $F(V)$.\n",
    "\n",
    "The more complex answer is that $\\lambda$ is *often* zero when the system is at a transition point between topologies, for example between a bistable and a monostable state. This is called a **bifurcation**:\n",
    "\n",
    "<img src=\"images/l11_i3.22.png\" alt=\"bifurcation\" style=\"width: 350px;\"/>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Saddle-node bifurcation\n",
    "\n",
    "The simplest bifurcation occurs when a stable and unstable equilibrium collide and vanish (or emerge from nothingness). Let's look at what happens around the resting potential of our simple 1-D model as we change the injected current:\n",
    "\n",
    "<img src=\"images/l11_i3.26.png\" alt=\"saddle-node\" style=\"width: 400px;\"/>\n",
    "\n",
    "The transition between a stable + unstable equilibrium and no equilibrium occurs when $F(V)$ is tangent to zero. This bifurcation is called a **saddle-node** or **fold** bifurcation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Saddle-node bifurcations are characterized by **slow transitions**. Even when the equilibria have vanished, as the system passes through their **ruins** it slows down, because $F(V)$ is still close to zero.\n",
    "\n",
    "In some neuron types, spike initiation near the **rheobase** (the minimum current needed to initiate spiking) exhibits behavior typical of saddle-node bifurcation:\n",
    "\n",
    "<img src=\"images/l11_i3.28.png\" alt=\"saddle-node spike initiation\" style=\"width: 600px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Bifurcation diagrams\n",
    "\n",
    "Another geometric tool used to analyze dynamical systems is the **bifurcation plot**, which shows how equilibria and other attractors are affected by changes in a **parameter** of the system.\n",
    "\n",
    "For example, here is the bifurcation digram for our simple 1D model:\n",
    "\n",
    "<img src=\"images/l11_i3.32.png\" alt=\"bifurcation\" style=\"width: 400px;\"/>\n",
    "\n",
    "Below the rheobase, there are a stable and unstable equilibrium; above the rheobase this equilibrium vanishes. We can tell that this is a saddle-node bifurcation because of the way the equilibria coalese (going to the right) or appear \"out of thin air\" (going to the left)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 2D Model\n",
    "\n",
    "Now we'll take what we learned about dynamical systems from the 1-D example and extend it to a model that can actually generate spikes. To do this, we need to add a mechanism that will repolarize the cell membrane - a potassium current. Because this current has slower kinetics than sodium, we need to add a dynamical variable, $n(t)$.\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",
    "The state of the system is now defined by two variables, $V(t)$ and $n(t)$, and its evolution is defined by two nonlinear ODEs:\n",
    "\n",
    "\\begin{align}\n",
    "\\dot{V} & = f(V,n) \\\\\n",
    "\\dot{n} & = g(V,n)\n",
    "\\end{align}\n",
    "\n",
    "As before, we're going to take a geometric approach to understanding this dynamical system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Vector fields\n",
    "\n",
    "Because the system is 2-dimensional, we now have to plot it in an x-y plane. Because the value at each $(V,n)$ coordinate is itself a vector $(\\dot{V}, \\dot{n})$, our plot is a **vector field**.\n",
    "\n",
    "Let's look at some very simple vector fields. Here's one with a constant gradient pointing to the right:\n",
    "\n",
    "\\begin{align}\n",
    "\\dot{V} & = 1 \\\\\n",
    "\\dot{n} & = 0\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "V, n = np.meshgrid(np.arange(-10, 10), np.arange(-10, 10))\n",
    "dV = np.ones_like(V)\n",
    "dn = np.zeros_like(n)\n",
    "q = plt.quiver(V, n, dV, dn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Here's one where all the arrows point toward the origin:\n",
    "\n",
    "\\begin{align}\n",
    "\\dot{V} & = -V \\\\\n",
    "\\dot{n} & = -n\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dV = -V\n",
    "dn = -n\n",
    "q = plt.quiver(V, n, dV, dn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "This is called a saddle field:\n",
    "\n",
    "\\begin{align}\n",
    "\\dot{V} & = -n \\\\\n",
    "\\dot{n} & = -V\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dV = -n\n",
    "dn = -V\n",
    "q = plt.quiver(V, n, dV, dn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Trajectories\n",
    "\n",
    "We can visualize the evolution of a system on a vector plot by numerically integrating from a set of initial conditions. Note that different initial conditions in a saddle field can lead to very divergent outcomes.\n",
    "\n",
    "With some practice, you can approximate numerical integration just by looking at the flow of the arrows in a vector field."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import integrate\n",
    "\n",
    "def grad(t, X):\n",
    "    V, n = X\n",
    "    return (-n, -V)\n",
    "\n",
    "T = 3\n",
    "t = np.arange(0, T, 0.001)\n",
    "X0 = (-7, -8)\n",
    "ret = integrate.solve_ivp(grad, (0, T), X0, t_eval=t)\n",
    "plt.quiver(V, n, dV, dn)\n",
    "plt.plot(*X0, 'bo')\n",
    "plt.plot(ret.y[0], ret.y[1])\n",
    "\n",
    "X0 = (1, 2)\n",
    "ret = integrate.solve_ivp(grad, (0, T), X0, t_eval=t)\n",
    "plt.plot(*X0, 'ro')\n",
    "plt.plot(ret.y[0], ret.y[1], 'r')\n",
    "\n",
    "X0 = (-4, -4)\n",
    "ret = integrate.solve_ivp(grad, (0, T), X0, t_eval=t)\n",
    "plt.plot(*X0, 'go')\n",
    "plt.plot(ret.y[0], ret.y[1], 'g')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Equilibria and Nullclines\n",
    "\n",
    "Just as we could understand a 1-D system from its phase portrait by plotting equilibrium points and arrows, we can annotate a 2D vector field with the equivalent.\n",
    "\n",
    "However, in a 2D system, equilibrium points are only found when both components of the gradient are zero, i.e. where $(\\dot{V}, \\dot{n}) = (0, 0)$.\n",
    "\n",
    "There may also be many other points in the state space where one component of the gradient is zero. These lines or manifolds are called **nullclines**. For a 2-D system, there will be one nullcline for each variable, and the point(s) where the nullclines meet are equilibria. Let's look at the nullclines for our 2D neuron 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",
    "<img src=\"images/l11_i4.4.png\" alt=\"phase plot for morris-lecar model\" style=\"width: 400px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Where are the nullclines in the saddle field? Where is the equilibrium? Is it a stable or unstable equilibrium?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dV = -n\n",
    "dn = -V\n",
    "q = plt.quiver(V, n, dV, dn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Limit cycles\n",
    "\n",
    "In 1-D systems, trajectories can only converge on a stable equilibrium or diverge to infinity. 2-D systems can do this as well, but also have the ability to form closed loops or **periodic trajectories**.\n",
    "\n",
    "If a periodic trajectory is isolated (i.e. surrounded by regions that are non-periodic), it is called a **limit cycle**.\n",
    "\n",
    "Like equilibria, limit cycles can be stable or unstable.\n",
    "\n",
    "<img src=\"images/l11_i4.9.png\" alt=\"limit cycles\" style=\"width: 400px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Local linearization\n",
    "\n",
    "Just as we saw with 1-D systems, we can understand the behavior of a 2-D dynamical system by linearizing it around an equilibrium. Consider a general nonlinear 2-D system:\n",
    "\n",
    "\\begin{align}\n",
    "\\dot{x} & = f(x, y) \\\\\n",
    "\\dot{y} & = g(x, y) \\\\\n",
    "\\end{align}\n",
    "\n",
    "The linear approximation around some equilibrium $(x_0,y_0)$ is\n",
    "\n",
    "\\begin{align}\n",
    "f(x, y) = a(x - x_0) + b(y - y_0)\\\\\n",
    "g(x, y) = c(x - x_0) + d(y - y_0)\\\\\n",
    "\\end{align}\n",
    "\n",
    "The slopes $a,b,c,d$ are given by the **partial derivatives** of $f$ and $g$ with respect to $x$ and $y$:\n",
    "\n",
    "\\begin{align}\n",
    "a & = \\frac{\\partial f}{\\partial x}\\Bigr|_{(x_0,y_0)}, b = \\frac{\\partial f}{\\partial y}\\Bigr|_{(x_0,y_0)} \\\\\n",
    "c & = \\frac{\\partial g}{\\partial x}\\Bigr|_{(x_0,y_0)}, b = \\frac{\\partial g}{\\partial y}\\Bigr|_{(x_0,y_0)} \\\\\n",
    "\\end{align}\n",
    "\n",
    "These derivatives form a matrix with the coefficents of the linearized system:\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "\\dot{u}\\\\\n",
    "\\dot{w}\\\\\n",
    "\\end{pmatrix}\n",
    "= \n",
    "\\begin{pmatrix}\n",
    "a & b\\\\\n",
    "c & d\\\\\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "u\\\\\n",
    "w\\\\\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "This matrix, which is just a set of 4 numbers (for 2-D systems), is called the **Jacobian** of the equilibrium."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Eigenvalues and eigenvectors\n",
    "\n",
    "If you've taken linear algebra, you know that matrices can be analyzed in terms of their eigenvalues and eigenvectors.  For a $2 \\times 2$ matrix, there are two eigenvectors $v_1$ and $v_2$ and two corresponding eigenvalues $\\lambda_1$ and $\\lambda_2$.\n",
    "\n",
    "If you're unfamiliar with these concepts, or need a refresher on how to find the eigenvalues of a matrix, take a look at Chapter 5 in *Elementary Linear Algebra* by Kirkwood and Kirkwood (available as an ebook through the UVA library) or the relevant chapter in any other introductory text.\n",
    "\n",
    "Fortunately, Python gives us access to linear algebra libraries that can find the eigenvalues of any matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import linalg\n",
    "L = np.array([[0, -1], [-1, 0]], dtype='d')\n",
    "w, v = linalg.eig(L)\n",
    "print(\"eigenvalues:\", w)\n",
    "print(\"eigenvectors:\", v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import linalg\n",
    "L = np.array([[3, 9], [-4, -3]], dtype='d')\n",
    "w, v = linalg.eig(L)\n",
    "print(\"eigenvalues:\", w)\n",
    "print(\"eigenvectors:\", v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What system does this particular matrix correspond to?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Types of 2-D equilibria\n",
    "\n",
    "Just as we were able to determine whether an equilibrium in a 1-D system was stable or unstable from the local slope $\\lambda$, we can determine the stability of 2-D equilibria from $\\lambda_1$ and $\\lambda_2$. However, the situation is more complex because $\\lambda_1$ and $\\lambda_2$ can have the same or different signs, and can also be purely real or complex. \n",
    "\n",
    "Node equilibria occur when both eigenvalues are real and have the same sign.\n",
    "\n",
    "<img src=\"images/l11_i4.16.png\" alt=\"nodes\" style=\"width: 400px;\"/>\n",
    "\n",
    "\n",
    "Saddle equilibria occur when both eigenvalues are real and have different signs.\n",
    "\n",
    "<img src=\"images/l11_i4.17.png\" alt=\"nodes\" style=\"width: 400px;\"/>\n",
    "\n",
    "\n",
    "Focus equilibria occur when both eigenvalues are complex. They are stable if the real part is positive, and the imaginary part determines the frequency of rotation. \n",
    "\n",
    "<img src=\"images/l11_i4.18.png\" alt=\"nodes\" style=\"width: 400px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Systems with multiple equilibria\n",
    "\n",
    "Nonlinear systems like our 2-D neuron model can have multiple equilibria. That means that there is a Jacobian for each equilibrium, and the behavior of the system near each equilibrium can be quite distinct.\n",
    "\n",
    "What are the three equilibria types in this phase portrait?\n",
    "\n",
    "<img src=\"images/l11_i4.19.png\" alt=\"nodes\" style=\"width: 400px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Bistability\n",
    "\n",
    "When systems have more than one stable equilibrium or attractor, they exhibit **bistability**. We saw this in the 1-D case where two stable equilibria were separated by an unstable equilibrium, which acted as a threshold between two attraction domains.\n",
    "\n",
    "The same can be true in 2-D systems, but the separation is defined by a manifold rather than a single point. This manifold is called a **separatrix**.\n",
    "\n",
    "<img src=\"images/l11_i4.22.png\" alt=\"bistability\" style=\"width: 400px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Attraction domains may also be centered around limit cycles:\n",
    "\n",
    "<img src=\"images/l11_i4.23.png\" alt=\"bistability limit cycle\" style=\"width: 400px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusions\n",
    "\n",
    "- Dynamical systems theory is the study of general principles of how systems of differential equations behave\n",
    "- Geometric analysis analyzes such systems in terms of phase portraits\n",
    "- Local equivalence allows us to understand how nonlinear dynamical systems behave around equilibrium points\n",
    "    - In 1-D systems, the slope of the system equation $F(V)$ determines whether equilibria are stable or not\n",
    "    - In 2-D systems, the \"slope\" is a matrix called the Jacobian\n",
    "- Changes to parameters of a system can cause qualitative shifts in behavior called **bifurcations**\n",
    "\n",
    "## Next week\n",
    "\n",
    "- 2-D bifurcations\n",
    "- Reduced dynamical neuron models"
   ]
  }
 ],
 "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
}