Skip to main content

sc_neurocore_engine/neurons/biophysical/
hodgkin_huxley.rs

1// SPDX-License-Identifier: AGPL-3.0-or-later
2// Commercial license available
3// © Concepts 1996–2026 Miroslav Šotek. All rights reserved.
4// © Code 2020–2026 Miroslav Šotek. All rights reserved.
5// ORCID: 0009-0009-3560-0851
6// Contact: www.anulum.li | protoscience@anulum.li
7// SC-NeuroCore — Hodgkin-Huxley Neuron Model
8
9//! Hodgkin-Huxley 1952 conductance-based neuron dynamics.
10
11use super::safe_rate;
12
13/// Hodgkin-Huxley 1952 — 4-ODE ion channel model.
14#[derive(Clone, Debug)]
15pub struct HodgkinHuxleyNeuron {
16    pub v: f64,
17    pub m: f64,
18    pub h: f64,
19    pub n: f64,
20    pub c_m: f64,
21    pub g_na: f64,
22    pub g_k: f64,
23    pub g_l: f64,
24    pub e_na: f64,
25    pub e_k: f64,
26    pub e_l: f64,
27    pub dt: f64,
28    pub v_threshold: f64,
29}
30
31impl HodgkinHuxleyNeuron {
32    pub fn new() -> Self {
33        Self {
34            v: -65.0,
35            m: 0.05,
36            h: 0.6,
37            n: 0.32,
38            c_m: 1.0,
39            g_na: 120.0,
40            g_k: 36.0,
41            g_l: 0.3,
42            e_na: 50.0,
43            e_k: -77.0,
44            e_l: -54.4,
45            dt: 0.01,
46            v_threshold: 0.0,
47        }
48    }
49    pub fn step(&mut self, current: f64) -> i32 {
50        let v_prev = self.v;
51        let steps = (1.0 / self.dt) as usize;
52        for _ in 0..steps {
53            let am = safe_rate(0.1, 40.0, self.v, 10.0, 1.0);
54            let bm = 4.0 * (-(self.v + 65.0) / 18.0).exp();
55            let ah = 0.07 * (-(self.v + 65.0) / 20.0).exp();
56            let bh = 1.0 / (1.0 + (-(self.v + 35.0) / 10.0).exp());
57            let an = safe_rate(0.01, 55.0, self.v, 10.0, 0.1);
58            let bn = 0.125 * (-(self.v + 65.0) / 80.0).exp();
59            self.m += (am * (1.0 - self.m) - bm * self.m) * self.dt;
60            self.h += (ah * (1.0 - self.h) - bh * self.h) * self.dt;
61            self.n += (an * (1.0 - self.n) - bn * self.n) * self.dt;
62            let i_na = self.g_na * self.m.powi(3) * self.h * (self.v - self.e_na);
63            let i_k = self.g_k * self.n.powi(4) * (self.v - self.e_k);
64            let i_l = self.g_l * (self.v - self.e_l);
65            self.v += (-i_na - i_k - i_l + current) / self.c_m * self.dt;
66        }
67        if self.v >= self.v_threshold && v_prev < self.v_threshold {
68            1
69        } else {
70            0
71        }
72    }
73    pub fn reset(&mut self) {
74        self.v = -65.0;
75        self.m = 0.05;
76        self.h = 0.6;
77        self.n = 0.32;
78    }
79}
80impl Default for HodgkinHuxleyNeuron {
81    fn default() -> Self {
82        Self::new()
83    }
84}
85
86#[cfg(test)]
87mod tests {
88    use super::*;
89
90    #[test]
91    fn default_matches_constructor_state() {
92        let default = HodgkinHuxleyNeuron::default();
93        let constructed = HodgkinHuxleyNeuron::new();
94        assert_eq!(default.v, constructed.v);
95    }
96
97    #[test]
98    fn hh_fires() {
99        let mut n = HodgkinHuxleyNeuron::new();
100        let t: i32 = (0..100).map(|_| n.step(10.0)).sum();
101        assert!(t > 0);
102    }
103
104    // -- HodgkinHuxley --
105    #[test]
106    fn hh_silent_without_input() {
107        let mut n = HodgkinHuxleyNeuron::new();
108        let t: i32 = (0..200).map(|_| n.step(0.0)).sum();
109        assert_eq!(t, 0);
110    }
111    #[test]
112    fn hh_reset_clears_state() {
113        let mut n = HodgkinHuxleyNeuron::new();
114        for _ in 0..100 {
115            n.step(10.0);
116        }
117        n.reset();
118        assert!((n.v - (-65.0)).abs() < 1e-10);
119        assert!((n.m - 0.05).abs() < 1e-10);
120        assert!((n.h - 0.6).abs() < 1e-10);
121        assert!((n.n - 0.32).abs() < 1e-10);
122    }
123    #[test]
124    fn hh_extreme_input_bounded() {
125        let mut n = HodgkinHuxleyNeuron::new();
126        for _ in 0..200 {
127            n.step(1e4);
128        }
129        assert!(n.v.is_finite());
130    }
131    #[test]
132    fn hh_gates_bounded() {
133        let mut n = HodgkinHuxleyNeuron::new();
134        for _ in 0..500 {
135            n.step(10.0);
136        }
137        assert!(n.m >= 0.0 && n.m <= 1.0, "m={}", n.m);
138        assert!(n.h >= 0.0 && n.h <= 1.0, "h={}", n.h);
139        assert!(n.n >= 0.0 && n.n <= 1.0, "n={}", n.n);
140    }
141    #[test]
142    fn hh_negative_input_no_crash() {
143        let mut n = HodgkinHuxleyNeuron::new();
144        for _ in 0..200 {
145            n.step(-20.0);
146        }
147        assert!(n.v.is_finite());
148    }
149    #[test]
150    fn hh_nan_input_no_panic() {
151        let mut n = HodgkinHuxleyNeuron::new();
152        n.step(f64::NAN);
153    }
154    #[test]
155    fn hh_sodium_potassium_opposition() {
156        // Na activation drives depolarisation, K drives repolarisation
157        let mut n = HodgkinHuxleyNeuron::new();
158        for _ in 0..50 {
159            n.step(10.0);
160        }
161        // After spiking, n (K activation) should have risen
162        assert!(n.n > 0.32, "K activation n should increase during spiking");
163    }
164}