{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "view-in-github"
},
"source": [
"<a href=\"https://colab.research.google.com/github/EugOT/compneuro/blob/main/notebooks/2023-12-03-DynamicNetworks-1.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a> <a href=\"https://kaggle.com/kernels/welcome?src=https://raw.githubusercontent.com/EugOT/compneuro/main/notebooks/2023-11-26-DynamicNetworks-1.ipynb\" target=\"_parent\"><img src=\"https://kaggle.com/static/images/open-in-kaggle.svg\" alt=\"Open in Kaggle\"/></a>"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "3QU6bQ9ylV4f"
},
"source": [
"# Tutorial 1: Neural Rate Models"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "Owt9JCmHlV4g"
},
"source": [
"---\n",
"# Tutorial Objectives\n",
"\n",
"The brain is a complex system, not because it is composed of a large number of diverse types of neurons, but mainly because of how neurons are connected to each other. The brain is indeed a network of highly specialized neuronal networks.\n",
"\n",
"The activity of a neural network constantly evolves in time. For this reason, neurons can be modeled as dynamical systems. The dynamical system approach is only one of the many modeling approaches that computational neuroscientists have developed (other points of view include information processing, statistical models, etc.).\n",
"\n",
"How the dynamics of neuronal networks affect the representation and processing of information in the brain is an open question. However, signatures of altered brain dynamics present in many brain diseases (e.g., in epilepsy or Parkinson's disease) tell us that it is crucial to study network activity dynamics if we want to understand the brain.\n",
"\n",
"In this tutorial, we will simulate and study one of the simplest models of biological neuronal networks. Instead of modeling and simulating individual excitatory neurons (e.g., LIF models that you implemented yesterday), we will treat them as a single homogeneous population and approximate their dynamics using a single one-dimensional equation describing the evolution of their average spiking rate in time.\n",
"\n",
"In this tutorial, we will learn how to build a firing rate model of a single population of excitatory neurons.\n",
"\n",
"<br>\n",
"\n",
"**Steps:**\n",
"- Write the equation for the firing rate dynamics of a 1D excitatory population.\n",
"- Visualize the response of the population as a function of parameters such as threshold level and gain, using the frequency-current (F-I) curve.\n",
"- Numerically simulate the dynamics of the excitatory population and find the fixed points of the system."
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "INdxjlYUlV4h"
},
"source": [
"---\n",
"# Setup"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"cellView": "both",
"execution": {},
"id": "HKzU_g8elV4h"
},
"outputs": [],
"source": [
"# Imports\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import scipy.optimize as opt # root-finding algorithm"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"cellView": "form",
"execution": {},
"id": "SQO7mBfelV4h"
},
"outputs": [],
"source": [
"# @title Figure Settings\n",
"import logging\n",
"logging.getLogger('matplotlib.font_manager').disabled = True\n",
"\n",
"import ipywidgets as widgets # interactive display\n",
"%config InlineBackend.figure_format = 'retina'\n",
"plt.style.use(\"https://raw.githubusercontent.com/NeuromatchAcademy/course-content/main/nma.mplstyle\")"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"cellView": "form",
"execution": {},
"id": "NqNEZiaelV4i"
},
"outputs": [],
"source": [
"# @title Plotting Functions\n",
"\n",
"def plot_fI(x, f):\n",
" plt.figure(figsize=(6, 4)) # plot the figure\n",
" plt.plot(x, f, 'k')\n",
" plt.xlabel('x (a.u.)', fontsize=14)\n",
" plt.ylabel('F(x)', fontsize=14)\n",
" plt.show()\n",
"\n",
"\n",
"def plot_dr_r(r, drdt, x_fps=None):\n",
" plt.figure()\n",
" plt.plot(r, drdt, 'k')\n",
" plt.plot(r, 0. * r, 'k--')\n",
" if x_fps is not None:\n",
" plt.plot(x_fps, np.zeros_like(x_fps), \"ko\", ms=12)\n",
" plt.xlabel(r'$r$')\n",
" plt.ylabel(r'$\\frac{dr}{dt}$', fontsize=20)\n",
" plt.ylim(-0.1, 0.1)\n",
" plt.show()\n",
"\n",
"\n",
"def plot_dFdt(x, dFdt):\n",
" plt.figure()\n",
" plt.plot(x, dFdt, 'r')\n",
" plt.xlabel('x (a.u.)', fontsize=14)\n",
" plt.ylabel('dF(x)', fontsize=14)\n",
" plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "mKdQuiI3lV4i"
},
"source": [
"---\n",
"# Section 1: Neuronal network dynamics"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "ACgrYXENlV4i"
},
"source": [
"## Section 1.1: Dynamics of a single excitatory population\n",
"\n",
"\n",
"<details>\n",
"<summary> <font color='blue'>Click here for text recap of relevant theory </font></summary>\n",
"\n",
"Individual neurons respond by spiking. When we average the spikes of neurons in a population, we can define the average firing activity of the population. In this model, we are interested in how the population-averaged firing varies as a function of time and network parameters. Mathematically, we can describe the firing rate dynamic of a feed-forward network as:\n",
"\n",
"\\begin{equation}\n",
"\\tau \\frac{dr}{dt} = -r + F(I_{\\text{ext}}) \\quad\\qquad (1)\n",
"\\end{equation}\n",
"\n",
"$r(t)$ represents the average firing rate of the excitatory population at time $t$, $\\tau$ controls the timescale of the evolution of the average firing rate, $I_{\\text{ext}}$ represents the external input, and the transfer function $F(\\cdot)$ (which can be related to f-I curve of individual neurons described in the next sections) represents the population activation function in response to all received inputs.\n",
"\n",
"</details>\n",
"\n",
"To start building the model, please execute the cell below to initialize the simulation parameters."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"cellView": "form",
"execution": {},
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "Y8bo260ClV4j",
"outputId": "14545790-feeb-4ed7-9860-2a023d4d0950"
},
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"{'tau': 1.0, 'a': 1.2, 'theta': 2.8, 'w': 0.0, 'I_ext': 0.0, 'T': 20.0, 'dt': 0.1, 'r_init': 0.2, 'range_t': array([ 0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ,\n",
" 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2. , 2.1,\n",
" 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3. , 3.1, 3.2,\n",
" 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4. , 4.1, 4.2, 4.3,\n",
" 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5. , 5.1, 5.2, 5.3, 5.4,\n",
" 5.5, 5.6, 5.7, 5.8, 5.9, 6. , 6.1, 6.2, 6.3, 6.4, 6.5,\n",
" 6.6, 6.7, 6.8, 6.9, 7. , 7.1, 7.2, 7.3, 7.4, 7.5, 7.6,\n",
" 7.7, 7.8, 7.9, 8. , 8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7,\n",
" 8.8, 8.9, 9. , 9.1, 9.2, 9.3, 9.4, 9.5, 9.6, 9.7, 9.8,\n",
" 9.9, 10. , 10.1, 10.2, 10.3, 10.4, 10.5, 10.6, 10.7, 10.8, 10.9,\n",
" 11. , 11.1, 11.2, 11.3, 11.4, 11.5, 11.6, 11.7, 11.8, 11.9, 12. ,\n",
" 12.1, 12.2, 12.3, 12.4, 12.5, 12.6, 12.7, 12.8, 12.9, 13. , 13.1,\n",
" 13.2, 13.3, 13.4, 13.5, 13.6, 13.7, 13.8, 13.9, 14. , 14.1, 14.2,\n",
" 14.3, 14.4, 14.5, 14.6, 14.7, 14.8, 14.9, 15. , 15.1, 15.2, 15.3,\n",
" 15.4, 15.5, 15.6, 15.7, 15.8, 15.9, 16. , 16.1, 16.2, 16.3, 16.4,\n",
" 16.5, 16.6, 16.7, 16.8, 16.9, 17. , 17.1, 17.2, 17.3, 17.4, 17.5,\n",
" 17.6, 17.7, 17.8, 17.9, 18. , 18.1, 18.2, 18.3, 18.4, 18.5, 18.6,\n",
" 18.7, 18.8, 18.9, 19. , 19.1, 19.2, 19.3, 19.4, 19.5, 19.6, 19.7,\n",
" 19.8, 19.9])}\n"
]
}
],
"source": [
"# @markdown *Execute this cell to set default parameters for a single excitatory population model*\n",
"\n",
"\n",
"def default_pars_single(**kwargs):\n",
" pars = {}\n",
"\n",
" # Excitatory parameters\n",
" pars['tau'] = 1. # Timescale of the E population [ms]\n",
" pars['a'] = 1.2 # Gain of the E population\n",
" pars['theta'] = 2.8 # Threshold of the E population\n",
"\n",
" # Connection strength\n",
" pars['w'] = 0. # E to E, we first set it to 0\n",
"\n",
" # External input\n",
" pars['I_ext'] = 0.\n",
"\n",
" # simulation parameters\n",
" pars['T'] = 20. # Total duration of simulation [ms]\n",
" pars['dt'] = .1 # Simulation time step [ms]\n",
" pars['r_init'] = 0.2 # Initial value of E\n",
"\n",
" # External parameters if any\n",
" pars.update(kwargs)\n",
"\n",
" # Vector of discretized time points [ms]\n",
" pars['range_t'] = np.arange(0, pars['T'], pars['dt'])\n",
"\n",
" return pars\n",
"\n",
"pars = default_pars_single()\n",
"print(pars)"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "z93VM81DlV4j"
},
"source": [
"You can now use:\n",
"- `pars = default_pars_single()` to get all the parameters.\n",
"- `pars = default_pars_single(T=T_sim, dt=time_step)` to set new simulation time and time step\n",
"- To update an existing parameter dictionary, use `pars['New_para'] = value`\n",
"\n",
"Because `pars` is a dictionary, it can be passed to a function that requires individual parameters as arguments using `my_func(**pars)` syntax."
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "J4A-wDvqlV4j"
},
"source": [
"## Section 1.2: F-I curves\n",
"\n",
"\n",
"<details>\n",
"<summary> <font color='blue'>Click here for text recap of relevant theory </font></summary>\n",
"\n",
"In electrophysiology, a neuron is often characterized by its spike rate output in response to input currents. This is often called the **F-I** curve, denoting the output spike frequency (**F**) in response to different injected currents (**I**). We estimated this for an LIF neuron in yesterday's tutorial.\n",
"\n",
"The transfer function $F(\\cdot)$ in Equation $1$ represents the gain of the population as a function of the total input. The gain is often modeled as a sigmoidal function, i.e., more input drive leads to a nonlinear increase in the population firing rate. The output firing rate will eventually saturate for high input values.\n",
"\n",
"A sigmoidal $F(\\cdot)$ is parameterized by its gain $a$ and threshold $\\theta$.\n",
"\n",
"$$ F(x;a,\\theta) = \\frac{1}{1+\\text{e}^{-a(x-\\theta)}} - \\frac{1}{1+\\text{e}^{a\\theta}} \\quad(2)$$\n",
"\n",
"The argument $x$ represents the input to the population. Note that the second term is chosen so that $F(0;a,\\theta)=0$.\n",
"\n",
"Many other transfer functions (generally monotonic) can be also used. Examples are the rectified linear function $ReLU(x)$ or the hyperbolic tangent $tanh(x)$."
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "RtblIbojlV4j"
},
"source": [
"### Coding Exercise 1.2: Implement F-I curve\n",
"\n",
"Let's first investigate the activation functions before simulating the dynamics of the entire population.\n",
"\n",
"In this exercise, you will implement a sigmoidal **F-I** curve or transfer function $F(x)$, with gain $a$ and threshold level $\\theta$ as parameters:\n",
"\n",
"\\begin{equation}\n",
"F(x;a,\\theta) = \\frac{1}{1+\\text{e}^{-a(x-\\theta)}} - \\frac{1}{1+\\text{e}^{a\\theta}}\n",
"\\end{equation}"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"execution": {},
"colab": {
"base_uri": "https://localhost:8080/",
"height": 332
},
"id": "tS0mnhfWlV4j",
"outputId": "09d7eea8-878c-498d-f0aa-9ea1699e6ffa"
},
"outputs": [
{
"output_type": "error",
"ename": "NotImplementedError",
"evalue": "Student exercise: implement the f-I function",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mNotImplementedError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-5-2a9bb0329f56>\u001b[0m in \u001b[0;36m<cell line: 30>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 28\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 29\u001b[0m \u001b[0;31m# Compute transfer function\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 30\u001b[0;31m \u001b[0mf\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mF\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpars\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'a'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpars\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'theta'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 31\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 32\u001b[0m \u001b[0;31m# Visualize\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m<ipython-input-5-2a9bb0329f56>\u001b[0m in \u001b[0;36mF\u001b[0;34m(x, a, theta)\u001b[0m\n\u001b[1;32m 14\u001b[0m \u001b[0;31m## TODO for students: compute f = F(x) ##\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 15\u001b[0m \u001b[0;31m# Fill out function and remove\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 16\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mNotImplementedError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Student exercise: implement the f-I function\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 17\u001b[0m \u001b[0;31m#################################################\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 18\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mNotImplementedError\u001b[0m: Student exercise: implement the f-I function"
]
}
],
"source": [
"def F(x, a, theta):\n",
" \"\"\"\n",
" Population activation function.\n",
"\n",
" Args:\n",
" x (float): the population input\n",
" a (float): the gain of the function\n",
" theta (float): the threshold of the function\n",
"\n",
" Returns:\n",
" float: the population activation response F(x) for input x\n",
" \"\"\"\n",
" #################################################\n",
" ## TODO for students: compute f = F(x) ##\n",
" # Fill out function and remove\n",
" raise NotImplementedError(\"Student exercise: implement the f-I function\")\n",
" #################################################\n",
"\n",
" # Define the sigmoidal transfer function f = F(x)\n",
" f = ...\n",
"\n",
" return f\n",
"\n",
"\n",
"# Set parameters\n",
"pars = default_pars_single() # get default parameters\n",
"x = np.arange(0, 10, .1) # set the range of input\n",
"\n",
"# Compute transfer function\n",
"f = F(x, pars['a'], pars['theta'])\n",
"\n",
"# Visualize\n",
"plot_fI(x, f)"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "5Xog0w7DlV4k"
},
"source": [
"### Interactive Demo 1.2 : Parameter exploration of F-I curve\n",
"\n",
"Here's an interactive demo that shows how the F-I curve changes for different values of the gain and threshold parameters.\n",
"\n",
"1. How does the gain parameter ($a$) affect the F-I curve?\n",
"2. How does the threshold parameter ($\\theta$) affect the F-I curve?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {},
"id": "UV5O9raUlV4k"
},
"outputs": [],
"source": [
"# @markdown Make sure you execute this cell to enable the widget!\n",
"\n",
"@widgets.interact(\n",
" a=widgets.FloatSlider(1.5, min=0.3, max=3., step=0.3),\n",
" theta=widgets.FloatSlider(3., min=2., max=4., step=0.2)\n",
")\n",
"\n",
"\n",
"def interactive_plot_FI(a, theta):\n",
" \"\"\"\n",
" Population activation function.\n",
"\n",
" Expecxts:\n",
" a : the gain of the function\n",
" theta : the threshold of the function\n",
"\n",
" Returns:\n",
" plot the F-I curve with give parameters\n",
" \"\"\"\n",
"\n",
" # set the range of input\n",
" x = np.arange(0, 10, .1)\n",
" plt.figure()\n",
" plt.plot(x, F(x, a, theta), 'k')\n",
" plt.xlabel('x (a.u.)', fontsize=14)\n",
" plt.ylabel('F(x)', fontsize=14)\n",
" plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "tweFbIiMlV4k"
},
"source": [
"## Section 1.3: Simulation scheme of E dynamics\n",
"\n",
"Because $F(\\cdot)$ is a nonlinear function, the exact solution of our differential equation of population activity can not be determined via analytical methods. As we have seen before, we can use numerical methods, specifically the Euler method, to find the solution (that is, simulate the population activity)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {},
"id": "_xTqFRIElV4k"
},
"outputs": [],
"source": [
"# @markdown Execute this cell to enable the single population rate model simulator: `simulate_single`\n",
"\n",
"\n",
"def simulate_single(pars):\n",
" \"\"\"\n",
" Simulate an excitatory population of neurons\n",
"\n",
" Args:\n",
" pars : Parameter dictionary\n",
"\n",
" Returns:\n",
" rE : Activity of excitatory population (array)\n",
"\n",
" Example:\n",
" pars = default_pars_single()\n",
" r = simulate_single(pars)\n",
" \"\"\"\n",
"\n",
" # Set parameters\n",
" tau, a, theta = pars['tau'], pars['a'], pars['theta']\n",
" w = pars['w']\n",
" I_ext = pars['I_ext']\n",
" r_init = pars['r_init']\n",
" dt, range_t = pars['dt'], pars['range_t']\n",
" Lt = range_t.size\n",
"\n",
" # Initialize activity\n",
" r = np.zeros(Lt)\n",
" r[0] = r_init\n",
" I_ext = I_ext * np.ones(Lt)\n",
"\n",
" # Update the E activity\n",
" for k in range(Lt - 1):\n",
" dr = dt / tau * (-r[k] + F(w * r[k] + I_ext[k], a, theta))\n",
" r[k+1] = r[k] + dr\n",
"\n",
" return r\n",
"\n",
"\n",
"help(simulate_single)"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "SVBuawFylV4k"
},
"source": [
"### Interactive Demo 1.3: Parameter Exploration of single population dynamics\n",
"\n",
"Explore these dynamics of the population activity in this interactive demo.\n",
"\n",
"1. How does $r_{\\text{sim}}(t)$ change with different $I_{\\text{ext}}$ values?\n",
"2. How does it change with different $\\tau$ values?\n",
"\n",
"<br>\n",
"\n",
"Note that, $r_{\\rm ana}(t)$ denotes the analytical solution - you will learn how this is computed in the next section."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {},
"id": "Ibo659HrlV4k"
},
"outputs": [],
"source": [
"# @title\n",
"\n",
"# @markdown Make sure you execute this cell to enable the widget!\n",
"\n",
"# get default parameters\n",
"pars = default_pars_single(T=20.)\n",
"\n",
"@widgets.interact(\n",
" I_ext=widgets.FloatSlider(5.0, min=0.0, max=10., step=1.),\n",
" tau=widgets.FloatSlider(3., min=1., max=5., step=0.2)\n",
")\n",
"\n",
"\n",
"def Myplot_E_diffI_difftau(I_ext, tau):\n",
" # set external input and time constant\n",
" pars['I_ext'] = I_ext\n",
" pars['tau'] = tau\n",
"\n",
" # simulation\n",
" r = simulate_single(pars)\n",
"\n",
" # Analytical Solution\n",
" r_ana = (pars['r_init']\n",
" + (F(I_ext, pars['a'], pars['theta'])\n",
" - pars['r_init']) * (1. - np.exp(-pars['range_t'] / pars['tau'])))\n",
"\n",
" # plot\n",
" plt.figure()\n",
" plt.plot(pars['range_t'], r, 'b', label=r'$r_{\\mathrm{sim}}$(t)', alpha=0.5,\n",
" zorder=1)\n",
" plt.plot(pars['range_t'], r_ana, 'b--', lw=5, dashes=(2, 2),\n",
" label=r'$r_{\\mathrm{ana}}$(t)', zorder=2)\n",
" plt.plot(pars['range_t'],\n",
" F(I_ext, pars['a'], pars['theta']) * np.ones(pars['range_t'].size),\n",
" 'k--', label=r'$F(I_{\\mathrm{ext}})$')\n",
" plt.xlabel('t (ms)', fontsize=16.)\n",
" plt.ylabel('Activity r(t)', fontsize=16.)\n",
" plt.legend(loc='best', fontsize=14.)\n",
" plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "kBGeRBvJlV4k"
},
"source": [
"## Think! 1.3: Finite activities\n",
"\n",
"Above, we have numerically solved a system driven by a positive input. Yet, $r_E(t)$ either decays to zero or reaches a fixed non-zero value.\n",
"\n",
"1. Why doesn't the solution of the system \"explode\" in a finite time? In other words, what guarantees that $r_E$(t) stays finite?\n",
"2. Which parameter would you change in order to increase the maximum value of the response?"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "yHww4G8UlV4k"
},
"source": [
"---\n",
"# Section 2: Fixed points of the single population system\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "BNbe9S7ElV4k"
},
"source": [
"## Section 2.1: Finding fixed points"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "kxIO_jW_lV4k"
},
"source": [
"This part introduces recurrent networks and how to derive their fixed points. It also introduces vector fields and phase planes in one dimension.\n",
"\n",
"<details>\n",
"<summary> <font color='blue'>Click here for text recap theory </font></summary>\n",
"\n",
"We can now extend our feed-forward network to a recurrent network, governed by the equation:\n",
"\n",
"\\begin{align}\n",
"\\tau \\frac{dr}{dt} &= -r + F(w\\cdot r + I_{\\text{ext}}) \\quad\\qquad (3)\n",
"\\end{align}\n",
"\n",
"where as before, $r(t)$ represents the average firing rate of the excitatory population at time $t$, $\\tau$ controls the timescale of the evolution of the average firing rate, $I_{\\text{ext}}$ represents the external input, and the transfer function $F(\\cdot)$ (which can be related to f-I curve of individual neurons described in the next sections) represents the population activation function in response to all received inputs. Now we also have $w$ which denotes the strength (synaptic weight) of the recurrent input to the population.\n",
"\n",
"As you varied the two parameters in the last Interactive Demo, you noticed that, while at first the system output quickly changes, with time, it reaches its maximum/minimum value and does not change anymore. The value eventually reached by the system is called the **steady state** of the system, or the **fixed point**. Essentially, in the steady states the derivative with respect to time of the activity ($r$) is zero, i.e. $\\displaystyle \\frac{dr}{dt}=0$.\n",
"\n",
"We can find that the steady state of the Equation. (1) by setting $\\displaystyle{\\frac{dr}{dt}=0}$ and solve for $r$:\n",
"\n",
"\\begin{equation}\n",
"-r_{\\text{steady}} + F(w\\cdot r_{\\text{steady}} + I_{\\text{ext}};a,\\theta) = 0 \\qquad (4)\n",
"\\end{equation}\n",
"\n",
"When it exists, the solution of Equation. (4) defines a **fixed point** of the dynamical system in Equation (3). Note that if $F(x)$ is nonlinear, it is not always possible to find an analytical solution, but the solution can be found via numerical simulations, as we will do later.\n",
"\n",
"From the Interactive Demo, one could also notice that the value of $\\tau$ influences how quickly the activity will converge to the steady state from its initial value.\n",
"\n",
"In the specific case of $w=0$, we can also analytically compute the solution of Equation (1) (i.e., the thick blue dashed line) and deduce the role of $\\tau$ in determining the convergence to the fixed point:\n",
"\n",
"\\begin{equation}\n",
"\\displaystyle{r(t) = \\big{[}F(I_{\\text{ext}};a,\\theta) -r(t=0)\\big{]} (1-\\text{e}^{-\\frac{t}{\\tau}})} + r(t=0)\n",
"\\end{equation}\n",
"\n",
"<br>\n",
"\n",
"We can now numerically calculate the fixed point with a root finding algorithm."
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "pOdpzLjWlV4k"
},
"source": [
"### Coding Exercise 2.1.1: Visualization of the fixed points\n",
"\n",
"When it is not possible to find the solution for Equation (3) analytically, a graphical approach can be taken. To that end, it is useful to plot $\\displaystyle{\\frac{dr}{dt}}$ as a function of $r$. The values of $r$ for which the plotted function crosses zero on the y axis correspond to fixed points.\n",
"\n",
"Here, let us, for example, set $w=5.0$ and $I^{\\text{ext}}=0.5$. From Equation (3), you can obtain\n",
"\n",
"\\begin{equation}\n",
"\\frac{dr}{dt} = [-r + F(w\\cdot r + I^{\\text{ext}})]\\,/\\,\\tau\n",
"\\end{equation}\n",
"\n",
"Then, plot the $dr/dt$ as a function of $r$, and check for the presence of fixed points."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"id": "wcTbINshlV4k"
},
"outputs": [],
"source": [
"def compute_drdt(r, I_ext, w, a, theta, tau, **other_pars):\n",
" \"\"\"Given parameters, compute dr/dt as a function of r.\n",
"\n",
" Args:\n",
" r (1D array) : Average firing rate of the excitatory population\n",
" I_ext, w, a, theta, tau (numbers): Simulation parameters to use\n",
" other_pars : Other simulation parameters are unused by this function\n",
"\n",
" Returns\n",
" drdt function for each value of r\n",
" \"\"\"\n",
" #########################################################################\n",
" # TODO compute drdt and disable the error\n",
" raise NotImplementedError(\"Finish the compute_drdt function\")\n",
" #########################################################################\n",
"\n",
" # Calculate drdt\n",
" drdt = ...\n",
"\n",
" return drdt\n",
"\n",
"\n",
"# Define a vector of r values and the simulation parameters\n",
"r = np.linspace(0, 1, 1000)\n",
"pars = default_pars_single(I_ext=0.5, w=5)\n",
"\n",
"# Compute dr/dt\n",
"drdt = compute_drdt(r, **pars)\n",
"\n",
"# Visualize\n",
"plot_dr_r(r, drdt)"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "l4najG0SlV4k"
},
"source": [
"### Coding Exercise 2.1.2: Numerical calculation of fixed points\n",
"\n",
"We will now find the fixed points numerically. To do so, we need to specify initial values ($r_{\\text{guess}}$) for the root-finding algorithm to start from. From the line $\\displaystyle{\\frac{dr}{dt}}$ plotted above in the last exercise, initial values can be chosen as a set of values close to where the line crosses zero on the y axis (real fixed point).\n",
"\n",
"The next cell defines three helper functions that we will use:\n",
"\n",
"- `my_fp_single(r_guess, **pars)` uses a root-finding algorithm to locate a fixed point near a given initial value\n",
"- `check_fp_single(x_fp, **pars)` verifies that the values of $r_{\\rm fp}$ for which $\\displaystyle{\\frac{dr}{dt}} = 0$ are the true fixed points\n",
"- `my_fp_finder(r_guess_vector, **pars)` accepts an array of initial values and finds the same number of fixed points, using the above two functions"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {},
"id": "hqzq0_lZlV4l"
},
"outputs": [],
"source": [
"# @markdown Execute this cell to enable the fixed point functions\n",
"\n",
"def my_fp_single(r_guess, a, theta, w, I_ext, **other_pars):\n",
" \"\"\"\n",
" Calculate the fixed point through drE/dt=0\n",
"\n",
" Args:\n",
" r_guess : Initial value used for scipy.optimize function\n",
" a, theta, w, I_ext : simulation parameters\n",
"\n",
" Returns:\n",
" x_fp : value of fixed point\n",
" \"\"\"\n",
" # define the right hand of E dynamics\n",
" def my_WCr(x):\n",
" r = x\n",
" drdt = (-r + F(w * r + I_ext, a, theta))\n",
" y = np.array(drdt)\n",
"\n",
" return y\n",
"\n",
" x0 = np.array(r_guess)\n",
" x_fp = opt.root(my_WCr, x0).x.item()\n",
"\n",
" return x_fp\n",
"\n",
"\n",
"def check_fp_single(x_fp, a, theta, w, I_ext, mytol=1e-4, **other_pars):\n",
" \"\"\"\n",
" Verify |dr/dt| < mytol\n",
"\n",
" Args:\n",
" fp : value of fixed point\n",
" a, theta, w, I_ext: simulation parameters\n",
" mytol : tolerance, default as 10^{-4}\n",
"\n",
" Returns :\n",
" Whether it is a correct fixed point: True/False\n",
" \"\"\"\n",
" # calculate Equation(3)\n",
" y = x_fp - F(w * x_fp + I_ext, a, theta)\n",
"\n",
" # Here we set tolerance as 10^{-4}\n",
" return np.abs(y) < mytol\n",
"\n",
"\n",
"def my_fp_finder(pars, r_guess_vector, mytol=1e-4):\n",
" \"\"\"\n",
" Calculate the fixed point(s) through drE/dt=0\n",
"\n",
" Args:\n",
" pars : Parameter dictionary\n",
" r_guess_vector : Initial values used for scipy.optimize function\n",
" mytol : tolerance for checking fixed point, default as 10^{-4}\n",
"\n",
" Returns:\n",
" x_fps : values of fixed points\n",
"\n",
" \"\"\"\n",
" x_fps = []\n",
" correct_fps = []\n",
" for r_guess in r_guess_vector:\n",
" x_fp = my_fp_single(r_guess, **pars)\n",
" if check_fp_single(x_fp, **pars, mytol=mytol):\n",
" x_fps.append(x_fp)\n",
"\n",
" return x_fps\n",
"\n",
"\n",
"help(my_fp_finder)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"id": "6s7M5vrFlV4l"
},
"outputs": [],
"source": [
"# Set parameters\n",
"r = np.linspace(0, 1, 1000)\n",
"pars = default_pars_single(I_ext=0.5, w=5)\n",
"\n",
"# Compute dr/dt\n",
"drdt = compute_drdt(r, **pars)\n",
"\n",
"#############################################################################\n",
"# TODO for students:\n",
"# Define initial values close to the intersections of drdt and y=0\n",
"# (How many initial values? Hint: How many times do the two lines intersect?)\n",
"# Calculate the fixed point with these initial values and plot them\n",
"raise NotImplementedError('student_exercise: find fixed points numerically')\n",
"#############################################################################\n",
"\n",
"# Initial guesses for fixed points\n",
"r_guess_vector = [...]\n",
"\n",
"# Find fixed point numerically\n",
"x_fps = my_fp_finder(pars, r_guess_vector)\n",
"\n",
"# Visualize\n",
"plot_dr_r(r, drdt, x_fps)"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "2NrhupwIlV4l"
},
"source": [
"### Interactive Demo 2.1: fixed points as a function of recurrent and external inputs.\n",
"\n",
"You can now explore how the previous plot changes when the recurrent coupling $w$ and the external input $I_{\\text{ext}}$ take different values. How does the number of fixed points change?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {},
"id": "p-zd9U-MlV4l"
},
"outputs": [],
"source": [
"# @markdown Make sure you execute this cell to enable the widget!\n",
"\n",
"@widgets.interact(\n",
" w=widgets.FloatSlider(4., min=1., max=7., step=0.2),\n",
" I_ext=widgets.FloatSlider(1.5, min=0., max=3., step=0.1)\n",
")\n",
"\n",
"\n",
"def plot_intersection_single(w, I_ext):\n",
" # set your parameters\n",
" pars = default_pars_single(w=w, I_ext=I_ext)\n",
"\n",
" # find fixed points\n",
" r_init_vector = [0, .4, .9]\n",
" x_fps = my_fp_finder(pars, r_init_vector)\n",
"\n",
" # plot\n",
" r = np.linspace(0, 1., 1000)\n",
" drdt = (-r + F(w * r + I_ext, pars['a'], pars['theta'])) / pars['tau']\n",
"\n",
" plot_dr_r(r, drdt, x_fps)"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "7ckt1wwGlV4l"
},
"source": [
"## Section 2.2: Relationship between trajectories & fixed points\n",
"\n",
"Let's examine the relationship between the population activity over time and the fixed points.\n",
"\n",
"Here, let us first set $w=5.0$ and $I_{\\text{ext}}=0.5$, and investigate the dynamics of $r(t)$ starting with different initial values $r(0) \\equiv r_{\\text{init}}$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {},
"id": "I2-ivVtnlV4l"
},
"outputs": [],
"source": [
"# @markdown Execute to visualize dr/dt\n",
"\n",
"def plot_intersection_single(w, I_ext):\n",
" # set your parameters\n",
" pars = default_pars_single(w=w, I_ext=I_ext)\n",
"\n",
" # find fixed points\n",
" r_init_vector = [0, .4, .9]\n",
" x_fps = my_fp_finder(pars, r_init_vector)\n",
"\n",
" # plot\n",
" r = np.linspace(0, 1., 1000)\n",
" drdt = (-r + F(w * r + I_ext, pars['a'], pars['theta'])) / pars['tau']\n",
"\n",
" plot_dr_r(r, drdt, x_fps)\n",
"\n",
"\n",
"plot_intersection_single(w = 5.0, I_ext = 0.5)"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "0I-SaIAElV4l"
},
"source": [
"### Interactive Demo 2.2: dynamics as a function of the initial value\n",
"\n",
"Let's now set $r_{\\rm init}$ to a value of your choice in this demo. How does the solution change? What do you observe? How does that relate to the previous plot of $\\frac{dr}{dt}$?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {},
"id": "d2R9ptPFlV4l"
},
"outputs": [],
"source": [
"# @markdown Make sure you execute this cell to enable the widget!\n",
"pars = default_pars_single(w=5.0, I_ext=0.5)\n",
"\n",
"@widgets.interact(\n",
" r_init=widgets.FloatSlider(0.5, min=0., max=1., step=.02)\n",
")\n",
"\n",
"def plot_single_diffEinit(r_init):\n",
" pars['r_init'] = r_init\n",
" r = simulate_single(pars)\n",
"\n",
" plt.figure()\n",
" plt.plot(pars['range_t'], r, 'b', zorder=1)\n",
" plt.plot(0, r[0], 'bo', alpha=0.7, zorder=2)\n",
" plt.xlabel('t (ms)', fontsize=16)\n",
" plt.ylabel(r'$r(t)$', fontsize=16)\n",
" plt.ylim(0, 1.0)\n",
" plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "Eb-UYoCTlV4l"
},
"source": [
"We will plot the trajectories of $r(t)$ with $r_{\\text{init}} = 0.0, 0.1, 0.2,..., 0.9$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {},
"id": "RAYbVBFVlV4l"
},
"outputs": [],
"source": [
"# @markdown Execute this cell to see the trajectories!\n",
"\n",
"pars = default_pars_single()\n",
"pars['w'] = 5.0\n",
"pars['I_ext'] = 0.5\n",
"\n",
"plt.figure(figsize=(8, 5))\n",
"for ie in range(10):\n",
" pars['r_init'] = 0.1 * ie # set the initial value\n",
" r = simulate_single(pars) # run the simulation\n",
"\n",
" # plot the activity with given initial\n",
" plt.plot(pars['range_t'], r, 'b', alpha=0.1 + 0.1 * ie,\n",
" label=r'r$_{\\mathrm{init}}$=%.1f' % (0.1 * ie))\n",
"\n",
"plt.xlabel('t (ms)')\n",
"plt.title('Two steady states?')\n",
"plt.ylabel(r'$r$(t)')\n",
"plt.legend(loc=[1.01, -0.06], fontsize=14)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "N_rrkoxnlV4l"
},
"source": [
"We have three fixed points but only two steady states showing up - what's happening?\n",
"\n",
"It turns out that the stability of the fixed points matters. If a fixed point is stable, a trajectory starting near that fixed point will stay close to that fixed point and converge to it (the steady state will equal the fixed point). If a fixed point is unstable, any trajectories starting close to it will diverge and go towards stable fixed points. In fact, the only way for a trajectory to reach a stable state at an unstable fixed point is if the initial value **exactly** equals the value of the fixed point."
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "Bpd3GcOilV4m"
},
"source": [
"### Think! 2.2: Stable vs unstable fixed points\n",
"\n",
"Which of the fixed points for the model we've been examining in this section are stable vs unstable?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {},
"id": "CIhNDYNnlV4p"
},
"outputs": [],
"source": [
"# @markdown Execute to print fixed point values\n",
"\n",
"# Initial guesses for fixed points\n",
"r_guess_vector = [0, .4, .9]\n",
"\n",
"# Find fixed point numerically\n",
"x_fps = my_fp_finder(pars, r_guess_vector)\n",
"\n",
"print(f'Our fixed points are {x_fps}')"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "ZpxD-76ulV4p"
},
"source": [
"We can simulate the trajectory if we start at the unstable fixed point: you can see that it remains at that fixed point (the red line below)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {},
"id": "kcKQfSnZlV4p"
},
"outputs": [],
"source": [
"# @markdown Execute to visualize trajectory starting at unstable fixed point\n",
"pars = default_pars_single()\n",
"pars['w'] = 5.0\n",
"pars['I_ext'] = 0.5\n",
"\n",
"plt.figure(figsize=(8, 5))\n",
"for ie in range(10):\n",
" pars['r_init'] = 0.1 * ie # set the initial value\n",
" r = simulate_single(pars) # run the simulation\n",
"\n",
" # plot the activity with given initial\n",
" plt.plot(pars['range_t'], r, 'b', alpha=0.1 + 0.1 * ie,\n",
" label=r'r$_{\\mathrm{init}}$=%.1f' % (0.1 * ie))\n",
"\n",
"pars['r_init'] = x_fps[1] # set the initial value\n",
"r = simulate_single(pars) # run the simulation\n",
"\n",
"# plot the activity with given initial\n",
"plt.plot(pars['range_t'], r, 'r', alpha=0.1 + 0.1 * ie,\n",
" label=r'r$_{\\mathrm{init}}$=%.4f' % (x_fps[1]))\n",
"\n",
"plt.xlabel('t (ms)')\n",
"plt.title('Two steady states?')\n",
"plt.ylabel(r'$r$(t)')\n",
"plt.legend(loc=[1.01, -0.06], fontsize=14)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "y07z9j-_lV4p"
},
"source": [
"See Bonus Section 1 to cover how to determine the stability of fixed points in a quantitative way."
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "XEDGyVLelV4p"
},
"source": [
"## Think! 2: Inhibitory populations\n",
"\n",
"Throughout the tutorial, we have assumed $w> 0 $, i.e., we considered a single population of **excitatory** neurons. What do you think will be the behavior of a population of inhibitory neurons, i.e., where $w> 0$ is replaced by $w< 0$?"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "FX76yXVQlV4p"
},
"source": [
"---\n",
"# Summary\n",
"\n",
"In this tutorial, we have investigated the dynamics of a rate-based single population of neurons.\n",
"\n",
"We learned about:\n",
"- The effect of the input parameters and the time constant of the network on the dynamics of the population.\n",
"- How to find the fixed point(s) of the system.\n",
"\n",
"We build on these concepts in the bonus material - check it out if you have time. You will learn:\n",
"\n",
"- How to determine the stability of a fixed point by linearizing the system.\n",
"- How to add realistic inputs to our model."
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "5xi4Ty6ZlV4p"
},
"source": [
"---\n",
"# Bonus"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "uEJZtsxllV4p"
},
"source": [
"## Bonus Section 1: Stability analysis via linearization of the dynamics"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "IwQ-3cErlV4p"
},
"source": [
"Here we will dive into the math of how to figure out the stability of a fixed point.\n",
"\n",
"\n",
"Just like in our equation for the feedforward network, a generic linear system\n",
"$\\frac{dx}{dt} = \\lambda (x - b)$, has a fixed point for $x=b$. The analytical solution of such a system can be found to be:\n",
"\n",
"\\begin{equation}\n",
"x(t) = b + \\big{(} x(0) - b \\big{)} \\text{e}^{\\lambda t}.\n",
"\\end{equation}\n",
"\n",
"Now consider a small perturbation of the activity around the fixed point: $x(0) = b + \\epsilon$, where $|\\epsilon| \\ll 1$. Will the perturbation $\\epsilon(t)$ grow with time or will it decay to the fixed point? The evolution of the perturbation with time can be written, using the analytical solution for $x(t)$, as:\n",
"\n",
"\\begin{equation}\n",
"\\epsilon (t) = x(t) - b = \\epsilon \\text{e}^{\\lambda t}\n",
"\\end{equation}\n",
"\n",
"<br>\n",
"\n",
"- if $\\lambda < 0$, $\\epsilon(t)$ decays to zero, $x(t)$ will still converge to $b$ and the fixed point is \"**stable**\".\n",
"\n",
"- if $\\lambda > 0$, $\\epsilon(t)$ grows with time, $x(t)$ will leave the fixed point $b$ exponentially, and the fixed point is, therefore, \"**unstable**\" ."
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "TlMgB1nwlV4p"
},
"source": [
"Similar to what we did in the linear system above, in order to determine the stability of a fixed point $r^{*}$ of the excitatory population dynamics, we perturb Equation (1) around $r^{*}$ by $\\epsilon$, i.e. $r = r^{*} + \\epsilon$. We can plug in Equation (1) and obtain the equation determining the time evolution of the perturbation $\\epsilon(t)$:\n",
"\n",
"\\begin{equation}\n",
"\\tau \\frac{d\\epsilon}{dt} \\approx -\\epsilon + w F'(w\\cdot r^{*} + I_{\\text{ext}};a,\\theta) \\epsilon\n",
"\\end{equation}\n",
"\n",
"where $F'(\\cdot)$ is the derivative of the transfer function $F(\\cdot)$. We can rewrite the above equation as:\n",
"\n",
"\\begin{equation}\n",
"\\frac{d\\epsilon}{dt} \\approx \\frac{\\epsilon}{\\tau }[-1 + w F'(w\\cdot r^* + I_{\\text{ext}};a,\\theta)]\n",
"\\end{equation}\n",
"\n",
"That is, as in the linear system above, the value of\n",
"\n",
"\\begin{equation}\n",
"\\lambda = [-1+ wF'(w\\cdot r^* + I_{\\text{ext}};a,\\theta)]/\\tau \\qquad (4)\n",
"\\end{equation}\n",
"\n",
"determines whether the perturbation will grow or decay to zero, i.e., $\\lambda$ defines the stability of the fixed point. This value is called the **eigenvalue** of the dynamical system."
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "gHL6o4L4lV4p"
},
"source": [
"The derivative of the sigmoid transfer function is:\n",
"\n",
"\\begin{align}\n",
"\\frac{dF}{dx} & = \\frac{d}{dx} (1+\\exp\\{-a(x-\\theta)\\})^{-1} \\\\\n",
"& = a\\exp\\{-a(x-\\theta)\\} (1+\\exp\\{-a(x-\\theta)\\})^{-2}. \\qquad (5)\n",
"\\end{align}\n",
"\n",
"We provide a helper function `dF` which computes this derivative."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {},
"id": "PPqhkiLilV4p"
},
"outputs": [],
"source": [
"# @markdown Execute this cell to enable helper function `dF` and visualize derivative\n",
"\n",
"def dF(x, a, theta):\n",
" \"\"\"\n",
" Population activation function.\n",
"\n",
" Args:\n",
" x : the population input\n",
" a : the gain of the function\n",
" theta : the threshold of the function\n",
"\n",
" Returns:\n",
" dFdx : the population activation response F(x) for input x\n",
" \"\"\"\n",
"\n",
" # Calculate the population activation\n",
" dFdx = a * np.exp(-a * (x - theta)) * (1 + np.exp(-a * (x - theta)))**-2\n",
"\n",
" return dFdx\n",
"\n",
"# Set parameters\n",
"pars = default_pars_single() # get default parameters\n",
"x = np.arange(0, 10, .1) # set the range of input\n",
"\n",
"# Compute derivative of transfer function\n",
"df = dF(x, pars['a'], pars['theta'])\n",
"\n",
"# Visualize\n",
"plot_dFdt(x, df)"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "uRYMdJHplV4p"
},
"source": [
"### Bonus Coding Exercise 1: Compute eigenvalues\n",
"\n",
"As discussed above, for the case with $w=5.0$ and $I_{\\text{ext}}=0.5$, the system displays **three** fixed points. However, when we simulated the dynamics and varied the initial conditions $r_{\\rm init}$, we could only obtain **two** steady states. In this exercise, we will now check the stability of each of the three fixed points by calculating the corresponding eigenvalues with the function `eig_single`. Check the sign of each eigenvalue (i.e., stability of each fixed point). How many of the fixed points are stable?\n",
"\n",
"Note that the expression of the eigenvalue at fixed point $r^*$:\n",
"\n",
"\\begin{equation}\n",
"\\lambda = [-1+ wF'(w\\cdot r^* + I_{\\text{ext}};a,\\theta)]/\\tau\n",
"\\end{equation}"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {},
"id": "vEerrUR5lV4p"
},
"outputs": [],
"source": [
"def eig_single(fp, tau, a, theta, w, I_ext, **other_pars):\n",
" \"\"\"\n",
" Args:\n",
" fp : fixed point r_fp\n",
" tau, a, theta, w, I_ext : Simulation parameters\n",
"\n",
" Returns:\n",
" eig : eigenvalue of the linearized system\n",
" \"\"\"\n",
" #####################################################################\n",
" ## TODO for students: compute eigenvalue and disable the error\n",
" raise NotImplementedError(\"Student exercise: compute the eigenvalue\")\n",
" ######################################################################\n",
" # Compute the eigenvalue\n",
" eig = ...\n",
"\n",
" return eig\n",
"\n",
"\n",
"# Find the eigenvalues for all fixed points\n",
"pars = default_pars_single(w=5, I_ext=.5)\n",
"r_guess_vector = [0, .4, .9]\n",
"x_fp = my_fp_finder(pars, r_guess_vector)\n",
"\n",
"\n",
"for fp in x_fp:\n",
" eig_fp = eig_single(fp, **pars)\n",
" print(f'Fixed point1 at {fp:.3f} with Eigenvalue={eig_fp:.3f}')"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "h8npq7hGlV4q"
},
"source": [
"**SAMPLE OUTPUT**\n",
"\n",
"```\n",
"Fixed point1 at 0.042 with Eigenvalue=-0.583\n",
"Fixed point2 at 0.447 with Eigenvalue=0.498\n",
"Fixed point3 at 0.900 with Eigenvalue=-0.626\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "LPp41FZTlV4q"
},
"source": [
"We can see that the first and third fixed points are stable (negative eigenvalues) and the second is unstable (positive eigenvalue) - as we expected!"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "sGMgVibJlV4q"
},
"source": [
"## Bonus Section 2: Noisy input drives the transition between two stable states"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "ZFhw3ov7lV4q"
},
"source": [
"As discussed in several previous tutorials, the Ornstein-Uhlenbeck (OU) process is usually used to generate a noisy input into the neuron. The OU input $\\eta(t)$ follows:\n",
"\n",
"\\begin{equation}\n",
"\\tau_\\eta \\frac{d}{dt}\\eta(t) = -\\eta (t) + \\sigma_\\eta\\sqrt{2\\tau_\\eta}\\xi(t)\n",
"\\end{equation}\n",
"\n",
"Execute the following function `my_OU(pars, sig, myseed=False)` to generate an OU process."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {},
"id": "HdXTrgB2lV4q"
},
"outputs": [],
"source": [
"# @markdown Execute to get helper function `my_OU` and visualize OU process\n",
"\n",
"\n",
"def my_OU(pars, sig, myseed=False):\n",
" \"\"\"\n",
" A functions that generates Ornstein-Uhlenback process\n",
"\n",
" Args:\n",
" pars : parameter dictionary\n",
" sig : noise amplitute\n",
" myseed : random seed. int or boolean\n",
"\n",
" Returns:\n",
" I : Ornstein-Uhlenbeck input current\n",
" \"\"\"\n",
"\n",
" # Retrieve simulation parameters\n",
" dt, range_t = pars['dt'], pars['range_t']\n",
" Lt = range_t.size\n",
" tau_ou = pars['tau_ou'] # [ms]\n",
"\n",
" # set random seed\n",
" if myseed:\n",
" np.random.seed(seed=myseed)\n",
" else:\n",
" np.random.seed()\n",
"\n",
" # Initialize\n",
" noise = np.random.randn(Lt)\n",
" I_ou = np.zeros(Lt)\n",
" I_ou[0] = noise[0] * sig\n",
"\n",
" # generate OU\n",
" for it in range(Lt - 1):\n",
" I_ou[it + 1] = (I_ou[it]\n",
" + dt / tau_ou * (0. - I_ou[it])\n",
" + np.sqrt(2 * dt / tau_ou) * sig * noise[it + 1])\n",
"\n",
" return I_ou\n",
"\n",
"\n",
"pars = default_pars_single(T=100)\n",
"pars['tau_ou'] = 1. # [ms]\n",
"sig_ou = 0.1\n",
"I_ou = my_OU(pars, sig=sig_ou, myseed=2020)\n",
"plt.figure(figsize=(10, 4))\n",
"plt.plot(pars['range_t'], I_ou, 'r')\n",
"plt.xlabel('t (ms)')\n",
"plt.ylabel(r'$I_{\\mathrm{OU}}$')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "vfIbRr6ZlV4q"
},
"source": [
"In the presence of two or more fixed points, noisy inputs can drive a transition between the fixed points! Here, we stimulate an E population for 1,000 ms applying OU inputs."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {},
"id": "uinhtExOlV4q"
},
"outputs": [],
"source": [
"# @markdown Execute this cell to simulate E population with OU inputs\n",
"\n",
"pars = default_pars_single(T=1000)\n",
"pars['w'] = 5.0\n",
"sig_ou = 0.7\n",
"pars['tau_ou'] = 1. # [ms]\n",
"pars['I_ext'] = 0.56 + my_OU(pars, sig=sig_ou, myseed=2020)\n",
"\n",
"r = simulate_single(pars)\n",
"\n",
"plt.figure(figsize=(10, 4))\n",
"plt.plot(pars['range_t'], r, 'b', alpha=0.8)\n",
"plt.xlabel('t (ms)')\n",
"plt.ylabel(r'$r(t)$')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {},
"id": "1j4VDl3ZlV4q"
},
"source": [
"You can see that the population activity is changing between fixed points (it goes up and down)!"
]
}
],
"metadata": {
"colab": {
"provenance": [],
"toc_visible": true
},
"kernel": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"kernelspec": {
"display_name": "Python 3",
"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.9.17"
}
},
"nbformat": 4,
"nbformat_minor": 0
}