{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Intrinsic currents subthreshold\n\nThis example illustrates how to record from a model with multiple\nintrinsic currents and visualize the results. This is illustrated\nusing the ``ht_neuron`` which has four intrinsic currents: ``I_NaP``,\n``I_KNa``, ``I_T``, and ``I_h``. It is a slightly simplified implementation of\nneuron model proposed in [1]_.\n\nThe neuron is driven by DC current, which is alternated\nbetween depolarizing and hyperpolarizing. Hyperpolarization\nintervals become increasingly longer.\n\n## References\n\n.. [1] Hill and Tononi (2005) Modeling Sleep and Wakefulness in the\n       Thalamocortical System J Neurophysiol 93:1671\n       http://dx.doi.org/10.1152/jn.00915.2004.\n\n## See Also\n\n:doc:`intrinsic_currents_spiking`\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We imported all necessary modules for simulation, analysis and plotting.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import nest\nimport matplotlib.pyplot as plt"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Additionally, we set the verbosity using ``set_verbosity`` to suppress info\nmessages. We also reset the kernel to be sure to start with a clean NEST.\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": [
        "We define simulation parameters:\n\n- The length of depolarization intervals\n- The length of hyperpolarization intervals\n- The amplitude for de- and hyperpolarizing currents\n- The end of the time window to plot\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "n_blocks = 5\nt_block = 20.\nt_dep = [t_block] * n_blocks\nt_hyp = [t_block * 2 ** n for n in range(n_blocks)]\nI_dep = 10.\nI_hyp = -5.\n\nt_end = 500."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We create the one neuron instance and the DC current generator and store\nthe returned handles.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nrn = nest.Create('ht_neuron')\ndc = nest.Create('dc_generator')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We create a multimeter to record\n\n- membrane potential ``V_m``\n- threshold value ``theta``\n- intrinsic currents ``I_NaP``, ``I_KNa``, ``I_T``, ``I_h``\n\nby passing these names in the ``record_from`` list.\n\nTo find out which quantities can be recorded from a given neuron,\nrun::\n\n  nest.GetDefaults('ht_neuron')['recordables']\n\nThe result will contain an entry like::\n\n  <SLILiteral: V_m>\n\nfor each recordable quantity. You need to pass the value of the\n``SLILiteral``, in this case ``V_m`` in the ``record_from`` list.\n\nWe want to record values with 0.1 ms resolution, so we set the\nrecording interval as well; the default recording resolution is 1 ms.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# create multimeter and configure it to record all information\n# we want at 0.1 ms resolution\nmm = nest.Create('multimeter',\n                 params={'interval': 0.1,\n                         'record_from': ['V_m', 'theta',\n                                         'I_NaP', 'I_KNa', 'I_T', 'I_h']}\n                 )"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We connect the DC generator and the multimeter to the neuron. Note that\nthe multimeter, just like the voltmeter is connected to the neuron,\nnot the neuron to the multimeter.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nest.Connect(dc, nrn)\nnest.Connect(mm, nrn)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We are ready to simulate. We alternate between driving the neuron with\ndepolarizing and hyperpolarizing currents. Before each simulation\ninterval, we set the amplitude of the DC generator to the correct value.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for t_sim_dep, t_sim_hyp in zip(t_dep, t_hyp):\n\n    dc.amplitude = I_dep\n    nest.Simulate(t_sim_dep)\n\n    dc.amplitude = I_hyp\n    nest.Simulate(t_sim_hyp)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We now fetch the data recorded by the multimeter. The data are returned as\na dictionary with entry ``times`` containing timestamps for all recorded\ndata, plus one entry per recorded quantity.\n\nAll data is contained in the ``events`` entry of the status dictionary\nreturned by the multimeter. Because all NEST function return arrays,\nwe need to pick out element `0` from the result of ``GetStatus``.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "data = mm.events\nt = data['times']"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The next step is to plot the results. We create a new figure, add a single\nsubplot and plot at first membrane potential and threshold.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig = plt.figure()\nVax = fig.add_subplot(111)\nVax.plot(t, data['V_m'], 'b-', lw=2, label=r'$V_m$')\nVax.plot(t, data['theta'], 'g-', lw=2, label=r'$\\Theta$')\nVax.set_ylim(-80., 0.)\nVax.set_ylabel('Voltageinf [mV]')\nVax.set_xlabel('Time [ms]')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To plot the input current, we need to create an input current trace. We\nconstruct it from the durations of the de- and hyperpolarizing inputs and\nadd the delay in the connection between DC generator and neuron:\n\n* We find the delay by checking the status of the dc->nrn connection.\n* We find the resolution of the simulation from the kernel status.\n* Each current interval begins one time step after the previous interval,\n  is delayed by the delay and effective for the given duration.\n* We build the time axis incrementally. We only add the delay when adding\n  the first time point after t=0. All subsequent points are then\n  automatically shifted by the delay.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "conns = nest.GetConnections(dc, nrn)\ndelay = conns.delay\ndt = nest.GetKernelStatus('resolution')\n\nt_dc, I_dc = [0], [0]\n\nfor td, th in zip(t_dep, t_hyp):\n    t_prev = t_dc[-1]\n    t_start_dep = t_prev + dt if t_prev > 0 else t_prev + dt + delay\n    t_end_dep = t_start_dep + td\n    t_start_hyp = t_end_dep + dt\n    t_end_hyp = t_start_hyp + th\n\n    t_dc.extend([t_start_dep, t_end_dep, t_start_hyp, t_end_hyp])\n    I_dc.extend([I_dep, I_dep, I_hyp, I_hyp])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The following function turns a name such as ``I_NaP`` into proper TeX code\n$I_{\\mathrm{NaP}}$ for a pretty label.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def texify_name(name):\n    return r'${}_{{\\mathrm{{{}}}}}$'.format(*name.split('_'))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next, we add a right vertical axis and plot the currents with respect to\nthat axis.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "Iax = Vax.twinx()\nIax.plot(t_dc, I_dc, 'k-', lw=2, label=texify_name('I_DC'))\n\nfor iname, color in (('I_h', 'maroon'), ('I_T', 'orange'),\n                     ('I_NaP', 'crimson'), ('I_KNa', 'aqua')):\n    Iax.plot(t, data[iname], color=color, lw=2, label=texify_name(iname))\n\nIax.set_xlim(0, t_end)\nIax.set_ylim(-10., 15.)\nIax.set_ylabel('Current [pA]')\nIax.set_title('ht_neuron driven by DC current')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We need to make a little extra effort to combine lines from the two axis\ninto one legend.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "lines_V, labels_V = Vax.get_legend_handles_labels()\nlines_I, labels_I = Iax.get_legend_handles_labels()\ntry:\n    Iax.legend(lines_V + lines_I, labels_V + labels_I, fontsize='small')\nexcept TypeError:\n    # work-around for older Matplotlib versions\n    Iax.legend(lines_V + lines_I, labels_V + labels_I)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Note that ``I_KNa`` is not activated in this example because the neuron does\nnot spike. ``I_T`` has only a very small amplitude.\n\n"
      ]
    }
  ],
  "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
}