8-Spiking-Neural-Networks.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Notebook 8: Spiking Neural Networks\n",
    "\n",
    "In Notebook 7, we used *Brian 2* to simulate individual neurons and explore their dynamical properties. This week we introduced spiking neural networks (SNNs) and how they can be used for biological simulation and machine learning problems. The goals of this notebook are:\n",
    "\n",
    "+ Learn additional features of *Brian 2* that are useful in making SNNs.\n",
    "+ How to create populations of neurons with heterogeneous dynamics.\n",
    "+ Create and manipulate synapses between neurons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install brian2\n",
    "!pip install brian2tools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from brian2 import *\n",
    "from brian2tools import *\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note**: Depending on the version of brian that gets installed, some of the code cells may generate a lot of log messages. You can safely\n",
    "ignore any of these that don't explicitly say there was an error."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulating Many Neurons\n",
    "\n",
    "It is easy to simulate large groups of neurons using *Brian 2*. We can add randomness to parameters of the neuron models to produce heterogeneous networks. Let's simulate 20 LIF neurons that have a random $\\tau$ parameter between 0-10 ms and constant input $I(t) = 2$. The **rand()** function used in a string in *Brian 2* will produce a random number between 0-1 for each neuron in the network. We can also set up monitors for the **v** state variables, the spike times of each neuron, and the overall population firing rate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_scope()\n",
    "\n",
    "# Network Parameters\n",
    "sim_length = 30*ms\n",
    "n_neurons = 20\n",
    "\n",
    "# Model Parameters\n",
    "I = 2\n",
    "theta = 1\n",
    "\n",
    "# Model Equations\n",
    "eqs = '''\n",
    "dv/dt = (-v+I)/tau:1\n",
    "tau:second #assigning tau as a state variable allows us to easily access it for introducing randomness\n",
    "'''\n",
    "\n",
    "# Create Network\n",
    "G = NeuronGroup(n_neurons,eqs, threshold = 'v > theta', reset = 'v = 0', method = 'euler')\n",
    "G.v = 0\n",
    "G.tau = '10*rand()*ms'\n",
    "\n",
    "# Monitor Network States\n",
    "M = StateMonitor(G,\"v\", record = True) #records state variable v for each neuron\n",
    "spkM = SpikeMonitor(G) #records spike times of each neuron in network\n",
    "Rates = PopulationRateMonitor(G) #records population firing rate\n",
    "\n",
    "run(sim_length)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can easily visualize the output of network by accessing our network monitors (the *brian_tools* package comes with many useful visualization features we will occasionally use)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=2, ncols=1, sharex=True, figsize=(9, 4))\n",
    "ax[0].plot(spkM.t/ms,spkM.i, '|',color = 'k')\n",
    "ax[0].set_title(\"Population Raster Plot\")\n",
    "ax[0].set_ylabel(\"Neuron Index\")\n",
    "\n",
    "brian_plot(Rates, axes=ax[1])\n",
    "ax[1].set_title(\"Population PSTH\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also look at the internal state variable for any of the neurons in the network. Here are the first 3 neurons:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_plot = 3\n",
    "fig, axes = plt.subplots(nrows = n_plot, ncols = 1, sharex=True, figsize=(9,6))\n",
    "for idx, ax in enumerate(axes):\n",
    "    ax.plot(M.t/ms,M.v[idx])\n",
    "    ax.set_ylabel(f\"V_{idx}\")\n",
    "ax.set_xlabel(\"Time (ms)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1\n",
    "\n",
    "**A**. Create a network with the following properties:\n",
    "\n",
    "+ 30 QIF neurons with a threshold of 40 and reset of 0.5 (look in the previous notebook for this model)\n",
    "+ random $\\tau$ between 5 and 10ms\n",
    "+ random initial **v** condition between 0 and 1\n",
    "+ constant input of 1\n",
    "\n",
    "**B**. Plot the PSTH, rasters, and some **v** state variables of the network."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Synapses\n",
    "\n",
    "Right now our network is totally unconnected: each neuron spikes without any dependence on the others. To make things more interesting, we need to create synapses between the neurons in our network. We will start out with a basic synapse. When the pre-synaptic neuron fires, the voltage of the post-synaptic neuron will jump by the synaptic weight value between the two neurons. Let's simulate a network identical to the one in the first example, but now with a synapse connectivity matrix with the following properties:\n",
    "\n",
    "+ the weights of each synapse will follow a Gaussian distribution with a mean of 0 (i.e., equal excitatory and inhibitory connections) and standard deviation of 0.1\n",
    "+ each neuron will have a probability of 0.5 to be synapsed to any other neuron\n",
    "+ no self-synapses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_scope()\n",
    "sim_length = 30*ms\n",
    "\n",
    "#Model Parameters\n",
    "n_neurons = 20\n",
    "I = 2\n",
    "theta = 1\n",
    "\n",
    "#Model Equations\n",
    "eqs = '''\n",
    "dv/dt = (-v+I)/tau:1\n",
    "tau:second\n",
    "'''\n",
    "\n",
    "#Create Network\n",
    "G = NeuronGroup(n_neurons,eqs, threshold = 'v > theta', reset = 'v = 0', method = 'euler')\n",
    "G.v = 0\n",
    "G.tau = '10*rand()*ms'\n",
    "\n",
    "#Create Synapses\n",
    "S = Synapses(G,G,\"w:1\",on_pre='v += w') #create synapses from neuron group G to neuron group G\n",
    "S.connect(p = 0.5, condition = 'i!=j') #connects synapses with optional conditons\n",
    "S.w = np.random.normal(0,.1,size = len(S)) #random matrix with size equal to number of synapses needed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the connectivity matrix of the network:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "brian_plot(S.w);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll run the simulation and plot the resulting activity:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Monitor Network State\n",
    "M = StateMonitor(G,\"v\", record = True)\n",
    "spkM = SpikeMonitor(G)\n",
    "Rates = PopulationRateMonitor(G)\n",
    "Syn_M = StateMonitor(S,'w',record = True) #monitor for synapses\n",
    "\n",
    "run(sim_length)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=2, ncols=1, sharex=True, figsize=(9, 4))\n",
    "ax[0].plot(spkM.t/ms,spkM.i, '|',color = 'k')\n",
    "ax[0].set_title(\"Population Raster Plot\")\n",
    "ax[0].set_ylabel(\"Neuron Index\")\n",
    "\n",
    "brian_plot(Rates, axes=ax[1])\n",
    "ax[1].set_title(\"Population PSTH\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_plot = 3\n",
    "fig, axes = plt.subplots(nrows = n_plot, ncols = 1, sharex=True, figsize=(9,6))\n",
    "for idx, ax in enumerate(axes):\n",
    "    ax.plot(M.t/ms,M.v[idx])\n",
    "    ax.set_ylabel(f\"V_{idx}\")\n",
    "ax.set_xlabel(\"Time (ms)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2\n",
    "\n",
    "**A**. Create a network that is composed of two \"areas\", i.e. you will make two neuron groups (G1, G2). \n",
    "\n",
    "Let G1 have the following properties:\n",
    "\n",
    "+ Composed of 20 LIF neurons with a constant input of 2.\n",
    "+ A random $\\tau$ between $0$ and $5ms$.\n",
    "+ Each neuron in G1 synapses with every neuron in G2.\n",
    "+ The weights of these synapses are normally distributed with a mean of 0 and a standard deviation of 2.\n",
    "\n",
    "Let G2 have the following properties:\n",
    "\n",
    "+ Composed of 40 LIF neurons (no direct input).\n",
    "+ A random $\\tau$ between $0$ and $10ms$.\n",
    "+ Each neuron in G2 synapses with every neuron in G2 with a probability of 0.5.\n",
    "+ Self-synapses are possible.\n",
    "+ The weights of these synapses are normally distributed with a mean of 0 and standard deviation of 1.\n",
    "\n",
    "**B**. Plot the PSTH, rasters, and synaptic connectivity matricies of G1 and G2.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spike-Timing-Dependent Plasticity\n",
    "\n",
    "One of the ways the brain implements learning is by modifying the strength of synapses. A common model of Hebbian learning is spike-timing-dependent plasticity (STDP). If the pre-synaptic neuron fires before the post-synaptic neuron, the strength of that synapses increases. If the opposite occurs, the synapse will decrease in strength. A simple way to model STDP is with the following equations:\n",
    "\n",
    "$\\Delta w = \\sum_{t_{pre}} \\sum_{t_{post}} W(t_{post} - t_{pre})$\n",
    "\n",
    "$ W(\\Delta t)=   \\left\\{\n",
    "\\begin{array}{ll}\n",
    "      A_{pre}exp(\\frac{-\\Delta t}{\\tau_{pre}}) & \\Delta t >0 \\\\\n",
    "      A_{post}exp(\\frac{-\\Delta t}{\\tau_{post}}) & \\Delta t <0\n",
    "\\end{array} \n",
    "\\right.$\n",
    "\n",
    "where\n",
    "\n",
    "+ $\\Delta w$ is the change in synaptic strength\n",
    "+ $t_{pre}$ and $t_{post}$ are the times of the pre- and post-synaptic spikes respectively\n",
    "+ $\\Delta t$ is $t_{post} - t_{pre}$\n",
    "+ $W(\\Delta t)$ is a function of the timing difference between the pre- and post-synatpic spikes, with parameters $A_{pre}, A_{post}, \\tau_{pre}, \\tau_{post}$\n",
    "\n",
    "The way to algorithmically implement this version of STDP in *Brian 2* is not obvious and as a bonus exercise you should try and figure out why these equations and the algorithm below are equivalent.\n",
    "\n",
    "Let's now simulate a network with STDP. We will use an additional parameter **wmax** to prevent our synapses from blowing up to arbitrarily high numbers. We will use a more complicated synaptic model than in the previous example. Each neuron will now have two inputs: $I_{ext}$ and $I_{syn}$. $I_{ext}$ will be an external input current and $I_{syn}$ will the the current due to synaptic inputs. A new state variable $I_{syn}$ will be added to the neuron equations and will be modeled as an exponential decay. A post-synaptic neuron receiving a spike will increase its $I_{syn}$ by an amount equal to the synapse weight, and this current will decay with parameter $\\tau_{syn}$. Each synapse will also have a delay associated with it (you can imagine this is modeling some notion of distance between neurons in the network).\n",
    "\n",
    "Additonally, let's demonstrate how to have a more complicated input structure into our network. The previous notebook gave you code for a *TimedArray* where you can give a neuron a time-varying input current. We can also use a *TimedArray* to give every neuron in the network a unique time-varying input. In the below code, each neuron receives a sine wave input with a random phase. Each sine wave has Gaussian random noise added to it. We chose a signal length of $5000$ since our simulation length is $500ms$ and our simulation is running at a default time resolution of $0.1ms$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_scope()\n",
    "\n",
    "# Simulation Parameters\n",
    "n_neurons = 20\n",
    "sim_length = 500*ms\n",
    "\n",
    "# STDP Parameters\n",
    "taupre = taupost = 5*ms\n",
    "wmax = 0.01\n",
    "Apre = 0.1\n",
    "Apost = -0.1\n",
    "\n",
    "#Input Signal\n",
    "I_ext = []\n",
    "for m in range(n_neurons):\n",
    "    random_phase = np.random.normal(0,2*pi,size = 1)\n",
    "    I_signal = 2*np.sin(np.arange(5000)*0.01+random_phase)\n",
    "    I_noise = np.random.normal(0,0.5,size = len(I_signal))\n",
    "    I_ext.append(I_signal+I_noise)\n",
    "I_ext = TimedArray(np.array(I_ext).T,dt = 0.1*ms) \n",
    "\n",
    "# Model Parameters\n",
    "theta = 1\n",
    "tau_syn = 2*ms\n",
    "\n",
    "# Model Equations\n",
    "eqs = '''\n",
    "dv/dt = (-v+I_ext(t,i)+I_syn)/tau:1 #notice new index i in the I_ext state variable, this iterates over each neuron\n",
    "dI_syn/dt = -I_syn/tau_syn:1\n",
    "tau:second\n",
    "'''\n",
    "\n",
    "# Create Network\n",
    "G = NeuronGroup(n_neurons,eqs, threshold = 'v > theta', reset = 'v = 0', method = 'euler')\n",
    "G.v = 0\n",
    "G.tau = '(rand()*5+5)*ms'\n",
    "\n",
    "# Create Synapses\n",
    "S = Synapses(G,G,\n",
    "            '''\n",
    "            w:1\n",
    "            dapre/dt = -apre/taupre:1 (event-driven)\n",
    "            dapost/dt = -apost/taupost:1 (event-driven)\n",
    "            ''',\n",
    "            on_pre= '''\n",
    "            I_syn += w\n",
    "            apre += Apre\n",
    "            w = clip(w+apost, -wmax, 0)\n",
    "            ''',\n",
    "            on_post = '''\n",
    "            apost += Apost\n",
    "            w = clip(w+apre, 0, wmax )\n",
    "            ''',\n",
    "            method = 'linear')\n",
    "S.connect(condition = 'i!=j')\n",
    "S.w = np.random.normal(0,.1,size = len(S))\n",
    "S.delay = '(rand()*2)*ms' #random delay in each synapse\n",
    "\n",
    "# Monitor Network States\n",
    "M = StateMonitor(G,\"v\", record = True)\n",
    "spkM = SpikeMonitor(G)\n",
    "Rates = PopulationRateMonitor(G)\n",
    "Syn_M = StateMonitor(S,'w',record = True)\n",
    "\n",
    "\n",
    "run(sim_length)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=2, ncols=1, sharex=True, figsize=(9, 4))\n",
    "ax[0].plot(spkM.t/ms,spkM.i, '|',color = 'k')\n",
    "ax[0].set_title(\"Population Raster Plot\")\n",
    "ax[0].set_ylabel(\"Neuron Index\")\n",
    "\n",
    "brian_plot(Rates, axes=ax[1])\n",
    "ax[1].set_title(\"Population PSTH\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at how the connectivity matrix changes as a result of the plasticity we added to the network. We can access the connectivity matrix at every time point by indexing the synapse monitor that we created. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We have to do some fancy indexing because Syn_M only contains the connected synapses (i.e. no self-synapses)\n",
    "W_int = np.zeros((n_neurons,n_neurons)) #initial connectivity matrix\n",
    "W_middle = np.zeros((n_neurons,n_neurons)) #halfway through the simulation\n",
    "W_final = np.zeros((n_neurons,n_neurons)) #final connectivity matrix\n",
    "\n",
    "W_int[S.i,S.j] = Syn_M.w[:,0]\n",
    "W_middle[S.i,S.j] = Syn_M.w[:,2500]\n",
    "W_final[S.i,S.j] = Syn_M.w[:,-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows = 1, ncols = 3, sharey=True, figsize=(15,5))\n",
    "ax[0].imshow(W_int.T,origin='lower')\n",
    "ax[0].set_title(\"Connectivity Matrix (0ms)\")\n",
    "ax[0].set_xlabel(\"Source Neuron Index\")\n",
    "ax[0].set_ylabel(\"Target Neuron Index\")\n",
    "\n",
    "ax[1].imshow(W_middle.T,origin='lower')\n",
    "ax[1].set_title(\"Connectivity Matrix (250ms)\")\n",
    "ax[1].set_xlabel(\"Source Neuron Index\")\n",
    "\n",
    "ax[2].imshow(W_final.T,origin='lower')\n",
    "ax[2].set_title(\"Connectivity Matrix (500ms)\")\n",
    "ax[2].set_xlabel(\"Source Neuron Index\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3\n",
    "\n",
    "**A**. Create a new network identical to the one above, but with the following changes:\n",
    "\n",
    "+ Increase the number of neurons to 40.\n",
    "+ Choose a new time-varying input function for the network. Add Gaussian noise so each neuron is getting a unique input. \n",
    "\n",
    "*Hint*: You may need to adjust the STDP parameters to prevent your network for blowing up or becoming silent.\n",
    "\n",
    "**B**. Plot the PSTH, rasters, and the synaptic connectivity matrices at $t=$ 0, 100, 300, and 500 ms."
   ]
  },
  {
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}