sc_neurocore_engine/neurons/biophysical/
golomb_fs.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 — Golomb Fast-Spiking Neuron Model
8
9//! Golomb Kv3-enabled fast-spiking interneuron dynamics.
10
11/// Golomb fast-spiking interneuron with Kv3. Golomb et al. 2007.
12#[derive(Clone, Debug)]
13pub struct GolombFSNeuron {
14    pub v: f64,
15    pub h: f64,
16    pub n: f64,
17    pub p: f64,
18    pub g_na: f64,
19    pub g_k: f64,
20    pub g_kv3: f64,
21    pub g_l: f64,
22    pub e_na: f64,
23    pub e_k: f64,
24    pub e_l: f64,
25    pub dt: f64,
26    pub v_threshold: f64,
27}
28
29impl GolombFSNeuron {
30    pub fn new() -> Self {
31        Self {
32            v: -65.0,
33            h: 0.9,
34            n: 0.1,
35            p: 0.0,
36            g_na: 112.5,
37            g_k: 225.0,
38            g_kv3: 150.0,
39            g_l: 0.25,
40            e_na: 50.0,
41            e_k: -90.0,
42            e_l: -70.0,
43            dt: 0.01,
44            v_threshold: -20.0,
45        }
46    }
47    /// Return `[dV, dh, dn, dp]` of the four-state Golomb-FS system at one
48    /// consistent state. The capacitance is unit-normalised (the membrane time
49    /// constant is folded into the conductances), matching the Python reference at
50    /// its default `c_m = 1`.
51    fn derivatives(&self, v: f64, h: f64, n: f64, p: f64, current: f64) -> [f64; 4] {
52        let m_inf = 1.0 / (1.0 + (-(v + 24.0) / 11.5).exp());
53        let h_inf = 1.0 / (1.0 + ((v + 58.3) / 6.7).exp());
54        let n_inf = 1.0 / (1.0 + (-(v + 12.4) / 6.8).exp());
55        let p_inf = 1.0 / (1.0 + (-(v + 3.0) / 8.0).exp());
56        let tau_h = 0.5 + 14.0 / (1.0 + ((v + 60.0) / 12.0).exp());
57        let tau_n = 0.087 + 11.4 / (1.0 + ((v + 14.6) / 8.6).exp());
58        let tau_p = 0.1 + 4.0 / (1.0 + ((v + 25.0) / 10.0).exp());
59        let dh = (h_inf - h) / tau_h;
60        let dn = (n_inf - n) / tau_n;
61        let dp = (p_inf - p) / tau_p;
62        let i_na = self.g_na * m_inf * m_inf * m_inf * h * (v - self.e_na);
63        let i_k = self.g_k * n * n * n * n * (v - self.e_k);
64        let i_kv3 = self.g_kv3 * p * p * (v - self.e_k);
65        let i_l = self.g_l * (v - self.e_l);
66        let dv = -i_na - i_k - i_kv3 - i_l + current;
67        [dv, dh, dn, dp]
68    }
69
70    /// Return one classical RK4 increment of `[V, h, n, p]`, holding `current`
71    /// constant across the four stages.
72    fn rk4_substep(&self, s: [f64; 4], current: f64) -> [f64; 4] {
73        let dt = self.dt;
74        let k1 = self.derivatives(s[0], s[1], s[2], s[3], current);
75        let k2 = self.derivatives(
76            s[0] + 0.5 * dt * k1[0],
77            s[1] + 0.5 * dt * k1[1],
78            s[2] + 0.5 * dt * k1[2],
79            s[3] + 0.5 * dt * k1[3],
80            current,
81        );
82        let k3 = self.derivatives(
83            s[0] + 0.5 * dt * k2[0],
84            s[1] + 0.5 * dt * k2[1],
85            s[2] + 0.5 * dt * k2[2],
86            s[3] + 0.5 * dt * k2[3],
87            current,
88        );
89        let k4 = self.derivatives(
90            s[0] + dt * k3[0],
91            s[1] + dt * k3[1],
92            s[2] + dt * k3[2],
93            s[3] + dt * k3[3],
94            current,
95        );
96        let mut out = [0.0_f64; 4];
97        for i in 0..4 {
98            out[i] = s[i] + dt * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]) / 6.0;
99        }
100        out
101    }
102
103    pub fn step(&mut self, current: f64) -> i32 {
104        let v_prev = self.v;
105        let mut s = [self.v, self.h, self.n, self.p];
106        for _ in 0..10 {
107            s = self.rk4_substep(s, current);
108        }
109        self.v = s[0];
110        self.h = s[1];
111        self.n = s[2];
112        self.p = s[3];
113        if self.v >= self.v_threshold && v_prev < self.v_threshold {
114            1
115        } else {
116            0
117        }
118    }
119    pub fn reset(&mut self) {
120        self.v = -65.0;
121        self.h = 0.9;
122        self.n = 0.1;
123        self.p = 0.0;
124    }
125}
126impl Default for GolombFSNeuron {
127    fn default() -> Self {
128        Self::new()
129    }
130}
131
132#[cfg(test)]
133mod tests {
134    use super::*;
135
136    #[test]
137    fn default_matches_constructor_state() {
138        let default = GolombFSNeuron::default();
139        let constructed = GolombFSNeuron::new();
140        assert_eq!(default.v, constructed.v);
141    }
142
143    #[test]
144    fn golomb_fires() {
145        let mut n = GolombFSNeuron::new();
146        let t: i32 = (0..2000).map(|_| n.step(200.0)).sum();
147        assert!(t > 0);
148    }
149
150    // -- GolombFS --
151    #[test]
152    fn golomb_silent_without_input() {
153        let mut n = GolombFSNeuron::new();
154        let t: i32 = (0..200).map(|_| n.step(0.0)).sum();
155        assert_eq!(t, 0);
156    }
157    #[test]
158    fn golomb_reset_clears_state() {
159        let mut n = GolombFSNeuron::new();
160        for _ in 0..100 {
161            n.step(200.0);
162        }
163        n.reset();
164        assert!((n.v - (-65.0)).abs() < 1e-10);
165    }
166    #[test]
167    fn golomb_extreme_bounded() {
168        // Golomb et al. 2007: n^4 kinetics diverge at extreme I; test at high but realistic drive
169        let mut n = GolombFSNeuron::new();
170        for _ in 0..200 {
171            n.step(200.0);
172        }
173        assert!(n.v.is_finite());
174    }
175    #[test]
176    fn golomb_kv3_enables_fast_spiking() {
177        // Kv3 current enables high-frequency firing
178        let mut n = GolombFSNeuron::new();
179        let t: i32 = (0..5000).map(|_| n.step(300.0)).sum();
180        assert!(t > 0, "Golomb FS should fire with high input, got {}", t);
181    }
182    #[test]
183    fn golomb_negative_no_crash() {
184        let mut n = GolombFSNeuron::new();
185        for _ in 0..200 {
186            n.step(-100.0);
187        }
188        assert!(n.v.is_finite());
189    }
190    #[test]
191    fn golomb_nan_no_panic() {
192        let mut n = GolombFSNeuron::new();
193        n.step(f64::NAN);
194    }
195}
sc_neurocore_engine/neurons/biophysical/golomb_fs.rs

sc_neurocore_engine/neurons/biophysical/
golomb_fs.rs
