{ "cells": [ { "cell_type": "markdown", "metadata": { "colab_type": "text", "execution": {}, "id": "view-in-github" }, "source": [ "\"Open   \"Open" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "# Tutorial 2: Wilson-Cowan Model\n", "\n", "**Week 2, Day 4: Dynamic Networks**\n", "\n", "**By Neuromatch Academy**\n", "\n", "**Content creators:** Qinglong Gu, Songtin Li, Arvind Kumar, John Murray, Julijana Gjorgjieva\n", "\n", "**Content reviewers:** Maryam Vaziri-Pashkam, Ella Batty, Lorenzo Fontolan, Richard Gao, Spiros Chavlis, Michael Waskom, Siddharth Suresh\n", "\n", "**Production editors:** Gagana B, Spiros Chavlis" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "---\n", "# Tutorial Objectives\n", "\n", "*Estimated timing of tutorial: 1 hour, 35 minutes*\n", "\n", "In the previous tutorial, you became familiar with a neuronal network consisting of only an excitatory population. Here, we extend the approach we used to include both excitatory and inhibitory neuronal populations in our network. A simple, yet powerful model to study the dynamics of two interacting populations of excitatory and inhibitory neurons, is the so-called **Wilson-Cowan** rate model, which will be the subject of this tutorial.\n", "\n", "The objectives of this tutorial are to:\n", "\n", "- Write the **Wilson-Cowan** equations for the firing rate dynamics of a 2D system composed of an excitatory (E) and an inhibitory (I) population of neurons\n", "- Simulate the dynamics of the system, i.e., Wilson-Cowan model.\n", "- Plot the frequency-current (F-I) curves for both populations (i.e., E and I).\n", "- Visualize and inspect the behavior of the system using **phase plane analysis**, **vector fields**, and **nullclines**.\n", "\n", "Bonus steps:\n", "\n", "- Find and plot the **fixed points** of the Wilson-Cowan model.\n", "- Investigate the stability of the Wilson-Cowan model by linearizing its dynamics and examining the **Jacobian matrix**.\n", "- Learn how the Wilson-Cowan model can reach an oscillatory state.\n", "\n", "Bonus steps (applications):\n", "- Visualize the behavior of an Inhibition-stabilized network.\n", "- Simulate working memory using the Wilson-Cowan model.\n", "\n", "
\n", "\n", "Reference paper:\n", "\n", "Wilson, H., and Cowan, J. (1972). Excitatory and inhibitory interactions in localized populations of model neurons. _Biophysical Journal_ **12**. doi: [10.1016/S0006-3495(72)86068-5](https://doi.org/10.1016/S0006-3495(72)86068-5)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "remove-input" ] }, "outputs": [], "source": [ "# @markdown\n", "from IPython.display import IFrame\n", "from ipywidgets import widgets\n", "out = widgets.Output()\n", "with out:\n", " print(f\"If you want to download the slides: https://osf.io/download/nvuty/\")\n", " display(IFrame(src=f\"https://mfr.ca-1.osf.io/render?url=https://osf.io/nvuty/?direct%26mode=render%26action=download%26mode=render\", width=730, height=410))\n", "display(out)" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "---\n", "# Setup" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Install and import feedback gadget\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# @title Install and import feedback gadget\n", "\n", "!pip3 install vibecheck datatops --quiet\n", "\n", "from vibecheck import DatatopsContentReviewContainer\n", "def content_review(notebook_section: str):\n", " return DatatopsContentReviewContainer(\n", " \"\", # No text prompt\n", " notebook_section,\n", " {\n", " \"url\": \"https://pmyvdlilci.execute-api.us-east-1.amazonaws.com/klab\",\n", " \"name\": \"neuromatch_cn\",\n", " \"user_key\": \"y1x3mpx5\",\n", " },\n", " ).render()\n", "\n", "\n", "feedback_prefix = \"W2D4_T2\"" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "both", "execution": {} }, "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": "markdown", "metadata": {}, "source": [ "## Figure Settings\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "hide-input" ] }, "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": "markdown", "metadata": {}, "source": [ "## Plotting Functions\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# @title Plotting Functions\n", "\n", "def plot_FI_inverse(x, a, theta):\n", " f, ax = plt.subplots()\n", " ax.plot(x, F_inv(x, a=a, theta=theta))\n", " ax.set(xlabel=\"$x$\", ylabel=\"$F^{-1}(x)$\")\n", "\n", "\n", "def plot_FI_EI(x, FI_exc, FI_inh):\n", " plt.figure()\n", " plt.plot(x, FI_exc, 'b', label='E population')\n", " plt.plot(x, FI_inh, 'r', label='I population')\n", " plt.legend(loc='lower right')\n", " plt.xlabel('x (a.u.)')\n", " plt.ylabel('F(x)')\n", " plt.show()\n", "\n", "\n", "def my_test_plot(t, rE1, rI1, rE2, rI2):\n", "\n", " plt.figure()\n", " ax1 = plt.subplot(211)\n", " ax1.plot(pars['range_t'], rE1, 'b', label='E population')\n", " ax1.plot(pars['range_t'], rI1, 'r', label='I population')\n", " ax1.set_ylabel('Activity')\n", " ax1.legend(loc='best')\n", "\n", " ax2 = plt.subplot(212, sharex=ax1, sharey=ax1)\n", " ax2.plot(pars['range_t'], rE2, 'b', label='E population')\n", " ax2.plot(pars['range_t'], rI2, 'r', label='I population')\n", " ax2.set_xlabel('t (ms)')\n", " ax2.set_ylabel('Activity')\n", " ax2.legend(loc='best')\n", "\n", " plt.tight_layout()\n", " plt.show()\n", "\n", "\n", "def plot_nullclines(Exc_null_rE, Exc_null_rI, Inh_null_rE, Inh_null_rI):\n", "\n", " plt.figure()\n", " plt.plot(Exc_null_rE, Exc_null_rI, 'b', label='E nullcline')\n", " plt.plot(Inh_null_rE, Inh_null_rI, 'r', label='I nullcline')\n", " plt.xlabel(r'$r_E$')\n", " plt.ylabel(r'$r_I$')\n", " plt.legend(loc='best')\n", " plt.show()\n", "\n", "\n", "def my_plot_nullcline(pars):\n", " Exc_null_rE = np.linspace(-0.01, 0.96, 100)\n", " Exc_null_rI = get_E_nullcline(Exc_null_rE, **pars)\n", " Inh_null_rI = np.linspace(-.01, 0.8, 100)\n", " Inh_null_rE = get_I_nullcline(Inh_null_rI, **pars)\n", "\n", " plt.plot(Exc_null_rE, Exc_null_rI, 'b', label='E nullcline')\n", " plt.plot(Inh_null_rE, Inh_null_rI, 'r', label='I nullcline')\n", " plt.xlabel(r'$r_E$')\n", " plt.ylabel(r'$r_I$')\n", " plt.legend(loc='best')\n", "\n", "\n", "def my_plot_vector(pars, my_n_skip=2, myscale=5):\n", " EI_grid = np.linspace(0., 1., 20)\n", " rE, rI = np.meshgrid(EI_grid, EI_grid)\n", " drEdt, drIdt = EIderivs(rE, rI, **pars)\n", "\n", " n_skip = my_n_skip\n", "\n", " plt.quiver(rE[::n_skip, ::n_skip], rI[::n_skip, ::n_skip],\n", " drEdt[::n_skip, ::n_skip], drIdt[::n_skip, ::n_skip],\n", " angles='xy', scale_units='xy', scale=myscale, facecolor='c')\n", "\n", " plt.xlabel(r'$r_E$')\n", " plt.ylabel(r'$r_I$')\n", "\n", "\n", "def my_plot_trajectory(pars, mycolor, x_init, mylabel):\n", " pars = pars.copy()\n", " pars['rE_init'], pars['rI_init'] = x_init[0], x_init[1]\n", " rE_tj, rI_tj = simulate_wc(**pars)\n", "\n", " plt.plot(rE_tj, rI_tj, color=mycolor, label=mylabel)\n", " plt.plot(x_init[0], x_init[1], 'o', color=mycolor, ms=8)\n", " plt.xlabel(r'$r_E$')\n", " plt.ylabel(r'$r_I$')\n", "\n", "\n", "def my_plot_trajectories(pars, dx, n, mylabel):\n", " \"\"\"\n", " Solve for I along the E_grid from dE/dt = 0.\n", "\n", " Expects:\n", " pars : Parameter dictionary\n", " dx : increment of initial values\n", " n : n*n trjectories\n", " mylabel : label for legend\n", "\n", " Returns:\n", " figure of trajectory\n", " \"\"\"\n", " pars = pars.copy()\n", " for ie in range(n):\n", " for ii in range(n):\n", " pars['rE_init'], pars['rI_init'] = dx * ie, dx * ii\n", " rE_tj, rI_tj = simulate_wc(**pars)\n", " if (ie == n-1) & (ii == n-1):\n", " plt.plot(rE_tj, rI_tj, 'gray', alpha=0.8, label=mylabel)\n", " else:\n", " plt.plot(rE_tj, rI_tj, 'gray', alpha=0.8)\n", "\n", " plt.xlabel(r'$r_E$')\n", " plt.ylabel(r'$r_I$')\n", "\n", "\n", "def plot_complete_analysis(pars):\n", " plt.figure(figsize=(7.7, 6.))\n", "\n", " # plot example trajectories\n", " my_plot_trajectories(pars, 0.2, 6,\n", " 'Sample trajectories \\nfor different init. conditions')\n", " my_plot_trajectory(pars, 'orange', [0.6, 0.8],\n", " 'Sample trajectory for \\nlow activity')\n", " my_plot_trajectory(pars, 'm', [0.6, 0.6],\n", " 'Sample trajectory for \\nhigh activity')\n", "\n", " # plot nullclines\n", " my_plot_nullcline(pars)\n", "\n", " # plot vector field\n", " EI_grid = np.linspace(0., 1., 20)\n", " rE, rI = np.meshgrid(EI_grid, EI_grid)\n", " drEdt, drIdt = EIderivs(rE, rI, **pars)\n", " n_skip = 2\n", " plt.quiver(rE[::n_skip, ::n_skip], rI[::n_skip, ::n_skip],\n", " drEdt[::n_skip, ::n_skip], drIdt[::n_skip, ::n_skip],\n", " angles='xy', scale_units='xy', scale=5., facecolor='c')\n", "\n", " plt.legend(loc=[1.02, 0.57], handlelength=1)\n", " plt.show()\n", "\n", "\n", "def plot_fp(x_fp, position=(0.02, 0.1), rotation=0):\n", " plt.plot(x_fp[0], x_fp[1], 'ko', ms=8)\n", " plt.text(x_fp[0] + position[0], x_fp[1] + position[1],\n", " f'Fixed Point1=\\n({x_fp[0]:.3f}, {x_fp[1]:.3f})',\n", " horizontalalignment='center', verticalalignment='bottom',\n", " rotation=rotation)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Helper Functions\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# @title Helper Functions\n", "\n", "def default_pars(**kwargs):\n", " pars = {}\n", "\n", " # Excitatory parameters\n", " pars['tau_E'] = 1. # Timescale of the E population [ms]\n", " pars['a_E'] = 1.2 # Gain of the E population\n", " pars['theta_E'] = 2.8 # Threshold of the E population\n", "\n", " # Inhibitory parameters\n", " pars['tau_I'] = 2.0 # Timescale of the I population [ms]\n", " pars['a_I'] = 1.0 # Gain of the I population\n", " pars['theta_I'] = 4.0 # Threshold of the I population\n", "\n", " # Connection strength\n", " pars['wEE'] = 9. # E to E\n", " pars['wEI'] = 4. # I to E\n", " pars['wIE'] = 13. # E to I\n", " pars['wII'] = 11. # I to I\n", "\n", " # External input\n", " pars['I_ext_E'] = 0.\n", " pars['I_ext_I'] = 0.\n", "\n", " # simulation parameters\n", " pars['T'] = 50. # Total duration of simulation [ms]\n", " pars['dt'] = .1 # Simulation time step [ms]\n", " pars['rE_init'] = 0.2 # Initial value of E\n", " pars['rI_init'] = 0.2 # Initial value of I\n", "\n", " # External parameters if any\n", " for k in kwargs:\n", " pars[k] = kwargs[k]\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", "\n", "def F(x, a, theta):\n", " \"\"\"\n", " Population activation function, F-I curve\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", " f : the population activation response f(x) for input x\n", " \"\"\"\n", "\n", " # add the expression of f = F(x)\n", " f = (1 + np.exp(-a * (x - theta)))**-1 - (1 + np.exp(a * theta))**-1\n", "\n", " return f\n", "\n", "\n", "def dF(x, a, theta):\n", " \"\"\"\n", " Derivative of the 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 : Derivative of the population activation function.\n", " \"\"\"\n", "\n", " dFdx = a * np.exp(-a * (x - theta)) * (1 + np.exp(-a * (x - theta)))**-2\n", "\n", " return dFdx" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "The helper functions included:\n", "\n", "- Parameter dictionary: `default_pars(**kwargs)`. You can use:\n", " - `pars = default_pars()` to get all the parameters, and then you can execute `print(pars)` to check these parameters.\n", " - `pars = default_pars(T=T_sim, dt=time_step)` to set a different simulation time and time step\n", " - After `pars = default_pars()`, use `par['New_para'] = value` to add a new parameter with its value\n", " - Pass to functions that accept individual parameters with `func(**pars)`\n", "- F-I curve: `F(x, a, theta)`\n", "- Derivative of the F-I curve: `dF(x, a, theta)`" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "---\n", "# Section 1: Wilson-Cowan model of excitatory and inhibitory populations\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Video 1: Phase analysis of the Wilson-Cowan E-I model\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "remove-input" ] }, "outputs": [], "source": [ "# @title Video 1: Phase analysis of the Wilson-Cowan E-I model\n", "from ipywidgets import widgets\n", "from IPython.display import YouTubeVideo\n", "from IPython.display import IFrame\n", "from IPython.display import display\n", "\n", "\n", "class PlayVideo(IFrame):\n", " def __init__(self, id, source, page=1, width=400, height=300, **kwargs):\n", " self.id = id\n", " if source == 'Bilibili':\n", " src = f'https://player.bilibili.com/player.html?bvid={id}&page={page}'\n", " elif source == 'Osf':\n", " src = f'https://mfr.ca-1.osf.io/render?url=https://osf.io/download/{id}/?direct%26mode=render'\n", " super(PlayVideo, self).__init__(src, width, height, **kwargs)\n", "\n", "\n", "def display_videos(video_ids, W=400, H=300, fs=1):\n", " tab_contents = []\n", " for i, video_id in enumerate(video_ids):\n", " out = widgets.Output()\n", " with out:\n", " if video_ids[i][0] == 'Youtube':\n", " video = YouTubeVideo(id=video_ids[i][1], width=W,\n", " height=H, fs=fs, rel=0)\n", " print(f'Video available at https://youtube.com/watch?v={video.id}')\n", " else:\n", " video = PlayVideo(id=video_ids[i][1], source=video_ids[i][0], width=W,\n", " height=H, fs=fs, autoplay=False)\n", " if video_ids[i][0] == 'Bilibili':\n", " print(f'Video available at https://www.bilibili.com/video/{video.id}')\n", " elif video_ids[i][0] == 'Osf':\n", " print(f'Video available at https://osf.io/{video.id}')\n", " display(video)\n", " tab_contents.append(out)\n", " return tab_contents\n", "\n", "\n", "video_ids = [('Youtube', 'GCpQmh45crM'), ('Bilibili', 'BV1CD4y1m7dK')]\n", "tab_contents = display_videos(video_ids, W=730, H=410)\n", "tabs = widgets.Tab()\n", "tabs.children = tab_contents\n", "for i in range(len(tab_contents)):\n", " tabs.set_title(i, video_ids[i][0])\n", "display(tabs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Submit your feedback\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# @title Submit your feedback\n", "content_review(f\"{feedback_prefix}_Phase_analysis_Video\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "This video explains how to model a network with interacting populations of excitatory and inhibitory neurons (the Wilson-Cowan model). It shows how to solve the network activity vs. time and introduces the phase plane in two dimensions." ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "## Section 1.1: Mathematical description of the WC model\n", "\n", "*Estimated timing to here from start of tutorial: 12 min*\n", "\n", "
\n", " Click here for text recap of relevant part of video \n", "\n", "Many of the rich dynamics recorded in the brain are generated by the interaction of excitatory and inhibitory subtype neurons. Here, similar to what we did in the previous tutorial, we will model two coupled populations of E and I neurons (**Wilson-Cowan** model). We can write two coupled differential equations, each representing the dynamics of the excitatory or inhibitory population:\n", "\n", "\\begin{align}\n", "\\tau_E \\frac{dr_E}{dt} &= -r_E + F_E(w_{EE}r_E -w_{EI}r_I + I^{\\text{ext}}_E;a_E,\\theta_E)\\\\\n", "\\tau_I \\frac{dr_I}{dt} &= -r_I + F_I(w_{IE}r_E -w_{II}r_I + I^{\\text{ext}}_I;a_I,\\theta_I) \\qquad (1)\n", "\\end{align}\n", "\n", "$r_E(t)$ represents the average activation (or firing rate) of the excitatory population at time $t$, and $r_I(t)$ the activation (or firing rate) of the inhibitory population. The parameters $\\tau_E$ and $\\tau_I$ control the timescales of the dynamics of each population. Connection strengths are given by: $w_{EE}$ (E $\\rightarrow$ E), $w_{EI}$ (I $\\rightarrow$ E), $w_{IE}$ (E $\\rightarrow$ I), and $w_{II}$ (I $\\rightarrow$ I). The terms $w_{EI}$ and $w_{IE}$ represent connections from inhibitory to excitatory population and vice versa, respectively. The transfer functions (or F-I curves) $F_E(x;a_E,\\theta_E)$ and $F_I(x;a_I,\\theta_I)$ can be different for the excitatory and the inhibitory populations.\n" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "### Coding Exercise 1.1: Plot out the F-I curves for the E and I populations\n", "\n", "Let's first plot out the F-I curves for the E and I populations using the helper function `F` with default parameter values." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {} }, "outputs": [], "source": [ "help(F)" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "execution": {} }, "source": [ "```python\n", "pars = default_pars()\n", "x = np.arange(0, 10, .1)\n", "\n", "print(pars['a_E'], pars['theta_E'])\n", "print(pars['a_I'], pars['theta_I'])\n", "\n", "###################################################################\n", "# TODO for students: compute and plot the F-I curve here #\n", "# Note: aE, thetaE, aI and theta_I are in the dictionary 'pairs' #\n", "raise NotImplementedError('student exercise: compute F-I curves of excitatory and inhibitory populations')\n", "###################################################################\n", "\n", "# Compute the F-I curve of the excitatory population\n", "FI_exc = ...\n", "\n", "# Compute the F-I curve of the inhibitory population\n", "FI_inh = ...\n", "\n", "# Visualize\n", "plot_FI_EI(x, FI_exc, FI_inh)\n", "\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {} }, "outputs": [], "source": [ "# to_remove solution\n", "pars = default_pars()\n", "x = np.arange(0, 10, .1)\n", "\n", "print(pars['a_E'], pars['theta_E'])\n", "print(pars['a_I'], pars['theta_I'])\n", "\n", "# Compute the F-I curve of the excitatory population\n", "FI_exc = F(x, pars['a_E'], pars['theta_E'])\n", "\n", "# Compute the F-I curve of the inhibitory population\n", "FI_inh = F(x, pars['a_I'], pars['theta_I'])\n", "\n", "# Visualize\n", "with plt.xkcd():\n", " plot_FI_EI(x, FI_exc, FI_inh)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Submit your feedback\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# @title Submit your feedback\n", "content_review(f\"{feedback_prefix}_Plot_FI_Exercise\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "## Section 1.2: Simulation scheme for the Wilson-Cowan model\n", "\n", "*Estimated timing to here from start of tutorial: 20 min*\n", "\n", "Once again, we can integrate our equations numerically. Using the Euler method, the dynamics of E and I populations can be simulated on a time-grid of stepsize $\\Delta t$. The updates for the activity of the excitatory and the inhibitory populations can be written as:\n", "\n", "\\begin{align}\n", "r_E[k+1] &= r_E[k] + \\Delta r_E[k]\\\\\n", "r_I[k+1] &= r_I[k] + \\Delta r_I[k]\n", "\\end{align}\n", "\n", "with the increments\n", "\n", "\\begin{align}\n", "\\Delta r_E[k] &= \\frac{\\Delta t}{\\tau_E}[-r_E[k] + F_E(w_{EE}r_E[k] -w_{EI}r_I[k] + I^{\\text{ext}}_E[k];a_E,\\theta_E)]\\\\\n", "\\Delta r_I[k] &= \\frac{\\Delta t}{\\tau_I}[-r_I[k] + F_I(w_{IE}r_E[k] -w_{II}r_I[k] + I^{\\text{ext}}_I[k];a_I,\\theta_I)]\n", "\\end{align}" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "### Coding Exercise 1.2: Numerically integrate the Wilson-Cowan equations\n", "\n", "We will implement this numerical simulation of our equations and visualize two simulations with similar initial points." ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "execution": {} }, "source": [ "```python\n", "def simulate_wc(tau_E, a_E, theta_E, tau_I, a_I, theta_I,\n", " wEE, wEI, wIE, wII, I_ext_E, I_ext_I,\n", " rE_init, rI_init, dt, range_t, **other_pars):\n", " \"\"\"\n", " Simulate the Wilson-Cowan equations\n", "\n", " Args:\n", " Parameters of the Wilson-Cowan model\n", "\n", " Returns:\n", " rE, rI (arrays) : Activity of excitatory and inhibitory populations\n", " \"\"\"\n", " # Initialize activity arrays\n", " Lt = range_t.size\n", " rE = np.append(rE_init, np.zeros(Lt - 1))\n", " rI = np.append(rI_init, np.zeros(Lt - 1))\n", " I_ext_E = I_ext_E * np.ones(Lt)\n", " I_ext_I = I_ext_I * np.ones(Lt)\n", "\n", " # Simulate the Wilson-Cowan equations\n", " for k in range(Lt - 1):\n", "\n", " ########################################################################\n", " # TODO for students: compute drE and drI and remove the error\n", " raise NotImplementedError(\"Student exercise: compute the change in E/I\")\n", " ########################################################################\n", "\n", " # Calculate the derivative of the E population\n", " drE = ...\n", "\n", " # Calculate the derivative of the I population\n", " drI = ...\n", "\n", " # Update using Euler's method\n", " rE[k + 1] = rE[k] + drE\n", " rI[k + 1] = rI[k] + drI\n", "\n", " return rE, rI\n", "\n", "\n", "pars = default_pars()\n", "\n", "# Simulate first trajectory\n", "rE1, rI1 = simulate_wc(**default_pars(rE_init=.32, rI_init=.15))\n", "\n", "# Simulate second trajectory\n", "rE2, rI2 = simulate_wc(**default_pars(rE_init=.33, rI_init=.15))\n", "\n", "# Visualize\n", "my_test_plot(pars['range_t'], rE1, rI1, rE2, rI2)\n", "\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "both", "execution": {} }, "outputs": [], "source": [ "# to_remove solution\n", "def simulate_wc(tau_E, a_E, theta_E, tau_I, a_I, theta_I,\n", " wEE, wEI, wIE, wII, I_ext_E, I_ext_I,\n", " rE_init, rI_init, dt, range_t, **other_pars):\n", " \"\"\"\n", " Simulate the Wilson-Cowan equations\n", "\n", " Args:\n", " Parameters of the Wilson-Cowan model\n", "\n", " Returns:\n", " rE, rI (arrays) : Activity of excitatory and inhibitory populations\n", " \"\"\"\n", " # Initialize activity arrays\n", " Lt = range_t.size\n", " rE = np.append(rE_init, np.zeros(Lt - 1))\n", " rI = np.append(rI_init, np.zeros(Lt - 1))\n", " I_ext_E = I_ext_E * np.ones(Lt)\n", " I_ext_I = I_ext_I * np.ones(Lt)\n", "\n", " # Simulate the Wilson-Cowan equations\n", " for k in range(Lt - 1):\n", "\n", " # Calculate the derivative of the E population\n", " drE = dt / tau_E * (-rE[k] + F(wEE * rE[k] - wEI * rI[k] + I_ext_E[k],\n", " a_E, theta_E))\n", "\n", " # Calculate the derivative of the I population\n", " drI = dt / tau_I * (-rI[k] + F(wIE * rE[k] - wII * rI[k] + I_ext_I[k],\n", " a_I, theta_I))\n", "\n", " # Update using Euler's method\n", " rE[k + 1] = rE[k] + drE\n", " rI[k + 1] = rI[k] + drI\n", "\n", " return rE, rI\n", "\n", "\n", "pars = default_pars()\n", "\n", "# Simulate first trajectory\n", "rE1, rI1 = simulate_wc(**default_pars(rE_init=.32, rI_init=.15))\n", "\n", "# Simulate second trajectory\n", "rE2, rI2 = simulate_wc(**default_pars(rE_init=.33, rI_init=.15))\n", "\n", "# Visualize\n", "with plt.xkcd():\n", " my_test_plot(pars['range_t'], rE1, rI1, rE2, rI2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Submit your feedback\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# @title Submit your feedback\n", "content_review(f\"{feedback_prefix}_Numerical_Integration_of_WE_model_Exercise\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "The two plots above show the temporal evolution of excitatory ($r_E$, blue) and inhibitory ($r_I$, red) activity for two different sets of initial conditions." ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "### Interactive Demo 1.2: population trajectories with different initial values\n", "\n", "In this interactive demo, we will simulate the Wilson-Cowan model and plot the trajectories of each population. We change the initial activity of the excitatory population.\n", "\n", "What happens to the E and I population trajectories with different initial conditions?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Make sure you execute this cell to enable the widget!\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# @markdown Make sure you execute this cell to enable the widget!\n", "@widgets.interact(\n", " rE_init=widgets.FloatSlider(0.32, min=0.30, max=0.35, step=.01)\n", ")\n", "\n", "\n", "def plot_EI_diffinitial(rE_init=0.0):\n", "\n", " pars = default_pars(rE_init=rE_init, rI_init=.15)\n", " rE, rI = simulate_wc(**pars)\n", "\n", " plt.figure()\n", " plt.plot(pars['range_t'], rE, 'b', label='E population')\n", " plt.plot(pars['range_t'], rI, 'r', label='I population')\n", " plt.xlabel('t (ms)')\n", " plt.ylabel('Activity')\n", " plt.legend(loc='best')\n", " plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {} }, "outputs": [], "source": [ "# to_remove explanation\n", "\"\"\"\n", "Sometimes both the excitatory and inhibitory population activities decay to 0 and\n", "sometimes they level out at different positive values.\n", "\"\"\";" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Submit your feedback\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# @title Submit your feedback\n", "content_review(f\"{feedback_prefix}_Population_trajectories_with_different_initial_values_Interactive_Demo_and_Discussion\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "### Think! 1.2\n", "It is evident that the steady states of the neuronal response can be different when different initial states are chosen. Why is that?\n", "\n", "We will discuss this in the next section but try to think about it first." ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "---\n", "# Section 2: Phase plane analysis\n", "\n", "*Estimated timing to here from start of tutorial: 45 min*\n", "\n", "Just like we used a graphical method to study the dynamics of a 1-D system in the previous tutorial, here we will learn a graphical approach called **phase plane analysis** to study the dynamics of a 2-D system like the Wilson-Cowan model. You have seen this before in the [pre-reqs calculus day](https://compneuro.neuromatch.io/tutorials/W0D4_Calculus/student/W0D4_Tutorial3.html#section-3-2-phase-plane-plot-and-nullcline) and on the [Linear Systems day](https://compneuro.neuromatch.io/tutorials/W2D2_LinearSystems/student/W2D2_Tutorial1.html#section-4-stream-plots)\n", "\n", "So far, we have plotted the activities of the two populations as a function of time, i.e., in the `Activity-t` plane, either the $(t, r_E(t))$ plane or the $(t, r_I(t))$ one. Instead, we can plot the two activities $r_E(t)$ and $r_I(t)$ against each other at any time point $t$. This characterization in the `rI-rE` plane $(r_I(t), r_E(t))$ is called the **phase plane**. Each line in the phase plane indicates how both $r_E$ and $r_I$ evolve with time." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Video 2: Nullclines and Vector Fields\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "remove-input" ] }, "outputs": [], "source": [ "# @title Video 2: Nullclines and Vector Fields\n", "from ipywidgets import widgets\n", "from IPython.display import YouTubeVideo\n", "from IPython.display import IFrame\n", "from IPython.display import display\n", "\n", "\n", "class PlayVideo(IFrame):\n", " def __init__(self, id, source, page=1, width=400, height=300, **kwargs):\n", " self.id = id\n", " if source == 'Bilibili':\n", " src = f'https://player.bilibili.com/player.html?bvid={id}&page={page}'\n", " elif source == 'Osf':\n", " src = f'https://mfr.ca-1.osf.io/render?url=https://osf.io/download/{id}/?direct%26mode=render'\n", " super(PlayVideo, self).__init__(src, width, height, **kwargs)\n", "\n", "\n", "def display_videos(video_ids, W=400, H=300, fs=1):\n", " tab_contents = []\n", " for i, video_id in enumerate(video_ids):\n", " out = widgets.Output()\n", " with out:\n", " if video_ids[i][0] == 'Youtube':\n", " video = YouTubeVideo(id=video_ids[i][1], width=W,\n", " height=H, fs=fs, rel=0)\n", " print(f'Video available at https://youtube.com/watch?v={video.id}')\n", " else:\n", " video = PlayVideo(id=video_ids[i][1], source=video_ids[i][0], width=W,\n", " height=H, fs=fs, autoplay=False)\n", " if video_ids[i][0] == 'Bilibili':\n", " print(f'Video available at https://www.bilibili.com/video/{video.id}')\n", " elif video_ids[i][0] == 'Osf':\n", " print(f'Video available at https://osf.io/{video.id}')\n", " display(video)\n", " tab_contents.append(out)\n", " return tab_contents\n", "\n", "\n", "video_ids = [('Youtube', 'V2SBAK2Xf8Y'), ('Bilibili', 'BV15k4y1m7Kt')]\n", "tab_contents = display_videos(video_ids, W=730, H=410)\n", "tabs = widgets.Tab()\n", "tabs.children = tab_contents\n", "for i in range(len(tab_contents)):\n", " tabs.set_title(i, video_ids[i][0])\n", "display(tabs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Submit your feedback\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# @title Submit your feedback\n", "content_review(f\"{feedback_prefix}_Nullclines_and_Vector_Fields_Video\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "## Interactive Demo 2: From the Activity - time plane to the **$r_I$ - $r_E$** phase plane\n", "\n", "In this demo, we will visualize the system dynamics using both the `Activity-time` and the `(rE, rI)` phase plane. The circles indicate the activities at a given time $t$, while the lines represent the evolution of the system for the entire duration of the simulation.\n", "\n", "Move the time slider to better understand how the top plot relates to the bottom plot. Does the bottom plot have explicit information about time? What information does it give us?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Make sure you execute this cell to enable the widget!\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# @markdown Make sure you execute this cell to enable the widget!\n", "\n", "pars = default_pars(T=10, rE_init=0.6, rI_init=0.8)\n", "rE, rI = simulate_wc(**pars)\n", "\n", "@widgets.interact(\n", " n_t=widgets.IntSlider(0, min=0, max=len(pars['range_t']) - 1, step=1)\n", ")\n", "\n", "\n", "def plot_activity_phase(n_t):\n", " plt.figure(figsize=(8, 5.5))\n", " plt.subplot(211)\n", " plt.plot(pars['range_t'], rE, 'b', label=r'$r_E$')\n", " plt.plot(pars['range_t'], rI, 'r', label=r'$r_I$')\n", " plt.plot(pars['range_t'][n_t], rE[n_t], 'bo')\n", " plt.plot(pars['range_t'][n_t], rI[n_t], 'ro')\n", " plt.axvline(pars['range_t'][n_t], 0, 1, color='k', ls='--')\n", " plt.xlabel('t (ms)', fontsize=14)\n", " plt.ylabel('Activity', fontsize=14)\n", " plt.legend(loc='best', fontsize=14)\n", "\n", " plt.subplot(212)\n", " plt.plot(rE, rI, 'k')\n", " plt.plot(rE[n_t], rI[n_t], 'ko')\n", " plt.xlabel(r'$r_E$', fontsize=18, color='b')\n", " plt.ylabel(r'$r_I$', fontsize=18, color='r')\n", "\n", " plt.tight_layout()\n", " plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {} }, "outputs": [], "source": [ "# to_remove explanation\n", "\n", "\"\"\"\n", "- Phase plane portraits allows us to visualize out of all possible states which\n", "states a system can take. For example in this demo you see that the among all\n", "possible pairs of values for r1 and r2, this system can take only a limited number\n", "of states (those that lie on the black line).\n", "\n", "- There are other things we can infer from the phase portraits e.g. fixed\n", "points, trajectory to the fixed points etc.\n", "\n", "- Explicit information about time is not visible in the phase portraits\n", "\"\"\";" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Submit your feedback\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# @title Submit your feedback\n", "content_review(f\"{feedback_prefix}_Time_plane_to_Phase_plane_Interactive_Demo_and_Discussion\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "## Section 2.1: Nullclines of the Wilson-Cowan Equations\n", "\n", "*Estimated timing to here from start of tutorial: 1 hour, 3 min*\n", "\n", "An important concept in the phase plane analysis is the \"nullcline\" which is defined as the set of points in the phase plane where the activity of one population (but not necessarily the other) does not change.\n", "\n", "In other words, the $E$ and $I$ nullclines of Equation $(1)$ are defined as the points where $\\displaystyle{\\frac{dr_E}{dt}}=0$, for the excitatory nullcline, or $\\displaystyle\\frac{dr_I}{dt}=0$ for the inhibitory nullcline. That is:\n", "\n", "\\begin{align}\n", "-r_E + F_E(w_{EE}r_E -w_{EI}r_I + I^{\\text{ext}}_E;a_E,\\theta_E) &= 0 \\qquad (2)\\\\[1mm]\n", "-r_I + F_I(w_{IE}r_E -w_{II}r_I + I^{\\text{ext}}_I;a_I,\\theta_I) &= 0 \\qquad (3)\n", "\\end{align}" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "### Coding Exercise 2.1: Compute the nullclines of the Wilson-Cowan model\n", "\n", "In the next exercise, we will compute and plot the nullclines of the E and I population." ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Along the nullcline of excitatory population Equation $2$, you can calculate the inhibitory activity by rewriting Equation $2$ into\n", "\n", "\\begin{align}\n", "r_I = \\frac{1}{w_{EI}}\\big{[}w_{EE}r_E - F_E^{-1}(r_E; a_E,\\theta_E) + I^{\\text{ext}}_E \\big{]}. \\qquad(4)\n", "\\end{align}\n", "\n", "where $F_E^{-1}(r_E; a_E,\\theta_E)$ is the inverse of the excitatory transfer function (defined below). Equation $4$ defines the $r_E$ nullcline." ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Along the nullcline of inhibitory population Equation $3$, you can calculate the excitatory activity by rewriting Equation $3$ into\n", "\\begin{align}\n", "r_E = \\frac{1}{w_{IE}} \\big{[} w_{II}r_I + F_I^{-1}(r_I;a_I,\\theta_I) - I^{\\text{ext}}_I \\big{]} \\qquad (5)\n", "\\end{align}\n", "\n", "shere $F_I^{-1}(r_I; a_I,\\theta_I)$ is the inverse of the inhibitory transfer function (defined below). Equation $5$ defines the $I$ nullcline." ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Note that, when computing the nullclines with Equations 4-5, we also need to calculate the inverse of the transfer functions.\n", "\n", "The inverse of the sigmoid shaped **f-I** function that we have been using is:\n", "\n", "\\begin{equation}\n", "F^{-1}(x; a, \\theta) = -\\frac{1}{a} \\ln \\left[ \\frac{1}{x + \\displaystyle \\frac{1}{1+\\text{e}^{a\\theta}}} - 1 \\right] + \\theta \\qquad (6)\n", "\\end{equation}\n", "\n", "The first step is to implement the inverse transfer function:" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "execution": {} }, "source": [ "```python\n", "def F_inv(x, a, theta):\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", " F_inverse : value of the inverse function\n", " \"\"\"\n", "\n", " #########################################################################\n", " # TODO for students: compute F_inverse\n", " raise NotImplementedError(\"Student exercise: compute the inverse of F(x)\")\n", " #########################################################################\n", "\n", " # Calculate Finverse (ln(x) can be calculated as np.log(x))\n", " F_inverse = ...\n", "\n", " return F_inverse\n", "\n", "\n", "# Set parameters\n", "pars = default_pars()\n", "x = np.linspace(1e-6, 1, 100)\n", "\n", "# Get inverse and visualize\n", "plot_FI_inverse(x, a=1, theta=3)\n", "\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {} }, "outputs": [], "source": [ "# to_remove solution\n", "def F_inv(x, a, theta):\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", " F_inverse : value of the inverse function\n", " \"\"\"\n", "\n", " # Calculate Finverse (ln(x) can be calculated as np.log(x))\n", " F_inverse = -1/a * np.log((x + (1 + np.exp(a * theta))**-1)**-1 - 1) + theta\n", "\n", " return F_inverse\n", "\n", "\n", "# Set parameters\n", "pars = default_pars()\n", "x = np.linspace(1e-6, 1, 100)\n", "\n", "# Get inverse and visualize\n", "with plt.xkcd():\n", " plot_FI_inverse(x, a=1, theta=3)" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Now you can compute the nullclines, using Equations 4-5 (repeated here for ease of access):\n", "\n", "\\begin{align}\n", "r_I &= \\frac{1}{w_{EI}}\\big{[}w_{EE}r_E - F_E^{-1}(r_E; a_E,\\theta_E) + I^{\\text{ext}}_E \\big{]} &\\qquad(4) \\\\\n", "r_E &= \\frac{1}{w_{IE}} \\big{[} w_{II}r_I + F_I^{-1}(r_I;a_I,\\theta_I) - I^{\\text{ext}}_I \\big{]} &\\qquad (5)\n", "\\end{align}" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "execution": {} }, "source": [ "```python\n", "def get_E_nullcline(rE, a_E, theta_E, wEE, wEI, I_ext_E, **other_pars):\n", " \"\"\"\n", " Solve for rI along the rE from drE/dt = 0.\n", "\n", " Args:\n", " rE : response of excitatory population\n", " a_E, theta_E, wEE, wEI, I_ext_E : Wilson-Cowan excitatory parameters\n", " Other parameters are ignored\n", "\n", " Returns:\n", " rI : values of inhibitory population along the nullcline on the rE\n", " \"\"\"\n", " #########################################################################\n", " # TODO for students: compute rI for rE nullcline and disable the error\n", " raise NotImplementedError(\"Student exercise: compute the E nullcline\")\n", " #########################################################################\n", "\n", " # calculate rI for E nullclines on rI\n", " rI = ...\n", "\n", " return rI\n", "\n", "\n", "def get_I_nullcline(rI, a_I, theta_I, wIE, wII, I_ext_I, **other_pars):\n", " \"\"\"\n", " Solve for E along the rI from dI/dt = 0.\n", "\n", " Args:\n", " rI : response of inhibitory population\n", " a_I, theta_I, wIE, wII, I_ext_I : Wilson-Cowan inhibitory parameters\n", " Other parameters are ignored\n", "\n", " Returns:\n", " rE : values of the excitatory population along the nullcline on the rI\n", " \"\"\"\n", " #########################################################################\n", " # TODO for students: compute rI for rE nullcline and disable the error\n", " raise NotImplementedError(\"Student exercise: compute the I nullcline\")\n", " #########################################################################\n", "\n", " # calculate rE for I nullclines on rI\n", " rE = ...\n", "\n", " return rE\n", "\n", "\n", "# Set parameters\n", "pars = default_pars()\n", "Exc_null_rE = np.linspace(-0.01, 0.96, 100)\n", "Inh_null_rI = np.linspace(-.01, 0.8, 100)\n", "\n", "# Compute nullclines\n", "Exc_null_rI = get_E_nullcline(Exc_null_rE, **pars)\n", "Inh_null_rE = get_I_nullcline(Inh_null_rI, **pars)\n", "\n", "# Visualize\n", "plot_nullclines(Exc_null_rE, Exc_null_rI, Inh_null_rE, Inh_null_rI)\n", "\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "both", "execution": {} }, "outputs": [], "source": [ "# to_remove solution\n", "def get_E_nullcline(rE, a_E, theta_E, wEE, wEI, I_ext_E, **other_pars):\n", " \"\"\"\n", " Solve for rI along the rE from drE/dt = 0.\n", "\n", " Args:\n", " rE : response of excitatory population\n", " a_E, theta_E, wEE, wEI, I_ext_E : Wilson-Cowan excitatory parameters\n", " Other parameters are ignored\n", "\n", " Returns:\n", " rI : values of inhibitory population along the nullcline on the rE\n", " \"\"\"\n", " # calculate rI for E nullclines on rI\n", " rI = 1 / wEI * (wEE * rE - F_inv(rE, a_E, theta_E) + I_ext_E)\n", "\n", " return rI\n", "\n", "\n", "def get_I_nullcline(rI, a_I, theta_I, wIE, wII, I_ext_I, **other_pars):\n", " \"\"\"\n", " Solve for E along the rI from dI/dt = 0.\n", "\n", " Args:\n", " rI : response of inhibitory population\n", " a_I, theta_I, wIE, wII, I_ext_I : Wilson-Cowan inhibitory parameters\n", " Other parameters are ignored\n", "\n", " Returns:\n", " rE : values of the excitatory population along the nullcline on the rI\n", " \"\"\"\n", " # calculate rE for I nullclines on rI\n", " rE = 1 / wIE * (wII * rI + F_inv(rI, a_I, theta_I) - I_ext_I)\n", "\n", " return rE\n", "\n", "\n", "# Set parameters\n", "pars = default_pars()\n", "Exc_null_rE = np.linspace(-0.01, 0.96, 100)\n", "Inh_null_rI = np.linspace(-.01, 0.8, 100)\n", "\n", "# Compute nullclines\n", "Exc_null_rI = get_E_nullcline(Exc_null_rE, **pars)\n", "Inh_null_rE = get_I_nullcline(Inh_null_rI, **pars)\n", "\n", "# Visualize\n", "with plt.xkcd():\n", " plot_nullclines(Exc_null_rE, Exc_null_rI, Inh_null_rE, Inh_null_rI)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Submit your feedback\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# @title Submit your feedback\n", "content_review(f\"{feedback_prefix}_Compute_the_nullclines_WE_Exercise\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Note that by definition along the blue line in the phase-plane spanned by $r_E, r_I$, $\\displaystyle{\\frac{dr_E(t)}{dt}} = 0$, therefore, it is called a nullcline.\n", "\n", "That is, the blue nullcline divides the phase-plane spanned by $r_E, r_I$ into two regions: on one side of the nullcline $\\displaystyle{\\frac{dr_E(t)}{dt}} > 0$ and on the other side $\\displaystyle{\\frac{dr_E(t)}{dt}} < 0$.\n", "\n", "The same is true for the red line along which $\\displaystyle{\\frac{dr_I(t)}{dt}} = 0$. That is, the red nullcline divides the phase-plane spanned by $r_E, r_I$ into two regions: on one side of the nullcline $\\displaystyle{\\frac{dr_I(t)}{dt}} > 0$ and on the other side $\\displaystyle{\\frac{dr_I(t)}{dt}} < 0$.\n" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "## Section 2.2: Vector field\n", "\n", "*Estimated timing to here from start of tutorial: 1 hour, 20 min*\n", "\n", "How can the phase plane and the nullcline curves help us understand the behavior of the Wilson-Cowan model?\n", "\n", "
\n", " Click here for text recap of relevant part of video \n", "\n", "The activities of the $E$ and $I$ populations $r_E(t)$ and $r_I(t)$ at each time point $t$ correspond to a single point in the phase plane, with coordinates $(r_E(t),r_I(t))$. Therefore, the time-dependent trajectory of the system can be described as a continuous curve in the phase plane, and the tangent vector to the trajectory, which is defined as the vector $\\bigg{(}\\displaystyle{\\frac{dr_E(t)}{dt},\\frac{dr_I(t)}{dt}}\\bigg{)}$, indicates the direction towards which the activity is evolving and how fast is the activity changing along each axis. In fact, for each point $(E,I)$ in the phase plane, we can compute the tangent vector $\\bigg{(}\\displaystyle{\\frac{dr_E}{dt},\\frac{dr_I}{dt}}\\bigg{)}$, which indicates the behavior of the system when it traverses that point.\n", "\n", "The map of tangent vectors in the phase plane is called **vector field**. The behavior of any trajectory in the phase plane is determined by i) the initial conditions $(r_E(0),r_I(0))$, and ii) the vector field $\\bigg{(}\\displaystyle{\\frac{dr_E(t)}{dt},\\frac{dr_I(t)}{dt}}\\bigg{)}$.\n", "\n", "In general, the value of the vector field at a particular point in the phase plane is represented by an arrow. The orientation and the size of the arrow reflect the direction and the norm of the vector, respectively." ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "### Coding Exercise 2.2: Compute and plot the vector field $\\displaystyle{\\Big{(}\\frac{dr_E}{dt}, \\frac{dr_I}{dt} \\Big{)}}$\n", "\n", "Note that\n", "\n", "\\begin{align}\n", "\\frac{dr_E}{dt} &= [-r_E + F_E(w_{EE}r_E -w_{EI}r_I + I^{\\text{ext}}_E;a_E,\\theta_E)]\\frac{1}{\\tau_E}\\\\\n", "\\frac{dr_I}{dt} &= [-r_I + F_I(w_{IE}r_E -w_{II}r_I + I^{\\text{ext}}_I;a_I,\\theta_I)]\\frac{1}{\\tau_I}\n", "\\end{align}" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "execution": {} }, "source": [ "```python\n", "def EIderivs(rE, rI,\n", " tau_E, a_E, theta_E, wEE, wEI, I_ext_E,\n", " tau_I, a_I, theta_I, wIE, wII, I_ext_I,\n", " **other_pars):\n", " \"\"\"Time derivatives for E/I variables (dE/dt, dI/dt).\"\"\"\n", " ######################################################################\n", " # TODO for students: compute drEdt and drIdt and disable the error\n", " raise NotImplementedError(\"Student exercise: compute the vector field\")\n", " ######################################################################\n", "\n", " # Compute the derivative of rE\n", " drEdt = ...\n", "\n", " # Compute the derivative of rI\n", " drIdt = ...\n", "\n", " return drEdt, drIdt\n", "\n", "\n", "# Create vector field using EIderivs\n", "plot_complete_analysis(default_pars())\n", "\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "both", "execution": {} }, "outputs": [], "source": [ "# to_remove solution\n", "def EIderivs(rE, rI,\n", " tau_E, a_E, theta_E, wEE, wEI, I_ext_E,\n", " tau_I, a_I, theta_I, wIE, wII, I_ext_I,\n", " **other_pars):\n", " \"\"\"Time derivatives for E/I variables (dE/dt, dI/dt).\"\"\"\n", "\n", " # Compute the derivative of rE\n", " drEdt = (-rE + F(wEE * rE - wEI * rI + I_ext_E, a_E, theta_E)) / tau_E\n", "\n", " # Compute the derivative of rI\n", " drIdt = (-rI + F(wIE * rE - wII * rI + I_ext_I, a_I, theta_I)) / tau_I\n", "\n", " return drEdt, drIdt\n", "\n", "\n", "# Create vector field using EIderivs\n", "with plt.xkcd():\n", " plot_complete_analysis(default_pars())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Submit your feedback\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# @title Submit your feedback\n", "content_review(f\"{feedback_prefix}_Compute_the_vector_field_Exercise\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "The last phase plane plot shows us that:\n", "\n", "- Trajectories seem to follow the direction of the vector field.\n", "- Different trajectories eventually always reach one of two points depending on the initial conditions.\n", "- The two points where the trajectories converge are the intersection of the two nullcline curves." ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "## Think! 2.2: Analyzing the vector field\n", "\n", "There are, in total, three intersection points, meaning that the system has three fixed points.\n", "\n", "1. One of the fixed points (the one in the middle) is never the final state of a trajectory. Why is that?\n", "2. Why do the arrows tend to get smaller as they approach the fixed points?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {} }, "outputs": [], "source": [ "# to_remove explanation\n", "\"\"\"\n", "1. Because the middle fixed point is unstable. Trajectories only point to stable\n", "fixed points.\n", "\n", "2. The slope of dr/dt determines the speed at which the system states will evolve.\n", "At the nullclines dr_e/dt = 0 and/or dr_i/dt = 0. That means as we move close to\n", "the nullclines the dr/dt becomes smaller and smaller. Therefore the system state\n", "evolves slower and slower as we approach the nullcline.\n", "\"\"\";" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Submit your feedback\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# @title Submit your feedback\n", "content_review(f\"{feedback_prefix}_Analyzing_the_vector_field_Discussion\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "---\n", "# Summary\n", "\n", "*Estimated timing of tutorial: 1 hour, 35 minutes*\n", "\n", "Congratulations! You have finished the fourth day of the second week of the neuromatch academy! Here, you learned how to simulate a rate based model consisting of excitatory and inhibitory population of neurons.\n", "\n", "In the last tutorial on dynamical neuronal networks you learned to:\n", "- Implement and simulate a 2D system composed of an E and an I population of neurons using the **Wilson-Cowan** model\n", "- Plot the frequency-current (F-I) curves for both populations\n", "- Examine the behavior of the system using phase **plane analysis**, **vector fields**, and **nullclines**.\n", "\n", "Do you have more time? Have you finished early? We have more fun material for you!\n", "\n", "In the bonus tutorial, there are some, more advanced concepts on dynamical systems:\n", "\n", "- You will learn how to find the fixed points on such a system, and to investigate its stability by linearizing its dynamics and examining the **Jacobian matrix**.\n", "- You will identify conditions under which the Wilson-Cowan model can exhibit oscillations.\n", "\n", "If you need even more, there are two applications of the Wilson-Cowan model:\n", "\n", "- Visualization of an Inhibition-stabilized network\n", "- Simulation of working memory" ] } ], "metadata": { "colab": { "collapsed_sections": [], "include_colab_link": true, "name": "W2D4_Tutorial2", "provenance": [], "toc_visible": true }, "kernel": { "display_name": "Python 3", "language": "python", "name": "python3" }, "kernelspec": { "display_name": "Python 3", "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.9.17" } }, "nbformat": 4, "nbformat_minor": 0 }