{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# The Unitary Events Analysis" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The executed version of this tutorial is at https://elephant.readthedocs.io/en/latest/tutorials/unitary_event_analysis.html\n", "\n", "The Unitary Events (UE) analysis \\[1\\] tool allows us to reliably detect correlated spiking activity that is not explained by the firing rates of the neurons alone. It was designed to detect coordinated spiking activity that occurs significantly more often than predicted by the firing rates of the neurons. The method allows one to analyze correlations not only between pairs of neurons but also between multiple neurons, by considering the various spike patterns across the neurons. In addition, the method allows one to extract the dynamics of correlation between the neurons by perform-ing the analysis in a time-resolved manner. This enables us to relate the occurrence of spike synchrony to behavior.\n", "\n", "The algorithm:\n", "\n", "1. Align trials, decide on width of analysis window.\n", "2. Decide on allowed coincidence width.\n", "3. Perform a sliding window analysis. In each window:\n", " 1. Detect and count coincidences.\n", " 2. Calculate expected number of coincidences.\n", " 3. Evaluate significance of detected coincidences.\n", " 4. If significant, the window contains Unitary Events.\n", "4. Explore behavioral relevance of UE epochs.\n", "\n", "References:\n", "\n", "1. Grün, S., Diesmann, M., Grammont, F., Riehle, A., & Aertsen, A. (1999). Detecting unitary events without discretization of time. Journal of neuroscience methods, 94(1), 67-79." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import random\n", "import string\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import quantities as pq\n", "import neo\n", "\n", "import elephant.unitary_event_analysis as ue\n", "from elephant.datasets import download_datasets\n", "\n", "# Fix random seed to guarantee fixed output\n", "random.seed(1224)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we download a data file containing spike train data from multiple trials of two neurons." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Download data\n", "repo_path='tutorials/tutorial_unitary_event_analysis/data/dataset-1.nix'\n", "filepath=download_datasets(repo_path)" ] }, { "cell_type": "markdown", "metadata": { "nbsphinx": "hidden" }, "source": [ "# Write a plotting function" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "nbsphinx": "hidden" }, "outputs": [], "source": [ "# borrowed from Viziphant\n", "\n", "plot_params_default = {\n", " # epochs to be marked on the time axis\n", " 'events': {},\n", " # figure size\n", " 'figsize': (10, 12),\n", " # right margin\n", " 'right': 0.9,\n", " # top margin\n", " 'top': 0.9,\n", " # bottom margin\n", " 'bottom': 0.1,\n", " # left margin\n", " 'left': 0.1,\n", " # horizontal white space between subplots\n", " 'hspace': 0.5,\n", " # width white space between subplots\n", " 'wspace': 0.5,\n", " # font size\n", " 'fsize': 12,\n", " # the actual unit ids from the experimental recording\n", " 'unit_real_ids': None,\n", " # line width\n", " 'lw': 2,\n", " # marker size for the UEs and coincidences\n", " 'ms': 5,\n", " # figure title\n", " 'suptitle': None,\n", "}\n", "\n", "\n", "def plot_ue(spiketrains, Js_dict, significance_level=0.05,\n", " **plot_params):\n", " \"\"\"\n", " Plots the results of pairwise unitary event analysis as a column of six\n", " subplots, comprised of raster plot, peri-stimulus time histogram,\n", " coincident event plot, coincidence rate plot, significance plot and\n", " unitary event plot, respectively.\n", "\n", " Parameters\n", " ----------\n", " spiketrains : list of list of neo.SpikeTrain\n", " A nested list of trials, neurons and their neo.SpikeTrain objects,\n", " respectively. This should be identical to the one used to generate\n", " Js_dict.\n", " Js_dict : dict\n", " The output of\n", " :func:`elephant.unitary_event_analysis.jointJ_window_analysis`\n", " function. The values of each key has the shape of:\n", "\n", " * different window --> 0-axis.\n", " * different pattern hash --> 1-axis;\n", "\n", " Dictionary keys:\n", "\n", " 'Js': list of float\n", " JointSurprise of different given patterns within each window.\n", " 'indices': list of list of int\n", " A list of indices of pattern within each window.\n", " 'n_emp': list of int\n", " The empirical number of each observed pattern.\n", " 'n_exp': list of float\n", " The expected number of each pattern.\n", " 'rate_avg': list of float\n", " The average firing rate of each neuron.\n", "\n", " significance_level : float\n", " The significance threshold used to determine which coincident events\n", " are classified as unitary events within a window.\n", " **plot_params\n", " User-defined plotting parameters used to update the default plotting\n", " parameter values. The valid keys:\n", "\n", " 'events' : list\n", " Epochs to be marked on the time axis.\n", " 'figsize' : tuple of int\n", " The dimensions for the figure size.\n", " 'right' : float\n", " The size of the right margin.\n", " 'top' : float\n", " The size of the top margin.\n", " 'bottom' : float\n", " The size of the bottom margin.\n", " 'left' : float\n", " The size of the left margin.\n", " 'hspace' : flaot\n", " The size of the horizontal white space between subplots.\n", " 'wspace' : float\n", " The width of the white space between subplots.\n", " 'fsize' : int\n", " The size of the font.\n", " 'unit_real_ids' : list of int\n", " The unit ids from the experimental recording.\n", " 'lw' : int\n", " The default line width.\n", " 'ms' : int\n", " The marker size for the unitary events and coincidences.\n", "\n", " Returns\n", " -------\n", " result : FigureUE\n", " The container for Axes objects generated by the function. Individual\n", " axes can be accessed using the following identifiers:\n", "\n", " * axes_spike_events : matplotlib.axes.Axes\n", "\n", " Contains the elements of the spike events subplot.\n", " * axes_spike_rates : matplotlib.axes.Axes\n", "\n", " Contains the elements of the spike rates subplot.\n", " * axes_coincident_events : matplotlib.axes.Axes\n", "\n", " Contains the elements of the coincident events subplot.\n", " * axes_coincidence_rates : matplotlib.axes.Axes\n", "\n", " Contains the elements of the coincidence rates subplot.\n", " * axes_significance : matplotlib.axes.Axes\n", "\n", " Contains the elements of the statistical significance subplot.\n", " * axes_unitary_events : matplotlib.axes.Axes\n", "\n", " Contains the elements of the unitary events subplot.\n", "\n", " Examples\n", " --------\n", " Unitary Events of homogenous Poisson random processes.\n", "\n", " Since we don't expect to find significant correlations in random processes,\n", " we show non-significant events (``significance_level=0.34``). Typically,\n", " in your analyses, the significant level threshold is ~0.05.\n", "\n", " .. plot::\n", " :include-source:\n", "\n", " import matplotlib.pyplot as plt\n", " import numpy as np\n", " import quantities as pq\n", "\n", " import viziphant\n", " from elephant.spike_train_generation import homogeneous_poisson_process\n", " from elephant.unitary_event_analysis import jointJ_window_analysis\n", "\n", " np.random.seed(10)\n", "\n", " spiketrains1 = [homogeneous_poisson_process(rate=20 * pq.Hz,\n", " t_stop=2 * pq.s) for _ in range(5)]\n", " spiketrains2 = [homogeneous_poisson_process(rate=50 * pq.Hz,\n", " t_stop=2 * pq.s) for _ in range(5)]\n", "\n", " spiketrains = np.stack((spiketrains1, spiketrains2), axis=1)\n", " ue_dict = jointJ_window_analysis(spiketrains,\n", " bin_size=5 * pq.ms,\n", " win_size=100 * pq.ms,\n", " win_step=10 * pq.ms)\n", " viziphant.unitary_event_analysis.plot_ue(spiketrains, Js_dict=ue_dict,\n", " significance_level=0.34,\n", " unit_real_ids=['1', '2'])\n", " plt.show()\n", "\n", " Refer to `UEA Tutorial `_ for real-case scenario.\n", " \"\"\"\n", " n_trials = len(spiketrains)\n", " n_neurons = len(spiketrains[0])\n", "\n", " input_parameters = Js_dict['input_parameters']\n", " t_start = input_parameters['t_start']\n", " t_stop = input_parameters['t_stop']\n", " bin_size = input_parameters['bin_size']\n", " win_size = input_parameters['win_size']\n", " win_step = input_parameters['win_step']\n", " pattern_hash = input_parameters['pattern_hash']\n", " if len(pattern_hash) > 1:\n", " raise ValueError(f\"To not clutter the plots, only one pattern hash is \"\n", " f\"required; got {pattern_hash}. You can call this \"\n", " f\"function multiple times for each hash at a time.\")\n", " for key in ['Js', 'n_emp', 'n_exp', 'rate_avg']:\n", " Js_dict[key] = Js_dict[key].squeeze()\n", " neurons_participated = ue.inverse_hash_from_pattern(pattern_hash,\n", " N=n_neurons).squeeze()\n", "\n", " t_winpos = ue._winpos(t_start=t_start, t_stop=t_stop, win_size=win_size,\n", " win_step=win_step)\n", " Js_sig = ue.jointJ(significance_level)\n", "\n", " # figure format\n", " plot_params_user = plot_params\n", " plot_params = plot_params_default.copy()\n", " plot_params.update(plot_params_user)\n", " if plot_params['unit_real_ids'] is None:\n", " plot_params['unit_real_ids'] = ['not specified'] * n_neurons\n", " if len(plot_params['unit_real_ids']) != n_neurons:\n", " raise ValueError('length of unit_ids should be' +\n", " 'equal to number of neurons!')\n", " plt.rcParams.update({'font.size': plot_params['fsize']})\n", " ls = '-'\n", " alpha = 0.5\n", "\n", " fig, axes = plt.subplots(nrows=6, sharex=True,\n", " figsize=plot_params['figsize'])\n", " axes[5].sharey(axes[0])\n", " axes[0].sharey(axes[2])\n", "\n", " for ax in (axes[0], axes[2], axes[5]):\n", " for n in range(n_neurons):\n", " for tr, data_tr in enumerate(spiketrains):\n", " ax.plot(data_tr[n].rescale('ms').magnitude,\n", " np.full_like(data_tr[n].magnitude,\n", " fill_value=n * n_trials + tr),\n", " '.', markersize=0.5, color='k')\n", " for n in range(1, n_neurons):\n", " # subtract 0.5 to separate the raster plots;\n", " # otherwise, the line crosses the raster spikes\n", " ax.axhline(n * n_trials - 0.5, lw=0.5, color='k')\n", " ymax = max(ax.get_ylim()[1], 2 * n_trials - 0.5)\n", " ax.set_ylim([-0.5, ymax])\n", " ax.set_yticks([n_trials - 0.5, 2 * n_trials - 0.5])\n", " ax.set_yticklabels([1, n_trials], fontsize=plot_params['fsize'])\n", " ax.set_ylabel('Trial', fontsize=plot_params['fsize'])\n", "\n", " for i, ax in enumerate(axes):\n", " ax.set_xlim([t_winpos[0], t_winpos[-1] + win_size])\n", " ax.text(-0.05, 1.1, string.ascii_uppercase[i],\n", " transform=ax.transAxes, size=plot_params['fsize'] + 5,\n", " weight='bold')\n", " for key in plot_params['events'].keys():\n", " for event_time in plot_params['events'][key]:\n", " ax.axvline(event_time, ls=ls, color='r', lw=plot_params['lw'],\n", " alpha=alpha)\n", "\n", " axes[0].set_title('Spike Events')\n", " axes[0].text(1.0, 1.0, f\"Unit {plot_params['unit_real_ids'][-1]}\",\n", " fontsize=plot_params['fsize'] // 2,\n", " horizontalalignment='right',\n", " verticalalignment='bottom',\n", " transform=axes[0].transAxes)\n", " axes[0].text(1.0, 0, f\"Unit {plot_params['unit_real_ids'][0]}\",\n", " fontsize=plot_params['fsize'] // 2,\n", " horizontalalignment='right',\n", " verticalalignment='top',\n", " transform=axes[0].transAxes)\n", "\n", " axes[1].set_title('Spike Rates')\n", " for n in range(n_neurons):\n", " axes[1].plot(t_winpos + win_size / 2.,\n", " Js_dict['rate_avg'][:, n].rescale('Hz'),\n", " label=f\"Unit {plot_params['unit_real_ids'][n]}\",\n", " lw=plot_params['lw'])\n", " axes[1].set_ylabel('Hz', fontsize=plot_params['fsize'])\n", " axes[1].legend(fontsize=plot_params['fsize'] // 2, loc='upper right')\n", " axes[1].locator_params(axis='y', tight=True, nbins=3)\n", "\n", " axes[2].set_title('Coincident Events')\n", " for n in range(n_neurons):\n", " if not neurons_participated[n]:\n", " continue\n", " for tr, data_tr in enumerate(spiketrains):\n", " indices = np.unique(Js_dict['indices'][f'trial{tr}'])\n", " axes[2].plot(indices * bin_size,\n", " np.full_like(indices, fill_value=n * n_trials + tr),\n", " ls='', ms=plot_params['ms'], marker='s',\n", " markerfacecolor='none',\n", " markeredgecolor='c')\n", " axes[2].set_ylabel('Trial', fontsize=plot_params['fsize'])\n", "\n", " axes[3].set_title('Coincidence Rates')\n", " axes[3].plot(t_winpos + win_size / 2.,\n", " Js_dict['n_emp'] / (\n", " win_size.rescale('s').magnitude * n_trials),\n", " label='Empirical', lw=plot_params['lw'], color='c')\n", " axes[3].plot(t_winpos + win_size / 2.,\n", " Js_dict['n_exp'] / (\n", " win_size.rescale('s').magnitude * n_trials),\n", " label='Expected', lw=plot_params['lw'], color='m')\n", " axes[3].set_ylabel('Hz', fontsize=plot_params['fsize'])\n", " axes[3].legend(fontsize=plot_params['fsize'] // 2, loc='upper right')\n", " axes[3].locator_params(axis='y', tight=True, nbins=3)\n", "\n", " axes[4].set_title('Statistical Significance')\n", " axes[4].plot(t_winpos + win_size / 2., Js_dict['Js'], lw=plot_params['lw'],\n", " color='k')\n", " axes[4].axhline(Js_sig, ls='-', color='r')\n", " axes[4].axhline(-Js_sig, ls='-', color='g')\n", " xlim_ax4 = axes[4].get_xlim()[1]\n", " alpha_pos_text = axes[4].text(xlim_ax4, Js_sig, r'$\\alpha +$', color='r',\n", " horizontalalignment='right',\n", " verticalalignment='bottom')\n", " alpha_neg_text = axes[4].text(xlim_ax4, -Js_sig, r'$\\alpha -$', color='g',\n", " horizontalalignment='right',\n", " verticalalignment='top')\n", " axes[4].set_yticks([ue.jointJ(1 - significance_level), ue.jointJ(0.5),\n", " ue.jointJ(significance_level)])\n", " # Try '1 - 0.34' to see the floating point errors\n", " axes[4].set_yticklabels(np.round([1 - significance_level, 0.5,\n", " significance_level], decimals=6))\n", "\n", " # autoscale fix to mind the text positions.\n", " # See https://stackoverflow.com/questions/11545062/\n", " # matplotlib-autoscale-axes-to-include-annotations\n", " plt.get_current_fig_manager().canvas.draw()\n", " for text_handle in (alpha_pos_text, alpha_neg_text):\n", " bbox = text_handle.get_window_extent()\n", " bbox_data = bbox.transformed(axes[4].transData.inverted())\n", " axes[4].update_datalim(bbox_data.corners(), updatex=False)\n", " axes[4].autoscale_view()\n", "\n", " mask_nonnan = ~np.isnan(Js_dict['Js'])\n", " significant_win_idx = np.nonzero(Js_dict['Js'][mask_nonnan] >= Js_sig)[0]\n", " t_winpos_significant = t_winpos[mask_nonnan][significant_win_idx]\n", " axes[5].set_title('Unitary Events')\n", " if len(t_winpos_significant) > 0:\n", " for n in range(n_neurons):\n", " if not neurons_participated[n]:\n", " continue\n", " for tr, data_tr in enumerate(spiketrains):\n", " indices = np.unique(Js_dict['indices'][f'trial{tr}'])\n", " indices_significant = []\n", " for t_sig in t_winpos_significant:\n", " mask = (indices * bin_size >= t_sig\n", " ) & (indices * bin_size < t_sig + win_size)\n", " indices_significant.append(indices[mask])\n", " indices_significant = np.hstack(indices_significant)\n", " indices_significant = np.unique(indices_significant)\n", " # does nothing if indices_significant is empty\n", " axes[5].plot(indices_significant * bin_size,\n", " np.full_like(indices_significant,\n", " fill_value=n * n_trials + tr),\n", " ms=plot_params['ms'], marker='s', ls='',\n", " mfc='none', mec='r')\n", " axes[5].set_xlabel(f'Time ({t_winpos.dimensionality})',\n", " fontsize=plot_params['fsize'])\n", " for key in plot_params['events'].keys():\n", " for event_time in plot_params['events'][key]:\n", " axes[5].text(event_time - 10 * pq.ms,\n", " axes[5].get_ylim()[0] - 35, key,\n", " fontsize=plot_params['fsize'], color='r')\n", "\n", " plt.suptitle(plot_params['suptitle'], fontsize=20)\n", " plt.subplots_adjust(top=plot_params['top'],\n", " right=plot_params['right'],\n", " left=plot_params['left'],\n", " bottom=plot_params['bottom'],\n", " hspace=plot_params['hspace'],\n", " wspace=plot_params['wspace'])\n", " \n", " return axes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Load data and extract spiketrains" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "io = neo.io.NixIO(f\"{filepath}\",'ro')\n", "block = io.read_block()\n", "\n", "spiketrains = []\n", "# each segment contains a single trial\n", "for ind in range(len(block.segments)):\n", " spiketrains.append (block.segments[ind].spiketrains)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Calculate Unitary Events" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "UE = ue.jointJ_window_analysis(\n", " spiketrains, bin_size=5*pq.ms, win_size=100*pq.ms, win_step=10*pq.ms, pattern_hash=[3])\n", "\n", "plot_ue(spiketrains, UE, significance_level=0.05)\n", "plt.show()" ] } ], "metadata": { "interpreter": { "hash": "623e048a0474aa032839f97d38ba0837cc9041adc49a14b480c72f2df8ea99e3" }, "kernelspec": { "display_name": "inm-elephant", "language": "python", "name": "inm-elephant" }, "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.8.10" }, "latex_envs": { "LaTeX_envs_menu_present": true, "autocomplete": true, "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 1, "hotkeys": { "equation": "Ctrl-E", "itemize": "Ctrl-I" }, "labels_anchors": false, "latex_user_defs": false, "report_style_numbering": false, "user_envs_cfg": false }, "toc": { "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }