sc_neurocore_engine/neurons/misc/
myelinated_axon.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 — Myelinated Axon Models
8
9//! Frankenhaeuser-Huxley and MRG nodal models for myelinated axons.
10
11// ═══════════════════════════════════════════════════════════════════
12// Frankenhaeuser-Huxley Axon
13// ═══════════════════════════════════════════════════════════════════
14
15/// Frankenhaeuser-Huxley 1964 — myelinated nerve fibre model.
16///
17/// Extension of HH for myelinated axons (Xenopus node of Ranvier).
18/// Uses Goldman-Hodgkin-Katz (GHK) permeability-based current equations
19/// with 4 gating variables: m (Na activation), h (Na inactivation),
20/// n (delayed rectifier K), p (slow non-specific current).
21///
22/// GHK current for monovalent ion:
23///   I = P * F²V/(RT) * (C_i - C_o * exp(-FV/RT)) / (1 - exp(-FV/RT))
24///
25/// Simplified (FH convention, V relative to rest, temperature factor absorbed):
26///   I_Na = P_Na * m² * h * ghk_drive(V, Na_i/Na_o)
27///   I_K  = P_K  * n² * ghk_drive(V, K_i/K_o)
28///   I_p  = P_p  * p² * ghk_drive(V, Na_i/Na_o)  [non-specific, Na-like]
29///   I_L  = g_L * (V - E_L)
30///
31/// C dV/dt = -(I_Na + I_K + I_p + I_L) + I_ext
32///
33/// The GHK driving force is the key distinction from conductance-based
34/// HH — it is nonlinear in V and depends on concentration ratios.
35///
36/// Frankenhaeuser & Huxley, J Physiol 171:302, 1964.
37/// Frankenhaeuser, J Physiol 160:46, 1962 (rate constants).
38#[derive(Clone, Debug)]
39pub struct FrankenhaeUserHuxleyAxon {
40    pub v: f64,    // Membrane potential (mV, relative to rest)
41    pub m: f64,    // Na activation
42    pub h: f64,    // Na inactivation
43    pub n: f64,    // K delayed rectifier
44    pub p: f64,    // Slow non-specific
45    pub c_m: f64,  // µF/cm²
46    pub p_na: f64, // Na permeability (10⁻³ cm/s, FH units)
47    pub p_k: f64,  // K permeability
48    pub p_p: f64,  // Slow current permeability
49    pub g_l: f64,  // Leak conductance (mS/cm²)
50    pub e_l: f64,  // Leak reversal (mV relative to rest)
51    // Concentration ratios (C_i / C_o) for GHK
52    pub na_ratio: f64, // [Na]_i / [Na]_o (~0.1 for frog)
53    pub k_ratio: f64,  // [K]_i / [K]_o (~30 for frog)
54    pub v_t: f64,      // RT/F thermal voltage (mV), ~25.3 at 20°C
55    pub dt: f64,
56    pub sub_steps: usize,
57    pub gain: f64,
58}
59
60impl Default for FrankenhaeUserHuxleyAxon {
61    fn default() -> Self {
62        Self::new()
63    }
64}
65
66impl FrankenhaeUserHuxleyAxon {
67    pub fn new() -> Self {
68        Self {
69            v: 0.0, // Relative to resting potential
70            m: 0.005,
71            h: 0.8,
72            n: 0.01,
73            p: 0.01,
74            c_m: 2.0, // µF/cm² (myelinated node, FH Table 4)
75            // Effective permeabilities: P_raw * F * [C]_o / 1000
76            // FH Table 4: P_Na=8e-3, P_K=1.2e-3, P_p=0.54e-3 cm/s
77            // [Na]_o=114.5, [K]_o=2.5 mM; F=96.485 C/mmol
78            p_na: 88.4,     // 8e-3 * 96.485 * 114.5 / 1000 (mA/cm² per unit gating)
79            p_k: 0.29,      // 1.2e-3 * 96.485 * 2.5 / 1000
80            p_p: 5.96,      // 0.54e-3 * 96.485 * 114.5 / 1000 (Na-like)
81            g_l: 30.3,      // Leak (FH Table 4: 30.3 mS/cm²)
82            e_l: 0.026,     // Leak reversal (FH Table 4: 0.026 mV)
83            na_ratio: 0.12, // [Na]_i/[Na]_o = 13.74/114.5 (FH 1962)
84            k_ratio: 48.0,  // [K]_i/[K]_o = 120/2.5 (FH 1962)
85            v_t: 25.3,      // RT/F at 20°C (293K)
86            dt: 0.5,        // External step (ms)
87            sub_steps: 50,  // dt_sub = 0.01 ms
88            gain: 1.0,
89        }
90    }
91
92    /// GHK current for monovalent ion (FH convention).
93    ///
94    /// Returns current density contribution in mA/cm² when P is in
95    /// FH-scaled units (absorbs F*[C_o]*1e-3 into P).
96    ///
97    /// I = P_eff * gates * V/V_T * (r - exp(-V/V_T)) / (1 - exp(-V/V_T))
98    ///
99    /// where r = [ion]_i / [ion]_o, V_T = RT/F.
100    /// At V→0: uses L'Hôpital limit = P_eff * (r - 1).
101    ///
102    /// The Faraday scaling: P_eff = P_raw * F * [C]_o / 1000
103    /// is absorbed into the permeability constant (FH Table 4 values
104    /// are already in these effective units: mA/cm² at unit gating).
105    #[inline]
106    fn ghk_current(v: f64, c_ratio: f64, v_t: f64) -> f64 {
107        if v.abs() < 0.01 {
108            // L'Hôpital limit
109            c_ratio - 1.0
110        } else {
111            let u = v / v_t;
112            let exp_neg_u = (-u).exp();
113            u * (c_ratio - exp_neg_u) / (1.0 - exp_neg_u)
114        }
115    }
116
117    pub fn step(&mut self, current: f64) -> i32 {
118        let input = self.gain * current;
119        let dt_sub = self.dt / self.sub_steps as f64;
120        let v_prev = self.v;
121
122        for _ in 0..self.sub_steps {
123            let v = self.v;
124
125            // FH alpha/beta rate functions (Frankenhaeuser 1962, Table 1)
126            // All rates in ms⁻¹, V in mV relative to rest.
127
128            // m gate (Na activation)
129            let am = if (v - 22.0).abs() < 0.1 {
130                1.87
131            } else {
132                0.36 * (v - 22.0) / (1.0 - (-(v - 22.0) / 3.0).exp())
133            };
134            let bm = if (v - 13.0).abs() < 0.1 {
135                1.87
136            } else {
137                0.4 * (13.0 - v) / (1.0 - ((v - 13.0) / 20.0).exp())
138            };
139
140            // h gate (Na inactivation)
141            let ah = if (v + 10.0).abs() < 0.1 {
142                0.08
143            } else {
144                0.1 * (-10.0 - v) / (1.0 - ((v + 10.0) / 6.0).exp())
145            };
146            let bh = 4.5 / (1.0 + ((45.0 - v) / 10.0).exp());
147
148            // n gate (K delayed rectifier)
149            let an = if (v - 13.0).abs() < 0.1 {
150                0.1
151            } else {
152                0.02 * (v - 13.0) / (1.0 - (-(v - 13.0) / 10.0).exp())
153            };
154            let bn = if (v - 23.0).abs() < 0.1 {
155                0.05
156            } else {
157                0.05 * (23.0 - v) / (1.0 - ((v - 23.0) / 10.0).exp())
158            };
159
160            // p gate (slow non-specific)
161            let ap = if (v - 21.0).abs() < 0.1 {
162                0.04
163            } else {
164                0.006 * (v - 21.0) / (1.0 - (-(v - 21.0) / 2.0).exp())
165            };
166            let bp = if (v + 4.0).abs() < 0.1 {
167                0.04
168            } else {
169                0.09 * (-4.0 - v) / (1.0 - ((v + 4.0) / 2.0).exp())
170            };
171
172            // Ensure rates are non-negative
173            let am = am.max(0.0);
174            let bm = bm.max(0.0);
175            let ah = ah.max(0.0);
176            let bh = bh.max(0.0);
177            let an = an.max(0.0);
178            let bn = bn.max(0.0);
179            let ap = ap.max(0.0);
180            let bp = bp.max(0.0);
181
182            // Gate updates
183            self.m += dt_sub * (am * (1.0 - self.m) - bm * self.m);
184            self.h += dt_sub * (ah * (1.0 - self.h) - bh * self.h);
185            self.n += dt_sub * (an * (1.0 - self.n) - bn * self.n);
186            self.p += dt_sub * (ap * (1.0 - self.p) - bp * self.p);
187
188            // Clamp gates
189            self.m = self.m.clamp(0.0, 1.0);
190            self.h = self.h.clamp(0.0, 1.0);
191            self.n = self.n.clamp(0.0, 1.0);
192            self.p = self.p.clamp(0.0, 1.0);
193
194            // GHK permeability-based currents (FH Table 4)
195            let i_na = self.p_na
196                * self.m
197                * self.m
198                * self.h
199                * Self::ghk_current(v, self.na_ratio, self.v_t);
200            let i_k = self.p_k * self.n * self.n * Self::ghk_current(v, self.k_ratio, self.v_t);
201            // p-current uses Na-like concentration ratio (non-specific cation)
202            let i_p = self.p_p * self.p * self.p * Self::ghk_current(v, self.na_ratio, self.v_t);
203            let i_l = self.g_l * (self.v - self.e_l);
204
205            let dv = (-(i_na + i_k + i_p + i_l) + input) / self.c_m;
206            self.v += dt_sub * dv;
207        }
208
209        // Safety
210        self.v = self.v.clamp(-50.0, 150.0);
211        if !self.v.is_finite() {
212            self.v = 0.0;
213        }
214        if !self.m.is_finite() {
215            self.m = 0.005;
216        }
217        if !self.h.is_finite() {
218            self.h = 0.8;
219        }
220        if !self.n.is_finite() {
221            self.n = 0.01;
222        }
223        if !self.p.is_finite() {
224            self.p = 0.01;
225        }
226
227        // Spike detection: V crosses 40 mV upward
228        if self.v >= 40.0 && v_prev < 40.0 {
229            1
230        } else {
231            0
232        }
233    }
234
235    pub fn reset(&mut self) {
236        *self = Self::new();
237    }
238}
239
240// ═══════════════════════════════════════════════════════════════════
241// Node of Ranvier (McIntyre-Richardson-Grill 2002)
242// ═══════════════════════════════════════════════════════════════════
243
244/// Node of Ranvier — McIntyre-Richardson-Grill 2002 model.
245///
246/// Gold-standard nodal model for mammalian myelinated axons. Includes
247/// the specific channel complement of nodes of Ranvier:
248///
249/// - **INaT** (transient Na, Nav1.6): m³h gating, fast activation
250/// - **INaP** (persistent Na, Nav1.6): p³ gating, subthreshold amplification
251/// - **IKs** (slow K, Kv7/KCNQ): s gating, membrane stabilisation
252/// - **IKf** (fast K, Kv3.1-like): fast repolarisation (n⁴ HH-style, optional)
253/// - **IL** (leak)
254///
255/// The persistent Na current (INaP) is critical — it provides subthreshold
256/// amplification and lowers the effective firing threshold, a key feature
257/// of Nav1.6-rich nodes that distinguishes them from generic HH models.
258///
259/// C dV/dt = -(INaT + INaP + IKs + IL) + I_ext
260///
261/// Gating uses Boltzmann steady-state + time-constant formulation
262/// (not alpha/beta) following MRG convention.
263///
264/// McIntyre, Richardson & Grill, J Neurophysiol 87:995, 2002.
265#[derive(Clone, Debug)]
266pub struct NodeOfRanvier {
267    pub v: f64,     // Membrane potential (mV)
268    pub m: f64,     // Nav1.6 transient activation
269    pub h: f64,     // Nav1.6 transient inactivation
270    pub p: f64,     // Nav1.6 persistent activation
271    pub s: f64,     // Kv7 slow K activation
272    pub c_m: f64,   // Nodal capacitance (µF/cm²)
273    pub g_nat: f64, // Transient Na conductance (mS/cm²)
274    pub g_nap: f64, // Persistent Na conductance
275    pub g_ks: f64,  // Slow K (Kv7) conductance
276    pub g_l: f64,   // Leak conductance
277    pub e_na: f64,  // Na reversal (mV)
278    pub e_k: f64,   // K reversal (mV)
279    pub e_l: f64,   // Leak reversal (mV)
280    pub dt: f64,    // External time step (ms)
281    pub sub_steps: usize,
282    pub gain: f64,
283}
284
285impl Default for NodeOfRanvier {
286    fn default() -> Self {
287        Self::new()
288    }
289}
290
291impl NodeOfRanvier {
292    pub fn new() -> Self {
293        Self {
294            v: -80.0,
295            m: 0.01,
296            h: 0.75,
297            p: 0.01,
298            s: 0.05,
299            c_m: 2.0,      // µF/cm² (MRG nodal value)
300            g_nat: 3000.0, // mS/cm² (high Nav1.6 density)
301            g_nap: 5.0,    // Persistent Na (small but critical)
302            g_ks: 80.0,    // Kv7/KCNQ slow K
303            g_l: 7.0,      // Nodal leak (higher than soma)
304            e_na: 50.0,
305            e_k: -90.0,
306            e_l: -90.0,    // MRG nodal resting ~-80 mV
307            dt: 0.5,       // External step (ms)
308            sub_steps: 20, // dt_sub = 0.025 ms
309            gain: 1.0,
310        }
311    }
312
313    /// Boltzmann steady-state: 1 / (1 + exp(-(V - V_half) / k))
314    #[inline]
315    fn boltz(v: f64, v_half: f64, k: f64) -> f64 {
316        1.0 / (1.0 + (-(v - v_half) / k).exp())
317    }
318
319    pub fn step(&mut self, current: f64) -> i32 {
320        let input = self.gain * current;
321        let dt_sub = self.dt / self.sub_steps as f64;
322        let v_prev_ext = self.v;
323
324        for _ in 0..self.sub_steps {
325            let v = self.v;
326
327            // Nav1.6 transient: m gate (fast activation)
328            // MRG: V_half = -26.8 mV, k = 9.2 mV
329            let m_inf = Self::boltz(v, -26.8, 9.2);
330            let tau_m = 0.025 + 0.14 / (1.0 + ((v + 25.0) / 10.0).powi(2)).max(0.01);
331            self.m += dt_sub * (m_inf - self.m) / tau_m;
332
333            // Nav1.6 transient: h gate (inactivation)
334            // MRG: V_half = -55.2 mV, k = -7.4 mV (negative slope)
335            let h_inf = Self::boltz(v, -55.2, -7.4);
336            let tau_h = 0.6 + 4.0 / (1.0 + ((v + 45.0) / 10.0).powi(2)).max(0.01);
337            self.h += dt_sub * (h_inf - self.h) / tau_h;
338
339            // Nav1.6 persistent: p gate (slow activation)
340            // MRG: V_half = -44.0 mV, k = 5.0 mV
341            let p_inf = Self::boltz(v, -44.0, 5.0);
342            let tau_p = 1.0 + 6.0 / (1.0 + ((v + 40.0) / 10.0).powi(2)).max(0.01);
343            self.p += dt_sub * (p_inf - self.p) / tau_p;
344
345            // Kv7 slow K: s gate
346            // MRG: V_half = -30.0 mV, k = 10.0 mV, slow
347            let s_inf = Self::boltz(v, -30.0, 10.0);
348            let tau_s = 20.0 + 60.0 / (1.0 + ((v + 30.0) / 15.0).powi(2)).max(0.01);
349            self.s += dt_sub * (s_inf - self.s) / tau_s;
350
351            // Clamp gates
352            self.m = self.m.clamp(0.0, 1.0);
353            self.h = self.h.clamp(0.0, 1.0);
354            self.p = self.p.clamp(0.0, 1.0);
355            self.s = self.s.clamp(0.0, 1.0);
356
357            // Currents
358            let i_nat = self.g_nat * self.m.powi(3) * self.h * (v - self.e_na);
359            let i_nap = self.g_nap * self.p.powi(3) * (v - self.e_na);
360            let i_ks = self.g_ks * self.s * (v - self.e_k);
361            let i_l = self.g_l * (v - self.e_l);
362
363            let dv = (-(i_nat + i_nap + i_ks + i_l) + input) / self.c_m;
364            self.v += dt_sub * dv;
365        }
366
367        // Safety bounds
368        self.v = self.v.clamp(-120.0, 60.0);
369        if !self.v.is_finite() {
370            self.v = -80.0;
371        }
372        if !self.m.is_finite() {
373            self.m = 0.01;
374        }
375        if !self.h.is_finite() {
376            self.h = 0.75;
377        }
378        if !self.p.is_finite() {
379            self.p = 0.01;
380        }
381        if !self.s.is_finite() {
382            self.s = 0.05;
383        }
384
385        // Spike: V crosses -10 mV upward
386        if self.v >= -10.0 && v_prev_ext < -10.0 {
387            1
388        } else {
389            0
390        }
391    }
392
393    pub fn reset(&mut self) {
394        *self = Self::new();
395    }
396}
397
398// ═══════════════════════════════════════════════════════════════════
399// Myelinated Axon (Saltatory Conduction Segment)
400// ═══════════════════════════════════════════════════════════════════
401
402/// Myelinated axon segment — node of Ranvier + internode cable.
403///
404/// Models a single saltatory conduction unit per the MRG 2002
405/// double-cable architecture: an active node (using the NodeOfRanvier
406/// model) coupled to a passive internode represented as a lumped
407/// RC cable.
408///
409/// The internode has:
410/// - Very low capacitance (~0.001 µF/cm², myelin layers)
411/// - Very high resistance (myelin sheath insulation)
412/// - Paranodal seal resistance (leakage at node-internode junction)
413///
414/// The node voltage drives current through the paranodal seal into
415/// the internode, and the internode voltage feeds back into the node
416/// via the return path. This bidirectional coupling determines the
417/// conduction velocity and safety factor.
418///
419/// The external input represents current arriving from the upstream
420/// node (saltatory propagation).
421///
422/// V_node equation: C_n dV_n/dt = I_ionic(node) + g_para*(V_i - V_n) + I_ext
423/// V_internode equation: C_i dV_i/dt = -g_l_myelin*(V_i - E_l) + g_para*(V_n - V_i)
424///
425/// McIntyre, Richardson & Grill, J Neurophysiol 87:995, 2002.
426/// Richardson et al., Clin Neurophysiol 111:2175, 2000.
427#[derive(Clone, Debug)]
428pub struct MyelinatedAxon {
429    // Node (active, MRG)
430    pub node: NodeOfRanvier,
431    // Internode (passive cable)
432    pub v_inter: f64,    // Internode voltage (mV)
433    pub c_inter: f64,    // Internode capacitance (µF/cm², very low)
434    pub g_l_myelin: f64, // Myelin leak conductance (very low)
435    pub e_l_myelin: f64, // Myelin leak reversal
436    pub g_para: f64,     // Paranodal seal conductance
437    pub dt: f64,
438    pub gain: f64,
439}
440
441impl Default for MyelinatedAxon {
442    fn default() -> Self {
443        Self::new()
444    }
445}
446
447impl MyelinatedAxon {
448    pub fn new() -> Self {
449        Self {
450            node: NodeOfRanvier::new(),
451            v_inter: -80.0,
452            c_inter: 0.001,    // Very low (myelin layers)
453            g_l_myelin: 0.001, // Very low leak through myelin
454            e_l_myelin: -80.0,
455            g_para: 0.01, // Paranodal seal conductance
456            dt: 0.5,
457            gain: 1.0,
458        }
459    }
460
461    pub fn step(&mut self, current: f64) -> i32 {
462        let input = self.gain * current;
463
464        // Paranodal coupling: node ↔ internode
465        let i_para_to_node = self.g_para * (self.v_inter - self.node.v);
466        let i_para_to_inter = self.g_para * (self.node.v - self.v_inter);
467
468        // Internode passive cable dynamics
469        let dv_inter =
470            (-self.g_l_myelin * (self.v_inter - self.e_l_myelin) + i_para_to_inter) / self.c_inter;
471        self.v_inter += self.node.dt / self.node.sub_steps as f64 * dv_inter;
472
473        // Safety bounds for internode
474        self.v_inter = self.v_inter.clamp(-120.0, 60.0);
475        if !self.v_inter.is_finite() {
476            self.v_inter = -80.0;
477        }
478
479        // Step the node with external input + paranodal current
480        // We modify the node's current to include paranodal coupling
481        let total_input = input + i_para_to_node * 100.0; // Scale for node's C_m
482        self.node.step(total_input)
483    }
484
485    /// Access the node membrane potential.
486    pub fn v(&self) -> f64 {
487        self.node.v
488    }
489
490    pub fn reset(&mut self) {
491        self.node.reset();
492        self.v_inter = -80.0;
493    }
494}
495
496#[cfg(test)]
497mod tests {
498    use super::*;
499
500    // -- Frankenhaeuser-Huxley Axon tests --
501
502    #[test]
503    fn fh_fires_with_input() {
504        // FH model uses µA/cm² — need ~1000+ for spiking (FH 1964 Fig 3)
505        let mut n = FrankenhaeUserHuxleyAxon::new();
506        let mut spikes = 0;
507        for _ in 0..2_000 {
508            spikes += n.step(2000.0);
509        }
510        assert!(
511            spikes > 0,
512            "FH axon must fire with strong input, got {spikes}"
513        );
514    }
515
516    #[test]
517    fn fh_silent_without_input() {
518        let mut n = FrankenhaeUserHuxleyAxon::new();
519        let mut spikes = 0;
520        for _ in 0..5_000 {
521            spikes += n.step(0.0);
522        }
523        assert_eq!(
524            spikes, 0,
525            "FH axon must be silent without input, got {spikes}"
526        );
527    }
528
529    #[test]
530    fn fh_action_potential_shape() {
531        // AP should depolarise well above 60 mV (spike threshold)
532        let mut n = FrankenhaeUserHuxleyAxon::new();
533        let mut v_max = -100.0_f64;
534        for _ in 0..500 {
535            n.step(2000.0);
536            v_max = v_max.max(n.v);
537        }
538        assert!(v_max > 40.0, "AP peak should exceed 40 mV, got {v_max:.1}");
539    }
540
541    #[test]
542    fn fh_gating_evolves() {
543        let mut n = FrankenhaeUserHuxleyAxon::new();
544        let m0 = n.m;
545        let h0 = n.h;
546        for _ in 0..100 {
547            n.step(2000.0);
548        }
549        assert!(n.m != m0 || n.h != h0, "Gating variables must evolve");
550    }
551
552    #[test]
553    fn fh_four_gates() {
554        // All 4 gates (m, h, n, p) must evolve during spiking
555        let mut n = FrankenhaeUserHuxleyAxon::new();
556        for _ in 0..200 {
557            n.step(2000.0);
558        }
559        // After spiking: m should have risen, h should have fallen
560        // n and p should have changed from initial
561        assert!(
562            n.m > 0.005 || n.h < 0.8 || n.n > 0.01 || n.p > 0.01,
563            "All gates must evolve: m={:.3}, h={:.3}, n={:.3}, p={:.3}",
564            n.m,
565            n.h,
566            n.n,
567            n.p
568        );
569    }
570
571    #[test]
572    fn fh_stronger_input_more_spikes() {
573        let mut weak = FrankenhaeUserHuxleyAxon::new();
574        let mut strong = FrankenhaeUserHuxleyAxon::new();
575        let (mut sw, mut ss) = (0, 0);
576        for _ in 0..2_000 {
577            sw += weak.step(1000.0);
578            ss += strong.step(3000.0);
579        }
580        assert!(
581            ss >= sw,
582            "Stronger input → more spikes: strong={ss} vs weak={sw}"
583        );
584    }
585
586    #[test]
587    fn fh_all_gates_bounded() {
588        let mut n = FrankenhaeUserHuxleyAxon::new();
589        for _ in 0..2_000 {
590            n.step(3000.0);
591        }
592        assert!(n.m >= 0.0 && n.m <= 1.0, "m out of bounds: {}", n.m);
593        assert!(n.h >= 0.0 && n.h <= 1.0, "h out of bounds: {}", n.h);
594        assert!(n.n >= 0.0 && n.n <= 1.0, "n out of bounds: {}", n.n);
595        assert!(n.p >= 0.0 && n.p <= 1.0, "p out of bounds: {}", n.p);
596    }
597
598    #[test]
599    fn fh_nan_input_stays_finite() {
600        let mut n = FrankenhaeUserHuxleyAxon::new();
601        n.step(f64::NAN);
602        assert!(n.v.is_finite());
603        assert!(n.m.is_finite());
604    }
605
606    #[test]
607    fn fh_reset_clears_state() {
608        let mut n = FrankenhaeUserHuxleyAxon::new();
609        for _ in 0..500 {
610            n.step(2000.0);
611        }
612        n.reset();
613        assert_eq!(n.v, 0.0);
614        assert_eq!(n.m, 0.005);
615        assert_eq!(n.h, 0.8);
616    }
617
618    #[test]
619    fn fh_performance_1k_steps() {
620        let start = std::time::Instant::now();
621        let mut n = FrankenhaeUserHuxleyAxon::new();
622        for _ in 0..1_000 {
623            std::hint::black_box(n.step(1500.0));
624        }
625        let elapsed = start.elapsed();
626        // 50 sub-steps per step → 50k total iterations
627        assert!(
628            elapsed.as_millis() < 100,
629            "1k steps must complete in <100ms"
630        );
631    }
632
633    // -- Node of Ranvier (MRG 2002) tests --
634
635    #[test]
636    fn nor_fires_with_input() {
637        let mut n = NodeOfRanvier::new();
638        let mut spikes = 0;
639        for _ in 0..2_000 {
640            spikes += n.step(500.0);
641        }
642        assert!(
643            spikes > 0,
644            "Node of Ranvier must fire with input, got {spikes}"
645        );
646    }
647
648    #[test]
649    fn nor_silent_without_input() {
650        let mut n = NodeOfRanvier::new();
651        let mut spikes = 0;
652        for _ in 0..5_000 {
653            spikes += n.step(0.0);
654        }
655        assert_eq!(spikes, 0, "Must be silent without input, got {spikes}");
656    }
657
658    #[test]
659    fn nor_high_nat_density() {
660        // Node of Ranvier has g_nat=3000 (much higher than standard HH ~120)
661        let n = NodeOfRanvier::new();
662        assert!(
663            n.g_nat > 1000.0,
664            "Nodal transient Na should be very high: g_nat={}",
665            n.g_nat
666        );
667    }
668
669    #[test]
670    fn nor_has_persistent_na() {
671        // MRG model must include persistent Na — distinguishes from generic HH
672        let n = NodeOfRanvier::new();
673        assert!(
674            n.g_nap > 0.0,
675            "MRG model must have persistent Na current: g_nap={}",
676            n.g_nap
677        );
678    }
679
680    #[test]
681    fn nor_has_kv7_slow_k() {
682        // Kv7 (KCNQ) is the dominant K channel at nodes, not Kv3 or Kv1
683        let n = NodeOfRanvier::new();
684        assert!(
685            n.g_ks > 0.0,
686            "MRG model must have slow K (Kv7): g_ks={}",
687            n.g_ks
688        );
689    }
690
691    #[test]
692    fn nor_persistent_na_lowers_threshold() {
693        // With persistent Na, less current is needed to fire
694        let mut with_nap = NodeOfRanvier::new();
695        let mut no_nap = NodeOfRanvier::new();
696        no_nap.g_nap = 0.0;
697        let (mut s_with, mut s_without) = (0, 0);
698        for _ in 0..2_000 {
699            s_with += with_nap.step(200.0);
700            s_without += no_nap.step(200.0);
701        }
702        assert!(
703            s_with >= s_without,
704            "Persistent Na should lower threshold: with={s_with} vs without={s_without}"
705        );
706    }
707
708    #[test]
709    fn nor_gating_evolves() {
710        let mut n = NodeOfRanvier::new();
711        let m0 = n.m;
712        let p0 = n.p;
713        for _ in 0..100 {
714            n.step(500.0);
715        }
716        assert!(
717            n.m != m0 || n.p != p0,
718            "Gating must evolve: m={:.3}, p={:.3}",
719            n.m,
720            n.p
721        );
722    }
723
724    #[test]
725    fn nor_nan_input_stays_finite() {
726        let mut n = NodeOfRanvier::new();
727        n.step(f64::NAN);
728        assert!(n.v.is_finite());
729        assert!(n.m.is_finite());
730        assert!(n.p.is_finite());
731        assert!(n.s.is_finite());
732    }
733
734    #[test]
735    fn nor_reset_clears_state() {
736        let mut n = NodeOfRanvier::new();
737        for _ in 0..500 {
738            n.step(500.0);
739        }
740        n.reset();
741        assert_eq!(n.v, -80.0);
742        assert_eq!(n.m, 0.01);
743        assert_eq!(n.p, 0.01);
744        assert_eq!(n.s, 0.05);
745    }
746
747    #[test]
748    fn nor_performance_1k_steps() {
749        let start = std::time::Instant::now();
750        let mut n = NodeOfRanvier::new();
751        for _ in 0..1_000 {
752            std::hint::black_box(n.step(500.0));
753        }
754        let elapsed = start.elapsed();
755        assert!(elapsed.as_millis() < 50, "1k steps must complete in <50ms");
756    }
757
758    // -- Myelinated Axon tests --
759
760    #[test]
761    fn myelin_fires_with_input() {
762        let mut n = MyelinatedAxon::new();
763        let mut spikes = 0;
764        for _ in 0..2_000 {
765            spikes += n.step(500.0);
766        }
767        assert!(spikes > 0, "Myelinated axon must fire, got {spikes}");
768    }
769
770    #[test]
771    fn myelin_silent_without_input() {
772        let mut n = MyelinatedAxon::new();
773        let mut spikes = 0;
774        for _ in 0..5_000 {
775            spikes += n.step(0.0);
776        }
777        assert_eq!(spikes, 0, "Must be silent without input, got {spikes}");
778    }
779
780    #[test]
781    fn myelin_internode_coupling() {
782        // Internode voltage should be affected by node spiking
783        let mut n = MyelinatedAxon::new();
784        let v_inter_0 = n.v_inter;
785        for _ in 0..500 {
786            n.step(500.0);
787        }
788        assert!(
789            (n.v_inter - v_inter_0).abs() > 0.001,
790            "Internode voltage should change with node activity: v_inter={}",
791            n.v_inter
792        );
793    }
794
795    #[test]
796    fn myelin_has_low_capacitance() {
797        let n = MyelinatedAxon::new();
798        assert!(
799            n.c_inter < 0.01,
800            "Myelin capacitance must be very low: {}",
801            n.c_inter
802        );
803    }
804
805    #[test]
806    fn myelin_has_low_myelin_leak() {
807        let n = MyelinatedAxon::new();
808        assert!(
809            n.g_l_myelin < 0.01,
810            "Myelin leak must be very low: {}",
811            n.g_l_myelin
812        );
813    }
814
815    #[test]
816    fn myelin_has_paranodal_seal() {
817        let n = MyelinatedAxon::new();
818        assert!(n.g_para > 0.0, "Must have paranodal seal conductance");
819    }
820
821    #[test]
822    fn myelin_stronger_input_more_spikes() {
823        let mut weak = MyelinatedAxon::new();
824        let mut strong = MyelinatedAxon::new();
825        let (mut sw, mut ss) = (0, 0);
826        for _ in 0..2_000 {
827            sw += weak.step(300.0);
828            ss += strong.step(1000.0);
829        }
830        assert!(ss >= sw, "Stronger → more spikes: strong={ss} vs weak={sw}");
831    }
832
833    #[test]
834    fn myelin_nan_input_stays_finite() {
835        let mut n = MyelinatedAxon::new();
836        n.step(f64::NAN);
837        assert!(n.v().is_finite());
838        assert!(n.v_inter.is_finite());
839    }
840
841    #[test]
842    fn myelin_reset_clears_state() {
843        let mut n = MyelinatedAxon::new();
844        for _ in 0..500 {
845            n.step(500.0);
846        }
847        n.reset();
848        assert_eq!(n.v_inter, -80.0);
849        assert_eq!(n.node.v, -80.0);
850    }
851
852    #[test]
853    fn myelin_performance_1k_steps() {
854        let start = std::time::Instant::now();
855        let mut n = MyelinatedAxon::new();
856        for _ in 0..1_000 {
857            std::hint::black_box(n.step(500.0));
858        }
859        let elapsed = start.elapsed();
860        assert!(elapsed.as_millis() < 50, "1k steps must complete in <50ms");
861    }
862
863    #[test]
864    fn myelinated_family_defaults_match_constructors() {
865        let fh_default = FrankenhaeUserHuxleyAxon::default();
866        let fh_constructed = FrankenhaeUserHuxleyAxon::new();
867        assert_eq!(fh_default.v, fh_constructed.v);
868        assert_eq!(fh_default.sub_steps, fh_constructed.sub_steps);
869
870        let node_default = NodeOfRanvier::default();
871        let node_constructed = NodeOfRanvier::new();
872        assert_eq!(node_default.v, node_constructed.v);
873        assert_eq!(node_default.g_nat, node_constructed.g_nat);
874
875        let axon_default = MyelinatedAxon::default();
876        let axon_constructed = MyelinatedAxon::new();
877        assert_eq!(axon_default.v_inter, axon_constructed.v_inter);
878        assert_eq!(axon_default.g_para, axon_constructed.g_para);
879    }
880
881    #[test]
882    fn fh_rate_singularities_use_finite_limits() {
883        for v in [22.0, 13.0, -10.0, 23.0, 21.0, -4.0] {
884            let mut n = FrankenhaeUserHuxleyAxon::new();
885            n.v = v;
886            n.dt = 0.0;
887            n.sub_steps = 1;
888            n.step(0.0);
889            assert!(n.v.is_finite(), "voltage must remain finite at V={v}");
890            assert!(n.m.is_finite(), "m gate must remain finite at V={v}");
891            assert!(n.h.is_finite(), "h gate must remain finite at V={v}");
892            assert!(n.n.is_finite(), "n gate must remain finite at V={v}");
893            assert!(n.p.is_finite(), "p gate must remain finite at V={v}");
894        }
895    }
896
897    #[test]
898    fn fh_nonfinite_gates_recover_to_initial_state() {
899        let mut n = FrankenhaeUserHuxleyAxon::new();
900        n.m = f64::NAN;
901        n.h = f64::NAN;
902        n.n = f64::NAN;
903        n.p = f64::NAN;
904        n.dt = 0.0;
905        n.sub_steps = 1;
906        n.step(0.0);
907        assert_eq!(n.m, 0.005);
908        assert_eq!(n.h, 0.8);
909        assert_eq!(n.n, 0.01);
910        assert_eq!(n.p, 0.01);
911    }
912
913    #[test]
914    fn myelin_nonfinite_internode_recovers_to_rest() {
915        let mut n = MyelinatedAxon::new();
916        n.v_inter = f64::NAN;
917        n.step(0.0);
918        assert_eq!(n.v_inter, -80.0);
919        assert!(n.node.v.is_finite());
920    }
921}
sc_neurocore_engine/neurons/misc/myelinated_axon.rs

sc_neurocore_engine/neurons/misc/
myelinated_axon.rs
