sc_neurocore_engine/neurons/interneurons/
sst_neuron.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 — SST interneuron model
8
9/// SST+ (somatostatin) low-threshold spiking interneuron.
10///
11/// Biophysics: Na+, K+, M-current (Kv7, slow K+ for adaptation),
12/// T-type Ca2+ (low-threshold burst), h-current (Ih, sag), leak.
13/// Key properties: spike frequency adaptation, rebound bursting,
14/// facilitating synapses, dendritic targeting.
15///
16/// Based on Pospischil et al. 2008 LTS parameterisation.
17#[derive(Clone, Debug)]
18pub struct SSTNeuron {
19    pub v: f64,
20    pub m: f64,
21    pub h: f64,
22    pub n: f64,
23    pub p: f64, // M-current activation
24    pub s: f64, // T-type Ca2+ inactivation
25    pub r: f64, // h-current activation
26    // Conductances
27    pub g_na: f64,
28    pub g_k: f64,
29    pub g_m: f64,
30    pub g_t: f64,
31    pub g_h: f64,
32    pub g_l: f64,
33    // Reversal potentials
34    pub e_na: f64,
35    pub e_k: f64,
36    pub e_ca: f64,
37    pub e_h: f64,
38    pub e_l: f64,
39    pub c_m: f64,
40    pub dt: f64,
41    pub v_threshold: f64,
42}
43
44impl SSTNeuron {
45    pub fn new() -> Self {
46        Self {
47            v: -65.0,
48            m: 0.02,
49            h: 0.8,
50            n: 0.2,
51            p: 0.0,
52            s: 0.9,
53            r: 0.1,
54            g_na: 50.0,
55            g_k: 5.0,
56            g_m: 0.12, // Strong M-current → adaptation
57            g_t: 0.01, // T-type Ca2+ for rebound (minimal window current)
58            g_h: 0.02, // Ih for sag
59            g_l: 0.05, // Leak for resting stability
60            e_na: 50.0,
61            e_k: -90.0,
62            e_ca: 120.0,
63            e_h: -40.0,
64            e_l: -65.0,
65            c_m: 1.0,
66            dt: 0.025,
67            v_threshold: -20.0,
68        }
69    }
70
71    /// Return `[dV, dm, dh, dn, dp, ds, dr]` of the seven-state SST system at one
72    /// consistent state. The Na/K activation rates use the L'Hôpital limit at the
73    /// removable Traub-Miles singularity; β_m carries the published `V - V_T - 40`
74    /// offset (an earlier `-17` offset drove the cell into depolarisation block).
75    fn derivatives(
76        &self,
77        v: f64,
78        m: f64,
79        h: f64,
80        n: f64,
81        p: f64,
82        s: f64,
83        r: f64,
84        current: f64,
85    ) -> [f64; 7] {
86        let dvt = v - (-56.2);
87        let asing = |num: f64, slope: f64, limit: f64| {
88            if num.abs() < 1e-6 {
89                limit
90            } else {
91                num / ((num / slope).exp() - 1.0)
92            }
93        };
94        let alpha_m = -0.32 * asing(dvt - 13.0, -4.0, -4.0);
95        let beta_m = 0.28 * asing(dvt - 40.0, 5.0, 5.0);
96        let alpha_h = 0.128 * (-(dvt - 17.0) / 18.0).exp();
97        let beta_h = 4.0 / (1.0 + (-(dvt - 40.0) / 5.0).exp());
98        let alpha_n = -0.032 * asing(dvt - 15.0, -5.0, -5.0);
99        let beta_n = 0.5 * (-(dvt - 10.0) / 40.0).exp();
100        let dm = alpha_m * (1.0 - m) - beta_m * m;
101        let dh = alpha_h * (1.0 - h) - beta_h * h;
102        let dn = alpha_n * (1.0 - n) - beta_n * n;
103        let p_inf = 1.0 / (1.0 + (-(v + 35.0) / 10.0).exp());
104        let tau_p = 400.0 / (3.3 * ((v + 35.0) / 20.0).exp() + (-(v + 35.0) / 20.0).exp());
105        let dp = (p_inf - p) / tau_p;
106        let m_t_inf = 1.0 / (1.0 + (-(v + 57.0) / 6.2).exp());
107        let s_inf = 1.0 / (1.0 + ((v + 81.0) / 4.0).exp());
108        let tau_s = 30.0 + 200.0 / (1.0 + ((v + 70.0) / 5.0).exp());
109        let ds = (s_inf - s) / tau_s;
110        let r_inf = 1.0 / (1.0 + ((v + 80.0) / 10.0).exp());
111        let tau_r = 100.0 + 500.0 / ((-(v + 70.0) / 20.0).exp() + ((v + 70.0) / 20.0).exp());
112        let dr = (r_inf - r) / tau_r;
113        let i_na = self.g_na * m * m * m * h * (v - self.e_na);
114        let i_k = self.g_k * n * n * n * n * (v - self.e_k);
115        let i_m = self.g_m * p * (v - self.e_k);
116        let i_t = self.g_t * m_t_inf * m_t_inf * s * (v - self.e_ca);
117        let i_h = self.g_h * r * (v - self.e_h);
118        let i_l = self.g_l * (v - self.e_l);
119        let dvdt = (-i_na - i_k - i_m - i_t - i_h - i_l + current) / self.c_m;
120        [dvdt, dm, dh, dn, dp, ds, dr]
121    }
122
123    /// Return one classical RK4 increment of `[V, m, h, n, p, s, r]`, holding
124    /// `current` constant across the four stages.
125    fn rk4_substep(&self, st: [f64; 7], current: f64) -> [f64; 7] {
126        let dt = self.dt;
127        let k1 = self.derivatives(st[0], st[1], st[2], st[3], st[4], st[5], st[6], current);
128        let mut a = [0.0_f64; 7];
129        for i in 0..7 {
130            a[i] = st[i] + 0.5 * dt * k1[i];
131        }
132        let k2 = self.derivatives(a[0], a[1], a[2], a[3], a[4], a[5], a[6], current);
133        for i in 0..7 {
134            a[i] = st[i] + 0.5 * dt * k2[i];
135        }
136        let k3 = self.derivatives(a[0], a[1], a[2], a[3], a[4], a[5], a[6], current);
137        for i in 0..7 {
138            a[i] = st[i] + dt * k3[i];
139        }
140        let k4 = self.derivatives(a[0], a[1], a[2], a[3], a[4], a[5], a[6], current);
141        let mut out = [0.0_f64; 7];
142        for i in 0..7 {
143            out[i] = st[i] + dt * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]) / 6.0;
144        }
145        out
146    }
147
148    pub fn step(&mut self, current: f64) -> i32 {
149        let v_prev = self.v;
150        let mut st = [self.v, self.m, self.h, self.n, self.p, self.s, self.r];
151        for _ in 0..4 {
152            st = self.rk4_substep(st, current);
153        }
154        self.v = st[0];
155        self.m = st[1];
156        self.h = st[2];
157        self.n = st[3];
158        self.p = st[4];
159        self.s = st[5];
160        self.r = st[6];
161        if self.v >= self.v_threshold && v_prev < self.v_threshold {
162            1
163        } else {
164            0
165        }
166    }
167
168    pub fn reset(&mut self) {
169        self.v = -65.0;
170        self.m = 0.02;
171        self.h = 0.8;
172        self.n = 0.2;
173        self.p = 0.0;
174        self.s = 0.9;
175        self.r = 0.1;
176    }
177}
178
179impl Default for SSTNeuron {
180    fn default() -> Self {
181        Self::new()
182    }
183}
184
185// ═══════════════════════════════════════════════════════════════════
186// VIP Irregular-Spiking Interneuron
187// ═══════════════════════════════════════════════════════════════════
188
189#[cfg(test)]
190mod tests {
191    use super::*;
192
193    #[test]
194    fn sst_fires_with_input() {
195        let mut n = SSTNeuron::new();
196        let spikes: i32 = (0..10000).map(|_| n.step(5.0)).sum();
197        assert!(spikes > 0, "SST+ must fire with sustained input");
198    }
199
200    #[test]
201    fn sst_no_fire_without_input() {
202        let mut n = SSTNeuron::new();
203        let spikes: i32 = (0..5000).map(|_| n.step(0.0)).sum();
204        assert_eq!(spikes, 0);
205    }
206
207    #[test]
208    fn sst_adaptation_reduces_rate() {
209        let mut n = SSTNeuron::new();
210        let first_half: i32 = (0..5000).map(|_| n.step(5.0)).sum();
211        let second_half: i32 = (0..5000).map(|_| n.step(5.0)).sum();
212        // M-current → spike frequency adaptation
213        assert!(
214            second_half <= first_half + 3,
215            "SST+ should adapt: first={first_half}, second={second_half}"
216        );
217    }
218
219    #[test]
220    fn sst_reset_roundtrip() {
221        let mut n = SSTNeuron::new();
222        for _ in 0..5000 {
223            n.step(5.0);
224        }
225        n.reset();
226        let mut fresh = SSTNeuron::new();
227        let r1: i32 = (0..2000).map(|_| n.step(5.0)).sum();
228        let r2: i32 = (0..2000).map(|_| fresh.step(5.0)).sum();
229        assert_eq!(r1, r2);
230    }
231
232    #[test]
233    fn sst_voltage_bounded() {
234        let mut n = SSTNeuron::new();
235        for _ in 0..20000 {
236            n.step(10.0);
237        }
238        assert!(n.v.is_finite());
239        assert!(n.p.is_finite());
240        assert!(n.s.is_finite());
241    }
242
243    #[test]
244    #[ignore = "wall-clock performance smoke; use Criterion benches for timing evidence"]
245    fn sst_performance_10k_steps() {
246        let mut n = SSTNeuron::new();
247        let start = std::time::Instant::now();
248        for _ in 0..10_000 {
249            n.step(5.0);
250        }
251        let elapsed = start.elapsed();
252        assert!(
253            elapsed.as_millis() < 500,
254            "10k SST steps took {:?}",
255            elapsed
256        );
257    }
258}
sc_neurocore_engine/neurons/interneurons/sst_neuron.rs

sc_neurocore_engine/neurons/interneurons/
sst_neuron.rs
