Network_of_LIF_neurons.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cb16c18d-a05d-4448-b537-4375182566e3",
   "metadata": {},
   "source": [
    "[![Open In Colab](./colab-badge.png)](https://colab.research.google.com/github/subhacom/moose-notebooks/blob/main/Network_of_LIF_neurons.ipynb) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c60110d4-b7fa-4499-8675-af40ef84bf33",
   "metadata": {},
   "source": [
    "If you are using `colab`, in a fresh runtime you need to run pip to install pymoose there. On the other hand, if you are running jupyter locally and have pymoose installed in that environment, skip the line below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "-3jqFhK9wPYq",
   "metadata": {
    "id": "-3jqFhK9wPYq"
   },
   "outputs": [],
   "source": [
    "# Only for colab\n",
    "!pip install pymoose --quiet"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99d7c0b4-0ddd-43ab-b543-c849de4683c4",
   "metadata": {
    "id": "99d7c0b4-0ddd-43ab-b543-c849de4683c4"
   },
   "source": [
    "# Network of LIF neurons\n",
    "**Credit: The original implementation of this example was created by Aditya Gilra in 2014 for the CAMP summer school at NCBS, Bangalore** based on `Ostojic, S. (2014). Two types of asynchronous activity in networks of excitatory and inhibitory spiking neurons. Nat Neurosci 17, 594-600.`\n",
    "\n",
    "*Look up `moose-examples` on `github` for his original object-oriented implementation. Here we shall implement this model in a simple procedural way.*\n",
    "\n",
    "Key parameters to play with synaptic coupling `J`, external input rate `nu_ext`, strength of inhibitory input relative to excitatory `ginh`.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "001cbf66-30eb-484b-b024-ada1f8a64cfb",
   "metadata": {
    "id": "001cbf66-30eb-484b-b024-ada1f8a64cfb"
   },
   "source": [
    "Rewritten by Subhasis Ray based on the article\n",
    "```\n",
    "N. Brunel, “Dynamics of Sparsely Connected Networks of Excitatory and Inhibitory Spiking Neurons,” J Comput Neurosci, vol. 8, no. 3, pp. 183–208, May 2000, doi: 10.1023/A:1008925309027.\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4042708-7255-4afe-9ea5-4653899d07e1",
   "metadata": {
    "id": "c4042708-7255-4afe-9ea5-4653899d07e1"
   },
   "outputs": [],
   "source": [
    "# Imports\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import moose"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1844f4c6-c6eb-4538-87de-ec0b7e1e453f",
   "metadata": {
    "id": "1844f4c6-c6eb-4538-87de-ec0b7e1e453f"
   },
   "outputs": [],
   "source": [
    "# Initialize numpy random number generator\n",
    "rng = np.random.default_rng(1)   # keep a constant seed for replicability\n",
    "\n",
    "N = 10\n",
    "N_E = int(0.8 * N)   # number of excitatory neurons\n",
    "N_I = N - N_E   # number of inhibitory neurons\n",
    "\n",
    "theta = 20e-3   # spiking threshold for all neurons\n",
    "Vr = 10e-3   # reset voltage 10 mV"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "608b9d2b-d1d4-4a26-8057-416645b16690",
   "metadata": {
    "id": "608b9d2b-d1d4-4a26-8057-416645b16690"
   },
   "source": [
    "## Model A\n",
    "Brunel calls the case where inhibitory and excitatory neurons have identical properties *model A*. Then considers the case where the two differ, and calls it *model B*.\n",
    "\n",
    "\n",
    "The model consists of $N_{E}$ excitatory neurons and $N_{I}$ inhibitory neurons.\n",
    "\n",
    "Each neuron receives $C$ random connections from other neurons.\n",
    "\n",
    "The number of connections from each type of neurons is proportional to the population size.\n",
    "\n",
    "So, the number of excitatory connections is $C_{E} = \\epsilon \\times N_{E}$\n",
    "\n",
    "and the number of inhibitory connections is $C_{I} = \\epsilon \\times N_{I}$.\n",
    "\n",
    "In addition, each neurons receives $C_{ext}$ excitatory connections from outside the network. Each of these external synapses transmit events from independent Poisson processes with rate $\\nu_{ext}$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7b3a9cb-5ddf-45f5-bea4-c91812cc9e18",
   "metadata": {
    "id": "f7b3a9cb-5ddf-45f5-bea4-c91812cc9e18"
   },
   "outputs": [],
   "source": [
    "# while J > 0 is PSP amplitude for all excitatory synapses,\n",
    "# PSP amplitude for inhibitory synapses is -gJ. We name g `ginh`\n",
    "J = 1e-3   # efficacy of excitatory synapse\n",
    "ginh = 1.0   # scale factor for inhibitory synapses\n",
    "# \"External synapses are activated by independent Poisson processes with rate nu_ext\"\n",
    "nu_ext = 100.0\n",
    "EPSILON = 1.0 # 0.1   # fraction of inputs from each type\n",
    "C_E = int(EPSILON * N_E)      # number of excitatory connections per neuron\n",
    "C_I = int(EPSILON * N_I)      # number of inhibitory connections per neuron\n",
    "\n",
    "C_EXT = C_E    # number of external excitatory connections per neuron\n",
    "\n",
    "DELAY = 1e-3    # synaptic delay\n",
    "tau_E = 20e-3   # time constant of excitatory neurons\n",
    "tau_rp = 2e-3   # refractory period\n",
    "\n",
    "Rm = tau_E # 20e6       # membrane resistance\n",
    "Cm = 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46a187e2-c293-473a-9676-8c781a9c4f1c",
   "metadata": {
    "id": "46a187e2-c293-473a-9676-8c781a9c4f1c"
   },
   "outputs": [],
   "source": [
    "# containers\n",
    "if moose.exists('/sim'):\n",
    "    moose.delete('/sim')\n",
    "\n",
    "sim = moose.Neutral('/sim')\n",
    "model = moose.Neutral(f'{sim.path}/model')\n",
    "data = moose.Neutral(f'{sim.path}/data')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e00e23df-2001-4aa7-aa26-9f405ad51055",
   "metadata": {
    "id": "e00e23df-2001-4aa7-aa26-9f405ad51055"
   },
   "source": [
    "Create an array of leaky-integrate fire neurons with `N` entries. All moose objects can be arrays if you specify the size as a second argument to the constructor.\n",
    "\n",
    "This will make `netA` an array / vector of LIF elements. While for single elements you access object attributes using dot notation like `{Object}.{field}`, to access the fields of a vector simultaneously, you have to do this via the `vec` attribute `{Object}.vec.{field}`.\n",
    "\n",
    "So, while for a single `LIF` object you would set the threshold as `LIF.thresh = value`, for the vector object it should be  `LIF.vec.thresh = value`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "434c6b94-55b0-4740-925e-14cb7e6cc22f",
   "metadata": {
    "id": "434c6b94-55b0-4740-925e-14cb7e6cc22f"
   },
   "outputs": [],
   "source": [
    "netA = moose.LIF(f'{model.path}/networkA', N)\n",
    "netA.vec.thresh = theta\n",
    "\n",
    "# Time constant of a parallel RC circuit is tau = Rm * Cm\n",
    "netA.vec.Cm = Cm\n",
    "netA.vec.Rm = Rm\n",
    "\n",
    "netA.vec.refractoryPeriod = tau_rp\n",
    "\n",
    "# set the voltage to reset potential when the similation is initialized with `moose.reinit()`.\n",
    "netA.vec.initVm = 0.0\n",
    "netA.vec.Em = 0.0\n",
    "netA.vec.vReset = Vr\n",
    "netA.vec.thresh = theta"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "AYG_4doiVxNX",
   "metadata": {
    "id": "AYG_4doiVxNX"
   },
   "source": [
    "Now we need to set up external inputs to the neurons. Each neuron receives $C_{ext} = C_{E}$ input connections from independent Poisson processes, each with $\\nu_{ext}$ rate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ftXovPiJYl3b",
   "metadata": {
    "id": "ftXovPiJYl3b"
   },
   "outputs": [],
   "source": [
    "ext_input = moose.RandSpike(f'{model.path}/ext_input', C_EXT * N)\n",
    "\n",
    "# The RandSpike object can generate random spikes with a given rate\n",
    "# approximating a Poisson process\n",
    "ext_input.vec.rate = nu_ext"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "EfaX-Tjf0cWt",
   "metadata": {
    "id": "EfaX-Tjf0cWt"
   },
   "source": [
    "Finally, we shall keep track of the spikes in the network using some a vector of `Table` objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "vwLtB5sm0xKo",
   "metadata": {
    "id": "vwLtB5sm0xKo"
   },
   "outputs": [],
   "source": [
    "ext_input_tab = moose.Table(f'{data.path}/ext_input_tab', ext_input.numData)\n",
    "moose.connect(ext_input, 'spikeOut', ext_input_tab, 'input', 'OneToOne')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "I6P4QR712GtT",
   "metadata": {
    "id": "I6P4QR712GtT"
   },
   "source": [
    "Let us break down the above. First, we are creating a `Table` array with the same number of elements as the `RandSpike` object `ext_input`. The `numData` field stores the number of elements in a MOOSE array object. A single object just has `numData = 1`.\n",
    "\n",
    "The `pymoose` interface also exposes the underlying vector as the attribute `vec` which is like a `sequence` (string, list, tuple, etc. are of the abstract data type sequence in Python, they all share some common properties). So, apart from `numData` attribute, you can also find the number of elements with `len(ext_input.vec)`.\n",
    "\n",
    "Then we are connecting `spikeOut` field of `ext_input` to the table's `input` field with a `OneToOne` message. `spikeOut` is a source field of `RandSpike` and  `input` is a destination field of `Table`. A `OneToOne` message connects these fields of corresponding elements in two vector objects, i.e., the `spikeOut` output value from the first element in `ext_input` will go into the `input` field of the first element of `ext_input_tab`, the output from the second element of `ext_input` will go into the `input` of the second element of `ext_input_tab`, and so on."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3q3pd94Z-Ihl",
   "metadata": {
    "id": "3q3pd94Z-Ihl"
   },
   "source": [
    "\n",
    "Next, we collect the spikes from all the elements in `netA` into another `Table` array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "MSlTL7jH4vof",
   "metadata": {
    "id": "MSlTL7jH4vof"
   },
   "outputs": [],
   "source": [
    "spike_tab = moose.Table(f'{data.path}/spike_tab', netA.numData)\n",
    "moose.connect(netA, 'spikeOut', spike_tab, 'input', 'OneToOne')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "jiSv-H3B5RnS",
   "metadata": {
    "id": "jiSv-H3B5RnS"
   },
   "source": [
    "Also, we want to monitor the state variable `Vm` of a few neurons. While spikes are discrete events in time, `Vm` is a continuous variable and is updated at every timestep of the simulation. However, the changes are so small from timestep to timestep, that we can sample its value at a much slower rate than it is computed. So instead of pushing the value at every timestep of `Vm` computation, we can get the `Table` to pull the value of `Vm` from the `LIF` object at a much slower rate. This will reduce both computer time and the amount of redundant data storage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "JqhM_I0W7KGY",
   "metadata": {
    "id": "JqhM_I0W7KGY"
   },
   "outputs": [],
   "source": [
    "vm_tab = moose.Table(f'{data.path}/vm_tab', 10)   # record from just 10 random LIF elements\n",
    "\n",
    "# Select random neurons from netA for recording Vm\n",
    "vm_idx = rng.choice(netA.numData, replace=False, size=vm_tab.numData)\n",
    "\n",
    "for ii, jj in enumerate(vm_idx):\n",
    "  moose.connect(vm_tab.vec[ii], 'requestOut', netA.vec[jj], 'getVm')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "kAdNGCGd0SbS",
   "metadata": {
    "id": "kAdNGCGd0SbS"
   },
   "source": [
    "Now we shall setup the connections in the network itself. First, we must create synapses on the neurons. We use a `SimpleSynHandler` which collects spike events arriving from all the presynaptic neurons in each timestep and sends out the combined value through `activationOut` field. This is passed on to the postsynaptic `LIF` neuron. Each event increments the `Vm` of the postsynaptic neuron by the `weight` value of the synapse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "O44T0OMF8yrp",
   "metadata": {
    "id": "O44T0OMF8yrp"
   },
   "outputs": [],
   "source": [
    "synh = moose.SimpleSynHandler(f'{model.path}/synh', netA.numData)\n",
    "moose.connect(synh, 'activationOut', netA, 'activation', 'OneToOne')\n",
    "# Each neuron gets\n",
    "# C_E synaptic inputs from excitatory neurons\n",
    "# C_I synaptic inputs from inhibitory neurons\n",
    "# and C_EXT = C_E sexternal excitatory synaptic inputs\n",
    "synh.vec.numSynapses = C_E + C_I + C_EXT\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "UVVxXhz-ARFF",
   "metadata": {
    "id": "UVVxXhz-ARFF"
   },
   "source": [
    "Out of the `N` elements in the `LIF` vector, we shall designate the first `N_E` as excitatory, and the rest `N_I` as inhibitory neurons. Whether a neuron is excitatory or inhibitory depends only on the synapse.\n",
    "\n",
    "- The first `C_E` connections to each synapse will be from randomly selected neurons out of the first `N_E`.\n",
    "- The next `C_I` connections will be from the next `N_I` inhibitory neurons.\n",
    "- Finally, the last `C_EXT` connections will come from the corresponding external Poisson spike generators in `ext_input`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "nyIs6yWiArBP",
   "metadata": {
    "id": "nyIs6yWiArBP"
   },
   "outputs": [],
   "source": [
    "for post_idx in range(N):\n",
    "  # Elements 0 to N_E of netA are the neurons we designated excitatory\n",
    "  # The first C_E synapses are excitatory, the rest C_I are inhibitory\n",
    "  syn_idx = 0\n",
    "\n",
    "  print(f'Connecting {syn_idx} to {syn_idx + C_E} excitatory inputs')\n",
    "  pre_indices = set(range(N_E))\n",
    "  # Avoid self connection.\n",
    "  # Unlike set.remove(), set.discard() does not raise an error if the element\n",
    "  # is not in the set.\n",
    "  pre_indices.discard(post_idx)\n",
    "  for pre_idx in rng.choice(list(pre_indices), size=C_E):\n",
    "    synapse = synh.vec[post_idx].synapse[syn_idx]\n",
    "    synapse.weight = J\n",
    "    synapse.delay = DELAY\n",
    "    moose.connect(netA.vec[pre_idx], 'spikeOut', synapse, 'addSpike')\n",
    "    syn_idx += 1\n",
    "  print(f'Connecting {syn_idx} to {syn_idx + C_I} inhibitory inputs')\n",
    "  # Elements N_E to N of netA are the N_I neurons we designated inhibitory\n",
    "  pre_indices = set(range(N_E, N))\n",
    "  pre_indices.discard(post_idx)  # avoid self connection\n",
    "  for pre_idx in rng.choice(list(pre_indices), size=C_I):\n",
    "    # Synapses C_E till the last are the C_I inhibitory synapses\n",
    "    synapse = synh.vec[post_idx].synapse[syn_idx]\n",
    "    synapse.weight = - ginh * J\n",
    "    synapse.delay = DELAY\n",
    "    moose.connect(netA.vec[pre_idx], 'spikeOut', synapse, 'addSpike')\n",
    "    syn_idx += 1\n",
    "\n",
    "  # Connect external inputs\n",
    "  print(f'Connecting {syn_idx} to {syn_idx + C_EXT} external inputs'\n",
    "  f' from {post_idx * C_EXT} to {(post_idx + 1) * C_EXT}')\n",
    "  for pre_idx in range(post_idx * C_EXT, (post_idx + 1) * C_EXT):\n",
    "    synapse = synh.vec[post_idx].synapse[syn_idx]\n",
    "    synapse.weight = J\n",
    "    moose.connect(ext_input.vec[pre_idx], 'spikeOut', synapse, 'addSpike')\n",
    "    syn_idx += 1\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "BCVsZ7gX1Qif",
   "metadata": {
    "id": "BCVsZ7gX1Qif"
   },
   "source": [
    "To check that model is correctly wired up, we can look at the connections using `moose.showmsg()` function.\n",
    "Comment these out when you have large number of neurons because it may stall the system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "JORukKBv1e7w",
   "metadata": {
    "id": "JORukKBv1e7w"
   },
   "outputs": [],
   "source": [
    "moose.showmsg(netA.vec[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03jJyDDxWo-C",
   "metadata": {
    "id": "03jJyDDxWo-C"
   },
   "outputs": [],
   "source": [
    "moose.showmsg(synh.vec[0].synapse[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "XGMJc1LIFWcg",
   "metadata": {
    "id": "XGMJc1LIFWcg"
   },
   "source": [
    "Now we can simulate the model after initializing it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "R0IchddmFZaG",
   "metadata": {
    "id": "R0IchddmFZaG"
   },
   "outputs": [],
   "source": [
    "#dt = 1e-5\n",
    "runtime = 1.0e-1\n",
    "# for ii in range(10):\n",
    "#   moose.setClock( ii, dt )\n",
    "moose.reinit()\n",
    "moose.start(runtime)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "Q7KBKifwIj5t",
   "metadata": {
    "id": "Q7KBKifwIj5t"
   },
   "source": [
    "First, for sanity check, we should see if the external spikes are being generated. For this we create a raster plot of the spike times of 1000 randomly selected external input source. Each row in the raster plot represents one external spike train, and the horizontal axis represents time. For each spike event in a spike train, we put a mark in its row at the time of the event."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9Z8KLIMxIr2l",
   "metadata": {
    "id": "9Z8KLIMxIr2l"
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "for ii, ext_idx in enumerate(rng.choice(ext_input.numData,\n",
    "                                        size=min(1000, ext_input.numData),\n",
    "                                        replace=False)):\n",
    "  spike_times = ext_input_tab.vec[ext_idx].vector\n",
    "  ax.plot(spike_times, [ii] * len(spike_times), 'k|')\n",
    "ax.set_xlabel('Time (s)')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "FEHhgWJJFie3",
   "metadata": {
    "id": "FEHhgWJJFie3"
   },
   "source": [
    "We can plot the `Vm` serieses we recorded."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "uGLA8a8_Hf_L",
   "metadata": {
    "id": "uGLA8a8_Hf_L"
   },
   "outputs": [],
   "source": [
    "vm_tab.vec[0].vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fEcNsmCUFmX_",
   "metadata": {
    "id": "fEcNsmCUFmX_"
   },
   "outputs": [],
   "source": [
    "t = np.arange(0, runtime + vm_tab.dt/2, vm_tab.dt)\n",
    "fig, ax = plt.subplots()\n",
    "for tab in vm_tab.vec:\n",
    "  ax.plot(t, tab.vector, label=tab.name)\n",
    "  # break\n",
    "ax.set_xlabel('Time (s)')\n",
    "ax.set_ylabel('Vm (mV)')\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "TP71wIzF6wyY",
   "metadata": {
    "id": "TP71wIzF6wyY"
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "for ii, idx in enumerate(rng.choice(netA.numData,\n",
    "                                        size=min(100, netA.numData),\n",
    "                                        replace=False)):\n",
    "  spike_times = spike_tab.vec[idx].vector\n",
    "  ax.plot(spike_times, [ii] * len(spike_times), 'k|')\n",
    "ax.set_xlabel('Time (s)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "XSn-1kukZeOj",
   "metadata": {
    "id": "XSn-1kukZeOj"
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "include_colab_link": true,
   "provenance": []
  },
  "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}