sc_neurocore_engine/neurons/biophysical/
prescott.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 — Prescott Neuron Model
8
9//! Prescott conductance model for Type I, II, and III excitability.
10
11/// Prescott 2008 — Type I/II/III excitability tuning via M-current.
12#[derive(Clone, Debug)]
13pub struct PrescottNeuron {
14    pub v: f64,
15    pub w: f64,
16    pub g_fast: f64,
17    pub g_slow: f64,
18    pub g_l: f64,
19    pub e_fast: f64,
20    pub e_slow: f64,
21    pub e_l: f64,
22    pub beta_w: f64,
23    pub gamma_w: f64,
24    pub tau_w: f64,
25    pub phi: f64,
26    pub dt: f64,
27    pub v_threshold: f64,
28}
29
30impl PrescottNeuron {
31    pub fn new() -> Self {
32        Self {
33            v: -65.0,
34            w: 0.0,
35            g_fast: 20.0,
36            g_slow: 20.0,
37            g_l: 2.0,
38            e_fast: 50.0,
39            e_slow: -100.0,
40            e_l: -70.0,
41            beta_w: -21.0,
42            gamma_w: 15.0,
43            tau_w: 100.0,
44            phi: 0.15,
45            dt: 0.1,
46            v_threshold: -20.0,
47        }
48    }
49
50    fn sigmoid(x: f64) -> f64 {
51        if x >= 0.0 {
52            let z = (-x).exp();
53            1.0 / (1.0 + z)
54        } else {
55            let z = x.exp();
56            z / (1.0 + z)
57        }
58    }
59
60    fn valid_state(v: f64, w: f64) -> bool {
61        v.is_finite() && w.is_finite() && (0.0..=1.0).contains(&w)
62    }
63
64    fn valid_runtime(&self) -> bool {
65        Self::valid_state(self.v, self.w)
66            && self.g_fast.is_finite()
67            && self.g_fast >= 0.0
68            && self.g_slow.is_finite()
69            && self.g_slow >= 0.0
70            && self.g_l.is_finite()
71            && self.g_l >= 0.0
72            && self.e_fast.is_finite()
73            && self.e_slow.is_finite()
74            && self.e_l.is_finite()
75            && self.beta_w.is_finite()
76            && self.gamma_w.is_finite()
77            && self.gamma_w > 0.0
78            && self.tau_w.is_finite()
79            && self.tau_w > 0.0
80            && self.phi.is_finite()
81            && self.phi >= 0.0
82            && self.dt.is_finite()
83            && self.dt > 0.0
84            && self.v_threshold.is_finite()
85    }
86
87    fn derivatives(&self, v: f64, w: f64, current: f64) -> Option<(f64, f64)> {
88        if !Self::valid_state(v, w) {
89            return None;
90        }
91        let m_inf = Self::sigmoid((v + 20.0) / 15.0);
92        let w_inf = Self::sigmoid((v - self.beta_w) / self.gamma_w);
93        let i_fast = self.g_fast * m_inf * (v - self.e_fast);
94        let i_slow = self.g_slow * w * (v - self.e_slow);
95        let i_l = self.g_l * (v - self.e_l);
96        let dv = -i_fast - i_slow - i_l + current;
97        let dw = self.phi * (w_inf - w) / self.tau_w;
98        if dv.is_finite() && dw.is_finite() {
99            Some((dv, dw))
100        } else {
101            None
102        }
103    }
104
105    fn rk4_step(&self, current: f64) -> Option<(f64, f64)> {
106        let dt = self.dt;
107        let (k1_v, k1_w) = self.derivatives(self.v, self.w, current)?;
108        let (k2_v, k2_w) =
109            self.derivatives(self.v + 0.5 * dt * k1_v, self.w + 0.5 * dt * k1_w, current)?;
110        let (k3_v, k3_w) =
111            self.derivatives(self.v + 0.5 * dt * k2_v, self.w + 0.5 * dt * k2_w, current)?;
112        let (k4_v, k4_w) = self.derivatives(self.v + dt * k3_v, self.w + dt * k3_w, current)?;
113        let next_v = self.v + dt * (k1_v + 2.0 * k2_v + 2.0 * k3_v + k4_v) / 6.0;
114        let next_w = self.w + dt * (k1_w + 2.0 * k2_w + 2.0 * k3_w + k4_w) / 6.0;
115        if Self::valid_state(next_v, next_w) {
116            Some((next_v, next_w))
117        } else {
118            None
119        }
120    }
121
122    pub fn step(&mut self, current: f64) -> i32 {
123        if !current.is_finite() || !self.valid_runtime() {
124            return 0;
125        }
126        let v_prev = self.v;
127        let Some((next_v, next_w)) = self.rk4_step(current) else {
128            return 0;
129        };
130        self.v = next_v;
131        self.w = next_w;
132        if self.v >= self.v_threshold && v_prev < self.v_threshold {
133            1
134        } else {
135            0
136        }
137    }
138
139    pub fn reset(&mut self) {
140        self.v = -65.0;
141        self.w = 0.0;
142    }
143}
144impl Default for PrescottNeuron {
145    fn default() -> Self {
146        Self::new()
147    }
148}
149
150#[cfg(test)]
151mod tests {
152    use super::*;
153
154    #[test]
155    fn default_matches_constructor_state() {
156        let default = PrescottNeuron::default();
157        let constructed = PrescottNeuron::new();
158        assert_eq!(default.v, constructed.v);
159    }
160
161    #[test]
162    fn derivative_rejects_invalid_and_nonfinite_candidates() {
163        let mut n = PrescottNeuron::new();
164        assert_eq!(n.derivatives(n.v, 2.0, 0.0), None);
165        n.g_fast = f64::MAX;
166        assert_eq!(n.derivatives(f64::MAX, 0.5, 0.0), None);
167    }
168
169    #[test]
170    fn invalid_rk4_candidate_preserves_state() {
171        let mut n = PrescottNeuron::new();
172        n.dt = 1.0e-300;
173        let before = (n.v, n.w);
174        assert_eq!(n.step(f64::MAX / 2.0), 0);
175        assert_eq!((n.v, n.w), before);
176    }
177
178    #[test]
179    fn prescott_fires() {
180        let mut n = PrescottNeuron::new();
181        let t: i32 = (0..500).map(|_| n.step(5.0)).sum();
182        assert!(t > 0);
183    }
184
185    // -- Prescott --
186    #[test]
187    fn prescott_zero_input_stable() {
188        let mut n = PrescottNeuron::new();
189        let _t: i32 = (0..500).map(|_| n.step(0.0)).sum();
190        // Prescott has fast Na conductance — may produce spontaneous activity
191        assert!(n.v.is_finite());
192    }
193    #[test]
194    fn prescott_reset_clears_state() {
195        let mut n = PrescottNeuron::new();
196        for _ in 0..100 {
197            n.step(5.0);
198        }
199        n.reset();
200        assert!((n.v - (-65.0)).abs() < 1e-10);
201        assert!((n.w - 0.0).abs() < 1e-10);
202    }
203    #[test]
204    fn prescott_rk4_reference_point() {
205        let mut n = PrescottNeuron::new();
206        assert_eq!(n.step(50.0), 0);
207        assert!((n.v - (-44.498914201492525)).abs() < 1e-12);
208        assert!((n.w - 1.4035864179018786e-05).abs() < 1e-17);
209    }
210    #[test]
211    fn prescott_extreme_bounded() {
212        let mut n = PrescottNeuron::new();
213        for _ in 0..200 {
214            n.step(1e4);
215        }
216        assert!(n.v.is_finite());
217    }
218    #[test]
219    fn prescott_slow_var_adapts() {
220        let mut n = PrescottNeuron::new();
221        for _ in 0..500 {
222            n.step(5.0);
223        }
224        assert!(n.w > 0.0, "slow variable w should activate during spiking");
225    }
226    #[test]
227    fn prescott_negative_no_crash() {
228        let mut n = PrescottNeuron::new();
229        for _ in 0..200 {
230            n.step(-10.0);
231        }
232        assert!(n.v.is_finite());
233    }
234    #[test]
235    fn prescott_nan_no_panic() {
236        let mut n = PrescottNeuron::new();
237        n.step(f64::NAN);
238    }
239}
sc_neurocore_engine/neurons/biophysical/prescott.rs

sc_neurocore_engine/neurons/biophysical/
prescott.rs
