Background.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ad6a89a2",
   "metadata": {},
   "source": [
    "# IIA project GG3: Neural Data Analysis\n",
    "\n",
    "Easter 2023<br>\n",
    "Project Leader: Yashar Ahmadian (ya311)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4e2031a-c438-4187-90bc-87cfc4ef7f4f",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Background"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43f7b50a-2dcb-4702-b1ab-38d7e58652a2",
   "metadata": {},
   "source": [
    "**Note:** The main point of this background handout is to introduce some terminology (which appear all in boldface),<br>\n",
    "and mathematical notation that will be used in the project notebook, beginning with the section,<br>\n",
    "\"What is the right model of LIP?\".<br> \n",
    "Deep understanding of this Background  section is not required for carrying out the project. But, apart from the<br>\n",
    "above reason, you are encourged to read it to understand the scientific motivations and significance of this problem."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aea7ff93-6ddb-4850-9508-13a815c36021",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "#### One-paragraph primer on neurons and spikes\n",
    "\n",
    "Neurons communicate and respond to stimuli by emitting action potentials or \"**spikes**\". You can <br>\n",
    "think of a spike as a sharp impulse of electric voltage that the neuron sends to its target neurons.<br>\n",
    "A sequence of spikes is called a **spike train**. We can represent a spike train, e.g. by providing the time<br>\n",
    "of different spikes. The time interval between two consecutive spikes is referred to as an inter-spike interval (ISI). \n",
    "\n",
    "Alternatively, we will represent a spike train using a binary (0/1) sequence or time-series, with each sequence<br>\n",
    "element corresponding to a narrow (e.g. 10 ms wide) **time bin**, with value 1 if a spike had occured in that time bin<br>\n",
    "and 0 otherwise. We can thus denote the spike train by a sequence $(n_t)_{t=1}^T$, of length $T$, where $n_t$ is the spike count in <br>\n",
    "in the $t$-th time bin (aka time-step). As we said, if the time bin is small enough the count $n_t$ can be either<br>\n",
    " 0 or 1. But  we will allow $n_t$ to acquire any (non-negative) integer value -- this is a real possibility<br>\n",
    "when the time bin is not so small.\n",
    "\n",
    "\n",
    "An important quantity characterising a neuron's response is its **firing rate**, namely the frequency or speed with which<br>\n",
    "the neuron fires spikes. It is quantified as spikes per unit of time, and measured in, *e.g.*, units of Hz.<br>\n",
    "The firing rate can and does change over time, so can be described by a function of time $r(t)$ or $r_t$.\n",
    "\n",
    "Experiments are normally carried out in **trials** where the animal, from whose brain we are recording neural activity,<br>\n",
    "performs a certain task over and over. The **neural data** we will be analysing in this project takes the form <br>\n",
    "of a collection of spike trains (form a simulated neuron) recorded over many trials. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f24f4770-e482-4324-99a2-cbc4d5890595",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "#### Random dot motion task\n",
    "The so-called random dot motion task is a standard psychophysical task commonly used in studies of animal and human *perceptual decision making*.<br> \n",
    "In each trial, the animal (e.g. a monkey) has to first fixate its eyes at the center of the monitor screen.<br>\n",
    "After fixation random dot motion (**RDM**) stimuli appear at a fixed location on the screen. <br>\n",
    "Based on the direction of motion of the random dots, in order to get reward, the monkey has to make a binary decision (within a limited time window)<br>\n",
    "that s/he announces by saccading to a left or right target.<br>\n",
    "Random dots can move with different levels of **coherence**, a parameter that controls the task difficulty, and can vary from trial to trial. <br>\n",
    "See the [first two videos here](https://elifesciences.org/articles/08825#media1) for examples of random dot motion stimui with high and low coherence, respectively.\n",
    "\n",
    "#### -----------------------   Figure 1   -----------------------\n",
    "<img src=\"figs/rdm_task.png\" width=500 height=500 />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84e9190e-1c13-4f14-9cce-68dc38cce6ae",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### Evidence accumulation by an ideal observer agent:\n",
    "\n",
    "To solve this task, the monkey has to monitor the motion of the dots during the period the stimulus is on,<br>\n",
    "and then judge the direction of motion, based on noisy (and possibly very incoherent) observations.<br> In principle,\n",
    "the longer you view the stimulus, the more evidence you have in favor one or the other alternative.<br>\n",
    "(But you also don't want to wait forever to get that reward.)\n",
    "\n",
    "This is similar to judging whether heads or tails are more probable, by observing the outcome of successive flips of a coin.<br>\n",
    "In this case the *probability of heads*, $p$,  or more precisely $|p - 0.5|$, mimicks coherence and controls the difficulty of <br>\n",
    "this decision making (or hypothesis testing) task.\n",
    "\n",
    "An ideal observer/decision-maker should accumulate evidence over all observed coin throws to make a decision. <br> \n",
    "Optimal judgment (assuming an indifferent prior between the two options) will be based on the likelihood ratio (i.e. the ratio  \n",
    "of the  likelihoods of the particular sequence of observed coin flips, under the two hypotheses),<br>\n",
    "or equivalently, based on the logarithm of this ratio. (See the notes in **Appendix A** for further details/derivation.)\n",
    "\n",
    "Log likelihood ratio (LLR) accumulates *additively* over independent coin flips, with increments that are stochastic<br>\n",
    " (since coin flip outcomes are).  However, for $p\\neq 0.5$ there will be a systematic drift; *e.g.*, for high $p$ <br>\n",
    "(high coherence in favor of the heads hypothesis) there will be a robust, systematic upward trend or \"drift\". <br>\n",
    "On the other hand, for $p = 0.5$, the LLR moves, or \"diffuses\", as a random walk. In general, the trajectory of LLR can<br>\n",
    " thus be described by a **drift-diffusion model** (see below). Finally, when LLR reaches a fixed positive or negative bound, <br>\n",
    "the decision is made -- in favor of heads or tails being more likely, respectively. A variable such as LLR which gets <br>\n",
    "updated dynamically based on observations and forms the basis of decision making is called a **decision variable**. \n",
    "\n",
    "Figure 2 below shows evidence accumulation for the coin flip task: <br> \n",
    "easier case with $p = 0.8$ (yellow), harder case with $p = 0.6$ (red).\n",
    "\n",
    "#### ----------------------------- Figure 2 -----------------------------\n",
    "<img src=\"figs/plot_coins_accumulation_4.png\" width=500 height=500 />\n",
    "\n",
    "#### --------------------------------------------------------------------\n",
    "\n",
    "This ideal observer model is extremely successful in explaining several systematic aspects of monkeys' (and humans')<br>\n",
    "decision making behaviour in this task.<br>\n",
    "In theoretical neuroscience a continuous-time version of the drift-diffusion model is used. <br>\n",
    "A sample trajectory of such a model is shown in Figure 3. As everything is discrete in computers, the ramp model simulator<br>\n",
    "provided to you at the outset of the project is an approximation to such a continuous-time model.<br> \n",
    "\n",
    "#### ----------------------------- Figure 3 -----------------------------\n",
    "<img src=\"figs/driftdiffusion_header.png\" width=500 height=500 />\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2410c7c2-c913-48e7-9ce1-9a7b4f38563a",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### Evidence accumulation in area LIP\n",
    "\n",
    "Area LIP (lateral intra-parietal cortex) of the cerebral cortex is involved in the control of saccadic eye movements.<br>\n",
    "A series of classic studies in systems neuroscience showed that many neurons in this area exhibit \"ramping activity\",<br>\n",
    "which in many respects mimicks the behavior of the decision variable (the LLR) in the ideal observer model.<br>\n",
    "For example, the firing rate ramps up faster in higher coherence trials. Moreover, the firing rate at the<br> \n",
    "onset of eye movement (saccade) is consistently at the same level (independent of the stimulus coherence),<br>  \n",
    "suggesting the decisionwas contingent on the firing rate of these neurons reaching a fixed, preset level.<br>\n",
    "These findings supported the hypothesis that area LIP and its neurons may be involved in the implementation of<br>\n",
    " some kind of evidence accumulation algorithm similar to and approximating the one used by the ideal-observer. <br>\n",
    "In particular, under this hypothesis, the firing rate of these LIP neurons is directly representing accumulated evidence. \n",
    "\n",
    "Figure 4 shows the PSTH of an LIP neuron for different coherence (and direction of motion, dashed lines) conditions.<br>\n",
    "The right side of the plot shows PSTH's after aligning trials based on the time when the marking made a saccade.<br>\n",
    "Remarkably, independent of trial coherence, the ramping firing rates have reached the same level at the decision time,<br>\n",
    "exactly as the decision variable in the drift diffusion model of the LLR in the coin flip example are supposed to behave.\n",
    "\n",
    "\n",
    "#### -----------------------   Figure 4   -----------------------\n",
    "<img src=\"figs/rdm_task_lip.png\" width=500 height=500 />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f18d919-a770-4f60-9f2f-224a9fec6c03",
   "metadata": {},
   "source": [
    "### The monkey and us\n",
    "\n",
    "As you will see, if Bayesian inference provides a good theoretical account of the Monkey's behaviour,<br>\n",
    "it is also a tool for us engineers and neuroscientists in studying the monkey's behaviour, the behaviour of <br>\n",
    "its neurons (or their simulators!). A fun fact about neural data analysis and computational neuroscience is that<br>\n",
    "we as engineers and scientists routinely using techniques and methods of inference that we hypothesise or believe<br>\n",
    "animal brains are employing in their day-to-day lives: to perceive the world, and to make good decisions -- that can<br>\n",
    "make the difference between death or survival till another day! \n",
    "\n",
    "In particular, you will be building techniques for accumulating evidence<br>\n",
    "in order to decide between two alternative proposed mechanisms underlying LIP neural activity. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27ca6c12-5c30-4cc2-bbb1-29fb05547262",
   "metadata": {},
   "source": [
    "# What is the right model of LIP?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6d0a6cd-7d16-43f4-9549-cb99481b2d01",
   "metadata": {},
   "source": [
    "### Alternative hypothesis: stepping model\n",
    "\n",
    "As we saw in the Background section, classic studies suggested that LIP neurons which exhibit ramping activity in<br>\n",
    "their trial-averaged PSTH's are involved in evidence accumulation. However, the story became more complicated,<br>\n",
    "when in 2015, [Latimer et al.](https://www.science.org/doi/10.1126/science.aaa4056) provided evidence that most LIP neurons are better modelled<br>\n",
    " as neurons with a \"stepping firing rate\". In this alternative model\n",
    " the rate does not continuously ramp up or down<br>  (albeit via a random walk) as in a drif-diffusion model.\n",
    "Rather, the rate is piece-wise constant:<br> it starts relatively low, but at some time point it jumps (\"steps\") up discontinuously to <br>\n",
    "a higher firing rate level. The jump point is random and varies from trial to trial, according to some distribution. <br>\n",
    "\n",
    "#### -------------------------------   Figure 5   -------------------------------\n",
    "<img src=\"figs/latimer-step-ramp.png\" width=600 height=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ff6abb3-eeb2-4630-8dfb-5ae866a8b81b",
   "metadata": {},
   "source": [
    "We will refer to these two competing hypotheses or models as the **ramping** and **stepping models**, respectively<br>\n",
    "(other common synonyms for the ramping model are \"the drif-diffusion model\", mentioned above, and \"the diffusion-to-bound model\";<br>\n",
    "we will also use **jump model** as synonymous with the stepping model.) \n",
    "\n",
    "In this project we aim to develop tools that allow us to reject or accept one of these hypotheses<br>\n",
    "based on observed spike trains. Understanding which of the two is a more accurate description of LIP activity <br>\n",
    "is scientifically significant. The ramping hypothesis suggests that LIP cortex is responsible<br>\n",
    "for accumulating evidence to inform and make decisions. On the other hand, the binary nature of the stepping model<br>\n",
    "suggests that LIP is downstream of the evidence accumulating area, and may simply reflect, in its activity, the decision already made<br>\n",
    "in an upstream area.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbad9546-c685-41f5-b86a-a0c2454ac6a9",
   "metadata": {},
   "source": [
    "__________________________________________________________________________________________________________________________________________________________________"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23b913af-6a16-4ba7-a9d4-8d05d749c1a9",
   "metadata": {},
   "source": [
    "# Appendix A: Bayesian decision making"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fac86e09-bc13-4fc5-a9f6-894a63379310",
   "metadata": {},
   "source": [
    "(Note that this appendix has close to no bearing on the project. It is provided mostly to\n",
    "satiate your curiousity<br> in case the stuff on log-likeihood based decision making in\n",
    "the Background section was completely new to you.) \n",
    "$\\newcommand{\\P}{\\mathbb{P}}$\n",
    "$\\newcommand{\\x}{x_{1:t}}$\n",
    "$\\newcommand{\\H}{\\mathrm{H}}$\n",
    "$\\newcommand{\\T}{\\mathrm{T}}$\n",
    "$\\newcommand{\\LR}{\\mathrm{LR}}$\n",
    "$\\newcommand{\\LLR}{\\mathrm{LLR}}$\n",
    "\n",
    "An ideal observer (such as you in this project!) should carry out its inferences based\n",
    "on the rules of probability theory.\n",
    "In other words, it has to carry out Bayesian inference. This means it has to compute \n",
    "the posterior probabilities of (discsrete) events or hypotheses (or posterior distributions\n",
    "for continuous latent variables) conditioned on observed data. In the coin flip example, the\n",
    "agent needs to decided between the following two alternative hypotheses: \n",
    "\n",
    "H: the coin is biased towards heads.<br>\n",
    "T: the coin is biased towards tails.\n",
    "\n",
    "Instead, here I will discuss the following simpler problem. Suppose the coherence $c = |p-0.5|$ \n",
    "is given to the ideal observer, and the observer's task is to (more) simply judge whether $p = 0.5 + c$ or \n",
    "$p = 0.5 - c$. I will still refer to these two alternatives by H and T, respectively.\n",
    "Let $x_t$ denote the outcome of the $t$-th coin flip. So $x_t$ can either be, say, 1 if coin lands heads,\n",
    "or 0 if it lands tails. I will show the set of all observed coin flips up to and including the $t$-th flip \n",
    "by $x_{1:t}$. By the Bayes' rule (itself a simple consequence of the product rule of probability theory)\n",
    "the posterior probability of H is given by \n",
    "\n",
    "$\n",
    "\\P(\\H | \\x) = \\frac{\\P(\\x| \\H) \\P(\\H)}{\\P(\\x)} = \\frac{\\P(\\x| \\H) \\P(\\H)}{\\P(\\x| \\H) \\P(\\H) + \\P(\\x| \\T) \\P(\\T)},\n",
    "$\n",
    "\n",
    "where $\\P(\\H)$ and $\\P(\\T)$ are the prior probabilities of H and T. If we assume the agent is *a prior* indifferent\n",
    "between these two alternatives, we have $\\P(\\H) = \\P(\\T) = 0.5$, and we can cancel them and obtain\n",
    "\n",
    "$\n",
    "\\P(\\H | \\x) =  \\frac{\\P(\\x| \\H)}{\\P(\\x| \\H) + \\P(\\x| T)} = \\frac{ \\LR_t}{1 + \\LR_t}, \\qquad\\qquad\\qquad\n",
    "$                  Eq. (1)\n",
    "\n",
    "where I defined $\\LR_t$ to be the likelihood ratio (at coin flip $t$):\n",
    "\n",
    "$\n",
    "\\LR_t \\equiv \\frac{\\P(\\x| \\H)}{\\P(\\x| \\T)}.\n",
    "$\n",
    "\n",
    "Finally, in a *decision making* task, simply calculating the posterior probability is not enough: ultimately you \n",
    "have to call the shots and decide H or T. The ideal agent has to do this when the posterior probability in favor of \n",
    "one or the other alternative exceeds some preset threshold. Thus the agent should decide H when and if the posterior probability \n",
    "$\\P(\\H | \\x)$ exceeds some preset threshold or bound, $B> 0$; and alternatively decided $T$ if $\\P(\\T | \\x) = 1 - \\P(\\H | \\x)$ exceeds\n",
    "that same bound (to be fair between the two). The latter is equivalent to $\\P(\\H | \\x)$ going below $1 - B$.\n",
    "\n",
    "Now since, according to Eq. (1), $\\P(\\H | \\x)$ depends monotonically on LR, thresholding the former is equivalent\n",
    "to thresholding the latter (using a different, but equivalent bound). Since the logarithm is also a monotonically increasing function,\n",
    "the same is true if we instead base the decision on the log likelihood ratio (LLR):\n",
    "\n",
    "$\n",
    "\\LLR_t = \\log \\LR_t = \\log\\frac{\\P(\\x| \\H)}{\\P(\\x| \\T)}.\n",
    "$\n",
    "\n",
    "In this case, we decide H as soon as $\\LLR_t  > B$ and decide T as soon as $\\LLR_t  < -B$ (the decision is made as soon as \n",
    "one or the other alternative happens). (Note that the $B$ here is different from, but equivalent to, the $B$ defined above as a bound for the \n",
    "posterior probability.)\n",
    "\n",
    "The nice thing about the $\\LLR_t$ is that (assuming different coin flips are independent events) it is *additive*. That's because\n",
    "probabilities of independent events multiply and so their logs add. In particular, since $\\P(x_{1:t}| \\H) = \\P(x_t| \\H)\\P(x_{1:t-1}| \\H)$, we have the update rule:\n",
    "\n",
    "$\n",
    "\\LLR_t =  \\log\\frac{\\P(x_{1:t-1}| \\H)}{\\P(x_{1:t-1}| \\T)} + \\log\\frac{\\P(x_t| \\H)}{\\P(x_t| \\T)} = \\LLR_{t-1} + \\log\\frac{\\P(x_t| \\H)}{\\P(x_t| \\T)}.\n",
    "$\n",
    "\n",
    "Note how the increment in $\\LLR_t$ depends only on the outcome of the current/last coin flip, $x_t$. Since we are assuming \n",
    "that the probability of heads is either $p = 0.5 + c$ (under the H hypotheses) or $p = 0.5 - c$ (under the T hypothesis), \n",
    "then the increment, $\\log\\frac{\\P(x_t| \\H)}{\\P(x_t| \\T)}$, is given either by \n",
    "\n",
    "$\\Delta \\equiv \\log\\frac{0.5 + c}{0.5 - c} > 0,$\n",
    "\n",
    "if the coin lands heads (i.e. $x_t = 1$) or by \n",
    "\n",
    "$\\log\\frac{0.5 - c}{0.5 + c} = - \\Delta.$\n",
    "\n",
    "\n",
    "While the agent doesn't know a priori whether H or T is the case, suppose we, as the experimentalist engineers observing her/him, <br>\n",
    "do know the truth. Let us suppose $H$ is true. Then heads have probability $0.5 + c$ and tails has probability $0.5 - c$. <br>\n",
    "This means that $\\LLR_t$ performs a **biased random walk** as follows:\n",
    "\n",
    "in every trial (and independently across trials) \n",
    "it takes a step of size $\\Delta$ either upwards, with probability $0.5 + c$, or downwards, with probability $0.5 - c$. <br>\n",
    "In the case, $c=0$ corresponding to a fair\n",
    "coin: $\\LLR_t$ performs a completley random walk, with equal probability of moving up or down in each trial.<br>\n",
    "\n",
    "It is not hard to see that (under the assumption that H is true) the *average* one-step change in $\\LLR_t$ is \n",
    "\n",
    "$2c \\Delta > 0,$ \n",
    "\n",
    "and the standard deviation (SD) is \n",
    "\n",
    "$\\Delta \\sqrt{1 - 4c^2}.$\n",
    "\n",
    "\n",
    "\n",
    "Finally, to derive a continuous-time drift-diffusion model out of this, assume now that the coherence is very small. \n",
    "In this case we can Taylor expand (using $log x \\approx x$ and $\\sqrt{1 + x}\\approx 1 + x/2$ for small $x$) to obtain $\\Delta \\approx 4c$, and \n",
    "the following expressions for the mean and SD of the incremement $d \\LLR_t$ in the log-likelihood ratio:\n",
    "\n",
    "$\\mathbb{E}[d \\LLR_t] = 8c^2$\n",
    "\n",
    "$\\mathrm{SD}[d \\LLR_t] = 4c$\n",
    "\n",
    "Intuitively, when coherence is low and the coin is very close to a fair coin, the agent has to observe many coin flips<br>\n",
    "in order to be able to make a decision with a decent degree of confidence. If we ignore the stochasticity of the increments\n",
    "and assume it changes by its average value $\\mathbb{E}[d \\LLR_t]$, the number of trials we need to observe to reach the same\n",
    "level of confidence scales like one over this average increment, and thus goes like $1/c^2$. Let us now suppose that each coin flip takes a \n",
    "small time, $d t$, to perform, and let us scale this so that decisions are always made at a time that is on the order of 1. \n",
    "The total time is $dt$ times the number of coin flips, while we have argued that the latter \n",
    "(for reasonably confident decisions) should grow like $1/c^2$ as coherence drops. Since we are requiring that total time stays order 1 and\n",
    "thus grow or shrink as $c$ grows or shrinks, then we need to scale $dt$ like $c^2$. Let me simply define $dt = 8\\beta^{-1} c^2$, for some (order 1) constant $\\beta$.<br>\n",
    "Then we find\n",
    "\n",
    "$\\mathbb{E}[d \\LLR_t] = \\beta dt$\n",
    "\n",
    "$\\mathrm{SD}[d \\LLR_t] = \\sigma \\sqrt{dt}$\n",
    "\n",
    "where $\\sigma^2 = 2 \\beta$. We can then approximate the increments in $\\LLR_t$ as\n",
    "\n",
    "$\n",
    "d\\LLR_t = \\beta dt + \\sigma \\sqrt{dt} \\epsilon_t\n",
    "$\n",
    "\n",
    "where $\\epsilon_t$ is a random variable with mean zero and standard deviation 1.<br> \n",
    "Importantly $\\epsilon_t$ at diffrent time steps, $t$, are independent.  We have thus arrive at equations describing a drift-diffusion stochastic process."
   ]
  }
 ],
 "metadata": {
  "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}