HODGKIN_HUXLEY

Overview

The HODGKIN_HUXLEY function numerically solves the Hodgkin–Huxley model, a landmark system of ordinary differential equations that describes how action potentials in neurons are initiated and propagated. Developed by Alan Hodgkin and Andrew Huxley in 1952—work that earned them the 1963 Nobel Prize in Physiology or Medicine—this model quantitatively explains the ionic mechanisms underlying nerve impulse transmission in excitable cells.

The model represents the neuronal membrane as an electrical circuit. The lipid bilayer acts as a capacitor (C_m), while voltage-gated ion channels for sodium (Na⁺) and potassium (K⁺), along with a passive leak current, are modeled as variable conductances. The membrane potential V evolves according to:

C_m \frac{dV}{dt} = I_{ext} - \bar{g}_{Na} m^3 h (V - E_{Na}) - \bar{g}_K n^4 (V - E_K) - g_L (V - E_L)

where I_{ext} is the injected current, \bar{g}_{Na}, \bar{g}_K, and g_L are the maximum conductances, and E_{Na}, E_K, and E_L are the reversal potentials for sodium, potassium, and leak channels respectively.

The gating variables m, h, and n represent the probability of ion channel subunits being in the open state. Each follows first-order kinetics governed by voltage-dependent rate constants \alpha and \beta:

\frac{dm}{dt} = \alpha_m(V)(1-m) - \beta_m(V)m

The sodium current depends on both an activation gate (m^3) and an inactivation gate (h), while the potassium current depends on an activation gate raised to the fourth power (n^4), reflecting the multi-subunit structure of these channels.

This implementation uses SciPy’s solve_ivp function to numerically integrate the system, with support for multiple integration methods including RK45 (default), RK23, DOP853, Radau, BDF, and LSODA. The function returns the time evolution of the membrane potential and all three gating variables.

This example function is provided as-is without any representation of accuracy.

Excel Usage

=HODGKIN_HUXLEY(v_initial, m_initial, h_initial, n_initial, i_ext, t_start, t_end, timesteps, solve_ivp_method)

v_initial (float, required): Initial membrane potential (mV)
m_initial (float, required): Initial sodium activation variable (0 to 1)
h_initial (float, required): Initial sodium inactivation variable (0 to 1)
n_initial (float, required): Initial potassium activation variable (0 to 1)
i_ext (float, required): External current (uA/cm^2)
t_start (float, required): Start time of integration (ms)
t_end (float, required): End time of integration (ms)
timesteps (int, optional, default: 10): Number of timesteps to solve for
solve_ivp_method (str, optional, default: “RK45”): Integration method for the ODE solver

Returns (list[list]): 2D list with header [t, V, m, h, n], or error message string.

Examples

Example 1: Demo case 1

Inputs:

v_initial	m_initial	h_initial	n_initial	i_ext	t_start	t_end
-65	0.05	0.6	0.32	10	0	10

Excel formula:

=HODGKIN_HUXLEY(-65, 0.05, 0.6, 0.32, 10, 0, 10)

Expected output:

t	V	m	h	n
0	-65	0.05	0.6	0.32
1.111	-54.96	0.1184	0.5745	0.3353
2.222	39.53	0.9394	0.3271	0.5346
3.333	-11.55	0.9682	0.1109	0.7602
4.444	-68.13	0.296	0.09	0.7419
5.556	-74.59	0.01608	0.1933	0.6451
6.667	-73.17	0.01892	0.2805	0.5661
7.778	-71.3	0.02368	0.3488	0.5043
8.889	-69.09	0.03085	0.3999	0.4577
10	-66.75	0.04077	0.4355	0.4248

Example 2: Demo case 2

Inputs:

v_initial	m_initial	h_initial	n_initial	i_ext	t_start	t_end	timesteps	solve_ivp_method
-65	0.05	0.6	0.32	10	0	10	10	RK23

Excel formula:

=HODGKIN_HUXLEY(-65, 0.05, 0.6, 0.32, 10, 0, 10, 10, "RK23")

Expected output:

t	V	m	h	n
0	-65	0.05	0.6	0.32
1.111	-54.96	0.1182	0.5745	0.3353
2.222	39.68	0.9384	0.3278	0.5335
3.333	-11.41	0.9684	0.1111	0.76
4.444	-68.14	0.3	0.08978	0.7421
5.556	-74.6	0.01597	0.1931	0.6453
6.667	-73.19	0.01887	0.2804	0.5663
7.778	-71.32	0.02363	0.3487	0.5044
8.889	-69.12	0.03065	0.3998	0.4578
10	-66.77	0.04065	0.4355	0.4248

Example 3: Demo case 3

Inputs:

v_initial	m_initial	h_initial	n_initial	i_ext	t_start	t_end	timesteps
-65	0.05	0.6	0.32	10	0	50	10

Excel formula:

=HODGKIN_HUXLEY(-65, 0.05, 0.6, 0.32, 10, 0, 50, 10)

Expected output:

t	V	m	h	n
0	-65	0.05	0.6	0.32
5.556	-74.59	0.01608	0.1933	0.6451
11.11	-64.46	0.05335	0.4575	0.4036
16.67	-29.86	0.3665	0.3555	0.4472
22.22	-71.44	0.02328	0.3288	0.5082
27.78	-60.48	0.08412	0.4594	0.3891
33.33	-42.81	0.7552	0.06882	0.7495
38.89	-67.31	0.03813	0.4151	0.4315
44.44	-56.15	0.1294	0.4353	0.4004
50	-73.81	0.01751	0.2268	0.5963

Example 4: Demo case 4

Inputs:

v_initial	m_initial	h_initial	n_initial	i_ext	t_start	t_end	timesteps
-65	0.05	0.6	0.32	15	0	10	10

Excel formula:

=HODGKIN_HUXLEY(-65, 0.05, 0.6, 0.32, 15, 0, 10, 10)

Expected output:

t	V	m	h	n
0	-65	0.05	0.6	0.32
1.111	-46.99	0.1805	0.5539	0.3465
2.222	24.03	0.9941	0.2209	0.6663
3.333	-29.09	0.893	0.08347	0.7727
4.444	-74.5	0.04235	0.1247	0.7081
5.556	-73.34	0.01851	0.222	0.6177
6.667	-71.36	0.02343	0.2984	0.5466
7.778	-68.91	0.03138	0.3555	0.4926
8.889	-66.21	0.04306	0.3951	0.4543
10	-63.52	0.05884	0.419	0.4296

Python Code

import math
from scipy.integrate import solve_ivp as scipy_solve_ivp
import numpy as np

def hodgkin_huxley(v_initial, m_initial, h_initial, n_initial, i_ext, t_start, t_end, timesteps=10, solve_ivp_method='RK45'):
    """
    Numerically solves the Hodgkin-Huxley system of ordinary differential equations for neuron action potentials.

    See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

    This example function is provided as-is without any representation of accuracy.

    Args:
        v_initial (float): Initial membrane potential (mV)
        m_initial (float): Initial sodium activation variable (0 to 1)
        h_initial (float): Initial sodium inactivation variable (0 to 1)
        n_initial (float): Initial potassium activation variable (0 to 1)
        i_ext (float): External current (uA/cm^2)
        t_start (float): Start time of integration (ms)
        t_end (float): End time of integration (ms)
        timesteps (int, optional): Number of timesteps to solve for Default is 10.
        solve_ivp_method (str, optional): Integration method for the ODE solver Valid options: RK45, RK23, DOP853, Radau, BDF, LSODA. Default is 'RK45'.

    Returns:
        list[list]: 2D list with header [t, V, m, h, n], or error message string.
    """
    # 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)
        ntp = int(timesteps)
    except Exception:
        return "Invalid input: All initial values, parameters, and timesteps must be numbers."
    if not (0 <= m0 <= 1):
        return "Invalid input: m_initial must be between 0 and 1."
    if not (0 <= h0 <= 1):
        return "Invalid input: h_initial must be between 0 and 1."
    if not (0 <= n0 <= 1):
        return "Invalid input: n_initial must be between 0 and 1."
    if t1 <= t0:
        return "Invalid input: t_end must be greater than t_start."
    if ntp <= 0:
        return "Invalid input: timesteps must be a positive integer."
    valid_methods = ['RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA']
    if solve_ivp_method not in valid_methods:
        return f"Invalid input: solve_ivp_method must be one of {', '.join(valid_methods)}."

    # 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 with singularity handling using L'Hopital's rule limits
    def alpha_m(V):
        x = (V + 40) / 10
        if abs(x) < 1e-7:
            return 1.0  # Limit as V -> -40
        return 0.1 * (V + 40) / (1 - math.exp(-x))

    def beta_m(V):
        return 4.0 * math.exp(-(V + 65) / 18)

    def alpha_h(V):
        return 0.07 * math.exp(-(V + 65) / 20)

    def beta_h(V):
        return 1.0 / (1.0 + math.exp(-(V + 35) / 10))

    def alpha_n(V):
        x = (V + 55) / 10
        if abs(x) < 1e-7:
            return 0.1  # Limit as V -> -55
        return 0.01 * (V + 55) / (1 - math.exp(-x))

    def beta_n(V):
        return 0.125 * math.exp(-(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_eval = np.linspace(t0, t1, ntp)
    y0 = [v0, m0, h0, n0]
    try:
        sol = scipy_solve_ivp(hh_ode, [t0, t1], y0, method=solve_ivp_method, t_eval=t_eval, vectorized=False)
        if not sol.success:
            return f"ODE solver failed: {sol.message}"
        result = [
            [float(sol.t[i]), float(sol.y[0][i]), float(sol.y[1][i]), float(sol.y[2][i]), float(sol.y[3][i])]
            for i in range(len(sol.t))
        ]
        return [["t", "V", "m", "h", "n"]] + result
    except Exception as e:
        return f"ODE solver error: {e}"

Overview

Excel Usage

Examples

Python Code

Online Calculator