{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Numerical phase-plane analysis of the Hodgkin-Huxley neuron\n\nhh_phaseplane makes a numerical phase-plane analysis of the Hodgkin-Huxley\nneuron (``hh_psc_alpha``). Dynamics is investigated in the V-n space (see remark\nbelow). A constant DC can be specified  and its influence on the nullclines\ncan be studied.\n\n## Remark\n\nTo make the two-dimensional analysis possible, the (four-dimensional)\nHodgkin-Huxley formalism needs to be artificially reduced to two dimensions,\nin this case by 'clamping' the two other variables, `m` and `h`, to\nconstant values (`m_eq` and `h_eq`).\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import nest\nimport numpy as np\nfrom matplotlib import pyplot as plt\n\n\namplitude = 100.  # Set externally applied current amplitude in pA\ndt = 0.1  # simulation step length [ms]\n\nv_min = -100.  # Min membrane potential\nv_max = 42.  # Max membrane potential\nn_min = 0.1  # Min inactivation variable\nn_max = 0.81  # Max inactivation variable\ndelta_v = 2.  # Membrane potential step length\ndelta_n = 0.01  # Inactivation variable step length\n\nV_vec = np.arange(v_min, v_max, delta_v)\nn_vec = np.arange(n_min, n_max, delta_n)\n\nnum_v_steps = len(V_vec)\nnum_n_steps = len(n_vec)\n\nnest.ResetKernel()\nnest.set_verbosity('M_ERROR')\nnest.resolution = dt\n\nneuron = nest.Create('hh_psc_alpha')\n\n# Numerically obtain equilibrium state\nnest.Simulate(1000)\n\nm_eq = neuron.Act_m\nh_eq = neuron.Inact_h\n\nneuron.I_e = amplitude  # Apply external current\n\n# Scan state space\nprint('Scanning phase space')\n\nV_matrix = np.zeros([num_n_steps, num_v_steps])\nn_matrix = np.zeros([num_n_steps, num_v_steps])\n\n# pp_data will contain the phase-plane data as a vector field\npp_data = np.zeros([num_n_steps * num_v_steps, 4])\n\ncount = 0\nfor i, V in enumerate(V_vec):\n    for j, n in enumerate(n_vec):\n        # Set V_m and n\n        neuron.set(V_m=V, Act_n=n, Act_m=m_eq, Inact_h=h_eq)\n        # Find state\n        V_m = neuron.V_m\n        Act_n = neuron.Act_n\n\n        # Simulate a short while\n        nest.Simulate(dt)\n\n        # Find difference between new state and old state\n        V_m_new = neuron.V_m - V\n        Act_n_new = neuron.Act_n - n\n\n        # Store in vector for later analysis\n        V_matrix[j, i] = abs(V_m_new)\n        n_matrix[j, i] = abs(Act_n_new)\n        pp_data[count] = np.array([V_m, Act_n, V_m_new, Act_n_new])\n\n        if count % 10 == 0:\n            # Write updated state next to old state\n            print('')\n            print('Vm:  \\t', V_m)\n            print('new Vm:\\t', V_m_new)\n            print('Act_n:', Act_n)\n            print('new Act_n:', Act_n_new)\n\n        count += 1\n\n# Set state for AP generation\nneuron.set(V_m=-34., Act_n=0.2, Act_m=m_eq, Inact_h=h_eq)\n\nprint('')\nprint('AP-trajectory')\n# ap will contain the trace of a single action potential as one possible\n# numerical solution in the vector field\nap = np.zeros([1000, 2])\nfor i in range(1000):\n    # Find state\n    V_m = neuron.V_m\n    Act_n = neuron.Act_n\n\n    if i % 10 == 0:\n        # Write new state next to old state\n        print('Vm: \\t', V_m)\n        print('Act_n:', Act_n)\n    ap[i] = np.array([V_m, Act_n])\n\n    # Simulate again\n    neuron.set(Act_m=m_eq, Inact_h=h_eq)\n    nest.Simulate(dt)\n\n# Make analysis\nprint('')\nprint('Plot analysis')\n\nnullcline_V = []\nnullcline_n = []\n\nprint('Searching nullclines')\nfor i in range(0, len(V_vec)):\n    index = np.nanargmin(V_matrix[:][i])\n    if index != 0 and index != len(n_vec):\n        nullcline_V.append([V_vec[i], n_vec[index]])\n\n    index = np.nanargmin(n_matrix[:][i])\n    if index != 0 and index != len(n_vec):\n        nullcline_n.append([V_vec[i], n_vec[index]])\n\nprint('Plotting vector field')\nfactor = 0.1\nfor i in range(0, np.shape(pp_data)[0], 3):\n    plt.plot([pp_data[i][0], pp_data[i][0] + factor * pp_data[i][2]],\n             [pp_data[i][1], pp_data[i][1] + factor * pp_data[i][3]],\n             color=[0.6, 0.6, 0.6])\n\nplt.plot(nullcline_V[:][0], nullcline_V[:][1], linewidth=2.0)\nplt.plot(nullcline_n[:][0], nullcline_n[:][1], linewidth=2.0)\n\nplt.xlim([V_vec[0], V_vec[-1]])\nplt.ylim([n_vec[0], n_vec[-1]])\n\nplt.plot(ap[:][0], ap[:][1], color='black', linewidth=1.0)\n\nplt.xlabel('Membrane potential V [mV]')\nplt.ylabel('Inactivation variable n')\nplt.title('Phase space of the Hodgkin-Huxley Neuron')\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.8.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}