Neurons.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 2. Neural Modeling and Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 2.1 Biophysical neuron models: \n",
    "    - HH model, IaF model, alpha function\n",
    "* 2.2 Abstract neuron models: f-I curve, sigmoid, MP neuron  \n",
    "* 2.3 Models of plasticity: Hebb rule, STDP  \n",
    "* 2.4 Theories of neural coding: rate, population, coincidence, synchrony, waves  \n",
    "* 2.5 Methods of neural decoding: PESH, PESC, STA, STC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from scipy.integrate import odeint, ode"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.1 Biophysical neuron models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Among variety of cells that compose animal body, what is striking about neurons is their complex shapes: branches of dendrites and long-projecting axons. \n",
    "For example, check [NeuroMorpho.org](http://neuromorpho.org) for 3D morphological data of thousansa of neurons.\n",
    "There are non-neural cells that show electric activities and chemical signaling, but neurons are specialized for collecting signals from thousands of different neurons and sending the output to target neurons far apart. \n",
    "\n",
    "Here we introduce basic mathematical models that capture biophysical properties of neurons, namely, electirc excitation, synaptic transmission, and dendritic integtion."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hodgkin-Huxley neuron models\n",
    "The Hodgkin-Huxley (HH) model considers a neuron as an electric circuit as depicted below.\n",
    "\n",
    ":::{figure-md} Fig_HH\n",
    "<img src=\"figures/hhmodel.jpg\" alt=\"HH model\" width=\"500px\">\n",
    "\n",
    "The electric circuit diagram for the Hodgkin-Huxley model.\n",
    ":::\n",
    "\n",
    "On the cellular membrane, there are *ionic channels* that pass specific type of ions. Sodium ions (Na$^+$) are scarce inside the cell, so that when sodium channel opens, positive charges flood into the cell to cause excitation. Potassium ions (K$^+$) are rich inside the cell, so that when potassium channel opens, positive charges flood out of the cell to cause inhibition. The HH model assumes a 'leak' current that put together all other ionic currents.\n",
    "\n",
    "The ingeniety of Hodgkin and Huxley is that they inferred from careful data analysis that a single sodium channel consists of three *activation* gates and one *inactivation* gate, and a single potassium channel consists of four activation gates. Such structures were later confirmed by genomics and imaging.\n",
    "\n",
    "The electric potential inside the neuron $V$ follows the following equation:\n",
    "\n",
    "$$C \\frac{dV}{dt} = g_{Na}m^3h(E_{Na}-V) + g_Kn^4(E_K-V) + g_L(E_L-V) + I$$\n",
    "\n",
    "Here, $m$, $h$, and $n$ represent the proportions of opening of sodium activation, sodium inactivation, and potassium activation gates, respectively. \n",
    "They follow the following differential equations with their rates of opening and closing, $\\alpha(V)$ and $\\beta(V)$, depending on the membrane voltage $V$.\n",
    "\n",
    "$$\\frac{dm}{dt} = \\alpha_m(V)(1-m) - \\beta_m(V)m$$\n",
    "\n",
    "$$\\frac{dh}{dt} = \\alpha_h(V)(1-h) - \\beta_h(V)h$$\n",
    "\n",
    "$$\\frac{dn}{dt} = \\alpha_n(V)(1-n) - \\beta_n(V)n$$\n",
    "\n",
    "These compose a system of four-dimensional non-linear differential equations. Another amazing thing about Hodgkin and Huxley is that they could simulate the solutions of these differential equations by a hand-powered computer.\n",
    "\n",
    "Below is a code to simulate the HH model by Python. Much easier!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# HH: Hodgkin-Huxley (1952) model\n",
    "C = 1.    # membrane capacitance (uF/cm^2)\n",
    "# maximum conductances (uS/cm^2)\n",
    "gna = 120.  # sodium\n",
    "gk = 36.    # potassium\n",
    "gl = 0.3   # leak\n",
    "# reversal potentials (mV)\n",
    "Ena = 50.   # sodium\n",
    "Ek = -77.   # potassium\n",
    "El = -54.4 # leak\n",
    "def hh(y, t, stim=0.):    \n",
    "    # state variables: potential and activation/inactivation\n",
    "    v, m, h, n = y\n",
    "    # membrane potential\n",
    "    if callable(stim):\n",
    "        I = stim(t)  # time-dependent\n",
    "    else:\n",
    "        I = stim  # constant\n",
    "    dv = (gna*m**3*h*(Ena-v) + gk*n**4*(Ek-v) + gl*(El-v) + I)/C\n",
    "    # sodium current activation\n",
    "    am = 0.1*(v+40)/(1-np.exp(-(v+40)/10))\n",
    "    bm = 4*np.exp(-(v+65)/18);\n",
    "    dm = am*(1-m) - bm*m\n",
    "    # sodium current inactivation\n",
    "    ah = 0.07*np.exp(-(v+65)/20)\n",
    "    bh = 1/(1+np.exp(-(v+35)/10))\n",
    "    dh = ah*(1-h) - bh*h\n",
    "    # potassium current activation\n",
    "    an = 0.01*(v+55)/(1-np.exp(-(v+55)/10))\n",
    "    bn = 0.125*np.exp(-(v+65)/80)\n",
    "    dn = an*(1-n) - bn*n\n",
    "    return [ dv, dm, dh, dn]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See the response to a ramp (gradually increasing) current."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# current stimulus I (uA/cm^2)\n",
    "def ramp(t):\n",
    "    return 0.1*np.array(t)  # ramp current"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run a simulation\n",
    "tt = np.arange(0, 300, 0.1)  # time to be simulated\n",
    "y0 = [ -70, 0.1, 0.6, 0.4]  # initial state: V, m, h, n\n",
    "yt = odeint(hh, y0, tt, args=(ramp,))   # simulated output\n",
    "# plot in separate rows\n",
    "plt.subplot(2, 1, 1)\n",
    "plt.plot(tt, ramp(tt));  # stim\n",
    "plt.plot(tt, yt[:,0]);   # Ie, V\n",
    "plt.ylabel(\"I, V\");\n",
    "plt.legend((\"I(t)\", \"V(t)\"), loc='upper left');\n",
    "plt.subplot(2, 1, 2)\n",
    "plt.plot(tt, yt[:,1:]);  # m, h, n\n",
    "plt.xlabel(\"t\")\n",
    "plt.ylabel(\"m, h, n\");\n",
    "plt.legend((\"m(t)\", \"h(t)\", \"n(t)\"), loc='upper left');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Integrate-and-fire neuron models\n",
    "\n",
    "The HH model explaines why a spike is generataed and followed by a refractory period based on the activation and inactionvation of sodium and potassium channels.\n",
    "We are, however, often interested in how a spike is triggered, rather than and electric mechanisms to create a spike. \n",
    "\n",
    "*Integrate-and-fire* (IaF) neurons models considers the dynamics of sub-threshold accumulation of input currents up to a spiking threshold, but take a spike as just an event and then reset the potential. From the HH equation, by removing the sodium and potassium current, we have\n",
    "\n",
    "$$C \\frac{dV(t)}{dt} = g_L(E_L-V(t)) + I$$\n",
    "\n",
    "By defining the membrane resistance $R=\\frac{1}{g_L}$ and the time constant $\\tau=RC$, we have the equation of *leaky integration*\n",
    "\n",
    "$$\\tau\\frac{dV(t)}{dt} = -V(t) + E_L + R I(t)$$\n",
    "\n",
    "When the membrane potential $V$ reaches to a threhold $V_\\theta$, a spike is generated and the membrane potential is reset to $V_r$.\n",
    "\n",
    "$$V(t) = V_r  \\quad\\mbox{ if } V(t) \\ge V_\\theta$$\n",
    "\n",
    "Below is an example of an IaF model. Here we use `ode()` function to detect reaching to the threshold by the `solout` option."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tau = 20\n",
    "R = 1\n",
    "El = -50  # resting potential\n",
    "Vth = -40  # threshold\n",
    "Vs = 40  # spike height\n",
    "Vr = -80  # reset\n",
    "def integ(t, v, stim=0.):  # integrate\n",
    "    if callable(stim):\n",
    "        I = stim(t)  # time-dependent\n",
    "    else:\n",
    "        I = stim  # constant\n",
    "    return (-v + El + R*I)/tau\n",
    "def fire(t, v):  # fire\n",
    "    return (-1 if v>=Vth else 0)  # stop if v>=V1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iaf = ode(integ)\n",
    "iaf.set_integrator('dopri5')  # Runke-Kutta with step size control\n",
    "iaf.set_solout(fire)  # stop when threshold is reached\n",
    "iaf.set_initial_value(Vr, 0)\n",
    "iaf.set_f_params(ramp)  # stimulus\n",
    "V = np.array([Vr])\n",
    "T = [0]\n",
    "tend = 300\n",
    "dt = 10\n",
    "while iaf.successful() and iaf.t <= tend:\n",
    "    v = iaf.integrate(iaf.t+dt)  # step\n",
    "    #V.append(v)       # record V\n",
    "    V = np.append(V, v)       # record V\n",
    "    T.append(iaf.t)   # record t\n",
    "    if v >= Vth:  # threshold reached\n",
    "        iaf.set_initial_value(Vr, iaf.t)  # reset V\n",
    "        V = np.append(V, Vs)  # spike\n",
    "        T.append(iaf.t)\n",
    "        V = np.append(V, Vr)  # reset\n",
    "        T.append(iaf.t)\n",
    "plt.plot(T, ramp(T))\n",
    "plt.plot(T, V);\n",
    "plt.xlabel(\"t\")\n",
    "plt.ylabel(\"I, V\");\n",
    "plt.legend((\"I(t)\", \"V(t)\"));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### F-I curve\n",
    "\n",
    "In the HH or IaF models, as we increase the input current $I$, the spike firing frequency $F$ increases. It is possible to characterize the property of a neuron by this $F-I$ relationship."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def iafspikes(Ic):  # spike times with constant current\n",
    "    iaf = ode(integ).set_integrator('dopri5')\n",
    "    iaf.set_solout(fire)  # stop when threshold is reached\n",
    "    iaf.set_f_params(Ic)  # stimulus\n",
    "    iaf.set_initial_value(Vr, 0)\n",
    "    tinit = 200  # initial transient\n",
    "    trun = 1000  # for 1000ms\n",
    "    dt = 1\n",
    "    tf = []  # spike timing\n",
    "    while iaf.successful() and iaf.t <= tinit+trun:\n",
    "        v = iaf.integrate(iaf.t+dt)  # step\n",
    "        if v >= Vth:  # threshold reached\n",
    "            iaf.set_initial_value(Vr, iaf.t)  # reset V\n",
    "            if iaf.t>=tinit:  # after transient\n",
    "                tf.append(iaf.t)\n",
    "    return np.array(tf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Is = np.arange(-20, 51, 5)\n",
    "Fs = []\n",
    "for c, Ic in enumerate(Is):\n",
    "    tf = iafspikes(Ic)\n",
    "    plt.plot(tf, tf*0+Ic, '.')\n",
    "    Fs.append(len(tf))\n",
    "plt.xlabel(\"t\")\n",
    "plt.ylabel(\"I\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(Is, Fs);\n",
    "plt.xlabel(\"I\")\n",
    "plt.ylabel(\"F (Hz)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Noisy integrate-and-fire model\n",
    "Try adding some noise to the neuron."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "In = 20  # noise size\n",
    "noise = np.random.normal(0, 1, 1500)  # normal gaussian noise every 1ms\n",
    "def Inoise(t):\n",
    "    i = int(np.floor(t))\n",
    "    return Ic + In*noise[i]\n",
    "Fs = []\n",
    "for c, Ic in enumerate(Is):\n",
    "    tf = iafspikes(Inoise)\n",
    "    #plt.plot(tf, tf*0+c, '.')\n",
    "    Fs.append(len(tf))\n",
    "plt.plot(Is, Fs);\n",
    "plt.xlabel(\"I\")\n",
    "plt.ylabel(\"F (Hz)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The firing threshold becomes smoother with noise."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear-nonlinear-Poisson models\n",
    "\n",
    "In some mathematical analysis, it is convenient to assume that a neuron collects synaptic inputs by linear weighted sum, sets its firing rate through a nonlinear function taking non-negative values, such as exponential and sigmoid function, and produces spikes stochastically according to its instantaneous firing rate.\n",
    "\n",
    "This is a *Poisson* process, for which the number of spiked generated in a time bin follows the Poisson distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Alpha function model of synaptic current\n",
    "\n",
    "When a spike travels through an axon and reaches to a synpase, it causes release of neurotransmitter molecules in the synaptic junction, which binds to the receptors on the postsynapic membrane, and cause either opening of ionic channels of the receptor or trigger molecular reactions in the postsynaptic cell.\n",
    "\n",
    "A simple model to approximate the dynamics of this complex cascades is the *alpha function*, defined as\n",
    "\n",
    "$$u(t) = \\frac{t}{\\tau_s} e^{-\\frac{t}{\\tau_s}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def alpha(t, taus):\n",
    "    return t/taus*np.exp(-t/taus)\n",
    "t = np.arange(0, 100, 1)\n",
    "taus = 10\n",
    "plt.plot(t, alpha(t, taus));\n",
    "plt.xlabel(\"t\")\n",
    "plt.ylabel(\"u(t)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The alpha function is a solution of a second-order dyamics to an impulse input at t=0\n",
    "\n",
    "$$\\tau_s\\frac{du_1(t)}{dt} = -u_1(t) + \\delta(t=0)$$\n",
    "\n",
    "$$\\tau_s\\frac{du_2(t)}{dt} = -u_2(t) + u_1(t)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def falpha(u, t):\n",
    "    # dynamics for alpha function. u = [u1, u2]\n",
    "    return np.array([-u[0], u[0] - u[1]])/taus\n",
    "ut = odeint(falpha, [1, 0], t)  # impulse input at t=0\n",
    "plt.plot(t, ut)\n",
    "plt.plot(t, alpha(t, taus), ':');  # analytic form in dotted line\n",
    "plt.xlabel(\"t\")\n",
    "plt.ylabel(\"u(t)\")\n",
    "plt.legend((\"u1\", \"u2\", \"u\"));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us make a network of IaF neurons with alpha function synapses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# IaF neuron network with alpha function synapses\n",
    "#I = 0   # bias current\n",
    "# Iaf neurons\n",
    "N = 10  # number of neurons\n",
    "tau = 50  # cellular time constant\n",
    "V0 = -45  # resting near threshold\n",
    "Vth = -40  # threshold\n",
    "Vs = 40  # spike height\n",
    "Vr = -80  # reset\n",
    "# Alpha synapses\n",
    "taus = 5  # synaptic time constant\n",
    "W = 40   # connection weight size; assume exponential distribution\n",
    "w = np.random.exponential(W, N**2).reshape((N,N))\n",
    "for i in range(N):\n",
    "    w[i,i] = 0  # remove self-connection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def integnet(t, vu):  # integrate\n",
    "    # vu = [v, u]\n",
    "    v = vu[:N]\n",
    "    u = vu[N:].reshape((N,N,2))\n",
    "    # synaptic potential\n",
    "    epsp = np.sum(w*u[:,:,1], axis=1)  # sum rows\n",
    "    # membrane dynamics\n",
    "    dv = (-v + V0 + epsp)/tau\n",
    "    # synaptic dynamics: for uniform taus, it can be reduced to N dim\n",
    "    du = np.stack((-u[:,:,0]/taus, (u[:,:,0]-u[:,:,1])/taus), axis=2)\n",
    "    return np.concatenate([dv, du.flatten()])\n",
    "def firenet(t, vu):  # fire\n",
    "    return (-1 if sum(vu[:N]>=Vth) else 0)  # stop if any of v>=V1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run a network simulation\n",
    "iafnet = ode(integnet).set_integrator('dopri5')\n",
    "iafnet.set_solout(firenet)  # stop when threshold is reached\n",
    "v0 = np.random.uniform(Vr, Vth, N)  # initial state\n",
    "u0 = np.random.uniform(0, 1, N*N*2)  # synapse state\n",
    "iafnet.set_initial_value(np.append(v0, u0, 0))\n",
    "V = [v0]\n",
    "T = [0]\n",
    "tend = 500\n",
    "dt = 1\n",
    "while iafnet.successful() and iafnet.t <= tend:\n",
    "    vu = iafnet.integrate(iafnet.t+dt)  # new v and u\n",
    "    V.append(vu[:N].copy())   # record V\n",
    "    T.append(iafnet.t)    # record t\n",
    "    if sum(vu[:N] >= Vth):  # any reached threshold\n",
    "        for i in range(N):   # check each neuron\n",
    "            if vu[i] >= Vth:  # reached threshold\n",
    "                vu[i] = Vr   # reset potential\n",
    "                vu[N:].reshape(N,N,2)[:,i,0] += 1  # synaptic input\n",
    "                V[-1][i] = Vs  # spike\n",
    "        iafnet.set_initial_value(vu, iafnet.t)  # reset state\n",
    "V = np.array(V)   # list to array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(T, V + np.arange(N)*100);  # shift traces vertically\n",
    "plt.xlabel(\"t\")\n",
    "plt.ylabel(\"V(t)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Abstract neuron models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mean-field models\n",
    "\n",
    "In the invertebrate neuvous system with small number of neurons, each single neuron has a distinct identity and signals specifit information. \n",
    "In the mammarian brain with millons to billions of neurons, on the other hand, neural responses tend to be redundant and distributed. \n",
    "For example, in the visual cortex, neurons in the same *column* have similar sensory tuning, such as presentation of edges in the same orientation.\n",
    "Moreover, responses of each neuron to the same stimulus can be variable across trials.\n",
    "This suggest that information is reliably represented by a population of neurons sharing the same tuning.\n",
    "\n",
    "*Mean-field* models, or neural *mass* model, capture the average firing rate of each population of neurons.\n",
    "\n",
    "The activity of a population of $N$ neurons is defined as\n",
    "\n",
    "$$A(t) = \\lim_{\\Delta t\\rightarrow 0} \\frac{n(t;t+\\Delta t)}{N\\Delta t}$$\n",
    "\n",
    "where $n(t;t+\\Delta t)$ is the number of spikes in the short time duration $\\Delta t$.\n",
    "\n",
    "For a population of identical IaF neurons with random homogeneous connections, the pupulation activity $A(t)$ and average membrane potential $h(t)$ can be apprximated by the following equations (Gerstner et al., 2014):\n",
    "\n",
    "$$A(t) = F(h(t))$$\n",
    "\n",
    "$$\\tau \\frac{dh(t)}{dt} = -h(t) + R I(t)$$\n",
    "\n",
    "The gain function $F$ is often approximated by a sigmoid function\n",
    "\n",
    "$$F(h) = \\frac{1}{1 + e^{-a(h-\\theta)}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Wilson-Cowan model\n",
    "A classic example of a mean-field model is the Wilson-Cowan model, which consists of excitatory and inhibitory populations connected to each other.\n",
    "Depending on the connection strengths, we can observe point attractor or oscillating behaviors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mean field network model\n",
    "N = 2  # number of populations\n",
    "w = np.array([[10, -10], [10, 0]])\n",
    "theta = np.array([-2.5, 2.5])\n",
    "tau = 10\n",
    "def gain(h):\n",
    "    h = np.array(h)  # to accept list input\n",
    "    return 1/(1 + np.exp(-h))\n",
    "def mfn(h, t):\n",
    "    return (-h + np.dot(w, gain(h)) - theta)/tau"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tt = np.arange(0, 200, 0.1)  # time to be simulated\n",
    "h0 = np.random.uniform(-1, 1, N)  # initial state\n",
    "ht = odeint(mfn, h0, tt)   # simulated output\n",
    "plt.plot(tt, gain(ht));   # A(t) = F(h(t))\n",
    "plt.xlabel(\"t\")\n",
    "plt.ylabel(\"A(t)\")\n",
    "plt.legend((\"excitatory\", \"inhibitory\"));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Artificial neural network models\n",
    "\n",
    "Abstraction of linear weighted input and nonlinear gain function from those neuron models brings us to *artificial neural networks* or *connectionist models*.\n",
    "\n",
    "$$x_i(t+1) = \\sum_{j=1}^n w_{ij} f(x_j(t))$$\n",
    "\n",
    "This is a recurrent neural network. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Illustration of a recurrent neural network\n",
    "n = 5  # number of neurons\n",
    "W = np.random.normal(0, 1, n*n).reshape(n,n)  # random weights\n",
    "i = np.arange(0,n)  # neuron index\n",
    "pos = np.array([np.sin(2*np.pi*i/n), np.cos(2*np.pi*i/n)]).T\n",
    "for i in range(n):  # target\n",
    "    for j in range(n):  # source\n",
    "        # edges from sources shifted inside, to target outside\n",
    "        plt.plot(pos[[j,i],0]*[.9,1.1], pos[[j,i],1]*[.9,1.1],\n",
    "                 lw=2*abs(W[i,j]), color=('r' if W[i,j]>0 else 'b'))\n",
    "plt.plot(pos[:,0], pos[:,1], 'oy', markersize=30)\n",
    "plt.axis('equal')\n",
    "plt.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Feed-forward neural networks\n",
    "\n",
    "A popular sub-class of artificial neural networks is *feed-forward* networks, which are organized by multiple layers and connections are uni-directional from lower to upper layers.\n",
    "\n",
    "$$x_i^{l+1} = \\sum_{j=1}^n w_{ij}^l f(x_j^l)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Illustration of a feed-forward neural network\n",
    "n = [4, 5, 3]  # numbers of neurons in layers\n",
    "W = []  # empty list\n",
    "for l in range(len(n)-1):  # weights from bottom to top layers\n",
    "    W.append(np.random.normal(0, 1, n[l+1]*n[l]).reshape(n[l+1],n[l]))\n",
    "    for i in range(n[l+1]):  # target\n",
    "        for j in range(n[l]):  # source\n",
    "            plt.plot([(j+1)/(n[l]+1), (i+1)/(n[l+1]+1)], [l, l+1],\n",
    "                 lw=2*abs(W[l][i,j]), color=('r' if W[l][i,j]>0 else 'b'))\n",
    "    plt.plot((np.arange(n[l])+1)/(n[l]+1), np.repeat(l,n[l]), 'oy', markersize=20)\n",
    "plt.plot((np.arange(n[-1])+1)/(n[-1]+1), np.repeat(len(n)-1,n[-1]), 'oy', markersize=20)\n",
    "plt.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Activation functions: Sigmoid, ReLU, and binary\n",
    "\n",
    "As the activation (or gain) function $f$, the most popular one is the (logistic) sigmoid function.\n",
    "\n",
    "$$f(x) = \\frac{1}{1 + e^{-x}}$$\n",
    "\n",
    "The use of *rectified linear unit* (ReLU) is also becoming popular for deep layered networks.\n",
    "\n",
    "$$f(x) = \\max(x, 0)$$\n",
    "\n",
    "The classic McCulloch-Pitts model took a binary activation function.\n",
    "\n",
    "$$f(x) = \\begin{cases} 1 & \\mbox{if } x \\ge 0 \\\\ 0 & \\mbox{if } x<0 \\end{cases}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Activation functions\n",
    "def sigmoid(x):\n",
    "    return 1/(1+np.exp(-x))\n",
    "def relu(x):\n",
    "    return np.maximum(x, 0)\n",
    "def binary(x):\n",
    "    return x>=0\n",
    "x = np.arange(-4, 4.1, 0.1)\n",
    "plt.plot(x, sigmoid(x))\n",
    "plt.plot(x, relu(x))\n",
    "plt.plot(x, binary(x))\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"f(x)\")\n",
    "plt.legend((\"sigmoid\", \"ReLU\", \"binary\"));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Softmax function\n",
    "When the output of a layer represents a proabaility distribution, it is common to use the *softmax* function whose outputs sum up to one.\n",
    "\n",
    "$$f_i(x) = \\frac{e^{x_i}}{\\sum_j e^{x_j}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Softmax function\n",
    "def softmax(x, beta=1):\n",
    "    \"\"\"softmax function: x is a vector, or column vectors\"\"\"\n",
    "    ex = np.exp(beta*np.array(x))\n",
    "    # return ex/np.sum(ex)   # for one vector\n",
    "    return ex/np.sum(ex, axis=0)  # for column vectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = np.linspace(-5, 5)\n",
    "x = np.stack( (-u, u*0+2, u))  # three components\n",
    "plt.subplot(2, 1, 1)\n",
    "plt.plot(u, x.T)\n",
    "plt.ylabel(\"x\");\n",
    "plt.subplot(2, 1, 2)\n",
    "plt.plot(u, softmax(x).T)   # constrained in [0, 1]\n",
    "plt.xlabel(\"u\")\n",
    "plt.ylabel(\"softmax(x)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Models of plasticity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hebb rule\n",
    "\"Cells fire together wire together\" is the basic concept proposed by Donald Hebb [Hebb1952]. More specifically, the Hebbian synaptic plasticity rule takes the form\n",
    "\n",
    "$$\\Delta w_{ij}(t) = \\alpha y_i(t) y_j(t)$$\n",
    "where $\\alpha$ is the learning rate parameter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spike timining dependent placticity (STDP)\n",
    "It has been observed in hippocampus, cortex and other networks that synaptic plasiticity is dependent on the timing of pre- and post-synapic spikes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.linspace(-100, 100)  # post-pre spike time difference\n",
    "tau = 20  # time constant\n",
    "dw = (t<0)*(-np.exp(t/tau)) + (t>0)*np.exp(-t/tau)  # depression if t<0\n",
    "plt.plot(t, dw)\n",
    "plt.grid(True)\n",
    "plt.xlabel(\"post-pre spike time difference\")\n",
    "plt.ylabel(\"synaptic weight change\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Theories of neural coding\n",
    "\n",
    "It is certain that neurons carry sensory, motor, or any cognitive information and perform computation by combining and transforming such information.\n",
    "But how exectly do they code a variety of information? \n",
    "This is not a trivial problem and there are many theories and debates.\n",
    "\n",
    "### Rate coding\n",
    "Each neuron encodes a certain variable by its firing rate. This is most evident in sensory receptor neurons and motor neurons.\n",
    "For example, the firing rate of a retinal photoreceptor is monotinically related to the strength of the light hitting the cell.\n",
    "\n",
    "### Population coding\n",
    "For motor control and cognitive processing, the brain has to combine multiple modalities of information, such as vision, hearing, touch and proprioception.\n",
    "To represent particular combinations of such information, some kind of multi-dimensional non-linear basis functions are required.\n",
    "\n",
    "*Population coding* is an idea in which a group of neurons with different response tuning functions represent information by their activity patterns.\n",
    "The recepient neurons can extract specific information by weighted sum of their activities.\n",
    "\n",
    "### Temporal coding\n",
    "In the rate coding framework, what matters is the frequency of spikes of neurons coding the same or related information.\n",
    "But there is a possiblity that not just the frequency, but the timing of each spike carry a certain information.\n",
    "\n",
    "For example, it is known that spikes of some auditory neurons are produced a certain phase of sound wave, which can be helpful for identification of sound source direction by detecting the time difference of the sound arriving in left and right ears.\n",
    "\n",
    "A more sophisticated idea is that even at the same firing frequencies, neurons represent related information by the synchrony of spikes.\n",
    "Experimental evidence suggests that visual cortex neurons that respond to line segments in a similar orientation spike synchronously for a single connected line, but asynchronously for separate line segments."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Methods for spike analysis/decoding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Peristimulus time histogram (PSTH)\n",
    "A neuron's repsonse to the same sensory stimulus can vary trial-by-trial.\n",
    "While primary sensory neurons respond reliably to sensory stimuli, sensory responses of cortical neurons are known to be highly variable across trials.\n",
    "\n",
    "To characterize the average response properties of a neuron, the most typical method is to align the spike trains from different trials at the stimulus onset and count the number of spikes in time bins of tens to hundreds of milliseconds.\n",
    "The original spike trains in multiple trials are often called *raster plot* and the average spike frequency around the time of stimulus onset is called *peristimulus time histogram (PSTH)*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spike-triggered average (STA)\n",
    "Characterizing the response property of a neuron the other way round is to see what happened before a spike is generated.\n",
    "That is *spike-triggered average (STA)* of sensory stimuli."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Neurons models\n",
    "\n",
    "* Hodgkin AL, Huxley AF (1952) A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol, 117(4):500–544. https://doi.org/10.1113/jphysiol.1952.sp004764\n",
    "\n",
    "* Gerstner W, Kistler WM, Naud R, Paninski L (2014) Neuronal Dynamics: From Single Neurons to Networks and Models of Cognition. [Cambridge University Press](https://www.cambridge.org/core/books/neuronal-dynamics/75375090046733765596191E23B2959D). \n",
    "([online version and Python exercise](http://neuronaldynamics.epfl.ch/))\n",
    "\n",
    "* Wilson HR, Cowan JD (1972) Excitatory and inhibitory interactions in localized populations of model neurons. Biophys. J., 12:1–24. https://doi.org/10.1016/S0006-3495(72)86068-5\n",
    "\n",
    "### Neural coding and analysis\n",
    "\n",
    "* Dayan P, Abott LF (2001) Theoretical Neuroscience: Compuational and\n",
    "Mathematical Modeling of Neural Systems. MIT Press. https://mitpress.mit.edu/9780262041997/theoretical-neuroscience/\n",
    "\n",
    "* Gray CM, Konig P, Engel AK, Singer W (1989) Oscillatory responses in cat visual cortex exhibit inter-columnar synchronization which reflects global stimulus properties. Nature, 338:334–337. https://doi.org/10.1038/338334a0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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": 2
}