{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Current-based generalized leaky integrate and fire (GLIF) neuron example\n\nSimple example of how to use the ``glif_psc`` neuron model for\nfive different levels of GLIF neurons.\n\nFour stimulation paradigms are illustrated for the GLIF model\nwith externally applied current and spikes impinging\n\nVoltage traces, current traces, threshold traces, and spikes are shown.\n\nKEYWORDS: glif_psc\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First, we import all necessary modules to simulate, analyze and plot this\nexample.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import nest\nimport matplotlib.pyplot as plt\nimport matplotlib.gridspec as gridspec"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We initialize NEST and set the simulation resolution.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.ResetKernel()\nresolution = 0.05\nnest.resolution = resolution"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We also pre-define the synapse time constant array, [2.0, 1.0] ms for\nthe two desired synaptic ports of the GLIF neurons. Note that the default\nsynapse time constant is [2.0] ms, which is for neuron with one port.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "syn_tau = [2.0, 1.0]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We create the five levels of GLIF model to be tested, i.e.,\n``lif``, ``lif_r``, ``lif_asc``, ``lif_r_asc``, ``lif_r_asc_a``.\nFor each level of GLIF model, we create  a ``glif_psc`` node. The node is\ncreated by setting relative model mechanism parameters and the time constant\nof the 2 synaptic ports as mentioned above. Other neuron parameters are set\nas default. The five ``glif_psc`` node handles were combined as a list.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "n_lif = nest.Create(\"glif_psc\",\n                    params={\"spike_dependent_threshold\": False,\n                            \"after_spike_currents\": False,\n                            \"adapting_threshold\": False,\n                            \"tau_syn\": syn_tau})\nn_lif_r = nest.Create(\"glif_psc\",\n                      params={\"spike_dependent_threshold\": True,\n                              \"after_spike_currents\": False,\n                              \"adapting_threshold\": False,\n                              \"tau_syn\": syn_tau})\nn_lif_asc = nest.Create(\"glif_psc\",\n                        params={\"spike_dependent_threshold\": False,\n                                \"after_spike_currents\": True,\n                                \"adapting_threshold\": False,\n                                \"tau_syn\": syn_tau})\nn_lif_r_asc = nest.Create(\"glif_psc\",\n                          params={\"spike_dependent_threshold\": True,\n                                  \"after_spike_currents\": True,\n                                  \"adapting_threshold\": False,\n                                  \"tau_syn\": syn_tau})\nn_lif_r_asc_a = nest.Create(\"glif_psc\",\n                            params={\"spike_dependent_threshold\": True,\n                                    \"after_spike_currents\": True,\n                                    \"adapting_threshold\": True,\n                                    \"tau_syn\": syn_tau})\n\nneurons = n_lif + n_lif_r + n_lif_asc + n_lif_r_asc + n_lif_r_asc_a"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For the stimulation input to the glif_psc neurons, we create one excitation\nspike generator and one inhibition spike generator, each of which generates\nthree spikes; we also create one step current generator and a Poisson\ngenerator, a parrot neuron (to be paired with the Poisson generator).\nThe three different injections are spread to three different time periods,\ni.e., 0 ms ~ 200 ms, 200 ms ~ 500 ms, 600 ms ~ 900 ms.\nEach of the excitation and inhibition spike generators generates three spikes\nat different time points. Configuration of the current generator includes the\ndefinition of the start and stop times and the amplitude of the injected\ncurrent. Configuration of the Poisson generator includes the definition of\nthe start and stop times and the rate of the injected spike train.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "espikes = nest.Create(\"spike_generator\",\n                      params={\"spike_times\": [10., 100., 150.],\n                              \"spike_weights\": [20.] * 3})\nispikes = nest.Create(\"spike_generator\",\n                      params={\"spike_times\": [15., 99., 150.],\n                              \"spike_weights\": [-20.] * 3})\ncg = nest.Create(\"step_current_generator\",\n                 params={\"amplitude_values\": [400., ],\n                         \"amplitude_times\": [200., ],\n                         \"start\": 200., \"stop\": 500.})\npg = nest.Create(\"poisson_generator\",\n                 params={\"rate\": 150000., \"start\": 600., \"stop\": 900.})\npn = nest.Create(\"parrot_neuron\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The generators are then connected to the neurons. Specification of\nthe ``receptor_type`` uniquely defines the target receptor.\nWe connect current generator, the spike generators, Poisson generator (via\nparrot neuron) to receptor 0, 1, and 2 of the GLIF neurons, respectively.\nNote that Poisson generator is connected to parrot neuron to transit the\nspikes to the glif_psc neuron.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.Connect(cg, neurons, syn_spec={\"delay\": resolution})\nnest.Connect(espikes, neurons,\n             syn_spec={\"delay\": resolution, \"receptor_type\": 1})\nnest.Connect(ispikes, neurons,\n             syn_spec={\"delay\": resolution, \"receptor_type\": 1})\nnest.Connect(pg, pn, syn_spec={\"delay\": resolution})\nnest.Connect(pn, neurons, syn_spec={\"delay\": resolution, \"receptor_type\": 2})"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "A ``multimeter`` is created and connected to the neurons. The parameters\nspecified for the multimeter include the list of quantities that should be\nrecorded and the time interval at which quantities are measured.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mm = nest.Create(\"multimeter\",\n                 params={\"interval\": resolution,\n                         \"record_from\": [\"V_m\", \"I\", \"I_syn\", \"threshold\",\n                                         \"threshold_spike\",\n                                         \"threshold_voltage\",\n                                         \"ASCurrents_sum\"]})\nnest.Connect(mm, neurons)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "A ``spike_recorder`` is created and connected to the neurons record the\nspikes generated by the glif_psc neurons.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sr = nest.Create(\"spike_recorder\")\nnest.Connect(neurons, sr)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Run the simulation for 1000 ms and retrieve recorded data from\nthe multimeter and spike recorder.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.Simulate(1000.)\n\ndata = mm.events\nsenders = data[\"senders\"]\n\nspike_data = sr.events\nspike_senders = spike_data[\"senders\"]\nspikes = spike_data[\"times\"]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We plot the time traces of the membrane potential (in blue) and\nthe overall threshold (in green), and the spikes (as red dots) in one panel;\nthe spike component of threshold (in yellow) and the voltage component of\nthreshold (in black) in another panel; the injected currents (in strong blue),\nthe sum of after spike currents (in cyan), and the synaptic currents (in\nmagenta) in responding to the spike inputs to the neurons in the third panel.\nWe plot all these three panels for each level of GLIF model in a separated\nfigure.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "glif_models = [\"lif\", \"lif_r\", \"lif_asc\", \"lif_r_asc\", \"lif_r_asc_a\"]\nfor i in range(len(glif_models)):\n\n    glif_model = glif_models[i]\n    node_id = neurons[i].global_id\n    plt.figure(glif_model)\n    gs = gridspec.GridSpec(3, 1, height_ratios=[2, 1, 1])\n    t = data[\"times\"][senders == 1]\n\n    ax1 = plt.subplot(gs[0])\n    plt.plot(t, data[\"V_m\"][senders == node_id], \"b\")\n    plt.plot(t, data[\"threshold\"][senders == node_id], \"g--\")\n    plt.plot(spikes[spike_senders == node_id],\n             [max(data[\"threshold\"][senders == node_id]) * 0.95] *\n             len(spikes[spike_senders == node_id]), \"r.\")\n    plt.legend([\"V_m\", \"threshold\", \"spike\"])\n    plt.ylabel(\"V (mV)\")\n    plt.title(\"Simulation of glif_psc neuron of \" + glif_model)\n\n    ax2 = plt.subplot(gs[1])\n    plt.plot(t, data[\"threshold_spike\"][senders == node_id], \"y\")\n    plt.plot(t, data[\"threshold_voltage\"][senders == node_id], \"k--\")\n    plt.legend([\"threshold_spike\", \"threshold_voltage\"])\n    plt.ylabel(\"V (mV)\")\n\n    ax3 = plt.subplot(gs[2])\n    plt.plot(t, data[\"I\"][senders == node_id], \"--\")\n    plt.plot(t, data[\"ASCurrents_sum\"][senders == node_id], \"c-.\")\n    plt.plot(t, data[\"I_syn\"][senders == node_id], \"m\")\n    plt.legend([\"I_e\", \"ASCurrents_sum\", \"I_syn\"])\n    plt.ylabel(\"I (pA)\")\n    plt.xlabel(\"t (ms)\")\n\nplt.show()"
      ]
    }
  ],
  "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.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}