{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Pulse packet example\n\nThis script compares the average and individual membrane potential excursions\nin response to a single pulse packet with an analytically acquired voltage\ntrace (see: Diesmann [1]_)\nA pulse packet is a transient spike volley with a Gaussian rate profile.\nThe user can specify the neural parameters, the parameters of the\npulse-packet and the number of trials.\n\n\n## References\n\n.. [1] Diesmann M. 2002. Dissertation. Conditions for stable propagation of\n       synchronous spiking in cortical neural networks: Single neuron dynamics\n       and network properties.\n       http://d-nb.info/968772781/34.\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First, we import all necessary modules for simulation, analysis and\nplotting.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import scipy.special as sp\nimport nest\nimport numpy\nimport matplotlib.pyplot as plt\n\n# Properties of pulse packet:\n\na = 100            # number of spikes in one pulse packet\nsdev = 10.         # width of pulse packet (ms)\nweight = 0.1       # PSP amplitude (mV)\npulsetime = 500.   # occurrence time (center) of pulse-packet (ms)\n\n\n# Network and neuron characteristics:\n\nn_neurons = 100    # number of neurons\ncm = 200.          # membrane capacitance (pF)\ntau_s = 0.5        # synaptic time constant (ms)\ntau_m = 20.        # membrane time constant (ms)\nV0 = 0.0           # resting potential (mV)\nVth = numpy.inf    # firing threshold, high value to avoid spiking\n\n\n# Simulation and analysis parameters:\n\nsimtime = 1000.               # how long we simulate (ms)\nsimulation_resolution = 0.1  # (ms)\nsampling_resolution = 1.   # for voltmeter (ms)\nconvolution_resolution = 1.   # for the analytics (ms)\n\n\n# Some parameters in base units.\n\nCm = cm * 1e-12            # convert to Farad\nWeight = weight * 1e-12    # convert to Ampere\nTau_s = tau_s * 1e-3       # convert to sec\nTau_m = tau_m * 1e-3       # convert to sec\nSdev = sdev * 1e-3         # convert to sec\nConvolution_resolution = convolution_resolution * 1e-3  # convert to sec"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This function calculates the membrane potential excursion in response\nto a single input spike (the equation is given for example in Diesmann [1]_,\neq.2.3).\nIt expects:\n\n* ``Time``: a time array or a single time point (in sec)\n* ``Tau_s`` and ``Tau_m``: the synaptic and the membrane time constant (in sec)\n* ``Cm``: the membrane capacity (in Farad)\n* ``Weight``: the synaptic weight (in Ampere)\n\nIt returns the provoked membrane potential (in mV)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def make_psp(Time, Tau_s, Tau_m, Cm, Weight):\n    term1 = (1 / Tau_s - 1 / Tau_m)\n    term2 = numpy.exp(-Time / Tau_s)\n    term3 = numpy.exp(-Time / Tau_m)\n    PSP = (Weight / Cm * numpy.exp(1) / Tau_s *\n           (((-Time * term2) / term1) + (term3 - term2) / term1 ** 2))\n    return PSP * 1e3"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This function finds the exact location of the maximum of the PSP caused by a\nsingle input spike. The location is obtained by setting the first derivative\nof the equation for the PSP (see ``make_psp()``) to zero. The resulting\nequation can be expressed in terms of a `LambertW function`.\nThis function expects:\n\n* ``Tau_s`` and ``Tau_m``: the synaptic and membrane time constant (in sec)\n\nIt returns the location of the maximum (in sec)\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 find_loc_pspmax(tau_s, tau_m):\n    var = tau_m / tau_s\n    lam = LambertWm1(-numpy.exp(-1 / var) / var)\n    t_maxpsp = (-var * lam - 1) / var / (1 / tau_s - 1 / tau_m) * 1e-3\n    return t_maxpsp"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First, we construct a Gaussian kernel for a given standard derivation\n(``sig``) and mean value (``mu``). In this case the standard derivation is\nthe width of the pulse packet (see [1]_).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sig = Sdev\nmu = 0.0\nx = numpy.arange(-4 * sig, 4 * sig, Convolution_resolution)\nterm1 = 1 / (sig * numpy.sqrt(2 * numpy.pi))\nterm2 = numpy.exp(-(x - mu) ** 2 / (sig ** 2 * 2))\ngauss = term1 * term2 * Convolution_resolution"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Second, we calculate the PSP of a neuron due to a single spiking input.\n(see Diesmann 2002, eq. 2.3).\nSince we do that in discrete time steps, we first construct an array\n(``t_psp``) that contains the time points we want to consider. Then, the\nfunction ``make_psp()`` (that creates the PSP) takes the time array as its\nfirst argument.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "t_psp = numpy.arange(0, 10 * (Tau_m + Tau_s), Convolution_resolution)\npsp = make_psp(t_psp, Tau_s, Tau_m, Cm, Weight)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now, we want to normalize the PSP amplitude to one. We therefore have to\ndivide the PSP by its maximum ([1]_ sec 6.1). The function\n``find_loc_pspmax()`` returns the exact time point (``t_pspmax``) when we\nexpect the maximum to occur. The function ``make_psp()`` calculates the\ncorresponding PSP value, which is our PSP amplitude (``psp_amp``).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "t_pspmax = find_loc_pspmax(Tau_s, Tau_m)\npsp_amp = make_psp(t_pspmax, Tau_s, Tau_m, Cm, Weight)\npsp_norm = psp / psp_amp"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we have all ingredients to compute the membrane potential excursion\n(`U`). This calculation implies a convolution of the Gaussian with the\nnormalized PSP (see [1]_, eq. 6.9). In order to avoid an offset in the\nconvolution, we need to add a pad of zeros on the left side of the\nnormalized PSP. Later on we want to compare our analytical results with the\nsimulation outcome. Therefore we need a time vector (`t_U`) with the correct\ntemporal resolution, which places the excursion of the potential at the\ncorrect time.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "psp_norm = numpy.pad(psp_norm, [len(psp_norm) - 1, 1])\nU = a * psp_amp * numpy.convolve(gauss, psp_norm)\nulen = len(U)\nt_U = (convolution_resolution * numpy.linspace(-ulen / 2., ulen / 2., ulen) +\n       pulsetime + 1.)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In this section we simulate a network of multiple neurons.\nAll these neurons receive an individual pulse packet that is drawn from a\nGaussian distribution.\n\nWe reset the Kernel, define the simulation resolution and set the\nverbosity using ``set_verbosity`` to suppress info messages.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.ResetKernel()\nnest.SetKernelStatus({'resolution': simulation_resolution})\nnest.set_verbosity(\"M_WARNING\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Afterwards we create several neurons, the same amount of\npulse-packet-generators and a voltmeter. All these nodes/devices\nhave specific properties that are specified in device specific\ndictionaries (here: `neuron_pars` for the neurons, `ppg_pars`\nfor the and pulse-packet-generators and `vm_pars` for the voltmeter).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "neuron_pars = {\n    'V_th': Vth,\n    'tau_m': tau_m,\n    'tau_syn_ex': tau_s,\n    'C_m': cm,\n    'E_L': V0,\n    'V_reset': V0,\n    'V_m': V0\n}\nneurons = nest.Create('iaf_psc_alpha', n_neurons, neuron_pars)\nppg_pars = {\n    'pulse_times': [pulsetime],\n    'activity': a,\n    'sdev': sdev\n}\nppgs = nest.Create('pulsepacket_generator', n_neurons, ppg_pars)\nvm_pars = {'interval': sampling_resolution}\nvm = nest.Create('voltmeter', 1, vm_pars)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now, we connect each pulse generator to one neuron via static synapses.\nWe want to keep all properties of the static synapse constant except the\nsynaptic weight. Therefore we change the weight with  the help of the command\n``SetDefaults``.\nThe command ``Connect`` connects all kinds of nodes/devices. Since multiple\nnodes/devices can be connected in different ways e.g., each source connects\nto all targets, each source connects to a subset of targets or each source\nconnects to exactly one target, we have to specify the connection. In our\ncase we use the ``one_to_one`` connection routine since we connect one pulse\ngenerator (source) to one neuron (target).\nIn addition we also connect the `voltmeter` to the `neurons`.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.SetDefaults('static_synapse', {'weight': weight})\nnest.Connect(ppgs, neurons, 'one_to_one')\nnest.Connect(vm, neurons)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In the next step we run the simulation for a given duration in ms.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.Simulate(simtime)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally, we record the membrane potential, when it occurred and to which\nneuron it belongs. The sender and the time point of a voltage\ndata point at position x in the voltage array (``V_m``), can be found at the\nsame position x in the sender (`senders`) and the time array (`times`).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "Vm = vm.get('events', 'V_m')\ntimes = vm.get('events', 'times')\nsenders = vm.get('events', 'senders')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here we plot the membrane potential derived from the theory and from the\nsimulation. Since we simulate multiple neurons that received slightly\ndifferent pulse packets, we plot the individual and the averaged membrane\npotentials.\n\nWe plot the analytical solution U (the resting potential V0 shifts the\nmembrane potential up or downwards).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.plot(t_U, U + V0, 'r', lw=2, zorder=3, label='analytical solution')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Then we plot all individual membrane potentials.\nThe time axes is the range of the simulation time in steps of ms.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "Vm_single = [Vm[senders == n.global_id] for n in neurons]\nsimtimes = numpy.arange(1, simtime)\nfor idn in range(n_neurons):\n    if idn == 0:\n        plt.plot(simtimes, Vm_single[idn], 'gray',\n                 zorder=1, label='single potentials')\n    else:\n        plt.plot(simtimes, Vm_single[idn], 'gray', zorder=1)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally, we plot the averaged membrane potential.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "Vm_average = numpy.mean(Vm_single, axis=0)\nplt.plot(simtimes, Vm_average, 'b', lw=4,\n         zorder=2, label='averaged potential')\nplt.legend()\nplt.xlabel('time (ms)')\nplt.ylabel('membrane potential (mV)')\nplt.xlim((-5 * (tau_m + tau_s) + pulsetime,\n          10 * (tau_m + tau_s) + pulsetime))\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
}