# Source code for brainmodels.neurons.HH

# -*- coding: utf-8 -*-

import brainpy as bp
import brainpy.math as bm
from .base import Neuron

__all__ = [
'HH'
]

[docs]class HH(Neuron):
r"""Hodgkin–Huxley neuron model.

**Model Descriptions**

The Hodgkin-Huxley (HH; Hodgkin & Huxley, 1952) model [1]_ for the generation of
the nerve action potential is one of the most successful mathematical models of
a complex biological process that has ever been formulated. The basic concepts
expressed in the model have proved a valid approach to the study of bio-electrical
activity from the most primitive single-celled organisms such as *Paramecium*,
right through to the neurons within our own brains.

Mathematically, the model is given by,

.. math::

C \frac {dV} {dt} = -(\bar{g}_{Na} m^3 h (V &-E_{Na})
+ \bar{g}_K n^4 (V-E_K) + g_{leak} (V - E_{leak})) + I(t)

\frac {dx} {dt} &= \alpha_x (1-x)  - \beta_x, \quad x\in {\rm{\{m, h, n\}}}

&\alpha_m(V) = \frac {0.1(V+40)}{1-\exp(\frac{-(V + 40)} {10})}

&\beta_m(V) = 4.0 \exp(\frac{-(V + 65)} {18})

&\alpha_h(V) = 0.07 \exp(\frac{-(V+65)}{20})

&\beta_h(V) = \frac 1 {1 + \exp(\frac{-(V + 35)} {10})}

&\alpha_n(V) = \frac {0.01(V+55)}{1-\exp(-(V+55)/10)}

&\beta_n(V) = 0.125 \exp(\frac{-(V + 65)} {80})

The illustrated example of HH neuron model please see this notebook <../neurons/HH_model.ipynb>_.

The Hodgkin–Huxley model can be thought of as a differential equation system with
four state variables, :math:V_{m}(t),n(t),m(t), and :math:h(t), that change
with respect to time :math:t. The system is difficult to study because it is a
nonlinear system and cannot be solved analytically. However, there are many numeric
methods available to analyze the system. Certain properties and general behaviors,
such as limit cycles, can be proven to exist.

*1. Center manifold*

Because there are four state variables, visualizing the path in phase space can
be difficult. Usually two variables are chosen, voltage :math:V_{m}(t) and the
potassium gating variable :math:n(t), allowing one to visualize the limit cycle.
However, one must be careful because this is an ad-hoc method of visualizing the
4-dimensional system. This does not prove the existence of the limit cycle.

.. image:: ../../images/Hodgkin_Huxley_Limit_Cycle.png
:align: center

A better projection can be constructed from a careful analysis of the Jacobian of
the system, evaluated at the equilibrium point. Specifically, the eigenvalues of
the Jacobian are indicative of the center manifold's existence. Likewise, the
eigenvectors of the Jacobian reveal the center manifold's orientation. The
Hodgkin–Huxley model has two negative eigenvalues and two complex eigenvalues
with slightly positive real parts. The eigenvectors associated with the two
negative eigenvalues will reduce to zero as time :math:t increases. The remaining
two complex eigenvectors define the center manifold. In other words, the
4-dimensional system collapses onto a 2-dimensional plane. Any solution
starting off the center manifold will decay towards the *center manifold*.
Furthermore, the limit cycle is contained on the center manifold.

*2. Bifurcations*

If the injected current :math:I were used as a bifurcation parameter, then the
Hodgkin–Huxley model undergoes a Hopf bifurcation. As with most neuronal models,
increasing the injected current will increase the firing rate of the neuron.
One consequence of the Hopf bifurcation is that there is a minimum firing rate.
This means that either the neuron is not firing at all (corresponding to zero
frequency), or firing at the minimum firing rate. Because of the all-or-none
principle, there is no smooth increase in action potential amplitude, but
rather there is a sudden "jump" in amplitude. The resulting transition is
known as a canard <http://www.scholarpedia.org/article/Canards>_.

.. image:: ../../images/Hodgkins_Huxley_bifurcation_by_I.gif
:align: center

The following image shows the bifurcation diagram of the Hodgkin–Huxley model
as a function of the external drive :math:I [3]_. The green lines show the amplitude
of a stable limit cycle and the blue lines indicate unstable limit-cycle behaviour,
both born from Hopf bifurcations. The solid red line shows the stable fixed point
and the black line shows the unstable fixed point.

.. image:: ../../images/Hodgkin_Huxley_bifurcation.png
:align: center

**Model Examples**

.. plot::
:include-source: True

>>> import brainpy as bp
>>> import brainmodels
>>> group = brainmodels.neurons.HH(2)
>>> runner = bp.StructRunner(group, monitors=['V'], inputs=('input', 10.))
>>> runner.run(200.)
>>> bp.visualize.line_plot(runner.mon.ts, runner.mon.V, show=True)

.. plot::
:include-source: True

>>> import brainpy as bp
>>> import brainmodels
>>> group = brainmodels.neurons.HH(2)
>>> runner = bp.ReportRunner(group, monitors=bp.Monitor(variables=['V'], intervals=[1.]),
>>>                          inputs=('input', 10.), jit=True)
>>> runner.run(200.)
>>> bp.visualize.line_plot(runner.mon['V.t'], runner.mon.V, show=True)

.. plot::
:include-source: True

>>> import brainpy as bp
>>> import brainmodels
>>> import matplotlib.pyplot as plt
>>>
>>> group = brainmodels.neurons.HH(2)
>>>
>>> I1 = bp.inputs.spike_input(sp_times=[500., 550., 1000, 1030, 1060, 1100, 1200], sp_lens=5, sp_sizes=5., duration=2000, )
>>> I2 = bp.inputs.spike_input(sp_times=[600.,       900, 950, 1500], sp_lens=5, sp_sizes=5., duration=2000, )
>>> I1 += bp.math.random.normal(0, 3, size=I1.shape)
>>> I2 += bp.math.random.normal(0, 3, size=I2.shape)
>>> I = bp.math.stack((I1, I2), axis=-1)
>>>
>>> runner = bp.StructRunner(group, monitors=['V'], inputs=('input', I, 'iter'))
>>> runner.run(2000.)
>>>
>>> fig, gs = bp.visualize.get_figure(1, 1, 3, 8)
>>> plt.plot(runner.mon.ts, runner.mon.V[:, 0])
>>> plt.plot(runner.mon.ts, runner.mon.V[:, 1] + 130)
>>> plt.xlim(10, 2000)
>>> plt.xticks([])
>>> plt.yticks([])
>>> plt.show()

**Model Parameters**

============= ============== ======== ====================================
**Parameter** **Init Value** **Unit** **Explanation**
------------- -------------- -------- ------------------------------------
V_th          20.            mV       the spike threshold.
E_Na          50.            mV       reversal potential of sodium.
E_K           -77.           mV       reversal potential of potassium.
E_leak        54.387         mV       reversal potential of unspecific.
g_Na          120.           msiemens conductance of sodium channel.
g_K           36.            msiemens conductance of potassium channel.
g_leak        .03            msiemens conductance of unspecific channels.
============= ============== ======== ====================================

**Model Variables**

================== ================= =========================================================
**Variables name** **Initial Value** **Explanation**
------------------ ----------------- ---------------------------------------------------------
V                        -65         Membrane potential.
m                        0.05        gating variable of the sodium ion channel.
n                        0.32        gating variable of the potassium ion channel.
h                        0.60        gating variable of the sodium ion channel.
input                     0          External and synaptic input current.
spike                    False       Flag to mark whether the neuron is spiking.
t_last_spike       -1e7               Last spike time stamp.
================== ================= =========================================================

**References**

.. [1] Hodgkin, Alan L., and Andrew F. Huxley. "A quantitative description
of membrane current and its application to conduction and excitation
in nerve." The Journal of physiology 117.4 (1952): 500.
.. [2] https://en.wikipedia.org/wiki/Hodgkin%E2%80%93Huxley_model
.. [3] Ashwin, Peter, Stephen Coombes, and Rachel Nicks. "Mathematical
frameworks for oscillatory network dynamics in neuroscience."
The Journal of Mathematical Neuroscience 6, no. 1 (2016): 1-92.
"""

[docs]  def __init__(self, size, ENa=50., gNa=120., EK=-77., gK=36., EL=-54.387, gL=0.03,
V_th=20., C=1.0, method='exp_auto', name=None):
# initialization
super(HH, self).__init__(size=size, method=method, name=name)

# parameters
self.ENa = ENa
self.EK = EK
self.EL = EL
self.gNa = gNa
self.gK = gK
self.gL = gL
self.C = C
self.V_th = V_th

# variables
self.m = bm.Variable(0.5 * bm.ones(self.num))
self.h = bm.Variable(0.6 * bm.ones(self.num))
self.n = bm.Variable(0.32 * bm.ones(self.num))

def dm(self, m, t, V):
alpha = 0.1 * (V + 40) / (1 - bm.exp(-(V + 40) / 10))
beta = 4.0 * bm.exp(-(V + 65) / 18)
dmdt = alpha * (1 - m) - beta * m
return dmdt

def dh(self, h, t, V):
alpha = 0.07 * bm.exp(-(V + 65) / 20.)
beta = 1 / (1 + bm.exp(-(V + 35) / 10))
dhdt = alpha * (1 - h) - beta * h
return dhdt

def dn(self, n, t, V):
alpha = 0.01 * (V + 55) / (1 - bm.exp(-(V + 55) / 10))
beta = 0.125 * bm.exp(-(V + 65) / 80)
dndt = alpha * (1 - n) - beta * n
return dndt

def dV(self, V, t, m, h, n, Iext):
I_Na = (self.gNa * m ** 3.0 * h) * (V - self.ENa)
I_K = (self.gK * n ** 4.0) * (V - self.EK)
I_leak = self.gL * (V - self.EL)
dVdt = (- I_Na - I_K - I_leak + Iext) / self.C
return dVdt

@property
def derivative(self):
return bp.JointEq([self.dV, self.dm, self.dh, self.dn])

[docs]  def update(self, _t, _dt):
V, m, h, n = self.integral(self.V, self.m, self.h, self.n, _t, self.input, dt=_dt)
self.spike.value = bm.logical_and(self.V < self.V_th, V >= self.V_th)
self.t_last_spike.value = bm.where(self.spike, _t, self.t_last_spike)
self.V.value = V
self.m.value = m
self.h.value = h
self.n.value = n
self.input[:] = 0.