{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spike Train Analysis: Exercise 2 - Unitary Event Analysis\n",
    "\n",
    "This notebook covers the following topics:\n",
    "\n",
    "* Unitary Event Analysis"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import quantities as pq\n",
    "import elephant.unitary_event_analysis as ue\n",
    "from viziphant.unitary_event_analysis import plot_ue\n",
    "\n",
    "import utils\n",
    "\n",
    "%matplotlib inline"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Task 2: Unitary Event Analysis\n",
    "The Unitary Event (UE) analysis is a statistical method to detect excess spike synchrony beyond chance coincidences. Here we use this method to investigate the correlation between given spike trains as a function of time. To understand how the method works, you first generate by yourself spike trains with desired correlation structures, such as which neurons are correlated, how the correlations evolve in time, and so on, and apply the method to those spike trains to see how the embedded correlation is reflected in the analysis result. Then you proceed to apply the method to real experimental data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.1 Dynamics of Synchronous Activity\n",
    "\n",
    "### Generating Data\n",
    "The function `utils.generate_spike_trains_with_coinc()` generates spike trains with the following features.\n",
    "\n",
    "- The spike trains for `num_trials` trials of `trial_duration` ms duration are generated.\n",
    "- The firing rate of each neuron is a stationary Poisson process with rate `rate_b` Hz.\n",
    "- The neurons are correlated with the coincidence rate `rate_c` Hz during `coinc_duration` ms in the middle of every trial and the rest of time they are independent.\n",
    "\n",
    "### Applying UE\n",
    "Once the data are generated, we are ready to apply the UE method on them! Choose values for the following parameters of the method:\n",
    "\n",
    "- `binsize`: size of the bin for discretizing each spike train\n",
    "- `winsize`: size of the window of analysis\n",
    "- `winstep`: size of the step for moving the window\n",
    "- `pattern`: binary list indicating which neurons participate in the synchrony pattern of interest, e.g., a list `[1, 1]` indicates a pattern of 2 neurons firing together, `[1, 0, 1]` a pattern of 3 neurons where neurons 1 and 3 fire together while neuron 2 doesn't.\n",
    "  - `pattern_hash`: use `ue.hash_from_pattern(pattern)` to get a unique hash number for a specific pattern\n",
    "- `significance_level`: significance level for evaluating p-values (or correspondingly surprise), e.g. 0.05 (5%)\n",
    "\n",
    "### 2.1.1 Analysis Window Width\n",
    "The time scales of the change in firing rate and the modulation of the rate of coincidence events can differ.\n",
    "The UE uses a sliding-window approach to capture such modulation of synchrony.\n",
    "The proper choice of the sliding window width required for detecting significant excess synchrony depends on two factors:\n",
    "the width of the time interval `coinc_duration` containing excess synchrony (“hot region”), and the rate of coincidences `rate_c` relative to the independent background rate `rate_b`.\n",
    "\n",
    "#### Exercise\n",
    "* Generate spike data by choosing the background rate `rate_b`, the coincidence rate `rate_c`, size of the hot region `coinc_duration` and number of trials `num_trials`\n",
    "* Apply the UE method on the generated data with a very large window size (e.g `winsize` = `trial_duration`). Do we observe any UEs?\n",
    "* Choose a smaler window size, test several different choices: `winsize = coinc_duration`, `winsize > coinc_duration` and `winsize < coinc_duration`. How does the surprise behave in each case around the hot region? How changing the number of trials `num_trials` changes the result?\n",
    "* Apply the UE method on data with various background rates `rate_b`, while keeping the coincidence rate `rate_c` fixed. Do we abserve UEs for high `rate_b`? Why?\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "scrolled": false
   },
   "source": [
    "# Spike train parameters\n",
    "num_trials = 40\n",
    "trial_duration = 1000 * pq.ms\n",
    "rate_b = 40 * pq.Hz\n",
    "rate_c = 2 * pq.Hz\n",
    "coinc_duration = 100 * pq.ms\n",
    "\n",
    "print('Generating data...')\n",
    "spiketrains = utils.generate_spike_trains_with_coinc(rate_b, rate_c, trial_duration, coinc_duration, num_trials)\n",
    "print('...done.')"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "# Analysis parameters\n",
    "binsize = 1 * pq.ms\n",
    "winsize = 500 * pq.ms\n",
    "winstep = 2 * pq.ms\n",
    "pattern = [1, 1]\n",
    "significance_level = 0.01\n",
    "\n",
    "print('Applying the UE method...')\n",
    "pattern_hash = [ue.hash_from_pattern(pattern)]\n",
    "ue_result = ue.jointJ_window_analysis(spiketrains, binsize, winsize, winstep, pattern_hash)\n",
    "print('...done.')"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "# Plot the result\n",
    "coinc_start = (trial_duration - coinc_duration).rescale('ms') / 2.\n",
    "coinc_stop = (coinc_start + coinc_duration).rescale('ms')\n",
    "plot_ue(spiketrains, ue_result, significance_level, unit_real_ids=[1, 2], events={'coinc\\nstart': [coinc_start], 'coinc\\nstop': [coinc_stop]})\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1.2 Pairwise or higher-order synchrony\n",
    "\n",
    "#### Exercise\n",
    "* Generate spike trains of three neurons as we have done in the previous excersice. Choose a subset of neurons (e.g. [1,2] neuron 1 and 2, or [1,2,3] neuron 1, 2 and 3) to inject coincidences with the `rate_c`.\n",
    "* Investigate the significance of triplet (pattern = [1,1,1]) or pair (pattern [1,1,0]) in each case. \n",
    "* Check if the triplet (pattern = [1,1,1]) is significant for the subset [[1,2],[2,3],[1,3]]. Does the fact that each neuron is correlated with the two other neurons mean that all the neurons are correlated? Discuss what we can conclude from the results."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "# Spike train parameters\n",
    "num_trials = 100\n",
    "trial_duration = 1000 * pq.ms\n",
    "rate_b = 20 * pq.Hz\n",
    "rate_c = 1 * pq.Hz\n",
    "coinc_duration = 100 * pq.ms\n",
    "num_sync_units = 3\n",
    "subsets = [[1, 2, 3],]  # or, e.g., [[1,2], [2,3], [1,3]]\n",
    "\n",
    "# Analysis parameters\n",
    "winsize = 100 * pq.ms\n",
    "binsize = 1 * pq.ms\n",
    "winstep = 2 * pq.ms\n",
    "pattern = [1, 1, 1]\n",
    "significance_level = 0.01\n",
    "\n",
    "#Generate correlated data\n",
    "print('Generating data...')\n",
    "spiketrains = utils.generate_spike_trains_with_coinc(rate_b, rate_c, trial_duration, coinc_duration, num_trials, num_sync_units, unit_ids_sync=subsets)\n",
    "print('...done.')\n",
    "\n",
    "print('Applying the UE method...')\n",
    "pattern_hash = [ue.hash_from_pattern(pattern)]\n",
    "ue_result = ue.jointJ_window_analysis(spiketrains, binsize, winsize, winstep, pattern_hash)\n",
    "print('...done.')"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "# Plot the result\n",
    "t_winpos = ue._winpos(0*pq.ms, trial_duration, winsize, winstep)\n",
    "times = (t_winpos + winsize/2.).rescale('ms').magnitude\n",
    "\n",
    "rate_avg = ue_result['rate_avg'][:, 0].rescale('Hz')\n",
    "\n",
    "coinc_start = (trial_duration - coinc_duration).rescale('ms').magnitude / 2.\n",
    "coinc_stop = coinc_start + coinc_duration.rescale('ms').magnitude\n",
    "rate_coinc = np.zeros_like(times)\n",
    "rate_coinc[(coinc_start <= times) & (times < coinc_stop)] = rate_c\n",
    "\n",
    "fig = plt.figure(figsize=(8, 6))\n",
    "\n",
    "ax1 = plt.subplot(2, 1, 1)\n",
    "ax1.set_ylabel('Rate (Hz)')\n",
    "ax1.plot(times, rate_avg[:, 0], label='Unit 1')\n",
    "ax1.plot(times, rate_avg[:, 1], label='Unit 2')\n",
    "ax1.plot(times, rate_avg[:, 2], label='Unit 3')\n",
    "ax1.plot(times, rate_coinc, label='Coincidence')\n",
    "ax1.axhline(rate_b, lw=0.5, color='black')\n",
    "ax1.legend()\n",
    "\n",
    "ax2 = plt.subplot(2, 1, 2, sharex=ax1)\n",
    "ax2.set_xlabel('Time (ms)')\n",
    "ax2.set_ylabel('Surprise')\n",
    "ax2.plot(times, ue_result['Js'][:, 0], color='k')\n",
    "ax2.axhline(0, ls='-', lw=0.5, color='black')\n",
    "ax2.axhline(ue.jointJ(significance_level), ls='-', lw=0.5, color='red')\n",
    "ax2.axhline(ue.jointJ(1 - significance_level), ls='-', lw=0.5, color='green')\n",
    "\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1.3 Non-stationary Firing Rate\n",
    "\n",
    "#### Generating Data\n",
    "The function `utils.generate_spike_trains_with_osc_coinc()` generates spike trains with the following features.\n",
    "\n",
    "- The firing rate of each individual spike train is modulated as a sinusoidal function, where `freq_bg` is the ordinary frequency (i.e. the number of oscillations that occur at the given time), `amp_bg` is the amplitude and `offset_bg` is a non-zero center amplitude. (As an example, if we choose `freq_bg` or `amp_bg` equal to zero, we will have a stationary firing rate equal to `offset_bg`). After specifying these parameters, `N` spike trains of length `T` will be generated independently with the same rate profile.\n",
    "- An extra sinusoidal rate profile with parameters `freq_coinc`, `amp_coinc` and `offset_coinc` generates a spike train which will be copied in all N spike trains to correlate them.\n",
    "- Spike trains for `nTrials` trials are generated.\n",
    "- `RateJitter` is the parameter to perturb the background firing rate of individual spike trains in each trial to introduce non-stationarity across trials.\n",
    "\n",
    "#### Exercise\n",
    "* Generate stationary data without firing rate modulation (i.e. `freq_bg=0*pq.Hz`, `freq_coinc=0*pq.Hz`) and vary the analysis parameters (`binsize`, `winsize`, etc)\n",
    "* Play with the values of `freq_bg` and `freq_coinc` and discuss the relation between modulation in the background firing rate and the coincidence rate with the expected and empirical coincidences.\n",
    "* Discuss the relation between surprise, firing rate modulation and coincidence modulation. Is the surprise correlated with firing rate? And with coincidences?\n",
    "* Discuss the relation between the size of window of analysis and the modulation of the firing rate. In the presence of modulation of firing rate what is the best choice of the window size?"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "scrolled": false
   },
   "source": [
    "# Spike train parameters\n",
    "num_trials = 100  # number of trials\n",
    "trial_duration = 1000 * pq.ms  # trial duration\n",
    "num_sync_units = 2  # number of neurons\n",
    "\n",
    "# - background rate\n",
    "freq_bg = 5 * pq.Hz  # frequency of oscillatory rate modulation\n",
    "amp_bg = 4. * pq.Hz  # modulation depth\n",
    "offset_bg = 30 * pq.Hz  # constant rate offset\n",
    "RateJitter = 0 * pq.Hz  # jitter of background firng rate across trials\n",
    "\n",
    "# - coincidence rate\n",
    "freq_coinc = 2 * pq.Hz\n",
    "amp_coinc = 1 * pq.Hz\n",
    "offset_coinc = 0 * pq.Hz\n",
    "\n",
    "# Analysis parameters\n",
    "winsize = 100 * pq.ms\n",
    "binsize = 1 * pq.ms\n",
    "winstep = 5 * pq.ms\n",
    "pattern = [1, 1]\n",
    "significance_level = 0.01\n",
    "\n",
    "print('Generating data ...')\n",
    "data = utils.generate_spike_trains_with_osc_coinc(num_trials, num_sync_units, trial_duration, freq_coinc, amp_coinc, offset_coinc, \n",
    "                                     freq_bg, amp_bg, offset_bg, RateJitter=RateJitter)\n",
    "spiketrains = data['st']\n",
    "print('..done')\n",
    "\n",
    "print('Applying the UE method...')\n",
    "pattern_hash = [ue.hash_from_pattern(pattern)]\n",
    "ue_result = ue.jointJ_window_analysis(spiketrains, binsize, winsize, winstep, pattern_hash)\n",
    "print('..done')"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "# Plot the result\n",
    "plot_ue(spiketrains, ue_result, significance_level, unit_real_ids=[1, 2])\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.2 The role of synchrony: behavior related occurrence of unitary events\n",
    "Here we apply the UE method on the experimental data recorded from motor cortex of a macaque monkey. The details of the experiment are explained in Riehle et al. (1997), but briefly on the task description:\n",
    "  - The monkey was involved in a delayed pointing task.\n",
    "  - The duration of the delay (from the preparatory signal (PS) to the reaction signal (RS)) was selected randomly from four possible durations (600, 900, 1200, 1500 ms) from trial to trial.\n",
    "  - 36 trials with the longest delay (1500 ms) were pooled in this example. Thus the monkey could expect the RS to occur at three successive moments (ES1, ES2, ES3) before it actually occurred at 1500 ms.\n",
    "\n",
    "Load `./data/Data14.npy` and `./data/Data15.npy`, which are the real spike train data of two simultaneously recorded neurons in the experiment of Riehle et al. (1997). Apply the UE analysis to them.\n",
    "\n",
    "* Is there a relation between UEs and behaviour?\n",
    "* Can the UEs be simply a reflection of rate modulation?"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "scrolled": false
   },
   "source": [
    "# Analysis parameters\n",
    "num_sync_units = 2\n",
    "data_id1, data_id2 = 14, 15\n",
    "winsize = 100 * pq.ms\n",
    "binsize = 5 * pq.ms\n",
    "winstep = 10 * pq.ms\n",
    "pattern = [1, 1]\n",
    "significance_level = 0.05\n",
    "method = 'analytic_TrialAverage'\n",
    "\n",
    "# Load the data\n",
    "data1 = np.load(f'./data/Data{data_id1}.npy', allow_pickle=True, encoding='latin1')\n",
    "data2 = np.load(f'./data/Data{data_id2}.npy', allow_pickle=True, encoding='latin1')\n",
    "spiketrains = [[x, y] for x, y in zip(data1.item()['st'], data2.item()['st'])]\n",
    "\n",
    "print('Applying the UE analysis...')\n",
    "pattern_hash = [ue.hash_from_pattern(pattern)]\n",
    "ue_result = ue.jointJ_window_analysis(spiketrains, binsize, winsize, winstep, pattern_hash, method=method)\n",
    "print('...done')"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "plot_ue(spiketrains, ue_result, significance_level, unit_real_ids=(data_id1, data_id2), events={'PS':[300*pq.ms],'ES1':[900*pq.ms], 'ES2':[1200*pq.ms], 'ES3':[1500*pq.ms],'RS': [1800*pq.ms]})\n",
    "plt.show()"
   ],
   "outputs": [],
   "execution_count": null
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}