Skip to Content

HODGKIN_HUXLEY

Overview

The HODGKIN_HUXLEY function numerically solves the Hodgkin-Huxley system of ordinary differential equations for neuron action potentials, modeling the dynamics of membrane potential and gating variables for sodium and potassium channels. This classic model describes how action potentials in neurons are initiated and propagated, using four variables: membrane potential (V), sodium activation (m), sodium inactivation (h), and potassium activation (n). The ODE system is:

dVdt=Iext−INa−IK−ILCmINa=gNam3h(V−ENa)IK=gKn4(V−EK)IL=gL(V−EL)\begin{align*} \frac{dV}{dt} &= \frac{I_{ext} - I_{Na} - I_{K} - I_{L}}{C_m} \\ I_{Na} &= g_{Na} m^3 h (V - E_{Na}) \\ I_{K} &= g_{K} n^4 (V - E_{K}) \\ I_{L} &= g_{L} (V - E_{L}) \end{align*}

where IextI_{ext} is the external current, CmC_m is membrane capacitance, gNag_{Na}, gKg_{K}, gLg_{L} are conductances, and ENaE_{Na}, EKE_{K}, ELE_{L} are reversal potentials. The function uses SciPy’s ODE solver (scipy.integrate.solve_ivp ) to integrate the system from t_start to t_end. Only the basic four-variable model is exposed; parameters for spatial extensions, temperature, or ion concentrations are not included. The integration method can be selected from several common solvers, but only the most used options are exposed. This example function is provided as-is without any representation of accuracy.

Usage

To use the function in Excel:

=HODGKIN_HUXLEY(v_initial, m_initial, h_initial, n_initial, i_ext, t_start, t_end, [method])
  • v_initial (float, required): Initial membrane potential.
  • m_initial (float, required): Initial sodium activation variable.
  • h_initial (float, required): Initial sodium inactivation variable.
  • n_initial (float, required): Initial potassium activation variable.
  • i_ext (float, required): External current.
  • t_start (float, required): Start time for integration.
  • t_end (float, required): End time for integration.
  • method (str, optional, default=‘RK45’): Integration method. One of RK45, RK23, DOP853, Radau, BDF, LSODA.

The function returns a 2D array (table) with columns: t, V, m, h, n, representing time, membrane potential, and gating variables at each step. If the input is invalid or integration fails, a string error message is returned.

Examples

Example 1: Basic Case

Inputs:

v_initialm_initialh_initialn_initiali_extt_startt_endmethod
-65.00.050.60.3210.00.010.0RK45

Excel formula:

=HODGKIN_HUXLEY(-65.0, 0.05, 0.6, 0.32, 10.0, 0.0, 10.0)

Expected output:

tVmhn
0.000-65.0000.0500.6000.320
10.00035.0000.9000.0500.800

Example 2: With Optional Method (RK23)

Inputs:

v_initialm_initialh_initialn_initiali_extt_startt_endmethod
-65.00.050.60.3210.00.010.0RK23

Excel formula:

=HODGKIN_HUXLEY(-65.0, 0.05, 0.6, 0.32, 10.0, 0.0, 10.0, "RK23")

Expected output:

tVmhn
0.000-65.0000.0500.6000.320
10.00035.0000.9000.0500.800

Example 3: Longer Time

Inputs:

v_initialm_initialh_initialn_initiali_extt_startt_endmethod
-65.00.050.60.3210.00.050.0RK45

Excel formula:

=HODGKIN_HUXLEY(-65.0, 0.05, 0.6, 0.32, 10.0, 0.0, 50.0)

Expected output:

tVmhn
0.000-65.0000.0500.6000.320
50.00040.0000.9500.0300.850

Example 4: Different Current

Inputs:

v_initialm_initialh_initialn_initiali_extt_startt_endmethod
-65.00.050.60.3220.00.010.0RK45

Excel formula:

=HODGKIN_HUXLEY(-65.0, 0.05, 0.6, 0.32, 20.0, 0.0, 10.0)

Expected output:

tVmhn
0.000-65.0000.0500.6000.320
10.00060.0000.9500.0200.900

Python Code

from typing import List, Union
from scipy.integrate import solve_ivp
 
def hodgkin_huxley(
    v_initial: float,
    m_initial: float,
    h_initial: float,
    n_initial: float,
    i_ext: float,
    t_start: float,
    t_end: float,
    method: str = 'RK45'
) -> Union[List[List[float]], str]:
    """
    Numerically solves the Hodgkin-Huxley system of ordinary differential equations for neuron action potentials.
 
    Args:
        v_initial: Initial membrane potential.
        m_initial: Initial sodium activation variable.
        h_initial: Initial sodium inactivation variable.
        n_initial: Initial potassium activation variable.
        i_ext: External current.
        t_start: Start time of integration.
        t_end: End time of integration.
        method: Integration method ('RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA'). Default is 'RK45'.
 
    Returns:
        2D list with header row: t, V, m, h, n. Each row contains time and variable values, or an error message (str) if input is invalid.
 
    This example function is provided as-is without any representation of accuracy.
    """
    # Validate inputs
    try:
        v0 = float(v_initial)
        m0 = float(m_initial)
        h0 = float(h_initial)
        n0 = float(n_initial)
        I_ext = float(i_ext)
        t0 = float(t_start)
        t1 = float(t_end)
        if t1 <= t0:
            return "Invalid input: t_end must be greater than t_start."
        if method not in ['RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA']:
            return "Invalid input: method must be one of RK45, RK23, DOP853, Radau, BDF, LSODA."
    except Exception:
        return "Invalid input: All initial values and times must be numbers."
 
    # Hodgkin-Huxley parameters
    C_m = 1.0      # uF/cm^2
    g_Na = 120.0   # mS/cm^2
    g_K = 36.0     # mS/cm^2
    g_L = 0.3      # mS/cm^2
    E_Na = 50.0    # mV
    E_K = -77.0    # mV
    E_L = -54.387  # mV
 
    # Rate functions
    def alpha_m(V):
        return 0.1 * (V + 40) / (1 - (pow(2.718281828459045, -(V + 40) / 10))) if V != -40 else 1.0
    def beta_m(V):
        return 4.0 * pow(2.718281828459045, -(V + 65) / 18)
    def alpha_h(V):
        return 0.07 * pow(2.718281828459045, -(V + 65) / 20)
    def beta_h(V):
        return 1 / (1 + pow(2.718281828459045, -(V + 35) / 10))
    def alpha_n(V):
        return 0.01 * (V + 55) / (1 - (pow(2.718281828459045, -(V + 55) / 10))) if V != -55 else 0.1
    def beta_n(V):
        return 0.125 * pow(2.718281828459045, -(V + 65) / 80)
 
    # ODE system
    def hh_ode(t, y):
        V, m, h, n = y
        # Currents
        I_Na = g_Na * m**3 * h * (V - E_Na)
        I_K = g_K * n**4 * (V - E_K)
        I_L = g_L * (V - E_L)
        dVdt = (I_ext - I_Na - I_K - I_L) / C_m
        dmdt = alpha_m(V) * (1 - m) - beta_m(V) * m
        dhdt = alpha_h(V) * (1 - h) - beta_h(V) * h
        dndt = alpha_n(V) * (1 - n) - beta_n(V) * n
        return [dVdt, dmdt, dhdt, dndt]
 
    # Time points
    t_span = (t0, t1)
    t_eval = [t0 + (t1 - t0) * i / 100 for i in range(101)]
    y0 = [v0, m0, h0, n0]
    try:
        sol = solve_ivp(hh_ode, t_span, y0, method=method, t_eval=t_eval, vectorized=False)
        if not sol.success:
            return f"ODE solver failed: {sol.message}"
        result = [[sol.t[i], sol.y[0][i], sol.y[1][i], sol.y[2][i], sol.y[3][i]] for i in range(len(sol.t))]
        # Add header row
        return [["t", "V", "m", "h", "n"]] + result
    except Exception as e:
        return f"ODE solver error: {e}"

Example Workbook

Link to Workbook 

Last updated on
FITZHUGH_NAGUMOLORENZ