sc_neurocore_engine/neurons/simple_spiking/
pernarowski.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 — Pernarowski Neuron Model
8
9//! Pernarowski beta-cell bursting dynamics.
10
11/// Pernarowski 1994 — simplified beta cell burster (3 ODE).
12#[derive(Clone, Debug)]
13pub struct PernarowskiNeuron {
14    pub v: f64,
15    pub w: f64,
16    pub z: f64,
17    pub alpha: f64,
18    pub beta: f64,
19    pub eps1: f64,
20    pub eps2: f64,
21    pub gamma: f64,
22    pub dt: f64,
23    pub v_threshold: f64,
24}
25
26impl PernarowskiNeuron {
27    pub fn new() -> Self {
28        Self {
29            v: -1.0,
30            w: 0.0,
31            z: 0.0,
32            alpha: 0.1,
33            beta: 0.5,
34            eps1: 0.1,
35            eps2: 0.001,
36            gamma: 0.5,
37            dt: 0.1,
38            v_threshold: 0.5,
39        }
40    }
41    fn is_valid(&self) -> bool {
42        self.v.is_finite()
43            && self.w.is_finite()
44            && self.z.is_finite()
45            && self.alpha.is_finite()
46            && self.beta.is_finite()
47            && self.eps1.is_finite()
48            && self.eps1 > 0.0
49            && self.eps2.is_finite()
50            && self.eps2 > 0.0
51            && self.gamma.is_finite()
52            && self.gamma > 0.0
53            && self.dt.is_finite()
54            && self.dt > 0.0
55            && self.v_threshold.is_finite()
56    }
57    fn derivatives(&self, v: f64, w: f64, z: f64, current: f64) -> Option<(f64, f64, f64)> {
58        if !(v.is_finite() && w.is_finite() && z.is_finite() && current.is_finite()) {
59            return None;
60        }
61        let dv = v - v.powi(3) / 3.0 - w - z + current;
62        let dw = self.eps1 * (v - self.gamma * w + self.alpha);
63        let dz = self.eps2 * (self.beta * (v + 0.7) - z);
64        if dv.is_finite() && dw.is_finite() && dz.is_finite() {
65            Some((dv, dw, dz))
66        } else {
67            None
68        }
69    }
70    fn rk4_candidate(&self, current: f64) -> Option<(f64, f64, f64)> {
71        let dt = self.dt;
72        let (k1v, k1w, k1z) = self.derivatives(self.v, self.w, self.z, current)?;
73        let (k2v, k2w, k2z) = self.derivatives(
74            self.v + 0.5 * dt * k1v,
75            self.w + 0.5 * dt * k1w,
76            self.z + 0.5 * dt * k1z,
77            current,
78        )?;
79        let (k3v, k3w, k3z) = self.derivatives(
80            self.v + 0.5 * dt * k2v,
81            self.w + 0.5 * dt * k2w,
82            self.z + 0.5 * dt * k2z,
83            current,
84        )?;
85        let (k4v, k4w, k4z) = self.derivatives(
86            self.v + dt * k3v,
87            self.w + dt * k3w,
88            self.z + dt * k3z,
89            current,
90        )?;
91        let v = self.v + dt * (k1v + 2.0 * k2v + 2.0 * k3v + k4v) / 6.0;
92        let w = self.w + dt * (k1w + 2.0 * k2w + 2.0 * k3w + k4w) / 6.0;
93        let z = self.z + dt * (k1z + 2.0 * k2z + 2.0 * k3z + k4z) / 6.0;
94        if v.is_finite() && w.is_finite() && z.is_finite() {
95            Some((v, w, z))
96        } else {
97            None
98        }
99    }
100    pub fn step(&mut self, current: f64) -> i32 {
101        if !self.is_valid() || !current.is_finite() {
102            return 0;
103        }
104        let v_prev = self.v;
105        let Some((v, w, z)) = self.rk4_candidate(current) else {
106            return 0;
107        };
108        self.v = v;
109        self.w = w;
110        self.z = z;
111        if self.v >= self.v_threshold && v_prev < self.v_threshold {
112            1
113        } else {
114            0
115        }
116    }
117    /// Run `n_steps` RK4 updates under a constant input, returning the `v` trace
118    /// and the upward-`v_threshold`-crossing spike count. Reuses `step`, so the
119    /// trace is bit-identical to the per-step path and to the Python reference
120    /// (the cubic uses `v.powi(3)` = `v*v*v`, matching the Python `v*v*v`; no
121    /// transcendental functions). The final state is left in `self.v` / `self.w`
122    /// / `self.z`.
123    pub fn simulate(&mut self, n_steps: usize, current: f64) -> (Vec<f64>, i64) {
124        let mut trace = Vec::with_capacity(n_steps);
125        let mut spikes: i64 = 0;
126        for _ in 0..n_steps {
127            let spiked = self.step(current);
128            trace.push(self.v);
129            spikes += spiked as i64;
130        }
131        (trace, spikes)
132    }
133    pub fn reset(&mut self) {
134        self.v = -1.0;
135        self.w = 0.0;
136        self.z = 0.0;
137    }
138}
139impl Default for PernarowskiNeuron {
140    fn default() -> Self {
141        Self::new()
142    }
143}
144
145#[cfg(test)]
146mod tests {
147    use super::*;
148
149    #[test]
150    fn default_matches_constructor_state() {
151        let default = PernarowskiNeuron::default();
152        let constructed = PernarowskiNeuron::new();
153        assert_eq!(default.v, constructed.v);
154    }
155
156    #[test]
157    fn simulate_matches_repeated_step() {
158        let mut simulated = PernarowskiNeuron::new();
159        let mut repeated = PernarowskiNeuron::new();
160        let (trace, spikes) = simulated.simulate(2_000, 1.0);
161        let mut expected_trace = Vec::with_capacity(2_000);
162        let mut expected_spikes = 0_i64;
163        for _ in 0..2_000 {
164            if repeated.step(1.0) == 1 {
165                expected_spikes += 1;
166            }
167            expected_trace.push(repeated.v);
168        }
169        assert_eq!(trace, expected_trace);
170        assert_eq!(spikes, expected_spikes);
171    }
172
173    #[test]
174    fn pernarowski_fires() {
175        let mut n = PernarowskiNeuron::new();
176        let t: i32 = (0..2000).map(|_| n.step(1.0)).sum();
177        assert!(t > 0);
178    }
179
180    #[test]
181    fn pernarowski_reset_clears_state() {
182        let mut n = PernarowskiNeuron::new();
183        for _ in 0..500 {
184            n.step(1.0);
185        }
186        n.reset();
187        assert!((n.v - (-1.0)).abs() < 1e-10);
188    }
189
190    #[test]
191    fn pernarowski_bounded() {
192        let mut n = PernarowskiNeuron::new();
193        for _ in 0..2000 {
194            n.step(50.0);
195        }
196        assert!(n.v.is_finite());
197    }
198
199    #[test]
200    fn pernarowski_slow_z() {
201        let mut n = PernarowskiNeuron::new();
202        let z0 = n.z;
203        for _ in 0..2000 {
204            n.step(1.0);
205        }
206        assert!((n.z - z0).abs() > 1e-6, "slow z should evolve");
207    }
208
209    #[test]
210    fn pernarowski_matches_rk4_candidate() {
211        let mut n = PernarowskiNeuron::new();
212        n.v = -0.8;
213        n.w = 0.2;
214        n.z = -0.1;
215        let current = 0.5;
216        let dt = n.dt;
217        let rhs = |v: f64, w: f64, z: f64| {
218            (
219                v - v.powi(3) / 3.0 - w - z + current,
220                n.eps1 * (v - n.gamma * w + n.alpha),
221                n.eps2 * (n.beta * (v + 0.7) - z),
222            )
223        };
224        let (k1v, k1w, k1z) = rhs(n.v, n.w, n.z);
225        let (k2v, k2w, k2z) = rhs(
226            n.v + 0.5 * dt * k1v,
227            n.w + 0.5 * dt * k1w,
228            n.z + 0.5 * dt * k1z,
229        );
230        let (k3v, k3w, k3z) = rhs(
231            n.v + 0.5 * dt * k2v,
232            n.w + 0.5 * dt * k2w,
233            n.z + 0.5 * dt * k2z,
234        );
235        let (k4v, k4w, k4z) = rhs(n.v + dt * k3v, n.w + dt * k3w, n.z + dt * k3z);
236        let expected = (
237            n.v + dt * (k1v + 2.0 * k2v + 2.0 * k3v + k4v) / 6.0,
238            n.w + dt * (k1w + 2.0 * k2w + 2.0 * k3w + k4w) / 6.0,
239            n.z + dt * (k1z + 2.0 * k2z + 2.0 * k3z + k4z) / 6.0,
240        );
241
242        assert_eq!(n.step(current), 0);
243        assert!((n.v - expected.0).abs() < 1e-14);
244        assert!((n.w - expected.1).abs() < 1e-14);
245        assert!((n.z - expected.2).abs() < 1e-14);
246    }
247
248    #[test]
249    fn pernarowski_invalid_input_preserves_state() {
250        let mut n = PernarowskiNeuron::new();
251        let before = (n.v, n.w, n.z);
252        assert_eq!(n.step(f64::NAN), 0);
253        assert_eq!((n.v, n.w, n.z), before);
254    }
255
256    #[test]
257    fn pernarowski_overflow_candidate_preserves_state() {
258        let mut n = PernarowskiNeuron::new();
259        n.v = 1.0e160;
260        let before = (n.v, n.w, n.z);
261        assert_eq!(n.step(0.5), 0);
262        assert_eq!((n.v, n.w, n.z), before);
263    }
264
265    #[test]
266    fn pernarowski_nan_no_panic() {
267        PernarowskiNeuron::new().step(f64::NAN);
268    }
269
270    #[test]
271    fn pernarowski_negative_no_crash() {
272        let mut n = PernarowskiNeuron::new();
273        for _ in 0..500 {
274            n.step(-5.0);
275        }
276        assert!(n.v.is_finite());
277    }
278}
sc_neurocore_engine/neurons/simple_spiking/pernarowski.rs

sc_neurocore_engine/neurons/simple_spiking/
pernarowski.rs
