sc_neurocore_engine/neurons/channels/
persistent_na.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 — Persistent sodium channel neuron
8
9use crate::neurons::biophysical::safe_rate;
10
11/// Persistent Na+ (INaP) neuron — WB base + non-inactivating Na+ current.
12///
13/// INaP activates at subthreshold voltages (-60 to -40 mV) and does not
14/// inactivate, providing a sustained depolarising drive. Key mechanism for:
15/// - Subthreshold membrane oscillations (entorhinal cortex, layer II stellate)
16/// - Plateau potentials and bistability (spinal motoneurons)
17/// - Amplification of synaptic inputs near threshold
18/// - Burst generation in respiratory neurons (pre-Bötzinger complex)
19///
20/// Biological context: Crill, Annu Rev Physiol 58:349, 1996; French et al.,
21/// J Gen Physiol 95:1139, 1990. The repository-specific WB+INaP recurrence
22/// is an experimental composite, not a publication-exact model.
23#[derive(Clone, Debug)]
24pub struct PersistentNaNeuron {
25    pub v: f64,
26    pub h: f64, // Transient Na+ inactivation
27    pub n: f64, // Kdr activation
28    pub p: f64, // INaP activation (slow)
29    // Conductances (mS/cm²)
30    pub g_na: f64,  // Transient Na+
31    pub g_nap: f64, // Persistent Na+
32    pub g_k: f64,   // Kdr
33    pub g_l: f64,
34    // Reversal potentials (mV)
35    pub e_na: f64,
36    pub e_k: f64,
37    pub e_l: f64,
38    pub c_m: f64,
39    pub phi: f64,
40    pub dt: f64,
41    pub v_threshold: f64,
42    pub gain: f64,
43}
44
45impl Default for PersistentNaNeuron {
46    fn default() -> Self {
47        Self::new()
48    }
49}
50
51impl PersistentNaNeuron {
52    pub fn new() -> Self {
53        Self {
54            v: -65.0,
55            h: 0.6,
56            n: 0.32,
57            p: 0.0,
58            g_na: 35.0,
59            g_nap: 0.15, // Persistent Na+ — small but significant
60            g_k: 9.0,
61            g_l: 0.3, // Higher leak to counteract INaP window current
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            dt: 0.5,
68            v_threshold: -20.0,
69            gain: 1.0,
70        }
71    }
72
73    fn valid(&self) -> bool {
74        let finite = [
75            self.v,
76            self.h,
77            self.n,
78            self.p,
79            self.g_na,
80            self.g_nap,
81            self.g_k,
82            self.g_l,
83            self.e_na,
84            self.e_k,
85            self.e_l,
86            self.c_m,
87            self.phi,
88            self.dt,
89            self.v_threshold,
90            self.gain,
91        ]
92        .into_iter()
93        .all(f64::is_finite);
94        finite
95            && (-100.0..=60.0).contains(&self.v)
96            && [self.h, self.n, self.p]
97                .into_iter()
98                .all(|gate| (0.0..=1.0).contains(&gate))
99            && (0.0..=200.0).contains(&self.g_na)
100            && (0.0..=20.0).contains(&self.g_nap)
101            && (0.0..=100.0).contains(&self.g_k)
102            && (0.0..=5.0).contains(&self.g_l)
103            && (30.0..=70.0).contains(&self.e_na)
104            && (-100.0..=-70.0).contains(&self.e_k)
105            && (-80.0..=-40.0).contains(&self.e_l)
106            && (0.5..=2.0).contains(&self.c_m)
107            && (0.5..=10.0).contains(&self.phi)
108            && self.dt > 0.0
109            && self.dt <= 1.0
110            && (-20.0..=20.0).contains(&self.v_threshold)
111            && (0.0..=10.0).contains(&self.gain)
112    }
113
114    pub fn try_step(&mut self, current: f64) -> Result<i32, &'static str> {
115        if !current.is_finite() {
116            return Err("current must be finite");
117        }
118        if !self.valid() {
119            return Err("PersistentNa state and parameters must satisfy the public bounds");
120        }
121
122        let mut candidate = self.clone();
123        let input = candidate.gain * current;
124        let sub_steps = 50;
125        let sub_dt = candidate.dt / sub_steps as f64;
126        let mut fired = 0i32;
127
128        for _ in 0..sub_steps {
129            let v = candidate.v;
130
131            // WB alpha/beta rates
132            let alpha_m = safe_rate(0.1, 35.0, v, 10.0, 1.0);
133            let beta_m = 4.0 * (-(v + 60.0) / 18.0).exp();
134            let m_inf = alpha_m / (alpha_m + beta_m);
135
136            let alpha_h = 0.07 * (-(v + 58.0) / 20.0).exp();
137            let beta_h = 1.0 / (1.0 + (-(v + 28.0) / 10.0).exp());
138
139            let alpha_n = safe_rate(0.01, 34.0, v, 10.0, 0.1);
140            let beta_n = 0.125 * (-(v + 44.0) / 80.0).exp();
141
142            // Persistent Na+ gating: slow activation, no inactivation
143            // Half-activation at -48 mV (subthreshold), tau ~10-50 ms
144            let p_inf = 1.0 / (1.0 + (-(v + 48.0) / 5.0).exp());
145            let tau_p = 10.0 + 40.0 / (1.0 + ((v + 48.0) / 10.0).powi(2));
146
147            // Gate updates
148            candidate.h +=
149                sub_dt * candidate.phi * (alpha_h * (1.0 - candidate.h) - beta_h * candidate.h);
150            candidate.n +=
151                sub_dt * candidate.phi * (alpha_n * (1.0 - candidate.n) - beta_n * candidate.n);
152            candidate.p += sub_dt * (p_inf - candidate.p) / tau_p;
153
154            // Currents
155            let i_na = candidate.g_na * m_inf.powi(3) * candidate.h * (v - candidate.e_na);
156            let i_nap = candidate.g_nap * candidate.p * (v - candidate.e_na);
157            let i_k = candidate.g_k * candidate.n.powi(4) * (v - candidate.e_k);
158            let i_l = candidate.g_l * (v - candidate.e_l);
159
160            let dv = (-i_na - i_nap - i_k - i_l + input) / candidate.c_m;
161            candidate.v += sub_dt * dv;
162            if ![candidate.v, candidate.h, candidate.n, candidate.p]
163                .into_iter()
164                .all(f64::is_finite)
165            {
166                return Err("PersistentNa candidate state became non-finite");
167            }
168
169            if candidate.v >= candidate.v_threshold {
170                fired = 1;
171                candidate.v = -65.0;
172            }
173        }
174
175        candidate.v = candidate.v.clamp(-100.0, 60.0);
176        candidate.h = candidate.h.clamp(0.0, 1.0);
177        candidate.n = candidate.n.clamp(0.0, 1.0);
178        candidate.p = candidate.p.clamp(0.0, 1.0);
179        *self = candidate;
180
181        Ok(fired)
182    }
183
184    pub fn step(&mut self, current: f64) -> i32 {
185        self.try_step(current).unwrap_or(0)
186    }
187
188    pub fn reset(&mut self) {
189        self.v = -65.0;
190        self.h = 0.6;
191        self.n = 0.32;
192        self.p = 0.0;
193    }
194}
195
196#[cfg(test)]
197mod tests {
198    use super::*;
199
200    // -- Persistent Na+ Neuron tests --
201
202    #[test]
203    fn nap_fires_with_input() {
204        let mut n = PersistentNaNeuron::new();
205        let mut spikes = 0;
206        for _ in 0..2_000 {
207            spikes += n.step(2.0);
208        }
209        assert!(spikes > 5, "NaP neuron must fire with input, got {spikes}");
210    }
211
212    #[test]
213    fn nap_subthreshold_oscillations() {
214        // INaP neurons often show subthreshold oscillations or low-rate
215        // spontaneous firing — this is a biological feature, not a bug.
216        // With negative input, INaP should be suppressed.
217        let mut n = PersistentNaNeuron::new();
218        let mut spikes_inhibited = 0;
219        for _ in 0..10_000 {
220            spikes_inhibited += n.step(-2.0);
221        }
222        assert_eq!(
223            spikes_inhibited, 0,
224            "INaP neuron must be silent with inhibitory input, got {spikes_inhibited}"
225        );
226    }
227
228    #[test]
229    fn nap_lowers_threshold() {
230        // INaP provides subthreshold depolarisation → lower effective threshold
231        let mut with_nap = PersistentNaNeuron::new();
232        let mut no_nap = PersistentNaNeuron::new();
233        no_nap.g_nap = 0.0;
234
235        // Use near-threshold input
236        let input = 1.0;
237        let mut spikes_nap = 0;
238        let mut spikes_no = 0;
239        for _ in 0..10_000 {
240            spikes_nap += with_nap.step(input);
241            spikes_no += no_nap.step(input);
242        }
243        assert!(
244            spikes_nap >= spikes_no,
245            "INaP must lower effective threshold: NaP={spikes_nap} vs none={spikes_no}"
246        );
247    }
248
249    #[test]
250    fn nap_p_gate_activates_at_subthreshold() {
251        // At -50 mV (subthreshold), p_inf should be significant
252        let mut n = PersistentNaNeuron::new();
253        n.v = -50.0;
254        // Step a few times for p to converge
255        for _ in 0..1000 {
256            // Hold at -50 mV artificially by resetting v each step
257            let _ = n.step(0.0);
258        }
259        // p_inf at -50 mV = 1/(1+exp(2/5)) = 1/(1+1.49) = 0.40
260        // After many steps p should approach p_inf
261        assert!(
262            n.p > 0.01,
263            "p gate must activate at subthreshold voltages, p={}",
264            n.p
265        );
266    }
267
268    #[test]
269    fn nap_increases_firing_rate() {
270        // Same input, higher g_nap → more spikes
271        let mut low = PersistentNaNeuron::new();
272        low.g_nap = 0.2;
273        let mut high = PersistentNaNeuron::new();
274        high.g_nap = 1.5;
275
276        let input = 1.5;
277        let mut spikes_low = 0;
278        let mut spikes_high = 0;
279        for _ in 0..10_000 {
280            spikes_low += low.step(input);
281            spikes_high += high.step(input);
282        }
283        assert!(
284            spikes_high >= spikes_low,
285            "Higher g_nap must increase firing: high={spikes_high} vs low={spikes_low}"
286        );
287    }
288
289    #[test]
290    fn nap_negative_input_no_crash() {
291        let mut n = PersistentNaNeuron::new();
292        for _ in 0..10_000 {
293            n.step(-100.0);
294        }
295        assert!(n.v.is_finite());
296        assert!(n.v >= -100.0);
297    }
298
299    #[test]
300    fn nap_nan_input_is_rejected_atomically() {
301        let mut n = PersistentNaNeuron::new();
302        let before = n.clone();
303        assert!(n.try_step(f64::NAN).is_err());
304        assert_eq!(n.v, before.v);
305        assert_eq!(n.h, before.h);
306        assert_eq!(n.n, before.n);
307        assert_eq!(n.p, before.p);
308    }
309
310    #[test]
311    fn nap_invalid_configuration_is_rejected_atomically() {
312        let mut n = PersistentNaNeuron::new();
313        n.c_m = 0.0;
314        let before = n.clone();
315        assert!(n.try_step(1.0).is_err());
316        assert_eq!(n.v, before.v);
317        assert_eq!(n.c_m, before.c_m);
318    }
319
320    #[test]
321    fn nap_extreme_input_bounded() {
322        let mut n = PersistentNaNeuron::new();
323        for _ in 0..1000 {
324            n.step(1e6);
325        }
326        assert!(n.v.is_finite() && n.v <= 60.0);
327    }
328
329    #[test]
330    fn nap_reset_clears_state() {
331        let mut n = PersistentNaNeuron::new();
332        n.g_nap = 0.3;
333        for _ in 0..1000 {
334            n.step(10.0);
335        }
336        n.reset();
337        assert_eq!(n.v, -65.0);
338        assert_eq!(n.p, 0.0);
339        assert_eq!(n.h, 0.6);
340        assert_eq!(n.g_nap, 0.3);
341    }
342
343    #[test]
344    fn nap_gates_bounded() {
345        let mut n = PersistentNaNeuron::new();
346        for _ in 0..10_000 {
347            n.step(10.0);
348        }
349        assert!(n.h >= 0.0 && n.h <= 1.0);
350        assert!(n.n >= 0.0 && n.n <= 1.0);
351        assert!(n.p >= 0.0 && n.p <= 1.0);
352    }
353
354    #[test]
355    fn nap_performance_1k_steps() {
356        let start = std::time::Instant::now();
357        let mut n = PersistentNaNeuron::new();
358        for _ in 0..1_000 {
359            std::hint::black_box(n.step(5.0));
360        }
361        let elapsed = start.elapsed();
362        assert!(
363            elapsed.as_millis() < 200,
364            "1k steps must complete in <200ms"
365        );
366    }
367}
sc_neurocore_engine/neurons/channels/persistent_na.rs

sc_neurocore_engine/neurons/channels/
persistent_na.rs
