{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Sensitivity to perturbation\n\nThis script simulates a network in two successive trials, which are identical\nexcept for one extra input spike in the second realisation (a small\nperturbation). The network consists of recurrent, randomly connected excitatory\nand inhibitory neurons. Its activity is driven by an external Poisson input\nprovided to all neurons independently. In order to ensure that the network is\nreset appropriately between the trials, we do the following steps:\n\n- resetting the network\n- resetting the random network generator\n- resetting the internal clock\n- deleting all entries in the spike recorder\n- introducing a hyperpolarisation phase between the trials\n  (in order to avoid that spikes remaining in the NEST memory\n  after the first simulation are fed into the second simulation)\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Importing all necessary modules for simulation, analysis and plotting.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy\nimport matplotlib.pyplot as plt\nimport nest"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here we define all parameters necessary for building and simulating the\nnetwork.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# We start with the global network parameters.\n\n\nNE = 1000      # number of excitatory neurons\nNI = 250       # number of inhibitory neurons\nN = NE + NI    # total number of neurons\nKE = 100       # excitatory in-degree\nKI = 25        # inhibitory in-degree"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Parameters specific for the neurons in the network. The  default values of\nthe reset potential ``E_L`` and the spiking threshold ``V_th`` are used to set\nthe limits of the initial potential of the neurons.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "neuron_model = 'iaf_psc_delta'\nneuron_params = nest.GetDefaults(neuron_model)\nVmin = neuron_params['E_L']   # minimum of initial potential distribution (mV)\nVmax = neuron_params['V_th']  # maximum of initial potential distribution (mV)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Synapse parameters. Changing the weights `J` in the network can lead to\nqualitatively different behaviors. If `J` is small (e.g. ``J = 0.1``), we\nare likely to observe a non-chaotic network behavior (after perturbation\nthe network returns to its original activity). Increasing `J`\n(e.g ``J = 5.5``) leads to rather chaotic activity. Given that in this\nexample the transition to chaos is probabilistic, we sometimes observe\nchaotic behavior for small weights (e.g. ``J = 0.5``) and non-chaotic\nbehavior for strong weights (e.g. ``J = 5.4``).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "J = 0.5                   # excitatory synaptic weight (mV)\ng = 6.                    # relative inhibitory weight\ndelay = 0.1               # spike transmission delay (ms)\n\n\n# External input parameters.\n\n\nJext = 0.2                # PSP amplitude for external Poisson input (mV)\nrate_ext = 6500.          # rate of the external Poisson input\n\n\n# Perturbation parameters.\n\n\nt_stim = 400.             # perturbation time (time of the extra spike)\nJstim = Jext              # perturbation amplitude (mV)\n\n\n# Simulation parameters.\n\n\nT = 1000.                 # simulation time per trial (ms)\nfade_out = 2. * delay     # fade out time (ms)\ndt = 0.01                 # simulation time resolution (ms)\nseed_NEST = 30            # seed of random number generator in Nest\nseed_numpy = 30           # seed of random number generator in numpy\n\nsenders = []\nspiketimes = []"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# we run the two simulations successively. After each simulation the\n# sender ids and spiketimes are stored in a list (``senders``, ``spiketimes``).\n\nfor trial in [0, 1]:\n\n    # Before we build the network, we reset the simulation kernel to ensure\n    # that previous NEST simulations in the Python shell will not disturb this\n    # simulation and set the simulation resolution (later defined\n    # synaptic delays cannot be smaller than the simulation resolution).\n    nest.ResetKernel()\n    nest.SetKernelStatus({\"resolution\": dt})\n\n    ###############################################################################\n    # Now we start building the network and create excitatory and inhibitory nodes\n    # and connect them. According to the connectivity specification, each neuron\n    # is assigned random KE synapses from the excitatory population and random KI\n    # synapses from the inhibitory population.\n\n    nodes_ex = nest.Create(neuron_model, NE)\n    nodes_in = nest.Create(neuron_model, NI)\n    allnodes = nodes_ex + nodes_in\n\n    nest.Connect(nodes_ex, allnodes,\n                 conn_spec={'rule': 'fixed_indegree', 'indegree': KE},\n                 syn_spec={'weight': J, 'delay': dt})\n    nest.Connect(nodes_in, allnodes,\n                 conn_spec={'rule': 'fixed_indegree', 'indegree': KI},\n                 syn_spec={'weight': -g * J, 'delay': dt})\n\n    ###############################################################################\n    # Afterwards we create a ``poisson_generator`` that provides spikes (the external\n    # input) to the neurons until time ``T`` is reached.\n    # Afterwards a ``dc_generator``, which is also connected to the whole population,\n    # provides a stong hyperpolarisation step for a short time period ``fade_out``.\n    #\n    # The ``fade_out`` period has to last at least twice as long as the simulation\n    # resolution to supress the neurons from firing.\n\n    ext = nest.Create(\"poisson_generator\",\n                      params={'rate': rate_ext, 'stop': T})\n    nest.Connect(ext, allnodes,\n                 syn_spec={'weight': Jext, 'delay': dt})\n\n    suppr = nest.Create(\"dc_generator\",\n                        params={'amplitude': -1e16, 'start': T,\n                                'stop': T + fade_out})\n    nest.Connect(suppr, allnodes)\n\n    spikerecorder = nest.Create(\"spike_recorder\")\n    nest.Connect(allnodes, spikerecorder)\n\n    ###############################################################################\n    # We then create the ``spike_generator``, which provides the extra spike\n    # (perturbation).\n\n    stimulus = nest.Create(\"spike_generator\")\n    stimulus.spike_times = []\n\n    ###############################################################################\n    # We need to reset the random number generator and the clock of\n    # the simulation Kernel. In addition, we ensure that there is no spike left in\n    # the spike recorder.\n\n    nest.SetKernelStatus({\"rng_seeds\": [seed_NEST], 'time': 0.0})\n    spikerecorder.n_events = 0\n\n    # We assign random initial membrane potentials to all neurons\n\n    numpy.random.seed(seed_numpy)\n    Vms = Vmin + (Vmax - Vmin) * numpy.random.rand(N)\n    allnodes.V_m = Vms\n\n    ##############################################################################\n    # In the second trial, we add an extra input spike at time ``t_stim`` to the\n    # neuron that fires first after perturbation time ``t_stim``. Thus, we make sure\n    # that the perturbation is transmitted to the network before it fades away in\n    # the perturbed neuron. (Single IAF-neurons are not chaotic.)\n\n    if trial == 1:\n        id_stim = [senders[0][spiketimes[0] > t_stim][0]]\n        nest.Connect(stimulus, nest.NodeCollection(id_stim),\n                     syn_spec={'weight': Jstim, 'delay': dt})\n        stimulus.spike_times = [t_stim]\n\n    # Now we simulate the network and add a fade out period to discard\n    # remaining spikes.\n\n    nest.Simulate(T)\n    nest.Simulate(fade_out)\n\n    # Storing the data.\n\n    senders += [spikerecorder.get('events', 'senders')]\n    spiketimes += [spikerecorder.get('events', 'times')]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We plot the spiking activity of the network (first trial in red, second trial\nin black).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.figure(1)\nplt.clf()\nplt.plot(spiketimes[0], senders[0], 'ro', ms=4.)\nplt.plot(spiketimes[1], senders[1], 'ko', ms=2.)\nplt.xlabel('time (ms)')\nplt.ylabel('neuron id')\nplt.xlim((0, T))\nplt.ylim((0, 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.6.12"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}