.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/hh_phaseplane.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        Click :ref:`here <sphx_glr_download_auto_examples_hh_phaseplane.py>`
        to download the full example code

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_hh_phaseplane.py:


Numerical phase-plane analysis of the Hodgkin-Huxley neuron
-----------------------------------------------------------

hh_phaseplane makes a numerical phase-plane analysis of the Hodgkin-Huxley
neuron (``hh_psc_alpha``). Dynamics is investigated in the V-n space (see remark
below). A constant DC can be specified  and its influence on the nullclines
can be studied.

Remark
~~~~~~

To make the two-dimensional analysis possible, the (four-dimensional)
Hodgkin-Huxley formalism needs to be artificially reduced to two dimensions,
in this case by 'clamping' the two other variables, `m` and `h`, to
constant values (`m_eq` and `h_eq`).

.. GENERATED FROM PYTHON SOURCE LINES 40-176

.. code-block:: default


    import nest
    import numpy as np
    from matplotlib import pyplot as plt


    amplitude = 100.  # Set externally applied current amplitude in pA
    dt = 0.1  # simulation step length [ms]

    v_min = -100.  # Min membrane potential
    v_max = 42.  # Max membrane potential
    n_min = 0.1  # Min inactivation variable
    n_max = 0.81  # Max inactivation variable
    delta_v = 2.  # Membrane potential step length
    delta_n = 0.01  # Inactivation variable step length

    V_vec = np.arange(v_min, v_max, delta_v)
    n_vec = np.arange(n_min, n_max, delta_n)

    num_v_steps = len(V_vec)
    num_n_steps = len(n_vec)

    nest.ResetKernel()
    nest.set_verbosity('M_ERROR')
    nest.resolution = dt

    neuron = nest.Create('hh_psc_alpha')

    # Numerically obtain equilibrium state
    nest.Simulate(1000)

    m_eq = neuron.Act_m
    h_eq = neuron.Inact_h

    neuron.I_e = amplitude  # Apply external current

    # Scan state space
    print('Scanning phase space')

    V_matrix = np.zeros([num_n_steps, num_v_steps])
    n_matrix = np.zeros([num_n_steps, num_v_steps])

    # pp_data will contain the phase-plane data as a vector field
    pp_data = np.zeros([num_n_steps * num_v_steps, 4])

    count = 0
    for i, V in enumerate(V_vec):
        for j, n in enumerate(n_vec):
            # Set V_m and n
            neuron.set(V_m=V, Act_n=n, Act_m=m_eq, Inact_h=h_eq)
            # Find state
            V_m = neuron.V_m
            Act_n = neuron.Act_n

            # Simulate a short while
            nest.Simulate(dt)

            # Find difference between new state and old state
            V_m_new = neuron.V_m - V
            Act_n_new = neuron.Act_n - n

            # Store in vector for later analysis
            V_matrix[j, i] = abs(V_m_new)
            n_matrix[j, i] = abs(Act_n_new)
            pp_data[count] = np.array([V_m, Act_n, V_m_new, Act_n_new])

            if count % 10 == 0:
                # Write updated state next to old state
                print('')
                print('Vm:  \t', V_m)
                print('new Vm:\t', V_m_new)
                print('Act_n:', Act_n)
                print('new Act_n:', Act_n_new)

            count += 1

    # Set state for AP generation
    neuron.set(V_m=-34., Act_n=0.2, Act_m=m_eq, Inact_h=h_eq)

    print('')
    print('AP-trajectory')
    # ap will contain the trace of a single action potential as one possible
    # numerical solution in the vector field
    ap = np.zeros([1000, 2])
    for i in range(1000):
        # Find state
        V_m = neuron.V_m
        Act_n = neuron.Act_n

        if i % 10 == 0:
            # Write new state next to old state
            print('Vm: \t', V_m)
            print('Act_n:', Act_n)
        ap[i] = np.array([V_m, Act_n])

        # Simulate again
        neuron.set(Act_m=m_eq, Inact_h=h_eq)
        nest.Simulate(dt)

    # Make analysis
    print('')
    print('Plot analysis')

    nullcline_V = []
    nullcline_n = []

    print('Searching nullclines')
    for i in range(0, len(V_vec)):
        index = np.nanargmin(V_matrix[:][i])
        if index != 0 and index != len(n_vec):
            nullcline_V.append([V_vec[i], n_vec[index]])

        index = np.nanargmin(n_matrix[:][i])
        if index != 0 and index != len(n_vec):
            nullcline_n.append([V_vec[i], n_vec[index]])

    print('Plotting vector field')
    factor = 0.1
    for i in range(0, np.shape(pp_data)[0], 3):
        plt.plot([pp_data[i][0], pp_data[i][0] + factor * pp_data[i][2]],
                 [pp_data[i][1], pp_data[i][1] + factor * pp_data[i][3]],
                 color=[0.6, 0.6, 0.6])

    plt.plot(nullcline_V[:][0], nullcline_V[:][1], linewidth=2.0)
    plt.plot(nullcline_n[:][0], nullcline_n[:][1], linewidth=2.0)

    plt.xlim([V_vec[0], V_vec[-1]])
    plt.ylim([n_vec[0], n_vec[-1]])

    plt.plot(ap[:][0], ap[:][1], color='black', linewidth=1.0)

    plt.xlabel('Membrane potential V [mV]')
    plt.ylabel('Inactivation variable n')
    plt.title('Phase space of the Hodgkin-Huxley Neuron')

    plt.show()


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** ( 0 minutes  0.000 seconds)


.. _sphx_glr_download_auto_examples_hh_phaseplane.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example


    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: hh_phaseplane.py <hh_phaseplane.py>`

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: hh_phaseplane.ipynb <hh_phaseplane.ipynb>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_