sc_neurocore_engine/neurons/channels/
bk.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 — BK calcium-activated potassium channel neuron
8
9use crate::neurons::biophysical::safe_rate;
10
11/// BK channel neuron — WB base + voltage- and Ca2+-dependent K+ current.
12///
13/// BK (big conductance, MaxiK) channels are activated by both membrane
14/// depolarisation and intracellular Ca2+. They have the largest single-
15/// channel conductance (~250 pS) of any K+ channel. During action
16/// potentials, Ca2+ influx through voltage-gated Ca2+ channels activates
17/// BK, producing fast repolarisation and a prominent fast AHP.
18///
19/// Key mechanism for:
20/// - Fast afterhyperpolarisation (fAHP): rapid spike repolarisation
21/// - Action potential narrowing: BK shortens AP duration
22/// - Burst termination: accumulated Ca2+ activates BK, ending burst
23/// - High-frequency firing: fast repolarisation enables rapid recovery
24///
25/// Bhatt & Storm, J Physiol 557:329, 2003; Faber & Bhatt, PNAS 100:2813, 2003.
26#[derive(Clone, Debug)]
27pub struct BKNeuron {
28    pub v: f64,
29    pub h: f64,  // Na+ inactivation
30    pub n: f64,  // Kdr activation
31    pub ca: f64, // Intracellular Ca2+ concentration
32    // Conductances (mS/cm²)
33    pub g_na: f64,
34    pub g_k: f64,
35    pub g_bk: f64, // BK conductance
36    pub g_l: f64,
37    // Reversal potentials (mV)
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, // Ca2+ decay time constant (ms)
44    pub dt: f64,
45    pub v_threshold: f64,
46    pub gain: f64,
47}
48
49impl Default for BKNeuron {
50    fn default() -> Self {
51        Self::new()
52    }
53}
54
55impl BKNeuron {
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_bk: 3.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: 50.0,
72            dt: 0.5,
73            v_threshold: -20.0,
74            gain: 1.0,
75        }
76    }
77
78    pub fn step(&mut self, current: f64) -> i32 {
79        let input = self.gain * current;
80        let sub_steps = 50;
81        let sub_dt = self.dt / sub_steps as f64;
82        let mut fired = 0i32;
83
84        for _ in 0..sub_steps {
85            let v = self.v;
86
87            let alpha_m = safe_rate(0.1, 35.0, v, 10.0, 1.0);
88            let beta_m = 4.0 * (-(v + 60.0) / 18.0).exp();
89            let m_inf = alpha_m / (alpha_m + beta_m);
90
91            let alpha_h = 0.07 * (-(v + 58.0) / 20.0).exp();
92            let beta_h = 1.0 / (1.0 + (-(v + 28.0) / 10.0).exp());
93
94            let alpha_n = safe_rate(0.01, 34.0, v, 10.0, 0.1);
95            let beta_n = 0.125 * (-(v + 44.0) / 80.0).exp();
96
97            // BK activation: joint voltage and Ca2+ dependence
98            // Half-activation shifts left (easier) with higher Ca2+
99            // BK half-activation: ~+10 mV without Ca2+, shifts to -20 mV with high Ca2+
100            let v_half_bk = 10.0 - 30.0 * (self.ca / (self.ca + 0.5));
101            let bk_inf = 1.0 / (1.0 + (-(v - v_half_bk) / 15.0).exp());
102
103            // Ca2+ dynamics: decay + spike-triggered influx
104            self.ca += sub_dt * (-self.ca / self.tau_ca);
105
106            self.h += sub_dt * self.phi * (alpha_h * (1.0 - self.h) - beta_h * self.h);
107            self.n += sub_dt * self.phi * (alpha_n * (1.0 - self.n) - beta_n * self.n);
108
109            let i_na = self.g_na * m_inf.powi(3) * self.h * (v - self.e_na);
110            let i_k = self.g_k * self.n.powi(4) * (v - self.e_k);
111            let i_bk = self.g_bk * bk_inf * (v - self.e_k);
112            let i_l = self.g_l * (v - self.e_l);
113
114            let dv = (-i_na - i_k - i_bk - i_l + input) / self.c_m;
115            self.v += sub_dt * dv;
116
117            if self.v >= self.v_threshold {
118                fired = 1;
119                self.v = -65.0;
120                self.ca += 0.3; // Ca2+ influx on spike
121            }
122        }
123
124        self.v = self.v.clamp(-100.0, 60.0);
125        if !self.v.is_finite() {
126            self.v = -65.0;
127            self.h = 0.6;
128            self.n = 0.32;
129        }
130        if !self.ca.is_finite() {
131            self.ca = 0.0;
132        }
133        self.h = self.h.clamp(0.0, 1.0);
134        self.n = self.n.clamp(0.0, 1.0);
135        self.ca = self.ca.max(0.0);
136
137        fired
138    }
139
140    pub fn reset(&mut self) {
141        *self = Self::new();
142    }
143}
144
145#[cfg(test)]
146mod tests {
147    use super::*;
148
149    // -- BK Neuron tests --
150
151    #[test]
152    fn bk_fires_with_input() {
153        let mut n = BKNeuron::new();
154        let mut spikes = 0;
155        for _ in 0..2_000 {
156            spikes += n.step(3.0);
157        }
158        assert!(spikes > 5, "BK neuron must fire with input, got {spikes}");
159    }
160
161    #[test]
162    fn bk_silent_without_input() {
163        let mut n = BKNeuron::new();
164        let mut spikes = 0;
165        for _ in 0..10_000 {
166            spikes += n.step(0.0);
167        }
168        assert_eq!(
169            spikes, 0,
170            "BK neuron must be silent without input, got {spikes}"
171        );
172    }
173
174    #[test]
175    fn bk_ca_accumulates_during_spiking() {
176        let mut n = BKNeuron::new();
177        assert_eq!(n.ca, 0.0);
178        for _ in 0..5000 {
179            n.step(5.0);
180        }
181        assert!(
182            n.ca > 0.0,
183            "Ca2+ must accumulate during spiking, ca={}",
184            n.ca
185        );
186    }
187
188    #[test]
189    fn bk_deepens_ahp() {
190        // BK should produce deeper AHP (more negative post-spike voltage)
191        // Compare with and without BK after a burst
192        let mut with_bk = BKNeuron::new();
193        let mut no_bk = BKNeuron::new();
194        no_bk.g_bk = 0.0;
195
196        // Drive both to spike, then check voltage
197        for _ in 0..2000 {
198            with_bk.step(5.0);
199            no_bk.step(5.0);
200        }
201        // After sustained spiking, BK with Ca2+ should keep voltage lower
202        // (stronger K+ current from BK)
203        // Test that BK neuron has non-zero Ca2+ (proves it's active)
204        assert!(with_bk.ca > 0.0, "BK neuron must have Ca2+ after spiking");
205    }
206
207    #[test]
208    fn bk_reduces_firing_rate() {
209        // BK should reduce firing rate via stronger repolarisation
210        let mut with_bk = BKNeuron::new();
211        let mut no_bk = BKNeuron::new();
212        no_bk.g_bk = 0.0;
213
214        let input = 3.0;
215        let mut spikes_bk = 0;
216        let mut spikes_no = 0;
217        for _ in 0..10_000 {
218            spikes_bk += with_bk.step(input);
219            spikes_no += no_bk.step(input);
220        }
221        // BK adds extra K+ → fewer spikes (or equal if Ca2+ builds slowly)
222        assert!(
223            spikes_no >= spikes_bk,
224            "BK should reduce firing: BK={spikes_bk} vs none={spikes_no}"
225        );
226    }
227
228    #[test]
229    fn bk_negative_input_no_crash() {
230        let mut n = BKNeuron::new();
231        for _ in 0..10_000 {
232            n.step(-100.0);
233        }
234        assert!(n.v.is_finite());
235    }
236
237    #[test]
238    fn bk_nan_input_stays_finite() {
239        let mut n = BKNeuron::new();
240        n.step(f64::NAN);
241        assert!(n.v.is_finite());
242    }
243
244    #[test]
245    fn bk_extreme_input_bounded() {
246        let mut n = BKNeuron::new();
247        for _ in 0..1000 {
248            n.step(1e6);
249        }
250        assert!(n.v.is_finite() && n.v <= 60.0);
251    }
252
253    #[test]
254    fn bk_reset_clears_state() {
255        let mut n = BKNeuron::new();
256        for _ in 0..1000 {
257            n.step(10.0);
258        }
259        n.reset();
260        assert_eq!(n.v, -65.0);
261        assert_eq!(n.ca, 0.0);
262    }
263
264    #[test]
265    fn bk_performance_1k_steps() {
266        let start = std::time::Instant::now();
267        let mut n = BKNeuron::new();
268        for _ in 0..1_000 {
269            std::hint::black_box(n.step(3.0));
270        }
271        let elapsed = start.elapsed();
272        assert!(
273            elapsed.as_millis() < 200,
274            "1k steps must complete in <200ms"
275        );
276    }
277}
sc_neurocore_engine/neurons/channels/bk.rs

sc_neurocore_engine/neurons/channels/
bk.rs
