{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"execution": {},
"id": "view-in-github"
},
"source": [
"
"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"# Tutorial 2: \"How\" models\n",
"\n",
"**Week 1, Day 1: Model Types**\n",
"\n",
"**By Neuromatch Academy**\n",
"\n",
"__Content creators:__ Matt Laporte, Byron Galbraith, Konrad Kording\n",
"\n",
"__Content reviewers:__ Dalin Guo, Aishwarya Balwani, Madineh Sarvestani, Maryam Vaziri-Pashkam, Michael Waskom, Ella Batty\n",
"\n",
"__Production editors:__ Gagana B, Spiros Chavlis"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"___\n",
"# Tutorial Objectives\n",
"\n",
"*Estimated timing of tutorial: 45 minutes*\n",
"\n",
"This is tutorial 2 of a 3-part series on different flavors of models used to understand neural data. In this tutorial we will explore models that can potentially explain *how* the spiking data we have observed is produced\n",
"\n",
"To understand the mechanisms that give rise to the neural data we save in Tutorial 1, we will build simple neuronal models and compare their spiking response to real data. We will:\n",
"- Write code to simulate a simple \"leaky integrate-and-fire\" neuron model\n",
"- Make the model more complicated — but also more realistic — by adding more physiologically-inspired details"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {}
},
"outputs": [],
"source": [
"\n",
"# @title Tutorial slides\n",
"# @markdown These are the slides for the videos in all tutorials today\n",
"from IPython.display import IFrame\n",
"link_id = \"6dxwe\"\n",
"print(f\"If you want to download the slides: https://osf.io/download/{link_id}/\")\n",
"IFrame(src=f\"https://mfr.ca-1.osf.io/render?url=https://osf.io/{link_id}/?direct%26mode=render%26action=download%26mode=render\", width=854, height=480)"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"---\n",
"# Setup"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {}
},
"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 = \"W1D1_T2\""
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "both",
"execution": {}
},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from scipy import stats"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {}
},
"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",
"%matplotlib inline\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": null,
"metadata": {
"cellView": "form",
"execution": {}
},
"outputs": [],
"source": [
"# @title Plotting Functions\n",
"\n",
"def histogram(counts, bins, vlines=(), ax=None, ax_args=None, **kwargs):\n",
" \"\"\"Plot a step histogram given counts over bins.\"\"\"\n",
" if ax is None:\n",
" _, ax = plt.subplots()\n",
"\n",
" # duplicate the first element of `counts` to match bin edges\n",
" counts = np.insert(counts, 0, counts[0])\n",
"\n",
" ax.fill_between(bins, counts, step=\"pre\", alpha=0.4, **kwargs) # area shading\n",
" ax.plot(bins, counts, drawstyle=\"steps\", **kwargs) # lines\n",
"\n",
" for x in vlines:\n",
" ax.axvline(x, color='r', linestyle='dotted') # vertical line\n",
"\n",
" if ax_args is None:\n",
" ax_args = {}\n",
"\n",
" # heuristically set max y to leave a bit of room\n",
" ymin, ymax = ax_args.get('ylim', [None, None])\n",
" if ymax is None:\n",
" ymax = np.max(counts)\n",
" if ax_args.get('yscale', 'linear') == 'log':\n",
" ymax *= 1.5\n",
" else:\n",
" ymax *= 1.1\n",
" if ymin is None:\n",
" ymin = 0\n",
"\n",
" if ymax == ymin:\n",
" ymax = None\n",
"\n",
" ax_args['ylim'] = [ymin, ymax]\n",
"\n",
" ax.set(**ax_args)\n",
" ax.autoscale(enable=False, axis='x', tight=True)\n",
"\n",
"def plot_neuron_stats(v, spike_times):\n",
" fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(12, 5))\n",
"\n",
" # membrane voltage trace\n",
" ax1.plot(v[0:100])\n",
" ax1.set(xlabel='Time', ylabel='Voltage')\n",
" # plot spike events\n",
" for x in spike_times:\n",
" if x >= 100:\n",
" break\n",
" ax1.axvline(x, color='red')\n",
"\n",
" # ISI distribution\n",
" if len(spike_times)>1:\n",
" isi = np.diff(spike_times)\n",
" n_bins = np.arange(isi.min(), isi.max() + 2) - .5\n",
" counts, bins = np.histogram(isi, n_bins)\n",
" vlines = []\n",
" if len(isi) > 0:\n",
" vlines = [np.mean(isi)]\n",
" xmax = max(20, int(bins[-1])+5)\n",
" histogram(counts, bins, vlines=vlines, ax=ax2, ax_args={\n",
" 'xlabel': 'Inter-spike interval',\n",
" 'ylabel': 'Number of intervals',\n",
" 'xlim': [0, xmax]\n",
" })\n",
" else:\n",
" ax2.set(xlabel='Inter-spike interval',\n",
" ylabel='Number of intervals')\n",
" plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"---\n",
"# \"How\" models"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {}
},
"outputs": [],
"source": [
"# @title Video 1: \"How\" models\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', 'PpnagITsb3E'), ('Bilibili', 'BV1yV41167Di')]\n",
"tab_contents = display_videos(video_ids, W=854, H=480)\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": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {}
},
"outputs": [],
"source": [
"# @title Submit your feedback\n",
"content_review(f\"{feedback_prefix}_How_models_Video\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"---\n",
"# Section 1: The Linear Integrate-and-Fire Neuron\n",
"\n",
"*Note that we'll be using some math to describe our leaky integrate-and-fire neurons. If you need to quickly check the notation (i.e. what a certain variable name stands for), you can look at the notation section after the summary!*\n",
"\n",
"How does a neuron spike?\n",
"\n",
"A neuron charges and discharges an electric field across its cell membrane. The state of this electric field can be described by the _membrane potential_. The membrane potential rises due to excitation of the neuron, and when it reaches a threshold a spike occurs. The potential resets, and must rise to a threshold again before the next spike occurs.\n",
"\n",
"One of the simplest models of spiking neuron behavior is the linear integrate-and-fire model neuron. In this model, the neuron increases its membrane potential $V_m$ over time in response to excitatory input currents $I$ scaled by some factor $\\alpha$:\n",
"\n",
"\\begin{equation}\n",
"dV_m = {\\alpha}I\n",
"\\end{equation}\n",
"\n",
"Once $V_m$ reaches a threshold value a spike is produced, $V_m$ is reset to a starting value, and the process continues.\n",
"\n",
"Here, we will take the starting and threshold potentials as $0$ and $1$, respectively. So, for example, if $\\alpha I=0.1$ is constant---that is, the input current is constant---then $dV_m=0.1$, and at each timestep the membrane potential $V_m$ increases by $0.1$ until after $(1-0)/0.1 = 10$ timesteps it reaches the threshold and resets to $V_m=0$, and so on.\n",
"\n",
"Note that we define the membrane potential $V_m$ as a scalar: a single real (or floating point) number. However, a biological neuron's membrane potential will not be exactly constant at all points on its cell membrane at a given time. We could capture this variation with a more complex model (e.g. with more numbers). Do we need to?\n",
"\n",
"The proposed model is a 1D simplification. There are many details we could add to it, to preserve different parts of the complex structure and dynamics of a real neuron. If we were interested in small or local changes in the membrane potential, our 1D simplification could be a problem. However, we'll assume an idealized \"point\" neuron model for our current purpose.\n",
"\n",
"#### Spiking Inputs\n",
"\n",
"Given our simplified model for the neuron dynamics, we still need to consider what form the input $I$ will take. How should we specify the firing behavior of the presynaptic neuron(s) providing the inputs to our model neuron?\n",
"\n",
"Unlike in the simple example above, where $\\alpha I=0.1$, the input current is generally not constant. Physical inputs tend to vary with time. We can describe this variation with a distribution.\n",
"\n",
"We'll assume the input current $I$ over a timestep is due to equal contributions from a non-negative ($\\ge 0$) integer number of input spikes arriving in that timestep. Our model neuron might integrate currents from 3 input spikes in one timestep, and 7 spikes in the next timestep. We should see similar behavior when sampling from our distribution.\n",
"\n",
"Given no other information about the input neurons, we will also assume that the distribution has a mean (i.e. mean rate, or number of spikes received per timestep), and that the spiking events of the input neuron(s) are independent in time. Are these reasonable assumptions in the context of real neurons?\n",
"\n",
"A suitable distribution given these assumptions is the Poisson distribution, which we'll use to model $I$:\n",
"\n",
"\\begin{equation}\n",
"I \\sim \\mathrm{Poisson}(\\lambda)\n",
"\\end{equation}\n",
"\n",
"where $\\lambda$ is the mean of the distribution: the average rate of spikes received per timestep."
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"## Coding Exercise 1: Compute $dV_m$\n",
"\n",
"For your first exercise, you will write the code to compute the change in voltage $dV_m$ (per timestep) of the linear integrate-and-fire model neuron. The rest of the code to handle numerical integration is provided for you, so you just need to fill in a definition for `dv` in the `lif_neuron` function below. The value of $\\lambda$ for the Poisson random variable is given by the function argument `rate`.\n",
"\n",
"\n",
"\n",
"The [`scipy.stats`](https://docs.scipy.org/doc/scipy/reference/stats.html) package is a great resource for working with and sampling from various probability distributions. We will use the `scipy.stats.poisson` class and its method `rvs` to produce Poisson-distributed random samples. In this tutorial, we have imported this package with the alias `stats`, so you should refer to it in your code as `stats.poisson`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {}
},
"outputs": [],
"source": [
"def lif_neuron(n_steps=1000, alpha=0.01, rate=10):\n",
" \"\"\" Simulate a linear integrate-and-fire neuron.\n",
"\n",
" Args:\n",
" n_steps (int): The number of time steps to simulate the neuron's activity.\n",
" alpha (float): The input scaling factor\n",
" rate (int): The mean rate of incoming spikes\n",
"\n",
" \"\"\"\n",
" # Precompute Poisson samples for speed\n",
" exc = stats.poisson(rate).rvs(n_steps)\n",
"\n",
" # Initialize voltage and spike storage\n",
" v = np.zeros(n_steps)\n",
" spike_times = []\n",
"\n",
" ################################################################################\n",
" # Students: compute dv, then comment out or remove the next line\n",
" raise NotImplementedError(\"Exercise: compute the change in membrane potential\")\n",
" ################################################################################\n",
"\n",
" # Loop over steps\n",
" for i in range(1, n_steps):\n",
"\n",
" # Update v\n",
" dv = ...\n",
" v[i] = v[i-1] + dv\n",
"\n",
" # If spike happens, reset voltage and record\n",
" if v[i] > 1:\n",
" spike_times.append(i)\n",
" v[i] = 0\n",
"\n",
" return v, spike_times\n",
"\n",
"# Set random seed (for reproducibility)\n",
"np.random.seed(12)\n",
"\n",
"# Model LIF neuron\n",
"v, spike_times = lif_neuron()\n",
"\n",
"# Visualize\n",
"plot_neuron_stats(v, spike_times)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "both",
"execution": {}
},
"outputs": [],
"source": [
"# to_remove solution\n",
"\n",
"def lif_neuron(n_steps=1000, alpha=0.01, rate=10):\n",
" \"\"\" Simulate a linear integrate-and-fire neuron.\n",
"\n",
" Args:\n",
" n_steps (int): The number of time steps to simulate the neuron's activity.\n",
" alpha (float): The input scaling factor\n",
" rate (int): The mean rate of incoming spikes\n",
"\n",
" \"\"\"\n",
" # Precompute Poisson samples for speed\n",
" exc = stats.poisson(rate).rvs(n_steps)\n",
"\n",
" # Initialize voltage and spike storage\n",
" v = np.zeros(n_steps)\n",
" spike_times = []\n",
"\n",
" # Loop over steps\n",
" for i in range(1, n_steps):\n",
"\n",
" # Update v\n",
" dv = alpha * exc[i]\n",
" v[i] = v[i-1] + dv\n",
"\n",
" # If spike happens, reset voltage and record\n",
" if v[i] > 1:\n",
" spike_times.append(i)\n",
" v[i] = 0\n",
"\n",
" return v, spike_times\n",
"\n",
"# Set random seed (for reproducibility)\n",
"np.random.seed(12)\n",
"\n",
"# Model LIF neuron\n",
"v, spike_times = lif_neuron()\n",
"\n",
"# Visualize\n",
"with plt.xkcd():\n",
" plot_neuron_stats(v, spike_times)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {}
},
"outputs": [],
"source": [
"# @title Submit your feedback\n",
"content_review(f\"{feedback_prefix}_Compute_dVm_Exercise\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"## Interactive Demo 1: Linear-IF neuron\n",
"Like last time, you can now explore how various parameters of the LIF model influence the ISI distribution. Specifically, you can vary `alpha`, which is input scaling factor, and `rate`, which is the mean rate of incoming spikes.\n",
"\n",
"\n",
"1. What is the spiking pattern of this model?\n",
"2. What effect does raising or lowering alpha have?\n",
"3. What effect does raising or lowering the rate have?\n",
"4. Does the distribution of ISIs ever look like what you observed in the data in Tutorial 1?\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {}
},
"outputs": [],
"source": [
"# @markdown You don't need to worry about how the code works – but you do need to **run the cell** to enable the sliders.\n",
"\n",
"def _lif_neuron(n_steps=1000, alpha=0.01, rate=10):\n",
" exc = stats.poisson(rate).rvs(n_steps)\n",
" v = np.zeros(n_steps)\n",
" spike_times = []\n",
" for i in range(1, n_steps):\n",
" dv = alpha * exc[i]\n",
" v[i] = v[i-1] + dv\n",
" if v[i] > 1:\n",
" spike_times.append(i)\n",
" v[i] = 0\n",
" return v, spike_times\n",
"\n",
"@widgets.interact(\n",
" alpha=widgets.FloatLogSlider(0.01, min=-2, max=-1),\n",
" rate=widgets.IntSlider(10, min=5, max=20)\n",
")\n",
"def plot_lif_neuron(alpha=0.01, rate=10):\n",
" v, spike_times = _lif_neuron(2000, alpha, rate)\n",
"\n",
" plot_neuron_stats(v, spike_times)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {}
},
"outputs": [],
"source": [
"# to_remove explanation\n",
"\n",
"\"\"\"\n",
"1. The spikes are pretty predictable - they occur at fairly regular intervals.\n",
"2. Raising alpha increases the rate of spikes (and thus decreases the inter-spike intervals)\n",
"3. Raising the rate increases the rate of spikes (and thus decreases the inter-spike intervals)\n",
"These two facts make sense because they both control the size of the input affecting the neuron - more\n",
"input means higher responses.\n",
"4. No, the ISIs look pretty symmetrical and Gaussian (not exponential)\n",
"\"\"\";"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {}
},
"outputs": [],
"source": [
"# @title Submit your feedback\n",
"content_review(f\"{feedback_prefix}_Linear_IF_neuron_Interactive_Demo_and_Discussion\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {}
},
"outputs": [],
"source": [
"# @title Video 2: Linear-IF models\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', 'QBD7kulhg4U'), ('Bilibili', 'BV1iZ4y1u7en')]\n",
"tab_contents = display_videos(video_ids, W=854, H=480)\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": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {}
},
"outputs": [],
"source": [
"# @title Submit your feedback\n",
"content_review(f\"{feedback_prefix}_Linear_IF_models_Video\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"---\n",
"# Section 2: Inhibitory signals\n",
"\n",
"*Estimated timing to here from start of tutorial: 20 min*\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"\n",
"Our linear integrate-and-fire neuron from the previous section was indeed able to produce spikes. However, our ISI histogram doesn't look much like empirical ISI histograms seen in Tutorial 1, which had an exponential-like shape. What is our model neuron missing, given that it doesn't behave like a real neuron?\n",
"\n",
"In the previous model we only considered excitatory behavior -- the only way the membrane potential could decrease was upon a spike event. We know, however, that there are other factors that can drive $V_m$ down. First is the natural tendency of the neuron to return to some steady state or resting potential. We can update our previous model as follows:\n",
"\n",
"\\begin{equation}\n",
"dV_m = -{\\beta}V_m + {\\alpha}I\n",
"\\end{equation}\n",
"\n",
"where $V_m$ is the current membrane potential and $\\beta$ is some leakage factor. This is a basic form of the popular Leaky Integrate-and-Fire model neuron (for a more detailed discussion of the LIF Neuron, see Bonus Section 2 and the Biological Neuron Models day later in this course).\n",
"\n",
"We also know that in addition to excitatory presynaptic neurons, we can have inhibitory presynaptic neurons as well. We can model these inhibitory neurons with another Poisson random variable:\n",
"\n",
"\\begin{align}\n",
"I &= I_{\\mathrm{exc}} - I_{\\mathrm{inh}} \\\\\n",
"I_{\\mathrm{exc}} &\\sim \\mathrm{Poisson}(\\lambda_{\\mathrm{exc}}) \\\\\n",
"I_{\\mathrm{inh}} &\\sim \\mathrm{Poisson}(\\lambda_{\\mathrm{inh}})\n",
"\\end{align}\n",
"\n",
"where $\\lambda_{\\mathrm{exc}}$ and $\\lambda_{\\mathrm{inh}}$ are the average spike rates (per timestep) of the excitatory and inhibitory presynaptic neurons, respectively."
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"### Coding Exercise 2: Compute $dV_m$ with inhibitory signals\n",
"\n",
"For your second exercise, you will again write the code to compute the change in voltage $dV_m$, though now of the LIF model neuron described above. Like last time, the rest of the code needed to handle the neuron dynamics are provided for you, so you just need to fill in a definition for `dv` below.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {}
},
"outputs": [],
"source": [
"def lif_neuron_inh(n_steps=1000, alpha=0.5, beta=0.1, exc_rate=10, inh_rate=10):\n",
" \"\"\" Simulate a simplified leaky integrate-and-fire neuron with both excitatory\n",
" and inhibitory inputs.\n",
"\n",
" Args:\n",
" n_steps (int): The number of time steps to simulate the neuron's activity.\n",
" alpha (float): The input scaling factor\n",
" beta (float): The membrane potential leakage factor\n",
" exc_rate (int): The mean rate of the incoming excitatory spikes\n",
" inh_rate (int): The mean rate of the incoming inhibitory spikes\n",
" \"\"\"\n",
"\n",
" # precompute Poisson samples for speed\n",
" exc = stats.poisson(exc_rate).rvs(n_steps)\n",
" inh = stats.poisson(inh_rate).rvs(n_steps)\n",
"\n",
" v = np.zeros(n_steps)\n",
" spike_times = []\n",
"\n",
" ###############################################################################\n",
" # Students: compute dv, then comment out or remove the next line\n",
" raise NotImplementedError(\"Exercise: compute the change in membrane potential\")\n",
" ################################################################################\n",
"\n",
" for i in range(1, n_steps):\n",
"\n",
" dv = ...\n",
"\n",
" v[i] = v[i-1] + dv\n",
" if v[i] > 1:\n",
" spike_times.append(i)\n",
" v[i] = 0\n",
"\n",
" return v, spike_times\n",
"\n",
"# Set random seed (for reproducibility)\n",
"np.random.seed(12)\n",
"\n",
"# Model LIF neuron\n",
"v, spike_times = lif_neuron_inh()\n",
"\n",
"# Visualize\n",
"plot_neuron_stats(v, spike_times)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "both",
"execution": {}
},
"outputs": [],
"source": [
"# to_remove solution\n",
"def lif_neuron_inh(n_steps=1000, alpha=0.5, beta=0.1, exc_rate=10, inh_rate=10):\n",
" \"\"\" Simulate a simplified leaky integrate-and-fire neuron with both excitatory\n",
" and inhibitory inputs.\n",
"\n",
" Args:\n",
" n_steps (int): The number of time steps to simulate the neuron's activity.\n",
" alpha (float): The input scaling factor\n",
" beta (float): The membrane potential leakage factor\n",
" exc_rate (int): The mean rate of the incoming excitatory spikes\n",
" inh_rate (int): The mean rate of the incoming inhibitory spikes\n",
" \"\"\"\n",
" # precompute Poisson samples for speed\n",
" exc = stats.poisson(exc_rate).rvs(n_steps)\n",
" inh = stats.poisson(inh_rate).rvs(n_steps)\n",
"\n",
" v = np.zeros(n_steps)\n",
" spike_times = []\n",
" for i in range(1, n_steps):\n",
"\n",
" dv = -beta * v[i-1] + alpha * (exc[i] - inh[i])\n",
"\n",
" v[i] = v[i-1] + dv\n",
" if v[i] > 1:\n",
" spike_times.append(i)\n",
" v[i] = 0\n",
"\n",
" return v, spike_times\n",
"\n",
"# Set random seed (for reproducibility)\n",
"np.random.seed(12)\n",
"\n",
"# Model LIF neuron\n",
"v, spike_times = lif_neuron_inh()\n",
"\n",
"# Visualize\n",
"with plt.xkcd():\n",
" plot_neuron_stats(v, spike_times)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {}
},
"outputs": [],
"source": [
"# @title Submit your feedback\n",
"content_review(f\"{feedback_prefix}_Compute_dVm_with_inhibitory_signals_Exercise\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"## Interactive Demo 2: LIF + inhibition neuron\n",
"\n",
"As in Interactive Demo 1, you can play with the parameters of the input to our LIF neuron and visualize what happens. Here, in addition to controlling `alpha`, which scales the input, you can also control `beta`, which is a leakage factor on the voltage, `exc_rate`, which is the mean rate of the excitatory presynaptic neurons, and `inh_rate`, which is the mean rate of the inhibitory presynaptic neurons.\n",
"\n",
"1. What effect does raising the excitatory rate have?\n",
"2. What effect does raising the inhibitory rate have?\n",
"3. What if you raise both the excitatory and inhibitory rate?\n",
"4. Does the distribution of ISIs ever look like what you observed in the data in Tutorial 1?\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {}
},
"outputs": [],
"source": [
"#@title\n",
"#@markdown **Run the cell** to enable the sliders.\n",
"def _lif_neuron_inh(n_steps=1000, alpha=0.5, beta=0.1, exc_rate=10, inh_rate=10):\n",
" \"\"\" Simulate a simplified leaky integrate-and-fire neuron with both excitatory\n",
" and inhibitory inputs.\n",
"\n",
" Args:\n",
" n_steps (int): The number of time steps to simulate the neuron's activity.\n",
" alpha (float): The input scaling factor\n",
" beta (float): The membrane potential leakage factor\n",
" exc_rate (int): The mean rate of the incoming excitatory spikes\n",
" inh_rate (int): The mean rate of the incoming inhibitory spikes\n",
" \"\"\"\n",
" # precompute Poisson samples for speed\n",
" exc = stats.poisson(exc_rate).rvs(n_steps)\n",
" inh = stats.poisson(inh_rate).rvs(n_steps)\n",
"\n",
" v = np.zeros(n_steps)\n",
" spike_times = []\n",
" for i in range(1, n_steps):\n",
" dv = -beta * v[i-1] + alpha * (exc[i] - inh[i])\n",
" v[i] = v[i-1] + dv\n",
" if v[i] > 1:\n",
" spike_times.append(i)\n",
" v[i] = 0\n",
"\n",
" return v, spike_times\n",
"\n",
"@widgets.interact(alpha=widgets.FloatLogSlider(0.5, min=-1, max=1),\n",
" beta=widgets.FloatLogSlider(0.1, min=-1, max=0),\n",
" exc_rate=widgets.IntSlider(12, min=10, max=20),\n",
" inh_rate=widgets.IntSlider(12, min=10, max=20))\n",
"def plot_lif_neuron(alpha=0.5, beta=0.1, exc_rate=10, inh_rate=10):\n",
" v, spike_times = _lif_neuron_inh(2000, alpha, beta, exc_rate, inh_rate)\n",
" plot_neuron_stats(v, spike_times)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {}
},
"outputs": [],
"source": [
"# to_remove explanation\n",
"\n",
"\"\"\"\n",
"1. Raising the excitatory rate while keeping inhibitory rate the same results in\n",
"more frequent firing and shorter ISIs. This makes intuitive sense - more excitatory input\n",
"means higher responses.\n",
"2. Raising the inhibitory rate while keeping the excitatory rate the same results\n",
"in less frequent firing and longer ISIs. This makes intuitive sense - more inhibitory input\n",
"means lower responses.\n",
"3. If the excitatory and inhibitory rates are equal, the average inter-spike interval\n",
"stays the same even if you raise both. This is because they balance each other out.\n",
"4. Yes, these ISIs look more exponential, like what we observed.\n",
"\"\"\";"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {}
},
"outputs": [],
"source": [
"# @title Submit your feedback\n",
"content_review(f\"{feedback_prefix}_LIF_and_inhibition_neuron_Interactive_Demo_and_Discussion\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {}
},
"outputs": [],
"source": [
"# @title Video 3: LIF + inhibition\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', 'Aq7JrxRkn2w'), ('Bilibili', 'BV1nV41167mS')]\n",
"tab_contents = display_videos(video_ids, W=854, H=480)\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": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {}
},
"outputs": [],
"source": [
"# @title Submit your feedback\n",
"content_review(f\"{feedback_prefix}_LIF_and_inhibition_Video\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"---\n",
"# Section 3: Reflecting on how models\n",
"\n",
"*Estimated timing to here from start of tutorial: 35 min*"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"## Think! 3: Reflecting on how models\n",
"\n",
"Please discuss the following questions for around 10 minutes with your group:\n",
"\n",
"- Have you seen how models before?\n",
"- Have you ever done one?\n",
"- Why are how models useful?\n",
"- When are they possible? Does your field have how models?\n",
"- What do we learn from constructing them?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"cellView": "form",
"execution": {}
},
"outputs": [],
"source": [
"# @title Submit your feedback\n",
"content_review(f\"{feedback_prefix}_Reflecting_on_how_models_Discussion\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"# Summary\n",
"\n",
"*Estimated timing of tutorial: 45 minutes*\n",
"\n",
"In this tutorial we gained some intuition for the mechanisms that produce the observed behavior in our real neural data. First, we built a simple neuron model with excitatory input and saw that its behavior, measured using the ISI distribution, did not match our real neurons. We then improved our model by adding leakiness and inhibitory input. The behavior of this balanced model was much closer to the real neural data."
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"---\n",
"# Notation\n",
"\n",
"\\begin{align}\n",
"V_m &\\quad \\text{membrane potential} \\\\\n",
"dV_m &\\quad \\text{change in membrane potential}\\\\\n",
"C_m &\\quad \\text{membrane capacitance}\\\\\n",
"I &\\quad \\text{input current}\\\\\n",
"R_m &\\quad \\text{membrane resistance}:\\\\\n",
"V_\\mathrm{rest} &\\quad \\text{resting potential}\\\\\n",
"\\alpha &\\quad \\text{scaling factor for input current}\\\\\n",
"\\beta &\\quad \\text{leakage factor}\\\\\n",
"\\lambda &\\quad \\text{average spike rate}\\\\\n",
"\\lambda_\\mathrm{exc} &\\quad \\text{average spike rate for excitatory neurons}\\\\\n",
"\\lambda_\\mathrm{inh} &\\quad \\text{average spike rate for inhibitory neurons}\\\\\n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"---\n",
"# Bonus"
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"## Bonus Section 1: Why do neurons spike?\n",
"\n",
"A neuron stores energy in an electric field across its cell membrane, by controlling the distribution of charges (ions) on either side of the membrane. This energy is rapidly discharged to generate a spike when the field potential (or membrane potential) crosses a threshold. The membrane potential may be driven toward or away from this threshold, depending on inputs from other neurons: excitatory or inhibitory, respectively. The membrane potential tends to revert to a resting potential, for example due to the leakage of ions across the membrane, so that reaching the spiking threshold depends not only on the amount of input ever received following the last spike, but also the timing of the inputs.\n",
"\n",
"The storage of energy by maintaining a field potential across an insulating membrane can be modeled by a capacitor. The leakage of charge across the membrane can be modeled by a resistor. This is the basis for the leaky integrate-and-fire neuron model."
]
},
{
"cell_type": "markdown",
"metadata": {
"execution": {}
},
"source": [
"## Bonus Section 2: The LIF Model Neuron\n",
"\n",
"The full equation for the LIF neuron is\n",
"\n",
"\\begin{equation}\n",
"C_{m}\\frac{dV_m}{dt} = -(V_m - V_{\\mathrm{rest}})/R_{m} + I\n",
"\\end{equation}\n",
"\n",
"where $C_m$ is the membrane capacitance, $R_m$ is the membrane resistance, $V_{\\mathrm{rest}}$ is the resting potential, and $I$ is some input current (from other neurons, an electrode, etc.).\n",
"\n",
"In our above examples we set many of these parameters to convenient values ($C_m = R_m = dt = 1$, $V_{\\mathrm{rest}} = 0$) to focus more on the general behavior of the model. However, these too can be manipulated to achieve different dynamics, or to ensure the dimensions of the problem are preserved between simulation units and experimental units (e.g. with $V_m$ given in millivolts, $R_m$ in megaohms, $t$ in milliseconds)."
]
}
],
"metadata": {
"colab": {
"collapsed_sections": [],
"include_colab_link": true,
"name": "W1D1_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
}