{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Example for the quantal_stp_synapse\n\nThe ``quantal_stp_synapse`` is a stochastic version of the Tsodys-Markram model\nfor synaptic short term plasticity (STP).\nThis script compares the two variants of the Tsodyks/Markram synapse in NEST.\n\nThis synapse model implements synaptic short-term depression and short-term\nfacilitation according to the quantal release model described by Fuhrmann et\nal. [1]_ and Loebel et al. [2]_.\n\nEach presynaptic spike will stochastically activate a fraction of the\navailable release sites.  This fraction is binomially distributed and the\nrelease probability per site is governed by the Fuhrmann et al. (2002) model.\nThe solution of the differential equations is taken from Maass and Markram\n2002 [3]_.\n\nThe connection weight is interpreted as the maximal weight that can be\nobtained if all n release sites are activated.\n\n## Parameters\n\nThe following parameters can be set in the status dictionary:\n\n* U         - Maximal fraction of available resources [0,1], default=0.5\n* u         - available fraction of resources [0,1], default=0.5\n* p         - probability that a vesicle is available, default = 1.0\n* n         - total number of release sites, default = 1\n* a         - number of available release sites, default = n\n* tau_rec   - time constant for depression in ms, default=800 ms\n* tau_rec   - time constant for facilitation in ms, default=0 (off)\n\n\n## References\n\n.. [1] Fuhrmann G, Segev I, Markram H, and Tsodyks MV. (2002). Coding of\n       temporal information by activity-dependent synapses. Journal of\n       Neurophysiology, 8. https://doi.org/10.1152/jn.00258.2001\n.. [2] Loebel, A., Silberberg, G., Helbig, D., Markram, H., Tsodyks,\n       M. V, & Richardson, M. J. E. (2009). Multiquantal release underlies\n       the distribution of synaptic efficacies in the neocortex. Frontiers\n       in Computational Neuroscience, 3:27. doi:10.3389/neuro.10.027.\n.. [3] Maass W, and Markram H. (2002). Synapses as dynamic memory buffers.\n       Neural Networks, 15(2), 155-161.\n       http://dx.doi.org/10.1016/S0893-6080(01)00144-7\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import nest\nimport numpy\nimport matplotlib.pyplot as plt"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "On average, the ``quantal_stp_synapse`` converges to the ``tsodyks2_synapse``,\nso we can compare the two by running multiple trials.\n\nFirst we define simulation time step and random seed\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "resolution = 0.1  # [ms]\nseed = 12345\n\n# We define the number of trials as well as the number of release sites.\n\nn_sites = 10.0  # number of synaptic release sites\nn_trials = 500  # number of measurement trials\n\n# The pre-synaptic neuron is driven by an injected current for a part of each\n# simulation cycle. We define here the parameters for this stimulation cycle.\n\nI_stim = 376.0   # [pA] stimulation current\nT_on = 500.0     # [ms] stimulation is on\nT_off = 1000.0   # [ms] stimulation is off\n\nT_cycle = T_on + T_off   # total duration of each stimulation cycle"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, we define parameter sets for facilitation and initial weight.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fac_params = {\"U\": 0.02,\n              \"u\": 0.02,\n              \"tau_fac\": 500.,\n              \"tau_rec\": 200.,\n              \"weight\": 1.}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Then, we assign the parameter set to the synapse models\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "tsyn_params = fac_params  # for tsodyks2_synapse\nqsyn_params = tsyn_params.copy()  # for quantal_stp_synapse\n\ntsyn_params[\"x\"] = tsyn_params[\"U\"]\nqsyn_params[\"n\"] = n_sites"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To make the responses comparable, we have to scale the weight by the\nnumber of release sites.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "qsyn_params[\"weight\"] = 1. / n_sites"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We reset NEST to have a well-defined starting point,\nmake NEST less verbose, and set some kernel attributes.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.ResetKernel()\nnest.set_verbosity(\"M_ERROR\")\nnest.resolution = resolution\nnest.rng_seed = seed"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We create three different neurons.\nNeuron one is the sender, the two other neurons receive the synapses.\nWe exploit Python's unpacking mechanism to assign the neurons to named\nvariables directly.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pre_neuron, tsyn_neuron, qsyn_neuron = nest.Create(\"iaf_psc_exp\",\n                                                   params={\"tau_syn_ex\": 3.},\n                                                   n=3)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We create two voltmeters, one for each of the postsynaptic neurons.\nWe start recording only after a first cycle, which is used for equilibration.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "tsyn_voltmeter, qsyn_voltmeter = nest.Create(\"voltmeter\",\n                                             params={\"start\": T_cycle,\n                                                     \"interval\": resolution},\n                                             n=2)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Connect one neuron with the deterministic tsodyks2 synapse and the other neuron\nwith the stochastic quantal stp synapse; then, connect a voltmeter to each neuron.\nHere, ``**tsyn_params`` inserts the content of the ``tsyn_params`` dict into the\ndict passed to ``syn_spec``.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.Connect(pre_neuron, tsyn_neuron,\n             syn_spec={\"synapse_model\": \"tsodyks2_synapse\", **tsyn_params})\n\n# For technical reasons, we currently must set the parameters of the\n# quantal_stp_synapse via default values. This will change in a future version\n# of NEST.\nnest.SetDefaults(\"quantal_stp_synapse\", qsyn_params)\nnest.Connect(pre_neuron, qsyn_neuron, syn_spec={\"synapse_model\": \"quantal_stp_synapse\"})\n\nnest.Connect(tsyn_voltmeter, tsyn_neuron)\nnest.Connect(qsyn_voltmeter, qsyn_neuron)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This loop runs over the `n_trials` trials and performs a standard protocol\nof a high-rate response, followed by a pause and then a recovery response.\n\nWe actually run over ``n_trials + 1`` rounds, since the first trial is for\nequilibration and is not recorded (see voltmeter parameters above).\n\nWe use the NEST ``:class:.RunManager`` to improve performance and call ``:func:.Run``\ninside for each part of the simulation.\n\nWe print a line of breadcrumbs to indicate progress.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(f\"Simulating {n_trials} times \", end=\"\", flush=True)\nwith nest.RunManager():\n    for t in range(n_trials + 1):\n        pre_neuron.I_e = I_stim\n        nest.Run(T_on)\n\n        pre_neuron.I_e = 0.0\n        nest.Run(T_off)\n\n        if t % 10 == 0:\n            print(\".\", end=\"\", flush=True)\nprint()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Simulate one additional time step. This ensures that the\nvoltage traces for all trials, including the last, have the full length, so we\ncan easily transform them into a matrix below.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.Simulate(nest.resolution)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Extract voltage traces and reshape the matrix with one column per trial\nand one row per time step. NEST returns results as NumPy arrays.\nWe extract times only once and keep only times for a single trial.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "vm_tsyn = tsyn_voltmeter.get(\"events\", \"V_m\")\nvm_qsyn = qsyn_voltmeter.get(\"events\", \"V_m\")\n\nsteps_per_trial = round(T_cycle / tsyn_voltmeter.get(\"interval\"))\nvm_tsyn.shape = (n_trials, steps_per_trial)\nvm_qsyn.shape = (n_trials, steps_per_trial)\n\nt_vm = tsyn_voltmeter.get(\"events\", \"times\")\nt_trial = t_vm[:steps_per_trial]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now compute the mean of all trials and plot against trials and references.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "vm_tsyn_mean = vm_tsyn.mean(axis=0)\nvm_qsyn_mean = vm_qsyn.mean(axis=0)\nrms_error = ((vm_tsyn_mean - vm_qsyn_mean) ** 2).mean() ** 0.5\n\nplt.plot(t_trial, vm_tsyn_mean, lw=2, alpha=0.7,\n         label=\"Tsodyks-2 synapse (deterministic)\")\nplt.plot(t_trial, vm_qsyn_mean, lw=2, alpha=0.7,\n         label=\"Quantal STP synapse (stochastic)\")\nplt.xlabel(\"Time [ms]\")\nplt.ylabel(\"Membrane potential [mV]\")\nplt.title(\"Comparison of deterministic and stochastic plasicity rules\")\nplt.text(0.95, 0.05, f\"RMS error: {rms_error:.3g}\",\n         horizontalalignment=\"right\", verticalalignment=\"bottom\",\n         transform=plt.gca().transAxes)  # relative coordinates for text placement\nplt.legend()\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
}