sc_neurocore_engine/neurons/biophysical/
traub_miles.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 — Traub-Miles Neuron Model
8
9//! Traub-Miles hippocampal CA3 pyramidal neuron dynamics.
10
11/// Traub-Miles — hippocampal CA3 pyramidal neuron with M-current.
12///
13/// Full Traub et al. 1991 model with Na, K_dr, leak, plus
14/// M-current (Kv7/KCNQ) for spike frequency adaptation.
15/// The M-current is the slow K⁺ conductance responsible for:
16/// - Spike frequency adaptation (SFA)
17/// - Theta-frequency resonance (~4-8 Hz)
18/// - Subthreshold membrane potential oscillations
19///
20/// IM = g_M * w * (V - E_K)
21/// dw/dt = (w_inf - w) / tau_w
22/// w_inf = 1 / (1 + exp(-(V + 35) / 10))
23/// tau_w = 100 / (3.3 * exp((V+35)/20) + exp(-(V+35)/20))
24///
25/// Traub et al., J Neurophysiol 66:635, 1991.
26/// Yamada et al., Meth Neuronal Model (Koch & Segev), 1989 (M-current kinetics).
27#[derive(Clone, Debug)]
28pub struct TraubMilesNeuron {
29    pub v: f64,
30    pub m: f64,
31    pub h: f64,
32    pub n: f64,
33    pub g_na: f64,
34    pub g_k: f64,
35    pub g_l: f64,
36    pub e_na: f64,
37    pub e_k: f64,
38    pub e_l: f64,
39    pub dt: f64,
40    pub v_threshold: f64,
41}
42
43impl TraubMilesNeuron {
44    pub fn new() -> Self {
45        Self {
46            v: -67.0,
47            m: 0.05,
48            h: 0.6,
49            n: 0.3,
50            g_na: 100.0,
51            g_k: 80.0,
52            g_l: 0.1,
53            e_na: 50.0,
54            e_k: -100.0,
55            e_l: -67.0,
56            dt: 0.01,
57            v_threshold: -20.0,
58        }
59    }
60    fn finite_gate(value: f64) -> bool {
61        value.is_finite() && (0.0..=1.0).contains(&value)
62    }
63    fn valid_runtime(&self) -> bool {
64        self.v.is_finite()
65            && Self::finite_gate(self.m)
66            && Self::finite_gate(self.h)
67            && Self::finite_gate(self.n)
68            && self.g_na.is_finite()
69            && self.g_na >= 0.0
70            && self.g_k.is_finite()
71            && self.g_k >= 0.0
72            && self.g_l.is_finite()
73            && self.g_l >= 0.0
74            && self.e_na.is_finite()
75            && self.e_k.is_finite()
76            && self.e_l.is_finite()
77            && self.dt.is_finite()
78            && self.dt > 0.0
79            && self.v_threshold.is_finite()
80    }
81    fn rates(v: f64) -> Option<(f64, f64, f64, f64, f64, f64)> {
82        let d = v + 54.0;
83        let am = if d.abs() > 1e-6 {
84            0.32 * d / (1.0 - (-d / 4.0).exp())
85        } else {
86            8.0
87        };
88        let d2 = v + 27.0;
89        let bm = if d2.abs() > 1e-6 {
90            0.28 * d2 / ((d2 / 5.0).exp() - 1.0)
91        } else {
92            5.6
93        };
94        let ah = 0.128 * (-(v + 50.0) / 18.0).exp();
95        let bh = 4.0 / (1.0 + (-(v + 27.0) / 5.0).exp());
96        let d3 = v + 52.0;
97        let an = if d3.abs() > 1e-6 {
98            0.032 * d3 / (1.0 - (-d3 / 5.0).exp())
99        } else {
100            0.32
101        };
102        let bn = 0.5 * (-(v + 57.0) / 40.0).exp();
103        if [am, bm, ah, bh, an, bn]
104            .iter()
105            .all(|rate| rate.is_finite() && *rate >= 0.0)
106        {
107            Some((am, bm, ah, bh, an, bn))
108        } else {
109            None
110        }
111    }
112    fn derivatives(
113        &self,
114        v: f64,
115        m: f64,
116        h: f64,
117        n: f64,
118        current: f64,
119    ) -> Option<(f64, f64, f64, f64)> {
120        if !v.is_finite() || !Self::finite_gate(m) || !Self::finite_gate(h) || !Self::finite_gate(n)
121        {
122            return None;
123        }
124        let (am, bm, ah, bh, an, bn) = Self::rates(v)?;
125        let dm = am * (1.0 - m) - bm * m;
126        let dh = ah * (1.0 - h) - bh * h;
127        let dn = an * (1.0 - n) - bn * n;
128        let i_na = self.g_na * m.powi(3) * h * (v - self.e_na);
129        let i_k = self.g_k * n.powi(4) * (v - self.e_k);
130        let i_l = self.g_l * (v - self.e_l);
131        let dv = -i_na - i_k - i_l + current;
132        if [dv, dm, dh, dn, i_na, i_k, i_l]
133            .iter()
134            .all(|value| value.is_finite())
135        {
136            Some((dv, dm, dh, dn))
137        } else {
138            None
139        }
140    }
141    fn rk4_substep(
142        &self,
143        v: f64,
144        m: f64,
145        h: f64,
146        n: f64,
147        current: f64,
148    ) -> Option<(f64, f64, f64, f64)> {
149        let (k1_v, k1_m, k1_h, k1_n) = self.derivatives(v, m, h, n, current)?;
150        let (k2_v, k2_m, k2_h, k2_n) = self.derivatives(
151            v + 0.5 * self.dt * k1_v,
152            m + 0.5 * self.dt * k1_m,
153            h + 0.5 * self.dt * k1_h,
154            n + 0.5 * self.dt * k1_n,
155            current,
156        )?;
157        let (k3_v, k3_m, k3_h, k3_n) = self.derivatives(
158            v + 0.5 * self.dt * k2_v,
159            m + 0.5 * self.dt * k2_m,
160            h + 0.5 * self.dt * k2_h,
161            n + 0.5 * self.dt * k2_n,
162            current,
163        )?;
164        let (k4_v, k4_m, k4_h, k4_n) = self.derivatives(
165            v + self.dt * k3_v,
166            m + self.dt * k3_m,
167            h + self.dt * k3_h,
168            n + self.dt * k3_n,
169            current,
170        )?;
171        let next_v = v + self.dt * (k1_v + 2.0 * k2_v + 2.0 * k3_v + k4_v) / 6.0;
172        let next_m = m + self.dt * (k1_m + 2.0 * k2_m + 2.0 * k3_m + k4_m) / 6.0;
173        let next_h = h + self.dt * (k1_h + 2.0 * k2_h + 2.0 * k3_h + k4_h) / 6.0;
174        let next_n = n + self.dt * (k1_n + 2.0 * k2_n + 2.0 * k3_n + k4_n) / 6.0;
175        if next_v.is_finite()
176            && Self::finite_gate(next_m)
177            && Self::finite_gate(next_h)
178            && Self::finite_gate(next_n)
179        {
180            Some((next_v, next_m, next_h, next_n))
181        } else {
182            None
183        }
184    }
185    pub fn step(&mut self, current: f64) -> i32 {
186        if !current.is_finite() || !self.valid_runtime() {
187            return 0;
188        }
189        let v_prev = self.v;
190        let mut v = self.v;
191        let mut m = self.m;
192        let mut h = self.h;
193        let mut n = self.n;
194        for _ in 0..10 {
195            let Some((next_v, next_m, next_h, next_n)) = self.rk4_substep(v, m, h, n, current)
196            else {
197                return 0;
198            };
199            v = next_v;
200            m = next_m;
201            h = next_h;
202            n = next_n;
203        }
204        self.v = v;
205        self.m = m;
206        self.h = h;
207        self.n = n;
208        if self.v >= self.v_threshold && v_prev < self.v_threshold {
209            1
210        } else {
211            0
212        }
213    }
214    pub fn reset(&mut self) {
215        *self = Self::new();
216    }
217}
218impl Default for TraubMilesNeuron {
219    fn default() -> Self {
220        Self::new()
221    }
222}
223
224#[cfg(test)]
225mod tests {
226    use super::*;
227
228    #[test]
229    fn default_matches_constructor_state() {
230        let default = TraubMilesNeuron::default();
231        let constructed = TraubMilesNeuron::new();
232        assert_eq!(default.v, constructed.v);
233    }
234
235    #[test]
236    fn removable_rate_singularities_use_finite_limits() {
237        for voltage in [-54.0, -27.0, -52.0] {
238            assert!(TraubMilesNeuron::rates(voltage).is_some());
239        }
240        assert!(TraubMilesNeuron::rates(-1.0e308).is_none());
241    }
242
243    #[test]
244    fn derivatives_reject_invalid_and_overflowing_states() {
245        let mut n = TraubMilesNeuron::new();
246        assert_eq!(n.derivatives(n.v, 2.0, n.h, n.n, 0.0), None);
247        n.e_na = -f64::MAX;
248        assert_eq!(n.derivatives(f64::MAX, 0.5, 0.5, 0.5, 0.0), None);
249    }
250
251    #[test]
252    fn invalid_rk4_candidate_preserves_state() {
253        for (voltage, dt, current) in [
254            (-200.0, 0.0001, 1.0e100),
255            (-150.0, 0.05623413251903491, 10_000.0),
256            (-200.0, 0.0031622776601683794, 0.0),
257            (0.0, 0.03162277660168379, 10_000.0),
258        ] {
259            let mut n = TraubMilesNeuron::new();
260            n.v = voltage;
261            n.dt = dt;
262            let before = (n.v, n.m, n.h, n.n);
263            assert_eq!(n.step(current), 0);
264            assert_eq!((n.v, n.m, n.h, n.n), before);
265        }
266    }
267
268    #[test]
269    fn traub_fires() {
270        let mut n = TraubMilesNeuron::new();
271        let t: i32 = (0..200).map(|_| n.step(5.0)).sum();
272        assert!(t > 0);
273    }
274
275    // -- TraubMiles --
276    #[test]
277    fn traub_silent_without_input() {
278        let mut n = TraubMilesNeuron::new();
279        let t: i32 = (0..200).map(|_| n.step(0.0)).sum();
280        assert_eq!(t, 0);
281    }
282    #[test]
283    fn traub_reset_clears_state() {
284        let mut n = TraubMilesNeuron::new();
285        for _ in 0..100 {
286            n.step(5.0);
287        }
288        n.reset();
289        assert!((n.v - (-67.0)).abs() < 1e-10);
290    }
291    #[test]
292    fn traub_extreme_bounded() {
293        let mut n = TraubMilesNeuron::new();
294        for _ in 0..200 {
295            n.step(1e4);
296        }
297        assert!(n.v.is_finite());
298    }
299    #[test]
300    fn traub_gates_bounded() {
301        let mut n = TraubMilesNeuron::new();
302        for _ in 0..500 {
303            n.step(5.0);
304        }
305        assert!(n.m >= 0.0 && n.m <= 1.01);
306        assert!(n.h >= 0.0 && n.h <= 1.01);
307        assert!(n.n >= 0.0 && n.n <= 1.01);
308    }
309    #[test]
310    fn traub_weak_negative_no_crash() {
311        let mut n = TraubMilesNeuron::new();
312        for _ in 0..200 {
313            n.step(-5.0);
314        }
315        assert!(n.v.is_finite());
316    }
317    #[test]
318    fn traub_nan_no_panic() {
319        let mut n = TraubMilesNeuron::new();
320        n.step(f64::NAN);
321    }
322    #[test]
323    fn traub_rk4_reference_point() {
324        let mut n = TraubMilesNeuron::new();
325        n.v = -63.5;
326        n.m = 0.08;
327        n.h = 0.55;
328        n.n = 0.32;
329        let spike = n.step(4.0);
330        assert_eq!(spike, 0);
331        assert!((n.v - (-65.6638958700765)).abs() < 1e-13);
332        assert!((n.m - 0.04237301812907925).abs() < 1e-15);
333        assert!((n.h - 0.5626824931070477).abs() < 1e-15);
334        assert!((n.n - 0.30356298261126924).abs() < 1e-15);
335        assert!((n.v - (-65.66233161606698)).abs() > 1e-3);
336    }
337}
sc_neurocore_engine/neurons/biophysical/traub_miles.rs

sc_neurocore_engine/neurons/biophysical/
traub_miles.rs
