{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Random balanced network (alpha synapses) connected with NEST\n\nThis script simulates an excitatory and an inhibitory population on\nthe basis of the network used in [1]_.\n\nIn contrast to ``brunel-alpha-numpy.py``, this variant uses NEST's builtin\nconnection routines to draw the random connections instead of NumPy.\n\nWhen connecting the network, customary synapse models are used, which\nallow for querying the number of created synapses. Using spike\nrecorders, the average firing rates of the neurons in the populations\nare established. The building as well as the simulation time of the\nnetwork are recorded.\n\n## References\n\n.. [1] Brunel N (2000). Dynamics of sparsely connected networks of excitatory and\n       inhibitory spiking neurons. Journal of Computational Neuroscience 8,\n       183-208.\n\n## See Also\n\n:doc:`brunel_alpha_numpy`\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Import all necessary modules for simulation, analysis and plotting. Scipy\nshould be imported before nest.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import time\nimport numpy as np\nimport scipy.special as sp\n\nimport nest\nimport nest.raster_plot\nimport matplotlib.pyplot as plt"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Definition of functions used in this example. First, define the `Lambert W`\nfunction implemented in SLI. The second function computes the maximum of\nthe postsynaptic potential for a synaptic input current of unit amplitude\n(1 pA) using the `Lambert W` function. Thus function will later be used to\ncalibrate the synaptic weights.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def LambertWm1(x):\n    # Using scipy to mimic the gsl_sf_lambert_Wm1 function.\n    return sp.lambertw(x, k=-1 if x < 0 else 0).real\n\n\ndef ComputePSPnorm(tauMem, CMem, tauSyn):\n    a = (tauMem / tauSyn)\n    b = (1.0 / tauSyn - 1.0 / tauMem)\n\n    # time of maximum\n    t_max = 1.0 / b * (-LambertWm1(-np.exp(-1.0 / a) / a) - 1.0 / a)\n\n    # maximum of PSP for current of unit amplitude\n    return (np.exp(1.0) / (tauSyn * CMem * b) *\n            ((np.exp(-t_max / tauMem) - np.exp(-t_max / tauSyn)) / b -\n             t_max * np.exp(-t_max / tauSyn)))\n\n\nnest.ResetKernel()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Assigning the current time to a variable in order to determine the build\ntime of the network.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "startbuild = time.time()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Assigning the simulation parameters to variables.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dt = 0.1    # the resolution in ms\nsimtime = 1000.0  # Simulation time in ms\ndelay = 1.5    # synaptic delay in ms"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Definition of the parameters crucial for asynchronous irregular firing of\nthe neurons.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "g = 5.0  # ratio inhibitory weight/excitatory weight\neta = 2.0  # external rate relative to threshold rate\nepsilon = 0.1  # connection probability"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Definition of the number of neurons in the network and the number of neurons\nrecorded from\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "order = 2500\nNE = 4 * order  # number of excitatory neurons\nNI = 1 * order  # number of inhibitory neurons\nN_neurons = NE + NI   # number of neurons in total\nN_rec = 50      # record from 50 neurons"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Definition of connectivity parameters\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "CE = int(epsilon * NE)  # number of excitatory synapses per neuron\nCI = int(epsilon * NI)  # number of inhibitory synapses per neuron\nC_tot = int(CI + CE)      # total number of synapses per neuron"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Initialization of the parameters of the integrate and fire neuron and the\nsynapses. The parameters of the neuron are stored in a dictionary. The\nsynaptic currents are normalized such that the amplitude of the PSP is J.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "tauSyn = 0.5  # synaptic time constant in ms\ntauMem = 20.0  # time constant of membrane potential in ms\nCMem = 250.0  # capacitance of membrane in in pF\ntheta = 20.0  # membrane threshold potential in mV\nneuron_params = {\"C_m\": CMem,\n                 \"tau_m\": tauMem,\n                 \"tau_syn_ex\": tauSyn,\n                 \"tau_syn_in\": tauSyn,\n                 \"t_ref\": 2.0,\n                 \"E_L\": 0.0,\n                 \"V_reset\": 0.0,\n                 \"V_m\": 0.0,\n                 \"V_th\": theta}\nJ = 0.1        # postsynaptic amplitude in mV\nJ_unit = ComputePSPnorm(tauMem, CMem, tauSyn)\nJ_ex = J / J_unit  # amplitude of excitatory postsynaptic current\nJ_in = -g * J_ex    # amplitude of inhibitory postsynaptic current"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Definition of threshold rate, which is the external rate needed to fix the\nmembrane potential around its threshold, the external firing rate and the\nrate of the poisson generator which is multiplied by the in-degree CE and\nconverted to Hz by multiplication by 1000.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nu_th = (theta * CMem) / (J_ex * CE * np.exp(1) * tauMem * tauSyn)\nnu_ex = eta * nu_th\np_rate = 1000.0 * nu_ex * CE"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Configuration of the simulation kernel by the previously defined time\nresolution used in the simulation. Setting ``print_time`` to `True` prints the\nalready processed simulation time as well as its percentage of the total\nsimulation time.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.resolution = dt\nnest.print_time = True\nnest.overwrite_files = True\n\nprint(\"Building network\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Creation of the nodes using ``Create``. We store the returned handles in\nvariables for later reference. Here the excitatory and inhibitory, as well\nas the poisson generator and two spike recorders. The spike recorders will\nlater be used to record excitatory and inhibitory spikes. Properties of the\nnodes are specified via ``params``, which expects a dictionary.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nodes_ex = nest.Create(\"iaf_psc_alpha\", NE, params=neuron_params)\nnodes_in = nest.Create(\"iaf_psc_alpha\", NI, params=neuron_params)\nnoise = nest.Create(\"poisson_generator\", params={\"rate\": p_rate})\nespikes = nest.Create(\"spike_recorder\")\nispikes = nest.Create(\"spike_recorder\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Configuration of the spike recorders recording excitatory and inhibitory\nspikes by sending parameter dictionaries to ``set``. Setting the property\n`record_to` to *\"ascii\"* ensures that the spikes will be recorded to a file,\nwhose name starts with the string assigned to the property `label`.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "espikes.set(label=\"brunel-py-ex\", record_to=\"ascii\")\nispikes.set(label=\"brunel-py-in\", record_to=\"ascii\")\n\nprint(\"Connecting devices\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Definition of a synapse using ``CopyModel``, which expects the model name of\na pre-defined synapse, the name of the customary synapse and an optional\nparameter dictionary. The parameters defined in the dictionary will be the\ndefault parameter for the customary synapse. Here we define one synapse for\nthe excitatory and one for the inhibitory connections giving the\npreviously defined weights and equal delays.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.CopyModel(\"static_synapse\", \"excitatory\",\n               {\"weight\": J_ex, \"delay\": delay})\nnest.CopyModel(\"static_synapse\", \"inhibitory\",\n               {\"weight\": J_in, \"delay\": delay})"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Connecting the previously defined poisson generator to the excitatory and\ninhibitory neurons using the excitatory synapse. Since the poisson\ngenerator is connected to all neurons in the population the default rule\n(``all_to_all``) of ``Connect`` is used. The synaptic properties are inserted\nvia ``syn_spec`` which expects a dictionary when defining multiple variables or\na string when simply using a pre-defined synapse.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.Connect(noise, nodes_ex, syn_spec=\"excitatory\")\nnest.Connect(noise, nodes_in, syn_spec=\"excitatory\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Connecting the first ``N_rec`` nodes of the excitatory and inhibitory\npopulation to the associated spike recorders using excitatory synapses.\nHere the same shortcut for the specification of the synapse as defined\nabove is used.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.Connect(nodes_ex[:N_rec], espikes, syn_spec=\"excitatory\")\nnest.Connect(nodes_in[:N_rec], ispikes, syn_spec=\"excitatory\")\n\nprint(\"Connecting network\")\n\nprint(\"Excitatory connections\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Connecting the excitatory population to all neurons using the pre-defined\nexcitatory synapse. Beforehand, the connection parameter are defined in a\ndictionary. Here we use the connection rule ``fixed_indegree``,\nwhich requires the definition of the indegree. Since the synapse\nspecification is reduced to assigning the pre-defined excitatory synapse it\nsuffices to insert a string.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "conn_params_ex = {'rule': 'fixed_indegree', 'indegree': CE}\nnest.Connect(nodes_ex, nodes_ex + nodes_in, conn_params_ex, \"excitatory\")\n\nprint(\"Inhibitory connections\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Connecting the inhibitory population to all neurons using the pre-defined\ninhibitory synapse. The connection parameter as well as the synapse\nparameter are defined analogously to the connection from the excitatory\npopulation defined above.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "conn_params_in = {'rule': 'fixed_indegree', 'indegree': CI}\nnest.Connect(nodes_in, nodes_ex + nodes_in, conn_params_in, \"inhibitory\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Storage of the time point after the buildup of the network in a variable.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "endbuild = time.time()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Simulation of the network.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(\"Simulating\")\n\nnest.Simulate(simtime)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Storage of the time point after the simulation of the network in a variable.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "endsimulate = time.time()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Reading out the total number of spikes received from the spike recorder\nconnected to the excitatory population and the inhibitory population.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "events_ex = espikes.n_events\nevents_in = ispikes.n_events"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Calculation of the average firing rate of the excitatory and the inhibitory\nneurons by dividing the total number of recorded spikes by the number of\nneurons recorded from and the simulation time. The multiplication by 1000.0\nconverts the unit 1/ms to 1/s=Hz.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rate_ex = events_ex / simtime * 1000.0 / N_rec\nrate_in = events_in / simtime * 1000.0 / N_rec"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Reading out the number of connections established using the excitatory and\ninhibitory synapse model. The numbers are summed up resulting in the total\nnumber of synapses.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "num_synapses = (nest.GetDefaults(\"excitatory\")[\"num_connections\"] +\n                nest.GetDefaults(\"inhibitory\")[\"num_connections\"])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Establishing the time it took to build and simulate the network by taking\nthe difference of the pre-defined time variables.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "build_time = endbuild - startbuild\nsim_time = endsimulate - endbuild"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Printing the network properties, firing rates and building times.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(\"Brunel network simulation (Python)\")\nprint(f\"Number of neurons : {N_neurons}\")\nprint(f\"Number of synapses: {num_synapses}\")\nprint(f\"       Exitatory  : {int(CE * N_neurons) + N_neurons}\")\nprint(f\"       Inhibitory : {int(CI * N_neurons)}\")\nprint(f\"Excitatory rate   : {rate_ex:.2f} Hz\")\nprint(f\"Inhibitory rate   : {rate_in:.2f} Hz\")\nprint(f\"Building time     : {build_time:.2f} s\")\nprint(f\"Simulation time   : {sim_time:.2f} s\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot a raster of the excitatory neurons and a histogram.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.raster_plot.from_device(espikes, hist=True)\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
}