{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Random balanced network (exp synapses, multiple time constants)\n\nThis script simulates an excitatory and an inhibitory population on\nthe basis of the network used in [1]_.\n\nThe example demonstrate the usage of the multisynapse neuron\nmodel. Each spike arriving at the neuron triggers an exponential\nPSP. The time constant associated with the PSP is defined in the\nreceptor type array tau_syn of each neuron. The receptor types of all\nconnections are uniformally distributed, resulting in uniformally\ndistributed time constants of the PSPs.\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_nest`\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Import all necessary modules for simulation, analysis and plotting.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import time\nimport nest\nimport nest.raster_plot\nimport matplotlib.pyplot as plt\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.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "tauMem = 20.0  # time constant of membrane potential in ms\ntheta = 20.0  # membrane threshold potential in mV\nJ = 0.1   # postsynaptic amplitude in mV\nnr_ports = 100  # number of receptor types\n# Create array of synaptic time constants for each neuron,\n# ranging from 0.1 to 1.09 ms.\ntau_syn = [0.1 + 0.01 * i for i in range(nr_ports)]\nneuron_params = {\"C_m\": 1.0,\n                 \"tau_m\": tauMem,\n                 \"t_ref\": 2.0,\n                 \"E_L\": 0.0,\n                 \"V_reset\": 0.0,\n                 \"V_m\": 0.0,\n                 \"V_th\": theta,\n                 \"tau_syn\": tau_syn}\nJ_ex = J       # 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 / (J * CE * tauMem)\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.SetKernelStatus({\"resolution\": dt, \"print_time\": True,\n                      \"overwrite_files\": True, 'local_num_threads': 4})\n\nprint(\"Building network\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Configuration of the model ``iaf_psc_exp_multisynapse`` and\n``poisson_generator`` using ``SetDefaults``. This function expects the model to\nbe the inserted as a string and the parameter to be specified in a\ndictionary. All instances of theses models created after this point will\nhave the properties specified in the dictionary by default.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.SetDefaults(\"iaf_psc_exp_multisynapse\", neuron_params)\nnest.SetDefaults(\"poisson_generator\", {\"rate\": p_rate})"
      ]
    },
    {
      "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.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nodes_ex = nest.Create(\"iaf_psc_exp_multisynapse\", NE)\nnodes_in = nest.Create(\"iaf_psc_exp_multisynapse\", NI)\nnoise = nest.Create(\"poisson_generator\")\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\npre-defined # in a dictionary and inserted via ``syn_spec``. As synaptic model\nthe pre-defined synapses \"excitatory\" and \"inhibitory\" are choosen,\nthus setting ``weight`` and ``delay``. The receptor type is drawn from a\ndistribution for each connection, which is specified in the synapse\nproperties by assigning a dictionary to the keyword ``receptor_type``,\nwhich includes the specification of the distribution and the associated\nparameter.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "syn_params_ex = {\"model\": \"excitatory\",\n                 \"receptor_type\": {\"distribution\": \"uniform_int\",\n                                   \"low\": 1, \"high\": nr_ports}\n                 }\nsyn_params_in = {\"model\": \"inhibitory\",\n                 \"receptor_type\": {\"distribution\": \"uniform_int\",\n                                   \"low\": 1, \"high\": nr_ports}\n                 }\n\nnest.Connect(noise, nodes_ex, syn_spec=syn_params_ex)\nnest.Connect(noise, nodes_in, syn_spec=syn_params_ex)"
      ]
    },
    {
      "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 while distribution the\nports. Here we use the previously defined parameter dictionary\n``syn_params_ex``. Beforehand, the connection parameter are defined in a\ndictionary. Here we use the connection rule ``fixed_indegree``,\nwhich requires the definition of the indegree.\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, syn_params_ex)\n\nprint(\"Inhibitory connections\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Connecting the inhibitory population to all neurons while distribution the\nports. Here we use the previously defined parameter dictionary\n``syn_params_in``.The connection parameter are defined analogously to the\nconnection from the excitatory population 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, syn_params_in)"
      ]
    },
    {
      "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(\"Number of neurons : {0}\".format(N_neurons))\nprint(\"Number of synapses: {0}\".format(num_synapses))\nprint(\"       Exitatory  : {0}\".format(int(CE * N_neurons) + N_neurons))\nprint(\"       Inhibitory : {0}\".format(int(CI * N_neurons)))\nprint(\"Excitatory rate   : %.2f Hz\" % rate_ex)\nprint(\"Inhibitory rate   : %.2f Hz\" % rate_in)\nprint(\"Building time     : %.2f s\" % build_time)\nprint(\"Simulation time   : %.2f s\" % sim_time)"
      ]
    },
    {
      "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.6.12"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}