{ "cells": [ { "cell_type": "markdown", "metadata": { "colab_type": "text", "execution": {}, "id": "view-in-github" }, "source": [ "\"Open   \"Open" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "# Tutorial 1: \"What\" 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\n", "\n", "
\n", "\n", "We would like to acknowledge [Steinmetz _et al._ (2019)](https://www.nature.com/articles/s41586-019-1787-x) for sharing their data, a subset of which is used here." ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "___\n", "# Tutorial Objectives\n", "\n", "*Estimated timing of tutorial: 50 minutes*\n", "\n", "This is tutorial 1 of a 3-part series on different flavors of models used to understand neural data. In this tutorial we will explore 'What' models, used to describe the data. To understand what our data looks like, we will visualize it in different ways. Then we will compare it to simple mathematical models. Specifically, we will:\n", "\n", "- Load a dataset with spiking activity from hundreds of neurons and understand how it is organized\n", "- Make plots to visualize characteristics of the spiking activity across the population\n", "- Compute the distribution of \"inter-spike intervals\" (ISIs) for a single neuron\n", "- Consider several formal models of this distribution's shape and fit them to the data \"by hand\"" ] }, { "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\n", "\n" ] }, { "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_T1\"" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Python requires you to explicitly \"import\" libraries before their functions are available to use. We will always specify our imports at the beginning of each notebook or script." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Tutorial notebooks typically begin with several set-up steps that are hidden from view by default.\n", "\n", "**Important:** Even though the code is hidden, you still need to run it so that the rest of the notebook can work properly. Step through each cell, either by pressing the play button in the upper-left-hand corner or with a keyboard shortcut (`Cmd-Return` on a Mac, `Ctrl-Enter` otherwise). A number will appear inside the brackets (e.g. `[3]`) to tell you that the cell was executed and what order that happened in.\n", "\n", "If you are curious to see what is going on inside each cell, you can double click to expand. Once expanded, double-click the white space to the right of the editor to collapse again." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "pycharm": { "name": "#%%\n" } }, "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 plot_isis(single_neuron_isis):\n", " plt.hist(single_neuron_isis, bins=50, histtype=\"stepfilled\")\n", " plt.axvline(single_neuron_isis.mean(), color=\"orange\", label=\"Mean ISI\")\n", " plt.xlabel(\"ISI duration (s)\")\n", " plt.ylabel(\"Number of spikes\")\n", " plt.legend()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "#@title Data retrieval\n", "#@markdown This cell downloads the example dataset that we will use in this tutorial.\n", "import io\n", "import requests\n", "r = requests.get('https://osf.io/sy5xt/download')\n", "if r.status_code != 200:\n", " print('Failed to download data')\n", "else:\n", " spike_times = np.load(io.BytesIO(r.content), allow_pickle=True)['spike_times']" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "---\n", "# \"What\" models" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @title Video 1: \"What\" 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', 'KgqR_jbjMQg'), ('Bilibili', 'BV1mz4y1X7ot')]\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}_What_models_Video\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "---\n", "\n", "# Section 1: Exploring the Steinmetz dataset\n", "\n", "In this tutorial we will explore the structure of a neuroscience dataset.\n", "\n", "We consider a subset of data from a study of [Steinmetz _et al._ (2019)](https://www.nature.com/articles/s41586-019-1787-x). In this study, Neuropixels probes were implanted in the brains of mice. Electrical potentials were measured by hundreds of electrodes along the length of each probe. Each electrode's measurements captured local variations in the electric field due to nearby spiking neurons. A spike sorting algorithm was used to infer spike times and cluster spikes according to common origin: a single cluster of sorted spikes is causally attributed to a single neuron.\n", "\n", "In particular, a single recording session of spike times and neuron assignments was loaded and assigned to `spike_times` in the preceding setup.\n", "\n", "Typically a dataset comes with some information about its structure. However, this information may be incomplete. You might also apply some transformations or \"pre-processing\" to create a working representation of the data of interest, which might go partly undocumented depending on the circumstances. In any case it is important to be able to use the available tools to investigate unfamiliar aspects of a data structure.\n", "\n", "Let's see what our data looks like..." ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "## Section 1.1: Warming up with `spike_times`" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "What is the Python type of our variable?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "type(spike_times)" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "You should see `numpy.ndarray`, which means that it's a normal NumPy array.\n", "\n", "If you see an error message, it probably means that you did not execute the set-up cells at the top of the notebook. So go ahead and make sure to do that.\n", "\n", "Once everything is running properly, we can ask the next question about the dataset: what's its shape?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "spike_times.shape" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "There are 734 entries in one dimension, and no other dimensions. What is the Python type of the first entry, and what is *its* shape?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "idx = 0\n", "print(\n", " type(spike_times[idx]),\n", " spike_times[idx].shape,\n", " sep=\"\\n\",\n", ")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "It's also a NumPy array with a 1D shape! Why didn't this show up as a second dimension in the shape of `spike_times`? That is, why not `spike_times.shape == (734, 826)`?\n", "\n", "To investigate, let's check another entry." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "idx = 321\n", "print(\n", " type(spike_times[idx]),\n", " spike_times[idx].shape,\n", " sep=\"\\n\",\n", ")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "It's also a 1D NumPy array, but it has a different shape. Checking the NumPy types of the values in these arrays, and their first few elements, we see they are composed of floating point numbers (not another level of `np.ndarray`):" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "i_neurons = [0, 321]\n", "i_print = slice(0, 5)\n", "\n", "for i in i_neurons:\n", " print(\n", " \"Neuron {}:\".format(i),\n", " spike_times[i].dtype,\n", " spike_times[i][i_print],\n", " \"\\n\",\n", " sep=\"\\n\"\n", " )" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Note that this time we've checked the NumPy `dtype` rather than the Python variable type. These two arrays contain floating point numbers (\"floats\") with 32 bits of precision.\n", "\n", "The basic picture is coming together:\n", "- `spike_times` is 1D, its entries are NumPy arrays, and its length is the number of neurons (734): by indexing it, we select a subset of neurons.\n", "- An array in `spike_times` is also 1D and corresponds to a single neuron; its entries are floating point numbers, and its length is the number of spikes attributed to that neuron. By indexing it, we select a subset of spike times for that neuron.\n", "\n", "Visually, you can think of the data structure as looking something like this:\n", "\n", "```\n", "| . . . . . |\n", "| . . . . . . . . |\n", "| . . . |\n", "| . . . . . . . |\n", "```\n", "\n", "Before moving on, we'll calculate and store the number of neurons in the dataset and the number of spikes per neuron:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "n_neurons = len(spike_times)\n", "total_spikes_per_neuron = [len(spike_times_i) for spike_times_i in spike_times]\n", "\n", "print(f\"Number of neurons: {n_neurons}\")\n", "print(f\"Number of spikes for first five neurons: {total_spikes_per_neuron[:5]}\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @title Video 2: Exploring the dataset\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', 'oHwYWUI_o1U'), ('Bilibili', 'BV1Hp4y1S7Au')]\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}_Exploring_the_dataset_Video\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "## Section 1.2: Getting warmer: counting and plotting total spike counts\n", "\n", "As we've seen, the number of spikes over the entire recording is variable between neurons. More generally, some neurons tend to spike more than others in a given period. Let's explore what the distribution of spiking looks like across all the neurons in the dataset." ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Are most neurons \"loud\" or \"quiet\", compared to the average? To see, we'll define bins of constant width in terms of total spikes and count the neurons that fall in each bin. This is known as a \"histogram\".\n", "\n", "You can plot a histogram with the matplotlib function `plt.hist`. If you just need to compute it, you can use the numpy function `np.histogram` instead." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "plt.hist(total_spikes_per_neuron, bins=50, histtype=\"stepfilled\")\n", "plt.xlabel(\"Total spikes per neuron\")\n", "plt.ylabel(\"Number of neurons\");" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Let's see what percentage of neurons have a below-average spike count:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "mean_spike_count = np.mean(total_spikes_per_neuron)\n", "frac_below_mean = (total_spikes_per_neuron < mean_spike_count).mean()\n", "print(f\"{frac_below_mean:2.1%} of neurons are below the mean\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "We can also see this by adding the average spike count to the histogram plot:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "plt.hist(total_spikes_per_neuron, bins=50, histtype=\"stepfilled\")\n", "plt.xlabel(\"Total spikes per neuron\")\n", "plt.ylabel(\"Number of neurons\")\n", "plt.axvline(mean_spike_count, color=\"orange\", label=\"Mean neuron\")\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "This shows that the majority of neurons are relatively \"quiet\" compared to the mean, while a small number of neurons are exceptionally \"loud\": they must have spiked more often to reach a large count.\n", "\n", "### Coding Exercise 1.2: Comparing mean and median neurons\n", "\n", "If the mean neuron is more active than 68% of the population, what does that imply about the relationship between the mean neuron and the median neuron?\n", "\n", "*Exercise objective:* Reproduce the plot above, but add the median neuron.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "#################################################################################\n", "## TODO for students:\n", "# Fill out function and remove\n", "raise NotImplementedError(\"Student exercise: complete histogram plotting with median\")\n", "#################################################################################\n", "\n", "# Compute median spike count\n", "median_spike_count = ... # Hint: Try the function np.median\n", "\n", "# Visualize median, mean, and histogram\n", "plt.hist(..., bins=50, histtype=\"stepfilled\")\n", "plt.axvline(..., color=\"limegreen\", label=\"Median neuron\")\n", "plt.axvline(mean_spike_count, color=\"orange\", label=\"Mean neuron\")\n", "plt.xlabel(\"Total spikes per neuron\")\n", "plt.ylabel(\"Number of neurons\")\n", "plt.legend()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "# to_remove solution\n", "\n", "# Compute median spike count\n", "median_spike_count = np.median(total_spikes_per_neuron) # Hint: Try the function np.median\n", "\n", "# Visualize median, mean, and histogram\n", "with plt.xkcd():\n", " plt.hist(total_spikes_per_neuron, bins=50, histtype=\"stepfilled\")\n", " plt.axvline(median_spike_count, color=\"limegreen\", label=\"Median neuron\")\n", " plt.axvline(mean_spike_count, color=\"orange\", label=\"Mean neuron\")\n", " plt.xlabel(\"Total spikes per neuron\")\n", " plt.ylabel(\"Number of neurons\")\n", " plt.legend()" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "\n", "*Bonus:* The median is the 50th percentile. What about other percentiles? Can you show the interquartile range on the histogram?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @title Submit your feedback\n", "content_review(f\"{feedback_prefix}_Comparing_mean_and_median_neurons_Exercise\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "---\n", "\n", "# Section 2: Visualizing neuronal spiking activity\n", "\n", "*Estimated timing to here from start of tutorial: 15 min*" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "## Section 2.1: Getting a subset of the data\n", "\n", "Now we'll visualize trains of spikes. Because the recordings are long, we will first define a short time interval and restrict the visualization to only the spikes in this interval. We defined a helper function, `restrict_spike_times`, to do this for you. If you call `help()` on the function, it will tell you a little bit about itself:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @markdown Execute this cell for helper function `restrict_spike_times`\n", "def restrict_spike_times(spike_times, interval):\n", " \"\"\"Given a spike_time dataset, restrict to spikes within given interval.\n", "\n", " Args:\n", " spike_times (sequence of np.ndarray): List or array of arrays,\n", " each inner array has spike times for a single neuron.\n", " interval (tuple): Min, max time values; keep min <= t < max.\n", "\n", " Returns:\n", " np.ndarray: like `spike_times`, but only within `interval`\n", " \"\"\"\n", " interval_spike_times = []\n", " for spikes in spike_times:\n", " interval_mask = (spikes >= interval[0]) & (spikes < interval[1])\n", " interval_spike_times.append(spikes[interval_mask])\n", " return np.array(interval_spike_times, object)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "help(restrict_spike_times)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "t_interval = (5, 15) # units are seconds after start of recording\n", "interval_spike_times = restrict_spike_times(spike_times, t_interval)" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Is this a representative interval? What fraction of the total spikes fall in this interval?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "original_counts = sum([len(spikes) for spikes in spike_times])\n", "interval_counts = sum([len(spikes) for spikes in interval_spike_times])\n", "frac_interval_spikes = interval_counts / original_counts\n", "print(f\"{frac_interval_spikes:.2%} of the total spikes are in the interval\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "How does this compare to the ratio between the interval duration and the experiment duration? (What fraction of the total time is in this interval?)\n", "\n", "We can approximate the experiment duration by taking the minimum and maximum spike time in the whole dataset. To do that, we \"concatenate\" all of the neurons into one array and then use `np.ptp` (\"peak-to-peak\") to get the difference between the maximum and minimum value:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "spike_times_flat = np.concatenate(spike_times)\n", "experiment_duration = np.ptp(spike_times_flat)\n", "interval_duration = t_interval[1] - t_interval[0]\n", "\n", "frac_interval_time = interval_duration / experiment_duration\n", "print(f\"{frac_interval_time:.2%} of the total time is in the interval\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "These two values—the fraction of total spikes and the fraction of total time—are similar. This suggests the average spike rate of the neuronal population is not very different in this interval compared to the entire recording.\n", "\n", "## Section 2.2: Plotting spike trains and rasters\n", "\n", "Now that we have a representative subset, we're ready to plot the spikes, using the matplotlib `plt.eventplot` function. Let's look at a single neuron first:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "neuron_idx = 1\n", "plt.eventplot(interval_spike_times[neuron_idx], color=\".2\")\n", "plt.xlabel(\"Time (s)\")\n", "plt.yticks([]);" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "We can also plot multiple neurons. Here are three:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "neuron_idx = [1, 11, 51]\n", "plt.eventplot(interval_spike_times[neuron_idx], color=\".2\")\n", "plt.xlabel(\"Time (s)\")\n", "plt.yticks([]);" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "This makes a \"raster\" plot, where the spikes from each neuron appear in a different row.\n", "\n", "Plotting a large number of neurons can give you a sense for the characteristics in the population. Let's show every 5th neuron that was recorded:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "neuron_idx = np.arange(0, len(spike_times), 5)\n", "plt.eventplot(interval_spike_times[neuron_idx], color=\".2\")\n", "plt.xlabel(\"Time (s)\")\n", "plt.yticks([]);" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "*Question*: How does the information in this plot relate to the histogram of total spike counts that you saw above?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {} }, "outputs": [], "source": [ "# to_remove explanation\n", "\n", "\"\"\"\n", "The above histogram is the distribution (across neurons) of the cumulative sum of\n", "all individual neurons' total spike count, i.e. sum across each row in the raster plot.\n", "\"\"\";" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @title Video 3: Visualizing activity\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', 'QGA5FCW7kkA'), ('Bilibili', 'BV1dt4y1Q7C5')]\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}_Visualizing_activity_Video\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "---\n", "\n", "# Section 3: Inter-spike intervals and their distributions\n", "\n", "*Estimated timing to here from start of tutorial: 25 min*" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Given the ordered arrays of spike times for each neuron in `spike_times`, which we've just visualized, what can we ask next?\n", "\n", "Scientific questions are informed by existing models. So, what knowledge do we already have that can inform questions about this data?\n", "\n", "We know that there are physical constraints on neuron spiking. Spiking costs energy, which the neuron's cellular machinery can only obtain at a finite rate. Therefore neurons should have a refractory period: they can only fire as quickly as their metabolic processes can support, and there is a minimum delay between consecutive spikes of the same neuron.\n", "\n", "More generally, we can ask \"how long does a neuron wait to spike again?\" or \"what is the longest a neuron will wait?\" Can we transform spike times into something else, to address questions like these more directly?\n", "\n", "We can consider the inter-spike times (or interspike intervals: ISIs). These are simply the time differences between consecutive spikes of the same neuron.\n", "\n", "### Exercise 3: Plot the distribution of ISIs for a single neuron\n", "\n", "*Exercise objective:* make a histogram, like we did for spike counts, to show the distribution of ISIs for one of the neurons in the dataset.\n", "\n", "Do this in three steps:\n", "\n", "1. Extract the spike times for one of the neurons\n", "2. Compute the ISIs (the amount of time between spikes, or equivalently, the difference between adjacent spike times)\n", "3. Plot a histogram with the array of individual ISIs" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "def compute_single_neuron_isis(spike_times, neuron_idx):\n", " \"\"\"Compute a vector of ISIs for a single neuron given spike times.\n", "\n", " Args:\n", " spike_times (list of 1D arrays): Spike time dataset, with the first\n", " dimension corresponding to different neurons.\n", " neuron_idx (int): Index of the unit to compute ISIs for.\n", "\n", " Returns:\n", " isis (1D array): Duration of time between each spike from one neuron.\n", " \"\"\"\n", " #############################################################################\n", " # Students: Fill in missing code (...) and comment or remove the next line\n", " raise NotImplementedError(\"Exercise: compute single neuron ISIs\")\n", " #############################################################################\n", "\n", " # Extract the spike times for the specified neuron\n", " single_neuron_spikes = ...\n", "\n", " # Compute the ISIs for this set of spikes\n", " # Hint: the function np.diff computes discrete differences along an array\n", " isis = ...\n", "\n", " return isis\n", "\n", "# Compute ISIs\n", "single_neuron_isis = compute_single_neuron_isis(spike_times, neuron_idx=283)\n", "\n", "# Visualize ISIs\n", "plot_isis(single_neuron_isis)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "# to_remove solution\n", "def compute_single_neuron_isis(spike_times, neuron_idx):\n", " \"\"\"Compute a vector of ISIs for a single neuron given spike times.\n", "\n", " Args:\n", " spike_times (list of 1D arrays): Spike time dataset, with the first\n", " dimension corresponding to different neurons.\n", " neuron_idx (int): Index of the unit to compute ISIs for.\n", "\n", " Returns:\n", " isis (1D array): Duration of time between each spike from one neuron.\n", " \"\"\"\n", " # Extract the spike times for the specified neuron\n", " single_neuron_spikes = spike_times[neuron_idx]\n", "\n", " # Compute the ISIs for this set of spikes\n", " # Hint: the function np.diff computes discrete differences along an array\n", " isis = np.diff(single_neuron_spikes)\n", "\n", " return isis\n", "\n", "# Compute ISIs\n", "single_neuron_isis = compute_single_neuron_isis(spike_times, neuron_idx=283)\n", "\n", "# Visualize ISIs\n", "with plt.xkcd():\n", " plot_isis(single_neuron_isis)" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "In general, the shorter ISIs are predominant, with counts decreasing rapidly (and smoothly, more or less) with increasing ISI. However, counts also rapidly decrease to zero with _decreasing_ ISI below the maximum of the distribution (8-11 ms). The absence of these very low ISIs agrees with the refractory period hypothesis: the neuron cannot fire quickly enough to populate this region of the ISI distribution.\n", "\n", "Check the distributions of some other neurons. To resolve various features of the distributions, you might need to play with the value of `bins` in the call to `plt.hist`. Using too few bins might smooth over interesting details, but if you use too many bins, the random variability will start to dominate.\n", "\n", "You might also want to restrict the range to see the shape of the distribution when focusing on relatively short or long ISIs. *Hint:* `plt.hist` takes a `range` argument" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @title Submit your feedback\n", "content_review(f\"{feedback_prefix}_ISIs_and_their_distributions_Exercise\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "---\n", "\n", "# Section 4: What is the functional form of an ISI distribution?\n", "\n", "*Estimated timing to here from start of tutorial: 35 min*" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @title Video 4: ISI distribution\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', 'DHhM80MOTe8'), ('Bilibili', 'BV1ov411B7Pm')]\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}_ISI_distribution_Video\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "The ISI histograms seem to follow continuous, monotonically decreasing functions above their maxima. The function is clearly non-linear. Could it belong to a single family of functions?\n", "\n", "To motivate the idea of using a mathematical function to explain physiological phenomena, let's define a few different function forms that we might expect the relationship to follow: exponential, inverse, and linear." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {}, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "def exponential(xs, scale, rate, x0):\n", " \"\"\"A simple parameterized exponential function, applied element-wise.\n", "\n", " Args:\n", " xs (np.ndarray or float): Input(s) to the function.\n", " scale (float): Linear scaling factor.\n", " rate (float): Exponential growth (positive) or decay (negative) rate.\n", " x0 (float): Horizontal offset.\n", "\n", " \"\"\"\n", " ys = scale * np.exp(rate * (xs - x0))\n", " return ys\n", "\n", "def inverse(xs, scale, x0):\n", " \"\"\"A simple parameterized inverse function (`1/x`), applied element-wise.\n", "\n", " Args:\n", " xs (np.ndarray or float): Input(s) to the function.\n", " scale (float): Linear scaling factor.\n", " x0 (float): Horizontal offset.\n", "\n", " \"\"\"\n", " ys = scale / (xs - x0)\n", " return ys\n", "\n", "def linear(xs, slope, y0):\n", " \"\"\"A simple linear function, applied element-wise.\n", "\n", " Args:\n", " xs (np.ndarray or float): Input(s) to the function.\n", " slope (float): Slope of the line.\n", " y0 (float): y-intercept of the line.\n", "\n", " \"\"\"\n", " ys = slope * xs + y0\n", " return ys" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "### Interactive Demo 4: ISI functions explorer\n", "\n", "Here is an interactive demo where you can vary the parameters of these functions and see how well the resulting outputs correspond to the data. Adjust the parameters by moving the sliders and see how close you can get the lines to follow the falling curve of the histogram. This will give you a taste of what you're trying to do when you *fit a model* to data.\n", "\n", "\"Interactive demo\" cells have hidden code that defines an interface where you can play with the parameters of some function using sliders. 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", "- Which type of function (exponential/inverse/linear) can you make match the data best?\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {}, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "#@markdown Be sure to run this cell to enable the demo\n", "# Don't worry about understanding this code! It's to setup an interactive plot.\n", "single_neuron_idx = 283\n", "single_neuron_spikes = spike_times[single_neuron_idx]\n", "single_neuron_isis = np.diff(single_neuron_spikes)\n", "\n", "counts, edges = np.histogram(\n", " single_neuron_isis,\n", " bins=50,\n", " range=(0, single_neuron_isis.max())\n", ")\n", "\n", "functions = dict(\n", " exponential=exponential,\n", " inverse=inverse,\n", " linear=linear,\n", ")\n", "\n", "colors = dict(\n", " exponential=\"C1\",\n", " inverse=\"C2\",\n", " linear=\"C4\",\n", ")\n", "\n", "@widgets.interact(\n", " exp_scale=widgets.FloatSlider(1000, min=0, max=20000, step=250),\n", " exp_rate=widgets.FloatSlider(-10, min=-200, max=50, step=1),\n", " exp_x0=widgets.FloatSlider(0.1, min=-0.5, max=0.5, step=0.005),\n", " inv_scale=widgets.FloatSlider(1000, min=0, max=3e2, step=10),\n", " inv_x0=widgets.FloatSlider(0, min=-0.2, max=0.2, step=0.01),\n", " lin_slope=widgets.FloatSlider(-1e5, min=-6e5, max=1e5, step=10000),\n", " lin_y0=widgets.FloatSlider(10000, min=0, max=4e4, step=1000),\n", ")\n", "def fit_plot(\n", " exp_scale=1000, exp_rate=-10, exp_x0=0.1,\n", " inv_scale=1000, inv_x0=0,\n", " lin_slope=-1e5, lin_y0=2000,\n", "):\n", " \"\"\"Helper function for plotting function fits with interactive sliders.\"\"\"\n", " func_params = dict(\n", " exponential=(exp_scale, exp_rate, exp_x0),\n", " inverse=(inv_scale, inv_x0),\n", " linear=(lin_slope, lin_y0),\n", " )\n", " f, ax = plt.subplots()\n", " ax.fill_between(edges[:-1], counts, step=\"post\", alpha=.5)\n", " xs = np.linspace(1e-10, edges.max())\n", " for name, function in functions.items():\n", " ys = function(xs, *func_params[name])\n", " ax.plot(xs, ys, lw=3, color=colors[name], label=name);\n", " ax.set(\n", " xlim=(edges.min(), edges.max()),\n", " ylim=(0, counts.max() * 1.1),\n", " xlabel=\"ISI (s)\",\n", " ylabel=\"Number of spikes\",\n", " )\n", " ax.legend()\n", " plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {} }, "outputs": [], "source": [ "# to_remove explanation\n", "\n", "\"\"\"\n", "The exponential function can be made to fit the data much better than the linear\n", "or inverse function.\n", "\"\"\";" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @title Submit your feedback\n", "content_review(f\"{feedback_prefix}_ISI_functions_explorer_Interactive_Demo_and_Discussion\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @title Video 5: Fitting models by hand\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', 'uW2HDk_4-wk'), ('Bilibili', 'BV1w54y1S7Eb')]\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}_Fitting_models_by_hand_Video\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "---\n", "# Section 5: Reflecting on what models\n", "\n", "*Estimated timing to here from start of tutorial: 40 min*" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "## Think! 5: Reflecting on what models\n", "\n", "Please discuss the following questions for around 10 minutes with your group:\n", "- Have you seen what models before?\n", "- Have you ever done one?\n", "- Why are what models useful?\n", "- When are they possible? Does your field have what 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_what_models_Discussion\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "---\n", "# Summary\n", "\n", "*Estimated timing of tutorial: 50 minutes*\n", "\n", "In this tutorial, we loaded some neural data and poked at it to understand how the dataset is organized. Then we made some basic plots to visualize (1) the average level of activity across the population and (2) the distribution of ISIs for an individual neuron. In the very last bit, we started to think about using mathematical formalisms to understand or explain some physiological phenomenon. All of this only allowed us to understand \"What\" the data looks like.\n", "\n", "This is the first step towards developing models that can tell us something about the brain. That's what we'll focus on in the next two tutorials." ] } ], "metadata": { "colab": { "collapsed_sections": [], "include_colab_link": true, "name": "W1D1_Tutorial1", "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 }