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/// Bhatt & Storm, J Physiol 557:329, 2003; Stocker, Nat Rev Neurosci 5:758, 2004.
24#[derive(Clone, Debug)]
25pub struct SKNeuron {
26    pub v: f64,
27    pub h: f64,
28    pub n: f64,
29    pub ca: f64,
30    pub g_na: f64,
31    pub g_k: f64,
32    pub g_sk: f64,
33    pub g_l: f64,
34    pub e_na: f64,
35    pub e_k: f64,
36    pub e_l: f64,
37    pub c_m: f64,
38    pub phi: f64,
39    pub tau_ca: f64,
40    pub dt: f64,
41    pub v_threshold: f64,
42    pub gain: f64,
43}
44
45impl Default for SKNeuron {
46    fn default() -> Self {
47        Self::new()
48    }
49}
50
51impl SKNeuron {
52    pub fn new() -> Self {
53        Self {
54            v: -65.0,
55            h: 0.6,
56            n: 0.32,
57            ca: 0.0,
58            g_na: 35.0,
59            g_k: 9.0,
60            g_sk: 2.0,
61            g_l: 0.1,
62            e_na: 55.0,
63            e_k: -90.0,
64            e_l: -65.0,
65            c_m: 1.0,
66            phi: 5.0,
67            tau_ca: 150.0, // Slower Ca2+ decay than BK → longer mAHP
68            dt: 0.5,
69            v_threshold: -20.0,
70            gain: 1.0,
71        }
72    }
73
74    pub fn step(&mut self, current: f64) -> i32 {
75        let input = self.gain * current;
76        let sub_steps = 50;
77        let sub_dt = self.dt / sub_steps as f64;
78        let mut fired = 0i32;
79
80        for _ in 0..sub_steps {
81            let v = self.v;
82
83            let alpha_m = safe_rate(0.1, 35.0, v, 10.0, 1.0);
84            let beta_m = 4.0 * (-(v + 60.0) / 18.0).exp();
85            let m_inf = alpha_m / (alpha_m + beta_m);
86
87            let alpha_h = 0.07 * (-(v + 58.0) / 20.0).exp();
88            let beta_h = 1.0 / (1.0 + (-(v + 28.0) / 10.0).exp());
89
90            let alpha_n = safe_rate(0.01, 34.0, v, 10.0, 0.1);
91            let beta_n = 0.125 * (-(v + 44.0) / 80.0).exp();
92
93            // SK activation: purely Ca2+-dependent (Hill function, n=2)
94            let ca2 = self.ca * self.ca;
95            let sk_inf = ca2 / (ca2 + 0.25); // Half-activation at [Ca2+]=0.5
96
97            // Ca2+ decay
98            self.ca += sub_dt * (-self.ca / self.tau_ca);
99
100            self.h += sub_dt * self.phi * (alpha_h * (1.0 - self.h) - beta_h * self.h);
101            self.n += sub_dt * self.phi * (alpha_n * (1.0 - self.n) - beta_n * self.n);
102
103            let i_na = self.g_na * m_inf.powi(3) * self.h * (v - self.e_na);
104            let i_k = self.g_k * self.n.powi(4) * (v - self.e_k);
105            let i_sk = self.g_sk * sk_inf * (v - self.e_k);
106            let i_l = self.g_l * (v - self.e_l);
107
108            let dv = (-i_na - i_k - i_sk - i_l + input) / self.c_m;
109            self.v += sub_dt * dv;
110
111            if self.v >= self.v_threshold {
112                fired = 1;
113                self.v = -65.0;
114                self.ca += 0.2;
115            }
116        }
117
118        self.v = self.v.clamp(-100.0, 60.0);
119        if !self.v.is_finite() {
120            self.v = -65.0;
121            self.h = 0.6;
122            self.n = 0.32;
123        }
124        if !self.ca.is_finite() {
125            self.ca = 0.0;
126        }
127        self.h = self.h.clamp(0.0, 1.0);
128        self.n = self.n.clamp(0.0, 1.0);
129        self.ca = self.ca.max(0.0);
130
131        fired
132    }
133
134    pub fn reset(&mut self) {
135        *self = Self::new();
136    }
137}
138
139#[cfg(test)]
140mod tests {
141    use super::*;
142
143    // -- SK Neuron tests --
144
145    #[test]
146    fn sk_fires_with_input() {
147        let mut n = SKNeuron::new();
148        let mut spikes = 0;
149        for _ in 0..2_000 {
150            spikes += n.step(2.0);
151        }
152        assert!(spikes > 5, "SK neuron must fire with input, got {spikes}");
153    }
154
155    #[test]
156    fn sk_silent_without_input() {
157        let mut n = SKNeuron::new();
158        let mut spikes = 0;
159        for _ in 0..10_000 {
160            spikes += n.step(0.0);
161        }
162        assert_eq!(
163            spikes, 0,
164            "SK neuron must be silent without input, got {spikes}"
165        );
166    }
167
168    #[test]
169    fn sk_adaptation() {
170        // SK causes spike frequency adaptation
171        let mut n = SKNeuron::new();
172        let input = 5.0;
173        let mut early = 0;
174        for _ in 0..2000 {
175            early += n.step(input);
176        }
177        let mut late = 0;
178        for _ in 0..2000 {
179            late += n.step(input);
180        }
181        assert!(
182            early >= late,
183            "SK should cause adaptation: early={early}, late={late}"
184        );
185    }
186
187    #[test]
188    fn sk_ca_dependent_only() {
189        // SK at rest (ca=0) should contribute zero current
190        let n = SKNeuron::new();
191        let ca2 = n.ca * n.ca;
192        let sk_inf = ca2 / (ca2 + 0.25);
193        assert!(
194            sk_inf < 0.001,
195            "SK must be inactive at ca=0, sk_inf={sk_inf}"
196        );
197    }
198
199    #[test]
200    fn sk_reduces_firing_rate() {
201        let mut with_sk = SKNeuron::new();
202        let mut no_sk = SKNeuron::new();
203        no_sk.g_sk = 0.0;
204
205        let input = 3.0;
206        let mut spikes_sk = 0;
207        let mut spikes_no = 0;
208        for _ in 0..10_000 {
209            spikes_sk += with_sk.step(input);
210            spikes_no += no_sk.step(input);
211        }
212        assert!(
213            spikes_no >= spikes_sk,
214            "SK should reduce firing: SK={spikes_sk} vs none={spikes_no}"
215        );
216    }
217
218    #[test]
219    fn sk_negative_input_no_crash() {
220        let mut n = SKNeuron::new();
221        for _ in 0..10_000 {
222            n.step(-100.0);
223        }
224        assert!(n.v.is_finite());
225    }
226
227    #[test]
228    fn sk_nan_input_stays_finite() {
229        let mut n = SKNeuron::new();
230        n.step(f64::NAN);
231        assert!(n.v.is_finite());
232    }
233
234    #[test]
235    fn sk_extreme_input_bounded() {
236        let mut n = SKNeuron::new();
237        for _ in 0..1000 {
238            n.step(1e6);
239        }
240        assert!(n.v.is_finite() && n.v <= 60.0);
241    }
242
243    #[test]
244    fn sk_reset_clears_state() {
245        let mut n = SKNeuron::new();
246        for _ in 0..1000 {
247            n.step(10.0);
248        }
249        n.reset();
250        assert_eq!(n.v, -65.0);
251        assert_eq!(n.ca, 0.0);
252    }
253
254    #[test]
255    fn sk_performance_1k_steps() {
256        let start = std::time::Instant::now();
257        let mut n = SKNeuron::new();
258        for _ in 0..1_000 {
259            std::hint::black_box(n.step(3.0));
260        }
261        let elapsed = start.elapsed();
262        assert!(
263            elapsed.as_millis() < 200,
264            "1k steps must complete in <200ms"
265        );
266    }
267}
sc_neurocore_engine/neurons/channels/sk.rs

sc_neurocore_engine/neurons/channels/
sk.rs
