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:

\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 $I_{ext}$ is the external current, $C_m$ is membrane capacitance, $g_{Na}$ , $g_{K}$ , $g_{L}$ are conductances, and $E_{Na}$ , $E_{K}$ , $E_{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, [timesteps], [solve_ivp_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.
timesteps (int, optional, default=10): Number of timesteps to solve for.
solve_ivp_method (string (enum), optional, default=‘RK45’): Integration method. Valid options: 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_initial	m_initial	h_initial	n_initial	i_ext	t_start	t_end	timesteps
-65.0	0.05	0.6	0.32	10.0	0.0	10.0	10

Excel formula:


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

Expected output:

t	V	m	h	n
0.000	-65.000	0.050	0.600	0.320
1.111	-54.960	0.118	0.575	0.335
2.222	39.530	0.939	0.327	0.535
3.333	-11.550	0.968	0.111	0.760
4.444	-68.130	0.296	0.090	0.742
5.556	-74.590	0.016	0.193	0.645
6.667	-73.170	0.019	0.281	0.566
7.778	-71.300	0.024	0.349	0.504
8.889	-69.090	0.031	0.400	0.458
10.000	-66.750	0.041	0.436	0.425

Example 2: With Optional Method (RK23)

Inputs:

v_initial	m_initial	h_initial	n_initial	i_ext	t_start	t_end	timesteps	method
-65.0	0.05	0.6	0.32	10.0	0.0	10.0	10	RK23

Excel formula:


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

Expected output:

t	V	m	h	n
0.000	-65.000	0.050	0.600	0.320
1.111	-54.960	0.118	0.575	0.335
2.222	39.680	0.938	0.328	0.534
3.333	-11.410	0.968	0.111	0.760
4.444	-68.140	0.300	0.090	0.742
5.556	-74.600	0.016	0.193	0.645
6.667	-73.190	0.019	0.280	0.566
7.778	-71.320	0.024	0.349	0.504
8.889	-69.120	0.031	0.400	0.458
10.000	-66.770	0.041	0.436	0.425

Example 3: Longer Time Span

Inputs:

v_initial	m_initial	h_initial	n_initial	i_ext	t_start	t_end	timesteps
-65.0	0.05	0.6	0.32	10.0	0.0	50.0	10

Excel formula:


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

Expected output:

t	V	m	h	n
0.000	-65.000	0.050	0.600	0.320
5.556	-74.590	0.016	0.193	0.645
11.111	-64.460	0.053	0.458	0.404
16.667	-29.860	0.367	0.355	0.447
22.222	-71.440	0.023	0.329	0.508
27.778	-60.480	0.084	0.459	0.389
33.333	-42.810	0.755	0.069	0.750
38.889	-67.310	0.038	0.415	0.432
44.444	-56.150	0.129	0.435	0.400
50.000	-73.810	0.018	0.227	0.596

Example 4: Different Current

Inputs:

v_initial	m_initial	h_initial	n_initial	i_ext	t_start	t_end	timesteps
-65.0	0.05	0.6	0.32	15.0	0.0	10.0	10

Excel formula:


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

Expected output:

t	V	m	h	n
0.000	-65.000	0.050	0.600	0.320
1.111	-46.990	0.181	0.554	0.347
2.222	24.030	0.994	0.221	0.666
3.333	-29.090	0.893	0.083	0.773
4.444	-74.500	0.042	0.125	0.708
5.556	-73.340	0.019	0.222	0.618
6.667	-71.360	0.023	0.298	0.547
7.778	-68.910	0.031	0.356	0.493
8.889	-66.210	0.043	0.395	0.454
10.000	-63.520	0.059	0.419	0.430

Python Code


import math
from scipy.integrate import 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.
 
    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.
        timesteps: Number of timesteps to solve for. Default is 10.
        solve_ivp_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)
        ntp = int(timesteps)
    except Exception:
        return "Invalid input: All initial values, parameters, and timesteps must be numbers."
    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."
    if solve_ivp_method not in ['RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA']:
        return "Invalid input: solve_ivp_method must be one of RK45, RK23, DOP853, Radau, BDF, LSODA."
 
    # 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):
        exp_val = math.exp(-(V + 40) / 10)
        if abs(1 - exp_val) < 1e-10:
            return 1.0
        return 0.1 * (V + 40) / (1 - exp_val)
    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 / (1 + math.exp(-(V + 35) / 10))
    def alpha_n(V):
        exp_val = math.exp(-(V + 55) / 10)
        if abs(1 - exp_val) < 1e-10:
            return 0.1
        return 0.01 * (V + 55) / (1 - exp_val)
    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 = 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))]
        # 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