{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Balanced neuron example\n\nThis script simulates a neuron driven by an excitatory and an\ninhibitory population of neurons firing Poisson spike trains. The aim\nis to find a firing rate for the inhibitory population that will make\nthe neuron fire at the same rate as the excitatory population.\n\nOptimization is performed using the ``bisection`` method from Scipy,\nsimulating the network repeatedly.\n\nThis example is also shown in the article [1]_\n\n## References\n\n.. [1] Eppler JM, Helias M, Mulller E, Diesmann M, Gewaltig MO (2009). PyNEST: A convenient interface to the NEST\n       simulator, Front. Neuroinform.\n       http://dx.doi.org/10.3389/neuro.11.012.2008\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First, we import all necessary modules for simulation, analysis and\nplotting. Scipy should be imported before nest.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from scipy.optimize import bisect\n\nimport nest\nimport nest.voltage_trace"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Additionally, we set the verbosity using ``set_verbosity`` to\nsuppress info messages.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.set_verbosity(\"M_WARNING\")\nnest.ResetKernel()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Second, the simulation parameters are assigned to variables.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "t_sim = 25000.0  # how long we simulate\nn_ex = 16000     # size of the excitatory population\nn_in = 4000      # size of the inhibitory population\nr_ex = 5.0       # mean rate of the excitatory population\nr_in = 20.5      # initial rate of the inhibitory population\nepsc = 45.0      # peak amplitude of excitatory synaptic currents\nipsc = -45.0     # peak amplitude of inhibitory synaptic currents\nd = 1.0          # synaptic delay\nlower = 15.0     # lower bound of the search interval\nupper = 25.0     # upper bound of the search interval\nprec = 0.01      # how close need the excitatory rates be"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Third, the nodes are created using ``Create``. We store the returned\nhandles in variables for later reference.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "neuron = nest.Create(\"iaf_psc_alpha\")\nnoise = nest.Create(\"poisson_generator\", 2)\nvoltmeter = nest.Create(\"voltmeter\")\nspikerecorder = nest.Create(\"spike_recorder\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Fourth, the ``poisson_generator`` (`noise`) is configured using ``set``.\nNote that we need not set parameters for the neuron, the spike recorder, and\nthe voltmeter, since they have satisfactory defaults.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "noise.rate = [n_ex * r_ex, n_in * r_in]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Fifth, the ``iaf_psc_alpha`` is connected to the ``spike_recorder`` and the\n``voltmeter``, as are the two Poisson generators to the neuron. The command\n``Connect`` has different variants. Plain `Connect` just takes the handles of\npre- and post-synaptic nodes and uses the default values for weight and\ndelay. It can also be called with a list of weights, as in the connection\nof the noise below.\nNote that the connection direction for the ``voltmeter`` is reversed compared\nto the ``spike_recorder``, because it observes the neuron instead of\nreceiving events from it. Thus, ``Connect`` reflects the direction of signal\nflow in the simulation kernel rather than the physical process of inserting\nan electrode into the neuron. The latter semantics is presently not\navailable in NEST.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.Connect(neuron, spikerecorder)\nnest.Connect(voltmeter, neuron)\nnest.Connect(noise, neuron, syn_spec={'weight': [[epsc, ipsc]], 'delay': 1.0})"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To determine the optimal rate of the neurons in the inhibitory population,\nthe network is simulated several times for different values of the\ninhibitory rate while measuring the rate of the target neuron. This is done\nby calling ``Simulate`` until the rate of the target neuron matches the rate\nof the neurons in the excitatory population with a certain accuracy. The\nalgorithm is implemented in two steps:\n\nFirst, the function ``output_rate`` is defined to measure the firing rate\nof the target neuron for a given rate of the inhibitory neurons.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def output_rate(guess):\n    print(\"Inhibitory rate estimate: %5.2f Hz\" % guess)\n    rate = float(abs(n_in * guess))\n    noise[1].rate = rate\n    spikerecorder.n_events = 0\n    nest.Simulate(t_sim)\n    out = spikerecorder.n_events * 1000.0 / t_sim\n    print(\"  -> Neuron rate: %6.2f Hz (goal: %4.2f Hz)\" % (out, r_ex))\n    return out"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The function takes the firing rate of the inhibitory neurons as an\nargument. It scales the rate with the size of the inhibitory population and\nconfigures the inhibitory Poisson generator (`noise[1]`) accordingly.\nThen, the spike counter of the ``spike_recorder`` is reset to zero. The\nnetwork is simulated using ``Simulate``, which takes the desired simulation\ntime in milliseconds and advances the network state by this amount of time.\nDuring simulation, the ``spike_recorder`` counts the spikes of the target\nneuron and the total number is read out at the end of the simulation\nperiod. The return value of ``output_rate()`` is the firing rate of the\ntarget neuron in Hz.\n\nSecond, the scipy function ``bisect`` is used to determine the optimal\nfiring rate of the neurons of the inhibitory population.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "in_rate = bisect(lambda x: output_rate(x) - r_ex, lower, upper, xtol=prec)\nprint(\"Optimal rate for the inhibitory population: %.2f Hz\" % in_rate)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The function ``bisect`` takes four arguments: first a function whose\nzero crossing is to be determined. Here, the firing rate of the target\nneuron should equal the firing rate of the neurons of the excitatory\npopulation. Thus we define an anonymous function (using `lambda`) that\nreturns the difference between the actual rate of the target neuron and the\nrate of the excitatory Poisson generator, given a rate for the inhibitory\nneurons. The next two arguments are the lower and upper bound of the\ninterval in which to search for the zero crossing. The fourth argument of\n``bisect`` is the desired relative precision of the zero crossing.\n\nFinally, we plot the target neuron's membrane potential as a function of\ntime.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.voltage_trace.from_device(voltmeter)\nnest.voltage_trace.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
}