sc_neurocore_engine/neurons/channels/
sk.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 — SK calcium-activated potassium channel neuron
8
9use crate::neurons::biophysical::safe_rate;
10
11/// SK channel neuron — WB base + Ca2+-only-dependent K+ current.
12///
13/// SK (KCa2.x) channels are activated solely by intracellular Ca2+
14/// (no voltage dependence). They have slower kinetics than BK and produce
15/// the medium afterhyperpolarisation (mAHP) lasting 50-200 ms after spikes.
16///
17/// Key mechanism for:
18/// - Medium AHP (mAHP): limits sustained firing rate
19/// - Spike frequency adaptation: Ca2+ builds → SK activates → firing slows
20/// - Rhythmic firing patterns: SK-mediated pauses create regular ISIs
21/// - Synaptic plasticity gating: SK in dendritic spines regulates NMDA currents
22///
23/// Stocker, Nat Rev Neurosci 5:758, 2004; Wang & Buzsáki, J Neurosci
24/// 16:6402, 1996. The threshold-reset event, the spike-triggered Ca2+
25/// increment, and the Hill constants are repository-specific
26/// specialisations of that review material, not a publication-exact
27/// recurrence.
28#[derive(Clone, Debug)]
29pub struct SKNeuron {
30    pub v: f64,
31    pub h: f64,
32    pub n: f64,
33    pub ca: f64,
34    pub g_na: f64,
35    pub g_k: f64,
36    pub g_sk: f64,
37    pub g_l: f64,
38    pub e_na: f64,
39    pub e_k: f64,
40    pub e_l: f64,
41    pub c_m: f64,
42    pub phi: f64,
43    pub tau_ca: f64,
44    pub dt: f64,
45    pub v_threshold: f64,
46    pub gain: f64,
47}
48
49impl Default for SKNeuron {
50    fn default() -> Self {
51        Self::new()
52    }
53}
54
55impl SKNeuron {
56    pub fn new() -> Self {
57        Self {
58            v: -65.0,
59            h: 0.6,
60            n: 0.32,
61            ca: 0.0,
62            g_na: 35.0,
63            g_k: 9.0,
64            g_sk: 2.0,
65            g_l: 0.1,
66            e_na: 55.0,
67            e_k: -90.0,
68            e_l: -65.0,
69            c_m: 1.0,
70            phi: 5.0,
71            tau_ca: 150.0, // Slower Ca2+ decay than BK → longer mAHP
72            dt: 0.5,
73            v_threshold: -20.0,
74            gain: 1.0,
75        }
76    }
77
78    fn valid(&self) -> bool {
79        let finite = [
80            self.v,
81            self.h,
82            self.n,
83            self.ca,
84            self.g_na,
85            self.g_k,
86            self.g_sk,
87            self.g_l,
88            self.e_na,
89            self.e_k,
90            self.e_l,
91            self.c_m,
92            self.phi,
93            self.tau_ca,
94            self.dt,
95            self.v_threshold,
96            self.gain,
97        ]
98        .into_iter()
99        .all(f64::is_finite);
100        finite
101            && (-100.0..=60.0).contains(&self.v)
102            && [self.h, self.n]
103                .into_iter()
104                .all(|gate| (0.0..=1.0).contains(&gate))
105            && self.ca >= 0.0
106            && (0.0..=200.0).contains(&self.g_na)
107            && (0.0..=100.0).contains(&self.g_k)
108            && (0.0..=50.0).contains(&self.g_sk)
109            && (0.0..=5.0).contains(&self.g_l)
110            && (30.0..=70.0).contains(&self.e_na)
111            && (-100.0..=-70.0).contains(&self.e_k)
112            && (-80.0..=-40.0).contains(&self.e_l)
113            && (0.5..=2.0).contains(&self.c_m)
114            && (0.5..=10.0).contains(&self.phi)
115            && (10.0..=2000.0).contains(&self.tau_ca)
116            && self.dt > 0.0
117            && self.dt <= 1.0
118            && (-20.0..=20.0).contains(&self.v_threshold)
119            && (0.0..=10.0).contains(&self.gain)
120    }
121
122    /// Advance one step after validating the drive and configuration.
123    ///
124    /// Computes the whole update on a candidate clone and commits only on
125    /// success: a non-finite `current`, a configuration outside the public
126    /// bounds, or a non-finite candidate returns `Err` with the pre-step
127    /// state preserved exactly.
128    pub fn try_step(&mut self, current: f64) -> Result<i32, &'static str> {
129        if !current.is_finite() {
130            return Err("current must be finite");
131        }
132        if !self.valid() {
133            return Err("SK state and parameters must satisfy the public bounds");
134        }
135
136        let mut candidate = self.clone();
137        let input = candidate.gain * current;
138        let sub_steps = 50;
139        let sub_dt = candidate.dt / sub_steps as f64;
140        let mut fired = 0i32;
141
142        for _ in 0..sub_steps {
143            let v = candidate.v;
144
145            let alpha_m = safe_rate(0.1, 35.0, v, 10.0, 1.0);
146            let beta_m = 4.0 * (-(v + 60.0) / 18.0).exp();
147            let m_inf = alpha_m / (alpha_m + beta_m);
148
149            let alpha_h = 0.07 * (-(v + 58.0) / 20.0).exp();
150            let beta_h = 1.0 / (1.0 + (-(v + 28.0) / 10.0).exp());
151
152            let alpha_n = safe_rate(0.01, 34.0, v, 10.0, 0.1);
153            let beta_n = 0.125 * (-(v + 44.0) / 80.0).exp();
154
155            // SK activation: purely Ca2+-dependent (Hill function, n=2)
156            let ca2 = candidate.ca * candidate.ca;
157            let sk_inf = ca2 / (ca2 + 0.25); // Half-activation at [Ca2+]=0.5
158
159            // Ca2+ decay
160            candidate.ca += sub_dt * (-candidate.ca / candidate.tau_ca);
161
162            candidate.h +=
163                sub_dt * candidate.phi * (alpha_h * (1.0 - candidate.h) - beta_h * candidate.h);
164            candidate.n +=
165                sub_dt * candidate.phi * (alpha_n * (1.0 - candidate.n) - beta_n * candidate.n);
166
167            let i_na = candidate.g_na * m_inf.powi(3) * candidate.h * (v - candidate.e_na);
168            let i_k = candidate.g_k * candidate.n.powi(4) * (v - candidate.e_k);
169            let i_sk = candidate.g_sk * sk_inf * (v - candidate.e_k);
170            let i_l = candidate.g_l * (v - candidate.e_l);
171
172            let dv = (-i_na - i_k - i_sk - i_l + input) / candidate.c_m;
173            candidate.v += sub_dt * dv;
174            if ![candidate.v, candidate.h, candidate.n, candidate.ca]
175                .into_iter()
176                .all(f64::is_finite)
177            {
178                return Err("SK candidate state became non-finite");
179            }
180
181            if candidate.v >= candidate.v_threshold {
182                fired = 1;
183                candidate.v = -65.0;
184                candidate.ca += 0.2;
185            }
186        }
187
188        candidate.v = candidate.v.clamp(-100.0, 60.0);
189        candidate.h = candidate.h.clamp(0.0, 1.0);
190        candidate.n = candidate.n.clamp(0.0, 1.0);
191        candidate.ca = candidate.ca.max(0.0);
192        *self = candidate;
193
194        Ok(fired)
195    }
196
197    /// Fail-closed wrapper for legacy callers: returns 0 on any rejected
198    /// input without mutating state.
199    pub fn step(&mut self, current: f64) -> i32 {
200        self.try_step(current).unwrap_or(0)
201    }
202
203    /// Restore the dynamic state to its initial values, preserving every
204    /// configuration parameter.
205    pub fn reset(&mut self) {
206        self.v = -65.0;
207        self.h = 0.6;
208        self.n = 0.32;
209        self.ca = 0.0;
210    }
211}
212
213#[cfg(test)]
214mod tests {
215    use super::*;
216
217    // -- SK Neuron tests --
218
219    #[test]
220    fn sk_fires_with_input() {
221        let mut n = SKNeuron::new();
222        let mut spikes = 0;
223        for _ in 0..2_000 {
224            spikes += n.step(2.0);
225        }
226        assert!(spikes > 5, "SK neuron must fire with input, got {spikes}");
227    }
228
229    #[test]
230    fn sk_silent_without_input() {
231        let mut n = SKNeuron::new();
232        let mut spikes = 0;
233        for _ in 0..10_000 {
234            spikes += n.step(0.0);
235        }
236        assert_eq!(
237            spikes, 0,
238            "SK neuron must be silent without input, got {spikes}"
239        );
240    }
241
242    #[test]
243    fn sk_nominal_step_matches_reference_anchor() {
244        let mut n = SKNeuron::new();
245        assert_eq!(n.try_step(5.0), Ok(0));
246        assert!((n.v - -63.180_064_213_072_19).abs() < 1.0e-12);
247        assert!((n.h - 0.648_122_835_749_998_1).abs() < 1.0e-12);
248        assert!((n.n - 0.237_186_365_946_861_5).abs() < 1.0e-12);
249        assert_eq!(n.ca, 0.0);
250    }
251
252    #[test]
253    fn sk_adaptation() {
254        // SK causes spike frequency adaptation
255        let mut n = SKNeuron::new();
256        let input = 5.0;
257        let mut early = 0;
258        for _ in 0..2000 {
259            early += n.step(input);
260        }
261        let mut late = 0;
262        for _ in 0..2000 {
263            late += n.step(input);
264        }
265        assert!(
266            early >= late,
267            "SK should cause adaptation: early={early}, late={late}"
268        );
269    }
270
271    #[test]
272    fn sk_ca_dependent_only() {
273        // SK at rest (ca=0) should contribute zero current
274        let n = SKNeuron::new();
275        let ca2 = n.ca * n.ca;
276        let sk_inf = ca2 / (ca2 + 0.25);
277        assert!(
278            sk_inf < 0.001,
279            "SK must be inactive at ca=0, sk_inf={sk_inf}"
280        );
281    }
282
283    #[test]
284    fn sk_reduces_firing_rate() {
285        let mut with_sk = SKNeuron::new();
286        let mut no_sk = SKNeuron::new();
287        no_sk.g_sk = 0.0;
288
289        let input = 3.0;
290        let mut spikes_sk = 0;
291        let mut spikes_no = 0;
292        for _ in 0..10_000 {
293            spikes_sk += with_sk.step(input);
294            spikes_no += no_sk.step(input);
295        }
296        assert!(
297            spikes_no >= spikes_sk,
298            "SK should reduce firing: SK={spikes_sk} vs none={spikes_no}"
299        );
300    }
301
302    #[test]
303    fn sk_negative_input_no_crash() {
304        let mut n = SKNeuron::new();
305        for _ in 0..10_000 {
306            n.step(-100.0);
307        }
308        assert!(n.v.is_finite());
309    }
310
311    #[test]
312    fn sk_nan_input_is_rejected_atomically() {
313        let mut n = SKNeuron::new();
314        let before = n.clone();
315        assert!(n.try_step(f64::NAN).is_err());
316        assert_eq!(n.v, before.v);
317        assert_eq!(n.h, before.h);
318        assert_eq!(n.n, before.n);
319        assert_eq!(n.ca, before.ca);
320    }
321
322    #[test]
323    fn sk_infinite_input_is_rejected_atomically() {
324        let mut n = SKNeuron::new();
325        let before = n.clone();
326        assert!(n.try_step(f64::INFINITY).is_err());
327        assert!(n.try_step(f64::NEG_INFINITY).is_err());
328        assert_eq!(n.v, before.v);
329        assert_eq!(n.ca, before.ca);
330    }
331
332    #[test]
333    fn sk_invalid_configuration_is_rejected_atomically() {
334        let mut n = SKNeuron::new();
335        n.c_m = 0.0;
336        let before = n.clone();
337        assert!(n.try_step(1.0).is_err());
338        assert_eq!(n.v, before.v);
339        assert_eq!(n.c_m, before.c_m);
340    }
341
342    #[test]
343    fn sk_extreme_input_bounded() {
344        let mut n = SKNeuron::new();
345        for _ in 0..1000 {
346            n.step(1e6);
347        }
348        assert!(n.v.is_finite() && n.v <= 60.0);
349    }
350
351    #[test]
352    fn sk_reset_clears_state() {
353        let mut n = SKNeuron::new();
354        for _ in 0..1000 {
355            n.step(10.0);
356        }
357        n.reset();
358        assert_eq!(n.v, -65.0);
359        assert_eq!(n.ca, 0.0);
360    }
361
362    #[test]
363    fn sk_reset_preserves_parameters() {
364        let mut n = SKNeuron::new();
365        n.g_sk = 4.0;
366        for _ in 0..100 {
367            n.step(5.0);
368        }
369        n.reset();
370        assert_eq!(n.v, -65.0);
371        assert_eq!(n.g_sk, 4.0);
372    }
373
374    #[test]
375    fn sk_performance_1k_steps() {
376        let start = std::time::Instant::now();
377        let mut n = SKNeuron::new();
378        for _ in 0..1_000 {
379            std::hint::black_box(n.step(3.0));
380        }
381        let elapsed = start.elapsed();
382        assert!(
383            elapsed.as_millis() < 200,
384            "1k steps must complete in <200ms"
385        );
386    }
387}
sc_neurocore_engine/neurons/channels/sk.rs

sc_neurocore_engine/neurons/channels/
sk.rs
