sc_neurocore_engine/neurons/biophysical/
pospischil.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 — Pospischil Neuron Model
8
9//! Pospischil minimal conductance model for cortical cell types.
10
11/// Pospischil — minimal HH for 5 cortical cell types. Pospischil et al. 2008.
12#[derive(Clone, Debug)]
13pub struct PospischilNeuron {
14    pub v: f64,
15    pub m: f64,
16    pub h: f64,
17    pub n: f64,
18    pub p: f64,
19    pub g_na: f64,
20    pub g_k: f64,
21    pub g_m: f64,
22    pub g_l: f64,
23    pub e_na: f64,
24    pub e_k: f64,
25    pub e_l: f64,
26    pub c_m: f64,
27    pub vt: f64,
28    pub dt: f64,
29    pub v_threshold: f64,
30}
31
32impl PospischilNeuron {
33    pub fn new() -> Self {
34        Self {
35            v: -70.0,
36            m: 0.05,
37            h: 0.6,
38            n: 0.3,
39            p: 0.0,
40            g_na: 50.0,
41            g_k: 5.0,
42            g_m: 0.07,
43            g_l: 0.1,
44            e_na: 50.0,
45            e_k: -90.0,
46            e_l: -70.0,
47            c_m: 1.0,
48            vt: -56.2,
49            dt: 0.025,
50            v_threshold: -20.0,
51        }
52    }
53    /// Return `[dV, dm, dh, dn, dp]` of the five-state system at one consistent
54    /// state. The Traub-Miles activation rates use the closed-form L'Hôpital limit
55    /// within `1e-6` of their `x/(exp(±x/k)-1)` removable singularities, matching
56    /// the Python/Julia/Go/Mojo kernels.
57    fn derivatives(&self, v: f64, m: f64, h: f64, n: f64, p: f64, current: f64) -> [f64; 5] {
58        let dv_vt = v - self.vt;
59        let x_m = dv_vt - 13.0;
60        let am = if x_m.abs() < 1e-6 {
61            1.28
62        } else {
63            -0.32 * x_m / ((-(x_m) / 4.0).exp() - 1.0)
64        };
65        let x_bm = dv_vt - 40.0;
66        let bm = if x_bm.abs() < 1e-6 {
67            1.4
68        } else {
69            0.28 * x_bm / ((x_bm / 5.0).exp() - 1.0)
70        };
71        let ah = 0.128 * (-(dv_vt - 17.0) / 18.0).exp();
72        let bh = 4.0 / (1.0 + (-(dv_vt - 40.0) / 5.0).exp());
73        let x_n = dv_vt - 15.0;
74        let an = if x_n.abs() < 1e-6 {
75            0.16
76        } else {
77            -0.032 * x_n / ((-(x_n) / 5.0).exp() - 1.0)
78        };
79        let bn = 0.5 * (-(dv_vt - 10.0) / 40.0).exp();
80        let p_inf = 1.0 / (1.0 + (-(v + 35.0) / 10.0).exp());
81        let tau_p = 608.0 / (3.3 * ((v + 35.0) / 20.0).exp() + (-(v + 35.0) / 20.0).exp());
82        let dm = am * (1.0 - m) - bm * m;
83        let dh = ah * (1.0 - h) - bh * h;
84        let dn = an * (1.0 - n) - bn * n;
85        let dp = (p_inf - p) / tau_p;
86        let i_na = self.g_na * m * m * m * h * (v - self.e_na);
87        let i_k = self.g_k * n * n * n * n * (v - self.e_k);
88        let i_m = self.g_m * p * (v - self.e_k);
89        let i_l = self.g_l * (v - self.e_l);
90        let dv = (-i_na - i_k - i_m - i_l + current) / self.c_m;
91        [dv, dm, dh, dn, dp]
92    }
93
94    /// Return one classical RK4 increment of `[V, m, h, n, p]`, holding `current`
95    /// constant across the four stages.
96    fn rk4_substep(&self, s: [f64; 5], current: f64) -> [f64; 5] {
97        let dt = self.dt;
98        let k1 = self.derivatives(s[0], s[1], s[2], s[3], s[4], current);
99        let k2 = self.derivatives(
100            s[0] + 0.5 * dt * k1[0],
101            s[1] + 0.5 * dt * k1[1],
102            s[2] + 0.5 * dt * k1[2],
103            s[3] + 0.5 * dt * k1[3],
104            s[4] + 0.5 * dt * k1[4],
105            current,
106        );
107        let k3 = self.derivatives(
108            s[0] + 0.5 * dt * k2[0],
109            s[1] + 0.5 * dt * k2[1],
110            s[2] + 0.5 * dt * k2[2],
111            s[3] + 0.5 * dt * k2[3],
112            s[4] + 0.5 * dt * k2[4],
113            current,
114        );
115        let k4 = self.derivatives(
116            s[0] + dt * k3[0],
117            s[1] + dt * k3[1],
118            s[2] + dt * k3[2],
119            s[3] + dt * k3[3],
120            s[4] + dt * k3[4],
121            current,
122        );
123        let mut out = [0.0_f64; 5];
124        for i in 0..5 {
125            out[i] = s[i] + dt * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]) / 6.0;
126        }
127        out
128    }
129
130    pub fn step(&mut self, current: f64) -> i32 {
131        let v_prev = self.v;
132        let mut s = [self.v, self.m, self.h, self.n, self.p];
133        for _ in 0..4 {
134            s = self.rk4_substep(s, current);
135        }
136        self.v = s[0];
137        self.m = s[1];
138        self.h = s[2];
139        self.n = s[3];
140        self.p = s[4];
141        if self.v >= self.v_threshold && v_prev < self.v_threshold {
142            1
143        } else {
144            0
145        }
146    }
147    pub fn reset(&mut self) {
148        self.v = -70.0;
149        self.m = 0.05;
150        self.h = 0.6;
151        self.n = 0.3;
152        self.p = 0.0;
153    }
154}
155impl Default for PospischilNeuron {
156    fn default() -> Self {
157        Self::new()
158    }
159}
160
161#[cfg(test)]
162mod tests {
163    use super::*;
164
165    #[test]
166    fn default_matches_constructor_state() {
167        let default = PospischilNeuron::default();
168        let constructed = PospischilNeuron::new();
169        assert_eq!(default.v, constructed.v);
170    }
171
172    #[test]
173    fn removable_rate_singularities_use_finite_limits() {
174        let n = PospischilNeuron::new();
175        for voltage in [n.vt + 13.0, n.vt + 40.0, n.vt + 15.0] {
176            assert!(n
177                .derivatives(voltage, n.m, n.h, n.n, n.p, 0.0)
178                .iter()
179                .all(|value| value.is_finite()));
180        }
181    }
182
183    #[test]
184    fn pospischil_fires() {
185        let mut n = PospischilNeuron::new();
186        let t: i32 = (0..200).map(|_| n.step(5.0)).sum();
187        assert!(t > 0);
188    }
189
190    // -- Pospischil --
191    #[test]
192    fn pospischil_silent_without_input() {
193        let mut n = PospischilNeuron::new();
194        let t: i32 = (0..200).map(|_| n.step(0.0)).sum();
195        assert_eq!(t, 0);
196    }
197    #[test]
198    fn pospischil_reset_clears_state() {
199        let mut n = PospischilNeuron::new();
200        for _ in 0..100 {
201            n.step(5.0);
202        }
203        n.reset();
204        assert!((n.v - (-70.0)).abs() < 1e-10);
205    }
206    #[test]
207    fn pospischil_moderate_input_stable() {
208        let mut n = PospischilNeuron::new();
209        for _ in 0..200 {
210            n.step(10.0);
211        }
212        assert!(n.v.is_finite());
213    }
214    #[test]
215    fn pospischil_m_current_present() {
216        let mut n = PospischilNeuron::new();
217        for _ in 0..200 {
218            n.step(5.0);
219        }
220        assert!(n.p > 0.0, "M-current (p) should activate during spiking");
221    }
222    #[test]
223    fn pospischil_negative_no_crash() {
224        let mut n = PospischilNeuron::new();
225        for _ in 0..200 {
226            n.step(-10.0);
227        }
228        assert!(n.v.is_finite());
229    }
230    #[test]
231    fn pospischil_nan_no_panic() {
232        let mut n = PospischilNeuron::new();
233        n.step(f64::NAN);
234    }
235}
sc_neurocore_engine/neurons/biophysical/pospischil.rs

sc_neurocore_engine/neurons/biophysical/
pospischil.rs
